MATLAB做遥感数字图像处理(1)

基础

编程tip

1. 能用向量不用循环

   i=1:100;
   mat=mat.*15;
   

2. BIL和BIP运算效率高 ，BSQ读取效率高

3. 数组下标是从1开始的！！！！！！！

4. Ctrl+C 强制退出运行

5. 像素处理一定将 uint8转换为double

   double(img);   
   img1=im2double(img);
   

6. 尽量少的用中间参数，减少空间复杂度。有些参数是最后不用的，直接覆盖

文件读写

文件读取

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'

遥感数字图像处理入门

图像统计

直方图和累计直方图

function LjorZft = ZhiFangTu_(image)
[rows,cols,bands]=size(image);
ljzft=zeros(bands+1,256)  %累计直方图每个波段一行，一行256列
ljzft(1,:)=0:255;
zft=ljzft;
%% 直方图
for band = 1:bands
	oneband=image(:,:,band);
	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

统计值

function st = statics(image)
[rows,cols,bands]=size(image);
%% 均值
average=zeros(1,bands);
for band =1:bands
    OneBand = image(:,:,band);
    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);
	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=zeros(1,256);
for band = 1:bands
	oneband=image(:,:,band);
	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);
	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

图像处理

线性拉伸

function xxls =  Xianxinglashen_(image,backmin,backmax)
[rows,cols,bands]=size(image);
for band = 1:bands
	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-maxnum)+backmax;
	images(:,:,band)=reshape(oneband,rows,cols);
end
xxls=images;
end

直方图均衡化

function jhhimage = Junhenghua_(image)
[rows,cols,bands] = size(image);
for band = 1:bands
	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;

直方图匹配

function ppimage = PiPei_(bimage,cimage)
[rowb,colb,bandsb]=size(bimage);
[rowc,colc,bandsc]=size(cimage);
for band =1:bandsc
    onebandb=bimage(:,:,band);
    onebandb=onebandb(:);
    zftsb=ZhiFangTu_(onebandb);
    onebandc=cimage(:,:,band);
    onebandc=onebandc(:);
    zftsc=ZhiFangTu_(onebandc);
    sumb=zftsb.ljzft(2,:)/(rowb*rowc);
    sumc=zftsc.ljzft(2,:)/(rowb*rowc);
    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)=j-1;
    end
    onebandc=match(onebandc+1);
    cimage(:,:,band)=reshape(onebandc,rowc,colc);
end
ppimage=cimage;