c语言数字图像处理（八）：噪声模型及均值滤波器

图像退化/复原过程模型

高斯噪声

PDF（概率密度函数）

生成高斯随机数序列

算法可参考<http://www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf>

代码实现

double gaussrand()
{
    static double V1, V2, S;
    static int phase = 0;
    double X;

    if(phase == 0) {
        do {
            double U1 = (double)rand() / RAND_MAX;
            double U2 = (double)rand() / RAND_MAX;

            V1 = 2 * U1 - 1;
            V2 = 2 * U2 - 1;
            S = V1 * V1 + V2 * V2;
            } while(S >= 1 || S == 0);

        X = V1 * sqrt(-2 * log(S) / S);
    } else
        X = V2 * sqrt(-2 * log(S) / S);

    phase = 1 - phase;

    return X * 50;
}

生成高斯噪声图及直方图

下面给一幅图添加高斯噪声

代码实现

void add_gaussian_noise(short** in_array, short** out_array, long height, long width)
{
    srand(time(NULL));
    for (int i = 0; i < height; i++){
        for (int j = 0; j < width; j++){
            out_array[i][j] = in_array[i][j] + (short)gaussrand();
            if (out_array[i][j] < 0x00)
                out_array[i][j] = 0x00;
            else if (out_array[i][j] > 0xff)
                out_array[i][j] = 0xff;
        }
    }
}

原图

添加高斯噪声

椒盐噪声

添加椒盐噪声（胡椒噪声和盐粒噪声概率分别为5%)

void add_salt_pepper_noise(short** in_array, short** out_array, long height, long width)
{
    srand(time(NULL));
    int noise_p;

    for (int i = 0; i < height; i++){
        for (int j = 0; j < width; j++){
            noise_p = rand() % 10;
            if (noise_p == 0){
                int temp = rand() % 2;
                if (temp)
                    out_array[i][j] = 0x00;
                else
                    out_array[i][j] = 0xff;
            }
            else
                out_array[i][j] = in_array[i][j];
        }
    }
}

均值滤波器

算术均值滤波器

代码实现

int is_in_array(short x, short y, short height, short width)
{
    if (x >= 0 && x < width && y >= 0 && y < height)
        return 1;
    else
        return 0;
}

/*
 * element
 * v0  v1  v2
 * v3  v4  v5
 * v6  v7  v8
 *
 */
void filtering(short** in_array, short** out_array, long height, long width)
{
    short value[9];

    for (int i = 0; i < height; i++){
        for (int j = 0; j < width; j++){
            value[0] = is_in_array(j-1, i-1, height, width) ? in_array[i-1][j-1] : 0;
            value[1] = is_in_array(j, i-1, height, width) ? in_array[i-1][j] : 0;
            value[2] = is_in_array(j+1, i-1, height, width) ? in_array[i-1][j+1] : 0;
            value[3] = is_in_array(j-1, i, height, width) ? in_array[i][j-1] : 0;
            value[4] = in_array[i][j];
            value[5] = is_in_array(j+1, i, height, width) ? in_array[i][j+1] : 0;
            value[6] = is_in_array(j-1, i+1, height, width) ? in_array[i+1][j-1] : 0;
            value[7] = is_in_array(j, i+1, height, width) ? in_array[i+1][j] : 0;
            value[8] = is_in_array(j+1, i+1, height, width) ? in_array[i+1][j+1] : 0;

            /* Arithmetic Mean Filter */
            for (int k = 0; k < ARRAY_SIZE*ARRAY_SIZE; k++)
                out_array[i][j] += value[k];
            out_array[i][j] /= ARRAY_SIZE*ARRAY_SIZE;

        }
    }
}

处理高斯噪声

处理椒盐噪声

结论：算术平均滤波对于高斯噪声和椒盐噪声都有一定的效果，但是同时会平滑图像

几何均值滤波器

实现

void filtering(short** in_array, short** out_array, long height, long width)
{
    short value[9];

    for (int i = 0; i < height; i++){
        for (int j = 0; j < width; j++){
            value[0] = is_in_array(j-1, i-1, height, width) ? in_array[i-1][j-1] : 0;
            value[1] = is_in_array(j, i-1, height, width) ? in_array[i-1][j] : 0;
            value[2] = is_in_array(j+1, i-1, height, width) ? in_array[i-1][j+1] : 0;
            value[3] = is_in_array(j-1, i, height, width) ? in_array[i][j-1] : 0;
            value[4] = in_array[i][j];
            value[5] = is_in_array(j+1, i, height, width) ? in_array[i][j+1] : 0;
            value[6] = is_in_array(j-1, i+1, height, width) ? in_array[i+1][j-1] : 0;
            value[7] = is_in_array(j, i+1, height, width) ? in_array[i+1][j] : 0;
            value[8] = is_in_array(j+1, i+1, height, width) ? in_array[i+1][j+1] : 0;

            /* Geometric Mean Filter */
            double product = 1.0;
            for (int k = 0; k < ARRAY_SIZE*ARRAY_SIZE; k++)
                product *= value[k];
            product = pow(product, 1.0 / 9.0);
            out_array[i][j] = (short)product;

            if (out_array[i][j] < 0x00)
                out_array[i][j] = 0x00;
            else if (out_array[i][j] > 0xff)
                out_array[i][j] = 0xff;
        }
    }
}

几何均值滤波器与算术均值滤波器相比，丢失的图像细节更少

谐波均值滤波器

实现

void filtering(short** in_array, short** out_array, long height, long width)
{
    short value[9];

    for (int i = 0; i < height; i++){
        for (int j = 0; j < width; j++){
            value[0] = is_in_array(j-1, i-1, height, width) ? in_array[i-1][j-1] : 0;
            value[1] = is_in_array(j, i-1, height, width) ? in_array[i-1][j] : 0;
            value[2] = is_in_array(j+1, i-1, height, width) ? in_array[i-1][j+1] : 0;
            value[3] = is_in_array(j-1, i, height, width) ? in_array[i][j-1] : 0;
            value[4] = in_array[i][j];
            value[5] = is_in_array(j+1, i, height, width) ? in_array[i][j+1] : 0;
            value[6] = is_in_array(j-1, i+1, height, width) ? in_array[i+1][j-1] : 0;
            value[7] = is_in_array(j, i+1, height, width) ? in_array[i+1][j] : 0;
            value[8] = is_in_array(j+1, i+1, height, width) ? in_array[i+1][j+1] : 0;

            /* Harmonic Mean Filter */
            double sum = 0;
            for (int k = 0; k < ARRAY_SIZE*ARRAY_SIZE; k++)
                sum += 1.0 / value[k];
            out_array[i][j] = (short)(9.0 / sum);

            if (out_array[i][j] < 0x00)
                out_array[i][j] = 0x00;
            else if (out_array[i][j] > 0xff)
                out_array[i][j] = 0xff;
        }
    }
}

处理高斯噪声

处理椒盐噪声

对盐粒噪声效果较好，不适用于胡椒噪声，善于处理高斯噪声

逆谐波均值滤波器

Q为滤波器的阶数，Q为正时，消除胡椒噪声，Q为负时消除盐粒噪声

Q=0为算术均值滤波器，Q=-1谐波均值滤波器

实现

void filtering(short** in_array, short** out_array, long height, long width)
{
    short value[9];

    for (int i = 0; i < height; i++){
        for (int j = 0; j < width; j++){
            value[0] = is_in_array(j-1, i-1, height, width) ? in_array[i-1][j-1] : 0;
            value[1] = is_in_array(j, i-1, height, width) ? in_array[i-1][j] : 0;
            value[2] = is_in_array(j+1, i-1, height, width) ? in_array[i-1][j+1] : 0;
            value[3] = is_in_array(j-1, i, height, width) ? in_array[i][j-1] : 0;
            value[4] = in_array[i][j];
            value[5] = is_in_array(j+1, i, height, width) ? in_array[i][j+1] : 0;
            value[6] = is_in_array(j-1, i+1, height, width) ? in_array[i+1][j-1] : 0;
            value[7] = is_in_array(j, i+1, height, width) ? in_array[i+1][j] : 0;
            value[8] = is_in_array(j+1, i+1, height, width) ? in_array[i+1][j+1] : 0;

            /* Contra-Harmonic Mean Filter */
            int Q = 2;
            double num = 0.0, den = 0.0;
            for (int k = 0; k < ARRAY_SIZE*ARRAY_SIZE; k++){
                num += pow(value[k], Q+1);
                den += pow(value[k], Q);
            }
            out_array[i][j] = (short)(num / den);

            if (out_array[i][j] < 0x00)
                out_array[i][j] = 0x00;
            else if (out_array[i][j] > 0xff)
                out_array[i][j] = 0xff;
        }
    }
}

Q = 2 消除胡椒噪声

Q = -2消除盐粒噪声

Q = -2消除盐粒噪声后的图像使用Q = 2消除胡椒噪声

再来一次

再来

此时椒盐噪声已经基本消除