基于小波分析和中值滤波 图像处理 matlab程序

程序是下面这样，但只能处理长宽一样的方形图像，灰度和彩祥宽谨色图像都可，要用其他图像只需把Lena.bmp改为其他图像，但图像要保存在m文件所在路径下。下面巧耐还有一个运行后的图像（之一），希望能对你有所帮助。

clearclc

%%%%%%%%%%测试图像只能是方形图像，长宽像素一样。

f=imread('Lena.bmp')%%读取图像数据，图像只能保存在m文件所在的路径下

d=size(f)

if length(d)>2

f=rgb2gray((f))%%%%%%%%如果是彩色图像则转化为灰度图

end

T=d(1)

SUB_T=T/2

% 2.进行二维小波分解

l=wfilters('db10','l') % db10（消失矩为10)低通分解滤波器冲击响应（长度为20）

L=T-length(l)

l_zeros=[l,zeros(1,L)] % 矩阵行数与输入图像一致，为2的整数幂

h=wfilters('db10','h') % db10（消失矩为10)高通分解滤波器冲击响应（长度为20）

h_zeros=[h,zeros(1,L)] % 矩阵行数与输入图像一致，为2的整数幂

for i=1:T % 列变换

row(1:SUB_T,i)=dyaddown( ifft( fft(l_zeros).*fft(f(:,i)') ) ).' % 圆周卷积<->FFT

row(SUB_T+1:T,i)=dyaddown( ifft( fft(h_zeros).*fft(f(:,i)') ) ).' % 圆周卷积<->FFT

end

for j=1:T % 行变换

line(j,1:SUB_T)=dyaddown( ifft( fft(l_zeros).*fft(row(j,:)) ) ) % 圆周卷积<->FFT

line(j,SUB_T+1:T)=dyaddown( ifft( fft(h_zeros).*fft(row(j,:)) ) ) % 圆周谨基卷积<->FFT

end

decompose_pic=line % 分解矩阵

% 图像分为四块

lt_pic=decompose_pic(1:SUB_T,1:SUB_T) % 在矩阵左上方为低频分量--fi(x)*fi(y)

rt_pic=decompose_pic(1:SUB_T,SUB_T+1:T) % 矩阵右上为--fi(x)*psi(y)

lb_pic=decompose_pic(SUB_T+1:T,1:SUB_T) % 矩阵左下为--psi(x)*fi(y)

rb_pic=decompose_pic(SUB_T+1:T,SUB_T+1:T) % 右下方为高频分量--psi(x)*psi(y)

% 3.分解结果显示

figure(1)

subplot(2,1,1)

imshow(f,[]) % 原始图像

title('original pic')

subplot(2,1,2)

image(abs(decompose_pic)) % 分解后图像

title('decomposed pic')

figure(2)

% colormap(map)

subplot(2,2,1)

imshow(abs(lt_pic),[]) % 左上方为低频分量--fi(x)*fi(y)

title('\Phi(x)*\Phi(y)')

subplot(2,2,2)

imshow(abs(rt_pic),[]) % 矩阵右上为--fi(x)*psi(y)

title('\Phi(x)*\Psi(y)')

subplot(2,2,3)

imshow(abs(lb_pic),[]) % 矩阵左下为--psi(x)*fi(y)

title('\Psi(x)*\Phi(y)')

subplot(2,2,4)

imshow(abs(rb_pic),[]) % 右下方为高频分量--psi(x)*psi(y)

title('\Psi(x)*\Psi(y)')

% 5.重构源图像及结果显示

% construct_pic=decompose_matrix'*decompose_pic*decompose_matrix

l_re=l_zeros(end:-1:1) % 重构低通滤波

l_r=circshift(l_re',1)' % 位置调整

h_re=h_zeros(end:-1:1) % 重构高通滤波

h_r=circshift(h_re',1)' % 位置调整

top_pic=[lt_pic,rt_pic] % 图像上半部分

t=0

for i=1:T % 行插值低频

if (mod(i,2)==0)

topll(i,:)=top_pic(t,:)% 偶数行保持

else

t=t+1

topll(i,:)=zeros(1,T) % 奇数行为零

end

end

for i=1:T % 列变换

topcl_re(:,i)=ifft( fft(l_r).*fft(topll(:,i)') )' % 圆周卷积<->FFT

end

bottom_pic=[lb_pic,rb_pic] % 图像下半部分

t=0

for i=1:T % 行插值高频

if (mod(i,2)==0)

bottomlh(i,:)=bottom_pic(t,:) % 偶数行保持

else

bottomlh(i,:)=zeros(1,T) % 奇数行为零

t=t+1

end

end

这个只是一级分解，matlab自带的函数可以实现多级分解，级数由编程者自己确定。

是的，是一样的。

%用于提升变换 --若血残阳

%TransImage.m 输入参数为图像文件贺蠢李的路径加文件名

function TransImage(filename)

f=imread(filename)

f=im2double(f) %将矩阵转化成double型

A=-1.586134342 %提升系数

B=-0.052980118

R=0.882911075

D=0.443506852

K=1.230174105

[m_orgHeight m_orgWidth]=size(f)

h=ones(m_orgHeight,m_orgWidth)

h=im2double(h)

h=h.*128

%第一次列提升开始------------------------------------------------------------------------------------------

f=im2double(f-h)%将矩阵的每个元素减128向直流靠近--

for i=2:2:(m_orgWidth-2) %偶数列提升（C++奇数角标列）

f(:,i)=f(:,i)+A*[f(:,i-1)+f(:,i+1)]

end

f(:,m_orgWidth)=f(:,m_orgWidth)+2*A*f(:,m_orgWidth-1)%最后一列用前一列的2倍来预测

for i=3:2:(m_orgWidth-1)%奇数列提升（C++偶数角标列）

f(:,i)=f(:,i)+B*[f(:,i-1)+f(:,i+1)]

end

f(:,1)=f(:,1)+2*B*f(:,2) %第一列用第二列的2倍来预测

%第一次列提升完毕--------------------------------------------------------------------------------------------

%第二次列提升开始--------------------------------------------------------------------------------------------

for i=2:2:(m_orgWidth-2) %偶数列提升

f(:,i)=f(:,i)+R*[f(:,i-1)+f(:,i+1)]

end

f(:,m_orgWidth)=f(:,m_orgWidth)+2*R*f(:,m_orgWidth-1)

for i=3:2:(m_orgWidth-1) %奇数列提升

f(:,i)=f(:,i)+D*[f(:,i-1)+f(:,i+1)]

end

f(:,1)=f(:,1)+2*D*f(:,2)

%第二次列提升完毕-------------------------------------------------------------------------------------------

for i=1:2:(m_orgWidth-1)%调整起奇、偶列值

f(:,i+1)=f(:,i+1)*(-K) %偶数列（C++奇数角标）

f(:,i)=f(:,i)/K %奇数列（C++偶数禅迟角标）

end

%第一次行提升档乱开始------------------------------------------------------------------------------------------

for i=2:2:(m_orgHeight-2)

f(i,:)=f(i,:)+A*[f(i-1,:)+f(i+1,:)]

end

f(m_orgHeight,:)=f(m_orgHeight,:)+2*A*f(m_orgHeight-1,:)

for i=3:2:(m_orgHeight-1)

f(i,:)=f(i,:)+B*[f(i-1,:)+f(i+1,:)]

end

f(1,:)=f(1,:)+2*B*f(2,:)

%第一次行提升完毕------------------------------------------------------------------------------------------

%第二次行提升开始------------------------------------------------------------------------------------------

for i=2:2:(m_orgHeight-2)

f(i,:)=f(i,:)+R*[f(i-1,:)+f(i+1,:)]

end

f(m_orgHeight,:)=f(m_orgHeight,:)+2*R*f(m_orgHeight-1,:)

for i=3:2:(m_orgHeight-1)

f(i,:)=f(i,:)+D*[f(i-1,:)+f(i+1,:)]

end

f(1,:)=f(1,:)+2*D*f(2,:)

%第二次行提升完毕------------------------------------------------------------------------------------------

for i=1:2:(m_orgHeight-1) %调整奇、偶行值

f(i+1,:)=f(i+1,:)*(-K)%偶数行（C++奇数角标）

f(i,:)=f(i,:)/K%奇数行（C++偶数角标）

end

for i=1:2:(m_orgWidth-1) %把f的奇数列（C语言中的偶数列[角标]）存放到h的前半部分

h(:,(i+1)/2)=f(:,i)

end

for i=2:2:m_orgWidth %把f的偶数列（c语言中的奇数列[角标]）存放到h的后半部分

h(:,m_orgWidth/2+i/2)=f(:,i)

end

for i=1:2:(m_orgHeight-1) %把h的奇数行（c语言中的偶数行[角标]）存放到f的上半部分

f((i+1)/2,:)=h(i,:)

end

for i=2:2:m_orgHeight %把h的偶数行 (c语言中的奇数行[角标])存放到f的下半部分

f(m_orgHeight/2+i/2,:)=h(i,:)

end

f=f.*(-10)%为看清楚高频部分 给矩阵每个元素同时乘以(-10)

h=ones(m_orgHeight/2,m_orgWidth/2) %定义与低频子带阶数相同的“1”阵

h=im2double(h)%将上述阵变换为double型

h=h.*128 %给阵的每个元素乘以128（估计是开始减去的128）

f(1:m_orgHeight/2,1:m_orgWidth/2)=f(1:m_orgHeight/2,1:m_orgWidth/2)./(-10)+h %给低频部分除以（-10）返回原来值，在加上“128”阵，看清楚低频轮廓

imshow(f) %显示变换后图片

生活中需要对一些图像进行处理，比如压缩，去噪，图像增强，图像锐化与钝化，图像融合，图像的分解等，以便对于图像的成分，边缘等细节信息有更加深刻的认识，小波分析由于其固有的时频特性，既可以对图像进行时域分析，也可以对图像进行频率分析，这使得小波分析在图像处理中得到了广泛的应用，本节对其中一些图像处理功能及函数进行讲解：

wavedec2函数用于对图像进行二维小波分解，其函数调用格式如下：

[c,l]=wavedec2(X,n,’wname’)

其中,X表示原始图像，n表示分解层数，wname表示小波函数，c表示各层系数，l表示各层系数对应的长度

ddencmp用于得到全局阀值，其调用格式如下：

[thr,sorh,keepapp]=ddencmp(‘cmp’,’wp’,X)

[thr,sorh,keepapp]=ddencmp(‘cmp’,’wv’,X)

其中cmp表示压缩，wp表示小波包，wv表示小波，X表示原始信号，thr表示阀值，sorh表示阀值类型，s表示软阀值，h表示硬阀值,keepapp=1表示保持近似系数不变

wdencmp用于对数据或图像进行阀值去噪或压缩，其调用格式如下：

[xcomp,c1,l1,perf0,perfl2]=wdencmp(‘gbl’,c,l,’wname’,n,thr,sorh,keepapp)

glb表示利用全局阀值，perf0表示恢复比，perfl2表示压缩比

示例：利用二维小波对图像进行压缩

编写对应的m文件如下：

clc

load woman

subplot(1,2,1)

imshow(X,map)

title('原始图像')

[c,l]=wavedec2(X,3,'sym4')

%%获取全指模局阀值%%

[thr,sorh,keepapp]=ddencmp('cmp','wp',X)

[xcmp,c1,l1,perf0,perfl2]=wdencmp('gbl',c,l,'sym4',3,thr,sorh,keepapp)

subplot(1,2,2)

imshow(xcmp,map)

title('压缩后图片')

程序运行结果如下图：

小波变换用与图像去噪，噪声会影响图像处理的输入，采集，处理的各个环节及输出结果等全过程，因此对于图像的噪声处理是一个不可忽略的重要的问题，去噪已经成为图像处理中不可或缺的一部分

示例：对图像进行二维小波去噪

编写对应的m文件如下：

load julia

%%产生噪声信号%%

init=3718025452

rand('seed',init)

xnoise=X+8*rand(size(X))

colormap(map)

subplot(1,3,1)

imshow(X,map)

title('原始信号')

subplot(1,3,2)

imshow(xnoise,map)

title('含有噪声的信号')

%%获取全局阀值%%

[thr,sorh,keepapp]=ddencmp('den','wp',xnoise)

[xden,c1,l1]=wdencmp('gbl',xnoise,'sym4'如消,3,thr,sorh,keepapp)

subplot(1,3,3)

imshow(xden,map)

title('去除噪声后信号')

程序运行结果如下图：

小波分析用于图像增强，图像增强是对图像进行一定处理，使图像比原图更加清晰，视觉渣逗知效果更好。

示例：利用小波分析对图像进行增强

编写对应的m文件如下：

clc

load facets

subplot(1,2,1)

imshow(X,map)

title('原始信号')

[c,l]=wavedec2(X,3,'sym4')

sizec=size(c)

fori=1:sizec(2)

if(c(i)>250)

c(i)=2*c(i)

else

c(i)=0.5*c(i)

end

end

y=waverec2(c,l,'sym4')

subplot(1,2,2)

imshow(y,map)

title('增强图像')

程序运行结果如下图：

图像钝化

图像的钝化可以在时域中，也可以在频域中，在时域中处理较为简单，只需要加一个平滑滤波器，使图像中每个点与其邻点做平滑处理即可，在此主要说明图像钝化在频域中的处理。图像钝化是为了突出低频信息，弱化高频信息。

示例：对图像进行频域钝化处理，

编写对应的m文件如下：

load chess

subplot(1,2,1)

imshow(X,map)

title('原始图像')

[c,l]=wavedec2(X,3,'db4')

sizec=size(c)

fori=1:sizec(2)

if(c(i)>280)

c(i)=c(i)*2

else

c(i)=c(i)*0.5

end

end

y=waverec2(c,l,'db4')

subplot(1,2,2)

imshow(y,map)

title('采用小波方法钝化图像')

程序运行结果如下图：

图像锐化，与图像钝化刚好相反，是为了突出高频信息，弱化低频信息，从快速变化的成分中分离出系统边界成分，以便进一步识别或者分割等 *** 作。

示例：对图像进行锐化处理

编写对应的m文件如下：

load chess

subplot(1,2,1)

imshow(X,map)

title('原始图像')

[c,l]=wavedec2(X,3,'db5')

sizec=size(c)

%%突出高频信息，弱化低频信息%%

fori=1:sizec(2)

if(abs(c(i))<280)

c(i)=c(i)*2

else

c(i)=c(i)*0.5

end

end

y=waverec2(c,l,'db5')

subplot(1,2,2)

imshow(y,map)

title('采用小波方法锐化图像')

程序运行结果如下图：

小波分析用于图像融合

图像融合是将同一图像的两个部分或者不同图像合成一张图，以便合成之后的图形比原来更容易理解。

示例：利用二维小波变换将两幅图像融合在一起

编写对应的m文件如下：

clear all

load bust

X1=X

map1=map

load woman

X2=X

map2=map

subplot(1,3,1)

imshow(X1,map1)

title('第一幅图像')

subplot(1,3,2)

imshow(X2,map2)

title('第二幅图像')

%%对第二幅图形低频部分和高频部分进行处理%%

fori=1:256

forj=1:256

if(X2(i,j)>120)

X2(i,j)=X2(i,j)*2

else

X2(i,j)=X2(i,j)*0.5

end

end

end

[c1,l1]=wavedec2(X1,3,'sym4')

[c2,l2]=wavedec2(X2,3,'sym4')

%%对图像进行融合%%

c=c1+c2

%%减少图像的亮度%%

c=c*0.5

y=waverec2(c,l1,'sym4')

subplot(1,3,3)

imshow(y,map2)

title('融合后图像')

程序运行结果如下图：

小波分析用于图像分解

对图像分解的目地在于可以更好的观察图像的细节，对图像做出更好的判断,swt2函数用于对图像进行分解，其调用格式如下：

[sa,sh,sv,sd]=swt2(X,N,’wname’)

其中sa，sh,sv,sd分别表示近似系数，水平系数，竖直系数，对角系数，x分解图像，N分解的层数，wname表示小波基名称

示例：对图像进行分解

编写对应的m文件如下：

clear all

load woman

[sa,sh,sv,sd]=swt2(X,3,'db3')

s=1

fori=1:3

subplot(3,4,s)

image(wcodemat(sa(:,:,i),192))

title(['第',num2str(i),'层近似系数'])

subplot(3,4,s+1)

image(wcodemat(sh(:,:,i),192))

title(['第',num2str(i),'层水平系数'])

subplot(3,4,s+2)

image(wcodemat(sv(:,:,i),192))

title(['第',num2str(i),'层竖直系数'])

subplot(3,4,s+3)

image(wcodemat(sd(:,:,i),192))

title(['第',num2str(i),'层对角系数'])

s=s+4

end

程序运行结果如下图：

---