使用 ARM NEON 的快速高斯模糊图像过滤器

Question

我正在尝试制作高斯模糊图像滤镜的移动快速版本。

我读过其他问题，例如：Fast Gaussian blur on unsigned char image- ARM Neon Intrinsics- iOS Dev

出于我的目的，我只需要一个固定大小 (7x7) 固定 sigma (2) 高斯滤波器。

因此，在针对 ARM NEON 进行优化之前，我正在 C++ 中实现一维高斯内核，并直接在移动环境（带有 NDK 的 Android）中将性能与 OpenCV GaussianBlur() 方法进行比较。这样一来，优化的代码就会简单得多。

但是结果是我的实现比 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;
}

Answer 1

根据 Google 文档，在 Android 设备上，使用 float/double 比使用 int/uchar 慢两倍。

您可能会在此 Android 文档中找到一些加快 C++ 代码速度的解决方案。 https://developer.android.com/training/articles/perf-tips

Answer 2

这是实现@Paul R 和@sh1 的所有建议后的代码，总结如下：

1) 只使用整数算术（精确到口味）

2）在进行乘法运算之前，将与mask中心相同距离的像素值相加，以减少乘法次数。

3) 仅应用水平过滤器以利用矩阵行的存储

4) 将边缘周围的循环与图像内部的循环分开，以免对反射函数进行不必要的调用。我完全移除了反射的功能，包括它们在沿边缘的循环内。

5) 此外，作为个人观察，为了在不调用（慢）函数“round”或“cvRound”的情况下改进舍入，我在临时和最终像素结果中添加了 0.5f（= 32768 整数精度) 以减少与 OpenCV 相比的错误/差异。

现在性能比 OpenCV 慢了大约 15 到大约 6 倍。

但是，生成的矩阵与使用 OpenCV 的高斯模糊获得的矩阵并不完全相同。这不是由于算术长度（足够）以及消除错误仍然存​​在。请注意，这是两个版本产生的矩阵之间像素强度的最小差异，介于 0 和 2（绝对值）之间。系数与 OpenCV 使用的相同，使用 getGaussianKernel 获得，大小和 sigma 相同。

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;
}

Answer 3

如果这是专门针对 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;
}

Answer 4

正如@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 真的不喜欢跨输入图像的多行获取少量数据。如果您收集一堆水平像素并过滤它们，它的工作效率会更高（因为缓存的工作方式）。如果临时缓冲区被转置，那么第二遍也收集一堆水平点（在原始方向上是垂直的）并再次转置其输出，因此它以正确的方式出现。

如果您优化以使您的工作集保持本地化，那么您可能不需要这种转置技巧，但值得了解一下，这样您就可以为自己设置一个健康的基准性能。不幸的是，像这样的本地化确实会迫使您返回到非最佳内存提取，但使用更广泛的数据类型可以减轻惩罚。