chatgpt 网站,apache php mysql wordpress,门户网站建设发展趋势,最近新闻热点事件2024专栏地址#xff1a;『youcans 的图像处理学习课』 文章目录#xff1a;『youcans 的图像处理学习课 - 总目录』 【youcans 的 OpenCV 学习课】8. 频率域图像滤波#xff08;上#xff09;
图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制#xff0… 专栏地址『youcans 的图像处理学习课』 文章目录『youcans 的图像处理学习课 - 总目录』 【youcans 的 OpenCV 学习课】8. 频率域图像滤波上
图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制是常用的图像预处理操作。 频率域图像处理先将图像进行傅里叶变换然后在变换域进行处理最后进行傅里叶反变换转换回空间域。 频率域滤波是滤波器传递函数与输入图像的傅里叶变换的对应像素相乘。频率域中的滤波概念更加直观滤波器设计也更容易。使用快速傅里叶变换算法比空间卷积运算速度快很多。 本文提供上述各种算法的完整例程和运行结果并通过一个应用案例示范综合使用多种图像增强方法。 文章目录【youcans 的 OpenCV 学习课】8. 频率域图像滤波上1. 傅里叶变换1.1 傅里叶级数与傅里叶变换例程 8.1连续周期信号的傅立叶级数例程 8.2连续非周期信号的傅立叶系数例程 8.3一维连续函数的傅里叶变换盒式函数1.2 连续函数的取样例程 8.4连续函数的取样1.3 一维离散傅里叶变换例程 8.6一维离散傅里叶变换2. 二维函数的傅里叶变换2.1 二维连续傅里叶变换例程 8.7二维连续函数的傅里叶变换二维盒式函数例程 8.8二维连续函数的傅里叶变换不同参数的盒式函数2.2 图像的混叠和重取样例程 8.9图像的抗混叠2.3 二维离散傅里叶变换DFT2.4 Numpy 实现图像傅里叶变换例程 8.10二维图像的离散傅里叶变换Numpy2.5 OpenCV 实现图像傅里叶变换cv.dft例程 8.11二维图像的离散傅里叶变换OpenCV例程 8.12OpenCV 快速傅里叶变换1. 傅里叶变换
滤波通常是指对图像中特定频率的分量进行过滤或抑制。图像滤波是在尽可能保留图像细节特征的条件下对目标图像的噪声进行抑制是常用的图像预处理操作。
图像滤波不仅可以在空间域进行还可以在频率域进行。空间滤波是图像与各种空间滤波器模板、核的卷积而空间卷积的傅里叶变换是频率域中相应变换的乘积因此频率域滤波是频率域滤波器传递函数与图像的傅里叶变换相乘。
频率域图像滤波先将图像进行傅里叶变换然后在变换域进行处理最后进行傅里叶反变换转换回空间域。
空间域滤波器和频率域滤波器形成一个傅里叶变换对 f(x,y)⊗h(x,y)⇔F(u,v)H(u,v)f(x,y)h(x,y)⇔F(u,v)⊗H(u,v)f(x,y) \otimes h(x,y) \Leftrightarrow F(u,v)H(u,v) \\f(x,y) h(x,y) \Leftrightarrow F(u,v) \otimes H(u,v) f(x,y)⊗h(x,y)⇔F(u,v)H(u,v)f(x,y)h(x,y)⇔F(u,v)⊗H(u,v) 也就是说空间域滤波器和频率域滤波器实际上是相互对应的有些空间域滤波器在频率域通过傅里叶变换实现会更方便、更快速。例如空间域的拉普拉斯滤波器就是频率域的高通滤波器。 1.1 傅里叶级数与傅里叶变换
傅里叶级数Fourier series在数论、组合数学、信号处理、概率论、统计学、密码学、声学、光学等领域都有着广泛的应用。
傅里叶级数公式指出任何周期函数都可以表示为不同频率的正弦函数和/或余弦函数的加权之和 f(t)A0∑n1∞Ansin(nωtψn)A0∑n1∞[ancos(nωt)bnsin(nωt)]\begin{aligned} f(t) A_0 \sum^{\infty}_{n1} A_n sin(n \omega t \psi _n)\\ A_0 \sum^{\infty}_{n1} [a_n cos(n \omega t) b_n sin(n \omega t)] \end{aligned} f(t)A0n1∑∞Ansin(nωtψn)A0n1∑∞[ancos(nωt)bnsin(nωt)] 这个和就是傅里叶级数。
进一步地任何非周期函数也可以表示为不同频率的正弦函数和/或余弦函数乘以加权函数的积分 F(ω)∫−∞∞f(t)e−jωtdtf(t)12π∫−∞∞F(ω)ejωtdω\begin{aligned} F(\omega) \int_{-\infty}^{\infty} f(t) e^{-j\omega t} dt\\ f(t) \frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) e^{j\omega t} d \omega \end{aligned} F(ω)f(t)∫−∞∞f(t)e−jωtdt2π1∫−∞∞F(ω)ejωtdω 这个公式就是傅里叶变换Fourier transform 和逆变换。
*傅里叶变换存在的充分条件是f(t) 的绝对值的积分是有限的在信号处理、图像处理领域这一条件都能满足。 例程 8.1连续周期信号的傅立叶级数 # 8.1连续周期信号的傅立叶级数from scipy import integratenf 30T 10tao 1.0y 1k np.arange(0, nf)an np.zeros(nf)bn np.zeros(nf)amp np.zeros(nf)pha np.zeros(nf)half0, err0 integrate.quad(lambda t: y, -tao/2, tao/2)an[0] 2 * half0 / Tfor n in range(1, nf):half1, err1 integrate.quad(lambda t: 2*y * np.cos(2.0/T * np.pi * n * t), -tao/2, tao/2)an[n] half1 / 10half2, err2 integrate.quad(lambda t: 2*y * np.sin(2.0/T * np.pi * n * t), -tao/2, tao/2)bn[n] half2 / 10amp[n] np.sqrt(an[n]**2 bn[n]**2)for i in range(0, nf):pha[i] 0.0 if an[i]0 else np.piplt.figure(figsize(9, 6))plt.subplot(211), plt.title(Amplitude spectrum), plt.stem(k, amp)plt.subplot(212), plt.title(Phase spectrum), plt.stem(k, pha)plt.show()例程 8.2连续非周期信号的傅立叶系数 # 8.2连续非周期信号的傅立叶系数dx 0.001x np.pi * np.arange(-1dx, 1dx, dx)n len(x)nquart int(np.floor(n/4))f np.zeros_like(x)f[nquart: 2*nquart] (4/n) * np.arange(1, nquart1)f[2*nquart: 3*nquart] np.ones(nquart) - (4/n) * np.arange(0, nquart)plt.figure(figsize(9, 6))plt.title(Fourier series of hat function)plt.plot(x, f, -, colork,labelOriginal)# Compute Fourier seriesA0 np.sum(f * np.ones_like(x)) * dxfFS A0 / 2orders 3A np.zeros(orders)B np.zeros(orders)for k in range(orders):A[k] np.sum(f * np.cos((k1) * x)) * dx # Inner productB[k] np.sum(f * np.sin((k1) * x)) * dxfFS fFS A[k] * np.cos((k 1) * x) B[k] * np.sin((k 1) * x)plt.plot(x, fFS, -, label{} order.format(k))print(Line , k, : A , A[k], B , B[k], fFS.shape)plt.legend(locupper right)plt.show()例程 8.3一维连续函数的傅里叶变换盒式函数
以一维盒式函数为例其傅里叶变换为
F(ω)∫−∞∞f(t)e−jωtdt∫−ππf(t)e−jωtdtsin(ω)ω\begin{aligned} F(\omega) \int_{-\infty}^{\infty} f(t) e^{-j\omega t} dt\\ \int_{-\pi}^{\pi} f(t) e^{-j\omega t} dt\\ \frac{sin(\omega)}{\omega} \end{aligned} F(ω)∫−∞∞f(t)e−jωtdt∫−ππf(t)e−jωtdtωsin(ω) # 8.3一维连续函数的傅里叶变换盒式函数# 盒式函数 (Box_function)x np.arange(-5, 5, 0.1)w 2 * np.pihalfW w / 2y np.where(x, x -halfW, 0)y np.where(x halfW, y, 0)plt.figure(figsize(9, 4))plt.subplot(131), plt.title(Box_function)plt.plot(x, y, -b)plt.subplot(132), plt.title(Fourier transform)fu np.sinc(x)plt.plot(x, fu, -b)plt.subplot(133), plt.title(Spectrum of box_func)fu abs(fu)plt.plot(x, fu, -b)plt.show()1.2 连续函数的取样
连续函数必须经过取样和量化转换为离散函数才能用计算机进行处理。
考虑一个连续函数 f(t)f(t)f(t)以自变量 t 的均匀间隔 ΔT\Delta TΔT 对函数取样取样序列中的任意取样值为
fk∫−∞∞f(t)δ(t−kΔT)dtf(kΔT)f_k \int_{-\infty}^{\infty} f(t) \delta (t-k \Delta T) dt f(k \Delta T) fk∫−∞∞f(t)δ(t−kΔT)dtf(kΔT)
取样后的函数的傅里叶变换为
F~(μ)(F⋆S)(μ)∫−∞∞F(τ)S(μ−τ)dτ1ΔT∑n−∞∞F(μ−nΔT)\begin{aligned} \tilde{F}(\mu) (F \star S) (\mu) \\ \int_{-\infty}^{\infty} F(\tau) S(\mu-\tau) d \tau\\ \frac{1}{\Delta T} \sum_{n-\infty}^{\infty}F(\mu-\frac{n}{\Delta T}) \end{aligned} F~(μ)(F⋆S)(μ)∫−∞∞F(τ)S(μ−τ)dτΔT1n−∞∑∞F(μ−ΔTn)
香农Shannon取样定理指出对于一个连续信号用大于信号最高频率 2倍的取样率来取样则不会丢失信号的有效信息。或者说以 1/ΔT1/\Delta T1/ΔT 的取样率对信号取样点的的最大频率是 μmax1/2ΔT\mu_{max}1/2\Delta Tμmax1/2ΔT 。 例程 8.4连续函数的取样 # 8.4连续函数的取样# 定义函数用于计算所有基矢的频率def gen_freq(N, fs):k np.arange(0, np.floor(N/2) 1, 1)return (k * fs) / NT 100# 定义多个不同频率的基频信号fk [2/T, 5/T, 12/T] # 频率A [7, 3, 2] # 振幅phi [np.pi, 2, 2*np.pi] # 初始相位n np.arange(T)s0 A[0] * np.sin(2 * np.pi * fk[0] * n phi[0])s1 A[1] * np.sin(2 * np.pi * fk[1] * n phi[1])s2 A[2] * np.sin(2 * np.pi * fk[2] * n phi[2])s s0 s1 s2 # 叠加生成混合信号g np.fft.rfft(s) # 傅里叶变换plt.figure(figsize(8, 6))plt.subplot(311)plt.plot(n, s0, n, s1, n, s2, :, marker, alpha0.5)plt.plot(n, s, r-, lw2)plt.title(Sampling of continuous functions)plt.subplot(312)fs 1 # 采样间隔为 1freq gen_freq(T, fsfs) # 计算频率序列ck np.abs(g) / T # 计算振幅plt.plot(freq, ck, .) # 频率-振幅图for f in fk:ck0 round(ck[np.where(freqf*fs)][0], 1)plt.annotate($({},{})$.format(f*fs, ck0), xy(f*fs, ck0), xytext(5, -10), textcoordsoffset points)plt.subplot(313)fs 10 # 采样间隔为 10freq gen_freq(T, fsfs) # 计算频率序列ck np.abs(g) / T # 计算振幅plt.plot(freq, ck, .) # 频率-振幅图for f in fk:ck0 round(ck[np.where(freqf*fs)][0], 1)plt.annotate($({},{})$.format(f*fs, ck0), xy(f*fs, ck0), xytext(5, -10), textcoordsoffset points)plt.show()1.3 一维离散傅里叶变换
数字信号和数字图像都是采样得到的离散变量。
对原函数的变换取样后的数据
F~(μ)∫−∞∞f~(t)e−j2πμtdt∫−∞∞∑n−∞∞f(t)δ(t−nΔT)e−j2πμtdt∑n−∞∞∫−∞∞f(t)δ(t−nΔT)e−j2πμtdt∑n−∞∞fne−j2πμnΔT\begin{aligned} \tilde{F}(\mu) \int_{-\infty}^{\infty} \tilde{f}(t) e^{-j 2 \pi \mu t} dt\\ \int_{-\infty}^{\infty} \sum_{n-\infty}^{\infty} f(t) \delta (t-n{\Delta T}) e^{-j 2 \pi \mu t} dt\\ \sum_{n-\infty}^{\infty} \int_{-\infty}^{\infty} f(t) \delta (t-n{\Delta T}) e^{-j 2 \pi \mu t} dt\\ \sum_{n-\infty}^{\infty} f_n e^{-j 2 \pi \mu n{\Delta T}}\\ \end{aligned} F~(μ)∫−∞∞f~(t)e−j2πμtdt∫−∞∞n−∞∑∞f(t)δ(t−nΔT)e−j2πμtdtn−∞∑∞∫−∞∞f(t)δ(t−nΔT)e−j2πμtdtn−∞∑∞fne−j2πμnΔT
当取样频率为 μm/MΔT\mu m/M \Delta Tμm/MΔT 时可以得到离散傅里叶变换DFT和逆变换公式为
Fm∑n0M−1fne−j2πμmn/M,m0,...M−1fn1M∑m0M−1Fmej2πμmn/M,n0,...M−1\begin{aligned} F_m \sum_{n0}^{M-1} f_n e^{-j\ 2\pi \mu mn/M}, m0,...M-1\\ f_n \frac{1}{M} \sum_{m0}^{M-1} F_m e^{j\ 2\pi \mu mn/M}, n0,...M-1 \end{aligned} Fmfnn0∑M−1fne−j 2πμmn/M,M1m0∑M−1Fmej 2πμmn/M,m0,...M−1n0,...M−1
由于任何周期或非周期函数都可以表示为不同频率的正弦函数和余弦函数之和的形式因此用傅里叶变换表示的函数特征完全可以通过傅里叶反变换来重建而且不会丢失任何信息。这是频率域图像处理的数学基础。
离散傅里叶变换 是将离散信号分解为一系列离散三角函数分量每个三角函数分量都有对应的幅值 A、频率 f 和相位 φ\varphiφ。通过所有分量叠加可以得到原离散信号。
在图像处理中通常使用 x,yx, yx,y 表示离散的空间域坐标变量用 u,vu,vu,v 表示离散的频率域变量。于是将一维离散傅里叶变换表示为
F(u)∑x0M−1f(x)e−j2πux/M,u0,...M−1f(x)1M∑u0M−1F(u)ej2πux/M,x0,...M−1\begin{aligned} F(u) \sum_{x0}^{M-1} f(x) e^{-j\ 2\pi u x/M}, u0,...M-1\\ f(x) \frac{1}{M} \sum_{u0}^{M-1} F(u) e^{j\ 2\pi u x/M}, x0,...M-1 \end{aligned} F(u)f(x)x0∑M−1f(x)e−j 2πux/M,M1u0∑M−1F(u)ej 2πux/M,u0,...M−1x0,...M−1 例程 8.6一维离散傅里叶变换 # 8.6一维离散傅里叶变换# 生成方波信号N 200t np.linspace(-10, 10, N)dt t[1] - t[0]sint np.sin(t)sig np.sign(sint)fig plt.figure(figsize(10, 4))plt.subplot(131), plt.title(source), plt.xticks([]), plt.yticks([])plt.plot(t, sig)# 离散傅里叶变换fft np.fft.fft(sig, N) # 离散傅里叶变换fftshift np.fft.fftshift(fft) # 对称平移amp abs(fftshift) / len(fft) # 幅值pha np.angle(fftshift) # 相位fre np.fft.fftshift(np.fft.fftfreq(ddt, nN)) # 频率plt.subplot(132), plt.title(DFT), plt.xticks([]), plt.yticks([])plt.plot(t, fft)# 信号恢复recover np.zeros(N)for a, p, f in zip(amp, pha, fre):sigCos a * np.cos(2 * np.pi * f * np.arange(N) * dt p) # 根据幅度,相位,频率构造三角函数recover sigCos # 把这些三角函数都加起来plt.subplot(133), plt.title(recover), plt.xticks([]), plt.yticks([])plt.plot(t, recover)plt.show()2. 二维函数的傅里叶变换
2.1 二维连续傅里叶变换
设 f(t,z)f(t,z)f(t,z) 是二维连续变量 t,zt, zt,z 的连续函数则其二维连续傅里叶变换和逆变换为
F(μ,ν)∫−∞∞∫−∞∞f(t,z)e−j2π(μtνz)dtdzf(t,z)∫−∞∞∫−∞∞F(μ,ν)ej2π(μtνz)dμdν\begin{aligned} F(\mu,\nu) \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(t,z) e^{-j 2\pi (\mu t \nu z)} dt \ dz\\ f(t,z) \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(\mu,\nu) e^{j 2\pi (\mu t \nu z)} d\mu \ d\nu \end{aligned} F(μ,ν)f(t,z)∫−∞∞∫−∞∞f(t,z)e−j2π(μtνz)dt dz∫−∞∞∫−∞∞F(μ,ν)ej2π(μtνz)dμ dν 对于图像处理式中的 t,zt, zt,z 表示连续空间变量 μ,ν\mu, \nuμ,ν 表示连续频率变量。 例程 8.7二维连续函数的傅里叶变换二维盒式函数
以盒式函数为例其傅里叶变换为
F(μ,ν)∫−∞∞∫−∞∞f(t,z)e−j2π(μtνz)dtdz∫−T/2T/2∫−Z/2Z/2Ae−j2π(μtνz)dtdzATZ[sin(πμT)πμT][sin(πνZ)πνZ]\begin{aligned} F(\mu,\nu) \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(t,z) e^{-j 2\pi (\mu t \nu z)} dt \ dz\\ \int_{-T/2}^{T/2} \int_{-Z/2}^{Z/2} A e^{-j 2\pi (\mu t \nu z)} dt \ dz\\ ATZ[\frac{sin(\pi \mu T)}{\pi \mu T}][\frac{sin(\pi \nu Z)}{\pi \nu Z}] \end{aligned} F(μ,ν)∫−∞∞∫−∞∞f(t,z)e−j2π(μtνz)dt dz∫−T/2T/2∫−Z/2Z/2Ae−j2π(μtνz)dt dzATZ[πμTsin(πμT)][πνZsin(πνZ)] # 8.7二维连续函数的傅里叶变换二维盒式函数# 二维盒式函数 (2D-Box_function)t np.arange(-5, 5, 0.1) # start, end, stepz np.arange(-5, 5, 0.1)height, width len(t), len(z)tt, zz np.meshgrid(t, z) # 将一维数组 xnew, ynew 转换为网格点集二维数组f np.zeros([len(t), len(z)]) # 初始化置零f[30:70, 30:70] 1 # 二维盒式函数fu np.sinc(t)fv np.sinc(z)fuv np.outer(np.sinc(t), np.sinc(z)) # 由公式计算连续傅里叶变换print(fu.shape, fv.shape, fuv.shape)fig plt.figure(figsize(10, 6))ax1 plt.subplot(121, projection3d)ax1.set_title(2D-Box function)ax1.plot_wireframe(tt, zz, f, rstride2, cstride2, linewidth0.5, colorc)ax1.set_xticks([]), ax1.set_yticks([]), ax1.set_zticks([])ax2 plt.subplot(122, projection3d)ax2.set_title(2D-Fourier transform)ax2.plot_wireframe(tt, zz, fuv, rstride2, cstride2, linewidth0.5, colorc)ax2.set_xticks([]), ax2.set_yticks([]), ax2.set_zticks([])plt.show()例程 8.8二维连续函数的傅里叶变换不同参数的盒式函数 # 8.8二维连续函数的傅里叶变换不同参数的盒式函数# 二维盒式函数 (2D-Box_function)height, width 128, 128m int((height - 1) / 2)n int((width - 1) / 2)fig plt.figure(figsize(9, 6))T [5, 10, 20]print(len(T))for i in range(len(T)):f np.zeros([height, width]) # 初始化置零f[m-T[i]:mT[i], n-T[i]:nT[i]] 1 # 不同尺寸的盒式函数fft np.fft.fft2(f)shift_fft np.fft.fftshift(fft)amp np.log(1 np.abs(shift_fft))plt.subplot(2,len(T),i1), plt.xticks([]), plt.yticks([])plt.imshow(f, gray), plt.title(Box Filter (T{}).format(T[i]))#plt.subplot(2,len(T),len(T)i1), plt.xticks([]), plt.yticks([])plt.imshow(amp, gray), plt.title(FFT Spectrum)plt.tight_layout()plt.show()2.2 图像的混叠和重取样
由于无法对一个函数无限地取样因此在数字图像中总是会出现混叠。
混叠分为空间混叠和时间混叠。空间混叠是由欠取样造成的在陌生重复的图像中较为明显。时间混叠与动态图像序列中图像的时间间隔相关如辐条倒转的“车轮效应”。
空间混叠的主要问题是会引入伪影即原始图像中未出现的线条锯齿、虚假高光和频率模式。
对图像进行缩放都会引入混叠。图像放大可以视为过取样使用双线性、双三次内插可以降低图像放大中的混叠。图像缩小可以视为欠取样混叠通常更为严重。
为了降低混叠在重取样前要使用低通滤波器来平滑以衰减图像的高频分量可以有效地抑制混叠但图像的清晰度也有所下降。 例程 8.9图像的抗混叠 # 8.9图像的抗混叠imgGray cv2.imread(../images/Fig0417a.tif, flags0) # flags0 读取为灰度图像height, width imgGray.shape[:2] # 图片的高度和宽度# 先缩小图像再放大图像imgZoom1 cv2.resize(imgGray, (int(0.33*width), int(0.33*height)))imgZoom2 cv2.resize(imgZoom1, None, fx3, fy3, interpolationcv2.INTER_AREA)# 先对原图像进行5x5的低通滤波再缩小图像再放大图像kSize (5, 5)imgBoxFilter cv2.boxFilter(imgGray, -1, kSize) # cv2.boxFilter 方法imgZoom3 cv2.resize(imgBoxFilter, (int(0.33*width), int(0.33*height)))imgZoom4 cv2.resize(imgZoom3, None, fx3, fy3, interpolationcv2.INTER_AREA)fig plt.figure(figsize(9,5))plt.subplot(131), plt.axis(off), plt.title(Origin)plt.imshow(imgGray, cmapgray)plt.subplot(132), plt.axis(off), plt.title(Resample)plt.imshow(imgZoom2, cmapgray)plt.subplot(133), plt.axis(off), plt.title(Box Resample)plt.imshow(imgZoom4, cmapgray)plt.tight_layout()plt.show()2.3 二维离散傅里叶变换DFT
对于二维图像处理通常使用 x,yx, yx,y 表示离散的空间域坐标变量用 u,vu,vu,v 表示离散的频率域变量。二维离散傅里叶变换DFT和反变换IDFT为
F(u,v)∑x0M−1∑y0N−1f(x,y)e−j2π(ux/Mvy/N)f(x,y)1MN∑u0M−1∑v0N−1F(u,v)ej2π(ux/Mvy/N)\begin{aligned} F(u,v) \sum_{x0}^{M-1} \sum_{y0}^{N-1} f(x,y) e^{-j 2\pi (ux/Mvy/N)}\\ f(x,y) \frac{1}{MN} \sum_{u0}^{M-1} \sum_{v0}^{N-1} F(u,v) e^{j 2\pi (ux/Mvy/N)} \end{aligned} F(u,v)f(x,y)x0∑M−1y0∑N−1f(x,y)e−j2π(ux/Mvy/N)MN1u0∑M−1v0∑N−1F(u,v)ej2π(ux/Mvy/N) 二维离散傅里叶变换也可以用极坐标表示 F(u,v)R(u,v)jI(u,v)∣F(u,v)∣ejϕ(u,v)F(u,v) R(u,v) j I(u,v) |F(u,v)| e^{j \phi (u,v)} F(u,v)R(u,v)jI(u,v)∣F(u,v)∣ejϕ(u,v) 傅里叶频谱Fourier spectrum为 ∣F(u,v)∣[R2(u,v)I2(u,v)]1/2|F(u,v)| [R^2(u,v) I^2(u,v)]^{1/2} ∣F(u,v)∣[R2(u,v)I2(u,v)]1/2 傅里叶相位谱Fourier phase spectrum为 ϕ(u,v)arctan[I(u,v)/R(u,v)]\phi (u,v) arctan[I(u,v)/R(u,v)] ϕ(u,v)arctan[I(u,v)/R(u,v)] 傅里叶功率谱Fourier power spectrum为 P(u,v)∣F(u,v)∣2R2(u,v)I2(u,v)P(u,v) |F(u,v)|^2 R^2(u,v) I^2(u,v) P(u,v)∣F(u,v)∣2R2(u,v)I2(u,v)
空间取样和频率间隔是相互对应的频率域所对应的离散变量间的间隔为Δu1/MΔTΔv1/NΔZ\Delta u 1/M \Delta T\Delta v 1/N \Delta ZΔu1/MΔTΔv1/NΔZ。即频域中样本之间的间隔与空间样本之间的间隔及样本数量的乘积成反比。
空间域滤波器和频率域滤波器也是相互对应的二维卷积定理是在空间域和频率域滤波之间建立等价关系的纽带 (f⋆h)(x,y)⇔(F⋅H)(u,v)(f \star h)(x,y) \Leftrightarrow (F \cdot H)(u,v) (f⋆h)(x,y)⇔(F⋅H)(u,v) 这表明 F 和 H 分别是 f 和 h 的傅里叶变换f 和 h 的空间卷积的傅里叶变换是它们的变换的乘积。
因此计算两个函数的空间卷积可以直接在空间域计算也可以在频率域计算先计算每个函数的傅里叶变换再对两个变换相乘最后进行傅里叶反变换转换回空间域。
也就是说空间域滤波器和频率域滤波器实际上是相互对应的有些空间域滤波器在频率域通过傅里叶变换实现会更方便、更快速。 2.4 Numpy 实现图像傅里叶变换
Numpy 中的 np.fft.fft2() 函数可以实现图像的傅里叶变换 。
函数说明 numpy.fft.fft2(a, sNone, axes(- 2, - 1), normNone) → out参数说明
a输入数组一维或多维数组可以是复数形式s输出数组的形状每个变换轴的长度对应于 fft(x,n) 中的 n可选项out输出数组复数形式的一维或多维数组complex ndarray
注意事项
用于二维图像傅里叶变换时输入数组 a 是 numpy 二维数组灰度图像或三维数组彩色图像。输出结果是复数形式 (Realj∗Imag)(Real j * Imag)(Realj∗Imag) 的数组不能直接用于显示图像。为了显示傅里叶变换结果的图像需要将数组的值调整到 [0,255] 的灰度空间内。
经过 np.fft.fft2() 函数实现傅里叶变换得到的图片频谱信息幅度谱的最大值低频分量在左上角 (0,0) 处。为了便于观察通常用 np.fft.fftshift() 函数将低频分量移动到频域图像的中心位置。
函数说明 numpy.fft.fftshift(x, axesNone) → y参数说明
x输入数组一维或多维数组axes整数或输入数组形状的元组用于指定移动的轴可选项y输出数组 例程 8.10二维图像的离散傅里叶变换Numpy # 8.10Numpy 实现二维离散傅里叶变换normalize lambda x: (x - x.min()) / (x.max() - x.min() 1e-6)imgGray cv2.imread(../images/Fig0424a.tif, flags0) # flags0 读取为灰度图像# imgGray cv2.imread(../images/imgBall.png, flags1) # flags0 读取为灰度图像# 傅里叶变换# fft np.fft.fft2(imgGray.astype(np.float32))fft np.fft.fft2(imgGray) # np.fft.fft2 实现傅里叶变换# 非中心化计算幅度谱和相位谱ampSpectrum np.sqrt(np.power(fft.real, 2) np.power(fft.imag, 2)) # 幅度谱print(ampSpectrum max{}, min{}.format(ampSpectrum.max(), ampSpectrum.min()))# phase np.arctan2(fft.imag, fft.real) # 计算相位角(弧度制)# phiSpectrum phase / np.pi*180 # 将相位角转换为 [-180, 180]phiSpectrum np.angle(fft)# 中心化将低频分量移动到频域图像的中心fftShift np.fft.fftshift(fft) # 将低频分量移动到频域图像的中心# 中心化后的幅度谱ampSpeShift np.sqrt(np.power(fftShift.real, 2) np.power(fftShift.imag, 2))ampShiftNorm np.uint8(normalize(ampSpeShift)*255) # 归一化为 [0,255]# 幅度谱做对数变换ampSpeLog np.log(1 ampSpeShift) # 幅度谱做对数变换以便于显示ampSpeLog np.uint8(normalize(ampSpeLog)*255) # 归一化为 [0,255]# np.fft.ifft2 实现图像的逆傅里叶变换invShift np.fft.ifftshift(fftShift) # 将低频逆转换回图像四角imgIfft np.fft.ifft2(invShift) # 逆傅里叶变换返回值是复数数组imgRebuild np.abs(imgIfft) # 将复数数组调整至灰度空间plt.figure(figsize(9, 6))plt.subplot(231), plt.title(Original image), plt.axis(off)plt.imshow(imgGray, cmapgray)plt.subplot(232), plt.title(FFT phase spectrum), plt.axis(off)plt.imshow(phiSpectrum, cmapgray)plt.subplot(233), plt.title(Rebuild image with IFFT), plt.axis(off)plt.imshow(imgRebuild, cmapgray)plt.subplot(234), plt.title(FFT amplitude spectrum), plt.axis(off)plt.imshow(ampSpectrum, cmapgray)plt.subplot(235), plt.title(FFT-shift amplitude), plt.axis(off)plt.imshow(ampSpeShift, cmapgray)plt.subplot(236), plt.title(Log-trans of FFT amp), plt.axis(off)plt.imshow(ampSpeLog, cmapgray)plt.tight_layout()plt.show()程序说明
图中未中心化的幅度谱FFT amp spe也并不是完全黑色在图像的四角位置都有微小的亮区域但是很难观察到这也是对幅度谱进行中心化处理fftShift的原因。 2.5 OpenCV 实现图像傅里叶变换cv.dft
使用 OpenCV 中的 cv.dft() 函数也可以实现图像的傅里叶变换cv.idft() 函数实现图像傅里叶逆变换。
函数说明 cv.dft(src[, dst[, flags[, nonzeroRows]]]) → dstcv.idft(src[, dst[, flags[, nonzeroRows]]]) → dst参数说明
src输入图像单通道灰度图像使用 np.float32 格式dst输出图像图像大小与 src 相同数据类型由 flag 决定flag转换标识符 cv.DFT_INVERSE用一维或二维逆变换取代默认的正向变换cv.DFT_SCALE缩放比例标识根据元素数量求出缩放结果常与DFT_INVERSE搭配使用cv.DFT_ROWS: 对输入矩阵的每行进行正向或反向的傅里叶变换常用于三维或高维变换等复杂操作cv.DFT_COMPLEX_OUTPUT对一维或二维实数数组进行正向变换默认方法结果是由 2个通道表示的复数阵列第一通道是实数部分第二通道是虚数部分cv.DFT_REAL_OUTPUT对一维或二维复数数组进行逆变换结果通常是一个尺寸相同的复数矩阵
注意事项
输入图像 src 是 np.float32 格式如图像使用 np.uint8 格式则必须先转换 np.float32 格式。默认方法 cv.DFT_COMPLEX_OUTPUT 时输入 src 是 np.float32 格式的单通道二维数组输出 dst 是 2个通道的二维数组第一通道 dft[:,:,0] 是实数部分第二通道 dft[:,:,1] 是虚数部分。不能直接用于显示图像。可以使用 cv.magnitude() 函数将傅里叶变换的结果转换到灰度 [0,255]。idft(src, dst, flags) 等价于 dft(src, dst, flagsDFT_INVERSE)。OpenCV 实现傅里叶变换计算速度比 Numpy 更快。
转换标识符为 cv.DFT_COMPLEX_OUTPUT 时cv.dft() 函数的输出是 2个通道的二维数组使用 cv.magnitude() 函数可以实现计算二维矢量的幅值 。
函数说明 cv.magnitude(x, y[, magnitude]) → dst参数说明
x一维或多维数组也表示复数的实部浮点型y一维或多维数组也表示复数的虚部浮点型数组大小必须与 x 相同dst输出数组数组大小和数据类型与 x 相同运算公式为
dst(I)x(I)2y(I)2dst(I) \sqrt{x(I)^2 y(I)^2} dst(I)x(I)2y(I)2
傅里叶变换及相关操作的取值范围可能不适于图像显示需要进行归一化处理。 OpenCV 中的 cv.normalize() 函数可以实现图像的归一化。
函数说明 cv.normalize(src, dst[, alpha[, beta[, norm_type[, dtype[, mask]]]]]) → dst参数说明
src输入图像dst输出结果与输入图像同尺寸同类型alpha归一化后的最小值可选项默认值为0beta归一化后的最大值可选项默认值为1norm_type归一化类型 NORM_INFLinf 范数绝对值的最大值NORM_L1L1 范数绝对值的和NORM_L2L2 范数欧几里德距离默认类型NORM_MINMAX线性缩放常用类型 dtype可选项默认值 -1表示输出矩阵与输入图像类型相同mask掩模遮罩可选项默认无遮罩
傅里叶变换在理论上需要 O(MN)2O(MN)^2O(MN)2 次运算非常耗时快速傅里叶变换只需要 O(MNlog(MN))O(MN log (MN))O(MNlog(MN)) 次运算就可以完成。
OpenCV 中的傅里叶变换函数 cv.dft() 对于行数和列数都可以分解为 2p∗3q∗5r2^p * 3^q * 5^r2p∗3q∗5r 的矩阵的计算性能最好。为了提高运算性能可以对原矩阵的右侧和下方补 0以满足该分解条件。OpenCV 中的 cv.getOptimalDFTSize() 函数可以实现图像的最优 DFT 尺寸扩充适用于 cv.dft() 和 np.fft.fft2()。
函数说明 cv.getOptimalDFTSize(versize) → retval参数说明
versize数组大小retvalDFT 扩充的最优数组大小 例程 8.11二维图像的离散傅里叶变换OpenCV # 8.11OpenCV 实现二维图像的离散傅里叶变换imgGray cv2.imread(../images/Fig0424a.tif, flags0) # flags0 读取为灰度图像# cv2.dft 实现图像的傅里叶变换imgFloat32 np.float32(imgGray) # 将图像转换成 float32dft cv2.dft(imgFloat32, flagscv2.DFT_COMPLEX_OUTPUT) # 傅里叶变换dftShift np.fft.fftshift(dft) # 将低频分量移动到频域图像的中心# 幅度谱# ampSpe np.sqrt(np.power(dft[:,:,0], 2) np.power(dftShift[:,:,1], 2))dftAmp cv2.magnitude(dft[:,:,0], dft[:,:,1]) # 幅度谱未中心化dftShiftAmp cv2.magnitude(dftShift[:,:,0], dftShift[:,:,1]) # 幅度谱中心化dftAmpLog np.log(1 dftShiftAmp) # 幅度谱对数变换以便于显示# 相位谱phase np.arctan2(dftShift[:,:,1], dftShift[:,:,0]) # 计算相位角(弧度制)dftPhi phase / np.pi*180 # 将相位角转换为 [-180, 180]print(dftMag max{}, min{}.format(dftAmp.max(), dftAmp.min()))print(dftPhi max{}, min{}.format(dftPhi.max(), dftPhi.min()))print(dftAmpLog max{}, min{}.format(dftAmpLog.max(), dftAmpLog.min()))# cv2.idft 实现图像的逆傅里叶变换invShift np.fft.ifftshift(dftShift) # 将低频逆转换回图像四角imgIdft cv2.idft(invShift) # 逆傅里叶变换imgRebuild cv2.magnitude(imgIdft[:,:,0], imgIdft[:,:,1]) # 重建图像plt.figure(figsize(9, 6))plt.subplot(231), plt.title(Original image), plt.axis(off)plt.imshow(imgGray, cmapgray)plt.subplot(232), plt.title(DFT Phase), plt.axis(off)plt.imshow(dftPhi, cmapgray)plt.subplot(233), plt.title(Rebuild image with IDFT), plt.axis(off)plt.imshow(imgRebuild, cmapgray)plt.subplot(234), plt.title(DFT amplitude spectrum), plt.axis(off)plt.imshow(dftAmp, cmapgray)plt.subplot(235), plt.title(DFT-shift amplitude), plt.axis(off)plt.imshow(dftShiftAmp, cmapgray)plt.subplot(236), plt.title(Log-trans of DFT amp), plt.axis(off)plt.imshow(dftAmpLog, cmapgray)plt.tight_layout()plt.show()例程 8.12OpenCV 快速傅里叶变换 # 8.12OpenCV 快速傅里叶变换imgGray cv2.imread(../images/Fig0429a.tif, flags0) # flags0 读取为灰度图像rows,cols imgGray.shape[:2] # 图像的行(高度)/列(宽度)# 快速傅里叶变换(要对原始图像进行矩阵扩充)rPad cv2.getOptimalDFTSize(rows) # 最优 DFT 扩充尺寸cPad cv2.getOptimalDFTSize(cols) # 用于快速傅里叶变换imgEx np.zeros((rPad,cPad,2),np.float32) # 对原始图像进行边缘扩充imgEx[:rows,:cols,0] imgGray # 边缘扩充下侧和右侧补0dftImgEx cv2.dft(imgEx, cv2.DFT_COMPLEX_OUTPUT) # 快速傅里叶变换# 傅里叶逆变换idftImg cv2.idft(dftImgEx) # 逆傅里叶变换idftMag cv2.magnitude(idftImg[:,:,0], idftImg[:,:,1]) # 逆傅里叶变换幅值# 矩阵裁剪得到恢复图像idftMagNorm np.uint8(cv2.normalize(idftMag, None, 0, 255, cv2.NORM_MINMAX)) # 归一化为 [0,255]imgRebuild np.copy(idftMagNorm[:rows, :cols])print(max(imgGray-imgRebuild) , np.max(imgGray-imgRebuild))print(imgGray:{}, dftImg:{}, idftImg:{}, imgRebuild:{}.format(imgGray.shape, dftImgEx.shape, idftImg.shape, imgRebuild.shape))plt.figure(figsize(9, 6))plt.subplot(131), plt.title(Original image), plt.axis(off)plt.imshow(imgGray, cmapgray)plt.subplot(132), plt.title(Log-trans of DFT amp), plt.axis(off)dftAmp cv2.magnitude(dftImgEx[:,:,0], dftImgEx[:,:,1]) # 幅度谱中心化dftAmpLog np.log(1 dftAmp) # 幅度谱对数变换以便于显示plt.imshow(cv2.normalize(dftAmpLog, None, 0, 255, cv2.NORM_MINMAX), cmapgray)plt.subplot(133), plt.title(Rebuild image with IDFT), plt.axis(off)plt.imshow(imgRebuild, cmapgray)plt.tight_layout()plt.show()版权声明
本文中部分原始图片来自 Rafael C. Gonzalez “Digital Image Processing, 4th.Ed.”特此致谢。 youcansxupt 原创作品转载必须标注原文链接(https://blog.csdn.net/youcans/article/details/120995650) Copyright 2022 youcans, XUPT 欢迎关注 『youcans 的 OpenCV 学习课』 系列持续更新 【youcans 的 OpenCV 学习课】1. 安装与环境配置 【youcans 的 OpenCV 学习课】2. 图像读取与显示 【youcans 的 OpenCV 学习课】3. 图像的创建与修改 【youcans 的 OpenCV 学习课】4. 图像的叠加与混合 【youcans 的 OpenCV 学习课】5. 图像的几何变换 【youcans 的 OpenCV 学习课】6. 灰度变换与直方图处理 【youcans 的 OpenCV 学习课】7. 空间域图像滤波 【youcans 的 OpenCV 学习课】8. 频率域图像滤波上 【youcans 的 OpenCV 学习课】9. 频率域图像滤波下
