8.图像变换技术

经过图像变换后，一方面能够更有效地反映图像自身的特征，另一方面也可使能量集中在少量数据上，更有利于图像的存储、传输和处理。

8.1 图像Radon变换

从检测器获取投影数据的过程，就是图像中的Radon变换。

8.1.1 Radon正变换

%对图像进行0°和45°方向上的Radon变换
clear all; close all;
I=zeros(200, 200); %建立图像
I(50:150, 50:150)=1;
[R, xp]=radon(I, [0, 45]); %Radon变换
figure;
subplot(131);imshow(I);
subplot(132);plot(xp, R(:, 1)); %0° Radon变换结果
subplot(133);plot(xp, R(:, 2)); %45° Radon变换结果

%对图像从0°到180°每隔10°做Radon变换
clear all; close all;
I=zeros(200, 200);
I(50:150, 50:150)=1;
theta=0:10:180; %角度值
[R, xp]=radon(I, theta); %Radon变换
figure;
subplot(121);imshow(I);
subplot(122);imagesc(theta, xp, R); %绘制各个角度变换结果
colormap(hot); %设置调色板
colorbar; %添加颜色条

%Radon变换检测直线
clear all; close all;
I=fitsread('solarspectra.fts');
J=mat2gray(I); %转换为灰度图像
BW=edge(J); %获取边缘
figure;subplot(121);imshow(J);
subplot(122);imshow(BW);
theta=0:179; %角度
[R, xp]=radon(BW, theta); %Radon变换
figure;imagesc(theta, xp, R); %绘制各个角度变换结果
colormap(hot); %设置调色板
colorbar; %添加颜色条
Rmax=max(max(R)) %获取最大值
[row, column]=find(R>=Rmax) %获取行和列值
x=xp(row) %获取位置
angel=theta(column) %获取角度

8.1.2 Radon反变换

%Radon反变换恢复图像
clear all; close all;
I=imread('circuit.tif');
theta=0:2:179;
[R, xp]=radon(I, theta); %Radon变换
J=iradon(R, theta); %Radon反变换
figure;
subplot(131);imshow(uint8(I));
subplot(132);imagesc(theta, xp, R);
axis normal;
subplot(133);imshow(uint8(J));

iradon()函数利用R矩阵各列的投影来构造图像I的近似值。投影数越多，重建后的图像越接近原始的图像。为了获取准确的图像，尽量使用较多的投影值。

%投影角度的多少对Radon变换和Radon反变换的影响
clear all; close all;
I=phantom(256); %读入图像
figure;
imshow(I);
theta1=0:10:170;
theta2=0:5:175;
theta3=0:2:178;
[R1, xp]=radon(I, theta1); 
[R2, xp]=radon(I, theta2);
[R3, xp]=radon(I, theta3);
figure;
imagesc(theta3, xp, R3);
colormap hot;
colorbar;
J1=iradon(R1, 10);
J2=iradon(R2, 5);
J3=iradon(R3, 2);
figure;
subplot(131);imshow(J1);
subplot(132);imshow(J2);
subplot(133);imshow(J3);

8.2 图像傅里叶变换

数学意义：将一个图像转换为一系列周期函数来处理

物理意义：傅里叶变换将图像从空间域转换到频率域，逆变换将图像从频率域转换到空间域

实际上对图像进行二维傅里叶变换得到频谱图，就是图像梯度的分布图。频谱图上看到的明暗不一的亮点，对应图像上某一点与邻域点差异的强弱，即梯度大小。如果频谱图中暗点数多，实际图像比较柔和，反之，如果频谱图中亮点多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大。

8.2.1 傅里叶变换的定义与性质

一维连续傅里叶变换：

 

 离散傅里叶变换：

 

 二维离散傅里叶变换的性质

 

8.2.2 傅里叶变换的MATLAB实现

%矩阵的二维离散傅里叶变换
clear all; close all;
I1=ones(4) %建立矩阵
I2=[2 2 2 2;1 1 1 1; 3 3 0 0; 0 0 0 0] 
J1=fft2(I1) %傅里叶变换
J2=fft2(I2)

%图像的二维离散傅里叶变换
clear all; close all;
I=imread('cameraman.tif'); %读入图像
J=fft2(I); %傅里叶变换
K=abs(J/256);
figure;
subplot(121);imshow(I);
subplot(122);imshow(uint8(K));

%通过函数ffthsift()进行平移
clear all; close all;
N=0:4
X=fftshift(N) %平移
Y=fftshift(fftshift(N)) %再平移
Z=ifftshift(fftshift(N)) %反平移

%图像进行傅里叶变换和平移
%图像的能量主要集中在低频
clear all; close all;
I=imread('peppers.png');
J=rgb2gray(I);
K=fft2(J); %傅里叶变换
K=fftshift(K); %平移
L=abs(K/256);
figure;
subplot(121);imshow(J);
subplot(122);imshow(uint8(L));

%图像变亮后进行傅里叶变换
%中央低频部分变大了，频谱图的中央低频部分代表了灰度图像的平均亮度
clear all; close all;
I=imread('peppers.png');
J=rgb2gray(I);
J=J*exp(1); %变亮
J(find(J>255))=255;
K=fft2(J);
K=fftshift(K);
L=abs(K/256);
figure;
subplot(121);imshow(J);
subplot(122);imshow(uint8(L));

%图像旋转后进行傅里叶变换
clear all; close all;
I=imread('peppers.png'); 
J=rgb2gray(I);
J=imrotate(J,45,'bilinear'); %图像旋转
K=fft2(J); %傅里叶变换
K=fftshift(K); %旋转
L=abs(K/256); 
figure;
subplot(121);imshow(J);
subplot(122);imshow(uint8(L));

%图像中添加高斯噪声后进行傅里叶变换
clear all; close all;
I=imread('peppers.png');
J=rgb2gray(I);
J=imnoise(J, 'gaussian', 0, 0.01);
K=fft2(J);
K=fftshift(K);
L=abs(K/256);
figure;
subplot(121);imshow(J);
subplot(122);imshow(uint8(L));

%灰度图像的傅里叶变换和反变换
clear all; close all;
I=imread('onion.png');
J=rgb2gray(I);
K=fft2(J); %傅里叶变换
L=fftshift(K); %平移
M=ifft2(K); %傅里叶反变换
figure;
subplot(121);imshow(uint8(abs(L)/198)); %显示频谱图
subplot(122);imshow(uint8(M)); %显示反变换得到的图像

%灰度图像的复制谱和相位谱
clear all; close all;
I=imread('peppers.png');
J=rgb2gray(I);
K=fft2(J);
L=fftshift(K);
fftr=real(L);
ffti=imag(L);
A=sqrt(fftr.^2+ffti.^2); %幅度谱
A=(A-min(min(A)))/(max(max(A))-min(min(A)))*255; %归一化到0-255
B=angle(K); %相位谱
figure;
subplot(121);imshow(A); %显示幅度谱
subplot(122);imshow(real(B)); %显示相位谱

%编程实现二维离散傅里叶变换
clear all; close all;
I=imread('onion.png');
J=rgb2gray(I);
J=double(J);
s=size(J);
M=s(1); N=s(2); %获取图像的行数和列数
for u=0:M-1
    for v=0:N-1
        k=0;
        for x=0:M-1
            for y=0:N-1
                k=J(x+1, y+1)*exp(-j*2*pi*(u*x/M+v*y/N))+k; %二维离散傅里叶变换公式
            end
        end
        F(u+1, v+1)=k; %傅里叶变换结果
    end
end
K=fft2(J); %采用了快速傅里叶变换算法，运算速度较快
figure;
subplot(121);imshow(K);%显示结果
subplot(122);imshow(F);

8.2.3 傅里叶变换的应用

%傅里叶变换识别图像中的字符
clear all; close all;
I=imread('text.png'); %读入图像
a=I(32:45, 88:98); %模板
figure;imshow(I); %显示原图像
figure;imshow(a); %显示模板
c=real(ifft2(fft2(I).*fft2(rot90(a, 2), 256, 256))); %傅里叶变换
figure;imshow(c, []); 显示卷积后的结果
max(c(:)) %取最大值
thresh=60; %设置阈值
figure;imshow(c>thresh) %显示识别结果

巴特沃斯滤波器

%对图像进行巴特沃斯低通滤波
clear all; close all;
I=imread('cameraman.tif');  
I=im2double(I);          
J=fftshift(fft2(I)); %傅里叶变换与平移
[x, y]=meshgrid(-128:127, -128:127); %产生离散数据
z=sqrt(x.^2+y.^2);
D1=10;  D2=30; %滤波器的截至频率
n=6; %滤波器的阶数
H1=1./(1+(z/D1).^(2*n)); %滤波器
H2=1./(1+(z/D2).^(2*n)); %滤波器
K1=J.*H1; %滤波
K2=J.*H2; %滤波
L1=ifft2(ifftshift(K1)); %傅里叶反变换
L2=ifft2(ifftshift(K2)); %傅里叶反变换
figure;
subplot(131);imshow(I); %原图像
subplot(132);imshow(real(L1)); %结果图像
subplot(133);imshow(real(L2)); %截止频率月底，图像变得越模糊

%对图像进行巴特沃斯高通滤波
clear all; close all;
I=imread('cameraman.tif');  
I=im2double(I);          
J=fftshift(fft2(I));    
[x, y]=meshgrid(-128:127, -128:127);
z=sqrt(x.^2+y.^2);
D1=10;  D2=40; %截止频率
n1=4;  n2=8; %滤波器阶数
H1=1./(1+(D1./z).^(2*n1)); 
H2=1./(1+(D2./z).^(2*n2));
K1=J.*H1; %滤波
K2=J.*H2;
L1=ifft2(ifftshift(K1)); %傅里叶反变换
L2=ifft2(ifftshift(K2));
figure;
subplot(131);imshow(I); 
subplot(132);imshow(real(L1));
subplot(133);imshow(real(L2)); %高通滤波保持了图像的边缘信息

8.3 图像离散余弦变换

8.3.1 离散余弦变换的定义

 

8.3.2 离散余弦变换的MATLAB实现

%对图像进行二维离散余弦变换
clear all; close all;
I=imread('coins.png');
I=im2double(I);
J=dct2(I); %二维离散余弦变换
figure;
subplot(121);imshow(I);
subplot(122);imshow(log(abs(J)), []); %显示变换系数，系数中的能量主要集中在左上角，其余大部分系数接近0

 

%dctmtx()生成离散余弦变换矩阵
clear all; close all;
A=[1 1 1 1; 2 2 2 2; 3 3 3 3] %建立矩阵
s=size(A); 
M=s(1); %矩阵的行数
N=s(2); %矩阵的列数
P=dctmtx(M) %离散余弦变换矩阵
Q=dctmtx(N) %离散余弦变换矩阵
B=P*A*Q' %离散余弦变换

%利用dctmtx()进行图像的离散余弦变换
clear all; close all;
I=imread('cameraman.tif');
I=im2double(I);
s=size(I);
M=s(1);
N=s(2);
P=dctmtx(M);
Q=dctmtx(N);
J=P*I*Q'; %离散余弦变换
K=dct2(I); %离散余弦变换
E=J-K; %变换系数之差
find(abs(E)>0.000001) %查找系数差大于0.000001
figure;
subplot(121);imshow(J); %显示离散余弦系数
subplot(122);imshow(K); %显示离散余弦系数

%图像的二维离散余弦反变换
clear all; close all;
I=imread('cameraman.tif');
I=im2double(I);
J=dct2(I); %二维离散余弦变换
J(abs(J)<0.1)=0; %绝对值小于0.1的系数设置为0
K=idct2(J); %二维离散余弦反变换
figure;
subplot(131);imshow(I); %显示原图像
subplot(132);imshow(J); %变换系数
subplot(133);imshow(K); %显示结果图像

8.3.3 离散余弦变换的应用

%通过函数blkproc()对图像进行块操作
clear all; close all;
I=imread('cameraman.tif');
fun1=@dct2; %函数句柄
J1=blkproc(I, [8 8], fun1); %块操作
fun2=@(x) std2(x)*ones(size(x)); %函数句柄
J2=blkproc(I, [8 8], fun2); %块操作
figure;
subplot(121);imagesc(J1);
subplot(122);imagesc(J2);
colormap gray; %设置调色板

 

%通过离散余弦变换进行图像压缩
clear all; close all;
I=imread('rice.png'); %读入图像
J=im2double(I);
T=dctmtx(8); %计算离散余弦变换矩阵
K=blkproc(J, [8 8], 'P1*x*P2', T, T');  %对每个小方块进行离散余弦变换
mask=[ 1  1  1  1  0  0  0  0
            1  1  1  0  0  0  0  0
            1  1  0  0  0  0  0  0
            1  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0 ];
K2=blkproc(K, [8 8], 'P1.*x', mask); %系数选择
L=blkproc(K2, [8 8], 'P1*x*P2', T', T); %对每个小方块进行离散余弦反变换
figure;
subplot(121);imshow(J);
subplot(122);imshow(L);

8.4 其他图像变换

8.4.1 Hadamard变换

 

%通过函数hadamard()产生Hadamard变换矩阵
clear all; close all;
H1=hadamard(2) %产生2阶Hadamard变换矩阵
H2=hadamard(4) %产生4阶Hadamard变换矩阵
H3=H2'*H2

%对图像进行Hadamard变换
clear all; close all;
I=imread('peppers.png'); %读入rgb图像
I=rgb2gray(I); %转换为灰度图像
I=im2double(I);
h1=size(I, 1); %行
h2=size(I, 2); %列
H1=hadamard(h1); %Hadamard变换矩阵
H2=hadamard(h2); %Hadamard变换矩阵
J=H1*I*H2/sqrt(h1*h2); %Hadamard变换
figure;
set(0,'defaultFigurePosition',[100,100,1000,500]);
set(0,'defaultFigureColor',[1 1 1])
subplot(121);imshow(I); %原图像
subplot(122);imshow(J); %Hadamard变换结果

8.4.2 Hough变换

%对图像进行Hough变换
clear all; close all;
I=imread('circuit.tif');
I=im2double(I);
BW=edge(I, 'canny'); %边缘检测
[H, Theta, Rho]=hough(BW, 'RhoResolution', 0.5, 'ThetaResolution', 0.5); %Hough变换
figure;
set(0,'defaultFigurePosition',[100,100,1000,500]);
set(0,'defaultFigureColor',[1 1 1])
subplot(121);imshow(BW);
subplot(122);imshow(imadjust(mat2gray(H)));
axis normal; %设置坐标轴
hold on;
colormap hot; %设置调色板

%通过图像的Hough变换检测直线
clear all; close all;
I=imread('gantrycrane.png');
I=rgb2gray(I);
BW=edge(I, 'canny'); %获取图像边缘
[H, Theta, Rho]=hough(BW, 'RhoResolution', 0.5, 'Theta', -90:0.5:89.5);%Hough变换
P=houghpeaks(H, 5, 'threshold', ceil(0.3*max(H(:)))); %获取5个最值点
x=Theta(P(:, 2)); %横坐标
y=Rho(P(:, 1)); %纵坐标
figure;
set(0,'defaultFigurePosition',[100,100,1000,500]);
set(0,'defaultFigureColor',[1 1 1])
subplot(121);
imshow(imadjust(mat2gray(H)), 'XData', Theta, 'YData', Rho,...
    'InitialMagnification', 'fit'); %绘制Hough变换结果
axis on;  %设置坐标轴
axis normal;
hold on;
plot(x, y, 's', 'color', 'white');
lines=houghlines(BW, Theta, Rho, P, 'FillGap', 5, 'MinLength', 7); %检测直线
subplot(122);imshow(I);
hold on;
maxlen=0;
for k=1:length(lines) %绘制多条直线
    xy=[lines(k).point1; lines(k).point2];
    plot(xy(:,1), xy(:, 2), 'linewidth', 2, 'color', 'green');
    plot(xy(1,1), xy(1, 2), 'linewidth', 2,  'color', 'yellow');
    plot(xy(2,1), xy(2, 2), 'linewidth', 2,  'color', 'red');
    len=norm(lines(k).point1-lines(k).point2);
    if (len>maxlen) %获取最长直线坐标
        maxlen=len;
        xylong=xy;
    end
end
hold on;
plot(xylong(:, 1), xylong(:, 2), 'color', 'blue'); %绘制最长的直线