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

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

《《数字图像处理》书中代码.docx》由会员分享，可在线阅读，更多相关《《数字图像处理》书中代码.docx（88页珍藏版）》请在冰豆网上搜索。

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

例1-1

m=uint8（zeros（128,128,1,27））;

fori=1:

27

[m（:

:

:

i）,map]=imread（'mri.tif',i）;

end

montage（m,map）

例2-1

a=imread（'cameraman.tif'）;%读入cameraman图像

figure

（1）;imshow（a）;

b1=a+50;%b1=a+45图像灰度值增加45

figure

（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=1:

255;

y=x+x.*（255-x）/255;

figure

（2）;plot（x,y）;%绘制

的曲线图

xlabel（'（b）函数

的曲线图'）;

b1=double（a）+0.006*double（a）.*（255-double（a））;

figure（3）;imshow（uint8（b1））;%显示非线性处理图像

xlabel（'（c）非线性处理效果'）;

例2-3

clear;

%%读入图像

a=imread（'cameraman.tif'）;

imhist（a）;

title（'原始cameraman图像的直方图'）;

%%b=f（a）

b1=1.25*double（a）+45;

figure

（2）;imhist（uint8（b1））;

title（'变换后的直方图'）;

例2-4

clear;

histgram=zeros（1,256）;%生成直方图数组

cdf=zeros（1,256）;

[cm,map]=imread（'cameraman.tif）;

[a,b]=size（cm）;

fori=1:

a

forj=1:

b

k=cm（i,j）;

histgram（k）=histgram（k）+1;

end

end%得到直方图

cdf

（1）=histgram

（1）;

fori=2:

256

cdf（i）=cdf（i-1）+histgram（i）;

end

fori=1:

a%点运算

fori=1:

b

k=cm（i,j）;

cm_equ（i,j）=cdf（k）*256/（a*b）;

end

end

imshow（uint8（cm_equ））;

figure

（2）;imhist（uint8（cm_equ））;

例2-5

I=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-6

a=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）;

fori=1:

s

（1）

forj=1:

s

（2）

fork=1:

s（3）

ifc（i,j,k）<0

c（i,j,k）=0;

end

ifc（i,j,k）>255

c（i,j,k）=255;

end

end

end

end

c=uint8（c）;

subplot（121）;imshow（a）;xlabel（'原始图像'）;

subplot（122）;imshow（c）;xlabel（'新图像'）;

例2-7

a=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-8

clear;

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-9

a=imread（'hill.jpg'）;

b=imread（'bom.jpg'）;

s=size（a）;

m=s

（1）;n=s

（2）;

b1=imresize（b,[mn]）;%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-10

a=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-2

clear;

loadwoman;%读入图像数据

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.ratio

whosdataunzippedzipped

%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（'inputargumentmustbeauint8vector'）;

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）;

whilelength（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）;

end

codeword=cell（256,1）;

codeword（symbols）=codeword_tmp;

len=0;

forindex=1:

length（vector）%得到整个图像所有比特数

len=len+length（codeword{double（vector（index））+1}）;

end

string=repmat（uint（0）,1,len）;

pointer=1;

forindex=1:

length（vector）%对输入图像进行编码

code=codeword{double（vector（index））+1};

len=length（code）;

string（pointer+（0:

len-1））=code;

pointer=pointer+len;

end

len=length（string）;

pad=8-mod（len,8）;%非8整数倍时,最后补pad个0

ifpad>0

string=[stringuint8（zeros（1,pad））];

end

codeword=codeword（symbols）;

codlen=zeros（size（codeword））;

weights=2.^（0:

23）;

maxcodelen=0;

forindex=1:

length（codeword）

len=length（codeword{index}）;

iflen>maxcodelen

maxcodelen=len;

end

iflen>0

code=sum（weights（codeword{index}==1））;

code=bitset（code,len+1）;

codeword{index}=code;

codelen（index）=len;

end

end

codeword=[codeword{:

}];

%计算压缩后的向量

cols=length（string）/8;

string=reshape（string,8,cols）;

weights=2.^（0:

7）;

zipped=uint8（weights*double（string））;

%码表存储到一个稀疏矩阵

huffcodes=sparse（1,1）;

forindex=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编码,返回解压后的图像数据

functionvector=huffdecode（zipped,info,image）

if~isa（zipped,'uint8'）

error（'inputargumentmustbeauint8vector'）;

end

%产生0,1序列,每位占一个字节

len=length（zipped）;

string=repmat（uint8（0）,1,len.*8）;

bitindex=1:

8;

forindex=1:

len

string（bitindex+8.*（index-1））=uint8（bitget（zipped（index）,bitindex））;

end

string=logical（string（:

）'）;

len=length（string）;

%开始解码

weights=2.^（0:

51）;

vector=repmat（uint8（0）,1,info.length）;

vectorindex=1;

codeindex=1;

code=0;

forindex=1:

len

code=bitset（code,codeindex,string（index））;

codeindex=codeindex+1;

byte=decode（bitset（code,codeindex）,info）;

ifbyte>0

vector（vectorindex）=byte-1;

codeindex=1;

code=0;

vectorindex=vectorindex+1;

end

end

vector=reshape（vector,info.rows,info.cols）;

%函数addnode添加节点

functioncodeword_new=addnode（codeword_old,item）

codeword_new=cell（size（codeword_old））;

forindex=1:

length（codeword_old）

codeword_new{index}=[itemcodeword_old{index}];

end

%函数frequency计算各符号出现的概率

functionf=frequency（vector）

if~isa（vector,'uint8'）

error（'inputaargumentmustbeauint8vector'）;

end

f=repmat（0,1,256）;

len=length（vector）;

forindex=0:

255

f（index+1）=sum（vector==uint8（index））;

end

f=f./len;

%函数decode返回码字对应的符号

functionbyte=decode（code,info）

byte=info.huffcodes（code）;

例3-5

clearall;

formatlonge;

symbol=['abcd'];

ps=[0.40.20.10.3];

inseq=（'dacab'）;

codeword=suanshubianma（symbol,ps,inseq）

outseq=suanshujiema（symbol,ps,codeword,length（inseq））

%算术编码函数suanshubianma

functionacode=suanshubianma（symbol,ps,inseq）

high_range=[];

fork=1:

length（ps）

high_range=[high_rangesum（ps（1:

k））];

end

low_range=[0high_range（1:

length（ps-1））];

sbidx=zeros（size（inseq））;

fori=1:

length（inseq）

sbidx（i）=find（symbol==inseq（i））;

end

low=0;

high=1;

fori=1:

length（inseq）

range=high-low;

high=low+range*high_range（sbidx（i））;

low=low+range*low_range（sbidx（i））;

end

acode=low;

%算术解码函数suanshujiema

functionsymbos=suanshujiema（symbol,ps,codeword,symlen）

formatlonge

high_range=[];

fork=1:

length（ps）

high_range=[high_rangesum（ps（1:

k））];

end

low_range=[0high_range（1:

length（ps）-1）];

psmin=min（ps）;

symbos=[];

fori=1:

symlen

idx=max（find（low_range<=codeword））;

codeword=codeword-low_range（idx）;

ifabs（codeword-ps（idx））<0.01*psmin

idx=idx+1;

codeword=0;

end

symbos=[symbossymbol（idx）];

codeword=codeword/ps（idx）;

ifabs（codeword）<0.01*psmin

i=symlen+1;

end

end

例3-6

clear;

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.ratio

whosI1Iunzippedzipped

%行程编码函数xingchengbianma

function[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;

forj=1:

L

if（vector（j）==c）

e（t1,2）=double（e（t1,2））+1;

else

c=vector（j）;

t1=t1+1;

e（t1,1）=c;

e（t1,2）=1;

end

end

zipped=e;

info.rows=m;

info.cols=n;

[m,n]=size（e）;

info.ratio=m*n/（info.rows*info.cols）;

%行程解码函数xingchengjiema

functionunzipped=xingchengjiema（zip,info）

zip=uint8（zip）;

[m,n]=size（zip）;

unzipped=[];

fori=1:

m

section=repmat（zip（i,1）,1,double（zip（i,2）））;

unzipped=[unzippedsection];

end

unzipped=reshape（unzipped,info.rows,info.cols）;

unzipped=double（unzipped）;

%计算均方根误差函数jfwucha

functionerms=jfwucha（f1,f2）

e=double（f1）-double（f2）;

[m,n]=size（e）;

erms=sqrt（sum（e（:

）.^2）/（m*n））;

iferms~=0

emax=max（abs（e（:

）））;

[h,x]=hist（e（:

））;

iflength（h）>=1

figure

（2）;

bar（x,h,'r'）;

e=mat2gray（e,[-emax,emax]）;

figure（3）;

imshow（e）;

end

end

例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应该为0

e=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,就是前值预测

functiony=LPCe