图像处理大作业

重庆科技学院

《图像处理》

课程大作业

院（系）:_电气与信息工程学院_专业班级:计科1004 学生姓名:谢子洋学号: 2010441747 设计地点（单位） I515

设计题目:数字图像处理大作业

完成日期： 2012年月日

指导教师评语: _______________________________________

___________________________________________________________________________

___________________________________________________________________________

成绩（五级记分制）:______ __________

指导教师（签字）:________ ________

1 大作业题目

1.1连续图像f(x.y)和数字图像I(r,c)中各分量的含义是什么？他们有什么联系和区别？取值范围在什么范围？

答：连续图像f（x，y）对应一个2-D数组，这里想x，y表示2D空间XY里面的坐标点的位置，f则表示图像在（x，y）的某种性质F的数值，其中图像的实际尺寸更具图像的实际大小所限制的。

数字图像I（r,c）也对应一个2-D数组，I代表离散化都的f，（r,c）代表离散化后的（x,y）,r代表图像的行，c代表图像的列，这里的I，r,c都是整数。

连续图像f（x,y）代表的连续实际图像，F是指图像在f的各个分量上的值。例如：二值灰度值图像，其中f在x，y分量上只有0 和1 两个值。这个也适用于rgb彩色图像。一个函数可由n阶泰特公式展开，因此，数字图像虽然没有完全刻画Original图像，利用计数机高速度计算，将这个整体的Original图像分割，以点的值存储到相应的坐标点上。

其中 f（x，y）的值是任意实数，而 I, c，r的值都是整数

1.2发光强度及亮度、照度各有什么不同？

答：发光强度只指光源的能量辐射强度；亮度是指扩展光源在某个方向上的单位投影面积的发光强度；照度是指光源照射到母体表面的光通量，它是光源对物体辐射的一种量度。答：1）发光强度，单位坎德拉，即cd。定义：光源在给定方向的单位立体角中发射的光通量定义为光源在该方向的发光强度。

解释：发光强度是针对点光源而言的，或者发光体的大小与照射距离相比比较小的场合。这个量是表明发光体在空间发射的会聚能力的。可以说，发光强度就是描述了光源到底有多“亮”，因为它是光功率与会聚能力的一个共同的描述。发光强度越大，光源看起来就越亮，同时在相同条件下被该光源照射后的物体也就越亮，因此，早些时候描述手电都用这个参数。

2）亮度，单位尼特，即nt。

定义：单位光源面积在法线方向上，单位立体角内所发出的光流。

解释：这个是最容易被误解的概念了。亮度是针对光源而言，而且不是对点光源，是对面光源而言的。无论是主动发光的还是被动（反射）发光的。亮度是一块比较小的面积看起来到底有多“亮”的意思。这个多“亮”，与取多少面积无关，但为了均匀，我们把面积取得比较小，因此才会出现“这一点的亮度”这样的说法。亮度不仅取决于光源的光通量，更取决于等价发光面积和发射的会聚程度。比如激光指示器，尽管其功率很小，但可会聚程度非常高，因此亮度非常

高。3）光照度，单位勒克斯，即lx。

定义：流明的光通量均匀分布在1平方米表面上所产生的光照度。解释：光照度是对被照地点而言的，但又与被照射物体无关。一个流明的光，均匀射到1平方米的物体上，照度就是1 lx。照度的测量，用照度表，或者叫勒克斯表、lux 表。常见照度：阳光直射（正午）下，110,000 lx，阴天室外，1000 lx，商场内，500 lx，阴天有窗室内，100 lx，普通房间灯光下，100 lx，满月照射下，

0.2 lx。

1.3试讨论连续卷积和离散卷积的不同。

答：连续卷和离散卷在定义域上不同，为了避免重叠失误，在离散卷中需要对被卷积系列进行扩展，即要将有限长的被卷积系列A和B扩展成长度为M=A+B-1的序列以得到以M为周期的结果序列。

1.4为什么一般情况下对离散图像的直方图均衡化并不能产生完全平坦的直方

答：因为同意灰度值的各个像素没有理由变换到不同灰度级，所以数字图像直方图均衡化的结果一般不能得到完全均匀分布的直方图，只有近视均衡的直方图。

1.5讨论用于空间滤波的平滑滤波器和锐化滤波器的相同点、不同点和联系。答：相同点：都能削弱或消除傅里叶空间的某些分量，而不影响或较少影响其它分量，从而达到增强某些频率分量的效果。

不同点：平滑滤波器削弱或消除了傅里叶空间的高频分量，所以达到了增强低频分量，平滑图像中细节的效果；锐化滤波器削弱或消除了傅里叶空间的低频分量，所以达到了增强高频分量，锐化图像中的细节的效果。

联系：两者效果相反，互为补充，从原始图像减去平滑滤波器的效果得到锐化滤波器的效果，如果减去锐化滤波器的效果则得到平滑滤波器的效果。

1.6参考教材公式(11.

2.1)，如果对数底数不为2，当对数底数为e时，得到的信息单位为奈特(nat)，当底数为10时，得到的信息单位成为哈特利(Hartley)，推导他们与比特的换算关系。

答：用不同的单位要表达相同数量的意思，根据定义得：

1bit =-log2P(E),

1nat = -lnP(E)= =-log2ebit

1hartley = -lgP(E)==-log210bit

2.1编程实现教材第11章习题11-13和11-16题。

答：11-13

function [output,table] = lzw2norm(vector)

%LZW2NORM LZW Data Compression (decoder)

% For vectors, LZW2NORM(X) is the uncompressed vector of X using the LZW algorithm.

% [...,T] = LZW2NORM(X) returns also the table that the algorithm produces.

%

% For matrices, X(:) is used as input.

%

% Input must be of uint16 type, while the output is a uint8.

% Table is a cell array, each element containig the corresponding code. %

% This is an implementation of the algorithm presented in the article % https://www.360docs.net/doc/609456718.html,/markn/articles/lzw/lzw.htm

%

% See also NORM2LZW

% $Author: Giuseppe Ridino' $

% $Revision: 1.0 $ $Date: 10-May-2004 14:16:08 $

% How it decodes:

%

% Read OLD_CODE

% output OLD_CODE

% CHARACTER = OLD_CODE

% WHILE there are still input characters DO

% Read NEW_CODE

% IF NEW_CODE is not in the translation table THEN

% STRING = get translation of OLD_CODE

% STRING = STRING+CHARACTER

% ELSE

% STRING = get translation of NEW_CODE

% END of IF

% output STRING

% CHARACTER = first character in STRING

% add translation of OLD_CODE + CHARACTER to the translation table % OLD_CODE = NEW_CODE

% END of WHILE

% ensure to handle uint8 input vector

if ~isa(vector,'uint16'),

error('input argument must be a uint16 vector')

end

% vector as a row

vector = vector(:)';

% initialize table (don't use cellstr because char(10) will be turned to empty!!!)

table = cell(1,256);

for index = 1:256,

table{index} = uint16(index-1);

end

% initialize output

output = uint8([]);

code = vector(1);

output(end+1) = code;

character = code;

for index=2:length(vector),

element = vector(index);

if (double(element)+1)>length(table),

% add it to the table

string = table{double(code)+1};

string = [string character];

else,

string = table{double(element)+1};

end

output = [output string];

character = string(1);

[table,code] = addcode(table,[table{double(code)+1} character]);

code = element;

end

% ###############################################

function code = getcodefor(substr,table)

code = uint16([]);

if length(substr)==1,

code = substr;

else, % this is to skip the first 256 known positions

for index=257:length(table),

if isequal(substr,table{index}),

code = uint16(index-1); % start from 0

break

end

end

end

% ###############################################

function [table,code] = addcode(table,substr)

code = length(table)+1; % start from 1

table{code} = substr;

code = uint16(code-1); % start from 0

11-16

function [output,table] = norm2lzw(vector)

%NORM2LZW LZW Data Compression (encoder)

% For vectors, NORM2LZW(X) is the compressed vector of X using the LZW algorithm.

% [...,T] = NORM2LZW(X) returns also the table that the algorithm produces.

%

% For matrices, X(:) is used as input.

%

% Input must be of uint8 type, while the output is a uint16.

% Table is a cell array, each element containig the corresponding code. %

% This is an implementation of the algorithm presented in the article % https://www.360docs.net/doc/609456718.html,/markn/articles/lzw/lzw.htm

%

% See also LZW2NORM

% $Author: Giuseppe Ridino' $

% $Revision: 1.0 $ $Date: 10-May-2004 14:16:08 $

% How it encodes:

%

% STRING = get input character

% WHILE there are still input characters DO

% CHARACTER = get input character

% IF STRING+CHARACTER is in the string table then

% STRING = STRING+character

% ELSE

% output the code for STRING

% add STRING+CHARACTER to the string table

% STRING = CHARACTER

% END of IF

% END of WHILE

% output the code for STRING

% ensure to handle uint8 input vector

if ~isa(vector,'uint8'),

error('input argument must be a uint8 vector')

end

% vector as uint16 row

vector = uint16(vector(:)');

% initialize table (don't use cellstr because char(10) will be turned to empty!!!)

table = cell(1,256);

for index = 1:256,

table{index} = uint16(index-1);

end

% initialize output

output = vector;

% main loop

outputindex = 1;

startindex = 1;

for index=2:length(vector),

element = vector(index);

substr = vector(startindex:(index-1));

code = getcodefor([substr element],table);

if isempty(code),

% add it to the table

output(outputindex) = getcodefor(substr,table);

[table,code] = addcode(table,[substr element]);

outputindex = outputindex+1;

startindex = index;

else,

% go on looping

end

end

substr = vector(startindex:index);

output(outputindex) = getcodefor(substr,table);

% remove not used positions

output((outputindex+1):end) = [];

% ############################################### function code = getcodefor(substr,table)

code = uint16([]);

if length(substr)==1,

code = substr;

else, % this is to skip the first 256 known positions for index=257:length(table),

if isequal(substr,table{index}),

code = uint16(index-1); % start from 0

break

end

end

end

% ###############################################

function [table,code] = addcode(table,substr)

code = length(table)+1; % start from 1

table{code} = substr;

code = uint16(code-1); % start from 0

2.2对“Fig0419(a)(chestXray_original).tif”的胸透图像进行空域、频域或混合方法增强，写出源代码并进行逐行注释，用两种以上的方法。

%% fftshift 对数变换,所应用的图片本身很简单，就只有黑白2种颜色

clc

clear

f = imread('.\images\Fig0419(a)(chestXray_original).tif \Fig0403(a) (image).tif');

imshow(f)

title('原始图像')

imfinfo('.\images\dipum_images_ch04\Fig0419(a)(chestXray_origin al).tif. ');%此处如果用Imfinfo(f)就会报错fft

%没有居中的傅里叶频谱

F=fft2(f);%进行二维快速傅里叶变换，其结果和DFT的一样，只是计算机的计算速度变快了而已，因而叫fft

S=abs(F);%求傅里叶变换后的幅值

figure,subplot(121),imshow(S,[])

title('傅里叶频谱图像1');%title函数一定要放在坐标显示的下一句才有效。

subplot(122),imshow(S),title('傅里叶频谱图像2')

%居中的傅里叶频谱

Fc=fftshift(F);%将频谱图像原点移至图像矩形中间

S1=abs(Fc);

figure,subplot(121),imshow(S1,[]);%加了第二个参数后显示的图像正常

subplot(122),imshow(S1);%当没有第二个参数时，显示的图像为竖线加一些孤立的黑点，why？

%使用对数后视觉增强后的傅里叶频谱

S2=log(1+S1);

imshow(S2,[]);

%%用fftshift对数变换显示稍复杂的图片

clc

clear

f=imread('.\images\dipum_images_ch04\Fig0409(a)(bld).tif ');

imshow(f);

f=double(f);%其实这句不要试过对后面的变换结果也没有影响

F=fft2(f);

Fc=fftshift(F);

S=abs(Fc);

S2=log(1+S);%如果没有这句的话，那么根本看不到细节的图，所以一定要用对数压缩增加对比度

figure,imshow(S2,[])

%%试一下ifft功能

clc

clear

f=imread('.\images\dipum_images_ch04\Fig0419(a)(chestXray_origi

nal).tif ');

imshow(f);

f=double(f);

f(1:8,1:8)

F=fft2(f);

f=ifft2(F);%取傅里叶变换后的反傅里叶变换,直接是整数，并不需要像某些代码一样取real部分

f(1:8,1:8)

%%先0扩充再滤波

clc

clear

f=imread('.\images\dipum_images_ch04\Fig0419(a)(chestXray_origi nal).tif ');

imshow(f);

[M N]=size(f);

F=fft2(f);%没有经过0扩充直接计算fft

sig=10;%高斯滤波参数

H=lpfilter('gaussian',M,N,sig);

G=H.*F;%加了.号的乘法表示对应每个元素都相乘，没加.号的乘法说明是真正的矩阵乘法

g=real(ifft2(G));

figure,imshow(g,[]);

PQ=paddedsize(size(f));%PQ为0扩展的面积，这里设置为与图像的大小相同？这里size(f)的大小为256*256

Fp=fft2(f,PQ(1),PQ(2));%如果fft2有3个参数的话，则是进行了0扩展了

Hp=lpfilter('gaussian',PQ(1),PQ(2),2*sig);%构造高斯滤波器

Gp=Hp.*Fp;%对0扩展后的图像进行高斯滤波

gp=real(ifft2(Gp));

figure,imshow(gp,[])'

gpc=gp(1:size(f,1),1:size(f,2));%取gp图像的1到f的第一维的行，

1到f第二维的列部分，即取与图像大小相同的部分

figure,imshow(gpc,[]);%此时显示的结果应该与g图像一样。

%直接使用空间域滤波

h=fspecial('gaussian',15,7);

gs=imfilter(f,h);

figure,imshow(gs,[])

%%

clc

clear

f=imread('.\images\dipum_images_ch04\Fig0419(a)(chestXray_origi nal).tif ');

f=double(f);

imshow(f);

PQ=paddedsize(size(f));%size(f)的大小为256*256

sig=10;

H = lpfilter('gaussian',PQ(1),PQ(2),2*sig); % PQ=[512512] figure,mesh(abs(H(1:10:512,1:10:512)));%画出滤波器的三维空间图

g=dftfilt(f,H);

figure,imshow(g,[])

%%空间滤波和频域滤波的比较

clc

clear

f=imread('.\images\dipum_images_ch04\Fig0419(a)(chestXray_origi nal).tif ');

imshow(f);

F=fft2(f);

S=fftshift(log(1+abs(F)));

S=gscale(S);%用gscale来进行将数据缩放到默认值0~255

figure,imshow(S);

h=fspecial('sobel');%构造sobel空间滤波器

figure,freqz2(h);%查看滤波器的图形

PQ=paddedsize(size(f));

H=freqz2(h,PQ(1),PQ(2));%将空间滤波器h转换成频域滤波器H，但是滤波器的原点在矩阵的中心处

H1=ifftshift(H);%原点迁移到左上角

figure,mesh(abs(H1(1:20:1200,1:20:1200)));

figure,subplot(121),imshow(abs(H),[]);

subplot(122),imshow(abs(H1),[]);

gs=imfilter(double(f),h);%空间滤波，其实就是边缘检测

gf=dftfilt(f,H);%频域滤波

subplot(121),imshow(gs,[]),subplot(122),imshow(gf,[]);%将2种滤波结果显示出来

figure,imshow(abs(gs)>0.2*abs(max(gs(:))));%只显示最大值的2 0%的边缘

figure,imshow(abs(gf)>0.2*abs(max(gf(:))));

d=abs(gs-gf);

a=max(d(:))

b=min(d(:))

%%构造低通滤波器

clc

clear

f=imread('.\images\dipum_images_ch04\Fig0419(a)(chestXray_origi nal).tif ');

imshow(f)

F1=fft2(f);

imshow(log(1+abs(fftshift(F1))),[]);%直接显示取对数后的傅里叶变换

PQ=paddedsize(size(f));

[U V]=dftuv(PQ(1),PQ(2));%这个函数还真不明白什么意思，help中解释的是计算网格频率矩阵U和V，这2个矩阵是用来频率滤波的

D0=0.05*PQ(2);

F=fft2(f,PQ(1),PQ(2));%用0扩充大小为PQ(1)*PQ(2)，然后fft变换

figure,imshow(log(1+abs(fftshift(F))),[]);%显示0扩充后的ff t变换图

H=exp(-(U.^2+V.^2)/(2*(D0^2)));%构造频域滤波算子

figure,imshow(fftshift(H),[]);%显示平移后滤波器算子

g=dftfilt(f,H)

figure,imshow(g,[]);

%%绘制一个低通滤波器的mesh图

clc

clear

H1=lpfilter('gaussian',500,500,50)%构造一个中心点在左上角的高斯低通滤波器,size为500*500

H2=fftshift(H1);%将中心点平移至矩阵中心

mesh(H1(1:10:500,1:10:500)),figure,mesh(H2(1:10:500,1:10: 500));%绘出2者的mesh图

axis([05005001]);%xy轴范围0~50，z轴范围0~1

figure,subplot(121),imshow(H1,[]),subplot(122),imshow(H2, []);%绘出2者的灰度图

%%绘制一个低通滤波器的mesh图

clc

clear

H=fftshift(hpfilter('ideal',500,500,100));%构造理想的高通滤波器

mesh(H(1:10:500,1:10:500));

axis([15015001]);

colormap([000]);%mesh网格全部用黑色画

axis off%关掉坐标轴显示

grid off%关掉网格显示

figure,imshow(H,[])

%%使用高通滤波器

clc

clear

f=imread('.\images\dipum_images_ch04\Fig0419(a)(chestXray_origi nal).tif ');

imshow(f)

title('原始图像')

PQ=paddedsize(size(f));

D0=0.05*PQ(1);

H=hpfilter('gaussian',PQ(1),PQ(2),D0);

g=dftfilt(f,H);

figure,imshow(g,[])

%%将高通滤波和直方图均衡结合起来

clc

clear

f=imread('.\images\dipum_images_ch04\Fig0419(a)(chestXray_origi nal).tif ');

imshow(f)

title('原始图像')

PQ = paddedsize(size(f));

D0 = 0.05*PQ(1);

HBW=hpfilter('btw',PQ(1),PQ(2),D0,2);%构造高通Butterworth 滤波器，截断频率为D0

H=0.5+2*HBW;

gdw=dftfilt(f,HBW);

gbw=gscale(gdw);%将灰度值自动缩放到0~255之间

figure,subplot(121),imshow(gdw,[]);

title('高通滤波后图像');

gbe=histeq(gbw,256);%直方图均衡化

subplot(122),imshow(gbe,[]);

title('高通滤波且直方图均衡化后图像');

ghf=dftfilt(f,H);%H是HBW的线性变换

ghf=gscale(ghf);

figure,subplot(121),imshow(ghf,[]); title('高频强调滤波后图像');

ghe=histeq(ghf,256);

subplot(122),imshow(ghe,[]);

title('高频强调且均衡化后图像');

%%fftshift和ifftshift的加深理解

clc

clear

A=[2001

0000

0000

3004]

B=fftshift(A)%中心点平移

C=fftshift(fftshift(A))%还原成A

D=ifftshift(fftshift(A))%也同样还原成A E=ifftshift(A)%平移结果和B一样