具有 FFT 互相关的 C++ 模式匹配（图片）

Question

每个人我都在尝试用 FFT 实现模式匹配，但我不确定结果应该是什么（我想我遗漏了一些东西，尽管阅读了很多关于这个问题的东西并尝试了很多不同的实现这个是迄今为止最好的）。这是我的 FFT 相关函数。

void fft2d(fftw_complex**& a, int rows, int cols, bool forward = true)
{
    fftw_plan p;
    for (int i = 0; i < rows; ++i)
    {
        p = fftw_plan_dft_1d(cols, a[i], a[i], forward ? FFTW_FORWARD : FFTW_BACKWARD, FFTW_ESTIMATE);
        fftw_execute(p);
    }

    fftw_complex* t = (fftw_complex*)fftw_malloc(rows * sizeof(fftw_complex));

    for (int j = 0; j < cols; ++j)
    {
        for (int i = 0; i < rows; ++i)
        {
            t[i][0] = a[i][j][0];
            t[i][1] = a[i][j][1];
        }

        p = fftw_plan_dft_1d(rows, t, t, forward ? FFTW_FORWARD : FFTW_BACKWARD, FFTW_ESTIMATE);
        fftw_execute(p);

        for (int i = 0; i < rows; ++i)
        {
            a[i][j][0] = t[i][0];
            a[i][j][1] = t[i][1];
        }
    }

    fftw_free(t);
}

int findCorrelation(int argc, char* argv[])
{
    BMP bigImage;
    BMP keyImage;
    BMP result;
    RGBApixel blackPixel = { 0, 0, 0, 1 };
    const bool swapQuadrants = (argc == 4);

    if (argc < 3 || argc > 4) {
        cout << "correlation img1.bmp img2.bmp" << endl;
        return 1;
    }

    if (!keyImage.ReadFromFile(argv[1])) {
        return 1;
    }

    if (!bigImage.ReadFromFile(argv[2])) {
        return 1;
}
    //Preparations
    const int maxWidth = std::max(bigImage.TellWidth(), keyImage.TellWidth());
    const int maxHeight = std::max(bigImage.TellHeight(), keyImage.TellHeight());

    const int rowsCount = maxHeight;
    const int colsCount = maxWidth;

    BMP bigTemp = bigImage;
    BMP keyTemp = keyImage;

    keyImage.SetSize(maxWidth, maxHeight);
    bigImage.SetSize(maxWidth, maxHeight);

    for (int i = 0; i < rowsCount; ++i)
        for (int j = 0; j < colsCount; ++j) {
            RGBApixel p1;

            if (i < bigTemp.TellHeight() && j < bigTemp.TellWidth()) {
                p1 = bigTemp.GetPixel(j, i);
            } else {
                p1 = blackPixel;
            }

            bigImage.SetPixel(j, i, p1);

            RGBApixel p2;

            if (i < keyTemp.TellHeight() && j < keyTemp.TellWidth()) {
                p2 = keyTemp.GetPixel(j, i);
            } else {
                p2 = blackPixel;
            }

            keyImage.SetPixel(j, i, p2);
        }
    //Here is where the transforms begin
    fftw_complex **a = (fftw_complex**)fftw_malloc(rowsCount * sizeof(fftw_complex*));
    fftw_complex **b = (fftw_complex**)fftw_malloc(rowsCount * sizeof(fftw_complex*));
    fftw_complex **c = (fftw_complex**)fftw_malloc(rowsCount * sizeof(fftw_complex*));

    for (int i = 0; i < rowsCount; ++i) {

        a[i] = (fftw_complex*)fftw_malloc(colsCount * sizeof(fftw_complex));
        b[i] = (fftw_complex*)fftw_malloc(colsCount * sizeof(fftw_complex));
        c[i] = (fftw_complex*)fftw_malloc(colsCount * sizeof(fftw_complex));

        for (int j = 0; j < colsCount; ++j) {
            RGBApixel p1;

            p1 = bigImage.GetPixel(j, i);

            a[i][j][0] = (0.299*p1.Red + 0.587*p1.Green + 0.114*p1.Blue);
            a[i][j][1] = 0.0;

            RGBApixel p2;

            p2 = keyImage.GetPixel(j, i);

            b[i][j][0] = (0.299*p2.Red + 0.587*p2.Green + 0.114*p2.Blue);
            b[i][j][1] = 0.0;
        }
    }

    fft2d(a, rowsCount, colsCount);
    fft2d(b, rowsCount, colsCount);

    result.SetSize(maxWidth, maxHeight);

    for (int i = 0; i < rowsCount; ++i)
        for (int j = 0; j < colsCount; ++j) {
            fftw_complex& y = a[i][j];
            fftw_complex& x = b[i][j];

            double u = x[0], v = x[1];
            double m = y[0], n = y[1];

            c[i][j][0] = u*m + n*v;
            c[i][j][1] = v*m - u*n;

            int fx = j;
            if (fx>(colsCount / 2)) fx -= colsCount;

            int fy = i;
            if (fy>(rowsCount / 2)) fy -= rowsCount;

            float r2 = (fx*fx + fy*fy);

            const double cuttoffCoef = (maxWidth * maxHeight) / 37992.;
            if (r2<128 * 128 * cuttoffCoef)
                c[i][j][0] = c[i][j][1] = 0;
        }

    fft2d(c, rowsCount, colsCount, false);

    const int halfCols = colsCount / 2;
    const int halfRows = rowsCount / 2;

    if (swapQuadrants) {
        for (int i = 0; i < halfRows; ++i)
            for (int j = 0; j < halfCols; ++j) {
                std::swap(c[i][j][0], c[i + halfRows][j + halfCols][0]);
                std::swap(c[i][j][1], c[i + halfRows][j + halfCols][1]);
            }

        for (int i = halfRows; i < rowsCount; ++i)
            for (int j = 0; j < halfCols; ++j) {
                std::swap(c[i][j][0], c[i - halfRows][j + halfCols][0]);
                std::swap(c[i][j][1], c[i - halfRows][j + halfCols][1]);
            }
    }

    for (int i = 0; i < rowsCount; ++i)
        for (int j = 0; j < colsCount; ++j) {
            const double& g = c[i][j][0];
            RGBApixel pixel;

            pixel.Alpha = 0;
            int gInt = 255 - static_cast<int>(std::floor(g + 0.5));
            pixel.Red = gInt;
            pixel.Green = gInt;
            pixel.Blue = gInt;

            result.SetPixel(j, i, pixel);
        }


    BMP res;
    res.SetSize(maxWidth, maxHeight);

    result.WriteToFile("result.bmp");

    return 0;
}

示例输出

Answer 1

这个问题可能更适合发布在另一个网站上，例如 cross validated（metaoptimize.com 曾经也是一个很好的问题，但它似乎已经消失了）

也就是说：

您可以使用 FFT 执行两种类似的操作：卷积和相关。卷积用于确定两个信号如何相互作用，而相关性可用于表示两个信号彼此的相似程度。确保您正在执行正确的操作，因为它们通常都是通过 DFT 实现的。

对于这种类型的 DFT 应用，您通常不会在傅里叶频谱中提取任何有用的信息，除非您正在寻找两个数据源或其他任何东西共有的频率（例如，如果您正在比较两个桥以查看是否他们的支持间隔相似）。

你的第三张图片看起来很像电源域；通常我看到相关输出完全是灰色的，除了发生重叠的地方。您的代码肯定看起来是在计算逆 DFT，所以除非我遗漏了一些东西，否则我对模糊外观提出的唯一其他解释可能是其中的一些 "fudge factor" 代码，例如：

if (r2<128 * 128 * cuttoffCoef)
  c[i][j][0] = c[i][j][1] = 0;

至于您应该期待什么：只要两幅图像之间有共同元素，您就会看到一个峰值。峰值越大，表示该区域附近的两幅图像越相似。

一些评论and/or 建议修改：

1) 卷积和相关不是尺度不变操作。换句话说，图案图像的大小会对输出产生重大影响。

2) 在相关之前对图像进行标准化。

当您为前向 DFT 传递准备好图像数据时：

a[i][j][0] = (0.299*p1.Red + 0.587*p1.Green + 0.114*p1.Blue);
a[i][j][1] = 0.0;
/* ... */

如何对图像进行灰度化是您的事（尽管我会选择 sqrt( r*r + b*b + g*g ) 之类的东西）。但是，我没有看到你做任何事情来规范化图像。

"normalize" 这个词在这种情况下可以有几种不同的含义。两种常见类型：

标准化 0.0 和 1.0 之间的值范围
标准化图像的 "whiteness"

3) 运行通过边缘增强滤镜处理您的图案图像。我个人使用了在线免费提供的 "dsp guide" 书中的 canny, sobel, and I think I messed with a few others. As I recall, canny was "quick'n dirty", sobel was more expensive, but I got comparable results when it came time to do correlation. See chapter 24。整本书都值得你花时间，但如果你时间不够，那么至少第 24 章会有很大帮助。

4) 在[0, 255]之间重新缩放输出图像；如果你想实现阈值，在这一步之后做，因为阈值步骤是有损的。

我对这个的记忆很模糊，但我记得（为清楚起见进行了编辑）：

您可以通过从整个功率谱中除掉最大的功率谱值来在 [-1.0, 1.0] 之间缩放最终图像像素（在重新缩放之前）
最大的功率谱值，足够方便，是功率谱中最中心的值（对应于最低频率）
如果你把它从功率谱中除掉，你最终会做两倍的功；由于 FFT 是线性的，因此您可以将除法延迟到逆 DFT 传递到您重新缩放 [0..255].
如果在重新调整后你的大部分值变得很黑你看不到它们，你可以使用 ODE y' = y(1 - y) 的解决方案（一个例子是 sigmoid f(x) = 1 / (1 + exp(-c*x) )，对于一些比例因子 c，可提供更好的层次）。这更多地是为了提高您以视觉方式解释结果的能力，而不是您可能用来以编程方式查找峰值的任何东西。

edit 我上面说的是[0, 255]。我建议您重新调整为 [128, 255] 或其他一些灰色而不是黑色的下限。

具有 FFT 互相关的 C++ 模式匹配（图片）

C++ Pattern Matching with FFT cross-correlation (Images)

c++

image-recognition

pattern-matching

fftw

cross-correlation