インテル® MKL クックブック

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

目的

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

ソリューション

表記規則:

仮定:

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

  1. p 1 次元 (1D) フーリエ変換 (j = 0 : p - 1 および r = -q : q) を評価します。

  2. g1[j, r] をラジアルグリッド (πr/q)(cosθj , sinθj) からデカルトグリッド (ξ, η) = (-q : q, -q : q) に補間し、f2(πξ/q, πη/q) を取得します。
  3. 逆 2 次元複素数-複素数 FFT を評価して、複素数値の復元イメージ f1 を取得します。

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

ステップ 1 と 3 の計算は、インテル® MKL FFT インターフェイスを呼び出します。ステップ 2 の計算は、インテル® MKL 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 が複素数値で復元される

ソースコード、イメージファイル、メイクファイル: サンプル (http://software.intel.com/en-us/mkl_cookbook_samples) の fft-ct フォルダーを参照してください。

説明

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

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

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

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

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

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

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

親トピック: インテル® MKL クックブック