HWC
← Note

No date

Free-Form Deformation Using Lower-Order B-spline for Nonrigid Image Registration

PaperImage Registration

快速解釋

這篇論文想解決一個很實際的問題：在非剛性影像配準裡，傳統的 B-spline free-form deformation（FFD）通常使用三次 B-spline，因為它夠平滑，但每次更新形變參數時都要處理較大的支撐區域，因此計算成本高。作者的做法不是直接放棄高階 spline 的平滑性，而是把優化階段改成較低階的 B-spline，並在每次隨機梯度更新時加入隨機擾動。這樣一來，每一步的計算更便宜，但在期望上仍可逼近高階 B-spline 的效果。等參數學完之後，再把同一組控制點係數放回高階 B-spline 中，得到平滑的最終形變場。作者把這件事稱為 postponed smoothing。

問題設定

在非剛性配準中，給定固定影像 F 與移動影像 M，目標是找到一個空間變換，使得變換後的移動影像與固定影像在某個相似度準則下最匹配。作者使用的整體最佳化問題可以寫成：

μ^=arg⁡min⁡μC(F,M(Tμ(x)))\hat{\mu} = \arg\min_{\mu} C\bigl(F, M(T_{\mu}(x))\bigr)

其中，mumu 是 FFD 的控制點參數，CC 是影像相似度或成本函數。

傳統 B-spline FFD 會在控制點格點上定義形變係數，再透過 tensor-product B-spline basis 對整個空間做插值。問題在於：B-spline 的階數越高，單一像素或體素會受到更多控制點影響，因此每次計算變形、梯度與更新時都更昂貴。在 D 維空間中，如果把 n 階 B-spline 換成 m 階 B-spline，且 mlenm le n，則相對計算成本大約變成：

(m+1n+1)D\left(\frac{m+1}{n+1}\right)^D

以 3D 為例，三次 B-spline 會讓一個體素同時受 43=644^3 = 64 個控制點影響；二次是 33=273^3 = 27 個；一次是 23=82^3 = 8 個；零次則只有 11 個。這正是高階平滑性與低階效率之間的核心張力。

核心想法

作者的方法可以拆成三個部分。

第一，優化時不用直接採用高階 B-spline，而是先用較低階的 basis 來表示形變。這會立刻縮小每個位置需要參與計算的控制點數量，因此每次隨機梯度更新都更快。

第二，作者指出高階 B-spline 可以看成低階 B-spline 與另一個核函數的卷積。對一維 B-spline 而言：

βn(x)=β0∗β0∗⋯∗β0⏟n+1 個 β0\beta_n(x) = \underbrace{\beta_0 * \beta_0 * \cdots * \beta_0}_{n+1\text{ 個 }\beta_0}

更關鍵的是，當 m≤nm \le n 時，可以把高階基底重寫成：

βn(x)=(βm∗βn−m−1)(x)\beta_n(x) = \bigl(\beta_m * \beta_{n-m-1}\bigr)(x)

也就是：

βn(x)=∫βm(x−t) βn−m−1(t) dt\beta_n(x) = \int \beta_m(x-t)\,\beta_{n-m-1}(t)\,dt

如果把

p(t)=βn−m−1(t)p(t) = \beta_{n-m-1}(t)

視為一個機率密度函數，那麼高階基底就等於「隨機位移後的低階基底」在期望下的平均效果。這是整篇論文最重要的觀察。

第三，作者把這個觀察搬到 FFD 優化中。標準的 n 階 FFD 可寫成：

Tμn(x)=x+∑ξ∈Ξcξ ΦDn ⁣(xη−ξ)T^n_{\mu}(x) = x + \sum_{\xi \in \Xi} c_{\xi}\,\Phi_D^n\!\left(\frac{x}{\eta} - \xi\right)

其中，XiXi 是控制點索引集合，cxic_{xi} 是控制點位移係數，etaeta 是控制點間距，PhiDnPhi_D^n 是 D 維的 tensor-product B-spline basis。

作者改成在每次更新時，先抽一個隨機擾動 t，再用較低階的 basis 形成 perturbed transform：

Tμm(x,t)=x+∑ξ∈Ξcξ ΦDm ⁣(xη−ξ−t)T^m_{\mu}(x,t) = x + \sum_{\xi \in \Xi} c_{\xi}\,\Phi_D^m\!\left(\frac{x}{\eta} - \xi - t\right)

因此每一步優化不再直接使用高階 basis，而是使用帶有隨機位移的低階 basis。若擾動分佈選得對，這個低階模型在期望上對應到高階模型。

最後，當參數收斂後，作者不直接把低階形變當成最終答案，而是把學到的同一組係數放回高階 spline 中，得到較平滑的最終變換：

Tμ^n(x)T^n_{\hat{\mu}}(x)

這就是 postponed smoothing。它的意義是：把高階 spline 的平滑性延後到輸出階段，而不是在整個最佳化過程中都承擔高階 basis 的成本。

關鍵公式

論文中的最佳化更新採用隨機梯度下降類方法，基本形式為：

μk+1=μk−ak g~(μk)\mu_{k+1} = \mu_k - a_k\,\tilde{g}(\mu_k)

其中，aka_k 是步長，tildeg(muk)tilde{g}(mu_k) 是以隨機樣本估計出的梯度。

在加入隨機擾動之後，更新式變成：

μk+1=μk−ak g~(μk,tk)\mu_{k+1} = \mu_k - a_k\,\tilde{g}(\mu_k, t_k)

其中，tkt_k 是第 k 次迭代抽到的擾動。這表示每次梯度估計其實都是在一個隨機平移過的低階 B-spline 網格上進行。

整篇方法的數學核心可以濃縮成下面這條鏈：

βn(x)=(βm∗βn−m−1)(x)=∫βm(x−t) βn−m−1(t) dt=∫βm(x−t) p(t) dt\beta_n(x) = \bigl(\beta_m * \beta_{n-m-1}\bigr)(x) = \int \beta_m(x-t)\,\beta_{n-m-1}(t)\,dt = \int \beta_m(x-t)\,p(t)\,dt

這條式子說明了高階基底可以由低階基底加上隨機擾動來表示。也就是說，作者不是單純以低階 spline 近似高階 spline，而是利用卷積結構，把「高階平滑」轉換成「低階基底在隨機位移下的期望」。

論文中也比較了不同的擾動分佈與階數組合，例如：

  1. Dirac-3-3：標準三次 B-spline FFD。
  2. Dirac-1-3：用一次 B-spline 做優化，但不加擾動，最後再做三次平滑。
  3. Uniform-1-3：一次 B-spline 加上均勻分佈擾動，再做三次平滑。
  4. Gaussian-1-3：一次 B-spline 加上高斯分佈擾動，再做三次平滑。

從這些設定可以看出，作者要驗證的不是「低階 spline 是否能完全取代高階 spline」，而是「低階 spline 加擾動與延後平滑，能否在更低成本下保留甚至改善高階方法的表現」。

模型結構

這篇論文不是深度學習模型，而是一個形變參數化與最佳化流程。整體結構可以整理成以下步驟：

  1. 在影像空間上建立控制點格點，並選定低階 B-spline 階數 m、最終輸出的高階階數 n，以及擾動分佈 p(t)p(t)。
  2. 用低階 B-spline basis 定義優化時使用的 FFD，讓每次成本函數與梯度估計只需依賴較少的控制點。
  3. 在每次隨機梯度更新前，先抽一個隨機擾動 tkt_k，以 perturbed transform 計算目前的配準誤差與梯度。
  4. 透過反覆的隨機更新估計控制點係數 hatmuhat{mu}。
  5. 最後不直接輸出低階變換，而是把 μ^\hat{\mu} 放入高階 B-spline FFD，得到平滑的最終形變場 Thatmun(x)T^n_{hat{mu}}(x)。

如果用一句話總結這個結構，就是：優化階段用便宜的低階 spline 做搜尋，輸出階段再用高階 spline 做重建。作者在肺部與腦部資料上的實驗顯示，這種設計通常能顯著降低運算時間，同時維持甚至提升配準精度與形變場的平滑性。