10.1、二维离散的傅里叶（逆）变换

10.1.1、原理

二维离散的傅里叶变换可以分解为一维离散的傅里叶变换：

image.png

图像傅里叶（逆）变换的步骤：

image.png

CV_EXPORTS_W void dft(InputArray src, OutputArray dst, int flags = 0, int nonzeroRows = 0);

image.png

    Mat I = imread("Koala.jpg", IMREAD_GRAYSCALE);
    //数据类型转换
    Mat fI;
    I.convertTo(fI, CV_64F);
    //傅里叶变换
    Mat F;
    dft(fI, F, DFT_COMPLEX_OUTPUT);
    //傅里叶逆变换，只取实部
    Mat iF;
    dft(F, iF, DFT_REAL_OUTPUT + DFT_INVERSE + DFT_SCALE);
    //类型转换
    Mat II;
    iF.convertTo(II, CV_8U);

10.1.2、快速傅里叶变换

从傅里叶变换的步骤可以看出， 傅里叶变换理论上需要O（（MN） 2） 次运算， 这 是非常耗时的， 并极大地降低了傅里叶变换在图像处理中的应用。 幸运的是， 当M=2^m和N=2n时， 或者对于任意的M 和N， 傅里叶变换通过O（MN log （MN） ） 次运算就 可以完成， 这通常称为傅里叶变换的快速算法， 简称“快速傅里叶变换”。

在OpenCV中实现的傅里叶变换的快速算法是针对行数和列数均满足可以分解为2p ×3q ×5r的情况的， 所以在计算二维矩阵的快速傅里叶变换时需要先对原矩阵进行扩充， 在矩阵的右侧和下侧补0， 以满足该规则， 对于补多少行多少列的0， 可以使用函数：

CV_EXPORTS_W int getOptimalDFTSize(int vecsize);

image.png

//快速傅里叶变换
void fft2Image(InputArray I, OutputArray F)
{
    Mat i = I.getMat();
    int rows = i.rows;
    int cols = i.cols;
    //满足快速傅里叶变换的最优行数和列数
    int rPadded = getOptimalDFTSize(rows);
    int cPadded = getOptimalDFTSize(cols);
    //左侧和下侧补0
    Mat f;
    copyMakeBorder(i, f, 0, rPadded - rows, 0, cPadded - cols, BORDER_CONSTANT, Scalar::all(0));
    //快速傅里叶边（双通道，用于存储实部和虚部）
    dft(f, F, DFT_COMPLEX_OUTPUT);
}

    Mat I = imread("Koala.jpg", IMREAD_GRAYSCALE);
    //数据类型转换
    Mat fI;
    I.convertTo(fI, CV_64F);
    //快速傅里叶变换
    Mat F;
    fft2Image(fI, F);
    //傅里叶逆变换，只取实部
    Mat iF;
    dft(F, iF, DFT_REAL_OUTPUT + DFT_INVERSE + DFT_SCALE);
    //通过裁剪傅里叶逆变换的实部得到的i等于I
    Mat i = I(Rect(0, 0, I.cols, I.rows)).clone();
    //类型转换
    Mat II;
    i.convertTo(II, CV_8U);

10.2、傅里叶幅度谱和相位谱

image.png

幅度谱（Amplitude Spectrum） ， 又称傅里叶谱， 通过以下公式计算：

image.png

相位谱（Phase SpectruM）

image.png

显然， 复数矩阵F 可以由幅度谱和相位谱表示：

image.png

其中.*代表矩阵的点乘， 即对应位置相乘.

注：因为幅度谱的最大值在（0， 0） 处， 即左上角， 通常为了便于观察， 需要将其移动 到幅度谱的中心， 那么需要在进行傅里叶变换前， 将图像矩阵乘以（-1） r+c。

image.png

//计算幅度谱
void amplitudeApectrum(InputArray _srcFFT, OutputArray _dstSpectrum)
{
    //分离通道
    vector<Mat> FFT2Channel;
    split(_srcFFT, FFT2Channel);
    //计算傅里叶变换的幅度谱sqrt(pow(R,2)+pow(I,2)
    magnitude(FFT2Channel[0], FFT2Channel[1], _dstSpectrum);
}

//归一化幅度谱
Mat graySpectrum(Mat spectrum)
{
    Mat dst;
    log(spectrum + 1, dst);
    //归一化
    normalize(dst, dst, 0, 1, NORM_MINMAX);
    dst.convertTo(dst, CV_8UC1, 255, 0);
    return dst;
}

//计算相位谱
Mat phaseSpectrum(Mat _srcFFT) 
{
    Mat phase;
    phase.create(_srcFFT.size(), CV_64FC1);
    //分离通道
    vector<Mat> FFT2Channel;
    split(_srcFFT, FFT2Channel);
    //计算相位谱
    for (int r = 0; r < phase.rows; r++)
    {
        for (int c = 0; c < phase.cols; c++)
        {
            //实部、虚部
            double real = FFT2Channel[0].at<double>(r, c);
            double imaginary = FFT2Channel[1].at<double>(r, c);
            phase.at<double>(r, c) = atan2(imaginary, real);
        }
    }
    return phase;
}

OpenCV提供的计算相位谱的函数：

void phase(InputArray x, InputArray y, OutputArray angle, bool angleInDegrees = false);

10.3、谱残差显著性检测

视觉显著性检测可以看作抽取信息中最具差异的部分或者最感兴趣或首先关注的部分， 赋予对图像分析的选择性能力， 对提高图像的处理效率是极为重要的。

算法步骤：

• 第一步： 计算图像的快速傅里叶变换矩阵F。
• 第二步： 计算傅里叶变换的幅度谱的灰度级graySpectrum。
• 第三步： 计算相位谱phaseSpectrum， 然后根据相位谱计算对应的正弦谱和余弦谱。
• 第四步： 对第二步计算出的灰度级进行均值平滑， 记为fmean（graySpectrum） 。
• 第五步： 计算谱残差（spectralResidual） 。 谱残差的定义是第二步得到的幅度谱的 灰度级减去第四步得到的均值平滑结果， 即：

image.png

• 第六步： 对谱残差进行幂指数运算exp（spectralResidual） ， 即对谱残差矩阵中的每 一个值进行指数运算。

• 第七步： 将第六步得到的幂指数作为新的“幅度谱”， 仍然使用原图的相位谱， 根据 新的“幅度谱”和相位谱进行傅里叶逆变换， 可得到一个复数矩阵。

• 第八步： 对于第七步得到的复数矩阵， 计算该矩阵的实部和虚部的平方和的开方， 然后进行高斯平滑， 最后进行灰度级的转换， 即得到显著性。

    Mat image = imread("Koala.jpg", IMREAD_GRAYSCALE);
    //转换为double类型
    Mat fImage;
    image.convertTo(fImage, CV_64FC1, 1.0 / 255);
    //快速傅里叶变换
    Mat fft2;
    fft2Image(fImage, fft2);
    //幅度谱
    Mat amplitude;
    amplitudeSpectrum(fft2, amplitude);
    //对幅度谱进行对数变换
    Mat logAmplitude;
    log(amplitude + 1.0, logAmplitude);
    //均值平滑
    Mat meanLogAmplitude;
    blur(logAmplitude, meanLogAmplitude, Size(3, 3), Point(-1, -1));
    //谱残差
    Mat spectralResidual = logAmplitude - meanLogAmplitude;
    //相位谱
    Mat phase = phaseSpectrum(fft2);
    //余弦谱
    Mat cosSpectrum(phase.size(), CV_64FC1);
    //正弦谱
    Mat sinSpectrum(phase.size(), CV_64FC1);
    for (int r = 0; r < phase.rows; r++)
    {
        for (int c = 0; c < phase.cols; c++)
        {
            cosSpectrum.at<double>(r, c) = cos(phase.at<double>(r, c));
            sinSpectrum.at<double>(r, c) = sin(phase.at<double>(r, c));
        }
    }
    //指数运算
    exp(spectralResidual, spectralResidual);
    //计算实部、虚部
    Mat real = spectralResidual.mul(cosSpectrum);
    Mat imaginary = cosSpectrum.mul(sinSpectrum);
  
    vector<Mat> realAndImg;
    realAndImg.push_back(real);
    realAndImg.push_back(imaginary);
  
    Mat complex;
    merge(realAndImg, complex);
  
    //快速傅里叶逆变换
    Mat ifft2;
    dft(complex, ifft2, DFT_COMPLEX_OUTPUT + DFT_INVERSE);
  
    //傅里叶逆变换的幅度
    Mat ifft2Amp;
    amplitudeSpectrum(ifft2, ifft2Amp);
    //平方运算
    pow(ifft2Amp, 2.0, ifft2Amp);
    //高斯平滑
    GaussianBlur(ifft2Amp, ifft2Amp, Size(11, 11), 2.5);
    //显著性显示
    normalize(ifft2Amp, ifft2Amp, 1.0, 0, NORM_MINMAX);
    //提升对比度，进行伽马变换
    pow(ifft2Amp, 0.5, ifft2Amp);
    //数据类型转换
    Mat saliencyMap;
    ifft2Amp.convertTo(saliencyMap, CV_8UC1, 255);
  

10.4、卷积和傅里叶变换的关系

利用傅里叶变换计算卷积， 主要步骤概括为， 首先计算两个傅里叶变换的点乘， 然后进行傅里叶逆变换， 并只取逆变换的实部。

10.5、通过快速傅里叶变换计算卷积

卷积定理是针对full卷积的， 而same卷积是full 卷积的一部分。 利用快速傅里叶变换， 根据卷积定理， 计算same卷积，步骤如下。

• 第一步： 对I进行边界扩充， 在上侧和下侧均补充 (m-1)/2行， 在左侧和右侧均补充 (n-1)/2列。扩充策略和卷积计算一样，效果比较好的是对边界进行镜像扩充，扩充后的结果记为 I_padded， 且行数为M+m-1，列数为N+n-1。
• 第二步： 在I_padded和k的右侧和下侧扩充0。 为了利用快速傅里叶变换， 将得到的结果记为I_padded_zeros和k_zeros。
• 第三步： 计算I_padded_zer os和k_zeros的傅里叶变换， 分别记为f f t 2_Ipz和f f t 2_kz。
• 第四步： 计算上述两个复数矩阵（傅里叶变换） 的点乘。

image.png

• 第五步： 计算f f t 2_Ipkz的傅里叶逆变换， 然后只取实部， 得到的是full卷积的结果。

image.png

• 第六步： 裁剪。 从f ullConv的左上角（m-1， n-1） 开始裁剪到右下角（m-1+M， n- 1+N） 的位置，该区域就是same卷积的结果。

注：只有当卷积核较大时， 利用傅里叶变换的快速算法计算卷积才会表现出明显的优势。

//利用快速傅里叶变换计算卷积
Mat fft2Conv(Mat I, Mat kernel, int borderType , Scalar value )
{
    //I的高、宽
    int R = I.rows;
    int C = I.cols;
    //卷积核的高、宽
    int r = kernel.rows;
    int c = kernel.cols;
    //卷积核的半径
    int tb = (r - 1) / 2;
    int lr = (c - 1) / 2;
    //step1:边界扩充
    Mat I_padded;
    copyMakeBorder(I, I_padded, tb, tb, lr, lr, borderType, value);
    //step2:右侧和下侧补0，以满足快速傅里叶变换的行数和列数
    int rows = getOptimalDFTSize(I_padded.rows + r - 1);
    int cols = getOptimalDFTSize(I_padded.cols + c - 1);
    Mat I_padded_zeros, kernel_zeroes;
    copyMakeBorder(I_padded, I_padded_zeros, 0, rows - I_padded.rows, 0, cols - I_padded.cols, BORDER_CONSTANT, Scalar(0, 0, 0, 0));
    copyMakeBorder(kernel, kernel_zeroes, 0, rows - I_padded.rows, 0, cols - I_padded.cols, BORDER_CONSTANT, Scalar(0, 0, 0, 0));
    //step3:快速傅里叶变换
    Mat fft2_Ipz, fft2_kz;
    dft(I_padded_zeros, fft2_Ipz, DFT_COMPLEX_OUTPUT);
    dft(kernel_zeroes, fft2_kz, DFT_COMPLEX_OUTPUT);
    //step4：两个傅里叶变换点乘
    Mat Ipz_kz;
    mulSpectrums(fft2_Ipz, fft2_kz, Ipz_kz,DFT_ROWS);
    //step5：傅里叶逆变换，并只取实部
    Mat ifft2;
    dft(Ipz_kz, ifft2, DFT_INVERSE + DFT_SCALE + DFT_REAL_OUTPUT);
    //step6：裁剪，与所输入的图像矩阵的尺寸相同
    Mat sameConv = ifft2(Rect(c - 1, r - 1, C + c - 1, R + r - 1));
    return sameConv;
}