傅立葉-梅林轉換

傅立葉-梅林轉換（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
```