图像处理入门4-傅里叶变换

图像处理4 傅里叶变换

  图像的空间域和频率域给我们提供了不同的视角，但所表达的内容是同一个，有点像“横看成岭侧成峰，远近高低各不同”。在空域中，自变量$(x,y)$表示图像中的位置，而函数$f(x,y)$表示图像在$(x,y)$处的像素值。换个角度，如果将图像看成二维信号，那么傅里叶变换是一种将时域信息转化为频域信息的一种数学变换，这样我们可以在频域上对信息进行提取、分析、操作。然后再反变换回空域中，完成想要的结果。关于傅里叶变换详细讲解可以参考：傅里叶变换

傅里叶变换原理

要理解傅里叶变换，必要的数学是必不可少的，做图像处理最需要的是理论基础和观察，傅里叶变换对于大学以上学过高等数学的人来说应该不算困难，困难在如何理解含义并将变换理论应用到实际工作中，但是实际上不同人群的数学基础差异较大，所以先简要说一下它的基础知识，有利于于掌握这个工具。

傅里叶级数

法国数学家在解热传导方程时，发现分离变量这种解法，这种方法给他带来了灵感。从此傅里叶级数在很多场合大放异彩。回到级数自身，它要求函数满足周期性和第一类迪利克雷间断点，才可以展开成这种级数，经过数学家添加延拓条件，对于非周期也可以使用。而且对于图像处理来说可以认为是“连续”的，所以可以直接使用。

周期为$T$的函数$f(x)$展开成三角函数形式：
$$f(x)=\frac{a_0}{2}+\sum _{k=1}^{+\infty }(a_{k}cos(n{\omega }_{0}x)+b_{k}sin(n{\omega }_{0}x))$$
其中：${\omega }_0=\frac{2\pi }{T}=2\pi u,u=\frac{1}{T}$是函数频率。$a_0,a_k,b_k$称为傅里叶系数。
傅里叶级数的思想是将一个周期函数用三角函数进行分解，然后对不同频率的函数做相关的处理。

傅里叶级数用一组正交三角函数作为无穷维基，正交基有许多良好性质。

傅里叶变换及其离散化

从傅里叶级数过渡到傅里叶变换，在定义域为时间轴的$(-\infty <t<\infty )$的非周期函数$f(t)$,无法通过周期延拓成周期函数，此种情况就用到了傅里叶变换。
变换：$F(u)={\int }_{-\infty }^{\infty }f(x)e^{-i2\pi ux}dx$
反变换：$f(x)={\int }_{-\infty }^{\infty }F(u)e^{-i2\pi ux}du$

图像是有限离散形式的，所以傅里叶变换需要换成离散的。二维离散傅里叶变换DFT可分离性的基本思想是DFT可分离为两次一维DFT。因此可以用通过计算两次一维的FFT来得到二维快速傅里叶FFT算法。根据快速傅里叶变换的计算要求，需要图像的行数、列数均满足2的n次方，如果不满足，在计算FFT之前先要对图像补零以满足2的n次。一维描述也方便易理解：

离散傅里叶变换：$$F(n)=\sum _{n=0}^{N-1}f(t)e^{-i\frac{2\pi }{N}nt}(t=0,1,2,...,N-1)$$
离散傅里叶逆变换：$$f(t)=\frac{1}{N}\sum _{n=0}^{N-1}F(n)e^{i\frac{2\pi }{N}nt}(n=0,1,2,...,N-1)$$
令$W_t^n=e^{-i\frac{2\pi }{N}nt}$，离散傅里叶变换写成矩阵形式
$$\left[\begin{array}{c}F(0)\\ F(1)\\ F2)\\ \cdot \\ \cdot \\ F(N-2)\\ F(N-1)\end{array}\right]=\left[\begin{array}{cccccc}{W}_0^0& W_1^0& \cdot & \cdot & \cdot & W_{N-1}^0\\ W_0^1& W_1^1& \cdot & \cdot & \cdot & W_{N-1}^1\\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ W_0^{N-1}& W_1^{N-1}& \cdot & \cdot & \cdot & W_{N-1}^{N-1}\end{array}\right]\left[\begin{array}{c}f(0)\\ f(1)\\ f(2)\\ \cdot \\ \cdot \\ f(N-2)\\ f(N-1)\end{array}\right]$$
反变换：

$$\left[\begin{array}{c}f(0)\\ f(1)\\ f2)\\ \cdot \\ \cdot \\ f(N-2)\\ f(N-1)\end{array}\right]=\left[\begin{array}{cccccc}{W}_0^0& W_1^0& \cdot & \cdot & \cdot & W_{N-1}^0\\ W_0^{-1}& W_1^{-1}& \cdot & \cdot & \cdot & W_{N-1}^{-1}\\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ W_0^{-(N-1)}& W_1^{-(N-1)}& \cdot & \cdot & \cdot & W_{N-1}^{-(N-1)}\end{array}\right]\left[\begin{array}{c}F(0)\\ F(1)\\ F(2)\\ \cdot \\ \cdot \\ F(N-2)\\ F(N-1)\end{array}\right]$$

快速傅里叶变换

直接计算离散傅里叶变换的话，时间复杂度$O(N^2)$，随着$N$增加，计算显著变慢。研究人员在离散傅里叶变换基础上又发展出了快速傅里叶变换，所谓快速就是计算复杂度降低，计算量变小。
考虑到离散傅里叶变换时系数矩阵存在对称性、奇偶性和周期性，简化矩阵计算，采用分治思想进行迭代。（想详细了解可以参考：1.快速傅里叶变换——线性代数视角算法剖析，2快速傅里叶变换-蝶形算法）

实例：

1.傅里叶滤波

*读取图像
read_image (Image, 'pic_scratch.bmp')
*反色
invert_image (Image, ImageInverted)
*获取图像尺寸
get_image_size (Image, Width, Height)

*设置频域滤波器
gen_sin_bandpass (ImageBandpass, 0.4, 'none', 'rft', Width, Height)
*快速傅里叶变换
rft_generic (ImageInverted, ImageFFT, 'to_freq', 'none', 'complex', Width)
*用滤波器对变换的频域滤波
convol_fft (ImageFFT, ImageBandpass, ImageConvol)
*反变换到时域
rft_generic (ImageConvol, Lines, 'from_freq', 'n', 'byte', Width)

*二值化提取区域
threshold (Lines, Region, 5, 255)
*分离不连通区域
connection (Region, ConnectedRegions)

*选择面积大于5小于5000的区域
select_shape (ConnectedRegions, SelectedRegions, 'area', 'and', 5, 5000)

*腐蚀扩展
dilation_circle (SelectedRegions, RegionDilation, 5.5)
*将区域合并
union1 (RegionDilation, RegionUnion)
reduce_domain (Image, RegionUnion, ImageReduced)
*从图像提取线段
lines_gauss (ImageReduced, LinesXLD, 0.8, 3, 5, 'dark', 'false', 'bar-shaped', 'false')
*合并线段
union_collinear_contours_xld (LinesXLD, UnionContours, 40, 3, 3, 0.2, 'attr_keep')
*选择线段长度在15到1000之间的有缺陷的线段
select_shape_xld (UnionContours, SelectedXLD, 'contlength', 'and', 15, 1000)
*生成区域
gen_region_contour_xld (SelectedXLD, RegionXLD, 'filled')
union1 (RegionXLD, RegionUnion)
dilation_circle (RegionUnion, RegionScratches, 10.5)

原图

滤波器　 　 　 　 　 　 　 　 　 　 　 　 　 　 　 变换后

划伤位置

2.傅里叶定位

除频域滤波之外，傅里叶可以用相位分析计算两个图像之间的变换，具体思路为：
1.读取图像，并设定平移参数，利用仿射变换将图像进行平移。这步是为了傅里叶相空间分析做准备；
2.利用二阶多项式创建一个灰色曲面，并将其与平移后的图像进行叠加
3.对原始图像以及叠加后图像进行傅里叶变换，都从从空间域转换到频率域
4.计算上述两幅频率域图像的相位关系
5.将相位相关性图像反变换到空间域，利用图像中心纠正和local_max_sub_pix算子计算出最后的偏移值
6.在窗口上进行相关显示

dev_update_off ()
read_image (Image, 'wafer/wafer_dies.png')
get_image_size (Image, Width, Height)
dev_close_window ()
dev_open_window (0, 0, Width, Height, 'black', WindowHandle)
set_display_font (WindowHandle, 16, 'mono', 'true', 'false')
dev_display (Image)
disp_message (WindowHandle, 'Original image', 'window', 12, 12, 'black', 'true')
disp_continue_message (WindowHandle, 'black', 'true')
stop ()
* Translate the image.
RowTrans := 92.4
ColumnTrans := 58.2
optimize_rft_speed (Width, Height, 'standard')  //优化傅里叶变换的运行时间
hom_mat2d_identity (HomMat2D)  //生成二维其次变换矩阵
hom_mat2d_translate (HomMat2D, RowTrans, ColumnTrans, HomMat2D) //将平移添加到齐次二维变换矩阵。
* We set 'init_new_image' to 'true' to ensure the translated image has
* defined values of 0 in the part that lies outside the original image.
*我们将“init_new_image”设置为“true”，以确保转换后的图像

*在原始图像之外的像素值为0。
get_system ('init_new_image', InitNewImage) //获取Halcon系统的当前值
set_system ('init_new_image', 'true')  //设置参数

affine_trans_image (Image, ImageTrans, HomMat2D, 'constant', 'false')  //仿射变换图像（将图像原点平移到(RowTrans ,ColumnTrans)）

full_domain (ImageTrans, ImageTranslated) //将图像的域扩展到最大。

set_system ('init_new_image', InitNewImage)
* Simulate a degradation of the image by an uneven illumination. 通过不均匀照明模拟图像的退化，(通过与二阶灰色曲面叠加进行实现)
*用二阶多项式创建一个灰色曲面。
gen_image_surface_second_order (ImageSurface, 'byte', 0.0005, -0.0008, 0, 0, 0, 128, Height / 2, Width / 2, Width, Height)
add_image (ImageTranslated, ImageSurface, ImageDegraded, 1, -128) //叠加两幅图像  g' := (g1 + g2) * Mult + Add
dev_display (ImageDegraded)
disp_message (WindowHandle, ['Image translated by (' + RowTrans$'4.1f' + ',' + ColumnTrans$'4.1f' + ')','and degraded by uneven illumination'], 'window', 12, 12, 'black', 'true')
disp_continue_message (WindowHandle, 'black', 'true')
stop ()
* Compute the phase correlation. 计算相位相关
rft_generic (Image, ImageFFT, 'to_freq', 'none', 'complex', Width) //原始图像  空间域到频率域  傅里叶变换
rft_generic (ImageDegraded, ImageTransFFT, 'to_freq', 'none', 'complex', Width) //模拟的退化后图像 空间域到频率域  傅里叶变换

phase_correlation_fft (ImageFFT, ImageTransFFT, ImagePhaseCorrelationFFT) //计算两幅图像在频域中的相位相关性

rft_generic (ImagePhaseCorrelationFFT, ImagePhaseCorrelation, 'from_freq', 'n', 'real', Width) //相位相关性图像  频率域到空间域 傅里叶变换
dev_display (ImagePhaseCorrelation)
gen_circle_contour_xld (Circle, RowTrans, ColumnTrans, 20, 0, 6.28318, 'positive', 1) //在平移后的原点上生成一个圆形轮廓
dev_set_color ('green')
dev_display (Circle)
*注意左上角的峰值
disp_message (WindowHandle, ['Phase correlation image:','note the peak in the upper left corner'], 'window', 12, 12, 'black', 'true')
disp_continue_message (WindowHandle, 'black', 'true')
stop ()
* Since the phase correlation is cyclic, negative translations result in
* peaks in lower or right part of the image.  If the translation is close
* 0 in one or two directions, the interpolation in local_max_sub_pix would
* therefore access wrong values because of its border treatment (which is
* not cyclic).  To obtain a translation that is correct in all cases, we
* shift the phase correlation cyclically so that a zero translation
* corresponds to the center of the image.  We then correct the coordinates
* returned by local_max_sub_pix.
*由于相位相关是周期性的，负平移将导致图像在下部或者右部出现峰值。
*如果平移在一个或两个方向上接近0，则利用local_max_sub_pix进行插值时将返回错误的值，因为边界问题(不是周期性的）。
*为了获得在任何情况下都是正确的平移信息，循环移动相位相关性，使零平移对应于图像的中心。然后我们修正local_max_sub_pix返回的坐标
RowOffset := Height / 2        //图像Row  中心
ColumnOffset := Width / 2      //图像Column  中心
cyclic_shift_image (ImagePhaseCorrelation, ImageCyclicShift, RowOffset, ColumnOffset)

local_max_sub_pix (ImageCyclicShift, 'facet', 1, 0.02, RowShifted, ColumnShifted)  //图像中局部极大值的亚像素精确检测
Row := RowShifted - RowOffset         //计算出的平移Row
Column := ColumnShifted - ColumnOffset   //计算出的平移Column
dev_set_part (0, 0, Height - 1, Width - 1)
dev_display (ImageDegraded)
disp_message (WindowHandle, 'Translation computed by the phase correlation:\n (' + Row$'5.2f' + ',' + Column$'5.2f' + ')', 'window', 12, 12, 'black', 'true')
gen_cross_contour_xld (Cross, Row, Column, 20, 0.0)   //显示平移后的原点位置
dev_display (Cross)

原图　 　 　 　 　 　 　 　图像平移  　 　 　 　 　 　 　 与二次曲面合成

 　 相位关系图　 　 　 　 　　 　平移位置在圆心  　 　　红十字计算出的偏移位置