MATLAB不使用histeq的直方图均衡与直方图匹配

XJTU机器视觉第一次作业

要求

1 将给定图像进行直方图均衡化处理，画出原图与直方图，以及直方图均衡化后的图像与直方图

2 对给定图像进行直方图匹配，画出原图，匹配图像，原图匹配之后的图像与各自的直方图

一

1. I = imread('001.png');   
2. [height, width] = size(I);  
3. N = numel(I);  % 求图像像素的总数  
4. Pr = imhist(I) / N;  % 原始图像的灰度分布频率  
5. PrA = cumsum(Pr); % 原始图像的累计灰度分布频率  
6. for i = 1:256  
7.     PrA(i,1) = round(PrA(i,1)*255); % 四舍五入  
8. end  
9. outImage = uint8(zeros(height, width));  % 预分配输出图片  
10. for i = 1 : height  
11.    for j = 1 : width  
12.       outImage(i,j) = PrA(I(i,j)); % 使用累积分布频率赋值  
13.    end  
14. end  
15. PrM = imhist(outImage)/N; % 获得新图像的灰度分布频率  
16. m = 0 : 255;   
17. figure;   
18. subplot(221), imshow(I); title('原图');  
19. subplot(222), bar(m, Pr);title('原图灰度分布频率');  
20. subplot(223), imshow(outImage);title('直方图均衡后');  
21. subplot(224),bar(m,PrM);title('直方图均衡后灰度分布频率');  

二

由于原图是rgb图片，故这里我分别对rgb三个通道进行了匹配，实际应用中可以结合图片自行决定

1. clc;  
2. clear;  
3. I1 = imread('original.bmp');  
4. I2 = imread('reference.jpg');  
5. [height1, width1] = size(I1);  
6. [height2, width2] = size(I2);  
7. N1 = numel(I1);  % 求图像像素的总数  
8. Pr1 = imhist(I1) / N1;  % 原始图像的灰度分布频率  
9. PrA1 = cumsum(Pr1); % 原始图像的累计灰度分布频率  
10. N2 = numel(I2);  % 求参考图像像素的总数  
11. N2 = N2/3;  
12. imRef1=I2(:,:,1);  
13. imRef2=I2(:,:,2);  
14. imRef3=I2(:,:,3);  
15. Pr2 = imhist(imRef1) / N2;  % 参考图像的灰度分布频率  
16. Pr3 = imhist(imRef2) / N2;  % 参考图像的灰度分布频率  
17. Pr4 = imhist(imRef3) / N2;  % 参考图像的灰度分布频率  
18. PrA2 = cumsum(Pr2); % 参考图像的累计灰度分布频率  
19. PrA3 = cumsum(Pr3); % 参考图像的累计灰度分布频率  
20. PrA4 = cumsum(Pr4); % 参考图像的累计灰度分布频率  
21. for i = 1:256  
22.     PrA1(i,1) = round(PrA1(i,1)*255); % 四舍五入  
23. end  
24. for i = 1:256  
25.     PrA2(i,1) = round(PrA2(i,1)*255); % 四舍五入  
26. end  
27. for i = 1:256  
28.     PrA3(i,1) = round(PrA3(i,1)*255); % 四舍五入  
29. end  
30. for i = 1:256  
31.     PrA4(i,1) = round(PrA4(i,1)*255); % 四舍五入  
32. end  
33. outImageA = uint8(zeros(height1, width1));  % 预分配输出图片  
34. for i = 1 : height1  
35.    for j = 1 : width1  
36.       outImageA(i,j) = PrA1(I1(i,j)); % 使用累积分布频率赋值  
37.    end  
38. end  
39. % 这里其实可以合在一起做，没必要将统计与变换分开  
40. PrmA1 = imhist(outImageA)/N1;  
41. PrmA2 = zeros(256,1); % 单通道匹配后灰度分布频率  
42. PrmA3 = zeros(256,1);  
43. PrmA4 = zeros(256,1);  
44. M2 = zeros(256,2); % 像素转换数组  
45. M3 = zeros(256,2);  
46. M4 = zeros(256,2);  
47. for i = 1:256  
48.     if PrmA1(i,1) ~= 0  
49.         k = abs(PrA2(1,1)-i);  
50.         h = 1;  
51.         for j = 1:256  
52.             if min(abs(PrA2(j,1)-i))<k  
53.                 k = min(abs(PrA2(j,1)-i));  
54.                 h = j;  
55.             end  
56.         end  
57.         PrmA2(h,1) = PrmA2(h,1) + PrmA1(i,1);  
58.         M2(i,1) = i;  
59.         M2(i,2) = h;  
60.     end  
61. end  
62. for i = 1:256  
63.     if PrmA1(i,1) ~= 0  
64.         k = abs(PrA3(1,1)-i);  
65.         h = 1;  
66.         for j = 1:256  
67.             if min(abs(PrA3(j,1)-i))<k  
68.                 k = min(abs(PrA3(j,1)-i));  
69.                 h = j;  
70.             end  
71.         end  
72.         PrmA3(h,1) = PrmA3(h,1) + PrmA1(i,1);  
73.         M3(i,1) = i;  
74.         M3(i,2) = h;  
75.     end  
76. end  
77. for i = 1:256  
78.     if PrmA1(i,1) ~= 0  
79.         k = abs(PrA4(1,1)-i);  
80.         h = 1;  
81.         for j = 1:256  
82.             if min(abs(PrA4(j,1)-i))<k  
83.                 k = min(abs(PrA4(j,1)-i));  
84.                 h = j;  
85.             end  
86.         end  
87.         PrmA4(h,1) = PrmA4(h,1) + PrmA1(i,1);  
88.         M4(i,1) = i;  
89.         M4(i,2) = h;  
90.     end  
91. end  
92. PrmA2f = cumsum(PrmA2)*255;  
93. PrmA3f = cumsum(PrmA3)*255;  
94. PrmA4f = cumsum(PrmA4)*255;  
95. outImageB = uint8(zeros(height1, width1));  % 预分配输出图片  
96. for i = 1 : height1  
97.    for j = 1 : width1  
98.        outImageB(i,j) = M2(outImageA(i,j)+1,2);  
99.    end  
100. end  
101. outImageC = uint8(zeros(height1, width1));  % 预分配输出图片  
102. for i = 1 : height1  
103.    for j = 1 : width1  
104.        outImageC(i,j) = M3(outImageA(i,j)+1,2);   
105.    end  
106. end  
107. outImageD = uint8(zeros(height1, width1));  % 预分配输出图片  
108. for i = 1 : height1  
109.    for j = 1 : width1  
110.       outImageD(i,j) = M4(outImageA(i,j)+1,2);   
111.    end  
112. end  
113. exp(:,:,1) = outImageB(:,:,1);  
114. exp(:,:,2) = outImageC(:,:,1);  
115. exp(:,:,3) = outImageD(:,:,1);  
116. Pr5 = imhist(exp)/N1;  
117. m = 0:255;  
118. subplot(541), imshow(I1);title('原图');  
119. subplot(542),bar(m,Pr1);title('原灰度分布频率');  
120. subplot(543),imshow(outImageA);title('直方图均衡化原图');  
121. subplot(544),bar(m,PrmA1);title('均衡化灰度分布频率');  
122. subplot(545), imshow(imRef1);title('参考图R通道');  
123. subplot(546),bar(m,PrA2);title('R-accumulation');  
124. subplot(547),bar(m,PrmA2f);title('processed-RA');  
125. subplot(548),bar(m,PrmA2);title('均衡化原图R通道匹配');  
126. subplot(549), imshow(imRef2);title('参考图G通道');  
127. subplot(5,4,10),bar(m,PrA3);title('G-accumulation');  
128. subplot(5,4,11),bar(m,PrmA3f);title('processed-GA');  
129. subplot(5,4,12),bar(m,PrmA3);title('均衡化原图G通道匹配');  
130. subplot(5,4,13), imshow(imRef3);title('参考图B通道');  
131. subplot(5,4,14),bar(m,PrA4);title('B-accumulation');  
132. subplot(5,4,15),bar(m,PrmA4f);title('processed-BA');  
133. subplot(5,4,16),bar(m,PrmA4);title('均衡化原图B通道匹配');  
134. subplot(5,4,17),imshow(I2);title('参考图');  
135. subplot(5,4,18),bar(m,Pr2);title('参考图灰度分布频率');  
136. subplot(5,4,19),imshow(exp);title('直方图匹配后结果');  
137. subplot(5,4,20),bar(m,Pr5);title('匹配后灰度分布频率');  
 

参考结果