《数字图像处理》书中代码

例1-1m=uint8(zeros(128,128,1,27));for i=1:27[m(:,:,:,i),map]=imread('mri.tif',i);endmontage(m,map)例2-1a=imread('cameraman.tif'); %读入cameraman图像figure(1);imshow(a);b1=a+50; %b1=a+45图像灰度值增加45figure(2);imshow(b1);b2=1.2*a; %b=1.2*a图像对比度增大figure(3);imshow(b2)b3=0.65*a; %b=0.65*a图像对比度减少figure(4);imshow(b3);b4=-double(a)+255; %b4=-1*a+255，图像求补，注意把a的类型转换为double figure(5);imshow(uint8(b4)); %再把double类型转换为unit8例2-2%读入lena图像a=imread('cameraman.tif'); %读取原始图像figure(1);imshow(a); %显示原始图像xlabel ('(a)原始图像');%显示函数⎰)(x的曲线图x=1:255;y=x+x.*(255-x)/255;figure(2);plot(x,y); %绘制⎰)(x的曲线图xlabel ('(b)函数⎰)(x的曲线图');b1=double(a)+0.006*double(a) .*(255-double(a));figure(3);imshow(uint8(b1)); %显示非线性处理图像xlabel('(c)非线性处理效果');例2-3clear;%%读入图像a=imread('cameraman.tif');imhist(a);title('原始cameraman图像的直方图');%%b=f(a)b1=1.25*double(a)+45;figure(2); imhist(uint8(b1));title ('变换后的直方图');例2-4clear;histgram=zeros(1,256); %生成直方图数组cdf=zeros(1,256);[cm,map]=imread('cameraman.tif);[a,b]=size(cm);for i=1:afor j=1:bk=cm(i,j);histgram(k)=histgram(k)+1;endend %得到直方图cdf(1)=histgram(1);for i=2:256cdf(i)=cdf(i-1)+histgram(i);endfor i=1:a %点运算for i=1:bk=cm(i,j);cm_equ(i,j)=cdf(k)*256/(a*b);endendimshow(uint8(cm_equ));figure(2);imhist(uint8(cm_equ));例2-5I=imread('tire.tif');J=histeq(I);H=adapthisteq(I);figure(1),imshow(I);xlabel('原始图像');figure(2),imshow(J);xlabel('histeq均衡化');figure(3),imshow(H);xlabel('adapthisteq均衡化');例2-6a=imread('hill.jpg');s=size(a);b=double(a);c(:,:,1)=b(:,:,1)+b(:,:,2);c(:,:,2)=b(:,:,2);c(:,:,3)=b(:,:,3)-b(:,:,2);for i=1:s(1)for j=1:s(2)for k=1:s(3)if c(i,j,k)<0c(i,j,k)=0;endif c(i,j,k)>255c(i,j,k)=255;endendendendc=uint8(c);subplot(121);imshow(a);xlabel('原始图像');subplot(122);imshow(c);xlabel('新图像');例2-7a=imread('eight.tif');a1=imnoise(a,'gaussian',0,0.006); %对原始图像加高斯噪声,共得到4幅图像a2=imnoise(a,'gaussian',0,0.006);a3=imnoise(a,'gaussian',0,0.006);a4=imnoise(a,'gaussian',0,0.006);k=imlincomb(0.25,a1,0.25,a2,0.25,a3,0.25,a4); %线性组合subplot(131);imshow(a);subplot(132);imshow(a1);subplot(133);imshow(k,[]);例2-8clear;a=imread('rice.png'); %读取图像figure(1);imshow(a); %显示原始图像background=imopen(a,strel('disk',15)); %在a上进行形态学运算;ap=imsubtract(a,background); %减法运算函数figure(2);imshow(background); %图像输出背景figure(3);imshow(ap,[]); %减法运算结果例2-9a=imread('hill.jpg');b=imread('bom.jpg');s=size(a);m=s(1);n=s(2);b1=imresize(b,[m n]); %MATLAB实现乘法运算函数a=double(a);c=double(b1);d=a.*c/128;d=uint8(d);subplot(131);imshow(a);subplot(132);imshow(b);subplot(133);imshow(d);例2-10a=imread('lena.bmp');background=imopen(a,strel('disk',15));a1=imdivide(a,background);subplot(131);imshow(a); %原始图像subplot(132);imshow(background); %Background结果subplot(133);imshow(a1,[]); %除法运算结果例3-2clear;load woman; %读入图像数据data=uint8(X);[zipped,info]=huffencode(data); %调用Huffman编码程序进行压缩unzipped=huffdecode(zipped,info); %调用Huffman解码程序进行解码%显示原始图像和经编码后的图像,显示压缩比,并计算均方根误差得erms=0,%表示是Huffman无失真编码subplot(121);imshow(data);subplot(122);imshow(unzipped);erms=compare(data(:),unzipped(:))cr=info.ratiowhos data unzipped zipped%Huffencode函数对输入矩阵vector进行Huffman编码,返回编码后的向量(压缩数据)及相关信息function [zippend,info]=huffencode(vector)%输入和输出都是uint8格式%info返回解码需要的结构信息%info.pad是添加的比特数%info.huffcodes是Huffman码字%info.rows是原始图像行数%info.cols是原始图像列数%info.length是最大码长if~isa(vector,'uint8')eror('input argument must be a uint8 vector');end[m,n]=size(vector);vector=vector(:)';f=frequency(vector); %计算各符号出现的概率symbols=find(f~=0);f=f(symbols);[f,sortindex]=sort(f); %将符号按照出现的概率大小排列symbols=symbols(sortindex);len=length(symbols);symbols_index=num2cell(1:len);codeword_tmp=cell(len,1);while length(f)>1 %生成Huffman树,得到码字编码表index1=symbols_index{1};index2=symbols_index{2};codeword_tmp(index1)=addnode(codeword_tmp(index1),uint8(0));codeword_tmp(index2)=addnode(codeword_tmp(index2),uint8(1));f=[sum(f(1:2)) f(3:end)];symbols_index=[{[index1,index2]} symbols_index(3:end)];[f,sortindex]=sort(f);symbols_index=symbols_index(sortindex);endcodeword=cell(256,1);codeword(symbols)=codeword_tmp;len=0;for index=1:length(vector) %得到整个图像所有比特数len=len+length(codeword{double(vector(index))+1});endstring=repmat(uint(0),1,len);pointer=1;for index=1:length(vector) %对输入图像进行编码code=codeword{double(vector(index))+1};len=length(code);string(pointer+(0:len-1))=code;pointer=pointer+len;endlen=length(string);pad=8-mod(len,8); %非8整数倍时,最后补pad个0if pad>0string=[string uint8(zeros(1,pad))];endcodeword=codeword(symbols);codlen=zeros(size(codeword));weights=2.^(0:23);maxcodelen=0;for index=1:length(codeword)len=length(codeword{index});if len>maxcodelenmaxcodelen=len;endif len>0code=sum(weights(codeword{index}==1));code=bitset(code,len+1);codeword{index}=code;codelen(index)=len;endendcodeword=[codeword{:}];%计算压缩后的向量cols=length(string)/8;string=reshape(string,8,cols);weights=2.^(0:7);zipped=uint8(weights*double(string));%码表存储到一个稀疏矩阵huffcodes=sparse(1,1);for index=1:nnz(codeword)huffcodes(codeword(index),1)=symbols(index);end%填写解码时所需的结构信息info.pad=pad;info.huffcodes=huffcodes;info.ratio=cols./length(vector);info.length=length(vector);info.maxcodelen=maxcodelen;info.rows=m;info.cols=n;%huffdecode函数对输入矩阵vector进行Huffman编码,返回解压后的图像数据function vector=huffdecode(zipped,info,image)if ~isa(zipped,'uint8')error('input argument must be a uint8 vector');end%产生0,1序列,每位占一个字节len=length(zipped);string=repmat(uint8(0),1,len.*8);bitindex=1:8;for index=1:lenstring(bitindex+8.*(index-1))=uint8(bitget(zipped(index),bitindex)); endstring=logical(string(:)');len=length(string);%开始解码weights=2.^(0:51);vector=repmat(uint8(0),1,info.length);vectorindex=1;codeindex=1;code=0;for index=1:lencode=bitset(code,codeindex,string(index));codeindex=codeindex+1;byte=decode(bitset(code,codeindex),info);if byte>0vector(vectorindex)=byte-1;codeindex=1;code=0;vectorindex=vectorindex+1;endendvector=reshape(vector,info.rows,info.cols);%函数addnode添加节点function codeword_new=addnode(codeword_old,item)codeword_new=cell(size(codeword_old));for index=1:length(codeword_old)codeword_new{index}=[item codeword_old{index}];end%函数frequency计算各符号出现的概率function f=frequency(vector)if ~isa(vector,'uint8')error('input a argument must be a uint8 vector');endf=repmat(0,1,256);len=length(vector);for index=0:255f(index+1)=sum(vector==uint8(index));endf=f./len;%函数decode返回码字对应的符号function byte=decode(code,info)byte=info.huffcodes(code);例3-5clear all;format long e;symbol=['abcd'];ps=[0.4 0.2 0.1 0.3];inseq=('dacab');codeword=suanshubianma(symbol,ps,inseq)outseq=suanshujiema(symbol,ps,codeword,length(inseq))%算术编码函数suanshubianmafunction acode=suanshubianma(symbol,ps,inseq)high_range=[];for k=1:length(ps)high_range=[high_range sum(ps(1:k))];endlow_range=[0 high_range(1:length(ps-1))];sbidx=zeros(size(inseq));for i=1:length(inseq)sbidx(i)=find(symbol==inseq(i));endlow=0;high=1;for i=1:length(inseq)range=high-low;high=low+range*high_range(sbidx(i));low=low+range*low_range(sbidx(i));endacode=low;%算术解码函数suanshujiemafunction symbos=suanshujiema(symbol,ps,codeword,symlen) format long ehigh_range=[];for k=1:length(ps)high_range=[high_range sum(ps(1:k))];endlow_range=[0 high_range(1:length(ps)-1)];psmin=min(ps);symbos=[];for i=1:symlenidx=max(find(low_range<=codeword));codeword=codeword-low_range(idx);if abs(codeword-ps(idx))<0.01*psminidx=idx+1;codeword=0;endsymbos=[symbos symbol(idx)];codeword=codeword/ps(idx);if abs(codeword)<0.01*psmini=symlen+1;endend例3-6clear;I1=imread('cameraman.tif');I=im2bw(I1,0.4);[zipped,info]=xingchengbianma(I); %调用xingchengbianma进行编码unzipped=xingchengjiema(zipped,info); %调用xingchengjiema进行解码subplot(131);imshow(I1); %显示原始图像xlabel('原始灰度图像');subplot(132);imshow(I);xlabel('二值图像'); %显示二值图像subplot(133);imshow(unzipped); %显示解码图像xlabel('解码图像');unzipped=uint8(unzipped);%计算均方根误差得erms=0,表示行程编码是无失真编码erms=jfwucha(I(:),unzipped(:))%显示压缩比cr=info.ratiowhos I1 I unzipped zipped%行程编码函数xingchengbianmafunction [zipped,info]=xingchengbianma(vector)[m,n]=size(vector);vector=vector(:)';vector=uint8(vector(:));L=length(vector);c=vector(1);e(1,1)=c;e(1,2)=0;t1=1;for j=1:Lif(vector(j)==c)e(t1,2)=double(e(t1,2))+1;elsec=vector(j);t1=t1+1;e(t1,1)=c;e(t1,2)=1;endendzipped=e;info.rows=m;info.cols=n;[m,n]=size(e);info.ratio=m*n/(info.rows*info.cols);%行程解码函数xingchengjiemafunction unzipped=xingchengjiema(zip,info)zip=uint8(zip);[m,n]=size(zip);unzipped=[];for i=1:msection=repmat(zip(i,1),1,double(zip(i,2)));unzipped=[unzipped section];endunzipped=reshape(unzipped,info.rows,info.cols); unzipped=double(unzipped);%计算均方根误差函数jfwuchafunction erms=jfwucha(f1,f2)e=double(f1)-double(f2);[m,n]=size(e);erms=sqrt(sum(e(:).^2)/(m*n));if erms~=0emax=max(abs(e(:)));[h,x]=hist(e(:));if length(h)>=1figure(2);bar(x,h,'r');e=mat2gray(e,[-emax,emax]);figure(3);imshow(e);endend例3-7%读入Lena图像,用LPCencode进行线性预测编码,用LPCdecode解码I=imread('lena.bmp');x=double(I);y=LPCencode(x);xx=LPCdecode(y);%显示线性误差图,如图3-12(a)所示figure(1);subplot(121);imshow(I);subplot(122);imshow(mat2gray(y));%计算均方根误差,因为是无损编码,那么erms应该为0e=double(x)-double(xx);[m,n]=size(e);erms=sqrt(sum(e(:).^2)/(m*n))%显示原图直方图figure(2);subplot(121);[h,f]=hist(x(:));bar(f,h,'k');%显示预测误差的直方图subplot(122);[h,f]=hist(y(:));bar(f,h,'k')%LPCencode函数用一维无损预测编码压缩图像x,a为预测系数,如果a为默认,则默认%a=1,就是前值预测function y=LPCencode(x,a)error(nargchk(1,2,nargin))if nargin<2a=1;endx=double(x);[m,n]=size(x);p=zeros(m,n); %存入预测值xs=x;zc=zeros(m,1);for i=1:length(a)xs=[zc xs(:,1:end-1)];p=p+a(i)*xs;endy=x-round(p);%LPCdecode函数是解码函数,与编码程序用的同一个预测器function x=LPCdecode(y,a)error(nargchk(1,2,nargin));if nargin<2a=1;enda=a(end:-1:1);[m,n]=size(y);order=length(a);a=repmat(a,m,1);x=zeros(m,n+order);for i=1:nii=i+order;x(:,ii)=y(:,i)+round(sum(a(:,order:-1:1).*x(:,(ii-1):-1:(ii-order)),2));endx=x(:,order+1:end);例3-8%设置压缩比crcr=0.5; %cr=0.5为2:1压缩;cr=0.125为8:1压缩I=imread('lena.bmp'); %图像的大小为512×512I1=double(I)/255; %图像为256级灰度图像,对图像进行归一化操作figure(1);imshow(I1); %显示原始图像%对图像进行FFTfftcoe=blkproc(I1,[8 8],'fft2(x)'); %将图像分割为8×8的子图像进行FFT coevar=im2col(fftcoe,[8 8],'distinct'); %将变换系数矩阵重新排列coe=coevar;[y,ind]=sort(coevar);[m,n]=size(coevar); %根据压缩比确定要变0的系数个数snum=64-64*cr;%舍去不重要的系数for i=1:ncoe(ind(1:snum),i)=0; %将最小的snum个变换系数清0endb2=col2im(coe,[8 8],[512 512],'distinct'); %重新排列系数矩阵%对子图像块进行FFT逆变换获得各个子图像的复原图像,并显示压缩图像I2=blkproc(b2,[8 8],'ifft2(x)'); %对截取后的变换系数进行FFT逆变换figure(2);imshow(I2);%计算均方根误差ermse=double(I1)-double(I2);[m,n]=size(e);erms=sqrt(sum(e(:).^2)/(m*n))例3-9%设置压缩比crcr=0.5; %cr=0.5为2:1压缩;cr=0.125为8:1压缩I=imread('lena.bmp'); %图像的大小为512×512I1=double(I)/255; %图像为256级灰度图像,对图像进行归一化操作figure(1);imshow(I1); %显示原始图像%对图像进行DCTt=dctmtx(8);dctcoe=blkproc(I1,[8 8],'P1*x*P2',t,t');coevar=im2col(dctcoe,[8 8],'distinct');coe=coevar;[y,ind]=sort(coevar);[m,n]=size(coevar); %根据压缩比确定要变0的系数个数%舍去不重要的系数snum=64-64*cr;for i=1:ncoe(ind(1:snum),i)=0; %将最小的snum个变换系数清0 endb2=col2im(coe,[8 8],[512 512],'distinct'); %重新排列系数矩阵%对截取后的变换系数进行DCT逆变换I2=blkproc(b2,[8 8],'P1*x*P2',t',t); %对截取后的变换系数进行DCT逆变换figure(2);imshow(I2);%计算均方根误差ermse=double(I1)-double(I2);[m,n]=size(e);erms=sqrt(sum(e(:).^2)/(m*n))例3-10%设置压缩比crcr=0.5; %cr=0.5为2:1压缩;cr=0.125为8:1压缩I=imread('lena.bmp'); %图像的大小为512×512I1=double(I)/255; %图像为256级灰度图像,对图像进行归一化操作figure(1);imshow(I1); %显示原始图像%对图像进行哈达玛变换t=hadamard(8);htcoe=blkproc(I1,[8 8],'P1*x*P2',t,t);coevar=im2col(htcoe,[8 8],'distinct');coe=coevar;[y,ind]=sort(coevar);[m,n]=size(coevar); %根据压缩比确定要变0的系数个数%舍去不重要的系数snum=64-64*cr;for i=1:ncoe(ind(1:snum),i)=0; %将最小的snum个变换系数清0endb2=col2im(coe,[8 8],[512 512],'distinct'); %重新排列系数矩阵%对截取后的变换系数进行哈达玛逆变换I2=blkproc(b2,[8 8],'P1*x*P2',t,t); %对截取后的变换系数进行哈达玛逆变换I2=I2./(8*8);figure(2);imshow(I2);%计算均方根误差ermse=double(I1)-double(I2);[m,n]=size(e);erms=sqrt(sum(e(:).^2)/(m*n))3.5.4节代码%jpegencode函数用来压缩图像,是一种近似的JPEG方法function y=jpegencode(x,quality)%x为输入图像%quality决定了截去的系数和压缩比error(nargchk(1,2,nargin)); %检查输入参数if nargin<=2quality=1; %默认时quality=1endx=double(x)-128; %像素层次移动-128[xm,xn]=size(x);t=dctmtx(8); %得到8*8 DCT矩阵%将图像分割成8*8子图像,进行DCT,然后进行量化y=blkproc(x,[8 8],'P1*x*P2',t,t');m=[16 11 10 16 24 40 51 6112 12 14 19 26 58 60 5514 13 16 24 40 57 69 5614 17 22 29 51 87 80 6218 22 67 56 68 109 103 7724 35 55 64 81 104 113 9249 64 78 87 103 121 120 9272 92 95 98 112 100 103 99]*quality;%用m量化步长短对变换矩阵进行量化,即根据式(3-37)量化yy=blkproc(y,[8 8],'round(x./P1)',m); %将图像块排列成向量y=im2col(yy,[8 8],'distinct'); %得到列数,也就是子图像个数xb=size(y,2); %变换系数排列次序order=[1 9 2 3 10 17 25 18 11 4 5 12 19 26 33...41 34 27 20 13 6 7 14 21 28 35 42 49 57 50...43 36 29 22 15 8 16 23 30 37 44 51 58 59 52...45 38 31 24 32 39 46 53 60 61 54 47 40 48 55 62 63 56 64];%用Z形扫描方式对变换系数重新排列y=y(order,:);eob=max(x(:))+1; %创建一个块结束符号num=numel(y)+size(y,2);r=zeros(num,1);count=0;% 将非零元素重新排列放到r中,-26-3 eob -25 1 eobfor j=1:xbi=max(find(y(:,j))); %每次对一列(即一块)进行操作if isempty(i)i=0;endp=count+1;q=p+i;r(p:q)=[y(1:i,j);eob]; %截去0并加上结束符号count=count+i+1;endr((count+1):end)=[]; %删除r的没有用的部分r=r+128;%保存编码信息y.size=uint16([xm,xn]);y.numblocks=uint16(xb);y.quality=uint16(quality*100);%对r进行Huffman编码[y.huffman ]=huffencode(uint8(r));%jpegdecode函数,jpegencode的解码程序function x=jpegdecode(y)error(nargchk(1,1,nargin)); %检查输入参数m=[16 11 10 16 24 40 51 6112 12 14 19 26 58 60 5514 13 16 24 40 57 69 5614 17 22 29 51 87 80 6218 22 67 56 68 109 103 7724 35 55 64 81 104 113 9249 64 78 87 103 121 120 9272 92 95 98 112 100 103 99];order=[1 9 2 3 10 17 25 18 11 4 5 12 19 26 33...41 34 27 20 13 6 7 14 21 28 35 42 49 57 50...43 36 29 22 15 8 16 23 30 37 44 51 58 59 52...45 38 31 24 32 39 46 53 60 61 54 47 40 48 55 62 63 56 64]; rev=order; %计算逆运算for k=1:length(order)rev(k)=find(order==k);end%ff=max(rev(:))+1;m=double(y.quality)/100*m;xb=double(y.numblocks); %得到图像的块数sz=double(y.size);xn=sz(1); %得到行数xm=sz(2); %得到列数x=huffdecode(y.huffman,); %Huffman解码x=double(x)-128;eob=max(x(:)); %得到块结束符z=zeros(64,xb);k=1;for i=1:xbfor j=1:64if x(k)==eobk=k+1;break;elsez(j,i)=x(k);k=k+1;endendendz=z(rev,:); %恢复次序x=col2im(z,[8 8],[xm xn],'distinct'); %重新排列成图像块x=blkproc(x,[8 8],'x.*P1',m); %逆量化t=dctmtx(8);x=blkproc(x,[8 8],'P1*x*P2',t',t); %DCT逆变换x=uint8(x+128); %进行位移例3-11clear all;x=imread('F:\教材各章图稿\第2章图\第2章图\lena.bmp'); %读入原始图像subplot(121); %显示原始图像imshow(x);y=jpegencode(x,5); %进行近似的JPEG编码X=jpegecode(y); %进行解码subplot(122);imshow(X); %显示压缩图像%计算均方根误差e=double(x)-double(X);[m,n]=size(e);erms=sqrt(sum(e(:).^2)/(m*n))%计算压缩比cr=imageratio(x,y)例4-1I=imread('fenj.jpg'); %读入清晰图像subplot(121);imshow(I); %显示原始图像len=30; %设置运动位移为30个像素theta=45; %设置运动角度为45度psf=fspecial('motion',len,theta); %建立二维仿真线性运动滤波器psf I1=imfilter(I,psf,'circular','conv'); %用psf产生退化图像subplot(122);imshow(I1); %显示运动后的图像例4-2[I1,map]=imread('fenj-mf.jpg'); %读入运动模糊图像figure(1);imshow(I1);len=30;theta=45;initpsf=fspecial('motion',len,theta); %建立复原点扩散函数[J,P]=deconvblind(I1,initpsf,30); %去卷积figure(2);imshow(J); %显示结果图像如图4-7c(所示)figure(3);imshow(P,[],'notruesize'); %显示复原点扩散函数如图4-7(b)所示例4-3I=imread('meihua.jpg'); %读入图像subplot(221);imshow(I);title('原始图像');H=fspecial('motion',30,45); %运动模糊PSFMotionBlur=imfilter(I,H); %卷积subplot(222);imshow(MotionBlur);title('运动模糊图像');H=fspecial('disk',10); %圆盘状模糊PSFbulrred=imfilter(I,H);subplot(223);imshow(bulrred);title('圆盘状模糊图像');H=fspecial('unsharp'); %钝化模糊PSFSharpened=imfilter(I,H);subplot(224);imshow(Sharpened);title('钝化模糊图像');例5-1A=imread('cat.jpg');figure;subplot(121);imshow(A);A=double(A);A_move=zeros(size(A));H=size(A);A_x=50;A_y=50;I_movesult(A_x+1:H(1),A_y+1:H(2),1:H(3))=I(1:H(1)-A_x,1:H(2)-A_y,1:H(3));subplot(122);imshow(uint8(A_move));例5-2I=imread('hua1.jpg');I1=imcrop(I,[80 60 50 50]);figure;subplot(121);imshow(I);subplot(122);imshow(I1);例5-3I=imread('fengj.jpg')n=50;m=50;figure;imshow(I);for i=1:50n=n-1;m=m+1;I1=imcrop(I,[n,n,m,m]);figure;imshow(I1);end例5-5I=imread('fish.jpg'); %I为原始图像figure;subplot(131);imshow(I); %显示原始图像I=double(I);I_en=imresize(I,4,'nearest'); %最近邻法标志函数nearest扩大4倍subplot(132);imshow(uint8(I_en)); %显示扩大4倍后的图像I_re=imresize(I,0.5,'nearest'); %缩小两倍subplot(133);imshow(uint8(I_re));%显示缩小2倍后的图像例5-6clear all;I=imread('che.jpg'); %本书用的图片都在当前目录下，故不用写路径subplot(121);imshow(I);I1=imrotate(I,30,'crop');subplot(122);imshow(I1);例5-7I=imread('che.jpg');for i=1:25I1=imrotate(I,15*i,'crop');imshow(I1);end例5-8I=imread('mao.jpg');figure;subplot(121);imshow(I);I=double(I);h=size(I);I1=zeros(h(1)+round(h(2)*tan(pi/6)),h(2),h(3));for m=1:h(1)for n=1:h(2)I1(m+round(n*tan(pi/6)),n,1:h(3))=I(m,n,1:h(3));endendI2=uint8(I1);subplot(122);imshow(I2);例5-9I=imread('ta.jpg');figure;subplot(221);imshow(I);I=double(I);h=size(I);I_fliplr(1:h(1),1:h(2),1:h(3))=I(1:h(1),h(2):-1:1,1:h(3)); %水平镜像变换I1=uint8(I_fliplr);subplot(222);imshow(I1);I_flipud(1:h(1),1:h(2),1:h(3))=I(h(1):-1:1,1:h(2),1:h(3)); %垂直镜像变换I2=uint8(I_flipud);subplot(223);imshow(I2);I_fliplr_flipud(1:h(1),1:h(2),1:h(3))=I(h(1):-1:1,h(2):-1:1,1:h(3)); %对角镜像变换I3=uint8(I_fliplr_flipud);subplot(224);imshow(I3);例5-10I=imread('moon.jpg'); %读取原始图像figure;subplot(121);imshow(I);I=double(I);I1=zeros(size(I))+255;h=size(I1);I1(50+1:h(1),50+1:h(2),1:h(3))=I(1:h(1)-50,1:h(2)-50,1:h(3)); %平移变换I2(1:h(1),1:h(2),1:h(3))=I1(h(1):-1:1,1:h(2),1:h(3)); %镜像变换I3=imrotate(I2,45,'nearest'); %旋转变换I4=imresize(I3,3,'nearest'); %扩大3倍I4=uint8(I4);subplot(122);imshow(I4);例5-11T=maketform('affine',[0.5 0 0;0.5 1 0;0 0 1]);I=imread('cameraman.tif');I1=imtransform(I,T);figure;subplot(121);imshow(I);subplot(122);imshow(I1);例5-12I=imread('yazi.jpg'); %读取原始图像I1=imresize(I,[60 60]);figure;subplot(121);imshow(I);T=maketform('projective',[1 1;31 1;31 31;1 31],[5 5;40 5;35 30;-10 30]);T1=makeresampler('cubic','circular');T2=imtransform(I1,T,T1,'size',[150 200],'XYScale',1);subplot(122);imshow(T2);例5-13I=imread('xibao.bmp'); %读取原始图像I1=double(I);I2=floor((I1(:,:,1)+I1(:,:,2)+I1(:,:,3))/3);A=[-1 -1 -1;-1 8 -1;-1 -1 -1];for i=2:15for j=2:15L=I2(i-1:i+1,j-1:j+1).*A;I3(i,j)=sum(sum(L));endendI3imshow(I3);例5-14I=imread('bao.bmp');I1=rgb2gray(I);figure;subplot(121);imshow(I1);fun=inline('median(x(:))');B=nlfilter(I1,[3 3],fun);subplot(122);imshow(B);例5-15I=imread(bao.bmp');I1=rgb2gray(I);figure;subplot(121);imshow(I1);fun=inline('mean(mean(x))');B=nlfilter(I1,[3 3],fun);subplot(122);image(B);axis off;例5-16I=imread('boy1.gif');figure;subplot(121);imshow(I);B=uint8(colfilt(I,[3 3],'sliding','mean'));subplot(122);imshow(B);例5-17I=imread('xibao.bmp');I1=rgb2gray(I);I2=double(I1);fun=inline('ones(16,1)*mean(x)');B=floor(colfilt(I2,[4 4],'distinct',fun))例5-18I=imread(xibao.bmp');I1=rgb2gray(I);fun=@dct2;I2=blkproc(I1,[8 8],fun);imagesc(I2);例5-19I=imread('fengja.jpg');subplot(121);imshow(I);C1=[10 40 90];C2=[30 98 54];I2=roipoly(I,C1,C2);subplot(122);imshow(I2);例5-20[I1,m]=imread('fengja.jpg');I2=rgb2gray(I1);I3=double(I2);figure;subplot(121);imshow(I3,m);C1=[10 40 90];C2=[30 98 54];BW=roipoly(I3,C1,C2);B=double(BW)*256;AB=(I3+B);subplot(122);imshow(AB,m);例5-21[I1,m]=imread('fengja.jpg');I2=rgb2gray(I1);I3=double(I2);figure;subplot(121);imshow(I3,m);C1=10:pi:180;C2=floor(20*cos(C1));BW=roipoly(I3,C1,C2);B=double(BW)*256;AB=(I3+B);subplot(122);imshow(AB,m);例6-2X=ones(100,100); %设定原始图像,大小为100×100 X(65:85,73:78)=0; %设定图像中间的黑带figure(1);imshow(X,'notruesize');F=fft2(X); %二维快速Fourier变换F1=abs(F);%结果输出figure(2);imshow(F1);colormap(jet);%阵元平移I1=[1,-1;-1,1];I2=repmat(I1,50,50);Y=X.*12;F2=abs(fft2(Y));figure(3);imshow(F2);colormap(jet);colorbar;例6-4figure(1);load imdemos saturn2; %将入MATLAB图像saturn2subplot(121);imshow(saturn2); %显示原始图像S=fftshift(fft2(saturn2)); %计算傅里叶变换并移位subplot(122);imshow(log(abs(S)),[]); %显示傅里叶变换谱例6-5I=imread('cat.jpg'); %读入真彩色figure(1);subplot(131);imshow(I);B=rgb2gray(I); %将真彩色转换为灰度图像subplot(132);imshow(B);C=fftshift(fft2(B)); %计算傅里叶变换并移位subplot(133);imshow(log(abs(C)),[]); %显示傅里叶变换谱例6-6RGB=imread('gantrycrane.png');I=rgb2gray(RGB); %转化成强度图像BW=edge(I,'canny'); %提取边界[H,T,R]=hough(BW,'RhoResolution',0.5,'ThetaResolution',0.5);%显示原始图像subplot(221);imshow(RGB);title('gantrycrane.png');%显示Hough矩阵subplot(212);imshow(imadjust(mat2gray(H)),'XData',T,'YData',R,'InitialMagnification','fit');title('Hough transform of gantrycrane.png');xlabel('\theta');ylabel('\rho');axis on;axis normal;hold on;colormap(hot);P=houghpeaks(H,5,'threshold',ceil(0.3*max(H(:))));lines=houghlines(BW,T,R,P,'FillGap',5,'MinLength',7);figure;imshow(BW);hold on;max_len=0;for k=1:length(lines)xy=[lines(k).point1;lines(k).point2];plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');%显示线段的开头和结尾plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');%检测最长线段的终点len=norm(lines(k).point1-lines(k).point2);if(len>max_len)max_len=len;xy_long=xy;endend%保存图像imwrite(I,'cha_1.bmp','bmp');imwrite(BW,'cha_2.bmp','bmp');H=H/(max(max(H)));imwrite(H,'cha_3.bmp','bmp');例6-8RGB=imread('gantrycrane.png');I=rgb2gray(RGB); %转化成强度图像BW=edge(I,'canny'); %提取边界figure(1);imshow(I);figure(2);imshow(BW);%计算theta_step=1;theta=0:theta_step:360;[R,xp]=radon(BW,theta);figure(3);imagesc(theta,xp,R);colormap(hot);xlabel('\theta(degrees)');ylabel('X\prime');title('R_{\theta}(X\prime)');colorbar;max_R=max(max(R)); %捕捉强累积点threshold=0.75;[II,JJ]=find(R>=(max_R*threshold));[line_n,d]=size(II);figure(4);imshow(BW);for k=1:line_nj=JJ(k);i=II(k);R_i=(j-1)*theta_step;xp_i=xp(i);[n,m]=size(BW); %计算直线在原始图像中的位置x_origin=m/2+(xp_i)*sin(R_i*pi/180);y_origin=n/2-(xp_i)*cos(R_i*pi/180)+1;x1=1;xe=m;y1=(y_origin-(x1-x_origin)*tan(((R_i)-90)*pi/180));ye=(y_origin-(xe-x_origin)*tan(((R_i)-90)*pi/180));xv=[x1 xe];yv=[y1,ye];hold on;line(xv,yv);plot(x_origin,y_origin,'.r');end例6-9I=zeros(100,100); %建立简单图像I(25:75,25:75)=1;figuresubplot(121);imshow(I); %显示建立的图像theta=0:180; %规定变换角度的范围[R,xp]=radon(I,theta); %计算radon变换subplot(122);imagesc(theta,xp,R); %以图像方式显示变换结果Rtitle('R_{\theta}(X\prime)'); %显示图像标题xlabel('\theta(degrees)'); %显示x坐标ylabel('X\prime'); %显示y坐标set(gca,'Xtick',0:20:180); %设置x坐标刻度colormap(hot); %设置调色板colorbar; %显示当前图像的调色板例6-10RGB=imread('hai1.jpg'); %读入真彩色图像GRAY=rgb2gray(RGB); %将真彩色图像转换为灰度图像figure;subplot(221);imshow(RGB); %显示真彩色图像subplot(222);imshow(GRAY); %显示灰度图像[R,xp]=radon(GRAY,[0 45]); %计算变换角度为0°和60°的Radon变换subplot(223);plot(xp,R(:,1));%显示0°方向上的Radon变换title('R_{0^o}(x\prime)');subplot(224);plot(xp,R(:,2));%显示60°方向上的Radon变换title('R_{60^o}(x\prime)');例6-11%产生测试图像并显示，如图6-30（a）所示I=phantom(256);figure;imshow(I);%计算映射数据，设定几何关系为“弧度”，然后分别给定光束密度为2弧度、1弧度%以及0.25弧度。