博客介绍了编程基础技巧，如文件读写、随机数等操作。重点讲解遥感数字图像处理，包括图像读取、统计描述、增强处理、融合处理、几何处理、特征提取和分类等内容，涉及多种算法和处理流程，如直方图均衡化、PCA影像融合等。
基础
编程tip
-
能用向量不用循环
i=1:100; mat=mat.*15; -
BIL和BIP运算效率高 ，BSQ读取效率高
-
数组下标是从1开始的！！！！！！！
-
Ctrl+C 强制退出运行
-
像素处理一定将 uint8转换为double
double(img); img1=im2double(img); -
尽量少的用中间参数，减少空间复杂度。有些参数是最后不用的，直接覆盖
文件读写
文件读取
f_id=fopen('filename','r');
num=fread(f_id,[m,n],'int8');
fclose(f_id);
文件存储
f_id=fopen('filenmae','a');
fwrite(f_id,num,'int8');
fclose(f_id);
| r | 行读操作（默打开文件进认形式） |
|---|---|
| r+ | 打开文件进行读和写操作 |
| w | 删除已存在文件中的内容或生成一个新文件，打开进行写操作 |
| w+ | 删除已存在文件中的内容或生成一个新文件，打开进行读和写操作 |
| a | 生成并打开一个新文件或打开一个已存在的文件，进行写操作，在文件末尾追加 |
| a+ | 生成并打开一个新文件或打开一个已存在的文件，进行读和写操作，在文件末尾追加 |
若在上表的字符后加 “b” ，则表示以二进制进行操作。
MAT文件读取
load('filename');
图像读取
imread
[A,map]=imread(filename,fmt)；%%将图像文件filename中的数据读入矩阵A中。
%%map：为彩色图像的调色板，它的值归一化到[0,1]
%fmt：指定了图像的格式
imwrite
imwrite(A,map,filename,fmt);
image与imshow
image(x);colormap(map); %%colormap 是调色板
imshow(I);
随机数
M = magic(3)
%生成一个n*n的矩阵，矩阵元素是由整数1到n^2组成的并且任何行任何列的和都相等，阶数n必须是大于等于3的标量。
M =
8 1 6
3 5 7
4 9 2
reshape
B = reshape(A,m,n) %将矩阵A的元素返回到一个m×n的矩阵B。如果A中没有m×n个元素则返回一个错误。
B=reshape(A,2,[]) %结果为两行
sprintf
将数据格式化为字符串
| 值类型 | 转换 | 详细信息 |
|---|---|---|
| 有符号整数 | %d 或 %i | 以 10 为基数 |
| 无符号整数 | %u | 以 10 为基数 |
%o | 以 8 为基数（八进制） | |
%x | 以 16 为基数（十六进制），小写字母 a–f | |
%X | 与 %x 相同，大写字母 A–F | |
| 浮点数 | %f | 定点记数法（使用精度操作符指定小数点后的位数。） |
%e | 指数记数法，例如 3.141593e+00（使用精度操作符指定小数点后的位数）。 | |
%E | 与 %e 相同，但为大写，例如 3.141593E+00（使用精度操作符指定小数点后的位数）。 | |
%g | 更紧凑的 %e 或 %f，不带尾随零（使用精度操作符指定有效数字位数。） | |
%G | 更紧凑的 %E 或 %f，不带尾随零（使用精度操作符指定有效数字位数。） | |
| 字符或字符串 | %c | 单个字符 |
%s | 字符向量或字符串数组。输出文本的类型与 formatSpec 的类型相同。 |
formatSpec = "The current time is: %d:%d %s";
A1 = 11;
A2 = 20;
A3 = 'a.m.';
str = sprintf(formatSpec,A1,A2,A3)
str =
"The current time is: 11:20 a.m."
str = sprintf('%d',round(pi)) % str='3'
遥感数字图像处理入门
程序涉及的方法、思路、计算、数据处理流程及编程
一、 图像读取
1.读取ENVI格式数据
读取文件三部曲：
fid=fopen(filename,‘r’);
a=fread(fid,[m,n],float);
fclose(fid);
存文件三部曲：
fid=fopen(filename);
fwrite(fid,a,’ ');
fclose(fid);
function image_m=readfile(fname)
%消除尾部空格
fname=deblank(fname);
elements={'samples ' 'lines ' 'bands ' 'data type '};
%将所有类型存入字典，之后可通过下标访问
datatype={'uint8' 'int16' 'int32' 'float32' 'float64' 'uint16' 'uint32' 'int64''uint64'};
if ~ischar(fname)
error("plese input str");
end
% 去掉文件名后缀
corename = strtok(fname,'.');
%搜索头文件
rfid=fopen(strcat(corename,'.hdr'),'r');
while 1
tline = fgetl(rfid);
if ~ischar(tline)
break;
end
%字符串按照=分割
[first,second]=strtok(tline , '=');
switch first
case elements(1)
[f,s]=strtok(second);
%将字符串改为数值，存储列数
sample=str2num(s);
case elements(2)
[f,s]=strtok(second);
%存储图像行数
line=str2num(s);
case elements(3)
[f,s]=strtok(second);
%存储波段数
band=str2num(s);
case elements(4)
[f,s]=strtok(second);
%提取数据格式
t=str2num(s);
switch t
case 1
t=datatype(1);
case 2
t=datatype(2);
case 3
t=datatype(3);
case 4
t=datatype(4);
case 5
t=datatype(5);
case 6
t=datatype(6);
case 7
t=datatype(7);
case 8
t=datatype(8);
case 9
t=datatype(9);
otherwise ,error('wrong');
end
end
end
fclose(rfid);
t=t{1,1};
disp([('打开：'),...
(num2str(sample)),('列×'),(num2str(line)),('行*'),(num2str(band)),('波段')...
('及数据类型是'),('"'),(t),('"'),('的ENVI image...')]);
%根据头文件的信息打开envi格式文件
fid=fopen(fname);
image=fread(fid,[sample,line],t);
image_m=uint8(image');
fclose(fid);
二、 统计与描述
1.均值、中值、众数、方差协方差、相关系数
function st = statics(image)
% 求取均值、中值、众数、方差协方差、相关系数
%求得的st是一个结构体，全部记录所有的统计信息
%利用size求取矩阵尺度
[rows,cols,bands]=size(image);
%% 均值
%返回的均值是一个向量，每个波段一个值
average=zeros(1,bands);
%分波段计算
for band =1:bands
%务必把二维矩阵变为一维
OneBand = image(:,:,band);
%务必把灰度值改为double才能保证累加正确
OneBand=double(OneBand(:));
%求和
for i = 1:rows*cols
average(band)=average(band)+OneBand(i);
end
%将求和值除以像元数
average(band)=average(band)/rows/cols;
end
%% 中值
midvalue=zeros(1,bands);
for band =1:bands
oneband = image(:,:,band);
%务必把灰度值改为double保证计算正确
oneband=double(oneband(:));
%通过排序，求取中间的位置
for i =1:rows*cols
for j = 1:rows*cols-1
if oneband(j)>oneband(j+1)
tp=oneband(j+1);
oneband(j+1)=oneband(j);
oneband(j)=tp;
end
end
end
%序号是中间的即为中值
midvalue(band)=oneband(round(rows*cols/2));
end
%% 众数
mode=zeros(1,bands);
%zft为直方图的意思
%众数求解即为先求取灰度直方图
%求取直方图最大值即为众数
zft=zeros(1,256);
for band = 1:bands
oneband=image(:,:,band);
%务必把灰度值改为double才能保证计算正确
oneband=double(oneband(:));
for i = 1:rows*cols
gray=oneband(i);
%根据下标计算可以有效代替写循环
%更高效地求解
zft(gray+1)=zft(gray+1)+1;
end
nummax=zft(1);maxindex=1;
%求解灰度直方图最大值以及其对应序列号
for i = 1:256
if zft(i)>nummax
nummax=zft(i);
maxindex=i;
end
end
mode(band)=maxindex-1;
end
%% 方差-协方差
%此处计算出的是一个方差-协方差矩阵
vars=zeros(bands,bands);
for iband = 1:bands
ioneband=image(:,:,iband);
%务必把灰度值改为double才能保证计算正确
ioneband=double(ioneband(:));
for jband = 1:bands
joneband=image(:,:,jband(:));
joneband=double(joneband);
re=0;
%根据方差求解公式遍历求和进行求解
for i =1:rows*cols
re=re+(ioneband(i)-average(iband))*(joneband(i)-average(jband));
end
vars(iband,jband)=re/rows/cols;
end
end
%% 相关系数
rou=zeros(bands,bands);
%相关系数即为对应的协方差除以各自方差，
%利用此求解
for i = 1:bands
for j = 1:bands
rou(i,j)=vars(i,j)/sqrt(vars(i,i)*vars(j,j));
end
end
%% 结构赋值
st.average=average;
st.midvalue=midvalue;
st.mode=mode;
st.vars=vars;
st.rou=rou;
end
2.排序
function array = my_sort(array)
%MYSORT 对二维数组进行排序
% 此处显示详细说明
[rows,cols]=size(array);
%将二维数组转为向量
array=array(:);
%运用冒泡排序排序
for i =1:rows*cols
for j = 1:rows*cols-1
if array(j)>array(j+1)
tp=array(j+1);
array(j+1)=array(j);
array(j)=tp;
end
end
end
3.直方图与累计直方图
function LjorZft = ZhiFangTu_(image)
% 计算灰度直方图与累计直方图
%% 初始化
[rows,cols,bands]=size(image);
%累计直方图每个波段一行，一行256列
ljzft=zeros(bands+1,256);
%第一行记录像元值
ljzft(1,:)=0:255;
zft=ljzft;
%% 直方图
for band = 1:bands
%务必把二维矩阵变为一维
oneband=image(:,:,band);
%务必把灰度值改为double保证计算正确
oneband=double(oneband(:));
for i = 1:rows*cols
value=oneband(i);
%根据下标计算可以有效代替写循环
%更高效地求解
zft(band+1,value+1)=zft(band+1,value+1)+1;
end
end
%% 累计直方图
for band = 1:bands
%对于累计直方图 第一个值是直方图的值，之后在前一个基础上加即可
ljzft(band+1,1)=zft(band+1,1);
for i = 2:256
%累计直方图第一个值等于灰度直方图第一个值
%从第二个开始，每个累计直方图的值等于前一个累计直方图的值加这一个灰度直方图的值
ljzft(band+1,i)=ljzft(band+1,i-1)+zft(band+1,i);
end
end
LjorZft.ljzft=ljzft;
LjorZft.zft=zft;
end
三、影像增强处理
1.直方图线性拉伸
变 换 前 亮 度 值 取 值 ： x a ∈ [ a 1 , a 2 ] ; 变 换 后 亮 度 值 取 值 ： x b ∈ [ b 1 , b 2 ] ; x b = b 2 − b 1 a 2 − a 1 ( x a − a 1 ) + b 1 变换前亮度值取值： x_a \in[a_1,a_2];\\ 变换后亮度值取值：x_b \in[b_1,b_2];\\ x_b=\frac{b_2-b_1}{a_2-a_1}(x_a-a_1)+b_1 变换前亮度值取值：xa∈[a1,a2];变换后亮度值取值：xb∈[b1,b2];xb=a2−a1b2−b1(xa−a1)+b1
function xxls = Xianxinglashen_(image,backmin,backmax)
% 线性拉伸计算
% 输入图像，以及要把图像拉伸到什么地步
%% 计算图像行列，波段数
[rows,cols,bands]=size(image);
%% 分波段拉伸
for band = 1:bands
%务必把二维矩阵变为一维，且转为double类型
oneband=image(:,:,band);
oneband=double(oneband(:));
maxnum=-1;minnum=1000;
% 找出最大最小值
images=image;
%通过一个循环，分别求取最大最小值
for i =1:rows*cols
if oneband(i)>maxnum
maxnum=oneband(i);
end
if oneband(i)<minnum
minnum=oneband(i);
end
end
%通过公式计算拉伸后矩阵
oneband = (backmax-backmin)/(maxnum-minnum)*(oneband-minnum)+backmin;
%通过reshape恢复二维矩阵
images(:,:,band)=reshape(oneband,rows,cols);
end
xxls=images;
end
2.直方图均衡化
x b k = L − 1 N ∑ j = 0 k h a ( x a j ) N 是 像 元 总 数 L 是 图 像 灰 度 级 （ 256 ） k = 0 , , , L − 1 h a ( x a j ) 为 原 图 像 第 j 灰 度 值 频 数 x_b^k=\frac{L-1}{N}\sum_{j=0}^{k}h_a(x_a^j)\\ N是像元总数 \qquad L是图像灰度级（256）\qquad k=0,,,L-1\\ h_a(x_a^j)为原图像第j灰度值频数 xbk=NL−1j=0∑kha(xaj)N是像元总数L是图像灰度级（256）k=0,,,L−1ha(xaj)为原图像第j灰度值频数
过程：
- 统计原图像每一灰度级的像元数和累积像元数；
- 计算原每一灰度级均衡化后的新值；
- 以新值代替原灰度值，形成均衡化后的新图像。
function jhhimage = Junhenghua_(image)
% 将图像均衡化处理
%% 计算图像行列，波段数
[rows,cols,bands] = size(image);
%% 循环分波段均衡化
for band = 1:bands
%务必把二维矩阵变为一维，且转为double类型
oneband=image(:,:,band);
oneband=double(oneband(:));
%引用函数求取累计直方图
myzft=ZhiFangTu_(oneband);
%根据公式计算变化后的表
bhh=myzft.ljzft(2,:)*255/rows/cols;
%将像素值四舍五入求整
bhh=round(bhh);
%将变化后的像素值替换原来图像
oneband=bhh(oneband+1); %%%%%这句非常快，但比较难想
image(:,:,band)=reshape(oneband,rows,cols);
end
jhhimage=image;
3.直方图匹配
function ppimage = PiPei_(bimage,cimage)
%直方图匹配
%% 分别计算图像行列，波段数
[rowb,colb,bandsb]=size(bimage);
[rowc,colc,bandsc]=size(cimage);
%循环各波段匹配
for band =1:bandsc
%务必把二维矩阵变为一维，且转为double类型
onebandb=bimage(:,:,band);
onebandb=onebandb(:);
%分别求取两幅图像的累计直方图
zftsb=ZhiFangTu_(onebandb);
onebandc=cimage(:,:,band);
onebandc=onebandc(:);
zftsc=ZhiFangTu_(onebandc);
%将累积直方图标准化
sumb=zftsb.ljzft(2,:)/(rowb*colb);
sumc=zftsc.ljzft(2,:)/(rowc*colc);
match=zeros(1,256);
%求取累积像元统计值与对应参考累积像元统计值距离的最小值
%距离最小值对应的参考累积像元统计值即为标准化后的像元值
for i = 1:256
mindist=5; % 本来就小于1
index=-1;
for j = 1:256
if abs(sumc(i)-sumb(j))<mindist
mindist=abs(sumc(i)-sumb(j));
index=j;
end
end
match(i)=index-1;
end
onebandc=match(onebandc+1);
%通过reshape恢复二维矩阵
cimage(:,:,band)=reshape(onebandc,rowc,colc);
end
ppimage=cimage;
4.均值滤波
function meanlb = MeanLB(image,n)
%均值滤波计算
%image为一个灰度图，n为自适应窗口大小
meanlb=image;
[rows,cols]=size(image);
%根据n计算窗口从中间到两边的距离
jian=(n-1)/2;
for row = 2:rows-1
for col = 2: cols-1
%提取处理的窗口
list=image(row-jian:row+jian,col-jian:col+jian);
%务必把二维矩阵变为一维，且转为double类型
list=double(list(:));
%求取均值，也可以用此前写过的函数计算
average=mean(list);
%把计算的值赋给中间的的像元
meanlb(row,col)=average;
end
end
end
5.中值滤波
function mlb= MidLB(image,n)
%中值滤波计算
%image为一个灰度图，n为自适应窗口大小
[rows,cols]=size(image);
mlb=image;
%根据n计算窗口从中间到两边的距离
jian=(n-1)/2;
for row = 2:rows-1
for col = 2:cols-1
%提取处理的窗口
list=image(row-jian:row+jian,col-jian:col+jian);
%务必把二维矩阵变为一维
list=list(:);
%排序，也可以用此前写过的函数排序
list=sort(list);
%求取向量的中间值
mid=list(floor(n*n/2));
%把计算的值赋给中间的的像元
mlb(row,col)=mid;
end
end
end
6.sobel算子锐化
function sbl = Sobel(image)
% 利用Sobel算子锐化
% image为灰度图
[rows,cols]=size(image);
%初始化sobel结果图，是一个0矩阵
sbl=zeros(rows,cols);
%此处对于边界采取赋0，即不处理的方式
for row = 2:rows-1
for col = 2:cols-1
%提取处理的窗口
list=image(row-1:row+1,col-1:col+1);
%务必把二维矩阵转为double类型
list=double(list);
%根据算子计算x y分量
x=list(3,1)+list(3,3)+list(3,2)*2-list(1,1)-list(1,2)*2-list(1,3);
y=list(1,3)+2*list(2,3)+list(3,3)-list(1,1)-2*list(2,1)-list(3,1);
%分别把x y求绝对值相加，即为该点处的梯度
xy=abs(x)+abs(y);
%把计算的值赋给中间的的像元
sbl(row,col)=xy;
end
end
end
7.Prewitt算子锐化
function prewitt = Prewitt(image)
% 利用Prewitt算子锐化
% image为灰度图
[rows,cols]=size(image);
%初始化Prewitt结果图，是一个0矩阵
prewitt=zeros(rows,cols);
%此处对于边界采取赋0，即不处理的方式
for row = 2:rows-1
for col = 2:cols-1
%提取处理的窗口
list=image(row-1:row+1,col-1:col+1);
%务必把二维矩阵转为double类型
list=double(list);
%根据算子计算x y分量
x=list(3,1)+list(3,3)+list(3,2)-list(1,1)-list(1,2)-list(1,3);
y=list(1,3)+list(2,3)+list(3,3)-list(1,1)-list(2,1)-list(3,1);
%分别把x y求绝对值相加，即为该点处的梯度
xy=abs(x)+abs(y);
%把计算的值赋给中间的的像元
prewitt(row,col)=xy;
end
end
end
8.RGB与HIS互相转换
function hou = rgb_ihs_f(data,judge)
%RGB HSI 互相转
% 其中judge若为0，则将数据转化为IHS
% judge若为1，则将数据转化为RGB
switch judge
%将RGB图像转为HSI
case 0
%将图像转换为0-1之间
sda=single(data)/255;
%取出行列数
[row,col]=size(sda(:,:,1));
%初始化返回图像大小
hou=zeros(row,col,3);
for i=1:row
for j=1:col
%取出三个波段中最大最小值
M=max(sda(i,j,:));
m=min(sda(i,j,:));
%根据公式进行计算
r=M-(M-m)*sda(i,j,1);
g=M-(M-m)*sda(i,j,2);
b=M-(M-m)*sda(i,j,3);
I=(M+m)/2;
%计算H
if M==m
H=0;
end
if M==sda(i,j,1)
H=60*(2+b-g);
end
if M==sda(i,j,2)
H=60*(4+r-b);
end
if M==sda(i,j,2)
H=60*(6+g-r);
end
if M==m
S=0;
end
%根据I的大小分情况计算S
if I<=0.5
S=(M-m)/(M+m);
end
if I>0.5
S=(M-m)/(2-M-m);
end
hou(i,j,1)=I;
hou(i,j,2)=H;
hou(i,j,3)=S;
end
end
case 1
%将HSI转为RGB
[row,col]=size(data(:,:,1));
hou=zeros(row,col,3);
for i=1:row
for j=1:col
%提取各波段
I=data(i,j,1);
H=data(i,j,2);
S=data(i,j,3);
%根据I的大小计算M
if I<=0.5
M=I*(1+S);
end
if I>0.5
M=I+S-I*S;
end
m=2*I-M;
%根据H的大小计算B
if H<60
B=M;
end
if 60<=H&&H<120
B=m+(M+m)*((120-H)/60);
end
if 120<=H&&H<240
B=m;
end
if 240<=H&&H<300
B=m+(M-m)*((H-240)/240);
end
if 300<=H&&H<=360
B=M;
end
%根据H的大小计算R
if H<60
R=m+(M-m)*(H/60);
end
if 60<=H&&H<180
R=M;
end
if 180<=H&&H<240
R=m+(M-m)*((240-H)/60);
end
if 240<=H&&H<=360
R=m;
end
%根据H的大小计算G
if H<120
G=m;
end
if 120<=H&&H<180
G=m+(M-m)*((H-120)/60);
end
if 180<=H&&H<300
G=M;
end
if 300<=H&&H<=360
G=m+(M-m)*((360-H)/60);
end
%将三个波段组合为一个图像
hou(i,j,1)=R;
hou(i,j,2)=G;
hou(i,j,3)=B;
end
end
end
end
9.PCA(K-L)变换
function pca = PCA(image,n)
% PCA计算
% image为原图像，n为最后图层数目
%计算原始图像的尺寸
[rows,cols,bands]=size(image);
%将原始图像每一个波段变为列向量，且为double型
imMatrix=double(reshape(image,[],bands));
%计算方差协方差矩阵，可以用之前写过的函数
varMat=cov(imMatrix);
%计算特征值特征向量
[X,B]=eig(varMat);
X=X';
%初始化输出图像大小
pca=zeros(rows,cols,n);
for newband=1:n
t=zeros(rows*cols,1);
for band=1:bands
%注意由于matlab计算的特征值是由小到大排列
%而pca要求的是从大到小
%所以我们用bands+1-newband对应起来
t=t+imMatrix(:,band)*X(bands+1-newband,band);
end
%重新把图像恢复为二维图像
%由于计算的结果有负值，故我们需要求绝对值
pca(:,:,newband)=uint8(abs(reshape(t,rows,cols)));
end
end
四、影像融合处理
1.Brovey变换融合算法
function brovimg = Brovey(mtl,pan)
% 计算brovey算法
% 输入 多光谱和全色波段
[rows,cols,bands]=size(pan);
brovimg=zeros(rows,cols,bands);
%由于此处提供的多光谱图像与
%全色波段图像大小不一样，所以
%此处将多光谱图像进行裁剪
mtl=double(mtl(:,1:cols,:));
%将全色波段的三个波段合并为一个
pan=(double(pan(:,:,1))+double(pan(:,:,2))+double(pan(:,:,3)))/3;
for band = 1:3
%利用brovey算法计算
oneband=(pan.*mtl(:,:,band))./(mtl(:,:,1)+mtl(:,:,2)+mtl(:,:,3));
brovimg(:,:,band)=oneband;
end
end
2.乘积变换融合算法
function cj = ChengJi(mtl,pan)
% 计算乘积融合算法
% 输入 多光谱和全色波段
[rows,cols,bands]=size(pan);
cj=zeros(rows,cols,bands);
%由于此处提供的多光谱图像与
%全色波段图像大小不一样，所以
%此处将多光谱图像进行裁剪
mtl=double(mtl(:,1:cols,:));
%将全色波段的三个波段合并为一个
pan=(double(pan(:,:,1))+double(pan(:,:,2))+double(pan(:,:,3)))/3;
for band = 1:bands
%利用乘积变换算法计算
oneband=pan.*mtl(:,:,band);
%由于计算出的乘积变换结果超出0-255
%所以需要拉伸到0-1之间方可显示
oneband=(oneband-min(oneband(:)))/(max(oneband(:))-min(oneband(:)));
cj(:,:,band)=(oneband);
end
3.PCA影像融合算法
function pcarh = pcaRH(mtl,pan)
% 计算PCA融合算法
% 输入 多光谱和全色波段
[rows,cols,bands]=size(pan);
%由于此处提供的多光谱图像与
%全色波段图像大小不一样，所以
%此处将多光谱图像进行裁剪
mtl=double(mtl(:,1:cols,:));
%将全色波段的三个波段合并为一个
pan=(double(pan(:,:,1))+double(pan(:,:,2))+double(pan(:,:,3)))/3;
%将原始图像每一个波段变为列向量，且为double型
imMatrix=reshape(mtl,[],bands);
%计算方差协方差矩阵，可以用之前写过的函数
varMatrix=cov(imMatrix);
%计算特征值特征向量
[X,B]=eig(varMatrix);
%初始化输出图像大小
pcarh=zeros(rows,cols,bands);
X=X';
for newband = 1:bands
t=zeros(rows*cols,1);
for i = 1:bands
%注意由于matlab计算的特征值是由小到大排列
%而pca要求的是从大到小
%所以我们用bands+1-newband对应起来
t=t+imMatrix(:,i)*X(bands+1-newband,i);
end
%重新把图像恢复为二维图像
%由于计算的结果有负值，故我们需要求绝对值
pcarh(:,:,newband)=abs(reshape(t,rows,cols));
end
%取出第一主成分层
pp1=pcarh(:,:,1);
%将第一主成分与全色波段线性拉伸至0-255
pp1=Xianxinglashen_(pp1,0,255);
pan=Xianxinglashen_(pan,0,255);
%将全色波段匹配至第一主成分
ppimage=PiPei_(uint8(pp1),uint8(pan));
%将第一主成分替换为全色波段
pcarh(:,:,1)=ppimage;
X=X';
%将主成分图像进行pca逆变换
X=inv(X);
X=X';
%方法类似正变换
for newband = 1:bands
t=zeros(rows*cols,1);
for i = 1:bands
t=t+imMatrix(:,i)*X(bands+1-newband,i);
end
pcarh(:,:,newband)=abs(reshape(t,rows,cols));
end
end
4.HSI变换融合算法
function hsirh = hsiRH(mtl,pan)
% 计算HSI融合算法
% 输入 多光谱和全色波段
[rows,cols,bands]=size(pan);
%由于此处提供的多光谱图像与
%全色波段图像大小不一样，所以
%此处将多光谱图像进行裁剪
%将图像做RGB转为HSI
ihs=rgb_ihs_f(mtl(:,1:cols,:),0);
%提取出HSI中的I图层并转为0-255
i=uint8(ihs(:,:,1)*255);
%将全色波段的三个波段合并为一个
pan=(double(pan(:,:,1))+double(pan(:,:,2))+double(pan(:,:,3)))/3;
%将全色波段直方图匹配到I图层
pp=double(PiPei_(i,pan))/255;
%并将其替换I图层
ihs(:,:,1)=pp;
%将HSI逆变换为RGB
hsirh=rgb_ihs_f(ihs,1);
end
五、影像几何处理
（本章程序属于课外练习，不作要求，但必须掌握思路和流程，了解是如何编程的）
1.影像手动选点配准(image to image)
步骤：
-
打开参考影像与待配准影像
-
选取几何校正的模型（二次多项式等）
-
手动选取多个同名点 若为自动配准则用NCC算法
-
同名点代入模型，计算模型参数
-
确定转换后影像边界
即原始影像四个角点在纠正后的位置
-
重新计算新影像分辨率：
新图像的行数 M＝(Y2-Y1)/△Y+1;
新图像的列数 N＝(X2-X1)/△X+1; -
根据模型逐个像元进行几何变化
-
重采样
如果求得的位置为整数，则该位置处的像元灰度就是新图像的灰度值。若不为整数，则需要用直接或间接法求取位置。
- 最近邻法：算法简单且保持原光谱信息不变；缺点是几何精度较差，图像灰度具有不连续性，边界出现锯齿状。
- 双线性插值：计算较简单，图像灰度具有连续性且采样精度比较精确；缺点是细节丧失。
- 三次卷积法：计算量大，图像灰度具有连续性且采样精度比较精确。
-
量化
六、影像特征提取
1.迭代阈值影像分割算法
function [d2,k,t1]=diedaiyuzhi(image)
%迭代阈值分割计算
%image是灰度图像
%d2是分割后的图像 k是迭代次数 t1是阈值
%获取行列数
[rows,cols]=size(image);
%将image变为一维向量
a=image(:);
%获取该波段影像均值
ave=mean(a);
%给定两个迭代初值
t1=ave;
t2=10000;
d2=zeros(rows*cols,1);
k=0;
%循环迭代求取
while(abs(t1-t2)>0.05)
%初始化两组像素集合
g1=double(0);g2=double(0);
%初始化两组像素集合数目
i1=double(0);i2=double(0);
for p=1:size(a)
%分别向左右两部分增加符合条件的值
if a(p,1)<t1
g1=g1+double(a(p,1));
i1=i1+1;
else
g2=g2+double(a(p,1));
i2=i2+1;
end
end
t2=t1;
%对G1和G2中所有像素计算平均灰度值
t1=(double(g1)/double(i1)+double(g2)/double(i2))/2;
%k记录迭代次数
k=k+1;
end
%将小于阈值的数换为255
%大于阈值的数换为0
for p=1:size(a)
if a(p,1)<t1
d2(p,1)=255;
else
d2(p,1)=0;
end
end
%将d2恢复二维形状并改为无符号8位
d2=reshape(d2,rows,cols);
d2=uint8(d2);
end
七、影像分类
1.K-means
function imageout= K_means(image,k)
% 功能：k-means算法
% image是一个灰度图像
% 输入2：k-分类数
% 输出2： imageout-分类后的图像
%% 将输入图像变为灰度图像
image=double(image);
%image=(image(:,:,1)+image(:,:,2)+image(:,:,3))/3;
image=image(:,:,1);
%一定要再转为整形，以便直方图计算
image=uint8(image);
% 求取图像尺寸
[width,high]=size(image);
ima=double(image);
copy=ima;
%将图像变为一维double型
ima=ima(:);
%求取最小值
mi=min(ima);
ima=ima-mi+1;
s=length(ima);
%% 计算图像灰度直方图
m=max(ima)+1;
%按灰度最大存储为一行
h=zeros(1,m);
%hc存储类别向量
hc=zeros(1,m);
for i=1:s
if(ima(i)>0)
h(ima(i))=h(ima(i))+1;%灰度直方图存储在h[]中
end
end
ind=find(h);
hl=length(ind);%统计灰度级的长度
% 初始化质心
mu=(1:k)*m/(k+1); %灰度按类别数分级
%% 开始分类
while(true)
oldmu=mu;
% 现有的分类
for i=1:hl
%计算各个灰度级距离哪一个分类中心最近
c=abs(ind(i)-mu);
cc=find(c==min(c));
hc(ind(i))=cc(1);
end
%重新计算均值
for i=1:k
a=find(hc==i);
mu(i)=sum(a.*h(a))/sum(h(a));
end
if(mu==oldmu)
break;
end
end
%% 计算分类后的图像
s=size(copy);
mask=zeros(s);
for i=1:s(1)
for j=1:s(2)
%分别判断各个灰度级距离哪一个分类中心最近
c=abs(copy(i,j)-mu);
a=find(c==min(c));
mask(i,j)=a(1);
end
end
%% 显示处理
imageout = ones([width,high]);
%此处为了有视觉效果
%我把灰度级数按照类别数分开
%每个类别一个灰度级数
jishu=255.0/(k+1);
for v=1:k
%找到各自类别所代表的位置
f=find(mask==v);
oneband=zeros(width,high);
oneband(f)=jishu*v;
%将各自类别叠加到图层上
imageout=imageout+oneband;
end
imshow(uint8(imageout)),title('分类后的图像');
end
2.ISODATA
不做要求。
3.最短距离分类器
function flag = zuiduanjuli(data,mc)
%马氏距离的最短距离监督分类
%data为待分类影像
%mc为元胞数组，用于储存已经选取好的分类样本，每一行是一个类别
%矩阵的列数必须与data的波段数相同
[m,n,p]=size(data); %获取图像大小
data=double(data); %转double型数据
[k,~]=size(mc); %获取样本个数
a=reshape(data,m*n,p); %重组图像成为二维数组
flag=zeros(m*n,1); %初始化
for i=1:m*n %遍历所有像元
%首先计算初始值
a1=a(i,:);
%计算第一类别的协方差
c1=inv((cov(mc{1,:})));
%计算第一类别的中心
ave=mean(mc{1,:});
%计算马氏距离
d1=sqrt((a1-ave)*c1*(a1-ave)');
p=1;
%遍历其他类别，找距离最近的类
for j=2:k
%计算类别的协方差
c1=inv((cov(mc{k,:})));
%计算第一类别的中心
ave=mean(mc{k,:});
%计算马氏距离
d=sqrt((a1-ave)*c1*(a1-ave)');
%找距离最近的类 如果先前样本距离较大，则替换成当前样本的分类
if d1>d
p=j;
d1=d;
end
end
flag(i,1)=p; %存放分类的标签结果
end
flag=reshape(flag,m,n);
end
4.最大似然
function flag= zuidasiran(image,mc )
%最大似然分类方法，data为待分类影响
%mc为已经选取好的分类样本,为n*1的元胞数组
[m,n,p]=size(image); %获取图像大小
[k,~]=size(mc); %获取类别个数
a=reshape(image,m*n,p); %重组图像成为二维数组
flag=zeros(m*n,1); %初始化标签结果数据集
ave=zeros(k,1); %初始化各类别的平均值
thi=zeros(k,1); %初始化各类别的协方差
for i=1:k
ave(i,1)=mean(mc{i,1}); %计算各个类别的平均值
thi(i,1)=cov(mc{i,1}); %计算各个类别的协方差
end
for i=1:m*n
a1=a(i,:);
%利用最大似然公式计算像元属于第一类别的概率
f1=exp(-(a1-ave(1,1))^2/(2*thi(1,1)))/sqrt(2*pi*thi(1,1));
p=1;
for j=2:k
%利用最大似然公式计算像元属于其他类别的概率
f2=exp(-(a1-ave(j,1))^2/(2*thi(j,1)))/sqrt(2*pi*thi(j,1));
if f1<f2 %如果下一个分类概率大
p=j; %将这个类别的类别号记录下来
f1=f2;
end
end
if f1<0.3 %如果概率小于0.3
flag(i,1)=k+1; %这一点可能不属于任何类别
else %如果概率大于0.3
flag(i,1)=p; %则为此分类的标签
end
end
flag=reshape(flag,m,n); %将标签结果转成图像形式大小的矩阵
imagesc(flag); %此处应用imagesc以全色图的方式显示
end
扫一扫
6972
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
