【youcans 的 OpenCV 例程200篇】232. 纹理特征之频谱方法

4.3 纹理特征之频谱方法

傅里叶谱可以描述图像中的周期性或半周期性二维模式的方向性，因此可以基于傅里叶变换对纹理进行频谱分析。

纹理与图像频谱中的高频分量密切相关，纹理模式在频谱图表现为高能量的爆发。

傅里叶谱具有三个特征可以描述纹理的性质：

（1）傅里叶频谱中的突出峰值对应于纹理模式的主方向；
（2）频域图中的峰值位置对应于纹理模式在空间上的基本周期；
（3）通过滤波消除周期分量，用统计方法描述剩下的非周期性信号。

把傅里叶幅度谱转换到极坐标中表示为函数 $S(r,θ)S(r,\theta)$ ，可以简化对频谱特性的解释。S 是频谱函数，r 和 $θ\theta$ 是极坐标系的半径和角度坐标轴。S 对于每一个方向 $θ\theta$ 可以简化为一维函数 $Sθ(r)S_{\theta}(r)$ ；S 对于每一个半径 r 也可以简化为一维函数 $Sr(θ)S_r(\theta)$ 。

分别对一维函数 $Sθ(r)S_{\theta}(r)$ 和 $Sr(θ)S_r(\theta)$ 积分，可以获得纹理频谱的全局描述：
$\sum_{\theta=0}^{\pi} S_{\theta}(r) \\ S(\theta) = \sum_{r=1}^{R} S_r (\theta)$
如果纹理具有空间的周期性或确定的方向性，则一维函数 $S (r)$ 和 $S(θ)S(\theta)$ 在对应的频率具有峰值。

例程 14.13：特征描述之傅里叶频谱纹理分析

本例选自冈萨雷斯《数字图像处理（第四版）》P617 例11.14。一幅包含随机分布目标的图像和另一幅周期性排列目标的图像，在傅里叶变换的幅度谱具有不同的表现形式。

对于随机摆放和整齐摆放的物体图像，使用傅里叶频谱可以显示不同的图像，但直接分析频谱图像比较困难。将其转换到极坐标系中，用一维函数 $S (r)$ 和 $S(θ)S(\theta)$ 表示就非常直观。

对于随机摆放的火柴图像，函数 $S (r)$ 曲线没有很强的周期分量；而对于整齐摆放的火柴图像，函数 $S (r)$ 曲线在r=15、r=25 附近具有很强的峰值，对应于频谱图中的亮区域的周期性水平重复。类似地，随机摆放火柴图像的 $Sr(θ)S_r(\theta)$ 曲线的分布混乱没有显著的峰值，而在整齐摆放的火柴图像中， $Sr(θ)S_r(\theta)$ 曲线在 $θ=90o\theta = 90^o$ 具有非常强烈的能量峰值，对应于频谱图中的水平方向的能量爆发。

    # 14.13 特征描述之纹理谱分析def halfcircle(radius, x0, y0):  # 计算圆心(x0,y0) 半径为 r 的半圆的整数坐标degree = np.arange(180, 360, 1)  # 因对称性可以用半圆 (180,)theta = np.float32(degree * np.pi / 180)  # 弧度，一维数组 (180,)xc = (x0 + radius * np.cos(theta)).astype(np.int)  # 计算直角坐标，整数yc = (y0 + radius * np.sin(theta)).astype(np.int)return xc, ycdef intline(x1, x2, y1, y2):  # 计算从(x1,y1)到(x2,y2)的线段上所有点的坐标dx, dy = np.abs(x2-x1), np.abs(y2-y1)  # x, y 的增量if dx==0 and dy==0:x, y = np.array([x1]), np.array([y1])return x, yif dx > dy:if x1>x2:x1, x2 = x2, x1y1, y2 = y2, y1m = (y2-y1) / (x2-x1)x = np.arange(x1, x2+1, 1)  #[x1,x2]y = (y1 + m*(x-x1)).astype(np.int)else:if y1>y2:x1, x2 = x2, x1y1, y2 = y2, y1m = (x2-x1) / (y2-y1)y = np.arange(y1, y2+1, 1)  # [y1,y2]x = (x1 + m*(y-y1)).astype(np.int)return x, ydef specxture(gray):# cv2.dft 实现图像的傅里叶变换height, width = gray.shapex0, y0 = int(height / 2), int(width / 2)  # x0=300, y0=300rmax = min(height, width) // 2 - 1  # rmax=299print(height, width, x0, y0, rmax)# FFT 变换 （youcans）gray32 = np.float32(gray)  # 将图像转换成 float32dft = cv2.dft(gray32, flags=cv2.DFT_COMPLEX_OUTPUT)  # 傅里叶变换，(600, 600, 2)dftShift = np.fft.fftshift(dft)  # 将低频分量移动到频域图像的中心sAmp = cv2.magnitude(dftShift[:, :, 0], dftShift[:, :, 1])  # 幅度谱，中心化 (600, 600)sAmpLog = np.log10(1 + np.abs(sAmp))  # 幅度谱对数变换 (600, 600)# FFT 频谱沿半径的分布函数sRad = np.zeros((rmax,))  # (299,)sRad[0] = sAmp[x0, y0]for r in range(1, rmax):xc, yc = halfcircle(r, x0, y0)  # 半径为 r 的圆的整数坐标 (360,)sRad[r] = sum(sAmp[xc[i], yc[i]] for i in range(xc.shape[0]))  # (360,)sRadLog = np.log10(1 + np.abs(sRad))  # 极坐标幅度谱 youcans 对数变换# FFT 频谱沿角度的分布函数xmax, ymax = halfcircle(rmax, x0, y0)  # 半径为 xupt 的圆的整数坐标 (360,)sAng = np.zeros((xmax.shape[0],))  # (360,)for a in range(xmax.shape[0]):  # xmax.shape[0]=(360,)xr, yr = intline(x0, xmax[a], y0, ymax[a])  # 从(x0,y0)到(xa,ya)线段所有点的坐标 (300,)sAng[a] = sum(sAmp[xr[i], yr[i]] for i in range(xr.shape[0]))  # (360,)return sAmpLog, sRadLog, sAng# 纹理的傅里叶频谱分析gray1 = cv2.imread("../images/Fig1135a.tif", flags=0)  # flags=0 读取为灰度图像gray2 = cv2.imread("../images/Fig1135b.tif", flags=0)sAmpLog1, sRadLog1, sAng1 = specxture(gray1)  # 图像纹理的频谱分析sAmpLog2, sRadLog2, sAng2 = specxture(gray2)print(sAmpLog1.shape, sRadLog1.shape, sAng1.shape)plt.figure(figsize=(9, 6))plt.subplot(241), plt.axis('off'), plt.title("Random matches"), plt.imshow(gray1, 'gray')plt.subplot(242), plt.axis('off'), plt.title("Amp spectrum"), plt.imshow(sAmpLog1, 'gray')plt.subplot(243), plt.axis('off'), plt.title("Arranged matches"), plt.imshow(gray2, 'gray')plt.subplot(244), plt.axis('off'), plt.title("Amp spectrum"), plt.imshow(sAmpLog2, 'gray')plt.subplot(245), plt.plot(sRadLog1), plt.title("S1 (radius)"), plt.xlim(0,300), plt.yticks([])plt.subplot(246), plt.plot(sAng1), plt.title("S1 (theta)"), plt.xlim(0,180), plt.yticks([])plt.subplot(247), plt.plot(sRadLog2), plt.title("S2 (radius)"), plt.xlim(0,300), plt.yticks([])plt.subplot(248), plt.plot(sAng2), plt.title("S2 (theta)"), plt.xlim(0,180), plt.yticks([])plt.tight_layout()plt.show()