還有辦法再加速嗎?

[DennisLiu-elogic]

Hi meiqua
統計了各段落與各層金字塔的耗時

使用OpenCL後已經讓OpenCV的函數耗時下降不少，關於加速ComputeResponseMap、Hysteresis與Linearlize這三段，不知道meiqua你還有沒有什麼想法?

目前實作過這個方法，ComputeResponseMap中Create 8個Mat，改為一次Create大張的Mat
(size = Size (src.cols * 8, src.rows))，response_maps[i] = BigSrc (Rect)指向大張的ROI，但實作後發現__m128i *map_data = response_maps[ori].ptr<__m128i>(); 這行似乎不正常，導致__m128i res1 = _mm_shuffle_epi8(lut[2 * ori + 0], lsb4_data[i]);出錯

[meiqua]

可以参考这里

[DennisLiu-elogic]

使用OpenCL後，sobel耗時占比約4%而已，這個github主要加速sobel應該沒太大用處?

還是ComputeResponseMap、Hysteresis與Linearlize這三段有辦法用裡面的東西加速

[meiqua]

opencl效果这么明显应该是用到并行了。
不过我说的这个不是单独针对某个算子，而是把所有算子fuse起来，这样只需要从主存中读一次数据。

[DennisLiu-elogic]

所有算子fuse起来->>是指sobel、phase、magnitude一起計算嗎?

[meiqua]

是的，因为都是局部读一块进去算，所以可以合并一下

[DennisLiu-elogic]

能提示下用halide如何實現嗎?

[meiqua]

可以参考halide github上的教程，讲的挺详细的

[zhengli233]

因为计算phase前总要计算gaussian和sobel，所以按照fuse的思路，我尝试把gaussian算子和sobel算子先算卷积，得到融合的卷积核，再用它对原图进行卷积。

Mat src = imread("Hydrangeas.jpg", CV_8UC1);
int g_k_size = 3;
Mat sobel_x, sobel_y, phase_;
Mat kernelX = (Mat_<float>(3, 3) << -1, 0, 1, -2, 0, 2, -1, 0, 1);
Mat kernelY = (Mat_<float>(3, 3) << -1, -2, -1, 0, 0, 0, 1, 2, 1);
Mat g_kernel = getGaussianKernel(g_k_size, 1, CV_32F);
Mat g_kernel_2d = g_kernel * g_kernel.t();
Mat fused_kernel;
Mat blur_;
src.convertTo(src, CV_32F);
GaussianBlur(kernelX, fused_kernel, Size(g_k_size, g_k_size), 0, 0, BORDER_REPLICATE);
filter2D(src, sobel_x, CV_32F, fused_kernel, Point(-1, -1), 0, BORDER_REPLICATE);
Mat cv_sobel_x, cv_sobel_y, cv_phase, cv_blur;
GaussianBlur(src, cv_blur, Size(g_k_size, g_k_size), 1, 0, BORDER_REPLICATE);
Sobel(cv_blur, cv_sobel_x, CV_32F, 1, 0, 3, 1.0, 0.0, BORDER_REPLICATE);
//Sobel(src, cv_sobel_y, CV_8U, 0, 1);
//phase(cv_sobel_x, cv_sobel_y, cv_phase, true);

for (int row = 0; row < src.rows; row++)
{
	for (int col = 0; col < src.cols; col++)
	{
		float fused = *sobel_x.ptr<float>(row, col);
		float cv = *cv_sobel_x.ptr<float>(row, col);
		if (fused != cv)
		{
			cout << fused << "    " << cv << endl;
		}
	}
}

理论上来说，卷积乘法是可以先算后面的，再算前面的，但是我这个结果和原结果总是有偏差，不知道为什么。
@meiqua 你能帮我看看吗

[meiqua]

GaussianBlur(kernelX, fused_kernel, Size(g_k_size, g_k_size), 0, 0, BORDER_REPLICATE);

这一句有些不对，想得到fused_kernel直接乘起来就好。
另外最好弄成seperatable filter的形式，两个维度分开算要快点。

[zhengli233]

直接乘起来是什么意思？

[meiqua]

把两个kernel逐元素相乘，用GaussianBlur在边界的地方有问题

[zhengli233]

逐元素相乘是指矩阵点乘吗？
但src * gaussian_kernel * sobel_kernel不都是卷积运算吗，而且两个kernel的大小不一定一样啊。

[meiqua]

不是点乘，直接对应位置相乘；
是的，所以能两个kernel相乘，大小不一样可以补0

[DennisLiu-elogic]

對應位置相乘不就是點乘嗎?

我改了一下@zhengli233的程式碼

1. 產生Gaussian (7x1)與(1*7)，矩陣相乘得到Gaussian(7x7)
2. 產生X方向Sobel，補0成7x7矩陣
3. convolve Sobel 和 Gaussian成fusedKernel
4. src 套用fusedKernel

還是有誤差 (max ~= 10)，有那裡錯誤嗎?懇請大神賜教 (抱歉不知道怎麼貼程式碼

    Mat matSrc = imread ("", IMREAD_GRAYSCALE);
Mat matSobelKernelX = (Mat_<float> (3, 3) << -1, 0, 1, -2, 0, 2, -1, 0, 1), matSobelFlip;
Mat matGaussianKernelX = getGaussianKernel (7, 1, CV_32F);//7x1 matrix
Mat matGaussianKennelY = matGaussianKernelX.t ();
Mat matGaussianKernel = matGaussianKernelX * matGaussianKennelY;// make 7x7 Gaussian

Mat matSobelBorder (7, 7, CV_32F, Scalar::all (0.0));
Rect rectCenter (2, 2, 3, 3);
matSobelKernelX.copyTo (matSobelBorder (rectCenter));//make 7x7 SobelX with 0 borders
flip (matSobelBorder, matSobelFlip, -1);//x, y flip

Mat matFusedKernel, matResult1;
//convolve matGaussianFlip and matSobelBorder
filter2D (matGaussianKernel, matFusedKernel, CV_32F, matSobelFlip , Point (-1, -1), 0, BORDER_REFLECT);
//Result
filter2D (matSrc, matResult1, CV_32F, matFusedKernel, Point (-1, -1), 0.0, BORDER_REFLECT);//convolve src and fused kernel

//opencv api result
Mat cv_sobel_x, cv_sobel_y, cv_phase, cv_blur;
GaussianBlur (matSrc, cv_blur, Size (7, 7), 1, 0, BORDER_REPLICATE);
Sobel (cv_blur, cv_sobel_x, CV_32F, 1, 0, 3, 1.0, 0.0, BORDER_REPLICATE);


Mat matDiff = matResult1 - cv_sobel_x;
matDiff.convertTo (matDiff, CV_8U);
imshow ("Diff", matDiff);
double dMax;
minMaxLoc (matDiff, 0, &dMax, 0, 0);

[zhengli233]

@DennisLiu-elogic 你好，
看了你的代码，发现和我的代码最大的差别是你有一个flip的操作（另外就是sobel算子边缘补全）。这个flip应该就是相当于矩阵转置的操作吧。这里我不太明白。。虽然从结果上来说，你的比我好得多，能说说为什么要这么做吗？
另外，误差的原因我想可能就是在边缘上。我从知乎上找了一张图，介绍怎么把卷积转换成矩阵乘法。

可以看到，如果按正常顺序计算卷积，我们需要

1. 把原图展开；
2. 与卷积核做点乘；
3. 把点乘结果展开；
4. 与下一个卷积核做点乘。

这里实际上对原图展开了2次，而如果先计算后面2个卷积核，虽然也有2次展开，但是其中一次是对卷积核做的展开，因此丢失了一部分数据。