高速フーリエ変換を使用したコンピューター・トモグラフィー・イメージの復元

目的

高速フーリエ変換 (FFT) 関数を使用してコンピューター・トモグラフィー (CT) データから元のイメージを復元する。

解決方法

表記規則:

インデックス範囲の表記には MATLAB* で使われている表記規則を採用しています。

例えば、k=-q : q は、k=-q, -q+1, -q+2,…, q-1, q を意味します。
f(x) は関数 f のポイント x の値を意味し、f[n] は離散データセット f の n 番目の要素の値を意味します。

仮定:

2 次元 (2D) イメージの密度 f(x, y) は単位円の外部で消えます。

x² + y² > 1 の場合 f = 0
CT データは、角度 θ_j = jπ/p で取得したイメージの p の投影からなります。ここで、j = 0 : p - 1 です。
各投影には、次の線に沿ったイメージの積分に近似する 2q + 1 密度値 g[j, l] = g(θ_j , s_l) が含まれます。

(x, y) = (-t sinθ_j+ s_l cosθ_j , t cosθ_j+ s_l sinθ_j),

ここで、l = -q : q、s_l= l/q、t は積分パラメーターです。

離散イメージ復元アルゴリズムは次のステップからなります。

p 1 次元 (1D) フーリエ変換 (j = 0 : p - 1 および r = -q : q) を評価します。
g₁[j, r] をラジアルグリッド (πr/q)(cosθ_j , sinθ_j) からデカルトグリッド (ξ, η) = (-q : q, -q : q) に補間し、f₂(πξ/q, πη/q) を取得します。
逆 2 次元複素数-複素数 FFT を評価して、複素数値の復元イメージ f₁ を取得します。

ここで、f(m/q, n/q) ≈f₁[m, n] (m = -q : q および n = -q : q) です。

ステップ 1 と 3 の計算は、oneMKL FFT インターフェイスを呼び出します。ステップ 2 の計算は、oneMKL FFT で使用されるデータレイアウトに合わせた、単純なバージョンの補間を実装します。

元の CT イメージの復元 (C/C++)

// 宣言
int Nq = 2*(q+1); // インプレースr2c FFT 用の空間
void       *gmem = mkl_malloc( sizeof(float)*p*Nq, 64 ); 
float       *g = gmem;   // g[j*Nq + ell+q] 
complex *g1 = gmem; // g1[j*Nq/2 + r+q]
 
// g をCT データで初期化
for (int j = 0; j < p; ++j) 
    for (int ell = 0; ell < 2*q+1; ++ell) { 
        g[j*Nq + ell+q] = get_g(theta_j,s_ell); 
    } 

// ステップ1: 1D 実数-複素数FFT を構成して計算 
DFTI_DESCRIPTOR_HANDLE h1 = NULL; 
DftiCreateDescriptor(&h1,DFTI_SINGLE,DFTI_REAL,1,(MKL_LONG)2*q); 
DftiSetValue(h1,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX); 
DftiSetValue(h1,DFTI_NUMBER_OF_TRANSFORMS,(MKL_LONG)p); 
DftiSetValue(h1,DFTI_INPUT_DISTANCE,(MKL_LONG)Nq); 
DftiSetValue(h1,DFTI_OUTPUT_DISTANCE,(MKL_LONG)Nq/2); 
DftiSetValue(h1,DFTI_FORWARD_SCALE,fscale); 
DftiCommitDescriptor(h1); DftiComputeForward(h1,g); // gmem にg1 が含まれる 

// ステップ2: g1 をf2 に補間- ここでは省略 
complex *f = mkl_malloc( sizeof(complex) * 2*q * 2*q, 64 ); 

// ステップ3: 2D 複素数-複素数FFT を構成して計算 
DFTI_DESCRIPTOR_HANDLE h3 = NULL; 
MKL_LONG sizes[2] = {2*q, 2*q}; DftiCreateDescriptor(&h3,DFTI_SINGLE,DFTI_COMPLEX,2,sizes); 
DftiCommitDescriptor(h3); 
DftiComputeBackward(h3,f); // f が複素数値で復元される

ソースコード、イメージファイル、メイクファイル: サンプル (https://www.intel.com/content/dam/develop/external/us/en/documents/mkl-cookbook-samples-120115.zip (英語)) の fft-ct フォルダーを参照してください。GitHub の block_lu_decomposition (英語) も参照してください。

説明

このコードは、最初に DftiComputeForward 関数の単一の呼び出しで 1 次元フーリエ変換のバッチを計算する oneMKL FFT ディスクリプターを構成した後、バッチ変換を計算します。複数の変換の距離は、対応する領域 (入力は実数、出力は複素数) の要素で設定されます。デフォルトでは、変換はインプレースです。

メモリーのフットプリントを小さくするため、FFT は「インプレース」で計算されます。つまり、計算の結果は入力データに上書きされます。インプレースの実数-複素数 FFT では、FFT の結果を格納するのに必要なメモリーが入力よりもやや多くなるため、入力配列は余分な空間を予約します。

ステップ 1 の入力で、配列 g には p x (2q+1)実数値のデータ要素 g(θ_j, s_l) が含まれます。このステップの出力の同じメモリーには、p x (q + 1) 複素数値の出力要素 g₁(θ_j, πr/q) が含まれます。結果の複素共役部分は格納されません。そのため、配列 g1 は、r の q + 1 値のみを指します。

g₁ から f₂ に補間するため、ステップ 3 の逆 FFT の複素数値のデータ f₂(ξ, η) と複素数値の出力 f₁(x, y) を格納する追加の配列 f が割り当てられます。補間ステップは oneMKL 関数を呼び出しません。C++ 実装は、この手法のソースコード (main.cpp) の step2_interpolation 関数を参照してください。補間の最も単純な実装は次のとおりです。

単位円の内部のすべての (ξ, η) で、最も近い (θ_j , πr/q) を見つけて、g₁(θ_j , πr/q) の値を f₂ に使用します。
単位円の外部のすべての (ξ, η) で、f₂ を 0 に設定します。
区間 -π < θ_j < 0 に対応する (ξ, η) の場合は、実数-複素数変換の結果の共役偶数プロパティー g₁(θ, ω)=conj(g(-θ, -ω)) を使用します。

ステップ 1 の FFT は、計算された出力が e^i(πr/q)q= (-1)^r で乗算される、表現区間の半分でオフセットされたデータに適用されることに注意してください。個別のパスで修正する代わりに、補間は乗数を考慮しています。

同様に、ステップ 3 の 2D FFT はイメージの中心を隅にシフトする出力を生成し、ステップ 2 はステップ 3 に入力をシフトする過程でこれを防いでいます。

ステップ 3 は、配列 f に含まれる補間されたデータについて 2 次元 (2q) x (2q) 複素数-複素数 FFT を計算します。この計算の後、複素数値のイメージ f₁ からビジュアルピクチャーへの変換が行われます。CT イメージ復元を実装する完全な C++ プログラムは、この手法のソースコード (main.cpp) を参照してください。