哈工大深研院数字图像处理第一次大作业：不调用Matlab现有库函数实现图像增强

题目是：

Image Enhancement
1. Histogram Equalization: Write a program for computing the histogram of the image (a) and implement the histogram equalization technique as we discussed in the class.

2. Remove the noise from the image (c). 

3. Filtering in the Frequency Domain: Implement the Gaussian highpass filters on test image (e).

Attention：
1. Don’t use any commercial tools such as matlab functions, Open CV or other library (except FFT functions). Please do it by c/c++, java or other languages.
2. Source codes and projects are required when you submit your report. Submitting names will be ‘Pro1_studentnumber_name.zip/rar’.
3. Each student should independently finish this work. 
4. Deadline: 2016-04-20.
5. All test images can be downloaded from the QQ group.

而本宝宝的报告是：

汇报人：张常乐 学号：15S158746 班级：电子与通信工程
程序清单：
1.Histogram_Equalization.m
2.Noise_Remove.m
3.Gaussian_Highpass_Filtering.m
程序说明：
1.依据灰度直方图均衡化处理原理，我首先得到了原始图像的灰度直方图Fig.1(a)，从图像中我们可以看到，图像灰度集中分布于灰度值较低的区域，表现在图像上Fig.2我们可以看出整幅图像颜色昏暗，对比度不高，为了解决这个问题，我们采用图像直方图均衡化的处理方法，分别计算出每个灰度出现的频率和概率，并利用离散积分函数：S_k=(L-1)/MN ∑_(j=0)^(k-1)▒n_j ,k=0,1,2,…,L-1(1)
计算出灰度的累积分布情况。利用这个分布，构造出原始灰度向均衡化灰度转换的映射关系，并利用这种映射关系来对灰度图像进行均衡化处理。最终可以得到灰度直方图分布如Fig.1(b)所示的灰度分布以及Fig.3所示的处理后结果。

2.针对Fig.4中存在的椒盐噪声，采用中值滤波可以达到很好的滤除效果。通过构建
1 1 1
1 1 1
1 1 1
这样的模板，将模板内的中值替代为模板中心的图像值，使这个模板扫过整幅图像便可以实现噪声滤除效果。同时，由于边界效应，我们将图像的边界进行了扩充。实践证明图像边界的扩充可以在处理后通过恰当的处理结果裁剪而消除影响，并且这种扩充对于图像来说将不会造成影响。这里有几点需要注意：
1)模板的维度必须是奇数维正方形，否则将会在模板中心选取时造成很大的困扰；
2)模板大小应当与需要处理的椒盐噪声的尺寸有关，同时，越大的模板将会对图像处理结果造成越模糊的影响。通过测试发现选用3*3的模板恰好可以满足本实验的要求；
3)针对中值选取，用于不允许使用Matlab自带的median以及sort函数，我采用了冒泡排序的方法，因为Matlab软件本身对于for循环处理的速度较慢，所以整个图像的处理速度显得较慢。

3.针对题目3中要求的的图像锐化，可以采用高斯高通滤波器进行锐化处理。通过构造图像Fig.6所对应的归一化频率分布，构造空间频率分布，利用高斯高通滤波器表达式(GHPF)：
H(u,v)=1-e^(-(D^2 (u,v))/(2〖D_0〗^2 ))           (2)
构造传递函数H(u,v)，对原始图像x(m,n)进行频谱对称化(fftshift)后的快速傅里叶变换(FFT)，再与传递函数相乘得到处理后的结果Y(u,v)，最终对该结果进行逆傅里叶变换，得到我们需要的锐化图像y(m,n)。如图Fig.7所示。

程序代码也给大家贴上来好了

1.

[plain]  view plain  copy

 

1. %Name:Histogram_Equalization  
2. %function:对直方图进行均衡化处理实现图像增强  
3. %Author:Changle Zhang  
4. %Student ID:15S158746  
5. clear all;  
6. close all;  
7. clc;                                  %Initialization  
8.   
9. fig = imread('a.jpg');                %import the image  
10. subplot(121)  
11. plot(fig);  
12. title('(a)Gray distribution of image befort process');  
13.   
14. [row,col] = size(fig);                %Image size  
15. PixelsNum = row * col;                %Number of Pixels  
16. freq = zeros(256,1);                  %record frequency of each gray level  
17. prob = zeros(256,1);                  %record probability of each gray level  
18. %calculate frequency and probability  
19. for i = 1:row  
20.     for j= 1:col  
21.         value = fig(i,j);  
22.         freq(value+1) = freq(value+1) + 1;  
23.         prob(value+1) = freq(value+1) / PixelsNum;  
24.     end  
25. end  
26.   
27. cums = zeros(256,1);                  %record the  cumulative distribution  
28. probc = zeros(256,1);                 %record the probability of each distribution  
29. output = zeros(256,1);                %gray level after H_E  
30. counter = 0;  
31. graylevel = 255;  
32. %calculate cumulative distribution and output gray level  
33. for i = 1:size(prob)  
34.     counter = counter + freq(i);  
35.     cums(i) = counter;  
36.     probc(i) = cums(i) / PixelsNum;  
37.     output(i) = round(probc(i) * graylevel);  
38. end  
39. %HE process  
40. outputimage = uint8(zeros(row,col));  %Final image  
41. for i = 1:row  
42.     for j = 1:col  
43.         outputimage(i,j) = output(fig(i,j)+1);  
44.     end  
45. end  
46. %outputing  
47. subplot(122)  
48. plot(outputimage);  
49. title('(b)Gray distribution of image after process');  
50. figure  
51. imshow(fig);  
52. title('Image before histogram equalization');  
53. figure  
54. imshow(outputimage);  
55. title('Image after histogram equalization');  

2.

[plain]  view plain  copy

 

1. %Name:Noise_Remove  
2. %function:利用3*3中值滤波模板去除含噪声图像中的椒盐噪声  
3. %Author:Changle Zhang  
4. %Student ID:15S158746  
5. clear all;  
6. close all;  
7. clc;                                                    %Initialization  
8.   
9. img=imread('c.jpg');                                    %Import image  
10. imshow(img);  
11. title('Image With Salt&pipe Noise');  
12. [m n]=size(img);  
13.   
14. Nmax=3;                                                 %Max filtering Radius(Nmax must be odd)  
15.   
16. imgn=zeros(m+2*Nmax+1,n+2*Nmax+1);  
17. imgn(Nmax+1:m+Nmax,Nmax+1:n+Nmax)=img;  
18.   
19. imgn(1:Nmax,Nmax+1:n+Nmax)=img(1:Nmax,1:n);                                     %Upper Expansion  
20. imgn(1:m+Nmax,n+Nmax+1:n+2*Nmax+1)=imgn(1:m+Nmax,n:n+Nmax);                     %Right Expansion  
21. imgn(m+Nmax+1:m+2*Nmax+1,Nmax+1:n+2*Nmax+1)=imgn(m:m+Nmax,Nmax+1:n+2*Nmax+1);   %Lower Expansion  
22. imgn(1:m+2*Nmax+1,1:Nmax)=imgn(1:m+2*Nmax+1,Nmax+1:2*Nmax);                     %Left Expansion  
23.   
24.   
25. [height, width]=size(imgn);                             %Size of Expansived Image  
26. x1=double(imgn);  
27. x2=x1;  
28. for i=1:height-Nmax+1  
29.     for j=1:height-Nmax+1  
30.         c=x1(i:i+(Nmax-1),j:j+(Nmax-1));                %get the initial filtering point  
31.         e=c(1,:);                                       %first row of matrix 'c'  
32.         for u=2:Nmax  
33.             e=[e,c(u,:)];                               %make 'c' a row vector 'e'     
34.         end  
35.         n1=Nmax*Nmax;                                   %Number of pixel in model  
36.         for l=1:n1-1                                    %Sort the row vector  
37.             for q=1:n1-1  
38.                 if e(q)>e(q+1)  
39.                     c1 = e(q);  
40.                     e(q) = e(q+1);  
41.                     e(q+1) = c1;  
42.                 end  
43.             end  
44.          end  
45.         mm = e((n1+1)/2);                               %mm is the median number  
46.         x2(i+(Nmax-1)/2,j+(Nmax-1)/2)=mm;               %assignment median value for the center point  
47.     end  
48. end   
49. %keep the unassignment points  
50. d=uint8(x2);  
51. figure;  
52. imshow(d(Nmax+1:m+Nmax,Nmax+1:n+Nmax),[]);  
53. title('Image After Noise Removement');  

3.

[plain]  view plain  copy

 

1. %Name:Gaussian_Highpass_Filtering  
2. %function:利用高斯高通滤波器锐化图像  
3. %Author:Changle Zhang  
4. %Student ID:15S158746  
5. clear all;  
6. close all;  
7. clc;                                                    %Initialization  
8.   
9. img=imread('e.tif');                                    %Import image  
10. imshow(img);  
11. title('Image Before Sharpen');  
12. [m n]=size(img);                                        %Get the image size m*n  
13. %Generate normalization uniformly-spaced point in frequency domain  
14. for i= 1:m  
15.     if mod(n,2)==0                                      %n is even  
16.         f1(i,:)=[-n+1:2:n-1]./(n-1);  
17.         f2(:,i)=f1(i,:)';  
18.     else                                                %n is odd  
19.         f1(i,:)=[-n:2:n-2]./(n-1);  
20.         f2(:,i)=f1(i,:)';  
21.     end  
22. end  
23.   
24. D = 0.2;                                                %Gaussian end frequency   
25. r = f1.^2+f2.^2;                                        %Filerting radius  
26. %Generate filering transformation H  
27. for i =1:m  
28.     for j = 1:n  
29.         t=r(i,j)/(D*D);  
30.         Hd(i,j)=1-exp(-t);  
31.     end  
32. end  
33. X=fftshift(fft2(double(img)));                          %FFT  
34. Y=X.*Hd;                                                %Apply high-pass filtering  
35. output=real(ifft2(ifftshift(Y)));                       %Ifft and get the real part  
36. figure;  
37. imshow(output,[]);  
38. title('Image After High-Pass filtering');