基于OpenCV的同态滤波

在4.0.1节中，我们已经介绍了一个简单的图像形成模型，即照射-反射模型。这个模型可以开发一种频率处理程序，该程序可以同时压缩灰度范围和增强对比度来改善一幅图像的表现。图像形成的照射-反射模型的表达式如下：

                                         6.4- 16

上式不能直接用于对照射和反射的频率分量进行操作，因为两个分量的傅里叶变换之积，不等于原图像的傅里叶变换：

                                6.4- 17

但是，我们可以定义如下变换关系：

                    6.4- 18

然后再对上式两边做傅里叶变换：

       6.4- 19

得到下面的傅里叶变换域的等式：

                                    6.4- 20

其中， 分别是 的傅里叶变换。

我们可以用一个滤波器 对 滤波，故有：

         6.4- 21

在空间域中，滤波后的图像为：

     6.4- 22

我们可以将上式改写成简略形式：

                                    6.4- 23

其中，

                                  6.4- 24

                                 6.4- 25

最后，因为 是通过取输入图像的自然对数形成的，我们可通过取滤波后的结果的指数这一反处理来形成输出图像：

                6.4- 26

                        6.4- 27

                                            6.4- 28

经滤波处理后的照射和反射分量。

以上推到的滤波方法流程图，如图6. 26所示。该滤波方法是针对特定成像系统，在此称为同态系统。在这种特殊应用中，方法的关键在于照射分量和反射分量的分离，其实现形式（6.4-20）所示。然后，如式（6.4-21）那样，用同态滤波器 对这些分量进行滤波操作。

图像的照射分量通常由慢的空间变化来表征，而反射分量往往引起突变，特别是在不同物体的连接部分。这些特性导致图像取对数后的傅里叶变换的低频成分与照射分量相联系，而高频成分与反射分量相联系。虽然这些联系只是粗略的近似，但它们用在图像滤波中是有益的。

使用同态滤波器可更好地控制照射分量和反射分量。这种控制需要指定一个滤波器函数 ，它可用不同的可控方法影响傅里叶变换的低频和高频分量。图6. 27显示了这种滤波器的剖面图。如果 和 选定， 且 ，那么图6. 27所示滤波器函数趋向于衰减低频（照射）的贡献。

图6. 27所示的函数形状可用高通滤波器的基本形式来近似。例如，采用形式稍微变化一下的高斯高通滤波器可得到函数：

                                                                        6.4- 29

其中 由式（6.3-2）定义，常数c控制函数坡度的锐利度，它在 和 之间过渡。这个滤波器类似于前面6.4.3节讨论的高频强调滤波器。

例6. 10 使用同态滤波增强图像

图6. 28（a）显示了一幅大小为1162×746像素的全身PET扫描图像。图像稍微有些模糊，并且由于主导显示动态范围的高灰度“热点”使得其许多低灰度特征很朦胧（这些“热点”是由脑部的肿瘤和肺部的肿瘤导致的）。

图6. 28（b）是图6. 28（a）经同态滤波湖得到的结果，使用式（6.4-8）中的滤波器， 。改滤波器的剖面看上去正好像图6. 27，只是坡度有点陡峭。

Mat homoMorphicFilter(const Mat& src,double D0,
                      double gammaH,double gammaL,double c)
{
    int rows=src.rows,cols=src.cols,channels=src.channels();
    Mat srcLog,src64d;
    src.convertTo(src64d,CV_64FC1);
    cv::log(src64d+1,srcLog);
    Mat srcLogFft=fft(srcLog);
    dftshift(srcLogFft);

    //创建同态滤波器
    Mat Huv(rows,cols,CV_64FC2,Scalar::all(0));//频域滤波器
    int cx=cols/2,cy=rows/2;
    double D2=D0*D0;
    for(int i=0;i<rows;i++)
    {
        Vec2d* p=Huv.ptr<Vec2d>(i);
        for(int j=0;j<cols;j++)
        {
            double x=j-cx,y=i-cy;
            p[j][0]=(gammaH-gammaL)*(1-c*exp(-(x*x+y*y)/D2))+gammaL;
            p[j][1]=p[j][0];
        }
    }
    Huv /=(rows*cols);

    //频域滤波，在频域与Huv逐点相乘
    Mat ftProduct=srcLogFft.mul(Huv);
    //逆傅里叶变换
    Mat ftInv=fft(ftProduct,DFT_INVERSE);
    //计算逆傅里叶变换的幅值
    Mat mag= myMagnitude(ftInv);
    Mat edst;
    cv::exp(mag,edst);
    normalize(edst,edst,255,0,NORM_MINMAX);
    edst.convertTo(edst,CV_8UC1);
    return edst;
}

下面是主程序：

int main()
{
    Mat img=imread("D:/CodeWork/MyImage/Fig0462a.tif",0);
    imshow("rogin img",img);
    Mat imgHomo=homoMorphicFilter(img,80,3,0.25,1.);
    imshow("homoMorphic filter",(imgHomo+img)/2);
    waitKey();
    return 0;
}

下面左图是经过同态滤波后的结果，右图是原图像：