傅里叶变换之工业缺陷检测应用

前言

目前，网上关于傅里叶变换的文章，主要是针对基础原理和一维信号滤波的讲解，对于二维图像处理则少有涉及。由于二维图像在频域处理的过程中，不仅需要做傅里叶变换，也会涉及到 幅度谱、相位谱、低频信号、高频信号、中心化操作、傅里叶反变换 等一系列概念，只有对这些概念和操作过程比较清楚，才能根据自己的需求得到较好的处理结果。本文就是针对这个问题进行创作的。

傅里叶变换在图像处理的应用主要包含一下几个方面：

1. 在图像高低通滤波和选择性滤波中的应用
2. 在图像压缩中的应用
3. 在图像增强中的应用

欢迎大家的评论，我们互相讨论学习，今后将针对这一些列问题继续发文介绍。

一、问题概述

在工业缺陷检测的场景中，我们经常会遇到一些难以直接通过空域方法进行检测的现象。如图所示，在进行织布赃污的检测时，由于脏污区域与周围像素点的灰度值相差较小，且脏污区域内的像素灰度值变化较大，因此很难选择合适的阈值来进行图像二值化，从而确定脏污区域。

此时，我们就要在频域内对图像进行处理。通过频域滤波，滤除掉频率较低的成分（无脏污部分）和频率较高的成分（织布纹理图案部分），即可得到脏污区域的图像。

二、二维图像的傅里叶变换

结合 欧拉公式 可知，二维傅里叶变换结果包含 实数 R(u,v) 和 虚数 I(u,v) 两部分，从而可以根据R(u,v) 和 I(u,v) 计算傅里叶变换结果的幅度谱和相位谱，如下图所示：

可得，可见，经过傅里叶变换后，得到了两个与原图像同样大小的二维图像，分别表示 实数 和 虚数 部分。再根据幅频谱和相位谱的计算公式，得到为 幅度谱（magnitude） 和 相位谱（phase）

C++ 结合 opencv 实现代码：

	/*对图像imgSrc进行傅里叶变换*/ 
    //扩展图像到合适尺寸以加快Fourier变化的速度
    Mat srcPadded;
    auto m = getOptimalDFTSize(imgSrc.rows);
    auto n = getOptimalDFTSize(imgSrc.cols);
    copyMakeBorder(imgSrc, srcPadded, 0, m - imgSrc.rows, 0, n - imgSrc.cols, BORDER_REPLICATE, Scalar(0));
    
    //提升像素精度并增加虚部通道以方便变换
    srcPadded.convertTo(srcPadded, CV_32FC1);
    Mat plane[] = { srcPadded, Mat::zeros(srcPadded.size(), srcPadded.type()) };
    Mat imgComplex;
    merge(plane, 2, imgComplex);
    
    //执行傅里叶变换
    dft(imgComplex, imgComplex);

结果如下图所示：

由于傅里叶变换后，图像中对应的低频成分在四角区域，高频成分在靠近中心区域（可以查阅相关资料进一步了解其计算原理）。为了便于滤波操作，通常需要对 幅度谱 和 相位谱 分别进行中心化处理，将低频成分集中在中心区域。

C++ 结合 opencv 实现代码：

 //频率中心化
 split(imgComplex, plane);  

 //实部操作
 Mat temp;
 int cx = plane[0].cols / 2;
 int cy = plane[0].rows / 2;
 Mat part1_r(plane[0], Rect(0, 0, cx, cy)); Mat part2_r(plane[0], Rect(cx, 0, cx, cy));
 Mat part3_r(plane[0], Rect(0, cy, cx, cy)); Mat part4_r(plane[0], Rect(cx, cy, cx, cy));
 part1_r.copyTo(temp); part4_r.copyTo(part1_r); temp.copyTo(part4_r);
 part2_r.copyTo(temp); part3_r.copyTo(part2_r); temp.copyTo(part3_r);
 //虚部操作
 Mat part1_i(plane[1], Rect(0, 0, cx, cy)); Mat part2_i(plane[1], Rect(cx, 0, cx, cy));
 Mat part3_i(plane[1], Rect(0, cy, cx, cy)); Mat part4_i(plane[1], Rect(cx, cy, cx, cy));
 part1_i.copyTo(temp); part4_i.copyTo(part1_i); temp.copyTo(part4_i);
 part2_i.copyTo(temp); part3_i.copyTo(part2_i); temp.copyTo(part3_i);

结果如下图所示：

接下来就可以进行滤波操作了。通常情况下，我们可以生成一个与原图同样大小的高斯核函数，与中心化处理的图像进行点乘操作，即可减少许多高频成分（即低通滤波）。再将实部和虚部的滤波结果合并，进行傅里叶反变换。

C++ 结合 opencv 实现代码：
GaussianLPFilter（…）等高斯滤波函数的实现方式见文末代码

    //根据输入滤波器类型生成相应的滤波器
    Mat filter;
    switch (filter_type) {
    case 1:
        filter = GaussianLPFilter(m, n, sigma1, false, rtype);    //
        break;
    case 2:
        filter = GaussianHPFilter(m, n, sigma1, false, rtype);
        break;
    case 3:
        filter = GaussianBPFilter(m, n, sigma2, sigma1, false, rtype);
        break;
    case 4:
        filter = GaussianBSFilter(m, n, sigma2, sigma1, false, rtype);
        break;
    default:
        filter = GaussianLPFilter(m, n, sigma1, false, rtype);
        break;
    }

    //实部、虚部分别和滤波器相乘
    Mat imgBlured_r, imgBlured_i;
    multiply(plane[0], filter, imgBlured_r);
    multiply(plane[1], filter, imgBlured_i);

    //实部、虚部合并
    Mat plane1[] = { imgBlured_r, imgBlured_i };
    Mat imgFreqBlured;
    merge(plane1, 2, imgFreqBlured);

    //反变换到空间域，求幅值并且转换到原图像素类型
    Mat result;    
    idft(imgFreqBlured, imgFreqBlured, DFT_SCALE);
    split(imgFreqBlured, plane);
    magnitude(plane[0], plane[1], result);
    result.convertTo(plane[0], rtype);     

结果如下图所示：

傅里叶反变换结果包含了实部和虚部两个二维图像，计算其幅度谱后，即可得到低通滤波结果。

为什么用反变换后的幅度谱来表示低通滤波结果，不考虑使用相位谱呢？幅度谱和相位谱分别代表什么含义呢？这个可以参考刘定生老师在《数字图像处理》中的ppt讲解（见下图）。这个涉及到较深的傅里叶变换原理的内容，能不能搞明白就看大家的能力了，我也没有很清楚

三、检测案例及代码分析

回到开头提到的织布脏污的检测，对于织布图像的纹理细节而言，其属于高频区域，对于没有脏污的区域而言，又属于低频区域。对图像分别进行低通，高通，带通滤波后，结果如下图所示：

低通滤波结果中，看不到脏污信息；
高通滤波结果中，能看到部分脏污信息，同时，纹理区域也被识别出来；
带通滤波结果中，能够明显地区分出脏污区域。

频域滤波函数代码如下，这里使用高斯核进行滤波操作。

namespace InspectCV
{

Mat GaussianLPFilter(int rows, int cols, double sigma_x, bool bSameSigma, int ktype)
{
    //高斯低通滤波器

    double sigma_y = 0.0;
    if (bSameSigma)
        sigma_y = sigma_x;
    else
        sigma_y = sigma_x * rows / cols;
    rows = getOptimalDFTSize(rows);
    cols = getOptimalDFTSize(cols);
    Mat kernelX = getGaussianKernel(cols, sigma_x, ktype);
    Mat kernelY = getGaussianKernel(rows, sigma_y, ktype);
    Mat kernel = kernelY * kernelX.t();
    normalize(kernel, kernel, 0, 1, NORM_MINMAX);
    return kernel;
}

Mat GaussianHPFilter(int rows, int cols, double sigma_x, bool bSameSigma, int ktype)
{
    //高斯高通滤波器

    Mat hpf = 1 - GaussianLPFilter(rows, cols, sigma_x, bSameSigma, ktype);
    //normalize(hpf, hpf, 0.0, 1.0, CV_MINMAX);
    return hpf;
}

Mat GaussianBPFilter(int rows, int cols, double sigma_x_high, double sigma_x_low, bool bSameSigma, int ktype)
{
    //高斯带通滤波器

    Mat filter_high = GaussianLPFilter(rows, cols, sigma_x_high, bSameSigma, ktype);
    Mat filter_low = GaussianLPFilter(rows, cols, sigma_x_low, bSameSigma, ktype);
    Mat bpf = filter_high - filter_low;
    normalize(bpf, bpf, 0.0, 1.0, NORM_MINMAX);
    return bpf;
}

Mat GaussianBSFilter(int rows, int cols, double sigma_x_high, double sigma_x_low, bool bSameSigma, int ktype)
{
    //高斯带阻滤波器

    return 1 - GaussianBPFilter(rows, cols, sigma_x_high, sigma_x_low, bSameSigma, ktype);
}

Mat GaussFreqBlur(const Mat& imgSrc, double sigma1, double sigma2, int filter_type, int rtype)
{
    /*
        filter_type == 1, 低通滤波；
        filter_type == 2, 高通滤波；
        filter_type == 3, 带通滤波；
        filter_type == 4, 带阻滤波；
    */
    
    //扩展图像到合适尺寸以加快Fourier变化的速度
    Mat srcPadded;
    auto m = getOptimalDFTSize(imgSrc.rows);
    auto n = getOptimalDFTSize(imgSrc.cols);
    copyMakeBorder(imgSrc, srcPadded, 0, m - imgSrc.rows, 0, n - imgSrc.cols, BORDER_REPLICATE, Scalar(0));
    
    //提升像素精度并增加虚部通道以方便变换
    srcPadded.convertTo(srcPadded, CV_32FC1);
    Mat plane[] = { srcPadded, Mat::zeros(srcPadded.size(), srcPadded.type()) };
    Mat imgComplex;
    merge(plane, 2, imgComplex);
    
    //正向变换
    dft(imgComplex, imgComplex);

    //频率中心化
    split(imgComplex, plane);  
   
    //实部操作
    Mat temp;
    int cx = plane[0].cols / 2;
    int cy = plane[0].rows / 2;
    Mat part1_r(plane[0], Rect(0, 0, cx, cy)); Mat part2_r(plane[0], Rect(cx, 0, cx, cy));
    Mat part3_r(plane[0], Rect(0, cy, cx, cy)); Mat part4_r(plane[0], Rect(cx, cy, cx, cy));
    part1_r.copyTo(temp); part4_r.copyTo(part1_r); temp.copyTo(part4_r);
    part2_r.copyTo(temp); part3_r.copyTo(part2_r); temp.copyTo(part3_r);
    //虚部操作
    Mat part1_i(plane[1], Rect(0, 0, cx, cy)); Mat part2_i(plane[1], Rect(cx, 0, cx, cy));
    Mat part3_i(plane[1], Rect(0, cy, cx, cy)); Mat part4_i(plane[1], Rect(cx, cy, cx, cy));
    part1_i.copyTo(temp); part4_i.copyTo(part1_i); temp.copyTo(part4_i);
    part2_i.copyTo(temp); part3_i.copyTo(part2_i); temp.copyTo(part3_i);

    //根据输入滤波器类型生成相应的滤波器
    Mat filter;
    switch (filter_type) {
    case 1:
        filter = GaussianLPFilter(m, n, sigma1, false, rtype);
        break;
    case 2:
        filter = GaussianHPFilter(m, n, sigma1, false, rtype);
        break;
    case 3:
        filter = GaussianBPFilter(m, n, sigma2, sigma1, false, rtype);
        break;
    case 4:
        filter = GaussianBSFilter(m, n, sigma2, sigma1, false, rtype);
        break;
    default:
        filter = GaussianLPFilter(m, n, sigma1, false, rtype);
        break;
    }

    //实部、虚部分别和滤波器相乘
    Mat imgBlured_r, imgBlured_i;
    multiply(plane[0], filter, imgBlured_r);
    multiply(plane[1], filter, imgBlured_i);

    //实部、虚部合并
    Mat plane1[] = { imgBlured_r, imgBlured_i };
    Mat imgFreqBlured;
    merge(plane1, 2, imgFreqBlured);

    //反变换到空间域，求幅值并且转换到原图像素类型
    idft(imgFreqBlured, imgFreqBlured, DFT_SCALE);
    split(imgFreqBlured, plane);
    magnitude(plane[0], plane[1], plane[0]);
    plane[0].convertTo(plane[0], rtype);

    //剪裁到输入图像尺寸大小并返回
    return plane[0](Rect(0, 0, imgSrc.cols, imgSrc.rows));
} 

}