%映射和重建
%產生測試圖像并顯示
P = phantom(256);
imshow(P)

%計算投影數據
%設定幾何關系為“線形”，分別給定光束數目為18，36和90。
theta1 = 0:10:170; 
[R1,xp] = radon(P,theta1); 
num_angles_R1 = size(R1,2)% num_angles_R1 = 18

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

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

%設定幾何關系為“弧度”，分別給定光束密度為2弧度、1弧度以及0.25弧度。
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);

%顯示投影數據
%幾何關系為“線形”
figure, imagesc(theta3,xp,R3)
colormap(hot)
colorbar
xlabel('Parallel Rotation Angle - \theta (degrees)'); 
ylabel('Parallel Sensor Position - x\prime (pixels)');

%幾何關系為“弧度”
figure, imagesc(fan_rot_angles3, sensor_pos3, F3)
colormap(hot); colorbar
xlabel('Fan Rotation Angle (degrees)')
ylabel('Fan Sensor Position (degrees)')

%重建圖像
%幾何關系為“線形”
dtheta1 = theta1(2) - theta1(1);
I1 = iradon(R1,dtheta1,output_size);
figure, imshow(I1)

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

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

%幾何關系為“弧度”
output_size = max(size(P));
Ifan1 = ifanbeam(F1,D, 'FanSensorSpacing',dsensor1,…
          'OutputSize',output_size);
figure, imshow(Ifan1)

Ifan2 = ifanbeam(F2,D, 'FanSensorSpacing',dsensor2,…
          'OutputSize',output_size);
figure, imshow(Ifan2)

Ifan3 = ifanbeam(F3,D, 'FanSensorSpacing',dsensor3,…
          'OutputSize',output_size);
figure, imshow(Ifan3)