图像的直方图均衡化和比特平面分层

（1）自定义一个函数，当输入为一幅图像(EXP3-1.tif)时，能输出该图像的直方图。

计算输入图像的直方图 getHist()

function H=getHist(pho)
ima=imread(pho);
ima=double(ima);
[r c]=size(ima);
H=zeros(1,256); % H(i)表示在灰度为i-1时图像中该灰度出现的次数
for k=1:256
    tempK=k-1; %计算灰度
    kPos=find(ima(:)==tempK); %寻找灰度为 tempK 所在图像的位置
    H(k)=length(kPos); %求得所有灰度为tempK的个数
end
H=H/(r*c); %对直方图进行归一化

按要求输入实验所给图像，输出原图像和对应的直方图：

function []=showSR()
H=getHist('EXP3-1.tif');
ima=imread('EXP3-1.tif');
subplot(121);
imshow(ima);
title('输入图像');
subplot(122);
bar(H,'b');
title('输入图像的直方图');

运行结果如下：

（2）自定义一个函数，将输入图像进行直方图均衡化处理，该函数输入参数为处理前的图像，输出参数为处理后的图像和直方图均衡化过程的变换函数。

计算图像直方图均衡化后的图像以及直方图均衡化过程的变换函数：

% im2表示均衡化的图像，T表示均衡化的变换函数
function [im2,T]=JunHengHua(pho)
im1=imread(pho);
[r c]=size(im1);
im1=double(im1(:)); %将其转化为行向量，易于下面对灰度级的查找
im2=zeros(r*c,1); %对应转化后的im1
L=256;
im1_H=zeros(1,L); %原图像的直方图
%im2_H=im1_H;
H_sum=0; %累计原图像的直方图
T=im1_H;
for k=1:L
    tempK=k-1;
    posK=find(im1==tempK); %寻找对应灰度值在原图像中的位置
    im1_H(k)=length(posK)/(r*c); %原图像的直方图，并进行归一化
    H_sum=H_sum+im1_H(k); %以求和代替积分，求累计原图像的直方图
    T(k)=(L-1)*H_sum; %求得原图像均衡化过程的变化函数
    if (T(k)-floor(T(k))) > 0.5 %进行四舍五入操作
        im2(posK)=ceil(T(k)); %得到变换后的灰度，赋值给原位置
    else
        im2(posK)=floor(T(k));
    end
end
im2=uint8(reshape(im2,r,c)); %将im2转为r行c列 无符号8位的矩阵

按要求输入实验所给图像，输出原图像、对应的直方图均衡化的图像以及其变换函数：

function []=showJunHengHua()
ima=imread('EXP3-1.tif');
subplot(131);
imshow(ima);
title('输入图像');
 
[im2 T]=JunHengHua('EXP3-1.tif');
subplot(132);
imshow(im2);
title('直方图均衡化的图像');
 
subplot(133);
plot(T);
title('直方图均衡化过程的变换函数');

运行结果如下：

（3）求图像的比特平面（该方法参考网上资料）

function res=get8BitPM(pho)
ima=imread(pho);
[r c]=size(ima);
n=8;
res=zeros(r,c,n); %各个比特平面，其中res(:,:,1)表示一个比特平面等
 
for i=1:n
    for j=1:r
        for k=1:c
            gray=bitand(ima(j,k),2^(i-1));  %进行按位与运算
            if gray==2^(i-1)
                res(j,k,i)=255;
            else
                res(j,k,i)=0;
            end
        end
    end
end

利用所得的比特平面，输出原图像和比特平面：

function []=showBitPM()
H=getHist('EXP3-2.tif');
ima=imread('EXP3-2.tif');
subplot(3,3,1);
imshow(ima);
title('输入图像');
 
%获取各比特平面
bitPM=get8BitPM('EXP3-2.tif');
 
%输出各比特平面
num=['1' '2' '3' '4' '5' '6' '7' '8']; %各比特平面的标号
for i=1:8
    subplot(3,3,i+1);
    imshow(uint8(bitPM(:,:,i)));
    title(['第',num(i),'个比特平面']);%设置该比特平面的标题
end

运行结果如下：

这里提供了一种我自己写的方法：

%图像的比特平面分层
function bitPM = get8BitPM2(pho)
ima=imread(pho);
subplot(331);
imshow(ima);
title('输入图像');
[r c] = size(ima);
im1 = dec2bin(ima(:),8); %将原图像的矩阵转化为列向量，并且将向量的各个元素转化为8位二进制数
A=zeros(r*c,1); %用于临时记录一个比特平面
bitPM=zeros(r,c,8); %用于存储所有的比特平面。如 bitPM(:,:,1)表示第1比特平面
 
% 8比特灰度图像，共有8个1比特平面
for n=1:8
    A=bin2dec(im1(:,8-n+1))*255; %其中 im1(:,8-n+1)表示取 im1中 所有行的元素的第 8-n+1 位，再将其转化为 十进制数。此时该数为0或1，将灰度进行拉伸，变换到0或255
    %即得第 8-n+1比特平面
    
    A=reshape(A,r,c); %将A转化为r行 c列的矩阵
    bitPM(:,:,8-n+1)=A; %保存该比特平面
    subplot(3,3,n+1);
    imshow(uint8(A)); %输出 第 8-n+1 比特平面
    title(['比特平面',num2str(n)]);
end

get8BitPM2(‘EXP3-2.tif’);

运行结果如下：