C++ - 本征 FFT 用于二维图像的反卷积
C++ - Eigen FFT use for Deconvolution of 2 dimensional image
我正在尝试对图像 I 执行反卷积,即 n x m。用于对其进行卷积的内核是 K,它也是 n x m。
现在我想通过执行反卷积找到原始图像 O。我知道我可以通过对 I 和 K 执行傅里叶变换并除法来检索图像 O :I / K,因为在傅立叶域中卷积是一个乘积。 (我从 here 那里得到了这些信息)。
我看到另一个关于如何使用 Eigen FFT 执行前向变换的 post,here。
我的正向傅立叶变换代码是:
I = 输入图像(时域)
O = 输出图像(频域)
tempFreq = 用于计算的临时矩阵(频域)
timevec1 = 浮点矢量
freqvec1、freqvec2、freqvec3 = 复向量
for (Int32 i = 0; i < I->uRowsCount; ++i)
{
for (Int32 j = 0; j < I->uColumnsCount; ++j)
{
timevec1.push_back((*I)(i, j));
}
fft.fwd(freqvec1, timevec1);
for (Int32 j = 0; j < I->uColumnsCount; ++j)
{
(tempFreq)(i, j) = freqvec1[j];
}
freqvec1.clear();
timevec1.clear();
}
freqvec1.clear();
timevec1.clear();
for (Int32 j = 0; j < I->uColumnsCount; ++j)
{
for (Int32 i = 0; i < I->uRowsCount; ++i)
{
freqvec2.push_back((tempFreq)(i, j));
}
fft.fwd(freqvec1, freqvec2);
for (Int32 i = 0; i < I->uRowsCount; ++i)
{
(O)(i, j) = freqvec1[i];
}
freqvec2.clear();
freqvec1.clear();
}
我的傅里叶逆变换代码是:
I = 输入图像(频域)
O = 输出图像(时域)
tempTime = 用于计算的临时矩阵(时域)
timevec1, timevec2 = 浮点矢量
freqvec1, freqvec2 = 复向量
for (Int32 j = 0; j < O->uColumnsCount; ++j)
{
for (Int32 i = 0; i < O->uRowsCount; ++i)
{
freqvec1.push_back((I)(i, j));
}
fft.inv(timevec1, freqvec1);
for (Int32 i = 0; i < O->uRowsCount; ++i)
{
(*tempTime)(i, j) = timevecCol[i];
}
freqvec1.clear();
timevec1.clear();
}
freqvec1.clear();
timevec1.clear();
for (Int32 i = 0; i < O->uRowsCount; ++i)
{
for (Int32 j = 0; j < O->uColumnsCount; ++j)
{
freqvec2.push_back((*tempTime)(i, j));
}
fft.inv(timevec2, freqvec2);
for (Int32 j = 0; j < O->uColumnsCount; ++j)
{
(*O)(i, j) = timevec2[j];
}
freqvec2.clear();
timevec2.clear();
}
对于反卷积,我将频域中的输入图像与频域中的内核分开:
freqDomainOutputImage = freqDomainInputImage.cwiseQuotient(freqDomainKernel);
为了获得原始图像,我对 freqDomainOutputImage 执行了逆傅立叶变换。
我得到的结果是:
我相信 FFT 将左上角镜像到另一边,但我不知道为什么?我没有使用 halfSpectrum。第二,为什么图像偏移了?如果我通过用这个替换最后一个循环将输出图像移动到中心:
(*O)((i + O->uRowsCount/2)%O->uRowsCount, (j + O->uColumnsCount/2)%O->uColumnsCount) = timevec2[j];
那么我的输出是:
(可以看到图像是从左上象限镜像的)。
最后,为什么我自己添加了没有噪点的模糊,但它似乎有噪点?
tempTime
需要很复杂。你似乎扔掉了那个虚构的部分,因此你得到了镜像效果。
从具有 complex-conjugate 对称性的复杂 frequency-domain 图像开始,应用一维 DFT(例如,沿行)。此转换的结果是一个复杂的图像,仅沿列具有 complex-conjugate 对称性。下一个 1D DFT 获取这些列并从中创建 real-values 信号,产生 real-valued 图像。
如果在第一步后丢弃虚部,则沿该维度删除 spatial-domain 图像的奇数部分,留下偶数(对称)图像。这就是您在结果中看到这种对称性的原因。
这样想:逆DFT反转了DFT的过程。如果 DFT 为 real->complex->complex
,则逆必须为 complex->complex->real
。中间图像具有空间和频率维度。频率分量总是需要很复杂。
您已经知道您需要交换象限。这是因为DFT在空间域和频域都以left-most像素为原点。因此,您的 spatial-domain 卷积核的原点必须位于 top-left 角,所有这些才能正常工作。
最后,您想知道噪音。这是由数值不稳定引起的。图像的某些频率分量只是具有非常小的值。对于 low-pass 过滤图像尤其如此。在那里你将非常小的值除以其他非常小的值,产生无意义。
反卷积在实践中总是需要正则化。维纳滤波器是向反卷积添加正则化的最简单方法。有关详细信息,请参阅 。
我正在尝试对图像 I 执行反卷积,即 n x m。用于对其进行卷积的内核是 K,它也是 n x m。 现在我想通过执行反卷积找到原始图像 O。我知道我可以通过对 I 和 K 执行傅里叶变换并除法来检索图像 O :I / K,因为在傅立叶域中卷积是一个乘积。 (我从 here 那里得到了这些信息)。
我看到另一个关于如何使用 Eigen FFT 执行前向变换的 post,here。
我的正向傅立叶变换代码是:
I = 输入图像(时域)
O = 输出图像(频域)
tempFreq = 用于计算的临时矩阵(频域)
timevec1 = 浮点矢量
freqvec1、freqvec2、freqvec3 = 复向量
for (Int32 i = 0; i < I->uRowsCount; ++i)
{
for (Int32 j = 0; j < I->uColumnsCount; ++j)
{
timevec1.push_back((*I)(i, j));
}
fft.fwd(freqvec1, timevec1);
for (Int32 j = 0; j < I->uColumnsCount; ++j)
{
(tempFreq)(i, j) = freqvec1[j];
}
freqvec1.clear();
timevec1.clear();
}
freqvec1.clear();
timevec1.clear();
for (Int32 j = 0; j < I->uColumnsCount; ++j)
{
for (Int32 i = 0; i < I->uRowsCount; ++i)
{
freqvec2.push_back((tempFreq)(i, j));
}
fft.fwd(freqvec1, freqvec2);
for (Int32 i = 0; i < I->uRowsCount; ++i)
{
(O)(i, j) = freqvec1[i];
}
freqvec2.clear();
freqvec1.clear();
}
我的傅里叶逆变换代码是:
I = 输入图像(频域)
O = 输出图像(时域)
tempTime = 用于计算的临时矩阵(时域)
timevec1, timevec2 = 浮点矢量
freqvec1, freqvec2 = 复向量
for (Int32 j = 0; j < O->uColumnsCount; ++j)
{
for (Int32 i = 0; i < O->uRowsCount; ++i)
{
freqvec1.push_back((I)(i, j));
}
fft.inv(timevec1, freqvec1);
for (Int32 i = 0; i < O->uRowsCount; ++i)
{
(*tempTime)(i, j) = timevecCol[i];
}
freqvec1.clear();
timevec1.clear();
}
freqvec1.clear();
timevec1.clear();
for (Int32 i = 0; i < O->uRowsCount; ++i)
{
for (Int32 j = 0; j < O->uColumnsCount; ++j)
{
freqvec2.push_back((*tempTime)(i, j));
}
fft.inv(timevec2, freqvec2);
for (Int32 j = 0; j < O->uColumnsCount; ++j)
{
(*O)(i, j) = timevec2[j];
}
freqvec2.clear();
timevec2.clear();
}
对于反卷积,我将频域中的输入图像与频域中的内核分开:
freqDomainOutputImage = freqDomainInputImage.cwiseQuotient(freqDomainKernel);
为了获得原始图像,我对 freqDomainOutputImage 执行了逆傅立叶变换。
我得到的结果是:
我相信 FFT 将左上角镜像到另一边,但我不知道为什么?我没有使用 halfSpectrum。第二,为什么图像偏移了?如果我通过用这个替换最后一个循环将输出图像移动到中心:
(*O)((i + O->uRowsCount/2)%O->uRowsCount, (j + O->uColumnsCount/2)%O->uColumnsCount) = timevec2[j];
那么我的输出是:
(可以看到图像是从左上象限镜像的)。
最后,为什么我自己添加了没有噪点的模糊,但它似乎有噪点?
tempTime
需要很复杂。你似乎扔掉了那个虚构的部分,因此你得到了镜像效果。
从具有 complex-conjugate 对称性的复杂 frequency-domain 图像开始,应用一维 DFT(例如,沿行)。此转换的结果是一个复杂的图像,仅沿列具有 complex-conjugate 对称性。下一个 1D DFT 获取这些列并从中创建 real-values 信号,产生 real-valued 图像。
如果在第一步后丢弃虚部,则沿该维度删除 spatial-domain 图像的奇数部分,留下偶数(对称)图像。这就是您在结果中看到这种对称性的原因。
这样想:逆DFT反转了DFT的过程。如果 DFT 为 real->complex->complex
,则逆必须为 complex->complex->real
。中间图像具有空间和频率维度。频率分量总是需要很复杂。
您已经知道您需要交换象限。这是因为DFT在空间域和频域都以left-most像素为原点。因此,您的 spatial-domain 卷积核的原点必须位于 top-left 角,所有这些才能正常工作。
最后,您想知道噪音。这是由数值不稳定引起的。图像的某些频率分量只是具有非常小的值。对于 low-pass 过滤图像尤其如此。在那里你将非常小的值除以其他非常小的值,产生无意义。
反卷积在实践中总是需要正则化。维纳滤波器是向反卷积添加正则化的最简单方法。有关详细信息,请参阅