手順 1: イメージの読み込み

この例では、同じ患者の頭部の 3 次元イメージを 2 つ使用します。レジストレーションの問題では、一方のイメージが固定イメージで他方が移動イメージであると見なします。レジストレーションの目標は、移動イメージを固定イメージに位置合わせすることです。この例では、固定イメージは T1 強調 MRI イメージです。レジストレーションを行う移動イメージは CT イメージです。データは、Retrospective Image Registration Evaluation (RIRE) プロジェクトで使用されるファイル形式で保存されています。multibandread を使用して、イメージ データを含むバイナリ ファイルを読み取ります。関数 helperReadHeaderRIRE を使用して、それぞれのイメージに関連付けられているメタデータを取得します。RIRE ファイル形式の詳細については、次のリンクを使用して確認してください。RIRE データ形式

fixedHeader  = helperReadHeaderRIRE('rirePatient007MRT1.header');
movingHeader = helperReadHeaderRIRE('rirePatient007CT.header');

fixedVolume  = multibandread('rirePatient007MRT1.bin',...
                            [fixedHeader.Rows, fixedHeader.Columns, fixedHeader.Slices],...
                            'int16=>single', 0, 'bsq', 'ieee-be' );
                        
movingVolume = multibandread('rirePatient007CT.bin',...
                            [movingHeader.Rows, movingHeader.Columns, movingHeader.Slices],...
                            'int16=>single', 0, 'bsq', 'ieee-be' );

関数 helperVolumeRegistration は、3 次元レジストレーションの結果の質を判断するために提供されている補助関数です。ビューは対話方式で回転でき、両軸は一致したまま動きます。

helperVolumeRegistration(fixedVolume,movingVolume);

また、imshowpair を使用して固定ボリュームと移動ボリュームの 1 つの平面を表示し、ボリュームの全体的な配置を把握することもできます。imshowpair で重ね合わせたイメージで、灰色の領域は類似した強度の領域に対応し、マゼンタおよび緑色の領域は片方のイメージが他方のイメージよりも明るい部分を示します。imshowpair を使用して、各ボリュームの中心を通る軸スライスに沿ってイメージ ボリュームの位置ずれを観察します。

centerFixed = size(fixedVolume)/2;
centerMoving = size(movingVolume)/2;
figure
imshowpair(movingVolume(:,:,centerMoving(3)), fixedVolume(:,:,centerFixed(3)));
title('Unregistered Axial Slice')

手順 2: 初期レジストレーションの設定

関数 imregconfig を使用すると、imregister での使用に適したオプティマイザーおよびメトリクスの構成を簡単に選択できます。これら 2 つのイメージは、MRI と CT という 2 つの異なる様式で撮画したものであるため、'multimodal' オプションが適しています。

[optimizer,metric] = imregconfig('multimodal');

imregister で使用されるアルゴリズムは、入力イメージの解像度や位置に関する空間参照情報が指定されると、より適切な結果に素早く収束します。この例では、CT と MRI のデータセットの解像度がイメージ メタデータで定義されています。そのメタデータを使用して、レジストレーションの入力引数として渡す空間参照オブジェクト imref3d を作成します。

Rfixed  = imref3d(size(fixedVolume),fixedHeader.PixelSize(2),fixedHeader.PixelSize(1),fixedHeader.SliceThickness);
Rmoving = imref3d(size(movingVolume),movingHeader.PixelSize(2),movingHeader.PixelSize(1),movingHeader.SliceThickness);

空間参照オブジェクトのプロパティは、関連付けられているイメージ ボリュームのワールド座標系における位置および各次元のピクセル範囲を定義します。Rmoving の XWorldLimits プロパティは、移動ボリュームの X 次元方向の位置を定義します。PixelExtentInWorld プロパティは、X 次元における (列に沿った) 各ピクセルのサイズをワールド単位で定義します。移動ボリュームはワールド座標系の X 軸の 0.3269 から 334.97 に広がっており、各ピクセルの範囲は 0.6536 mm です。単位がミリメートルなのは、空間参照の作成に使用されたヘッダー情報がミリメートル単位であったためです。ImageExtentInWorldX プロパティによって、移動イメージ ボリュームの全範囲がワールド単位で決められます。

Rmoving.XWorldLimits

ans = 1×2

    0.3268  334.9674

Rmoving.PixelExtentInWorldX

ans = 0.6536

Rmoving.ImageExtentInWorldX

ans = 334.6406

2 つのボリューム間の位置のずれには、平行移動、拡大縮小、回転が関わっています。イメージのレジストレーションには相似変換を使用します。

まず、レジストレーション結果の質を確認するため、imregister を使用して、直接表示と観察のできるレジストレーション済み出力イメージのボリュームを取得します。

レジストレーション結果で収束の結果を向上させるために、オプティマイザーの InitialRadius プロパティに既定以外の設定を指定します。

optimizer.InitialRadius = 0.004;
movingRegisteredVolume = imregister(movingVolume,Rmoving, fixedVolume,Rfixed, 'rigid', optimizer, metric);

再び imshowpair を使用して、レジストレーション済みボリュームの中心を通る軸スライスの配置を調べるプロセスを繰り返し、レジストレーションがどの程度適切に実行されているかを把握します。

figure
imshowpair(movingRegisteredVolume(:,:,centerFixed(3)), fixedVolume(:,:,centerFixed(3)));
title('Axial Slice of Registered Volume')

上記の軸スライスからは、レジストレーションが適切に実行されているように見えます。再び helperVolumeRegistration を使用してレジストレーションされたボリュームを表示し、レジストレーションが適切かどうかを引き続き判断します。

helperVolumeRegistration(fixedVolume,movingRegisteredVolume);

手順 3: 移動イメージを固定イメージに位置合わせする 3 次元幾何学的変換の取得

関数 imregtform は、幾何学的変換の推定を行う場合に使用できます。この推定は、レジストレーションされた出力イメージを形成するために imregister で使用されます。imregtform は imregister と同じアルゴリズムを使用し、imregister と同じ入力引数をとります。imregister の結果として出力されたボリュームを視覚的に調べると、レジストレーションが適切に実行されたことを示しているため、このレジストレーション結果に関連付けられている幾何学的変換を取得するために同じ入力引数を用いて imregtform を呼び出すことができます。

geomtform = imregtform(movingVolume,Rmoving, fixedVolume,Rfixed, 'rigid', optimizer, metric)

geomtform = 
  affine3d with properties:

                 T: [4x4 double]
    Dimensionality: 3

imregtform を実行すると、幾何学的変換オブジェクトが得られます。このオブジェクトには、3 次元アフィン変換行列を定義するプロパティ T が含まれています。

geomtform.T

ans = 4×4

    0.9704   -0.0143   -0.2410         0
    0.0228    0.9992    0.0324         0
    0.2404   -0.0369    0.9700         0
  -15.8648  -17.5692   29.1830    1.0000

幾何学的変換の transformPointsForward メソッドを使用して、レジストレーションの結果として移動イメージの点 [u,v,w] がマッピングされる位置を判定することができます。imregtform に対し空間参照入力が指定されたので、幾何学的変換によって、移動イメージから固定イメージへとワールド座標系の点がマッピングされます。以下では、transformPointsForward メソッドを使用して、ワールド座標系における移動イメージの変換後の中心位置を判定しています。

centerXWorld = mean(Rmoving.XWorldLimits);
centerYWorld = mean(Rmoving.YWorldLimits);
centerZWorld = mean(Rmoving.ZWorldLimits);
[xWorld,yWorld,zWorld] = transformPointsForward(geomtform,centerXWorld,centerYWorld,centerZWorld);

Rfixed の worldToSubscript メソッドを使用して、移動ボリュームの中心と位置合わせする固定ボリュームの要素を決定することができます。

[r,c,p] = worldToSubscript(Rfixed,xWorld,yWorld,zWorld)

r = 116

c = 132

p = 13

手順 4: 移動イメージ ボリュームへの幾何学的変換の推定の適用

関数 imwarp を使用して、imregtform で推定された幾何学的変換を 3 次元ボリュームに適用することができます。'OutputView' の名前/値を使用して、リサンプリングされた出力イメージのワールド範囲と解像度を決める空間参照引数を定義します。固定イメージに関連付けられている空間参照オブジェクトを使用することで、imregister と同じ結果を得ることができます。これによって、ワールド範囲と解像度が固定イメージと移動イメージで等しい出力ボリュームが作成されます。両方のボリュームでワールド範囲と解像度が同じになれば、移動ボリュームと固定ボリュームの各サンプル間でピクセル同士が対応するようになります。

movingRegisteredVolume = imwarp(movingVolume,Rmoving,geomtform,'bicubic','OutputView',Rfixed);

再び imshowpair を使用して、imwarp で作成されたレジストレーション済みボリュームの中心を通る軸スライスを表示します。

figure 
imshowpair(movingRegisteredVolume(:,:,centerFixed(3)), fixedVolume(:,:,centerFixed(3)));
title('Axial Slice of Registered Volume')