matlab直方图拉伸、均衡化和匹配

问题分析

内容 1：图像增强显示

• 编写一个灰度拉伸函数，实现功能：将原始图像的灰度值以线性方式拉伸到 0-255 之间。

• 读取./data/flower.tif 文件，显示原始图像、原始图像的直方图，然后分别采用上述步骤中编写的函数、imadjust 函数以实现图像增强，并显示图像增强后的效果以及相应的图像直方图。

• 读取./data/mars_moon.tif 文件，显示原始图像、原始图像的直方图，阅读imadjust 函数，理解 gamma 参数的意义，然后设置不同的 gamma 值进行图像增强，并显示图像和增强后的图像直方图。

内容 2：直方图均衡化

• 根据直方图均衡化的原理，实现直方图均衡化功能，编写直方图均衡化函数。

• 读取./data/flower.tif 文件，分别利用你编写的直方图均衡化函数和 histeq。

• 函数对其进行直方图均衡化操作，比较显示效果，以及直方图形态。

内容 3：直方图匹配

• 读取./data/bhtmref.img 文件，以./data/bhtmref_reference.tif 文件为参考图像，对所有波段图像进行直方图匹配，并使用 writeenvi 函数输出匹配后的图像为 img 格式。

实验数据

• flower.tif，多波段图像文件；

• mars_moon.tif，灰度图像；

• bhtmref.img，遥感图像；

• bhtmref_reference.tif ，遥感图像直方图匹配的参考。

操作步骤

直方图拉伸

1 灰度拉伸函数

%%灰度拉伸函数
function  Img2 = enhance( Img )
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
    low=min(min(Img));
    hi=max(max(Img));
    cha=hi-low;
    Img2=(Img-low).*(255/cha);
    return 
end

2 图像拉伸

%% 图像拉伸
% 读取./data/flower.tif 文件，显示原始图像、原始图像的直方图，然后分别
% 采用上述步骤中编写的函数、imadjust 函数以实现图像增强，并显示图像增强后
% 的效果以及相应的图像直方图。

img=imread('flower.tif');
img2=enhance(img);
img3=imadjust(img);

figure(1)
subplot(3,2,1)
imshow(img)
title('原图')
subplot(3,2,2)
imhist(img,64)
subplot(3,2,3)
imshow(img2)
subplot(3,2,4)
imhist(img2,64)
subplot(3,2,5)
imshow(img3)
subplot(3,2,6)
imhist(img3,64)

4 GAMMA参数的调整

%% gamma参数的调整
% 读取./data/mars_moon.tif 文件，显示原始图像、原始图像的直方图，阅读
% imadjust 函数，理解 gamma 参数的意义，然后设置不同的 gama 值进行图像增强，
% 并显示图像和增强后的图像直方图

% GAMMA指定描述关系的曲线的形状
% 在I和J中的值之间。如果GAMMA小于1，则映射为
% 向更高（更亮）的输出值加权。如果GAMMA大于
% 如果大于1，则映射将向较低（较暗）的输出值进行加权。如果
%如果省略参数，GAMMA默认为1（线性映射）。

img4=imread('mars_moon.tif');
img5=imadjust(img4,[],[],0.2);
img6=imadjust(img4,[],[],5);

figure(2)
subplot(3,2,1)
title('原图')
imshow(img4)
subplot(3,2,2)
imhist(img,64)
subplot(3,2,3)
imshow(img5)
subplot(3,2,4)
imhist(img2,64)
subplot(3,2,5)
imshow(img6)
subplot(3,2,6)
imhist(img3,64)

直方图均衡化

1 均衡化函数

function  outImage  = junheng( 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

% 接下来，将频数值算为频率
ProbPixel = zeros(1,256); % 统计各灰度级出现的频率
for i = 1 : 256
    ProbPixel(i) = NumPixel(i) / (height * width);
end


% 再用函数cumsum()来计算累积分布函数（CDF），并将频率（取值范围是0~1）映射到0~255的无符号整数
CumPixel = cumsum(ProbPixel);  % 这里的数组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));
   end
end

end

2 均衡化应用

%% 直方图均衡化

img=imread('flower.tif');
img2=junheng(img);
img3=histeq(img);

figure(1)
subplot(3,2,1)
imshow(img)
title('原图')
subplot(3,2,2)
imhist(img,64)

subplot(3,2,3)
imshow(img2)
title('均衡化')
subplot(3,2,4)
imhist(img2,64)
subplot(3,2,5)
imshow(img3)
title('均衡化')
subplot(3,2,6)
imhist(img3,64)

直方图匹配

%% 直方图匹配
% 读取./data/bhtmref.img 文件，以./data/bhtmref_reference.tif 文件为参考图
% 像，对所有波段图像进行直方图匹配，并使用 fwriteenvi 函数输出匹配后的图像
% 为 img 格式
ref=imread('bhtmref_reference.tif');
img=freadenvi('bhtmref.img');
HistRef=imhist(ref);
img=uint8(img);
[m,n,p]=size(img);
newImg=[];
figure;
imshow(ref)
for i=1:p
    matched = histeq(img(:,:,i),HistRef);
    newImg=cat(3,newImg,matched);
end
figure(2);
for i=1:p
    subplot(2,3,i)
    imshow(newImg(:,:,i))
end
fwriteenvi('bhtmref_new.img',newImg);