c语言数字图像处理(八):噪声模型及均值滤波器
图像退化/复原过程模型
高斯噪声
PDF(概率密度函数)
生成高斯随机数序列
算法可参考<http://www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf>
代码实现
1 double gaussrand() 2 { 3 static double V1, V2, S; 4 static int phase = 0; 5 double X; 6 7 if(phase == 0) { 8 do { 9 double U1 = (double)rand() / RAND_MAX; 10 double U2 = (double)rand() / RAND_MAX; 11 12 V1 = 2 * U1 - 1; 13 V2 = 2 * U2 - 1; 14 S = V1 * V1 + V2 * V2; 15 } while(S >= 1 || S == 0); 16 17 X = V1 * sqrt(-2 * log(S) / S); 18 } else 19 X = V2 * sqrt(-2 * log(S) / S); 20 21 phase = 1 - phase; 22 23 return X * 50; 24 }
生成高斯噪声图及直方图
下面给一幅图添加高斯噪声
代码实现
1 void add_gaussian_noise(short** in_array, short** out_array, long height, long width) 2 { 3 srand(time(NULL)); 4 for (int i = 0; i < height; i++){ 5 for (int j = 0; j < width; j++){ 6 out_array[i][j] = in_array[i][j] + (short)gaussrand(); 7 if (out_array[i][j] < 0x00) 8 out_array[i][j] = 0x00; 9 else if (out_array[i][j] > 0xff) 10 out_array[i][j] = 0xff; 11 } 12 } 13 }
原图
添加高斯噪声
椒盐噪声
添加椒盐噪声(胡椒噪声和盐粒噪声概率分别为5%)
1 void add_salt_pepper_noise(short** in_array, short** out_array, long height, long width) 2 { 3 srand(time(NULL)); 4 int noise_p; 5 6 for (int i = 0; i < height; i++){ 7 for (int j = 0; j < width; j++){ 8 noise_p = rand() % 10; 9 if (noise_p == 0){ 10 int temp = rand() % 2; 11 if (temp) 12 out_array[i][j] = 0x00; 13 else 14 out_array[i][j] = 0xff; 15 } 16 else 17 out_array[i][j] = in_array[i][j]; 18 } 19 } 20 }
均值滤波器
算术均值滤波器
代码实现
1 int is_in_array(short x, short y, short height, short width) 2 { 3 if (x >= 0 && x < width && y >= 0 && y < height) 4 return 1; 5 else 6 return 0; 7 } 8 9 /* 10 * element 11 * v0 v1 v2 12 * v3 v4 v5 13 * v6 v7 v8 14 * 15 */ 16 void filtering(short** in_array, short** out_array, long height, long width) 17 { 18 short value[9]; 19 20 for (int i = 0; i < height; i++){ 21 for (int j = 0; j < width; j++){ 22 value[0] = is_in_array(j-1, i-1, height, width) ? in_array[i-1][j-1] : 0; 23 value[1] = is_in_array(j, i-1, height, width) ? in_array[i-1][j] : 0; 24 value[2] = is_in_array(j+1, i-1, height, width) ? in_array[i-1][j+1] : 0; 25 value[3] = is_in_array(j-1, i, height, width) ? in_array[i][j-1] : 0; 26 value[4] = in_array[i][j]; 27 value[5] = is_in_array(j+1, i, height, width) ? in_array[i][j+1] : 0; 28 value[6] = is_in_array(j-1, i+1, height, width) ? in_array[i+1][j-1] : 0; 29 value[7] = is_in_array(j, i+1, height, width) ? in_array[i+1][j] : 0; 30 value[8] = is_in_array(j+1, i+1, height, width) ? in_array[i+1][j+1] : 0; 31 32 /* Arithmetic Mean Filter */ 33 for (int k = 0; k < ARRAY_SIZE*ARRAY_SIZE; k++) 34 out_array[i][j] += value[k]; 35 out_array[i][j] /= ARRAY_SIZE*ARRAY_SIZE; 36 37 } 38 } 39 }
处理高斯噪声
处理椒盐噪声
结论:算术平均滤波对于高斯噪声和椒盐噪声都有一定的效果,但是同时会平滑图像
几何均值滤波器
实现
1 void filtering(short** in_array, short** out_array, long height, long width) 2 { 3 short value[9]; 4 5 for (int i = 0; i < height; i++){ 6 for (int j = 0; j < width; j++){ 7 value[0] = is_in_array(j-1, i-1, height, width) ? in_array[i-1][j-1] : 0; 8 value[1] = is_in_array(j, i-1, height, width) ? in_array[i-1][j] : 0; 9 value[2] = is_in_array(j+1, i-1, height, width) ? in_array[i-1][j+1] : 0; 10 value[3] = is_in_array(j-1, i, height, width) ? in_array[i][j-1] : 0; 11 value[4] = in_array[i][j]; 12 value[5] = is_in_array(j+1, i, height, width) ? in_array[i][j+1] : 0; 13 value[6] = is_in_array(j-1, i+1, height, width) ? in_array[i+1][j-1] : 0; 14 value[7] = is_in_array(j, i+1, height, width) ? in_array[i+1][j] : 0; 15 value[8] = is_in_array(j+1, i+1, height, width) ? in_array[i+1][j+1] : 0; 16 17 /* Geometric Mean Filter */ 18 double product = 1.0; 19 for (int k = 0; k < ARRAY_SIZE*ARRAY_SIZE; k++) 20 product *= value[k]; 21 product = pow(product, 1.0 / 9.0); 22 out_array[i][j] = (short)product; 23 24 if (out_array[i][j] < 0x00) 25 out_array[i][j] = 0x00; 26 else if (out_array[i][j] > 0xff) 27 out_array[i][j] = 0xff; 28 } 29 } 30 }
几何均值滤波器与算术均值滤波器相比,丢失的图像细节更少
谐波均值滤波器
实现
1 void filtering(short** in_array, short** out_array, long height, long width) 2 { 3 short value[9]; 4 5 for (int i = 0; i < height; i++){ 6 for (int j = 0; j < width; j++){ 7 value[0] = is_in_array(j-1, i-1, height, width) ? in_array[i-1][j-1] : 0; 8 value[1] = is_in_array(j, i-1, height, width) ? in_array[i-1][j] : 0; 9 value[2] = is_in_array(j+1, i-1, height, width) ? in_array[i-1][j+1] : 0; 10 value[3] = is_in_array(j-1, i, height, width) ? in_array[i][j-1] : 0; 11 value[4] = in_array[i][j]; 12 value[5] = is_in_array(j+1, i, height, width) ? in_array[i][j+1] : 0; 13 value[6] = is_in_array(j-1, i+1, height, width) ? in_array[i+1][j-1] : 0; 14 value[7] = is_in_array(j, i+1, height, width) ? in_array[i+1][j] : 0; 15 value[8] = is_in_array(j+1, i+1, height, width) ? in_array[i+1][j+1] : 0; 16 17 /* Harmonic Mean Filter */ 18 double sum = 0; 19 for (int k = 0; k < ARRAY_SIZE*ARRAY_SIZE; k++) 20 sum += 1.0 / value[k]; 21 out_array[i][j] = (short)(9.0 / sum); 22 23 if (out_array[i][j] < 0x00) 24 out_array[i][j] = 0x00; 25 else if (out_array[i][j] > 0xff) 26 out_array[i][j] = 0xff; 27 } 28 } 29 }
处理高斯噪声
处理椒盐噪声
对盐粒噪声效果较好,不适用于胡椒噪声,善于处理高斯噪声
逆谐波均值滤波器
Q为滤波器的阶数,Q为正时,消除胡椒噪声,Q为负时消除盐粒噪声
Q=0为算术均值滤波器,Q=-1谐波均值滤波器
实现
1 void filtering(short** in_array, short** out_array, long height, long width) 2 { 3 short value[9]; 4 5 for (int i = 0; i < height; i++){ 6 for (int j = 0; j < width; j++){ 7 value[0] = is_in_array(j-1, i-1, height, width) ? in_array[i-1][j-1] : 0; 8 value[1] = is_in_array(j, i-1, height, width) ? in_array[i-1][j] : 0; 9 value[2] = is_in_array(j+1, i-1, height, width) ? in_array[i-1][j+1] : 0; 10 value[3] = is_in_array(j-1, i, height, width) ? in_array[i][j-1] : 0; 11 value[4] = in_array[i][j]; 12 value[5] = is_in_array(j+1, i, height, width) ? in_array[i][j+1] : 0; 13 value[6] = is_in_array(j-1, i+1, height, width) ? in_array[i+1][j-1] : 0; 14 value[7] = is_in_array(j, i+1, height, width) ? in_array[i+1][j] : 0; 15 value[8] = is_in_array(j+1, i+1, height, width) ? in_array[i+1][j+1] : 0; 16 17 /* Contra-Harmonic Mean Filter */ 18 int Q = 2; 19 double num = 0.0, den = 0.0; 20 for (int k = 0; k < ARRAY_SIZE*ARRAY_SIZE; k++){ 21 num += pow(value[k], Q+1); 22 den += pow(value[k], Q); 23 } 24 out_array[i][j] = (short)(num / den); 25 26 if (out_array[i][j] < 0x00) 27 out_array[i][j] = 0x00; 28 else if (out_array[i][j] > 0xff) 29 out_array[i][j] = 0xff; 30 } 31 } 32 }
Q = 2 消除胡椒噪声
Q = -2消除盐粒噪声
Q = -2消除盐粒噪声后的图像使用Q = 2消除胡椒噪声
再来一次
再来
此时椒盐噪声已经基本消除