logo标签列表

使用ARM NEON实现的快速高斯模糊图像滤波器

opencvimage-processinggaussianneonimagefilter
4
我正在尝试制作高斯模糊图像滤镜的移动快速版本。我已经阅读了其他问题，例如：Fast Gaussian blur on unsigned char image- ARM Neon Intrinsics- iOS Dev。为了我的目的，我只需要一个固定大小（7x7）和固定sigma（2）的高斯滤波器。因此，在优化ARM NEON之前，我正在使用C++实现1D高斯核，并将其与OpenCV GaussianBlur()方法在移动环境（带有NDK的Android）中直接进行性能比较。这样可以得到一个简单得多的代码来优化。然而，结果是我的实现比OpenCV4Android版本慢10倍。我已经阅读了OpenCV4 Tegra有优化的GaussianBlur实现，但我不认为标准的OpenCV4Android有这种优化，那么我的代码为什么这么慢呢？以下是我的实现（注意：在应用滤镜靠近边界时，使用reflect101进行像素反射）：
Mat myGaussianBlur(Mat src){
    Mat dst(src.rows, src.cols, CV_8UC1);
    Mat temp(src.rows, src.cols, CV_8UC1);
    float sum, x1, y1;

    // coefficients of 1D gaussian kernel with sigma = 2
    double coeffs[] = {0.06475879783, 0.1209853623, 0.1760326634, 0.1994711402, 0.1760326634, 0.1209853623, 0.06475879783};
    //Normalize coeffs
    float coeffs_sum = 0.9230247873f;
    for (int i = 0; i < 7; i++){
        coeffs[i] /= coeffs_sum;
    }

    // filter vertically
    for(int y = 0; y < src.rows; y++){
        for(int x = 0; x < src.cols; x++){
            sum = 0.0;
            for(int i = -3; i <= 3; i++){
                y1 = reflect101(src.rows, y - i);
                sum += coeffs[i + 3]*src.at<uchar>(y1, x);
            }
            temp.at<uchar>(y,x) = sum;
        }
    }

    // filter horizontally
    for(int y = 0; y < src.rows; y++){
        for(int x = 0; x < src.cols; x++){
            sum = 0.0;
            for(int i = -3; i <= 3; i++){
                x1 = reflect101(src.rows, x - i);
                sum += coeffs[i + 3]*temp.at<uchar>(y, x1);
            }
            dst.at<uchar>(y,x) = sum;
        }
    }

    return dst;
}
编辑：reflect101 可能有一点错误，它们的参数应该是 (src.cols，y+i) 和 (src.rows，x+i)。 - Alessandro Gaietta
4个回答
6
问题的一个重要部分在于算法过于精确，正如@PaulR所指出的那样。通常最好保持您的系数表不会比您的数据更精确。在这种情况下，因为您似乎正在处理uchar数据，所以您将使用大约8位系数表。
在NEON实现中，保持这些权重小尤其重要，因为算术越窄，您可以一次处理更多的通道。
除此之外，第一个主要的减速突出显示的是在主循环中具有图像边缘反射代码。这会使得大部分工作效率降低，因为通常在这种情况下不需要做任何特殊的事情。
如果您在边缘使用特殊版本的循环，然后当您安全时，使用简化的内部循环，不调用该reflect101()函数，可能会更好一些。
第二个问题(与原型代码更相关)是，在应用加权函数之前，可以将窗口的两侧翼相加，因为表格在两侧包含相同的系数。
sum = src.at<uchar>(y1, x) * coeffs[3];
for(int i = -3; i < 0; i++) {
    int tmp = src.at<uchar>(y + i, x) + src.at<uchar>(y - i, x);
    sum += coeffs[i + 3] * tmp;
}

这样每像素可以少进行6次乘法运算，也为控制溢出条件做了一步优化。

接下来有几个与内存系统相关的问题。

两次通过的方法在原则上是不错的，因为它可以避免大量的重新计算。但不幸的是，它可能会将有用的数据推出L1高速缓存，从而使所有操作变得更慢。它还意味着当您将结果写入内存时，您正在量化中间和，这可能会降低精度。

将此代码转换为NEON代码时，您需要关注的一点是尝试将工作集保持在寄存器文件中，但不要在完全利用之前丢弃计算。

当人们使用两次通过时，通常会将中间数据转置--即输入的一列变成输出的一行。

这是因为CPU不喜欢跨多行获取少量数据。如果将一堆水平像素收集在一起并进行过滤，则效率更高（因为缓存的工作方式）。如果临时缓冲区被转置，那么第二次通过也会将一堆水平点（在原始方向上为垂直）收集在一起，并再次转置其输出，以便正确输出。

如果您优化以使工作集局部化，则可能不需要此转置技巧，但最好知道它以便设置健康的基准性能。不幸的是，像这样本地化确实迫使您返回非最优内存获取，但使用更宽的数据类型可以缓解这种惩罚。

真的是非常好的建议！在应用加权函数之前将窗口的翼部相加以获得加速效果非常有趣。然而，考虑到转换到ARM NEON，我不知道是否值得这样做，因为如果可以在单个操作中完成所有乘法，则可能没有必要。 - Alessandro Gaietta
如果您一次完成，可以将图像的行折叠到彼此上，并过滤这些行以获取一个行的多个过滤像素，然后可以在寄存器文件中使用旋转窗口（大量的“VEXT”操作），并使用未折叠的方法在另一个方向上进行过滤。另一种方法可能是过滤图像的8x8或4x4块，使用相同的技术但在寄存器文件内部转置数据...但我不知道它会怎样，因为我从未尝试过。 - sh1
我的意思是我可以通过1个NEON操作一起完成7个乘法，因此在乘法之前添加翅膀是无用的。 - Alessandro Gaietta
2

如果这是针对8位图像的特定内容，那么您真的不希望使用浮点系数，尤其是双精度。此外，您也不希望对x1、y1使用浮点数。您应该只使用整数来表示坐标，并且可以使用定点数（即整数）来表示系数，以保持所有滤波算术在整数域中进行，例如：

Mat myGaussianBlur(Mat src){
    Mat dst(src.rows, src.cols, CV_8UC1);
    Mat temp(src.rows, src.cols, CV_16UC1); // <<<
    int sum, x1, y1;  // <<<

    // coefficients of 1D gaussian kernel with sigma = 2
    double coeffs[] = {0.06475879783, 0.1209853623, 0.1760326634, 0.1994711402, 0.1760326634, 0.1209853623, 0.06475879783};
    int coeffs_i[7] = { 0 }; // <<<
    //Normalize coeffs
    float coeffs_sum = 0.9230247873f;
    for (int i = 0; i < 7; i++){
        coeffs_i[i] = (int)(coeffs[i] / coeffs_sum * 256); // <<<
    }

    // filter vertically
    for(int y = 0; y < src.rows; y++){
        for(int x = 0; x < src.cols; x++){
            sum = 0; // <<<
            for(int i = -3; i <= 3; i++){
                y1 = reflect101(src.rows, y - i);
                sum += coeffs_i[i + 3]*src.at<uchar>(y1, x); // <<<
            }
            temp.at<uchar>(y,x) = sum;
        }
    }

    // filter horizontally
    for(int y = 0; y < src.rows; y++){
        for(int x = 0; x < src.cols; x++){
            sum = 0; // <<<
            for(int i = -3; i <= 3; i++){
                x1 = reflect101(src.rows, x - i);
                sum += coeffs_i[i + 3]*temp.at<uchar>(y, x1); // <<<
            }
            dst.at<uchar>(y,x) = sum / (256 * 256); // <<<
        }
    }

    return dst;
}
这个近似并不是很好。我的算法与OpenCV GaussianBlur结果像素的平均差为1.83，而这个近似的平均差大于5。然而，这个想法显然是好的，可以将计算时间减少25-30％，然后用65536替换256也可以保持相当高的精度。 - Alessandro Gaietta
是的，您可以在固定点系数的准确性与头房、计算成本等之间进行权衡。您还可以使用舍入而不是截断来提高准确性。更进一步的改进是将临时数据保持在比输入/输出像素更高的分辨率，并在最后处理缩放。 - Paul R
请注意，我现在已经修改了上面的代码，将temp保持在16位，并在最后组合X/Y缩放因子-这应该显着减少总误差。 - Paul R
2
这是在实施@Paul R和@sh1的所有建议后的代码，总结如下：
1）仅使用整数算术（精度自定）。
2）在应用乘法之前，将距离掩模中心相同距离的像素值相加，以减少乘法的数量。
3）仅应用水平滤波器，以利用矩阵按行存储的优势。
4）将边缘周围的循环与图像内部的循环分开，以避免不必要地调用反射函数。我完全删除了反射功能，并将其包含在沿边缘的循环中。
5）此外，作为个人观察，为了改善舍入而不调用“round”或“cvRound”（慢速）函数，我向临时和最终像素结果添加了0.5f（= 32768在整数精度中），以减少误差/与OpenCV的差异。
现在性能比OpenCV慢大约15倍到6倍。
然而，生成的矩阵与OpenCV的高斯模糊得到的矩阵并不完全相同。这不是由于算术长度（足够长）以及误差的去除。请注意，这是两个版本产生的矩阵之间的最小差异，在绝对值上介于0和2之间。系数与OpenCV使用的相同，使用相同大小和sigma获得的getGaussianKernel。
Mat myGaussianBlur(Mat src){

Mat dst(src.rows, src.cols, CV_8UC1);
Mat temp(src.rows, src.cols, CV_8UC1);
int sum;
int x1;

double coeffs[] = {0.070159, 0.131075, 0.190713, 0.216106, 0.190713, 0.131075, 0.070159};
int coeffs_i[7] = { 0 };
for (int i = 0; i < 7; i++){
        coeffs_i[i] = (int)(coeffs[i] * 65536); //65536
}

// filter horizontally - inside the image
for(int y = 0; y < src.rows; y++){
    uchar *ptr = src.ptr<uchar>(y);
    for(int x = 3; x < (src.cols - 3); x++){
        sum = ptr[x] * coeffs_i[3];
        for(int i = -3; i < 0; i++){
            int tmp = ptr[x+i] + ptr[x-i];
            sum += coeffs_i[i + 3]*tmp;
        }
        temp.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}
// filter horizontally - edges - needs reflect
for(int y = 0; y < src.rows; y++){
    uchar *ptr = src.ptr<uchar>(y);
    for(int x = 0; x <= 2; x++){
        sum = 0;
        for(int i = -3; i <= 3; i++){
            x1 = x + i;
            if(x1 < 0){
                x1 = -x1;
            }
            sum += coeffs_i[i + 3]*ptr[x1];
        }
        temp.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}
for(int y = 0; y < src.rows; y++){
    uchar *ptr = src.ptr<uchar>(y);
    for(int x = (src.cols - 3); x < src.cols; x++){
        sum = 0;
        for(int i = -3; i <= 3; i++){
            x1 = x + i;
            if(x1 >= src.cols){
                x1 = 2*src.cols - x1 - 2;
            }
            sum += coeffs_i[i + 3]*ptr[x1];
        }
        temp.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}

// transpose to apply again horizontal filter - better cache data locality
transpose(temp, temp);

// filter horizontally - inside the image
for(int y = 0; y < src.rows; y++){
    uchar *ptr = temp.ptr<uchar>(y);
    for(int x = 3; x < (src.cols - 3); x++){
        sum = ptr[x] * coeffs_i[3];
        for(int i = -3; i < 0; i++){
            int tmp = ptr[x+i] + ptr[x-i];
            sum += coeffs_i[i + 3]*tmp;
        }
        dst.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}
// filter horizontally - edges - needs reflect
for(int y = 0; y < src.rows; y++){
    uchar *ptr = temp.ptr<uchar>(y);
    for(int x = 0; x <= 2; x++){
        sum = 0;
        for(int i = -3; i <= 3; i++){
            x1 = x + i;
            if(x1 < 0){
                x1 = -x1;
            }
            sum += coeffs_i[i + 3]*ptr[x1];
        }
        dst.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}
for(int y = 0; y < src.rows; y++){
    uchar *ptr = temp.ptr<uchar>(y);
    for(int x = (src.cols - 3); x < src.cols; x++){
        sum = 0;
        for(int i = -3; i <= 3; i++){
            x1 = x + i;
            if(x1 >= src.cols){
                x1 = 2*src.cols - x1 - 2;
            }
            sum += coeffs_i[i + 3]*ptr[x1];
        }
        dst.at<uchar>(y,x) = (sum + 32768) / 65536;
    }
}

transpose(dst, dst);

return dst;
}
0
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接