Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

マルチモーダル 3 次元医用画像のレジストレーション

この例では、強度ベースのレジストレーションを使用して 2 つのボリューム イメージを自動的に位置合わせする方法を示します。この例では、磁気共鳴法 (MRI) とコンピューター断層撮影 (CT) を使用して取得されたマルチモーダル イメージをレジストレーションします。マルチモーダル イメージは、患者の位置 (並進または回転) やピクセル サイズ (スケーリング) の違いにより位置がずれる可能性があります。

イメージのレジストレーションには "固定イメージ""移動イメージ" があります。固定イメージの位置は変わりません。移動イメージには幾何学的変換が適用され、固定イメージと揃えられます。"強度ベース" のイメージ レジストレーション手法は、ピクセル強度に関する情報を使用して幾何学的変換を計算します。多くの場合、強度ベースのレジストレーションは医用画像および遠隔測定イメージに適しています。特徴ベースのレジストレーションやコントロール ポイント レジストレーションなど、代替のレジストレーション手法の詳細については、イメージ レジストレーション方法の選択を参照してください。この例では、3 次元イメージのレジストレーションに次の 2 つの強度ベースの方法を使用します。

  • imregisterを使用してイメージ レジストレーションを直接行う。

  • imregtformを使用して、移動イメージを固定イメージにマッピングするために必要な幾何学的変換を推定し、imwarpを使用してその変換を適用する。

Medical Imaging Toolbox™ は、このワークフローを簡素化し、患者座標系での空間参照を自動的に管理するオブジェクトと関数を提供します。開始するには、Medical Image Registration (Medical Imaging Toolbox)を参照してください。

イメージの読み込み

この例では、同じ患者から収集した CT イメージと T1 強調 MRI イメージを使用します。この 3 次元 CT と MRI のデータ セットは、Michael Fitzpatrick 博士によって提供された The Retrospective Image Registration Evaluation (RIRE) Dataset からのものです。詳細については、RIRE プロジェクトのホームページを参照してください。

この例では、MRI スキャンが固定イメージで、CT イメージが移動イメージです。イメージは、RIRE プロジェクトで使用されるファイル形式で保存されています。補助関数 helperReadHeaderRIRE を使用して、それぞれのイメージに関連付けられているメタデータを取得します。補助関数は、この例にサポート ファイルとして添付されています。

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

multibandread を使用して、イメージ データを含むバイナリ ファイルを読み取ります。

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");

レジストレーションされていないイメージの表示

imshowpairを使用して各ボリュームの中央のトランスバース スライスを表示することで、レジストレーションされていないボリュームの位置関係を判断します。MRI スライスはマゼンタ、CT スライスは緑です。スケーリングの違いは、イメージのピクセル間隔が異なっていることを示しています。

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

ボリュームを 3 次元オブジェクトとして表示することもできます。複数のボリュームを表示できる viewer3d オブジェクトを作成します。

viewerUnregistered = viewer3d(BackgroundColor="black",BackgroundGradient="off");

関数volshowを使用して、medicalVolume オブジェクトを 3 次元等値面として表示します。ビューアー ウィンドウ内をクリックしてドラッグすることで、ボリュームを回転させることができます。3 次元等値面を固有座標でプロットすると、スライスの厚さがワールド座標での面内のピクセル間隔と異なるため、垂直方向に歪んで表示されます。

volshow(fixedVolume,Parent=viewerUnregistered,RenderingStyle="Isosurface", ...
    Colormap=[1 0 1],Alphamap=1);
volshow(movingVolume,Parent=viewerUnregistered,RenderingStyle="Isosurface", ...
    Colormap=[0 1 0],Alphamap=1);

空間参照の定義

空間参照を使用すると、ピクセル間隔が異なるイメージの表示とレジストレーション結果を改善できます。このデータ セットでは、ヘッダー ファイルによって CT および MRI のイメージのピクセル間隔が定義されています。このメタデータとボクセル単位の各ボリュームのサイズを使用して、空間参照オブジェクトimref3dを作成します。imref3d オブジェクトのプロパティは、ワールド座標系でのイメージ ボリュームの位置と各次元でのピクセル間隔を定義します。たとえば、固定および移動 imref3d オブジェクトの PixelExtentInWorldX プロパティは、固定ボリュームおよび移動ボリュームの X 軸に沿ったピクセル間隔がそれぞれ 1.25 mm および 0.6536 mm であることを示します。空間参照の構築に使用されるヘッダー情報はミリメートルであるため、単位はミリメートルです。

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]

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

関数imregisterはレジストレーションを行い、レジストレーションされた移動イメージ ボリュームを返します。この方法は、レジストレーションされたイメージ データが必要かつレジストレーションを実行するために使用される幾何学的変換が不要な場合に使用します。

関数 imregister は、レジストレーションされるボリュームに追加する塗りつぶしピクセルの値を 0 に設定します。レジストレーション結果の表示を改善するには、塗りつぶしの値がイメージ データ範囲の最小値と等しくなるように、CT 強度を範囲 [0, 1] にスケーリングします。

rescaledMovingVolume = rescale(movingVolume);

関数 imregister では、レジストレーション計算用のオプティマイザーとメトリクス構成を指定する必要があります。関数imregconfigを使用して構成を定義します。CT イメージと MRI イメージは異なるモダリティからのものであるため、モダリティを "multimodal" として指定します。

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

レジストレーション時の収束性能を向上させるために、オプティマイザーの InitialRadius プロパティの値を変更します。

optimizer.InitialRadius = 0.004;

imregister を使用してイメージを位置合わせします。空間参照情報を指定すると、関数がより適切な結果に素早く収束します。2 つの空間参照ボリューム間の位置ずれには並進と回転が含まれるため、剛体変換を使用してイメージをレジストレーションします。

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

レジストレーション結果を評価するには、固定ボリュームとレジストレーションされた移動ボリュームの中央のトランスバース スライスを表示します。イメージがより整列して表示されます。また、関数 imregister によって移動イメージの空間参照が固定イメージと一致するように調整されるため、イメージはスケーリングの違いがなく、各次元で同じ数のボクセルをもつようになります。

figure
imshowpair(movingRegisteredVolume(:,:,centerMoving),fixedVolume(:,:,centerFixed))
title("Registered Transverse Slice - imregister")

volshow を使用して、レジストレーションされたボリュームを 3 次元等値面として表示します。表示をズームしたり回転させたりして、レジストレーション結果を評価できます。

viewerRegistered1 = viewer3d(BackgroundColor="black",BackgroundGradient="off");
volFixed1 = volshow(fixedVolume,Parent=viewerRegistered1,RenderingStyle="Isosurface", ...
    Colormap=[1 0 1],Alphamap=1);
volRegistered1 = volshow(movingRegisteredVolume,Parent=viewerRegistered1,RenderingStyle="Isosurface", ...
    Colormap=[0 1 0],Alphamap=1);

オプションで、volshow によって作成された Volume オブジェクトの Transformation プロパティを設定することにより、レジストレーションされたボリュームを正しいボクセル間隔でワールド座標に表示できます。固定ボリューム ヘッダー ファイルからのボクセル間隔に基づいてボリュームをスケーリングするように、Transformation の値を指定します。レジストレーションされた移動ボリュームは固定ボリュームのボクセル グリッドでリサンプリングされているため、両方のボリュームに固定ボリューム ファイルの間隔を使用します。

sx = fixedHeader.PixelSize(1);
sy= fixedHeader.PixelSize(2);
sz = fixedHeader.SliceThickness;

A = [sx 0 0 0; 0 sy 0 0; 0 0 sz 0; 0 0 0 1];
tformVol = affinetform3d(A);

volFixed1.Transformation = tformVol;
volRegistered1.Transformation = tformVol;

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

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

tform = imregtform(rescaledMovingVolume,Rmoving,fixedVolume,Rfixed, ...
    "rigid",optimizer,metric)
tform = 
  rigidtform3d with properties:

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

tformA プロパティは、移動イメージを固定イメージに位置合わせするために使用される 3 次元アフィン変換行列を指定します。imregtform への入力は空間参照であるため、幾何学的変換によって、移動イメージから固定イメージへとワールド座標系の点がマッピングされます。

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

関数imwarpを使用して、imregtform から移動イメージ ボリュームに幾何学的変換を適用できます。移動ボリューム、移動ボリュームの空間参照、および imregtform からの変換を指定します。名前と値の引数 OutputView は、変換された出力イメージ ボリュームの空間参照を指定します。imregister と同じ結果を生成するには、固定イメージの imref3d オブジェクトとして OutputView を指定します。これにより、固定ボリュームと同じ空間参照をもつレジストレーション イメージ ボリュームが作成されます。

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

レジストレーション結果を評価するには、固定ボリュームと変換された移動ボリュームの中央のトランスバース スライスを表示します。結果は imregister の結果と同じです。

figure 
imshowpair(movingRegisteredVolume(:,:,centerFixed), ...
    fixedVolume(:,:,centerFixed))
title("Registered Transverse Slice - imregtform")

volshow を使用して、レジストレーションされたボリュームを 3 次元等値面として表示します。Transformation プロパティを名前と値の引数として設定して、ボリュームをワールド座標で表示します。

viewerRegistered2 = viewer3d(BackgroundColor="black",BackgroundGradient="off");
volshow(fixedVolume,Parent=viewerRegistered2,RenderingStyle="Isosurface", ...
    Colormap=[1 0 1],Alphamap=1,Transformation=tformVol);
volshow(movingRegisteredVolume,Parent=viewerRegistered2,RenderingStyle="Isosurface", ...
    Colormap=[0 1 0],Alphamap=1,Transformation=tformVol);

参考

| | | | | |

関連する例

詳細