Enter your search text in the box above

Select a result to preview

MRI 背景與 Iterative Reconstruction 問題建構

#application-case-study #mri-reconstruction #k-space #conjugate-gradient #algorithm-selection

重點總覽 (Overview)

項目 內容 重點
MRI 兩階段 acquisition (scan) → reconstruction 掃描在 k-space 取樣;重建把樣本變成影像 m(r)
k-space spatial frequency / Fourier transform domain 掃描沿著預定 trajectory 在頻域取樣
三方衝突 scan time ↔ resolution ↔ SNR 改善其一通常犧牲其他;平行運算同時突破三者
Cartesian 軌跡 均勻網格 → W(k) 為常數 直接 inverse FFT,極快
Non-Cartesian spiral / radial / rosette 對動作不敏感、降低硬體需求 → 需 gridding 或 iterative solver
逆問題公式 ρ = (FᴴF + λWᴴW)⁻¹ FᴴD F 模型化物理、D 為掃描資料、W 為先驗(解剖)約束
為何 iterative F 太大(128³≈2M 欄)→ 直接 Gaussian elimination 不可行 改用 Conjugate Gradient (CG) 迭代解
平行化目標 FHD(= FᴴD)與 Q(加速 FᴴF) 本章聚焦 FHD;CG solver 暫不平行化
Important

本筆記只建立「問題」:MRI 物理背景 → k-space 軌跡 → 為何要 iterative linear solver → FHD/Q 是計算瓶頸。實際 kernel 平行化、記憶體優化、硬體三角函數 在 17-Iterative-MRI-Reconstruction/02-FHD-Kernel-Parallelism-Structure 等姊妹筆記。

MRI 兩階段與 k-space (MRI Acquisition / Reconstruction & k-space)

MRI 安全、非侵入地探測生物組織結構與功能,流程分兩階段:

   Acquisition (scan)                      Reconstruction
   ┌─────────────────────┐                 ┌──────────────────────┐
   │ scanner 沿 trajectory │   k-space       │  把樣本轉換成影像      │     image
   │ 在 k-space 取樣        │ ── 資料 s(k) ──▶ │  (估計組織形狀/紋理)   │ ───▶  m(r)
   └─────────────────────┘                 └──────────────────────┘
        (空間頻率域 / Fourier transform 域)
三方指標互相衝突

短掃描時間、高解析度、高 SNR(signal-to-noise ratio) 通常無法兼得 — 改善一項往往犧牲另一或兩項。大規模平行運算正是讓三者同時提升的突破口。

重建方程式 (Reconstruction equation, Eq. 17.1)

m(r)=∑jW(k_j)s(k_j)ei2πk_j⋅r
符號 意義
m(r) 重建後的影像(voxel 值)
s(k) 量測到的 k-space 資料
W(k) weighting function:補償非均勻取樣密度;亦可當 apodization filter 降噪/降 artifacts

Cartesian vs Non-Cartesian Trajectory

  (A) Cartesian 均勻網格           (C) Spiral 非均勻 (non-Cartesian)
   ky                              ky
   │  • • • • •                     │       .-‾‾‾-.
   │  • • • • •                     │     /   ___  \
   │  • • • • •                     │    |  /    \  |
   │  • • • • •                     │     \  \__/  /
   │  • • • • •  ── kx              │       '----'   ── kx
   指數項均勻間隔 → inverse FFT      指數項非均勻 → gridding 或 linear solver
面向 Cartesian Non-Cartesian(spiral / radial / rosette)
取樣點 均勻 Cartesian 網格 螺旋、放射線(projection imaging)、玫瑰線
W(k) 常數,可移出加總 隨密度變化,需保留
指數項間隔 均勻 → 加總即 FFT 非均勻,不再是 FFT
重建法 inverse FFT(極快) gridding+FFT,或 iterative linear solver
對患者動作 較敏感 較不敏感(motion artifacts 少)
掃描硬體需求 較高 較低;可自校正 field inhomogeneity
典型用途 臨床標準 高 SNR 應用,如 sodium imaging(追蹤中風/癌症組織健康)
Cartesian 為何能用 FFT

均勻網格下 W(k) 是常數可提出加總,且指數項在 k-space 均勻間隔 → 整個加總恰好是 inverse FFT,計算極有效率,因此臨床廣泛採用。

Non-Cartesian 的 gridding 折衷

一個常見做法是 gridding:先把非均勻樣本內插到均勻 Cartesian 網格,再用 FFT(見 Ch.7 的 convolution gridding mask)。但簡單 gridding/iFFT 會有嚴重 artifacts(PSNR 較低);本章的 iterative 法精度更高。FHD 也可用 gridding 近似在數秒內完成,但會降低最終影像品質。

為何要 Iterative Linear-Solver Reconstruction (Why Iterative Reconstruction)

Haldar 與 Liang 提出的 linear-solver-based iterative reconstruction(Stone et al., 2008)可顯式建模掃描器資料擷取的物理,進而降低 artifacts;代價是計算昂貴。

準貝氏估計問題 (quasi-Bayesian formulation, Fig. 17.3)

ρ=(FHF+λWHW)−1FHD
符號 意義
ρ 重建影像的 voxel 向量(要解的未知數)
F 模型化成像物理的矩陣
D 掃描器資料樣本向量
W 納入先驗資訊(如解剖約束)的矩陣;由高解析水分子掃描導出
Fᴴ, Wᴴ Hermitian / conjugate transpose(轉置後再取每項複共軛,a+ib → a−ib)
λ 正則化權重

表面上只是:矩陣相乘相加 (FᴴF+λWᴴW)、矩陣–向量乘 FᴴD、矩陣求逆、再相乘。問題出在矩陣尺寸。

為何不能直接求逆

方法 直接解 (Gaussian elimination) 迭代解 (Conjugate Gradient, CG)
矩陣求逆 顯式計算 (·)⁻¹ 迭代逼近,不顯式求逆
對巨大 F 實務上不可行(intractable) 偏好做法
每次迭代主成本 — matrix-vector mult (FᴴF+λWᴴW)·ρ
矩陣有多大

即使只是中等解析度的 128³ voxel 重建,F 就有 128³ = 2,097,152 欄,每欄有 N 個元素(N = 使用的 k-space 樣本數 = D 的大小)。F 極其龐大 → 直接 Gaussian elimination 不可行,必須改用 CG 迭代法。這類巨大維度在 big data analytics 用 iterative solver 估計大量雜訊觀測資料時很常見。

CG 演算法:每次迭代更新目前影像估計 ρ,以改進準貝氏 cost function;其效率主要由每次迭代的 (FᴴF+λWᴴW)·ρ matrix-vector 乘法決定。

FHF / FHD / Q 計算目標 (The Computation Target)

兩個關鍵結構使 CG 可行:W 通常稀疏 → WᴴW 可高效實作;FᴴF 是 Toeplitz → 可用 FFT 做高效 matrix-vector 乘。Stone et al. 用 GPU 算出 Q 這個資料結構,讓我們不必真的算出 FᴴF 就能快速做含 FᴴF 的 matrix-vector 乘。

   scanner 資料 D ─┐
                   ▼
   FHD = Fᴴ·D    matrix-VECTOR  (~3h CPU @128³, 每次掃描都要)  ◀── 本章平行化目標
                   │
                   ▼
   find ρ : CG 迭代  ── 每步用 Q 快速算 (FᴴF)·ρ
                   │     Q : matrix-MATRIX (數天 CPU,每軌跡只算一次,可重用)
                   ▼
   ρ = 重建影像 (voxels)
FHD (= FᴴD) Q (加速 FᴴF)
運算型態 matrix–vector 乘 matrix–matrix 乘(計算量大得多)
核心計算結構 與 Q 相同 與 FHD 相同
CPU 時間 ~3 小時(128³) 數天
計算頻率 每次影像擷取都要(D 每次不同) 每個 scanner+trajectory 只算一次,可重用
為何只討論 FHD

FHD 與 Q 的核心計算結構完全相同,Q 只是量更大(matrix-matrix)。因此從平行化角度討論其中一個就夠;選 FHD 因為它每次掃描都要跑。

為何把 CG solver 排除在平行化之外

因為 Q 已預先算好,"find ρ" 的 CG 步驟 在序列 CPU 上只占整體重建 < 1% 的時間。所以本章把 CG solver 留在平行化範圍外,專注於 FHD。
但要注意 Amdahl 式反轉:FHD 大幅加速(全章約 10× speedup)後,原本 < 1% 的 CG solver 反而可能變成約 50% 的執行時間 — 加速成功的階段會讓其他階段成為新瓶頸(詳見 17-Iterative-MRI-Reconstruction/04-Hardware-Trig-Accuracy-and-Performance-Tuning)。

考試/面試重點 (Exam / Test Patterns)

情境 / 關鍵字 答案 / 技巧
「k-space 是什麼」 空間頻率域(spatial frequency / Fourier transform 域);掃描器沿 trajectory 在此取樣
「MRI 哪兩階段」 acquisition (scan) + reconstruction
「Cartesian 軌跡用哪種重建」 inverse FFT(W 為常數、指數項均勻間隔)
「為何用 non-Cartesian」 對動作較不敏感、自校正 field inhomogeneity、降低硬體需求;支援高 SNR(如 sodium imaging)
「non-Cartesian 為何不能直接 FFT」 指數項不再均勻間隔,加總不是 FFT 形式 → 需 gridding 或 iterative solver
「W(k) 角色」 補償非均勻取樣密度 + apodization(降噪/降 artifacts)
「為何不直接矩陣求逆」 F 太大(128³≈2M 欄)→ Gaussian elimination intractable → 用 CG
「CG 每次迭代主成本」 matrix-vector 乘 (FᴴF+λWᴴW)·ρ
「Q 是什麼」 預算的資料結構,讓含 FᴴF(Toeplitz)的 matrix-vector 乘無需真的形成 FᴴF;每軌跡算一次
「FHD vs Q 差異」 核心結構相同;Q 是 matrix-matrix(每軌跡一次),FHD 是 matrix-vector(每次掃描)
「為何專注 FHD 而非 CG」 CG 占 CPU 時間 < 1%(Q 已預算);FHD ~3h 且每次擷取都要
「Hermitian transpose 是什麼」 轉置後再取每項複共軛(a+ib → a−ib)