C语言数字图像处理（七）：图像噪声&降噪算法，同态滤波器和带阻滤波器

0. 完整仓库 & 教程：

这一章节的完整代码在：Chapter 7. Image Noise and Noise Reduction

如果你喜欢这个系列的文章或者感觉对你有帮助，请给我的仓库一个⭐️。

1. 同态滤波器

算法：

同态滤波器的过程与第六章(这部分可以参考：6：低通滤波器和高通滤波器)中的高通/低通滤波器类似：

首先对图像进行二维离散傅里叶变换 DFT (这部分原理需要理解： 5: 二维离散傅里叶变换和重建)，得到每个像素的复数数组。 每个复数的实部和虚部分别乘以 𝐻(𝑢, 𝑣)。 𝐻(𝑢, 𝑣) 的公式为：

其中 𝛾𝐻 > 1， 𝛾𝐿 < 1，并且 表示通带的半径。 𝐷(𝑢, 𝑣)的计算方法也是两点之间的距离，通过公式：

其中 𝑃 和𝑄 是图像的长和宽。在后续的代码处理中，设置𝛾𝐻 = 1.5, 𝛾𝐿 = 0.75,  = 30, 𝑐 = 1。

图像和结果对比（bridge, goldhill）：

结果分析：

同态滤波本质上是频域滤波器（frequency domain filter），它结合了频率滤波和空间灰度变换。 通常，图像的亮部（illumination）是缓慢变化的部分，集中在低频部分。 反射部分（reflection）属于快速变化的部分，往往集中在高频部分。

这里选用模型的形状类似于一个逆高斯滤波器（见章节6），可以减弱低频分量（illumination）的影响并强调高频分量（reflection）。 最终的结果是压缩了动态范围和增强了对比度，使纹理更加突出。

代码实现(完整代码见顶部GitHub)：

for(int i = 0; i < height; i++) {
    for(int j = 0; j < width; j++) {
        float dis = sqrt(pow((float)i - height/2, 2) + pow((float)j - width/2, 2));
        float H = (1.5 - 0.75) * (1 - pow(M_E, pow(dis / 30, 2) * -1)) + 0.75;
        real[(int)(i*height + j)] = real_array[(int)(i*height + j)] * H;
        imaginary[(int)(i*height + j)] = imaginary_array[(int)(i*height + j)] * H;
    }
}
printf("Homomorphic Finished!\n");
outimage = iDFT(image, real, imaginary);

这里再把第5章里傅里叶变换和重建的代码搬过来下：

void DFT(Image *image, float *real_array, float *imaginary_array) {
    unsigned char *tempin, *tempout;
    float real, imaginary;
    Image *outimage;
    outimage = CreateNewImage(image, (char*)"#testing function");
    tempin = image->data;
    tempout = outimage->data;
    
    printf("DFT algorithm is executing...\n");
    for(int i = 0; i < image->Height; i++) {
        for(int j = 0; j < image->Width; j++) {
            real = 0;
            imaginary = 0;
            
            for(int m = 0; m < image->Height; m++) {
                for(int n = 0; n < image->Width; n++) {
                    float temp = (float)i * m / (float)image->Height + (float)j * n / (float)image->Width;
                    int offset = (m + n) % 2 == 0 ? 1 : -1;
                    real += tempin[image->Width * m + n] * cos(-2 * pi * temp) * offset;
                    imaginary += tempin[image->Width * m + n] * sin(-2 * pi * temp) * offset;
                }
            }
            real_array[image->Width * i + j] = real;
            imaginary_array[image->Width * i + j] = imaginary;

            int temp = (int)(sqrt(real*real + imaginary*imaginary) / sqrt((float)image->Height*(float)image->Width));
            if(temp < 0) temp = 0;
            if(temp > 255) temp = 255;
            tempout[image->Width * i + j] = temp;
        }
    }
    printf("DFT Finished!\n");
    SavePNMImage(outimage, (char*)"DFT.pgm");
}

Image *iDFT(Image *image, float *real_array, float *imaginary_array) {
    unsigned char *tempin, *tempout;
    float real;
    Image *outimage;
    outimage = CreateNewImage(image, (char*)"#testing fun