数字图像处理是一门集计算机科学、光学、数学、物理学等多学科的综合科学。随着计算机科学的发展，数字图像处理技术取得了巨大的进展，呈现出强大的生命力，已经在多种领域取得了大量的应用，推动了社会的发展。其中，遥感领域中，对于影像数据的处理均基于数字图像处理的技术。而遥感影像数据作为地理信息科学的重要数据源，如何从中获取有用的信息，是地理信息数据处理中重要的内容。
 
MATLAB作为数学领域应用最广泛的一种软件，集成了对于图片处理的函数和功能，成为了处理数字图像问题的佼佼者。其出众的计算能力和简便的绘图能力可以有效进行数字图像的变换和操作。本文探究了MATLAB(R2014a)软件下数字图像处理部分简单内容方法的操作实现，包括图像变换、图像增强等。
 
由于时间所迫和水平不足，对于本文中的问题请指出。
 
在这里首先将介绍基本的操作。
 
1　　数字图像处理基础
 
1.1            MATLAB图像处理基本操作
 
本文中对于大多数的操作，是对数字图像处理领域中最为著名的“lena”图片进行操作的。原图如下(Figure 1)：
 

 
Figure 1
 
首先，在MATLAB中显示这幅图片：
 
1
 
2
 
3I=imread('lena.jpg');
 
imfinfo('lena.jpg')
 
imshow(I);
 
其中，imread()可以读取图片，imfinfo()可以获取图片的信息(Figure 2)，imshow()可以显示图片。
 

 
Figure 2
 
从Figure 2中可以看出，该图片格式为png，长宽均为512px，颜色类型为真彩色。
 
处理后的图片还需要进行保存：
 
1imwrite(I,'lenaSave.jpg');
 
1.2            图像数字化
 
使用imread()读取图像后，可以看到，读入的图片I是以一个512*512*3的矩阵进行保存的，即分别是RGB颜色的数字图片。为了方便起见，在本文中，笔者将lena图进行了转换，转换为灰度图，即保存为一个512*512*1的矩阵。
 
为了完成这样一个步骤，MATLAB中提供了相应的函数：
 
1
 
2I1=rgb2gray(I);
 
imshow(I1);
 
I1以512*512*1的矩阵形式进行存储。效果如下图(Figure 3)：
 

 
Figure 3
 
1.3            直方图
 
为了显示图像灰度的分布情况，还需要绘制灰度直方图。可以使用如下代码：
 
1figure,imhist(I1);
 
效果如下(Figure 4)：
 

 
Figure 4
 
1.4            图像噪声的添加
 
为了完成多种图像处理的操作和试验，还可以对图片添加噪声。所用函数为imnoise (I, type)，该函数中的type可以为5种噪声参数，分别为：'gaussian'(高斯白噪声)，'localvar'(与图象灰度值有关的零均值高斯白噪声)，'poisson'(泊松噪声)，'salt & pepper'(椒盐噪声)和'speckle'(斑点噪声)。
 
以下为高斯噪声(Figure 5)和椒盐噪声(Figure 6)的代码和效果：
 
1
 
2I2=imnoise(I1,'gaussian');
 
figure,imshow(I2);
 

 
Figure 5
 
1
 
2I3=imnoise(I1,'salt & pepper');
 
figure,imshow(I3);
 

 
Figure 6
 
2         图像变换
 
2.1            傅立叶变换
 
傅立叶变换可以将图像从空间域转换到频率域，然后再进行相应的处理。MATLAB中有傅立叶变换的函数。
 
原图：
 

 
Figure 7
 
傅立叶变换：
 
1
 
2
 
3F=fft2(I1);
 
S=abs(F);
 
figure,imshow(S,[]);
 
效果：
 

 
Figure 8
 
看上去是一片黑，事实上在全图的左上角可以看出有一个白点(即左上角缺失了)。
 
平移：
 
1
 
2Fc=fftshift(F);
 
figure,imshow(abs(Fc),[]);
 
效果：
 

 
Figure 9
 
将频率平移到中间。
 
频谱图：
 
1figure,imshow(uint8(abs(Fc/256)));
 
效果：
 

 
Figure 10
 
之所以在本例中没有使用lena的图片，是因为该图傅立叶频谱图并无特点。
 
傅立叶逆变换：
 
1
 
2
 
3f=real(ifft2(F)/255);
 
figure,imshow(f);
 
效果：
 

 
Figure 11
 
经过逆变换后，仍然为原图。

 
 
 
   
   
  
注：本系列来自于图像处理课程实验，用Matlab实现最基本的图像处理算法
 
 


1.Fourier变换
 
（1）频域增强
 
除了在空间域内可以加工处理图像以外，我们还可以将图像变换到其他空间后进行处理，这些方法称为变换域方法，最常见的变换域是频域。
 
使用Fourier变换把图像从空间域变换到频域，在频域内做相应增强处理，再从频域变换到空间域得到处理后的图像。
 
我们这里主要学习Fourier变换和FFT变换的算法，没有学过通信原理，我对信号、时域分析也不是很清楚。
 

2.FFT算法
 
（1）离散Fourier变换，DFT
 
函数的计算式为：
 

的FFT算法。 
（2）快速Fourier变换，FFT
 
快速傅立叶变换FFT是利用单位复数根的特殊性质（消去引理和折半引理，见《算法导论》（第3版中文版）P532详细论述），在时间内计算出DFT的一种快速算法。
 
FFT有两种基本实现方式：
 
递归FFT
迭代FFT，也叫FFT蝶式运算
递归FFT由于递归栈开销大且容量有限等缺陷（但理解容易），一般计算采取迭代FFT实现。
 
（3）迭代FFT
 
直接给出算法导论版本的迭代FFT算法：
 
 
其中求逆序数拷贝的函数为：
 
 
FFT采用折半迭代的思想，因此速度能降低到。
 
（4）迭代FFTMatlab实现
 
Matlab有fft函数，我们也可以自己实现：
 
function [ fft_m ] = IterativeFFT( vec )
    clear i;

    n = length(vec);

    fft_m = BitReverseCopy(vec);
    for s = 1 : log2(n)
        m = power(2, s);
        wm = exp(- 2 * pi * i / m);

        for k = 0 : m : n - 1
            w = 1;
            for j = 0 : m / 2 - 1
                t = w * fft_m(k + j + m / 2 + 1);
                u = fft_m(k + j + 1);
                fft_m(k + j + 1) = u + t;
                fft_m(k + j + m / 2 + 1) = u - t;
                w = w * wm;
            end
        end
    end
end
 

 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 
BitReverseCopy函数如下：
 
function [ copy ] = BitReverseCopy( vec )
    n = length(vec);
    copy = zeros(1, n);

    bitn = log2(n);

    for i = 0 : n - 1
        revi = bin2dec(fliplr(dec2bin(i, bitn)));
        copy(revi + 1) = vec(i + 1);
    end
end
 

 
1
2
3
4
5
6
7
8
9
10
11
 
需要特别注意的是:
 
一般给出的FFT算法伪代码都是基于下标从零开始的数组，写在Matlab需要考虑映射关系
clear i是为了怕之前有变量i和复数记号i混淆，清楚Matlab workspace中的缓存
默认vec是double类型的！因为中间采用许多double类型运算

3.图像的二维Fourier变换
 
（1）二维DFT
 
二维DFT定义公式为： 

的时间，效率很低。 
（2）将二维DFT分解为两个一维DFT
 
Fourier变换的变换核（公式中和
再对变换后的矩阵的每一行做一维DFT：
 

的时间完成二维傅立叶变换。 
（3）用一维FFT实现二维FFT
 
同样的我们可以用两个一维FFT实现二维FFT，时间复杂度：
 
function [ mfft2 ] = JCGuoFFT2( data )
    h = size(data, 1);
    w = size(data, 2);
    mfft2 = data;

    if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w
        disp('JCGuoFFT2 exit: h and w must be the power of 2!')
    else
        for i = 1 : h
            mfft2(i, :) = IterativeFFT(mfft2(i, :));
        end

        for j = 1 : w
            mfft2(:, j) = IterativeFFT(mfft2(:, j));
        end
    end
end

 

 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
代码很简单，先做FFT行变换再做FFT列变换。之前忘记提到，我这里实现的都是长度必须是2的次方的Fourier变换，因此有时候会做一些长度规范检查。
 
（4）变换结果
 
经过JCGuoFFT2的二维傅立叶变换，我们可以得到复平面内各个点的复数值，那么显示在图像中需要先求出幅值：
 
pic1_fft = JCGuoFFT2(double(pic1));
pic1_fft_amp = abs(pic1_fft);
 

 
1
2
 
在对幅值做一次log变换得到较好的频域图像：
 
pic1_fft_amp_log = log(1 + pic1_fft_amp);
 

 
1
 
绘制结果如下图：
 
 
（5）低频信号移到图像中心点
 
由于复数运算的周期特性，图像的Fourier变换在复平面上是完全对称的，可以想象平面是由无限多个上图（右）频域图像拼接而成的二维阵列。一般研究频域图像是把低频部分，也就是变换后的边缘部分移到图像中心点，Matlab提供fftshift函数完成平移。
 
平移的思路有两个
 
通过Fourier变换平移定理先把原始图像做变换再做FFT
先做FFT后再根据频域图像的对称性做对称变换
查阅Matlab文档发现它是采用第二种方法，对图像做以下子矩阵交换：
 
 
那么我们可以很容易的写出自己的fftshift：
 
function [ after ] = FFTShift( before )
    h = size(before, 1);
    w = size(before, 2);

    after = before;

    if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w
        disp('FFTShift exit: h and w must be the power of 2!')
    else
        hh = h / 2;
        hw = w / 2;
        after(1 : hh, 1 : hw) = before(hh + 1 : h, hw + 1 : w);
        after(1 : hh, hw + 1 : w) = before(hh + 1 : h, 1 : hw);
        after(hh + 1 : h, 1 : hw) = before(1 : hh, hw + 1 : w);
        after(hh + 1 : h, hw + 1 : w) = before(1 : hh, 1 : hw);
    end
end

 

 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
将低频部分平移到中心点后结果为：
 
 

4.图像的二维反Fourier变换
 
（1）一维逆DFT和一维逆FFT
 
一维离散傅立叶变换的逆变换是将e的指数部分变号，然后整体除以长度N（Fourier变换与逆变换关于符号、系数有很多种组合的定义，但他们都是等价且对称的。这里的定义配合Matlab做fft实际效果。）：
 


同样我们可以根据iDFT的定义推演iFFT的算法，其迭代版本的Matlab实现如下：
 
function [ ifft_m ] = IterativeIFFT( vec )
    clear i;
    n = length(vec);
    ifft_m = BitReverseCopy(vec);
    for s = 1 : log2(n)
        m = power(2, s);
        wm = exp(2 * pi * i / m);

        for k = 0 : m : n - 1
            w = 1;
            for j = 0 : m / 2 - 1
                t = w * ifft_m(k + j + m / 2 + 1);
                u = ifft_m(k + j + 1);
                ifft_m(k + j + 1) = u + t;
                ifft_m(k + j + m / 2 + 1) = u - t;
                w = w * wm;
            end
        end
    end
    ifft_m = ifft_m ./ n;
end

 

 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 
（2）二维逆FFT
 
二维逆FFT和二维FFT的思路一致，对所有行和列分别作一次iFFT即可：
 
function [ mifft2 ] = JCGuoIFFT2( data )
    h = size(data, 1);
    w = size(data, 2);
    mifft2 = data;

    if power(2, log2(h)) ~= h || power(2, log2(w)) ~= w
        disp('JCGuoIFFT2 exit: h and w must be the power of 2!')
    else
        for i = 1 : h
            mifft2(i, :) = IterativeIFFT(mifft2(i, :));
        end

        for j = 1 : w
            mifft2(:, j) = IterativeIFFT(mifft2(:, j));
        end
    end
end

 

 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
（3）逆FFT结果
 
对之前Rect1.bmp用JCGuoFFT2变换的到的Fourier变换（非shift和log之后、非仅幅度部分）作FFT2逆变换可以直接得到原图像：
 
 
这幅图有多个结果，可以看title知道每个结果的含义，图2-1是用JCGuoIFFT2做傅立叶反变换的结果，得到的图像和原图像、和Matlab ifft2反变换后的图像都是一致的。
 

5.幅频特性与相频特性
 
（1）对振幅和相位单独做逆FFT
 
如果我们只把图像Fourier变换的振幅部分和相位部分单独做二维逆FFT，可以直观的感受人眼对图像幅频特性和相频特性的敏感度。
 
复数的幅度/振幅定义为：
 


相位角和相位定义为：
 


对相位反变换需要在乘以一个系数（我是调出来的，针对图像，肯定可以自动的做均衡化）：
 
pic2_fft_angle = angle(pic2_fft);
clear i;
tmp = 10000 * exp(i * pic2_fft_angle);
pic2_ifft_angle = uint8(JCGuoIFFT2(tmp));

 

 
1
2
3
4
5
 
对振幅和相位单独做逆FFT结果如下：
 
 
（2）人眼敏感度
 
观察上图，幅频特性主要涵盖了图像颜色的分布，相频特性主要刻画了图像的边界信息。人眼对图像的相频特性更加敏感，看相频特性能够大概知道图像内容。
 

6.Fourier变换的旋转定理
 
（1）Fourier变换旋转定理
 


（2）结果
 
Rect2.bmp是Rect1.bmp旋转45度的示意图（注：原Rect2.bmp是二值的，做了预处理，但其本身边界不平滑，导致效果不太好，但不影响观察旋转定理）：
 
 
我们可以看到幅度FFT正变换和相位FFT你变换都是旋转了45度，但是幅度FFT逆变换区别较大。
 
 

 

 
实验四 图像的傅立叶变换与频域滤波
 
一、 实验目的
 
1了解图像变换的意义和手段；
 
2熟悉傅里叶变换的基本性质；
 
3熟练掌握FFT 方法的应用；
 
4通过实验了解二维频谱的分布特点；
 
5通过本实验掌握利用MATLAB 编程实现数字图像的傅立叶变换。
 
6、掌握怎样利用傅立叶变换进行频域滤波
 
7、掌握频域滤波的概念及方法
 
8、熟练掌握频域空间的各类滤波器
 
9、利用MATLAB 程序进行频域滤波
 
二、 实验原理
 
1应用傅立叶变换进行图像处理
 
傅里叶变换是线性系统分析的一个有力工具，它能够定量地分析诸如数字化系统、采样点、电子放大器、卷积滤波器、噪音和显示点等的作用。通过实验培养这项技能，将有助于解决大多数图像处理问题。对任何想在工作中有效应用数字图像处理技术的人来说，把时间用在学习和掌握博里叶变换上是很有必要的。
 
2傅立叶(Fourier )变换的定义
 
对于二维信号，二维Fourier 变换定义为 :
 
??∞
 
∞-+-==dxdy e y x f v u F y x f F vy ux j )(2),(),()},({π 二维离散傅立叶变换为：
 
∑
 
∑-=+--==10)(21
 
01
 
),(),(N y N y u M x u j M x MN e y x f v u F π
 
图像的傅立叶变换与一维信号的傅立叶变换变换一样，有快速算法，具体参见参考书目，有关傅立叶变换的快速算法的程序不难找到。实际上，现在有

一、获取代码方式
 
获取代码方式1：
 完整代码已上传我的资源：【图像几何】基于matlab GUI傅立叶+DCT图像变换【含Matlab源码 844期】
 
获取代码方式2：
 通过紫极神光博客主页开通CSDN会员，凭支付凭证，私信博主，可获得此代码。
 
获取代码方式3：
 通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。
 
备注：开通CSDN会员，仅只能免费获得1份代码（有效期为开通日起，三天内有效）；
 订阅紫极神光博客付费专栏，可免费获得2份代码（有效期为订阅日起，三天内有效）；
 
二、傅立叶+DCT简介
 
FFT(Fast Fourier transform)：快速傅里叶变换，是DFT的工程化实现方法。
 DFT直接求解太过于复杂，FFT方法根据DFT求解过程中旋转因子的性质并引入分治算法思想，大大简化计算过程，被广泛应用在频谱分析的工程实践中，如matlab，C，C++，CUDA等底层实现
 
1 DFT简介
 频谱分析是信号处理中的重要环节，从傅里叶变换FT，到拉普拉斯变换LT，离散时间傅里叶变换DTFT，Z变换ZT，到我们所讲的离散傅里叶变换DFT（他们之间的联系和区别见我的其他博客）。
 相比于其他变换，DFT被广泛应用的原因是其输入的时域信号是离散的，输出的频域结果也是离散的。这就极大方便了我们进行基于计算机的频谱计算，存储和分析，没办法数字信号处理是大趋势。
 ​ DFT变换的公式为：
 
 
 但为了分析方便，在FFT的计算过程中，我们依然使用k = 0 ∼ N − 1 k = 0 \sim N-1k=0∼N−1的选取策略。也即，如下：
 
 2 旋转因子WWW的性质
 
 3 FFT蝶形计算证明
 
 
 
 
 
 
 4 FFT计算过程
 
 
三、部分源代码
 
function varargout = fushidctbianhuan(varargin)
% FUSHIDCTBIANHUAN M-file for fushidctbianhuan.fig
%      FUSHIDCTBIANHUAN, by itself, creates a new FUSHIDCTBIANHUAN or raises the existing
%      singleton*.
%
%      H = FUSHIDCTBIANHUAN returns the handle to a new FUSHIDCTBIANHUAN or the handle to
%      the existing singleton*.
%
%      FUSHIDCTBIANHUAN('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in FUSHIDCTBIANHUAN.M with the given input arguments.
%
%      FUSHIDCTBIANHUAN('Property','Value',...) creates a new FUSHIDCTBIANHUAN or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before fushidctbianhuan_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to fushidctbianhuan_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help fushidctbianhuan

% Last Modified by GUIDE v2.5 06-Apr-2014 17:16:14

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @fushidctbianhuan_OpeningFcn, ...
                   'gui_OutputFcn',  @fushidctbianhuan_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


% --- Executes just before fushidctbianhuan is made visible.
function fushidctbianhuan_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to fushidctbianhuan (see VARARGIN)
axes(handles.axes1);
img=imread('lena.bmp');
imshow(img);
axes(handles.axes2);
imgf=double(img);
imgs=fftshift(fft2(imgf));
imshow(log(1+abs(imgs)),[]);
set(handles.dct_radiobutton,'value',0);
set(handles.fft_radiobutton,'value',1);
% Choose default command line output for fushidctbianhuan
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes fushidctbianhuan wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = fushidctbianhuan_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in fft_radiobutton.
function fft_radiobutton_Callback(hObject, eventdata, handles)
% hObject    handle to fft_radiobutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.dct_radiobutton,'value',0);
set(handles.fft_radiobutton,'value',1);
img = getimage(handles.axes1);
axes(handles.axes2);
imgf=double(img);
imgs=fftshift(fft2(imgf));
imshow(log(1+abs(imgs)),[]);
 set(handles.text2,'string','傅立叶变换谱');
 set(handles.text2,'fontsize',12);
 set(handles.text2,'fontname','宋体');
 set(handles.text2,'fontweight','bold');

% Hint: get(hObject,'Value') returns toggle state of fft_radiobutton


% --- Executes on button press in dct_radiobutton.
function dct_radiobutton_Callback(hObject, eventdata, handles)
% hObject    handle to dct_radiobutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
set(handles.dct_radiobutton,'value',1);
set(handles.fft_radiobutton,'value',0);
img = getimage(handles.axes1);
axes(handles.axes2);
imgd=img;
imgs=dct2(imgd);
imshow(log(abs(imgs)));
set(handles.text2,'string','DCT变换系数');
set(handles.text2,'fontsize',12);
set(handles.text2,'fontname','宋体');
set(handles.text2,'fontweight','bold');
% Hint: get(hObject,'Value') returns toggle state of dct_radiobutton


% --- Executes on button press in close_button.
function close_button_Callback(hObject, eventdata, handles)
% hObject    handle to close_button (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
close(fushidctbianhuan);

% --- Executes on selection change in fft_dct_pop.
function fft_dct_pop_Callback(hObject, eventdata, handles)
% hObject    handle to fft_dct_pop (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
val=get(hObject,'value');
str=get(hObject,'string');
switch str{val};
    case 'lena'
        lena = [];
        lena = imread('lena.bmp');
        img = lena;
    case 'Saturn',
        saturn = [];
        load imdemos saturn
        img = saturn;
    case 'Tire',
        tire = [];
        load imdemos tire
        img = tire;
    case 'Pout',
        pout = [];
        load imdemos pout
        img = pout;
    case 'Trees',
        trees = [];
        load imdemos trees
        img = trees;
    case 'Quarter',
        quarter = [];
        load imdemos quarter
        img = quarter;
    case 'Circuit',
        circuit = [];
        load imdemos circuit
        img = circuit;
    case 'Rice',
        rice = [];
        load imdemos rice
        img = rice;
    case 'Fingerprint'
        fingerprint = [];
        fingerprint = imread('fingerprint.jpg');
        img = fingerprint;
    case 'Licenceplate'
        licenceplate = [];
        licenceplate = imread('licenceplate.jpg');
        img = licenceplate;
    case 'Haze'
        haze = [];
        haze = imread('haze.jpg');
        img = haze;
    case 'Cloudy'
        cloudy = [];
        cloudy = imread('cloudy.tif');
        img = cloudy;
end
axes(handles.axes1);
imshow(img);
control=get(handles.fft_radiobutton,'value');
switch control;
 
四、运行结果
 
 
五、matlab版本及参考文献
 
1 matlab版本
 2014a
 
2 参考文献
 [1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
 [2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
 [3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
 [4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.
 [5]陈浩,方勇,朱大洲,王成,陈子龙.基于蚁群算法的玉米植株热红外图像边缘检测[J].农机化研究. 2015,37(06)