イメージの読み込み

この例では、同じ患者について別々の時期に取得した CT イメージと T1 強調 MRI イメージを使用します。この 3 次元 CT と MRI のデータセットは、The Retrospective Image Registration Evaluation (RIRE) Dataset の一部として Michael Fitzpatrick 博士から提供されたものです。

この例では、MRI イメージを固定イメージとして、CT イメージを移動イメージとして指定します。データは、Retrospective Image Registration Evaluation (RIRE) プロジェクトで使用されるファイル形式で保存されています。multibandread を使用して、イメージ データを含むバイナリ ファイルを読み取ります。関数 helperReadHeaderRIRE を使用して、それぞれのイメージに関連付けられているメタデータを取得します。

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

figure
helperVolumeRegistration(fixedVolume,movingVolume);

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

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

空間参照情報を組み込むことによって、表示とレジストレーション結果を改善できます。このデータの場合、CT と MRI のデータセットの解像度はイメージ メタデータで定義されています。そのメタデータを使用して、空間参照オブジェクト imref3d を作成します。

Rfixed  = imref3d(size(fixedVolume),fixedHeader.PixelSize(2), ...
    fixedHeader.PixelSize(1),fixedHeader.SliceThickness)

Rfixed = 
  imref3d with properties:

           XWorldLimits: [0.6250 320.6250]
           YWorldLimits: [0.6250 320.6250]
           ZWorldLimits: [2 106]
              ImageSize: [256 256 26]
    PixelExtentInWorldX: 1.2500
    PixelExtentInWorldY: 1.2500
    PixelExtentInWorldZ: 4
    ImageExtentInWorldX: 320
    ImageExtentInWorldY: 320
    ImageExtentInWorldZ: 104
       XIntrinsicLimits: [0.5000 256.5000]
       YIntrinsicLimits: [0.5000 256.5000]
       ZIntrinsicLimits: [0.5000 26.5000]

Rmoving = imref3d(size(movingVolume),movingHeader.PixelSize(2), ...
    movingHeader.PixelSize(1),movingHeader.SliceThickness)

Rmoving = 
  imref3d with properties:

           XWorldLimits: [0.3268 334.9674]
           YWorldLimits: [0.3268 334.9674]
           ZWorldLimits: [2 114]
              ImageSize: [512 512 28]
    PixelExtentInWorldX: 0.6536
    PixelExtentInWorldY: 0.6536
    PixelExtentInWorldZ: 4
    ImageExtentInWorldX: 334.6406
    ImageExtentInWorldY: 334.6406
    ImageExtentInWorldZ: 112
       XIntrinsicLimits: [0.5000 512.5000]
       YIntrinsicLimits: [0.5000 512.5000]
       ZIntrinsicLimits: [0.5000 28.5000]

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

方法 1: imregister を使用してイメージ レジストレーションを行う

関数 imregister では、レジストレーション結果の質を直接表示して観察することができるレジストレーション済み出力イメージのボリュームを取得できます。

関数 imregconfig を使用して、imregister で使用するための正しいオプティマイザーとメトリクス構成を選択します。今回の 2 つのイメージは、MRI と CT という 2 つの異なるモダリティで撮像したものであるため、"multimodal" オプションが適しています。レジストレーション結果で収束性能を向上させるために、オプティマイザーの InitialRadius プロパティの値を変更します。

[optimizer,metric] = imregconfig("multimodal");
optimizer.InitialRadius = 0.004;

2 つのボリューム間の位置のずれには、平行移動と回転が関わっています。このため、剛体変換を使用してイメージのレジストレーションを行います。空間参照情報を指定すると、imregister で使用されるアルゴリズムがより適切な結果に素早く収束できるようになります。

movingRegisteredVolume = imregister(movingVolume,Rmoving, ...
    fixedVolume,Rfixed,"rigid",optimizer,metric);

レジストレーション済みボリュームの中心を通る軸スライスを表示して、レジストレーションがどの程度適切に実行されているかを把握します。レジストレーション プロセスは、イメージの位置合わせに加えて、移動イメージの空間参照情報を調整することによって、移動イメージが固定イメージと一致するようにします。これで、イメージのサイズが同じになり、適切に位置合わせが行われます。

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

再び helperVolumeRegistration を使用してレジストレーションされたボリュームを表示し、レジストレーションが適切かどうかを引き続き判断します。

helperVolumeRegistration(fixedVolume,movingRegisteredVolume);

方法 2: 3 次元幾何学的変換を推定して適用する

関数 imregister は、イメージのレジストレーションを行いますが、移動イメージに適用された幾何学的変換に関する情報は返しません。推定された幾何学的変換を行う場合は、関数 imregtform を使用することで、変換に関する情報を格納している幾何学的変換オブジェクトを取得できます。imregtform は imregister と同じアルゴリズムを使用し、imregister と同じ入力引数をとります。

tform = imregtform(movingVolume,Rmoving,fixedVolume,Rfixed, ...
    "rigid",optimizer,metric)

tform = 
  rigidtform3d with properties:

    Dimensionality: 3
                 R: [3×3 double]
       Translation: [-15.8648 -17.5692 29.1830]
                 A: [4×4 double]

プロパティ A は、移動イメージを固定イメージに位置合わせするために使用される 3 次元アフィン変換行列を定義します。

tform.A

ans = 4×4

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

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

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

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

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

r = 116

c = 132

p = 13

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

movingRegisteredVolume = imwarp(movingVolume,Rmoving,tform, ...
    "bicubic",OutputView=Rfixed);

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

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