傅立叶-梅林转换

傅立叶-梅林转换（Fourier–Mellin Transform, FMT）是一种结合傅立叶变换与梅林变换概念的积分变换方法，主要用于影像处理中获得对旋转与尺度变化不敏感的特征表示。傅立叶-梅林转换常用于解决影像间的旋转不变性与尺度不变性问题，例如在模板匹配、模式识别、影像注册（配准）以及数位浮水印等应用中，都利用其对旋转与缩放的稳健特性。需要注意的是，傅立叶-梅林转换本身并未直接提供平移不变性；影像的平移（位移差异）通常需要透过额外的相位相关（phase correlation）技术来估计处理。

## 定义

傅立叶–梅林转换（Fourier–Mellin Transform, FMT）是对以极座标表示之函数进行傅立叶与梅林变换的积分变换。令函数以极座标表示为 $f_p(r,\theta)$（其中 $r\ge0$ 为径向座标，$\theta$为角度），则傅立叶–梅林转换定义为：

$$M_f(u,v)=\frac{1}{2\pi}\int_0^{\infty}\int_0^{2\pi}f_p(r,\theta)\,r^{-ju}\,e^{-jv\theta}\,d\theta\,\frac{dr}{r}\,.$$

令$$\rho=\ln r,$$ 则有$r=e^{\rho}$ 且 $\frac{dr}{r}=d\rho$。因此，上式可改写为：

$$M_f(u,v)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\int_0^{2\pi}f_p(e^{\rho},\theta)\,e^{-ju\rho}\,e^{-jv\theta}\,d\theta\,d\rho\,.$$

这也就是说对于一个在笛卡尔座标系统的函数$f_c(x,y)$，只要经过指数-极座标转换(log-polar transformation) 后，就可用傅立叶转换来得到$f_c(x,y)$的傅立叶-梅林转换。

## 性质

### 旋转不变性

若令$$f'_p(r,\theta)=f_p(r,\theta-\Delta\theta),$$ 即原函数旋转 $\Delta\theta$ 后，则傅立叶–梅林转换为$$M_{f'}(u,v)=e^{-jv\Delta\theta}\,M_f(u,v)\,.$$ 这表示傅立叶–梅林谱的震幅保持不变，仅在角度方向上乘上一个相位因子 $e^{-ik\Delta\theta}$。

### 缩放不变性

若令$$f'_p(r,\theta)=f_p(ar,\theta),\quad a>0,$$ 则可推导出$$M_{f'}(u,v)=a^{-ju}\,M_f(u,v) =\;e^{-ju\ln a}\,M_f(u,v)\,.$$ 也就是说，尺度变化仅使傅立叶–梅林谱在径向（对数）方向上乘上一个相位因子 $e^{-ju\ln a}$，而震幅不变。

## 原理与实现流程

在傅立叶-梅林转换的实作中，通常采用以下两个步骤来将原始函数（或影像震幅）转换到傅立叶-梅林域：

- 指数-极座标转换

先将以笛卡尔座标表示的函数 $f_c(x,y)$ 转换成以极座标表示的形式。令影像中心为 $(c_x,c_y)$，对于任意像素 $(x,y)$，计算其径向距离与角度：

$$r=\sqrt{(x-c_x)^2+(y-c_y)^2},\quad \theta=\arctan\frac{y-c_y}{\,x-c_x\,}.$$ 接着，将径向座标取对数，令$$\rho=\ln r,$$ 则原函数可写为 $f_p(e^{\rho},\theta)$。在这个新表示中，尺度（缩放）变化会转化为 $\rho$ 轴上的平移，而旋转则表现为 $\theta$ 轴上的平移。

- 傅立叶转换

对经过指数-极座标转换后的函数 $f(e^{\rho},\theta)$，在 $\rho$ 与 $\theta$ 方向上进行二维傅立叶变换，即可得到傅立叶-梅林表示：

$$M_f(u,v)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\int_0^{2\pi}f_p(e^{\rho},\theta)\,e^{-ju\rho}\,e^{-jv\theta}\,d\theta\,d\rho\,.$$

【延伸说明】

在实际影像应用中，常将原始影像先经傅立叶变换取得震幅图(忽略相位移)，再对该震幅图执行指数-极座标转换，最后利用快速傅立叶转换(FFT)计算傅立叶-梅林表示；另外，针对平移不变性部分，通常会结合相位相关法来估计平移量，此处不在本定义范畴内。

## 应用

傅立叶-梅林转换因其对旋转与尺度变化具有不变性，故在下列领域中广泛应用：

- 影像配准：例如 Reddy 与 Chatterji (1996) 提出的方法利用傅立叶-梅林转换校正影像间的旋转与尺度差异，并搭配相位相关法估计平移。
- 视觉里程计：Kazik 与 Göktoğan (2011) 采用傅立叶-梅林转换估计单目影像序列中的旋转与尺度变化，以计算机器人运动参数。
- 光流估计：Ho 与 Goecke (2008) 探讨利用傅立叶-梅林转换估计带旋转运动的影像块运动。
- 数位浮水印：O' Ruanaidh 与 Pun (1997) 利用傅立叶-梅林表示实现对几何变换具有鲁棒性的水印嵌入。

## Python 实现范例

以下示例程式码展示如何以 Python 实现傅立叶-梅林转换的正向与反向过程。此处假设输入影像已为震幅图，流程为：将震幅图映射至对数极座标 → 对映射结果执行二维傅立叶变换；反向转换则为逆过程。

```python
import cv2
import numpy as np

def fmt_forward(image):
    """
    計算影像的傅立葉-梅林正轉換表示。
    輸入: 灰階影像（震幅圖，numpy array）。
    輸出: 複數型傅立葉-梅林表示 M_f(u,v)。
    """
    # 1. 對數極座標映射：將震幅圖映射至對數極座標系統
    M, N = image.shape
    center = (N/2, M/2)
    maxRadius = np.hypot(center[0], center[1])
    log_polar = cv2.warpPolar(image.astype(np.float32), (N, M), center, maxRadius, cv2.WARP_POLAR_LOG)
    
    # 2. 對映射後的圖像執行二維傅立葉變換
    T = np.fft.fft2(log_polar)
    return T

def fmt_inverse(T, image_shape):
    """
    將傅立葉-梅林表示反轉換為對數極座標域中的圖像。
    再進行逆對數極座標映射以重建原始震幅圖。
    輸入: T - 傅立葉-梅林表示；image_shape - 原始影像尺寸 (M, N)。
    輸出: 重建的震幅圖 (numpy array)。
    """
    M, N = image_shape
    center = (N/2, M/2)
    maxRadius = np.hypot(center[0], center[1])
    
    # 1. 對 T 進行逆傅立葉變換，得到對數極座標域中的圖像
    log_polar_recon = np.fft.ifft2(T)
    log_polar_recon = np.abs(log_polar_recon)
    
    # 2. 從對數極座標映射回笛卡兒座標（使用 warpPolar 的逆映射）
    image_recon = cv2.warpPolar(log_polar_recon.astype(np.float32), (N, M), center, maxRadius,
                                flags=cv2.WARP_INVERSE_MAP + cv2.WARP_POLAR_LOG).astype(np.uint8)
    return image_recon
```