数字图像处理（二）傅里叶变换

目录

• 数字图像处理（二）之傅里叶变换
  • 二维傅里叶变换和逆变换
  • 图像2的整数次幂填充
• 参考

数字图像处理（二）之傅里叶变换

使用python实现数字图像处理中如下功能：

1. 二维傅里叶变换
2. 图像二维傅里叶逆变换
3. 图像中心化傅里叶变换和谱图像
4. 图像2的整数次幂填充

代码链接：xiaohuiduan/image_process (github.com)

二维傅里叶变换和逆变换

图像二维傅里叶变换的公式为，(M,N)为图像的大小，(x,y)为图像的坐标。

[egin{aligned} &F(u, v)=sum_{x=0}^{M-1} sum_{y=0}^{N-1} f(x, y) e^{-j 2 pileft(frac{u x}{M}+frac{v y}{N} ight)} \ &f(x, y)=frac{1}{M N} sum_{u=0}^{M-1} sum_{v=0}^{N-1} F(u, v) e^{j 2 pileft(frac{u x}{M}+frac{v y}{N} ight)} end{aligned} ]

其中，图像的傅里叶变换公式可以等价于：

[egin{aligned} F(u, v) &=sum_{x=0}^{M-1} e^{-j 2 pi frac{u x}{M}} sum_{y=0}^{N-1} f(x, y) e^{-j 2 pi frac{v y}{N}} \ &=sum_{x=0}^{M-1} F(x, v) e^{-j 2 pi frac{u x}{M}} \ F(x, v) &=sum_{y=0}^{N-1} f(x, y) e^{-j 2 pi frac{v y}{N}} end{aligned} ]

也就是说，计算二维傅里叶变换可以先计算行的一维DFT然后再计算列的一维DFT。而一维的傅里叶变换可以使用快速傅里叶变换的方式提高运算速度。因此二维傅里叶变换这样做可以提高计算的速度。

在逆变换中，可以进行如下的变换，将逆变换变成傅立叶正变换的形式。

[egin{aligned} &f(x, y)=frac{1}{M N} sum_{u=0}^{M-1} sum_{v=0}^{N-1} F(u, v) e^{j 2 pileft(frac{u x}{M}+frac{v y}{N} ight)} \ &M N f^{*}(x, y)=sum_{u=0}^{M-1} sum_{v=0}^{N-1} F^{*}(u, v) e^{-j 2 pileft(frac{u x}{M}+frac{v y}{N} ight)} end{aligned} ]

其中(F^*)代表共轭。(f^* = f)（因为(f)是实数，共轭保持不变。）因此：

[f(x, y)=frac{sum_{u=0}^{M-1} sum_{v=0}^{N-1} F^{*}(u, v) e^{-j 2 pileft(frac{u x}{M}+frac{v y}{N} ight)}}{M N} ]

也就是说，可以将傅里叶变换后的结果进行共轭，然后再次进行傅里叶变换并除以(MN)，便可以得到IDFT后的结果。

进行一维傅里叶变换，如果使用公式(F(u)=sum_{x=0}^{M-1} f(x) e^{-j 2 pileft(frac{u x}{M} ight)})进行结算，则时间复杂度为(M^2)，而使用快速傅里叶变换，则时间复杂度为(log_2M)。

一维快速傅里叶变换可以借助numpy中的函数进行实现。下面定义dft2D函数，对图像顺序进行x轴y轴进行一维傅里叶变换，便可以得到二维傅里叶变换。

在图像的傅里叶变换中，如果对傅里叶变换结果进行展示，会发现图像的亮度部分分散在边缘。于是可以对傅里叶变换进行中心化，这样人眼看起来就比较明显。（图像傅里叶的中心化等价于对图像进行image_copy[i][j] = image_copy[i][j]*((-1)**(i+j))）

def dft2D(self, image: np.array, shift: bool = False):
    """二维傅里叶变换，先对X轴进行傅里叶变换，然后再对Y轴进行傅里叶变换。

    Args:
        image (np.array): [图像ORkernel]
        shift (bool, optional): [是否中心化]. Defaults to False.

    Returns:
        [type]: [返回傅里叶变换的值：a+bj]
    """
    image_copy = image.copy()
    if shift:
        for i in range(image_copy.shape[0]):
            for j in range(image_copy.shape[1]):
                image_copy[i][j] = image_copy[i][j]*((-1)**(i+j))

    dft_row = np.fft.fft(image_copy, axis=0)
    dft_2 = np.fft.fft(dft_row, axis=1)
    return dft_2

同时定义IDFT变换函数如下所示，对于IDFT变换我们只需要按照(f(x, y)=frac{sum_{u=0}^{M-1} sum_{v=0}^{N-1} F^{*}(u, v) e^{-j 2 pileft(frac{u x}{M}+frac{v y}{N} ight)}}{M N})来写代码，就很简单了。

def idft2D(self, dft_2: np.array, shift: bool = False):
    """使用dft进行idft变换，F -> F*(共轭变换)

    Args:
        dft_2 (np.array): [傅里叶变换]
        shift (bool,optional): [是否反中心化]. Defaults to False.
    Returns:
        [type]: [image进行反傅里叶变换，可能会产生j虚值。返回幅值]
    """
    dft_2_copy = dft_2.copy()
    idft_row = np.fft.fft(dft_2_copy.conjugate(), axis=0)
    image = np.fft.fft(idft_row, axis=1)
    image = image/(image.shape[0]*image.shape[1])
    if shift:
        for i in range(image.shape[0]):
            for j in range(image.shape[1]):
                image[i][j] = image[i][j]*(-1)**(i+j)
    return abs(image)

下面简单的定义两个画图函数，分别对傅里叶变换后的图像进行画图。

def show_dft2D_in_2D(self, title: str, dft2: np.array, set_log: bool = True):
    """在2维平面上展示傅里叶变换，幅值

    Args:
        title (str): [标题]
        dtf2 (np.array): [傅里叶变换的图像]
        set_log (bool): [对傅里叶变换后的结果取log]
    """
    dft2_copy = dft2.copy()
    dft2_copy = abs(dft2_copy)
    if set_log:
        dft2_copy = np.log2(dft2_copy+1)
    self.show_img(title, [dft2_copy], ['gray'])
    return dft2_copy

def show_dft2D_in_3D(self, title: str, image: np.array, set_log: bool = True):
    """在3维平面上展示傅里叶变换

    Args:
        title (str): [标题]
        image (np.array): [傅里叶变换的图像]
        set_log (bool): [对傅里叶变换后的结果取log]
    """
    image = abs(image)
    if set_log:
        image = np.log10(image+1)
    fig = plt.figure()
    plt.title(title)
    ax3 = plt.axes(projection='3d')

    xx = np.arange(0, image.shape[0])
    yy = np.arange(0, image.shape[1])
    X, Y = np.meshgrid(xx, yy)
    ax3.plot_surface(X, Y, image.astype('uint8'), cmap='rainbow')
    plt.show()

比如说画出(sigma=1)的高斯傅里叶变化的图：

画出rose原图，idft还原的图以及两者差值

图像2的整数次幂填充

在快速傅里叶变换中，假定(N)为2的整数次方。因此如果图像的shape不为2的整数次幂，就需要对其进行填充。

下面是numpy fft函数的一些说明

FFT (Fast Fourier Transform) refers to a way the discrete Fourier Transform (DFT) can be calculated efficiently, by using symmetries in the calculated terms. The symmetry is highest when n is a power of 2, and the transform is therefore most efficient for these sizes.

对于某个数，找到最近的2的整数次幂可以用如下的算法（算法来源于Java HashMap）：

def find2power(self,cap):
    """找到距离cap最近的2的整数次幂,activate by the hashmap in the Java 

    Args:
        cap ([type]): [cap > 0]

    Returns:
        [type]: [2的整数次幂]
    """
    n = cap - 1
    n |= n >> 1	
    n |= n >> 2
    n |= n >> 4
    n |= n >> 8
    n |= n >> 1
    return n+1

因此，对某张图片（灰度图片）进行2的整数次幂的填充操作可以使用如下操作。

image_w = image.shape[0]
image_h = image.shape[1]
fit_w = ImageP.find2power(image_w)
fit_h = ImageP.find2power(image_h)
padding_w = int((fit_w-image_w)/2)
padding_h = int((fit_h - image_h)/2)
image_pad = np.zeros((fit_w,fit_h,1))
image_pad[padding_w:image_w+padding_w,padding_h:image_h+padding_h,:] = image

参考

1. 数字图像处理（第三版）