傅里叶空间应用高斯核时出现图像阴影/旋转/翻转问题排查

前置声明：我知道可以用OpenCV实现图像平滑/模糊，但本次作业要求底层实现以掌握完整流程。我刚接触图像处理和DFT相关知识，目前遇到了相关难题。

问题描述 #

我通过FFTW库和OpenCV进行图像处理，将图像读取为OpenCV::Mat对象并转换为grayscale灰度空间。为满足FFTW的FFT函数输入要求，我将数据转换为double类型，存储在名为imageDouble的Mat对象中，对应指针为pImageDouble。

std::string filename = "<path_to_your_image>";
Mat image = imread(filename, IMREAD_GRAYSCALE);

/***********************************************************************
*** creating Mat object of type *doube* and copying image contents to it
*** then creating a pointer to it.
***********************************************************************/
Mat imageDouble;
image.convertTo(imageDouble, CV_64FC1);
double* pImageDouble = imageDouble.ptr<double>(0);

我同时将高斯核读入Mat对象，先执行循环移位操作，将高斯核的中心（数值最高点）移位到核的左上角(0,0)位置，随后对核进行零填充至与输入原图尺寸一致，结果（double类型）存储在名为paddedGKernel的Mat对象中，对应指针为pPaddedGKernel。

Mat paddedGKernel = padGKernelWithZeros(nRows, nCols, rolledGKernel);
double* pPaddedGKernel = paddedGKernel.ptr<double>(0);

我初始化了fftw_complex对象存储FFT输出结果、fftw_plan对象，完成内存分配后执行变换。我调用FFTW的fftw_plan_dft_r2c_2d()函数分别对输入图像、移位并填充后的高斯核执行二维FFT，随后在傅里叶空间做点乘操作，以对原图施加高斯滤波。

fftw_complex* out;               //for result of FFT for original image
fftw_complex* outGaussian;       //for result of FFT for Gaussian kernel
fftw_plan p;

/*************************************************
*** Allocating memory for the fftw_complex objects
*************************************************/
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * nRows * nCols);
outGaussian = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * nRows * nCols);


/****************************************
*** FFT on imageDouble, outputting to out
****************************************/
p = fftw_plan_dft_r2c_2d(imageDouble.rows, imageDouble.cols, pImageDouble, out, FFTW_ESTIMATE);
fftw_execute(p);


/**************************************************
*** FFT on paddedGKernel, outputting to outGuassian
**************************************************/
p = fftw_plan_dft_r2c_2d(paddedGKernel.rows, paddedGKernel.cols, pPaddedGKernel, outGaussian, FFTW_ESTIMATE);
fftw_execute(p);


/****************************************************
*** Pointwise multiplication to apply Gaussian kernel
****************************************************/
for (size_t i = 0; i < nCCols * nRows; i++)
{
    out[i][0] = out[i][0] * (1 / outGaussian[i][0]);
}

之后我调用fftw_plan_dft_c2r_2d()执行逆FFT，将点乘结果输出到OpenCV::Mat对象imageCloneDouble中，再将数据转换为uchar类型存储到OpenCV::Mat imageBack。

/***********************************************************************
*** Making Mat object (type *double*) to put data into and pointer to it
***********************************************************************/
Mat imageCloneDouble = Mat::zeros(nRows, nCols, CV_64FC1);
double* pImageCloneDouble = imageCloneDouble.ptr<double>(0);

/*****************************************************************************
*** Making and executing plan for inverse FFT, outputting to imageCloneDouble,
*** then normalizing since this function puts out unnormalized values.
*****************************************************************************/
fftw_plan pp = fftw_plan_dft_c2r_2d(imageCloneDouble.rows, imageCloneDouble.cols, out, pImageCloneDouble, FFTW_BACKWARD | FFTW_ESTIMATE);
fftw_execute(pp);
imageCloneDouble = imageCloneDouble / (nCols * nRows *2);

/***************************************************************************
*** Making Mat object (type *uchar*) and copying data inverse FFT data to it
***************************************************************************/
Mat imageBack;
imageCloneDouble.convertTo(imageBack, CV_8UC1);

我期望输出imageBack为施加高斯模糊后的原图，但实际输出的图像上叠加了疑似旋转、修改后的原图副本，还存在阴影问题。我无法确定翻转/旋转及叠加效果的产生原因，推测问题出在傅里叶空间点乘环节，或是遗漏了点乘结果的必要处理步骤。请问我需要对傅里叶空间的数据做什么处理，才能得到正确的滤波结果？

你代码中存在5个核心错误，逐一修正即可得到正常的高斯模糊结果：

• 1. 频域运算逻辑错误
     你当前写的实部相除操作是反卷积逻辑，而非高斯模糊需要的复数逐点相乘。傅里叶域的卷积对应两个复数完整相乘，不能只操作实部，修正后的循环代码如下：

// 注意r2c输出的复数个数是nRows*(nCols/2 +1)
for (size_t i = 0; i < nRows * (nCols / 2 + 1); i++)
{
    double real = out[i][0] * outGaussian[i][0] - out[i][1] * outGaussian[i][1];
    double imag = out[i][0] * outGaussian[i][1] + out[i][1] * outGaussian[i][0];
    out[i][0] = real;
    out[i][1] = imag;
}

• 2. FFTW内存分配与循环边界错误
     fftw_plan_dft_r2c_2d实数转复数的2D变换利用了共轭对称性，输出的复数数组长度为nRows * (nCols / 2 + 1)而非nRows * nCols，你之前的内存分配不仅冗余，还可能出现边界访问越界，对应修改内存分配代码：

out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * nRows * (nCols / 2 + 1));
outGaussian = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * nRows * (nCols / 2 + 1));

• 3. 逆变换归一化系数错误
     FFTW的c2r逆变换输出结果仅需要除以总像素数nRows * nCols即可完成归一化，你额外除以2会导致图像整体亮度偏低、灰度截断出现异常阴影，修正为：

imageCloneDouble = imageCloneDouble / (double)(nCols * nRows);

• 4. 逆变换后缺失循环移位操作
     你预处理高斯核时将核中心移到了左上角，逆变换输出的空域结果默认是经过循环移位的状态，需要再次执行ifftshift（将左上角的有效内容移回图像中心）才能得到正常的滤波结果，在转uchar之前添加象限交换代码：

int cx = imageCloneDouble.cols / 2;
int cy = imageCloneDouble.rows / 2;
Mat tmp;
// 拆分四个象限
Mat q0(imageCloneDouble, Rect(0, 0, cx, cy));
Mat q1(imageCloneDouble, Rect(cx, 0, imageCloneDouble.cols - cx, cy));
Mat q2(imageCloneDouble, Rect(0, cy, cx, imageCloneDouble.rows - cy));
Mat q3(imageCloneDouble, Rect(cx, cy, imageCloneDouble.cols - cx, imageCloneDouble.rows - cy));
// 交换象限完成ifftshift
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);

• 5. 高斯核未归一化
     你填充后的paddedGKernel所有像素值之和需要为1，否则会导致输出图像整体亮度异常，若未做归一化可以在移位填充后添加代码：

paddedGKernel /= sum(paddedGKernel)[0];

内容的提问来源于stack exchange，提问作者T-Bo