彩色图像的空间域滤波

（1） RGB彩色空间向HSI彩色空间的转换：

自定义一个函数，实现RGB彩色空间向HSI彩色空间的转换，要求该函数的输入参数为RGB彩色图像，输出参数为HSI彩色图像。

根据RGB彩色空间到HSI彩色空间的转换公式，编写函数 RGBtoHSI(pho)

% pho 表示输入图像路径（包括图像名+后缀）
function hsi = RGBtoHSI(pho)
ima=imread(pho); %读取输入图像
rgb=double(ima); %将图像转为双精度类型
rgb=rgb/255; %归一化
R=rgb(:,:,1); %提取第1个通道
G=rgb(:,:,2); %提取第2个通道
B=rgb(:,:,3); %提取第3个通道
 
angle=(acos(1/2*((R-G)+(R-B))./sqrt((R-G).^2+(R-B).*(G-B))))*180/pi;
H=angle; %计算H（色调）通道
H(B>G)=360-H(B>G); %将 B>G 对应位置上的值赋值为 2*pi - angle
S=1-3./(R+G+B+eps).*min(min(R,G),B); %计算S（色饱和度）通道
I=1/3*(R+G+B+eps); %计算I（亮度）通道
 
hsi=cat(3,H/360,S,I); %合成3个通道，为HSI图像

输入原图像并进行显示，将原图像从RGB彩色空间向HSI彩色空间转换，显示HSI彩色空间的图像 useRGBtoHSI('EXP6.tif')

% pho 表示输入图像路径（包括图像名+后缀）
function [] = useRGBtoHSI(pho)
ima=imread(pho); %读入输入图像
subplot(121); %设置输入图像的显示位置（1行2列的第1个位置）
imshow(ima); %显示输入图像
title('RGB彩色图像'); %设置输入图像的标题
hsi=RGBtoHSI(pho); %将输入图像由RGB空间向HSI空间转化
subplot(122); %设置输出图像的显示位置（1行2列的第2个位置）
imshow(hsi); %显示输出图像
title('HSI彩色图像'); %设置输出图像的标题
%imwrite(hsi,'EXP6_HSI.tif','tif');

输出结果如下：

注意：

在编写程序的过程中运行程序曾出现如下结果：

解决：

HSI彩色图像输出结果与正确结果的颜色不一样。思考发现，H通道处并未做归一化处理，即 H=H/360才是正确值。

（2）HSI彩色空间向RGB彩色空间的转换：

自定义一个函数，实现HSI彩色空间向RGB彩色空间的转换，要求该函数的输入参数为HSI彩色图像，输出参数为RGB彩色图像。

根据HSI彩色空间到RGB彩色空间的转换公式，编写函数 HSItoRGB(hsi)

% hsi 表示HSI空间图像
function rgb = HSItoRGB(hsi)
%ima=imread(pho);
hsi=double(hsi); %转化为双精度类型
[r c m]=size(hsi); %计算hsi的行列和维度
H=hsi(:,:,1); %提取H分量
S=hsi(:,:,2); %提取S分量
I=hsi(:,:,3); %提取I分量
H=H*360; %将色调值变为原来的范围[0°,360°]
R=zeros(r,c); %R分量
G=zeros(r,c); %G分量
B=zeros(r,c); %B分量
 
inv=0<=H & H<120; %RG扇区（inv是H中值在RG扇区的位置）
B(inv)=I(inv).*(1-S(inv)); %根据在RG扇区时的计算公式计算B分量
R(inv)=I(inv).*(1+S(inv).*cosd(H(inv))./cosd(60-H(inv))); %根据在RG扇区时的计算公式计算R分量
G(inv)=3*I(inv)-(R(inv)+B(inv)); %根据在RG扇区时的计算公式计算G分量
 
inv=120<=H & H<240; %GB扇区（inv是H中值在GB扇区的位置）
H(inv)=H(inv)-120;
R(inv)=I(inv).*(1-S(inv)); %根据在GB扇区时的计算公式计算R分量
G(inv)=I(inv).*(1+S(inv).*cosd(H(inv))./cosd(60-H(inv))); %根据在GB扇区时的计算公式计算G分量
B(inv)=3*I(inv)-(R(inv)+G(inv)); %根据在GB扇区时的计算公式计算B分量
 
inv=240<=H & H<360; %BR扇区（inv是H中值在BR扇区的位置）
H(inv)=H(inv)-240;
G(inv)=I(inv).*(1-S(inv)); %根据在BR扇区时的计算公式计算G分量
B(inv)=I(inv).*(1+S(inv).*cosd(H(inv))./cosd(60-H(inv))); %根据在BR扇区时的计算公式计算B分量
R(inv)=3*I(inv)-(G(inv)+B(inv)); %根据在BR扇区时的计算公式计算R分量
 
%将各个分量的灰度级拉伸到256个灰度级上
R=R*255;
G=G*255;
B=B*255;
rgb=cat(3,R,G,B); %合成彩色图像
rgb=uint8(rgb); %转为8位无符号整型

输入原图像，将原图像从RGB彩色空间向HSI彩色空间转换，显示HSI彩色空间的图像，并将该图像向RGB彩色空间转化，得到结果进行显示  useHSItoRGB('EXP6.tif')

% pho 表示输入图像路径（包括图像名+后缀）
function [] = useHSItoRGB(pho)
hsi=RGBtoHSI(pho); %将输入图像转化为HSI空间彩色图像
subplot(121); %设置HSI空间彩色图像的显示位置（1行2列的第1个位置）
imshow(hsi); %显示HSI空间彩色图像
title('HSI彩色图像'); %设置HSI空间彩色图像的标题
rgb=HSItoRGB(hsi); %由HSI空间转化为RGB空间
subplot(122); %设置输出图像的显示位置（1行2列的第2个位置）
imshow(rgb); %显示输出图像
title('RGB彩色图像'); %设置输出图像的标题

输出结果如下：

注意：

在编写程序的过程中运行程序曾出现如下结果：

不仅如此，其中还出现一些“奇特”的图片，观察原图像转为HSI彩色图像后中H通道的数值，再看看作为实验2小题输入的HSI彩色图像的H通道的数值，就可以发现这两个通道的数值本应该一样却不一样了。我能想到的原因应该是作为实验2小题的输入图像——HSI彩色图像是我由RGB转HSI得到的（输出保存），其中保存的过程中HSI彩色图像H通道内的值转化为了256灰度级来呈现（保存）。所以第2小题的输入图像不能由第1小题的输入结果保存得到，而是直接使用第1小题的结果作为输入参数。

解决方法：对原程序进行修改——输入参数改为转化后的HSI。

还出现过如下结果：

该结果是经上述问题后修改的程序运行的结果，再次查看原程序，会发现H还未乘以360。

解决方法：对原程序进行修改——H分量乘以360。

接下来说的是：原图像与由原图像转化HSI彩色图像再转为RGB彩色图像的对比，查看两者差值，判断RGB——>HSI和HSI——>RGB两个函数的正确性。

errorAnaly('EXP6.tif')

% pho 表示输入图像路径（包括图像名+后缀）
% 该函数显示输入图像和其HSI空间图像、由HSI空间转化来的RGB图像、以及原图像与由HSI空间转化来的RGB图像的差值图像
% 观察在转化过程是否出现过大误差
function []=errorAnaly(pho)
ima=imread(pho); %读入输入图像
subplot(221); %设置输入图像的显示位置（2行2列的第1个位置）
imshow(ima); %显示输入图像
title('原图像'); %设置输入图像的标题
 
hsi=RGBtoHSI(pho); %将输入图像转化为HSI空间彩色图像
subplot(222); %设置HSI空间彩色图像的显示位置（2行2列的第2个位置）
imshow(hsi); %显示HSI空间彩色图像
title('HSI彩色图像'); %设置HSI空间彩色图像的标题
 
rgb=HSItoRGB(hsi); %由HSI空间转化为RGB空间
subplot(223); %设置输出图像的显示位置（2行2列的第3个位置）
imshow(rgb); %显示输出图像
title('RGB彩色图像'); %设置输出图像的标题
 
erra=ima-rgb; %原图像与由HSI空间转化来的RGB图像的差值图像
subplot(224); %设置差值图像的显示位置（2行2列的第4个位置）
imshow(erra); %显示差值图像
title('原图像与RGB彩色图像的差值'); %设置差值图像的标题

输出结果如下：

由差值图像可见，原图像与RGB彩色图像的“一致性”，说明两个函数的正确性。

（3）要求在HSI彩色空间中，实现对含噪图像EXP6.tif的空间域平滑处理，滤波器为5×5高斯均值滤波器（滤波器模板如下图所示），并在RGB彩色空间中显示处理前后的结果。

输入图像在HSI彩色空间中的空间域平滑处理：

使用上次我们介绍到的自定义的空间滤波器（ http://blog.csdn.net/qq_15096707/article/details/50061003 ），加入5x5高斯均值滤波器，对原函数进行修改补充，如下：

%filter_type 所选择的滤波器类型。有：3x3高斯均值滤波器、5x5高斯均值滤波器、中值滤波器、最大值滤波器、最小值滤波器
%filter_size 定义滤波器大小，如 filter_size=[3 3];
function ima2 = mySpatialFilter(ima, filter_type, filter_size)
%ima=imread(pho); %根据路径读取原图像
ima=double(ima); %转化为双精度类型
[r c]=size(ima); %读取原图像的行数和列数
 
%设置高斯均值滤波器只有3x3或5x5大小的
if strcmp(filter_type, '3x3高斯均值滤波器')
    filter_size=[3 3];
end
 
if strcmp(filter_type, '5x5高斯均值滤波器')
    filter_size=[5 5];
end
m=filter_size(1); %读取设置的滤波器模板的大小
n=filter_size(2);
 
rs=r+m-1; %计算边界填充后的图像大小（填充模板大小）
cs=c+n-1;
tIma=zeros(rs,cs); %创建一个边界填充后的矩阵
tIma((m-1)/2+1:rs-(m-1)/2, (n-1)/2+1:cs-(n-1)/2)=ima; %矩阵中间填充原图像灰度级
 
%以最邻近的边界作为边界填充的灰度值
%先填充上边界
for i=1:(m-1)/2
    tIma(i,(n-1)/2+1:cs-(n-1)/2)=ima(1);
end
%填充下边界
for i=rs-(m-1)/2+1:rs
    tIma(i,(n-1)/2+1:cs-(n-1)/2)=ima(r);
end
%填充左边界
for i=1:(n-1)/2
    tIma(:,i)=tIma(:,(n-1)/2+1);
end
%填充右边界
for i=cs-(n-1)/2+1:cs
    tIma(:,i)=tIma(:,cs-(n-1)/2);
end
 
rIma=zeros(rs,cs);
%遍历该矩阵，套用滤波模板
for i=1:rs-m+1
    for j=1:cs-n+1
        %读取模板覆盖处的图像
        tempIma=tIma(i:i+m-1, j:j+n-1);
        %计算该模板覆盖的中心位置
        centerR=i+(m-1)/2;
        centerC=j+(m-1)/2;
        switch filter_type
            case '3x3高斯均值滤波器'
                %定义3 x 3的高斯均值滤波器模板
                module=[1 2 1; 2 4 2; 1 2 1];
                module=module(:)';
                %为邻域以该模板为系数进行加权求平均，得出的值赋值得邻域中心
                rIma(centerR,centerC)=module*tempIma(:)/sum(module);
            case '5x5高斯均值滤波器'
                %定义5 x 5的高斯均值滤波器模板
                module=[1 2 3 2 1; 2 5 6 5 2; 3 6 8 6 3; 2 5 6 5 2; 1 2 3 2 1];
                module=module(:)';
                %为邻域以该模板为系数进行加权求平均，得出的值赋值得邻域中心
                rIma(centerR,centerC)=module*tempIma(:)/sum(module);
            case '中值滤波器'
                %为邻域内的值排序，求得中值，将该值赋值给邻域中心
                tempIma=sort(tempIma(:));
                rIma(centerR,centerC)=tempIma((length(tempIma)+1)/2,1);
            case '最大值滤波器'
                %求邻域中的最大灰度值，将该值赋值给邻域中心
                rIma(centerR,centerC)=max(tempIma(:));
            case '最小值滤波器'
                %求邻域中的最小灰度值，将该值赋值给邻域中心
                rIma(centerR,centerC)=min(tempIma(:));
            otherwise
                error('不存在该滤波器');
        end
    end
end
%去除原先填充的边界，得出最终结果
ima2=rIma((m-1)/2+1:rs-(m-1)/2, (n-1)/2+1:cs-(n-1)/2);
ima2=uint8(ima2);

将原图像转为HSI空间图像，在HSI空间中对I分量使用5x5高斯均值滤波器，再将HSI空间向RGB空间转化，显示结果：

useHSIGS('EXP6.tif')

% pho 表示输入图像路径（包括图像名+后缀）
% 该函数输入原图像，将原图像转化为HSI彩色图像，在HSI空间使用5x5高斯均值滤波器，再将结果转为RGB输出
function []=useHSIGS(pho)
ima=imread(pho); %读取输入图像
subplot(121); %设置输入图像的显示位置（1行2列的第1个位置）
imshow(ima); %显示输入图像
title('输入图像'); %设置输入图像的标题
 
hsi=RGBtoHSI(pho); %将输入图像由RGB空间向HSI空间转化
%在使用5x5高斯均值滤波器前，需要将I分量拉伸到[0 255]灰度区间
%为什么只对I分量进行操作呢？I分量与图像的彩色信息无关。
%注意差异性：
%RGB空间中的平均是不同彩色的平均
%HIS空间中仅仅是强度的平均，色调H和饱和度S均保持不变
hsi(:,:,3)=mySpatialFilter(hsi(:,:,3)*255,'5x5高斯均值滤波器',[5 5]);
hsi(:,:,3)=hsi(:,:,3)/255; %将I分量再转回[0 1]区间
rgb=HSItoRGB(hsi); %由HSI空间再转为RGB空间
 
subplot(122); %显示输出图像
imshow(rgb); %显示输出图像
title('使用5x5高斯均值滤波器的输出图像'); %设置输出图像的标题

输出结果如下：

注意：

为什么在HSI空间使用高斯均值滤波器时只对I分量进行处理？

I分量与图像的彩色信息无关，H&S分量与人感受颜色的方式紧密相连。

注意差异性：

RGB空间中的锐化是不同彩色的平均；

HSI空间中仅仅是强度的锐化，色调H和饱和度S均保持不变。