一、实验目的

1. 掌握RGB彩色模型和HSI彩色模型之间的转换方法，以及彩色图像的平滑方法。
2. 掌握彩色图像边缘检测的方法。
3. 掌握一维和二维信号的单尺度和多尺度小波变换。
4. 掌握小波包分解的方法。

二、实验内容

1. 彩色图像平滑。(课本P310 例6.12)

图1显示了一幅RGB图像，请进行以下操作：
(1) 显示图1的红、绿、蓝分量图像。
(2) 显示图1的H、S、I分量图像。
(3) 用大小为 $5\times 5$ 的平均核分别对图1的R、G、B分量进行平滑，然后混合这些独立平滑后的图像以形成平滑后的全彩色图像。
(4) 保持图1的H分量和S分量不变，仅用大小为 $5\times 5$ 的平均核对I分量进行平滑，然后把最终的HSI转换成RGB图像以便显示。
(5) 求(3)和(4)中所得到的两幅图像的差值。
(6) 分析实验结果。

答：
(1)
(2)

%色调分量
numerator=0.5*((R-G)+(R-B));
denominator=sqrt((R-G).^2+(R-B).*(G-B));
theta=acos(numerator./(denominator+eps));

H = zeros(size(f(:,:,1)));
H(B > G) = 2*pi - theta(B > G);
H(B <= G) = theta(B <= G);
H=H/(2*pi);

%饱和度分量
min_val=min(R,min(G,B));
S=1-3./(R+G+B).*min_val;

%亮度分量
I=1/3*(R+G+B);
subplot(1,3,1);imshow(H);title('色调分量','FontSize',15);
subplot(1,3,2),imshow(S);title('饱和度分量','FontSize',15);
subplot(1,3,3),imshow(I);title('亮度分量','FontSize',15);

(3)

%用大小为5×5的平均核分别对图1的R、G、B分量进行平滑

% 定义平均核
kernel = ones(5, 5) / 25;

% 对每个通道应用平滑核
red_channel = imfilter(R, kernel, 'replicate');
green_channel = imfilter(G, kernel, 'replicate');
blue_channel = imfilter(B, kernel, 'replicate');

% 合并三个通道
smoothed_img = cat(3, red_channel, green_channel, blue_channel);

% 显示结果
imshow(smoothed_img);

(4)

(5)

%没有定标时，图片接近0，不是特别清晰
smoothed_img=im2uint8(smoothed_img);
diff=smoothed_img-rgb;
subplot(1,2,1);imshow(diff);title('未调整I分量之后的diff图','FontSize',13);

%将图片重新定标，使得图像更加清晰

% 转换数据类型为double
diff = im2double(diff);
% 仅对I通道执行调整
hsi=rgb_to_hsi(diff);
I=hsi(:,:,3);
I = imadjust(I);
hsi(:,:,3)=I;

diff=hsi_to_rgb(hsi);
subplot(1,2,2);imshow(diff);title('调制I分量之后的diff图','FontSize',13);


function rgb = hsi_to_rgb(hsi)
hsi=im2double(hsi);
H=hsi(:,:,1)*2*pi;
S=hsi(:,:,2);
I=hsi(:,:,3);
R=zeros(size(hsi,1),size(hsi,2));
G=zeros(size(hsi,1),size(hsi,2));
B=zeros(size(hsi,1),size(hsi,2));
%RG 扇区 (0<=H<2*pi/3)
c=find((0<=H) & (H<2*pi/3));
B(c)=I(c).*(1-S(c));
R(c)=I(c).*(1+S(c).*cos(H(c))./cos((pi/3)-H(c)));
G(c)=3*I(c)-(R(c)+B(c));
%BG 扇区 (2*pi/3<=H<4*pi/3)
c=find((2*pi/3<=H) & (H<4*pi/3));
R(c)=I(c).*(1-S(c));
G(c)=I(c).*(1+S(c).*cos(H(c)-2*pi/3)./cos(pi-H(c)));
B(c)=3*I(c)-(R(c)+G(c));
%BR 扇区 (4*pi/3<=H<2*pi)
c=find((4*pi/3<=H) & (H<2*pi));
G(c)=I(c).*(1-S(c));
B(c)=I(c).*(1+S(c).*cos(H(c)-4*pi/3)./cos((5*pi/3)-H(c)));
R(c)=3*I(c)-(B(c)+G(c));

rgb=cat(3,R,G,B);
rgb=im2uint8(rgb);
end

function hsi=rgb_to_hsi(rgb)
R=rgb(:,:,1);
G=rgb(:,:,2);
B=rgb(:,:,3);

%色调分量
numerator=0.5*((R-G)+(R-B));
denominator=sqrt((R-G).^2+(R-B).*(G-B));
theta=acos(numerator./(denominator+eps));

H = zeros(size(rgb(:,:,1)));
H(B > G) = 2*pi - theta(B > G);
H(B <= G) = theta(B <= G);
H=H/(2*pi);

%饱和度分量
min_val=min(R,min(G,B));
S=1-3./(R+G+B).*min_val;

%亮度分量
I=1/3*(R+G+B);


hsi=cat(3,H,S,I);
hsi=im2uint8(hsi);
end

(6)
可以看到，在RGB色彩模型中对彩色图像做平滑滤波，每个像素的颜色都是邻域中像素的平均颜色。而HSI彩色模型中仅对亮度分量做平均，解除了亮度和色彩信息的联系。两种平滑方法的区别会随着核的增大而变得更为明显。

2. 彩色边缘检测。(课本P318 例6.16)

(1) 请在RGB向量空间中计算图1的梯度。
(2) 请求图1的3幅分量图像的梯度，并通过在每一坐标点 处叠加相应的3幅分量图像的梯度值，合成一幅梯度图像。
(3) 请求(1)和(2)中所得的两幅梯度图像的差的绝对值。

答：

close all; clear all; clc;
f=imread('1.Lena.tif');
f=im2double(f);
% Sobel 算子
sobel = fspecial('sobel');
sv = sobel';
sh = sobel; 
% 计算 X、Y 方向梯度
Rx = imfilter(f, sv, 'replicate');
Ry = imfilter(f, sh, 'replicate');
Gx = imfilter(f, sv', 'replicate');
Gy = imfilter(f, sh', 'replicate');
Bx = imfilter(f, sv, 'replicate') + imfilter(f, sv', 'replicate');
By = imfilter(f, sh, 'replicate') + imfilter(f, sh', 'replicate');
% 计算 gxx、gyy 和 gxy
gxx = sum(cat(3, Rx, Gy, Bx).^2, 3);
gyy = sum(cat(3, Ry, Gy, By).^2, 3);
gxy = sum(cat(3, Rx, Gy, Bx) .* cat(3, Ry, Gy, By), 3);
% 计算 theta
theta = 0.5 * atan(2*gxy./(gxx - gyy + eps));
%计算 F
F1=sqrt(1/2*((gxx+gyy)+(gxx-gyy).*cos(2*theta)+2*gxy.*sin(2*theta)));
theta1=theta+pi/2;
F2=sqrt(1/2*((gxx+gyy)+(gxx-gyy).*cos(2*theta1)+2*gxy.*sin(2*theta1)));
F = max(F1, F2);
F = mat2gray(F);
figure; imshow(F);
imwrite(F,'边缘检测图像.tif');

(2)

%%
figure;
%R、 G、 B 分量梯度图像
RG=sqrt(Rx.^2+Ry.^2);
subplot(131); imshow(RG(:,:,1),[]);title('R分量图像梯度');
GG=sqrt(Gx.^2+Gy.^2);
subplot(132); imshow(GG(:,:,2),[]);title('G分量图像梯度');
BG=sqrt(Bx.^2+By.^2);
subplot(133); imshow(BG(:,:,3),[]);title('B分量图像梯度');
%叠加成一幅梯度图像
F3=RG + GG + BG;
F3=mat2gray(F3);
figure; imshow(F3);
imwrite(F3,'3 幅分量图像叠加的边缘检测图像.tif');

(3)

%两幅梯度图像做差
Fd=F-F3;
Fd=abs(Fd);
Fd=mat2gray(Fd);
figure; imshow(Fd);

3. 一维小波变换。(课本P364 例7.20)

有一离散函数 f ( x ) = { 1 , 4 , − 3 , 0 } f(x)=\{1,4,-3,0\} f(x)={1,4,−3,0} 请根据式(1)所示的哈尔尺度函数系数和式(2)所示的小波函数，请用实验讲义中的两种方法对 进行二尺度小波变换。

$$h_{\phi }(n)=\{\begin{array}{lc}1/\sqrt{2},& n=0,1\\ 0,& \text{ 其他 }\end{array}\text{\hspace{1em}(1)}$$

$$h_{\psi }(n)=\{\begin{array}{lc}1/\sqrt{2},& n=0\\ -1/\sqrt{2},& n=1\\ 0,& \text{ 其他 }\end{array}\text{\hspace{1em}(2)}$$

答：

方法一：

fn = [1,4,-3,0];
h_phi = [1/sqrt(2),1/sqrt(2)];
h_psi = [-1/sqrt(2),1/sqrt(2)];

% 一级小波分解
wf1 = conv(fn,h_psi);%高通细节分量
wf1 = wf1(2:2:end-mod(length(wf1),2));
sca1 = conv(fn,h_phi);%低通尺度分量
sca1 = sca1(2:2:end-mod(length(sca1),2))

% 二级小波分解,对尺度空间V_{J-1}进行分解
wf2 = conv(sca1,h_psi);
wf2 = wf2(2:2:end-mod(length(wf11),2));
sca2 = conv(sca1,h_phi);
sca2 = sca2(2:2:end-mod(length(sca2),2));

方法二：

fn=[1,4,-3,0];
[ca,cd] = dwt(fn,'haar');
[ca1, cd1]=dwt(ca,'haar')

4. 二维小波变换。(课本P369 例7.22)

图2给出了一幅512×512的花瓶和窗台图像，请对其进行以下操作：

(1) 请求图2的一尺度小波分解后图像。
(2) 请求图2的二尺度小波分解后图像。

答：
(1)

f=imread('4.一幅512×512的图像.tif'); 
f=im2double(f); 
[cA,cH,cV,cD]=dwt2(f, 'haar'); 
cA=mat2gray(cA); 
cH=mat2gray(cH); 
cV=mat2gray(cV); 
cD=mat2gray(cD); 
w=[cA,cH;cV,cD];
figure; imshow(w); 

(2)

[cA2,cH2,cV2,cD2]=dwt2(cA, 'haar'); 
cA2=mat2gray(cA2); 
cH2=mat2gray(cH2); 
cV2=mat2gray(cV2); 
cD2=mat2gray(cD2);
w1=[cA2,cH2;cV2,cD2];
w2=[w1,cH;cV,cD];
figure;imshow(w2);

5. 小波包分解。(课本P376 例7.24)

图3给出了一幅扫描的指纹图像，请对其进行二尺度、全小波包分解。

答：

clear all;close all;clc;
f = imread('5.扫描的指纹图像.tif');
f=im2double(f);
[cA1,cH1,cV1,cD1]=dwt2(f, 'haar');
cA1=mat2gray(cA1);
cH1=mat2gray(cH1);
cV1=mat2gray(cV1);
cD1=mat2gray(cD1);

[cA2,cH2,cV2,cD2]=dwt2(cA1,'haar');
cA2=mat2gray(cA2);
cH2=mat2gray(cH2);
cV2=mat2gray(cV2);
cD2=mat2gray(cD2);
w2=[cA2,cH2;cV2,cD2];

[cA2,cH2,cV2,cD2]=dwt2(cH1,'haar');
cA3=mat2gray(cA2);
cH3=mat2gray(cH2);
cV3=mat2gray(cV2);
cD3=mat2gray(cD2);
w3=[cA3,cH3;cV3,cD3];

[cA2,cH2,cV2,cD2]=dwt2(cV1,'haar');
cA4=mat2gray(cA2);
cH4=mat2gray(cH2);
cV4=mat2gray(cV2);
cD4=mat2gray(cD2);
w4=[cA4,cH4;cV4,cD4];

[cA2,cH2,cV2,cD2]=dwt2(cD1,'haar');
cA5=mat2gray(cA2);
cH5=mat2gray(cH2);
cV5=mat2gray(cV2);
cD5=mat2gray(cD2);
w5=[cA5,cH5;cV5,cD5];

%二尺度分解
wt1=[w2,cH1;cV1,cD1];
figure;imshow(wt1);
%全分解
wall=[w2,w3;w4,w5];
figure;imshow(wall)