我该如何优化这个侵蚀过滤器代码？

Question

我已经实现了将侵蚀滤镜应用于图像的函数

void applyErosionFilter(QImage &input, int matrixSize)
{
    int filterOffset = (matrixSize - 1) / 2;
    int byteOffset = 0;

    uchar red, green, blue = 0;
    uchar morphResetValue = 255;

    uchar *data = input.bits();
    int stride = input.bytesPerLine();
    uchar *newdata = new uchar[stride * input.height()];

    int i = 0;
    for (int y = filterOffset; y < input.height() - filterOffset; y++)
    {
        for (int x = filterOffset; x < input.width() - filterOffset; x++)
        {
            byteOffset = y * stride + x * 4;

            red = morphResetValue;
            green = morphResetValue;
            blue = morphResetValue;

            for (int filterY = -filterOffset; filterY <= filterOffset; filterY++)
            {
                for (int filterX = -filterOffset; filterX <= filterOffset; filterX++)
                {
                    i = byteOffset + (filterX * 4) + (filterY * stride);

                    if (data[i] < red)
                        red = data[i];

                    if (data[i + 1] < green)
                        green = data[i + 1];

                    if (data[i + 2] < blue)
                        blue = data[i + 2];
                }
            }

            newdata[byteOffset] = red;
            newdata[byteOffset + 1] = green;
            newdata[byteOffset + 2] = blue;
            newdata[byteOffset + 3] = 255;
        }
    }

    input = input.fromImage(QImage(newdata, input.width(), input.height(), QImage::Format::Format_ARGB32));
    delete [] newdata;
}

它工作得很好，但我一直想知道是否有另一种方法可以更有效地做到这一点，也许有一种方法可以使用 openGL 在 GPU 上执行计算。

Answer 1

作为一种快速改进，我建议使用多个线程并行计算多行。您也可以使用 OpenCL 或 Cuda 在 GPU 上实现此目的，但这需要大量样板文件。

我已修改您的代码以使用多线程，但我尚未对此进行测试，因为我目前没有在此设备上安装 Qt。但这至少可以提示您从哪里开始。

（顺便说一句，uchar red, green, blue = 0; 仅将 blue 初始化为 0，而 red 和 green 保持未初始化）

#include <thread>
#include <vector>
void applyErosionFilter(QImage &input, int matrixSize)
{
    int filterOffset = (matrixSize - 1) / 2;
    int byteOffset = 0;

    uchar morphResetValue = 255;

    uchar *data = input.bits();
    int stride = input.bytesPerLine();
    uchar *newdata = new uchar[stride * input.height()];

    unsigned num_threads = std::thread::hardware_concurrency();
    if (num_threads == 0)
        num_threads = 1;
    std::vector<std::thread> threads;
    int i = 0;
    for (unsigned i = 0; i < num_threads; ++i)
    {
        int start = (input.height() - 2 * filterOffset) * i / num_threads;
        int end = (input.height() - 2 * filterOffset) * (i+1) / num_threads;
        threads.emplace_back([&](int start_line, int end_line){
            unsigned char red = 0;
            unsigned char green = 0;
            unsigned char blue = 0;
            for (int y = start_line; y < end_line; y++)
            {
                for (int x = filterOffset; x < input.width() - filterOffset; x++)
                {
                    byteOffset = y * stride + x * 4;

                    red = morphResetValue;
                    green = morphResetValue;
                    blue = morphResetValue;

                    for (int filterY = -filterOffset; filterY <= filterOffset; filterY++)
                    {
                        for (int filterX = -filterOffset; filterX <= filterOffset; filterX++)
                        {
                            i = byteOffset + (filterX * 4) + (filterY * stride);

                            if (data[i] < red)
                                red = data[i];

                            if (data[i + 1] < green)
                                green = data[i + 1];

                            if (data[i + 2] < blue)
                                blue = data[i + 2];
                        }
                    }

                    newdata[byteOffset] = red;
                    newdata[byteOffset + 1] = green;
                    newdata[byteOffset + 2] = blue;
                    newdata[byteOffset + 3] = 255;
                }
            }
        },
            start, end);
    }
    for (auto &thread : threads)
        thread.join();

    input = input.fromImage(QImage(newdata, input.width(), input.height(), QImage::Format::Format_ARGB32));
    delete [] newdata;
}

Answer 2

1. 使用多线程：将图像水平切割成带状，并分别处理。我是用 OpenMP 做的，它真的很困难。
2. 大小为 NxN 的平方结构元素可以分解为大小为 1xN 和 Nx1 的两段（水平和垂直）。因此，您无需对每个像素进行 NxN 测试，而是进行 2xN：N=3 9 vs 6、n=5 25 vs 10 等等。速度要快得多！
3. 使用已在以下库中实现的算法：SMIL，Matthieu Faessel（基于行比较的 C++ 自动矢量化代码，最快！！！），libmorpho，Marc Van Droogenbroeck（C++ 但有限到 8 位编码，如果我是对的），或者 Christophe Clienti 的 Fulguro（非常适合 SIMD 优化的朗伯算法）。正如您在这些库中看到的那样，它们使用智能方法/算法/架构来快速获得结果。您开发的是课程中讲授的基础知识，以便于理解，但它是最慢的！

Answer 3

谢谢大家。 我找到了我想要的东西，我想和你分享。

QOpenGLTexture *m_texImageInput;
QOpenGLShaderProgram *m_shaderComputeH;

void initiateShader()
{
    if (m_texImageInput)
    {
        delete m_texImageInput;
        m_texImageInput = nullptr;
    }

    QImage img(":/image.png");
    m_texImageInput = new QOpenGLTexture(img.convertToFormat(QImage::Format_RGBA8888).mirrored());

    if (m_shaderComputeH)
    {
        delete m_shaderComputeH;
        m_shaderComputeH = nullptr;
    }

    m_shaderComputeH = new QOpenGLShaderProgram;
    m_shaderComputeH->addShaderFromSourceFile(QOpenGLShader::Compute, ":/csErosionFilter.fsh");
    m_shaderComputeH->link();
}

QSize getWorkGroups(int workGroupSize, const QSize &imageSize)
{
    int x = imageSize.width();
    x = (x % workGroupSize) ? (x / workGroupSize) + 1 : (x / workGroupSize);
    int y = imageSize.height();
    y = (y % workGroupSize) ? (y / workGroupSize) + 1 : (y / workGroupSize);

    return QSize(x, y);
}

void executeFilter(int radius)
{
    QOpenGLExtraFunctions *f = QOpenGLContext::currentContext()->extraFunctions();

    // Process input image
    QSize workGroups = getWorkGroups(32, QSize(m_texImageInput->width(), m_texImageInput->height()));
    // Pass 2
    f->glBindImageTexture(0, m_texImageInput->textureId(), 0, 0, 0,  GL_READ_WRITE, GL_RGBA8);
    f->glBindImageTexture(1, m_texImageProcessed->textureId(), 0, 0, 0,  GL_READ_WRITE, GL_RGBA8);
    m_shaderComputeH->bind();
    m_shaderComputeH->setUniformValue("radius", radius);
    f->glDispatchCompute(workGroups.width(), workGroups.height(), 1);
    f->glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);
    m_shaderComputeH->release();
    // Compute cleanup
    f->glBindImageTexture(0, 0, 0, 0, 0, GL_READ_WRITE, GL_RGBA8);
    f->glBindImageTexture(1, 0, 0, 0, 0, GL_READ_WRITE, GL_RGBA8);
}

这是计算着色器

   /// csErosionFilter.fsh
    #version 430 core
    #define COMPUTEPATCHSIZE 32
    layout (local_size_x = COMPUTEPATCHSIZE, local_size_y = COMPUTEPATCHSIZE) in;
    layout(binding=0, rgba8) uniform readonly highp image2D inputImage;
    layout(binding=1, rgba8) uniform writeonly highp image2D resultImage;
    uniform int radius;
    void main()
    {
        ivec2 imgSize = imageSize(resultImage);
        int x = int(gl_GlobalInvocationID.x);
        int y = int(gl_GlobalInvocationID.y);
        if ((x >= imgSize.x) || (y >= imgSize.y))
            return;

        vec4 newValue = vec4(1);
        int left   = clamp(x - radius, 0, imgSize.x - 1);
        int right  = clamp(x + radius, 0, imgSize.x - 1);
        int top    = clamp(y - radius, 0, imgSize.y - 1);
        int bottom = clamp(y + radius, 0, imgSize.y - 1);
        for (int iX = left; iX <= right; iX++)
        {
            for (int iY = top; iY <= bottom; iY++)
            {
                vec4 value = imageLoad(inputImage, ivec2(iX, iY));
                if(value.x < newValue.x)
                    newValue.x = value.x;
                if(value.y < newValue.y)
                    newValue.y = value.y;
                if(value.z < newValue.z)
                    newValue.z = value.z;
            }
        }

        imageStore(resultImage, ivec2(x,y), newValue);
    }