OpenCV-Python教程:27.图像转换

理论

傅里叶变换用来分析多种过滤器的频率特征。对于图片,2D离散傅里叶变换(DFT)用来找频率范围。一个快速算法叫快速傅里叶变换(FFT)用来计算DFT。

对于正弦信号,x(t) = Asin(2πft).我们可以说f是信号的频率,如果频率范围给定,我们可以看到f的峰值。如果信号是从离散信号采样,我们还是一样的频率范围,但是是周期范围[-π, π] 或[0, 2π]。你可以认为一个图像是从两个方向采样的信号。所以在X和Y方向做傅里叶变换得到图像的频率表现。

对于正弦信号,如果振幅变化很快,你可以说它是高频信号。如果变化很慢,就是低频信号。这个也可以扩展到图像领域,图像里的振幅变化,在边缘,或者噪点变化最大。所以我们可以说,边缘和噪点是图像里的高频内容。如果在振幅上没特别大变化,就是低频内容

我们来看怎么做傅里叶变换

Numpy 里的傅里叶变换

首先我们来看怎么在Numpy里找傅里叶变换。Numpy有一个FFT包来做这个。np.fft.fft2()让我们做频率变换。第一个参数是输入图像,是灰度的,第二个参数是可选的,决定输出数组的大小。如果比输入图像大的话,输入图像会补全0然后再做转换。如果比输入图像小,输入图像会被裁切,如果不传,输出图像大小和输入一样

当你得到了结果,0频率内容会在左上角,如果你想拿到中间来,你需要N/2在两个方向上来移动结果。这个是用函数np.fft.fftshift().当你找到了频率转换,你可以找到振幅谱线。

import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread('messi5.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

你可以看到白色区域在中心,表示低频内容更多。

现在你达到了频率转换,就可以在频率范围内做更多运算了,比如高通过滤和重建图像,比如反向离散傅里叶变换。你只需要通过一个60x60的矩形窗口去掉低频内容,然后用np.fft.ifftshift()做反向变换,0频率内容又在左上角了。然后做反向快速傅里叶变换 np.ifft2()。结果是复数。

rows, cols = img.shape
crow,ccol = rows/2 , cols/2
fshift[crow-30:crow+30, ccol-30:ccol+30] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.abs(img_back)

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back, cmap ='gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])

plt.show()

结果是

结果显示高通滤波器是一个边缘检测运算,这个我们在图片地图的章节看过,这还显示了大多数凸显该数据是在光谱的低频区域。

如果你仔细看结果,特别是JET色下的图像,你可以看到一些(红色箭头指向的位置),这些波纹状的结构叫做ringing effects(振铃效应)。这是由于我们用来掩图的矩形窗口导致的。所以矩形窗口不用来做过滤,更好的选择是高斯窗口

OpenCV里的傅里叶变换

OpenCV提供了函数cv2.dft()和cv2.idft()。它返回和前面同样的结果但是是双通道的。第一个通道包含了结果的实数部分,第二个通道包含虚数部分,输入图像应该先被转换成np.float32。我们来看看:

import numpy as np
import cv2
from matplotlib import pyplot as plt

img = cv2.imread('messi5.jpg',0)

dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

magnitude_spectrum = 20*np.log(cv2.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

你也可以用cv2.cartToPolar()同时返回magnitude和phase

现在我们改做反向DFT了,在前面的例子里我们创建了HPF,这次我们看如何用LPF来去除图像里的高频内容,它实际上会模糊图片。我们用高值(1)在低频创建一个掩图,传入低频内容,在高频区域用0.

rows, cols = img.shape
crow,ccol = rows/2 , cols/2

# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1

# apply mask and inverse DFT
fshift = dft_shift*mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

注意:


和平常一样,OpenCV函数cv2.dft()和cv2.idft()比Numpy要快,但是Numpy函数更友好。

DFT性能优化

DFT计算的性能在数组大小在有些情况下会好一些,当数组大小是2的乘方的时候是最快的。大小是2,3,5的城际是也处理的比较有效。所以如果你担心你的代码的性能,你可以修改数组的大小到某些优化的大小(通过补0),之后再做DFT,对于OpenCV,你必须手动补0,对于Numpy,你指定FFT计算的大小,它会自动补0.

我们怎么找优化大小呢?OpenCV提供了函数cv2.getOptimalDFTSize()。它可以给cv2.dft()和np.fft.fft2()用,我们可以用IPython魔法命令%timeit

In [16]: img = cv2.imread('messi5.jpg',0)
In [17]: rows,cols = img.shape
In [18]: print rows,cols
342 548

In [19]: nrows = cv2.getOptimalDFTSize(rows)
In [20]: ncols = cv2.getOptimalDFTSize(cols)
In [21]: print nrows, ncols
360 576

(342,548)被修改成了(360,576)。现在用0补全(对于OpenCV)然后看计算DFT的性能。你可以创建一个新的大的全0的数组,然后吧数据拷贝进去,使用cv2.copyMakeBorder()。

nimg=np.zeros((nrows,ncols))
nimg[:rows,:cols]=img

或者

right = ncols - cols
bottom = nrows - rows
bordertype = cv2.BORDER_CONSTANT    #just to avoid line breakup in PDF file
nimg = cv2.copyMakeBorder(img,0,bottom,0,right,bordertype,value=0)

现在我们计算DFT性能和Numpy函数对比

In [22]: %timeit fft1 = np.fft.fft2(img)
10 loops, best of 3: 40.9 ms per loop
In [23]: %timeit fft2 = np.fft.fft2(img,[nrows,ncols])
100 loops, best of 3: 10.4 ms per loop

显示有4倍的性能差别,现在我们用OpenCV的函数

In [24]: %timeit dft1= cv2.dft(np.float32(img),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 13.5 ms per loop
In [27]: %timeit dft2= cv2.dft(np.float32(nimg),flags=cv2.DFT_COMPLEX_OUTPUT)
100 loops, best of 3: 3.11 ms per loop

也是4倍的差别。你还可以看到OpenCV函数比Numpy的函数快3倍。在反向FFT也是一样。

为什么拉普拉斯是高通滤波器?

有很多人提类似的问题,为什么拉普拉斯是一个高通滤波器?为什么Sobel是高通滤波器等,第一个回答是通过傅里叶变换。对于一些高的FFT取拉普拉斯的傅里叶变换,分析它

import cv2
import numpy as np
from matplotlib import pyplot as plt

# simple averaging filter without scaling parameter
mean_filter = np.ones((3,3))

# creating a guassian filter
x = cv2.getGaussianKernel(5,10)
gaussian = x*x.T

# different edge detecting filters
# scharr in x-direction
scharr = np.array([[-3, 0, 3],
[-10,0,10],
[-3, 0, 3]])

# sobel in x direction
sobel_x= np.array([[-1, 0, 1],
[-2, 0, 2],
[-1, 0, 1]])

# sobel in y direction
sobel_y= np.array([[-1,-2,-1],
[0, 0, 0],
[1, 2, 1]])

# laplacian
laplacian=np.array([[0, 1, 0],
[1,-4, 1],
[0, 1, 0]])

filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
filter_name = ['mean_filter', 'gaussian','laplacian', 'sobel_x', \
'sobel_y', 'scharr_x']

fft_filters = [np.fft.fft2(x) for x in filters]
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
mag_spectrum = [np.log(np.abs(z)+1) for z in fft_shift]

for i in xrange(6):
    plt.subplot(2,3,i+1),plt.imshow(mag_spectrum[i],cmap = 'gray')
    plt.title(filter_name[i]), plt.xticks([]), plt.yticks([])

plt.show()

从图像里你能看到每个kernel的屏蔽的频率区域,什么区域被通过了。从这个信息我们可以知道哪个kernel是HPF哪个是LPF。

推荐阅读更多精彩内容