基于MATLAB的常见图像处理技术--图像变换技术

图像变换技术介绍

图像Radon变换

Radon变换介绍

Radon正变换

对图像进行0°和45°方向上的Radon变换

close all; clear all; clc;
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变换

close all; clear all; clc;
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);   %绘制各个角度的Radon变换结果
colormap(hot);           %设置调色板
colorbar;                %添加颜色条

通过radon变换来检测直线

close all; clear all; clc;
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)        %获取角度

Radon反变换

通过Radon反变换来恢复图像

close all; clear all; clc;
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));    %显示图像

投影角度的多少对Radon变换和反变换的影响

close all; clear all; clc;
I=phantom(256);   %读入图像
figure;
imshow(I);         %显示图像
theta1=0:10:170; 
theta2=0:5:175; 
theta3=0:2:178;             
[R1,xp]=radon(I,theta1);   %Radon变换
[R2,xp]=radon(I,theta2);
[R3,xp]=radon(I,theta3);
figure; 
imagesc(theta3,xp,R3);     %显示结果
colormap hot;              %设置调色板
colorbar;                  %添加颜色条
J1=iradon(R1,10);         %Radon反变换
J2=iradon(R2,5); 
J3=iradon(R3,2); 
figure;
subplot(131);imshow(J1);    %显示结果图像   
subplot(132);imshow(J2);
subplot(133);imshow(J3);    

图像傅里叶变换

傅里叶变换的物理意义

傅里叶变换的物理意义是：将图像的灰度分布函数变换为图像的频率分布函数；傅里叶逆变换：将图像的频率扥不函数变换为灰度分布函数

对图像进行二维傅里叶变换得到频谱图，就是图像梯度的分布图，傅里叶频谱图上看到的明暗不易的亮点，是图像上某一点与邻域点差异的强弱，即梯度的大小，即该点的频率大小。

傅里叶变换的定义及性质

1. 连续傅里叶变换
2. 离散傅里叶变换
3. 二维离散傅里叶变换的性质
   1）可分离性
   2）线性
   3）共轭对称性
   4）位移性
   5）尺度变换性
   6）旋转不变性
   7）卷积性
   8）DC系数

傅里叶变换的MATLAB实现

fft()快速FFT算法：
矩阵的二维离散傅里叶变换

close all; clear all; clc;
I1=ones(4)          
I2=[2 2 2 2;1 1 1 1;3 3 0 0;0 0 0 0]
J1=fft(I1)
J2=fft(I2)

图像的二维离散傅里叶变换

close all; clear all; clc;
I=imread('cameraman.tif'); 
J=fft2(I);
K=abs(J/256);
figure;
subplot(121);imshow(I);
subplot(122);imshow(uint8(K));

注：通过函数fft2()得到的频谱，坐标原点位于左上角。可以通过函数fftshift()将变换后的坐标原点移到频谱图窗口中央，坐标原点是低频，向外是高频。

通过函数fftshift()进行平移：

close all; clear all; clc;
N=0:4
X=fftshift(N)            %平移
Y=fftshift(fftshift(N))  %平移后再进行平移
Z=ifftshift(fftshift(N)) %平移后再进行反平移

图像通过傅里叶变换和平移

close all; clear all; clc;
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));

图像变亮后进行傅里叶变换

close all; clear all; clc;
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));

图像旋转后进行傅里叶变换

close all; clear all; clc;
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));

图像中添加高斯噪声后进行傅里叶变换

close all; clear all; clc;
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));

灰度图像的傅里叶变换和反变换

close all; clear all; clc;
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));

灰度图像的幅值谱和相位谱

close all; clear all; clc;
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;  %归一化
B=angle(K);                   %相位谱
figure;
subplot(121);imshow(A);
subplot(122);imshow(real(B));

编程实现二维离散傅里叶变换

close all; clear all; clc;
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);   %采用函数fft2()
figure;
subplot(121);imshow(K);
subplot(122);imshow(F);

傅里叶变换的应用

通过傅里叶变换识别图像中的字符

close all; clear all; clc;
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)

对图像进行巴特沃斯低通滤波

close all; clear all; clc;
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));

对图像进行巴特沃斯高通滤波

close all; clear all; clc;
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));

图像离散余弦变换

离散余弦变换（DCT）是以一组不同频率和幅值的余弦函数和来近似一幅图像，实际上是傅里叶变换的实数部分。离散余弦变换经常用于图像压缩。

离散余弦变换的定义

离散余弦变换的MATLAB实现

函数dct()：一维离散余弦变换
函数idct()：一维离散余弦反变换

对图像进行二维离散余弦变换

close all; clear all; clc;
I=imread('coins.png'); 
I=im2double(I);
J=dct2(I);   %二维离散余弦变换
figure;
subplot(121);imshow(I);
subplot(122);imshow(log(abs(J)),[]);

通过函数dctmtx()生成离散余弦变换矩阵

close all; clear all; clc;
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()进行图像的离散余弦变换

close all; clear all; clc;
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);

图像的二维离散反余弦变换

close all; clear all; clc;
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);

离散余弦变换的应用

通过函数blkproc()对图像进行块操作

close all; clear all; clc;
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;             %设置调色板

通过离散余弦操作进行图像压缩

close all; clear all; clc;
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              %只选择左上角的10个系数
    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);

其他图像变换

Hadamard变换

Hadamard变换：相当于在原来的图像矩阵左右分别乘以一个矩阵，这两个矩阵都是正交矩阵。Hadamard变换矩阵中所有的元素都是+1或-1。

通过函数hadamard()产生Hadamard变换矩阵

close all; clear all; clc;
H1=hadamard(2)   %产生2阶Hadamard变换矩阵
H2=hadamard(4)   %产生4阶Hadamard变换矩阵
H3=H2'*H2        %验证

对图像进行Hadamard变换

close all; clear all; clc;
I=imread('peppers.png');
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);

Hough变换

Hough变换是图像处理中从图像中识别几何形状的基本方法之一。
基本原理：利用点与线的对偶性，将原始图像空间给定的曲线通过曲线表达形式变为参数空间的一个点。把原始图像中给定曲线的检测问题转化为寻找参数空间中的峰值问题。

函数[H, Theta, Rho]=hough(BW)：对二值图像BW进行Hough变换，返回值H为Hough变换矩阵，theta为变换角度，rho为变换半径r。
对图像进行Hough变换

close all; clear all; clc;
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)));   %显示Hough变换结果
axis normal;    %设置坐标轴
hold on;
colormap hot;   %设置调色板

函数houghpeaks()：用于在Hough变换后的矩阵中寻找最值，该最值可以用于定位直线段。
函数houghlines()：用于绘制找到的直线段

通过图像的Hough变换检测直线

close all; clear all; clc;
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'); %绘制最长的直线

注：首先通过函数hough()对图像进行Hough变换，然后采用函数houghpeaks()获取矩阵中较大的5个点，通过函数houghlines()获取线段的端点。通过循环语句，绘制多条直线，同时将最长的直线设置为蓝色。