数字图像处理程序合集中(matlab实现)


前言

接上篇文章,这是这个合集的中集,主要包含影像的增强处理影像的融合处理两大部分内容,请大家理性食用,欢迎大家批评指正!

其中用到了很多上篇中的函数,在用到的地方我会点出!

1.影像增强处理

1.1 直方图线性拉伸

function B=xxls(A) % 线性拉伸的函数
[row,col,band]=size(A);
samples=row*col;
A1=reshape(A,[samples,band]);
D=A1;
maxA=max(A1);
minA=min(A1);
for i=1:band
    for j=1:samples
        D(j,i)=255*(A1(j,i)-minA(i))/(maxA(i)-minA(i))+0;
    end
end
B=reshape(D,[row,col,band]);
figure(1)
subplot(1,2,1)
imshow(uint8(A))
% histogram(A)
title('原始图像')
subplot(1,2,2)
imshow(uint8(B))
% histogram(B)
title('线性拉伸后的图像')

1.2 直方图均衡化

function B=zftjh(A) % 直方图均衡函数
[row,col,band]=size(A);
samples=row*col;
A1=reshape(A,[samples,band]);
A2=cal_hist(A);
A3=zeros(samples,band);
for i=1:band
    for j=1:samples
        A3(j,i)=round(255*A2(A1(j,i)+1,2,i)); % 注意要向下取整,离散的灰度值只能为整数
    end
end
B=reshape(A3,[row,col,band]);
figure(1)
subplot(1,2,1)
imshow(uint8(A))
title('原始图像')
subplot(1,2,2)
imshow(uint8(B))
title('直方图均衡化后的图像')
end

用到了上篇中计算直方图函数cal-hist()


1.3 直方图匹配

function C=zftpp(A,B) % A是待匹配图像,B是参考图像,C是输出的直方图匹配后的图像
[row,col,band]=size(A);
C=A;
hist1=cal_hist(A);
hist2=cal_hist(B);
P=zeros(256,band);
for k=1:band
   for i=1:256 % 256是hist矩阵的行数
     [~,index] = min(abs(hist1(i,2,k)-hist2(:,2,k)));
     P(i,k)=index-1;
   end
end
for i=1:band
    for j=1:row
        for k=1:col
            C(j,k,i)=P(A(j,k,i)+1,i);
        end
    end
end
figure(1);
subplot(1,3,1);
imshow(uint8(A));
title('原始图像');
subplot(1,3,2);
imshow(uint8(B));
title('参照图像');
subplot(1,3,3);
imshow(uint8(C));
title('直方图匹配化后的图像')
end

用到了上篇中计算直方图函数cal-hist()


1.4 均值滤波

function B=jzlb(A) % 计算均值滤波的函数
[m,n,band]=size(A);
format long
B=A;
mode=[0.111,0.111,0.111;0.111,0.111,0.111;0.111,0.111,0.111];
for i=1:band
    for j=2:m-1
        for k=2:n-1
            linyu=[A(j-1,k-1,i),A(j,k-1,i),A(j+1,k-1,i);A(j-1,k,i),A(j,k,i),A(j+1,k,i);A(j-1,k+1,i),A(j,k+1,i),A(j+1,k+1,i)];
            B(j,k,i)=sum(dot(linyu,mode),2); % dot函数用来求内积,x1*y1+x2*y2+xn*yn
        end
    end
end
figure(1)
subplot(1,2,1)
% imshow(uint8(A))
histogram(A)
title('原始图像')
subplot(1,2,2)
% imshow(uint8(B))
histogram(B)
title('均值滤波后的图像')
end

1.5 中值滤波

function B=zzlb(A,n1)% 中值滤波
[m,n,band]=size(A);
format long % 设置输出格式
B=zeros(m,n,band);
mode=[0.111,0.111,0.111;0.111,0.111,0.111;0.111,0.111,0.111];%不知道怎么输入1/9,所以用0.111来替代了
for i=1:band
    for j=2:m-1
        for k=2:n-1
             linyu=[A(j-1,k-1,i),A(j,k-1,i),A(j+1,k-1,i);A(j-1,k,i),A(j,k,i),A(j+1,k,i);A(j-1,k+1,i),A(j,k+1,i),A(j+1,k+1,i)];
             sunzi=linyu.*mode;
             sunzi=sort(reshape(sunzi,[1,n1]),2);
             B(j,k,i)=sunzi(round(n1/2)); % round函数是向上取整
        end
    end
end
figure(1)
subplot(1,2,1)
imshow(uint8(A))
title('原始图像')
subplot(1,2,2)
imshow(uint8(B))
title('中值滤波后的图像')
end

1.6 sobel算子

function B=sobel(A) %sobel锐化算子
[m,n,band]=size(A);
format long
mode1=[-1,-2,-1;0,0,0;1,2,1];
mode2=[-1,0,1;-2,0,2;-1,0,1];
B=A;
for i=1:band
    for j=2:m-1
        for k=2:n-1
            linyu=[A(j-1,k-1,i),A(j,k-1,i),A(j+1,k-1,i);A(j-1,k,i),A(j,k,i),A(j+1,k,i);A(j-1,k+1,i),A(j,k+1,i),A(j+1,k+1,i)];
            B(j,k,i)=abs(sum(sum(linyu.*mode1),2))+abs(sum(sum(linyu.*mode2),2));
        end
    end
end
figure(1)
subplot(1,2,1)
imshow(uint8(A))
title('原始图像')
subplot(1,2,2)
imshow(uint8(B))
title('sobel算子锐化后的图像')
end

1.7 prewitt算子

function B=prewitt(A) % prewitt算子锐化
[m,n,band]=size(A);
format long
mode1=[-1,-1,-1;0,0,0;1,1,1];
mode2=[-1,0,1;-1,0,1;-1,0,1];
B=A;
for i=1:band
    for j=2:m-1
        for k=2:n-1
            linyu=[A(j-1,k-1,i),A(j,k-1,i),A(j+1,k-1,i);A(j-1,k,i),A(j,k,i),A(j+1,k,i);A(j-1,k+1,i),A(j,k+1,i),A(j+1,k+1,i)];
            B(j,k,i)=abs(sum(sum(linyu.*mode1),2))+abs(sum(sum(linyu.*mode2),2));
        end
    end
end
figure(1)
subplot(1,2,1)
imshow(uint8(A))
title('原始图像')
subplot(1,2,2)
imshow(uint8(B))
title('prewitt算子锐化后的图像')
end

1.8 RGB-HIS转化

function B=rgb_his(A) % rgb-his,A为输入矩阵,已经转为double类型
[row,col,band]=size(A);
% 如果输入图像的波段数!=3,取前三个波段
if band~=3
    disp('输入的不是rgb图像,请选择三个波段');
    A=A(:,:,1:3);
end
% 初始化输出矩阵
B=zeros(row,col,3);
for i =1:row
    for j=1:col
        % 得到中间变量sita角,需要分情形讨论
        if (A(i,j,1)==A(i,j,2))&&(A(i,j,2)==A(i,j,3))
            sita=0;
        else
            fenzi=0.5*(2*A(i,j,1)-A(i,j,2)-A(i,j,3));
            fenmu=pow2(pow2(A(i,j,1)-A(i,j,2),2)+(A(i,j,1)-A(i,j,3)*(A(i,j,2)-A(i,j,3))),0.5);
            sita=180*acos(fenzi/fenmu)/pi;
        end
        % 求解i分量
        B(i,j,3)=(A(i,j,1)+A(i,j,2)+A(i,j,3))/3; 
        % 求解s分量
        B(i,j,2)=1-3/(A(i,j,1)+A(i,j,2)+A(i,j,3))*min([A(i,j,1);A(i,j,2);A(i,j,3)]);
        % 求解h分量
        if A(i,j,2)>=A(i,j,3)
            B(i,j,1)=sita;
        else
            B(i,j,1)=360-sita;
        end
    end
end
% 绘图展示进行对比
figure(1)
subplot(1,2,1)
imshow(uint8(A))
title('RGB图像')
subplot(1,2,2)
imshow(uint8(B))
title('HIS图像')
end

1.9 HIS-RGB转化

function B=his_rgb(A)% HIS转化为RGB
[row,col,~]=size(A);
% 初始化输出矩阵
B=A;
%% 以下是分角度进行讨论,代码就是课件上公式的复刻
%% 重点注意各分量代表的含义
% H对应1,S对应2,I对应3
% R对应1,G对应2,B对应3
for i=1:row
    for j=1:col
        if A(i,j,1)>=0&&A(i,j,1)<120
            B(i,j,3)=A(i,j,3)*(1-A(i,j,2));
            B(i,j,1)=A(i,j,3)*(1+(A(i,j,2)*cos(A(i,j,1))/cos(60-A(i,j,1))));
            B(i,j,2)=3*A(i,j,3)-(B(i,j,1)+B(i,j,3));
        elseif A(i,j,1)<240&&A(i,j,1)>=120
            B(i,j,1)=A(i,j,3)*(1-A(i,j,2));
            B(i,j,2)=A(i,j,3)*(1+(A(i,j,3)*cos(A(i,j,1)-120)/cos(180-A(i,j,1))));
            B(i,j,3)=3*A(i,j,3)-(B(i,j,1)+B(i,j,2));
        else
            B(i,j,2)=A(i,j,3)*(1-A(i,j,2));
            B(i,j,3)=A(i,j,3)*(1+(A(i,j,2)*cos(A(i,j,1)-240)/cos(300-A(i,j,1))));
            B(i,j,1)=3*A(i,j,3)-(B(i,j,2)+B(i,j,3));
        end
    end
end
% 将原始图像与结果图像展示比对
figure(1)
subplot(1,2,1)
imshow(uint8(A))
title('HIS图像')
subplot(1,2,2)
imshow(uint8(B))
title('RGB图像')
end

1.10 PCA变换

function [TZXL,B] = cal_pca(C) % 将A进行主成分分析得到矩阵B,并且将主成分展示出来
[lines,samples,bands] = size(C);
JZ = zeros(bands,lines*samples);
for i = 1:bands
   JZ(i,:) = reshape(C(:,:,i),[1,lines*samples]); % 将三维矩阵变为二维矩阵
end
COV = cal_cov(C); %得到每个波段的协方差矩阵
[TZXL,TZZ] = eigs(COV);  %得到特征值和特征向量
[~,index] = sort(diag(TZZ),'descend');   %将特征值按降序排列
TZXL = TZXL(:,index); %特征值为bands列,特征向量矩阵的每一列代表对应特征值的特征向量
A = TZXL';
B = A * JZ; %B的每一行代表第n各主成分,如果波段数大于3,只要取前三个主成分就可以了
B=reshape(B',[lines,samples,bands]);%将二维矩阵重新变成三维矩阵,便于图像的显示
figure(1)
subplot(1,2,1)
imshow(uint8(C))
title('原图像')
subplot(1,2,2) 
imshow(uint8(B))
title('k-l图像')
end

%% 1.cal-cov计算协方差函数
%% 2.eigs计算特征值和特征向量


2.影像融合处理

2.1 Brovey变换融合算法

function C=Brovey(A,B) % A为全色影像,B为高光谱影像,C为融合后的影像
% 分别得到A与B矩阵的大小便于内插
[m1,n1]=size(A);
[m2,n2,band]=size(B);
% 构造梯度向量F,并且设置参数将高光谱影像大小调整为与全色影像一致
F = griddedInterpolant(double(B));
xq = linspace(0,m2,m1)';
yq = linspace(0,n2,n1)';
zq = (1:band)';
F.Method = 'linear';
vq = F({xq,yq,zq});
% 初始化D矩阵用来B矩阵每个位置光谱向量的和
D=zeros(m1,n1);
C=zeros(m1,n1,band);
% 循环累加得到D矩阵中每个位置的值
for i=1:m1
    for j=1:n1
        D(i,j)=sum(vq(i,j,:));
    end
end
% brovey融合的核心公式:C=A.*B(i)/sum(B(i))
for i=1:band
    C(:,:,i)=A.*vq(:,:,i)./D; %得到每个单波段矩阵
end
% 将输入图像与融合图像展示
figure(1)
subplot(1,3,1)
imshow(uint8(A))
title('全色图像')
subplot(1,3,2) 
imshow(uint8(B(:,:,2:4)))
title('多光谱图像')
subplot(1,3,3) 
imshow(uint8(C(:,:,2:4)))
title('融合后图像')

2.2 乘积变换融合算法

function C=chengji(A,B) % A为全色影像,B为高光谱影像,C为融合后的影像
[m1,n1]=size(A);
[m2,n2,band]=size(B);
%构造梯度向量F,并且设置参数将高光谱影像大小调整为与全色影像一致
F = griddedInterpolant(double(B));
xq = linspace(0,m2,m1)';
yq = linspace(0,n2,n1)';
zq = (1:band)';
F.Method = 'linear';
vq = F({xq,yq,zq});
C=zeros(m1,n1,band);
% 乘积运算的核心公式 C(i)=A.*B(i)
for i=1:band
    C(:,:,i)=A.*vq(:,:,i); %得到每个单波段矩阵
end
% 将融合后的C影像的灰度值变为0-255之间
for j=1:band
C(:,:,j) = round((C(:,:,j) - min(min(C(:,:,j))))*255./(max(max(C(:,:,j))) - min(min(C(:,:,j)))));% 将灰度值复原回0-255
end
% 分三个图框将A,B,以及最终的融合结果输出对比
figure(1)
subplot(1,3,1)
imshow(uint8(A))
title('全色图像')
subplot(1,3,2) 
imshow(uint8(B(:,:,2:4)))
title('多光谱图像')
subplot(1,3,3) 
imshow(uint8(C(:,:,2:4)))
title('融合后图像')

2.3 PCA影像融合算法

function C=pca_ronghe(A,B) % A为全色影像,B为高光谱影像,C为融合后的影像
% 分别得到A与B矩阵的大小便于内插
[m1,n1]=size(A);
[m2,n2,band]=size(B);
% 构造梯度向量F,并且设置参数将高光谱影像大小调整为与全色影像一致
F = griddedInterpolant(double(B));%创建一个重采样对象
xq = linspace(0,m2,m1)';
yq = linspace(0,n2,n1)';
zq = (1:band)';% 重新分割序列,便于插值
F.Method = 'linear'; % 使用线性方式进行插值
vq = round(F({xq,yq,zq})); % F是一个梯度变换向量
% 调用之前编写的函数得到B主成分分析之后的矩阵
[TZXL,vq]=cal_pca(vq);
% 对灰度值取整
vq=round(vq);
% 以pca变换后的高光谱影像的第一主成分为参考对全色影像进行直方图匹配
pipei=zftpp(A(:,:,1),vq(:,:,1)); 
% 匹配结果替代原来的第一主成分
vq(:,:,1)=pipei;
% pca逆变换得到结果
C=vq*TZXL;
% 将输入与输出展示对比
figure(1)
subplot(1,3,1)
imshow(uint8(A))
title('全色图像')
subplot(1,3,2) 
imshow(uint8(B))
title('多光谱图像')
subplot(1,3,3) 
imshow(uint8(A))
title('融合后图像')
end

%%1.cal-pca,pca变换 %%2.zftpp,直方图匹配函数


2.4 HSI变换融合算法

function C=hsi_ronghe(A,B)% his融合,A为全色影像,B为多波段影像
% 获取影像大小
[m1,n1]=size(A);
[m2,n2,band]=size(B);
%创建一个重采样对象
F = griddedInterpolant(double(B));
% 重新分割序列,便于插值
xq = linspace(0,m2,m1)';
yq = linspace(0,n2,n1)';
zq = (1:band)';
% 使用线性方式进行插值
F.Method = 'linear'; 
% F是一个梯度变换向量
vq = round(F({xq,yq,zq})); 
% 调用rgb影像转化为his影像的函数
vq=rgb_his(vq);
% 将数据类型改为可以代数运算的形式
vq=double(vq);
% 以pca变换后的高光谱影像的第一主成分为参考对全色影像进行直方图匹配
pipei=zftpp(A(:,:,1),vq(:,:,1)); 
% 将高光谱影像的第一主成分进行替换
vq(:,:,1)=pipei;
% 将灰度值化为整数
C=round(his_rgb(vq));
% 画图展示比对
figure(1)
subplot(1,3,1)   
imshow(uint8(A))
title('全色图像')
subplot(1,3,2) 
imshow(uint8(B))
title('多光谱图像')
subplot(1,3,3) 
imshow(uint8(C))
title('融合后图像')
end

%% 1.rgb-his,rgb格式转化为his格式图像 %%2.zftpp,直方图匹配函数


评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

会飞的神里绫华

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值