還有辦法再加速嗎?

[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次展开，但是其中一次是对卷积核做的展开，因此丢失了一部分数据。

[meiqua]

哦，之前我理解的有问题，fusion说的不是把卷积核都合成一个，因为这样会越来越大，比如两个3x3的合起来就变成5x5，虽然方法不明，但肯定不是直接一个filter2D另一个。
fusion说的是局部读一块进去，一起算完再写出来，关键是提高了局部性，而且kernel不仅不应该合起来，还应该尽量拆开算，比如3x3算两次有mn*18次乘法，5x5一次就mn*25次。这也是为什么可分离卷积分开算要快点。

[DennisLiu-elogic]

@zhengli233
根據這篇的回覆
做兩次卷積src * Sobel * Gaussian = src * (Sobel * Gaussain) = src * (Sobel convolve Gassian)
opencv Filter2D做的是correlation，所以其中一個kernel要XY反轉再帶入Filter2D

但根據版主開示"提高局部性"，這種方法應該不適用於此情況，然後我的算法在邊緣上還是有一些誤差，我猜是opencv Sobel的kernel可能不是 -1, 0, 1, -2, 0, 2, -1, 0, 1
用getDerivKenels拿Mat出來看的確是 -1, 0, 1

[meiqua]

卷积应该不满足交换律，从感受野来看，7x7 3x3卷积两次一定会扩大，变成(1+3)*2+1=9*9。
但这引入了一个有意思的问题，可能是高斯核边界的地方本来值就小，所以忽略结合后内核最外面一圈影响不大。也就是说这样近似反而减小了9/(49+9)的运算量，是不是直接内核卷积近似最好还需要推导。

[DennisLiu-elogic]

google查了前幾筆資料，都說捲積是有結合率的(a * b) * c = a * (b * c)

[meiqua]

推导了下，还真满足结合律，不过扩大成9x9也是应该的，改了下一点误差都没有了：

    Mat matSrc = imread ("", IMREAD_GRAYSCALE);

    // modification 1: double type, otherwise GaussianBlur will output CV_8U type
    matSrc.convertTo(matSrc, CV_64F);

    Mat matSobelKernelX = (Mat_<double> (3, 3) << -1, 0, 1, -2, 0, 2, -1, 0, 1), matSobelFlip;
    Mat matGaussianKernelX = getGaussianKernel (7, 1, CV_64F);//7x1 matrix
    Mat matGaussianKennelY = matGaussianKernelX.t ();
    Mat matGaussianKernel = matGaussianKernelX * matGaussianKennelY;// make 7x7 Gaussian

    // modification 2: enlarge fused kernel
    Mat matGaussianKernel9(9, 9, CV_64F, Scalar::all (0.0));
    matGaussianKernel.copyTo(matGaussianKernel9(Rect(1, 1, 7, 7)));

    flip(matSobelKernelX, matSobelFlip, -1);//x, y flip

    Mat matFusedKernel, matResult1;

    // modification 3: BORDER_CONSTANT, padding by 0.
    // Padding sobel kernel happens to have the same effect.
    filter2D (matGaussianKernel9, matFusedKernel, CV_64F, matSobelFlip , Point (-1, -1), 0, BORDER_CONSTANT);

    //Result
    filter2D (matSrc, matResult1, CV_64F, matFusedKernel, Point (-1, -1), 0.0, BORDER_REPLICATE);

    //opencv api result
    Mat cv_sobel_x, cv_blur;

    // BORDER_REPLICATE twice make it different from fused kernel around border
    GaussianBlur(matSrc, cv_blur, Size (7, 7), 1, 0, BORDER_REPLICATE);
    Sobel (cv_blur, cv_sobel_x, CV_64F, 1, 0, 3, 1.0, 0.0, BORDER_REPLICATE);


    Mat matDiff = cv::abs(matResult1 - cv_sobel_x);
//    matDiff.convertTo (matDiff, CV_8U);
    imshow ("Diff", matDiff > 1e-12);
    // no difference at all except border (limits of precision: float is 1e-3, double is 1e-12)
    // original case has small difference because it's not 9x9 kernel,
    // which means border of fused kernel are set to 0

    double dMax;
    minMaxLoc (matDiff, 0, &dMax, 0, 0);
    std::cout << "dmax: " << dMax << endl;

[zhengli233]

哦，之前我理解的有问题，fusion说的不是把卷积核都合成一个，因为这样会越来越大，比如两个3x3的合起来就变成5x5，虽然方法不明，但肯定不是直接一个filter2D另一个。
fusion说的是局部读一块进去，一起算完再写出来，关键是提高了局部性，而且kernel不仅不应该合起来，还应该尽量拆开算，比如3x3算两次有mn18次乘法，5x5一次就mn25次。这也是为什么可分离卷积分开算要快点。

@meiqua 好像确实是这么个道理，kernel不应该合起来。

第二段中你说的提高局部性是指节省了写内存的开销吗？在计算量没有减少的情况下，写内存的开销与计算比起来应该占比很小吧。

推导了下，还真满足结合律，不过扩大成9x9也是应该的，改了下一点误差都没有了：

    Mat matSrc = imread ("", IMREAD_GRAYSCALE);

    // modification 1: double type, otherwise GaussianBlur will output CV_8U type
    matSrc.convertTo(matSrc, CV_64F);

    Mat matSobelKernelX = (Mat_<double> (3, 3) << -1, 0, 1, -2, 0, 2, -1, 0, 1), matSobelFlip;
    Mat matGaussianKernelX = getGaussianKernel (7, 1, CV_64F);//7x1 matrix
    Mat matGaussianKennelY = matGaussianKernelX.t ();
    Mat matGaussianKernel = matGaussianKernelX * matGaussianKennelY;// make 7x7 Gaussian

    // modification 2: enlarge fused kernel
    Mat matGaussianKernel9(9, 9, CV_64F, Scalar::all (0.0));
    matGaussianKernel.copyTo(matGaussianKernel9(Rect(1, 1, 7, 7)));

    flip(matSobelKernelX, matSobelFlip, -1);//x, y flip

    Mat matFusedKernel, matResult1;

    // modification 3: BORDER_CONSTANT, padding by 0.
    // Padding sobel kernel happens to have the same effect.
    filter2D (matGaussianKernel9, matFusedKernel, CV_64F, matSobelFlip , Point (-1, -1), 0, BORDER_CONSTANT);

    //Result
    filter2D (matSrc, matResult1, CV_64F, matFusedKernel, Point (-1, -1), 0.0, BORDER_REPLICATE);

    //opencv api result
    Mat cv_sobel_x, cv_blur;

    // BORDER_REPLICATE twice make it different from fused kernel around border
    GaussianBlur(matSrc, cv_blur, Size (7, 7), 1, 0, BORDER_REPLICATE);
    Sobel (cv_blur, cv_sobel_x, CV_64F, 1, 0, 3, 1.0, 0.0, BORDER_REPLICATE);


    Mat matDiff = cv::abs(matResult1 - cv_sobel_x);
//    matDiff.convertTo (matDiff, CV_8U);
    imshow ("Diff", matDiff > 1e-12);
    // no difference at all except border (limits of precision: float is 1e-3, double is 1e-12)
    // original case has small difference because it's not 9x9 kernel,
    // which means border of fused kernel are set to 0

    double dMax;
    minMaxLoc (matDiff, 0, &dMax, 0, 0);
    std::cout << "dmax: " << dMax << endl;

我跑了一下这个代码，并不是没有偏差啊，dMax最大有74多，我读的win7示例图片中八仙花。

[xinsuinizhuan]

@meiqua 好像确实是这么个道理，kernel不应该合起来。

第二段中你说的提高局部性是指节省了写内存的开销吗？在计算量没有减少的情况下，写内存的开销与计算比起来应该占比很小吧。

您这是在做pyrDown的融合吗？结果怎么样？有没有实质性的效果？

[meiqua]

@zhengli233 是吗，我跑的图片边界地方是4，imshow出来中间全黑，说明差异比1e-12还小。边界地方特别大也是有可能的，因为BORDER_REPLICATE两次肯定跟一次不一样

[zhengli233]

@meiqua
没有修改你的代码，只是读了我本机路径的图片，结果是这样的

可以看到，偏差还是挺大的 : (

[zhengli233]

@xinsuinizhuan 不是在做pyrDown的融合，是高斯和sobel的融合

[xinsuinizhuan]

@xinsuinizhuan 不是在做pyrDown的融合，是高斯和sobel的融合

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

[zhengli233]

@xinsuinizhuan phase, magnitude不是卷积计算，所以也就没有算子，谈不上融合

[xinsuinizhuan]

@zhengli233 见笑了。只是现在的速度黑白图也只有90ms左右，达不到要求，在寻求加速。看到说融合后可以达到halcon的水平，您又在做类似融合的东西。

[meiqua]

@zhengli233 图片贴一下我试试？我看有可能是1e-12的阈值太小了

[zhengli233]

@meiqua 就是这张图，你试试看

[meiqua]

@zhengli233 试过了，中间还是全黑呀。dMax确实有74.2986

[zhengli233]

@meiqua
我把1e-12的阈值调大了点到1e-8，这个精度应该足够了，确实中间是全黑的了。
然后我测试了一下速度，99的fusedKernel比77 + 3*3的kernels慢了接近3倍。。怀疑opencv的高斯和sobel函数内部都是将kernel拆开算的。
所以回到最初所讨论的：

局部一块读进去算

它是怎么缩短运行时间的？

[meiqua]

对，opencv里面算高斯 sobel是用可分离卷积算的。
fusion这个问题halide talk里面讲的非常好。

[zhengli233]

@meiqua
halide talk里讲的方法我总结一下，请你看看大致是不是这样：

1. locality。实际上和linearize的思想差不多，就是将内存的读取做成cache-friendly，并且减少重复读取。具体的做法就是用当前读取到的内存做尽可能多的计算，这应该也就是fusion的意思。
2. parallel。这个和使用cpu向量指令集和openmp是一样的，但难点在于同时要照顾到locality。

halide会通过它的方法帮我们实现locality和parallel的优化，但代价就是要使用它的语法将计算组成pipeline。

但说实话，这个代价还挺大的，因为它的资料感觉挺少的，即使是opencv中常用kernel的组合，网上都找不到现成的例子。我自己尝试了halide的官方的例子，连编译都通不过。。。

另外，opencv似乎想要实现一个叫mini-halide的东西，但还只是在draft阶段。

[xinsuinizhuan]

@meiqua
halide talk里讲的方法我总结一下，请你看看大致是不是这样：

1. locality。实际上和linearize的思想差不多，就是将内存的读取做成cache-friendly，并且减少重复读取。具体的做法就是用当前读取到的内存做尽可能多的计算，这应该也就是fusion的意思。

2. parallel。这个和使用cpu向量指令集和openmp是一样的，但难点在于同时要照顾到locality。

halide会通过它的方法帮我们实现locality和parallel的优化，但代价就是要使用它的语法将计算组成pipeline。

但说实话，这个代价还挺大的，因为它的资料感觉挺少的，即使是opencv中常用kernel的组合，网上都找不到现成的例子。我自己尝试了halide的官方的例子，连编译都通不过。。。

另外，opencv似乎想要实现一个叫mini-halide的东西，但还只是在draft阶段。

可以参考我的思路走一走，进行加速。使用SSE优化soble gaussblur等函数进行优化，这个加速也是很明显的，只是不太懂SSE。参考博客：https://www.cnblogs.com/Imageshop/

[meiqua]

@zhengli233
1 应该说linearize就是locality的思想
确实，直接用halide会提高使用门槛，最好的方式是参考思想简单实现一个，因为这里也不用到halide这么复杂的程度。

[meiqua]

最近简单实现了整个fusion的过程，大家有兴趣可以跑跑看。

[DennisLiu1993]

@xinsuinizhuan
@zhengli233
兩位可以參考我的github，shape matching的替代方案，可以應用於許多場域
https://github.com/DennisLiu1993/Fastest_Image_Pattern_Matching

[xinsuinizhuan]

shape matching

好的。我有空测试一下您的方案。