顧險峰教授：解讀醫學影像配準的基本演算法

本文來自微信公眾號「老顧談幾何」，作者顧險峰教授，雷鋒網經顧險峰授權轉載。

顧險峰教授，美國紐約州立大學石溪分校計算機系和應用數學系的終身教授，也是清華大學丘成桐數學科學中心訪問教授。曾獲美國國家自然科學基金CAREER獎，中國國家自然科學基金海外傑出青年獎（與胡事民教授合作），「華人菲爾茨獎」：晨興應用數學金獎。

丘成桐先生和顧險峰博士團隊，將微分幾何，代數拓撲，黎曼面理論，偏微分方程與計算機科學相結合，創立跨領域學科「計算共形幾何」，並廣泛應用於計算機圖形學，計算機視覺，幾何建模，無線感測器網路，醫學圖像等領域。目前已經發表二百篇余篇國際論文，學術專著包括「Computational Conformal Geometry」（計算共形幾何）， 「Ricci Flow for Surface Registration and Shape Analysis」等。

以下是顧險峰教授著作的原文：

紐約長島有著得天獨厚的自然條件，這裡氣候溫潤，海岸線漫長。美麗的石溪大學坐落在長島中部，坐擁一片私人海灘。這片海灘不為人知，空曠寂寥，深沉靜謐。學校師生經常來這裡沐浴海風，遠眺大洋，特別在黃昏時分，在落日餘暉中參省冥思，體悟自然真理。

1970年代初期，年輕的化學系助理教授Paul Lauterbur博士經常帶著他的女兒來這裡散步思考。有一天，三歲的小女孩在岸邊撿到了一隻微小的黑色蛤蜊，Paul為這隻蛤蜊拍攝了人類歷史上第一張活體斷層圖像。很可惜，短視的石溪大學拒絕為Paul的發明申請專利，專家們一致認為「這項發明可能帶來的轉讓費不會彌補專利申請費。」最終，石溪失去了Lauterbur博士，也失去了經濟上騰飛的一次機遇。歷史是公正的，因為發明了核磁共振斷層掃描技術，Paul Lauterbur博士於2003年獲得諾貝爾醫學獎。

數十年來，核磁共振技術和CT斷層掃描技術徹底地革命了醫學，醫學影像技術使得醫生可以直接看到病人體內，從而精準地進行診斷，制定治療方案，檢驗治療效果。很多病變都會誘發器官組織的變形，或者由器官變形所誘發，例如大腦皮層的萎縮退化誘導老年失智，各種腫瘤會在器官表面形成凸起，骨質流失會引起骨骼的變形。因此，醫生可以通過精確比對器官的幾何形狀，來判斷臟器是否反常；通過分析腫瘤的幾何特徵，來判斷腫瘤的良性惡性。這些可以歸結為醫學影像的配准（registration）和分析（analysis）問題。

醫學影像是一個非常龐大的學科，有種類繁多的計算方法。圖像配準的方法也是異常豐富，理論嚴謹，富有實效。從基礎的方法論角度而言，比較普適的主要有基於微分幾何的方法，基於流體力學的方法和基於概率測度理論的方法。雖然，這些途徑中問題提法、數學工具、計算機演算法非常不同，但是最終都可以歸結為幾何偏微分方程，進而轉換為變分優化問題。這三種方法各有優缺點，相輔相成，都有很大的實用價值。每種方法都在迅猛發展，分支眾多。下面，我們簡略討論這三種方法的核心思想和最簡單的演算法。每一種方法深入下去，都是一片浩瀚的海洋。

微分幾何方法

首先，我們從醫學圖像中提取我們感興趣的器官，表示成二維曲面或者三維實體。例如，我們希望研究大腦皮層曲面。首先得到整個顱部的核磁共振斷層掃描圖像；然後對每個斷層的圖像進行分割（segmentation），將不同的組織進行區分，將骨骼、腦灰質、腦白質分開；其次，提取我們感興趣的組織輪廓（contour），主要是腦灰質；然後將每個斷層圖像中的輪廓線摞在一起，組合成封閉曲面。當然，實際演算法遠比理想描述複雜得多，每一步都會引入大量的幾何、拓撲誤差，需要精細而嚴密的方法以提高穩定性和精度。

圖1. 從醫學圖像中重建的大腦皮層曲面（王雅琳）

例如，我們比較同一個人不同時期掃描得到的大腦皮層曲面，來定量分析大腦萎縮情況。從根本層面上而言，我們希望將第一個大腦皮層上的每一個細胞映到第二個大腦皮層上相同的細胞位置。當然，在目前技術現實中，這是無法真正做到的，我們可以通過設計各種能量，加上一些限制條件來逼近理想情形。大腦皮層上面有解剖學特徵，例如溝回皺褶，我們要求微分同胚將主要的溝回映射到相應的溝回上面。同時，我們要求曲面映射的幾何畸變最小。如果我們假設大腦皮層曲面是具有彈性的，那麼我們要求映射帶來的彈性形變勢能最小。彈性勢能一般用調和能量來表述。根據不同的醫學應用，我們可以設計不同的能量。假設我們已經得到同一器官的不同曲面，或者兩個三維實體，註冊問題的微分幾何提法如下：給定三維歐氏空間中的兩個嵌入流形，，給定初始映射 ，求兩個流形間的微分同胚，這個微分同胚和初始微分同胚同倫，，滿足一定的限制條件，例如將特徵點映到相應的特徵點，；同時優化某種能量，例如調和能量。形式上，我們可以把註冊問題表示成：

從理論層面來講，我們需要明確知曉這一問題解的存在性、唯一性、穩定性和正則性。從演算法角度而言，我們需要求出能量的一階變分和二階變分。理論上的困難有很多，通常情況下，我們考慮的泛函空間-流形間的微分同胚空間本身並沒有緊性，即便我們找到了能量極小化序列，其極限不見得是微分同胚，因此解的存在性證明比較困難；如果能量沒有凸性，我們很難保證解的唯一性，能量的凸性很多時候依賴於流形本身的幾何特性，例如高斯曲率；各種直觀的幾何限制條件，例如特徵點的對應，很難轉換成解析形式，往往用偏微分方程來表述，或者用變分的局部線性逼近來顯示錶達。演算法角度的困難也很多，計算機所能處理的計算都是離散的，如何將微積分表述的概念來離散逼近，如何加速收斂，如何用單純復形的組合結構（離散多面體）來描述流形結構，如何尋找適合數值偏微分方程的三角剖分或者樣條表示，這些問題也都具有本質困難。

基於微分幾何的註冊方法也有多種，比較常用的是基於共形幾何（conformal geometry）的內蘊方法，這種方法只依賴曲面的黎曼度量，不需要曲面的嵌入。如圖2所示，其核心方法如下：首先，我們用單值化定理（uniformization）將曲面映射到標準空間中，；其次，我們在標準空間中尋找最優微分同胚，；最後，曲面間的微分同胚由複合映射給出， 。

圖2. 基於共形幾何的曲面配准方法

第一步歸結為尋找一個函數，變換曲面本身的黎曼度量，使得新的黎曼度量誘導常值高斯曲率。根據曲面的拓撲，高斯曲率有三種選擇+1,0，-1，對應的標準空間為單位球面，歐氏平面和雙曲平面，如圖3所示。曲面單值化可以用曲面黎奇流方法得到。這樣，我們將三維空間中曲面間的映射歸結為標準空間的自同胚，減小了搜索範圍，降低了計算難度。

圖3. 曲面單值化

第二步，我們首先將標準空間的自同胚用擬共形映射（quasi-conformal map）來表述。所有的K-擬共形映射構成的空間具有緊性，可以保證解的存在性。每個自映射可以用Beltrami微分來表示，因此我們在Beltrami微分空間中對相應能量進行變分，構造能量極小序列，從而求得極限。這裡問題的關鍵在於如何將限制條件用Beltrami微分來表示，和如何求得歐拉-拉格朗日方程。如果我們要求解是同胚，則Beltrami 微分的模小於1，如果我們要求解保持共形結構，則Beltrami微分和源曲面上的全純二次微分作用為0，等等。很多時候，解的唯一性和正則性依賴於目標曲面上的高斯曲率。例如，曲面間度為1的調和映射，如果目標曲面上的高斯曲率處處為負，則必為微分同胚。再如，曲面間映射的任意同倫類中存在唯一的Teichmuller映射，使得曲面共形結構畸變最小。調和映照和Teichmuller映射在醫學圖像配准問題中被經常使用。在求解這些映射過程中，曲面的拓撲、黎曼度量和共形結構經常被靈活變化，從而達到最優的目的。

圖4. 腰椎骨質流失測量，共形幾何方法（雷諾明）

流場方法

大形變微分同胚度量映射（LDDMM）方法考慮空間中的流場，每個粒子在空中流動，每一剎那的流場用粒子的速度向量場來描述。固定時刻，粒子從起始位置到達終點位置，這給出了空間到自身的微分同胚。

流場的速度場可以被視作是有限支集的光滑映射，所有的流速場記為。給定時變流速場，，給出了空間中每一點每一刻運動的速度向量，我們得到依賴於時間的微分同胚滿足微分方程：

如果我們固定一點，，給出了此點的運動軌跡。我們定義矢量場的模為

歐氏空間的微分同胚變換群記為：

微分同胚群中的右不變距離為：

可以證明，是的一個微分同胚群，在度量下是完備的，任意兩個微分同胚之間存在一條極短測地線。

給定兩幅醫學圖像，記為和，我們尋找一個微分同胚，使得圖像和圖像盡量接近，同時和恆同映射盡量接近。由此，我們得到能量：

,

等價地，我們可以用流速場來表示同樣的能量：

,

圖像註冊問題歸結為極小化能量，求得流速場。

我們可以用梯度下降法來優化能量。微分同胚對路徑非線性依賴，但是的變分可以表示成：

,

這裡。由此，我們可以得到能量關於流速場的梯度，從而可以進行優化。

圖5.孕期嬰兒大腦皮層36周到43周的生長情況（Laurent Risser， LDDMM 方法）

概率測度方法

從根本層面上而言，CT圖像和MRI圖像反映的是人體內部組織的密度。由此，我們可以將醫學圖像的灰度值看成是概率測度。兩個圖像之間的匹配應該滿足如下條件：密度相近的點彼此對應，每個原像點和其像點的距離盡量接近。從這個角度出發，人們提出了基於最優傳輸理論（optimal transportation theory）的圖像匹配方法。

例如我們匹配CT體圖像，設是三維歐氏空間中的凸區域，（例如單位立方體），兩幅圖像代表定義其上的概率測度，，滿足總測度相同，

,

我們尋找自同胚，將推前到，這意味著對於一切波萊爾集,

,

記成，對應的微分方程為雅克比方程：

。

雅克比方程的解有無窮多個。每個映射的傳輸代價定義為：

我們尋找最優傳輸映射，極小化傳輸代價，。當p為2時，最優傳輸映射是某個凸函數的梯度映射，，。代入雅克比方程，

我們的蒙日-安培方程：

.

圖6. 大腦白質圖像匹配（Allen Tannenbaum，最優傳輸方法）

小結

這三類方法基於不同的物理和幾何的理解，所用的理論工具相距甚遠，風格迥異，各有千秋。

基於微分幾何的方法有著強烈的幾何直觀，只用黎曼度量，不需要流形在歐氏空間的嵌入，洗鍊簡潔，嚴密普適，計算效率較高。但是需要比較抽象的黎曼幾何概念，需要先從圖像中提取器官曲面。基於流場的方法直接靈活，適用於各種數據類型，但是需要流形的嵌入，形變過程中曲面不允許自交。同時為了二維曲面匹配，需要計算三維背景空間的自同胚，計算量較大。基於最優傳輸的方法兼顧物理直覺，可向高維直接推廣，但是所解方程高度非線性。另一方面，這幾種方法具有內在的密切聯繫：最優傳輸的蒙日-安培方程和微分幾何中的Alexandrov定理等價；最優傳輸存在流體力學的解釋。總之，所有方法都是在微分同胚群中進行優化，問題具有本質的難度。

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！