影像增强处理

影像增强处理(matlab)

本文将介绍遥感数字图像处理中的影像增强处理部分。将按照以下9个方面内容的基本原理和代码实现进行介绍。
分别为直方图线性拉伸、直方图均衡化、直方图匹配、均值滤波、中值滤波、Sobel算子锐化、Prewitt算子锐化、RGB->HIS变换&&HIS->RGB变换、PCA(K-L)变换。
所有程序均用matlab实现。

直方图线性拉伸

基本原理：
对直方图进行线性拉伸是为了改善图像的对比度，采用线性或分段线性的函数关系称为线性变换。线性变换公式如下所示：

首先找出图像的最大、最小灰度值，然后依据上述公式进行线性变换
代码实现如下所示：

function image_=xxls(image,bmax,amax)
[rows,cols,bands]=size(image);
for band=1:bands
	oneband=image(:,:,band);
	oneband=double(oneband(:));
	%寻找该波段灰度值的最大值和最小值
	images=image;
	amax=oneband(1);
	amin=oneband(1);
	for i=2:rows*cols
		if(oneband(i)>max)
			max=oneband(i);
		end
		if(oneband(i)<min)
			min=oneband(i);
		end
	end
	%根据公式做线性变换
	oneband=(bmax-bmin)*(oneband-amin)/(amax-amin)+bmin;
	%通过reshape恢复图像
	images(:,:,band)=reshape(oneband,rows,cols);
end
image_=images;
end

直方图均衡化

1.基本原理：是将一已知灰度概率密度分布的图像，经过某种变换，变成一幅具有均匀灰度概率密度分布的新图像，其结果是扩展了像元取值的动态范围，从而达到增强图像整体对比度的效果。
2.步骤：
1>计算图像的累计直方图；
2>计算原每一灰度级均衡化后的新值（取整）；
3>以新值代替原灰度值，形成均衡化后的新图像；
3.计算方法

4.程序实现

function image_=equalization(image)
[rows,cols,bands]=size(image);
for band=1:bands
	oneband=image(:,:,band);
	oneband=double(oneband(:));
	%调用static函数得到累计直方图
	myst=static(image);
	myljzft=myst.ljzft;
	%寻找该波段图像灰度最大值
	max=oneband(1);
	for i=2:rows*cols
		if(oneband(i)>max)
			max=oneband(i);
		end
	end
	%计算均衡化后的新值
	b=ljzft(oneband+1)/rows/cols*max;
	%取整
	b=round(b);
	%以新值代替原来的旧值
	oneband=b;
	%调用reshape恢复图像
	image_(:,:,band)=reshape(oneband,rows,cols);
end
end

直方图匹配

1.基本原理：直方图匹配是指使一副图像的直方图变成规定形状的直方图而进行的图像增强的方法。通过变换，是两幅图像的亮度变化规律尽可能地接近。
2.步骤：
1>计算原图像和参考图像的累计直方图
2>将两个累计直方图标准化
3>计算匹配后的新灰度值
4>新值代替旧值
3.代码实现

function image_=matching(rimage,oimage)
%rimage为参考图像、oimage为原图像
[orows,ocols,obands]=size(oimage);
[rrows,rcols,rbands]=size(rimage);
for band=1:obands
	oneoband=oimage(:,:,band);
	oneoband=double(oneoband(:));
	onerband=rimage(:,:,band);
	onerband=double(onerband(:));
	%调用static函数计算两个图像的累计直方图
	myst1=static(oimage);
	myst2=static(rimage);
	oljzft=myst1.ljzft;
	rljzft=myst2.ljzft;
	%将两个图像的累计直方图标准化
	oljzft(band+1,:)=oljzft(band+1,:)/orows/ocols;
	rljzft(band+1,:)=rljzft(band+1,:)/rrows/rcols;
	%初始化匹配后的灰度值
	mat=zeros(1,256);
	%计算距离最小值
	for i=1:256
		mindist=1;
		index=-1;
		for j=1:256
			if(abs(oneoband(i)-onerband(j)<mindist)
				mindist=abs(oneoband(i)-onerband(j);
				index=j;
			end
		end
		mat(i)=index-1;
	end
	%新值代替旧值
	oneoband=mat(oneoband+1);
	oimage(:,:,band)=reshape(oneoband,orows,ocols);
end
image_=oimage;
end

均值滤波

基本原理：将每个像元在以其为中心的区域内，取一定像元的平均值来代替该像元的值，从而平滑图像细节，降低图像反差。该算法简单，计算速度块，再去掉尖锐噪声的同时，造成图像模糊，对特别是图像的边缘和细节消弱较多，且随着邻域范围的扩大，去噪能力增强的同时，模糊程度越严重。
代码实现：

function image_=MeanLB(n,image)
%n为模板大小
[rows,cols]=size(image);
%初始化输出图像
image_=image;
%根据n计算模板中间到两边的距离
dist=(n-1)/2;
%边缘像元不做处理
for i=2:rows-1
	for j=2:cols-1
		%提取模板窗口
		board=image(i-dist:i+dist,j-dist:j+dist);
		%转换为一维向量并转换为double类型
		board=double(board(:));
		%调用mean函数计算均值
		average=mean(board);
		%把计算的值赋给中间像元
		image_(i,j)=average;
	end
end
end

中值滤波

基本原理：中值滤波的基本原理和均值滤波类似，只是中间像元的灰度值不是均值，而是中值。
代码实现：

function image_=Median(n,image)
%n为模板大小
[rows,cols]=size(image);
%初始化输出图像
image_=image;
%根据n计算模板中间到两边的距离
dist=(n-1)/2;
%边缘像元不做处理
for i=2:rows-1
	for j=2:cols-1
		%提取模板窗口
		board=image(i-dist:i+dist,j-dist:j+dist);
		%转换为一维向量并转换为double类型
		board=double(board(:));
		%调用sort函数排序
		board=mean(board);
		%求中值
		mid=board(floot(n*n/2));
		%把计算的值赋给中间像元
		image_(i,j)=mid;
	end
end
end

Sobel算子锐化

基本原理：图像的边缘部分不做处理，利用Sobel算子计算水平、垂直分量，并将二者相加得到新的灰度值。
代码实现：

function image_ = Sobel(image)
%Sobel算子锐化
%image为灰度图

%用size函数获取图像的行列数
[rows,cols]=size(image);
%初始化image_图像，其为一张灰度为0的图像
image_=zeros(rows,cols);
%对原始图像边缘部分不做处理，设置灰度值为0
for row=2:rows-1
    for col=2:cols-1
        %提取模板的大小
        board=image(row-1:row+1,col-1:col+1);
        %将二维矩阵转换成double类型
        board=double(board);
        %根据Sobel算子计算水平、垂直分量
        x=board(3,1)+2*board(3,2)+board(3,3)-board(1,1)-board(1,3)-2*board(1,2);
        y=board(1,3)+2*board(2,3)+board(3,3)-board(1,1)-board(3,1)-2*board(3,1);
        %把水平、垂直分量的绝对值相加
        ima=abs(x)+abs(y);
        image_(row,col)=ima;
    end
end
end

Prewitt算子锐化

基本原理：Prewitt算子锐化的基本原理和Sobel算子锐化的基本原理类似，只是算子不同。
代码实现：

function image_= Prewitt(image)
%Prewitt算子锐化
%image为原始图像，image_是经锐化后的图像

%用size函数获取原图像的行列数
[rows,cols]=size(image);
%初始化image_,灰度值全为0
image_=zero(rows,cols);
%边缘不做处理，设置灰度值为0
for row=2:rows-1
    for col=2:cols-1
        %提取模板大小
        board=image(row-1:row+1,col-1:col+1);
        %转换为double类型
        board=double(board);
        %根据prewitt算子计算水平、垂直分量
        x=board(3,1)+board(3,2)+board(3,3)-board(1,1)-board(1,2)-board(1,3);
        y=board(1,3)+board(2,3)+board(3,3)-board(1,1)-board(2,1)-board(3,1);
        ima=abs(x)+abs(y);
        image_(row,col)=ima;
    end
end
end

RGB->HIS变换&&HIS->RGB变换

基本原理：RGB彩色模型到HSI彩色模型的转换

HSI彩色模型到RGB彩色模型的转换



程序实现如下：

function image_ = RGB_HIS(image,k)
%RGB和HIS的相互变换
%image为原始图像,image_为经过变换后的输出图像
% 若k=0则执行从RGB到HIS的变换，若k=1则执行从HIS到RGB的变换

switch k
    case 0
        %将图像归一化
        data=single(image)/255;
        %取出行列数
        [rows,cols]=size(data(:,:,1));
        %初始化输出图像
        image_=zeros(rows,cols,3);
        for row=1:rows
            for col=1:cols
                %取出红黄蓝波段各值
                R=data(row,col,1);
                G=data(row,col,2);
                B=data(row,col,3);
                m=min(R,G,B);
                theta=acos(((R-G)+(R-B))/(2*sqrt((R-G)*(R-G)+(R-B)*(G-B))));
                I=(R+G+B)/3;
                S=1-3*m/(R+G+B);
                if (G>=B)
                    H=theta;
                end
                if (G<B)
                    H=360-theta;
                end
                %将三个波段组合成一个图像
                image_(row,col,1)=I;
                image_(row,col,2)=H;
                image_(row,col,3)=S;
            end
        end
    case 1
        [rows,cols]=size(image);
        image_=zeros(rows,cols,3);
        for row=1:rows
            for col=1:cols
                %取出H,I,S
                I=image(row,col,1);
                H=image(row,col,2);
                S=image(row,col,3);
                %判断H处于那个区间
                if (H>=0&&H<120)
                    B=1/(1-S);
                    R=1/(1+S*cos(H)/(cos(60-H)));
                    G=-3/(B+R);
                end
                if (H>=120&&H<240)
                    R=I*(1-S);
                    G=1/(1+(S*cos(H-120)/(cos(180-H))));
                    B=-3/(R+G);
                end
                if (H>=240&&H<360)
                    G=1/(1-S);
                    B=1/(1+(S*cos(H-240))/cos(300-H));
                    R=-3/(G+B);
                end
                image_(row,col,1)=R;
                image_(row,col,2)=G;
                image_(row,col,3)=B;
            end
        end
end
end

PCA变换（K-L变换）

代码实现：

function image_=PCA(image,n)
%n为最后图层数目
[rows,cols,bands]=size(image);
%将每个波段转换为一个列向量，且为double型
imMatrix=double(reshape(image,[],bands));
%调用static函数求协方差矩阵
myst=static(image);
var=myst.var;
%初始化输出图像
image_=(rows,cols,n);
%调用eig函数求取特征函数和特征向量
(X,B)=eig(var);
X=X';
for iband=1:n
	t=zeros(rows,cols,n);
	for band=1:bands
		t=t+imMatrix(:,band)*X(bands+1-iband,band);
	end
end
%有可能出现负数，取绝对值
image_=uint8(abs(reshape(t,rows,cols);
end