MATLAB编写图像直方图均衡化函数

clc;
% 首先读入灰度图像，并提取图像的高度和宽度
A = imread('lena.jpg'); 
subplot(2,2,1);
imshow(A);
image=rgb2gray(A);%设置为灰度图像
subplot(2,2,2);
imhist(uint8(image));%图像直方图
[height, width] = size(image); 
% 然后统计每个灰度的像素值的累计数目
NumPixel = zeros(1,256);  % 建立一个256列的行向量，以统计各灰度级的像素个数
for i = 1 : height
   for j = 1 : width
       k = image(i,j);  % k是像素点(i,j)的灰度值
       % 因为NumPixel数组的下标是从1开始的，但是图像像素的取值范围是0~255
       % 所以用NumPixel(k+1)
       NumPixel(k+1) = NumPixel(k+1) + 1;  % 对应灰度值像素点数量加1 
   end
end
% 将频数值算为频率
pinlv = zeros(1,256); % 统计各灰度级出现的频率
for i = 1 : 256
    pinlv(i) = NumPixel(i) / (height * width);
end
 
% 再用函数cumsum()来计算累积分布函数（CDF），并将频率（取值范围是0~1）映射到0~255的无符号整数
CumPixel = cumsum(pinlv);  % 这里的数组CumPixel大小也是1×256
CumPixel = uint8((256-1) .* CumPixel + 0.5);  
% 在下列用作直方图均衡化实现的赋值语句右端，image(i,j)被用来作为CumPixel的索引
% 例如，image(i,j)=120，则从CumPixel中取出第120个值作为image(i,j)的新像素值
outImage = uint8(zeros(height, width));  % 预分配数组
for i = 1 : height
   for j = 1 : width
      outImage(i,j) = CumPixel(image(i,j)+1);%因为MATLAB数组下标从1开始，而灰度级从0开始
   end
end
subplot(2,2,3);
imshow(outImage);
subplot(2,2,4);
imhist(uint8(outImage));%图像直方图