Main Content

投影データからイメージを復元

ライブ スクリプトを開く

この例では、頭部ファントム イメージからパラレル ビーム投影およびファン ビーム投影を行う方法と、ラドン変換およびファンビーム変換を使用してイメージを復元する方法を示します。

関数 radon と関数 iradon が投影にパラレル ビーム形状を使用するのに対し、fanbeam と ifanbeam はファンビーム形状を使用します。下記の例は、パラレル ビームとファンビームの形状を比較するために、各形状に対する合成投影を作成します。次に、これらの合成投影を使用して元のイメージを復元します。

イメージの復元を必要とする現実世界での用途は、X 線による吸収断層の撮影です。物理的な標本を通過する放射線の減衰を、さまざまな角度から測定することで投影が形成されます。元のイメージは標本の断面であり、ここでの強度は標本の密度を表すと考えることができます。映像は特殊な医療撮像装置で収集され、標本の内部イメージは iradon または ifanbeam を使用して復元されます。

関数 iradon は、パラレル ビーム投影によりイメージを復元します。パラレル ビーム形状では、各投影は、指定した角度でイメージを通った一連の線積分を組み合わせて形成されます。関数 ifanbeam は、ファンビーム投影からイメージを復元します。ファンビーム投影には 1 つのエミッターと複数のセンサーがあります。

頭部ファントムの作成

テスト イメージは Shepp-Logan の頭部ファントムであり、関数 phantom を使用して生成できます。ファントム イメージは、現実世界での人体頭部実際の断層を撮影したイメージで見られるさまざまな要素を示します。外側の明るい楕円の外殻は頭蓋骨と考えられ、内部の多くの楕円は脳の特徴または腫瘍と考えられます。 

P = phantom(256);
imshow(P)

パラレル ビーム - 合成投影の計算

パラレル ビーム形状を使用して合成投影を計算し、投影角度をさまざまに変化させます。radon の各呼び出しでは、配列が出力され、各列はいずれかの角度で対応する theta のラドン変換になります。 

theta1 = 0:10:170; 
[R1,~] = radon(P,theta1); 
num_angles_R1 = size(R1,2)
num_angles_R1 = 18

theta2 = 0:5:175;  
[R2,~] = radon(P,theta2);
num_angles_R2 = size(R2,2)
num_angles_R2 = 36

theta3 = 0:2:178;  
[R3,xp] = radon(P,theta3); 
num_angles_R3 = size(R3,2)
num_angles_R3 = 90

各角度で、投影は xp 軸沿いの N ポイントで計算されます。N は、可能なすべての投影角度で各ピクセルが投影されるようなイメージの対角距離に基づく定数です。

N_R1 = size(R1,1)
N_R1 = 367

N_R2 = size(R2,1)
N_R2 = 367

N_R3 = size(R3,1)
N_R3 = 367

そのため、小さな頭部ファントムを使用する場合は、xp 軸のポイントを少なくして計算する必要があります。

P_128 = phantom(128);
[R_128,xp_128] = radon(P_128,theta1);
N_128 = size(R_128,1)
N_128 = 185

投影データ R3 を表示します。元のファントム イメージの一部の特徴が R3 のイメージで見られます。R3 の最初の列は 0 度での投影に対応し、垂直方向に積分します。真ん中の列は 90 度での投影に対応し、水平方向に積分します。90 度での投影は、ファントムの最も外側の楕円の垂直半軸が長いため、0 度での投影に比べてプロファイルの幅が広くなります。

imagesc(theta3,xp,R3)
colormap(hot)
colorbar
xlabel('Parallel Rotation Angle - \theta (degrees)'); 
ylabel('Parallel Sensor Position - x\prime (pixels)');

パラレル ビーム - 投影データから頭部ファントムを復元

各復元のパラレル回転のインクリメント dtheta を上記合成投影で使用したインクリメントと一致させます。実際には、トランスミッターとセンサーの形状はわかりますが、ソース イメージ P の形状はわかりません。

次の 3 つの復元 (I1、I2、および I3) では、投影を行う角度を変えたときの影響を示します。I1 と I2 では、元のファントムで見られた特徴が不明瞭です。特に、各イメージの底部にある 3 つの楕円を見てください。I3 の結果は、元のイメージ P によく似ています。

I1 と I2 にはかなり不自然な結果が含まれていることに注意してください。これらの不自然な影響を回避するには、角度のステップを増やします。

% Constrain the output size of each reconstruction to be 
% the same as the size of the original image, |P|.
output_size = max(size(P));

dtheta1 = theta1(2) - theta1(1);
I1 = iradon(R1,dtheta1,output_size);

dtheta2 = theta2(2) - theta2(1);
I2 = iradon(R2,dtheta2,output_size);

dtheta3 = theta3(2) - theta3(1);
I3 = iradon(R3,dtheta3,output_size);

figure
montage({I1,I2,I3},'Size',[1 3])
title(['Reconstruction from Parallel Beam Projection ' ...
    'with 18, 24, and 90 Projection Angles'])

ファン ビーム - 合成投影の計算

ファンビーム形状を使用して合成投影を計算し、'FanSensorSpacing' をさまざまに変化させます。

D = 250; 
dsensor1 = 2;
F1 = fanbeam(P,D,'FanSensorSpacing',dsensor1);

dsensor2 = 1;
F2 = fanbeam(P,D,'FanSensorSpacing',dsensor2);

dsensor3 = 0.25;
[F3,sensor_pos3,fan_rot_angles3] = fanbeam(P,D, ...
    'FanSensorSpacing',dsensor3);

投影データ F3 を表示します。ファン回転角の範囲は 0 ～ 360 度で、同じ特徴が両側からサンプリングされているため、180 度のオフセットの位置で同じパターンが得られます。このファンビーム投影イメージの特徴を上記のパラレル ビーム投影イメージの同じ特徴と相互に関連付けることができます。 

imagesc(fan_rot_angles3,sensor_pos3,F3)
colormap(hot)
colorbar
xlabel('Fan Rotation Angle (degrees)')
ylabel('Fan Sensor Position (degrees)')

ファン ビーム - 投影データから頭部ファントムを復元

各復元のファンとセンサーの間隔を、合成投影の作成時に使用した間隔と一致させます。実際には、トランスミッターとセンサーの形状はわかりますが、ソース イメージ P の形状はわかりません。

'FanSensorSpacing' の値を変えると、各回転角で使用するセンサーの数が効率的に変更されます。これらのファンビーム復元では、同じ回転角を使用します。これは、その都度異なる回転角を使用するパラレル ビーム復元とは対照的です。

'FanSensorSpacing' は、fanbeam と ifanbeam を使用するときに制御できる複数のパラメーターの中の 1 つにすぎないことに注意してください。また、関数 fan2para と para2fan を使用して、パラレル ビーム投影データとファンビーム投影データを切り替えることもできます。

Ifan1 = ifanbeam(F1,D,'FanSensorSpacing',dsensor1, ...
    'OutputSize',output_size);
Ifan2 = ifanbeam(F2,D,'FanSensorSpacing',dsensor2, ...
    'OutputSize',output_size);
Ifan3 = ifanbeam(F3,D,'FanSensorSpacing',dsensor3, ...
    'OutputSize',output_size);

montage({Ifan1,Ifan2,Ifan3},'Size',[1 3])
title(['Reconstruction from Fan Beam Projection with '...
    '18, 24, and 90 Projection Angles'])