程
MATLAB 源码：
序生成128的流图：
clc;
clear all; close all;
% load mri %载入mri数据，是matlab自带库 % ph = squeeze(D); %压缩载入的数据D，并赋值给ph 输入图片数字ph = phantom3d(128);
prompt={'Enter the Piece num(1 to 128):'}; %提示信息“输入1到27的片的数字” name='Input number'; %弹出框名称 defaultanswer={'1'}; %默认数字 numInput=inputdlg(prompt,name,1,defaultanswer) %弹出框，并得到用户的输入信息 对图片信息进行预处P= squeeze(ph(:,:,str2num(cell2mat(numInput))));%将用户的输入信息转换成数字，并从ph中得到相应的片信息P
imshow(P) %展示图片P
D = 250; %将D赋值为250,是从扇束顶点到旋转中心的像素距离。
用函数fanbeam进行映dsensor1 = 2; %正实数指定扇束传感器的间距2 F1 = fanbeam(P,D,'FanSensorSpacing',dsensor1); %通过P，D等计算扇束的数据值
射，得到扇束的数据，并dsensor2 = 1; %正实数指定扇束传感器的间距1 F2 = fanbeam(P,D,'FanSensorSpacing',dsensor2); %通过P，D等计算扇束的数据值 dsensor3 = 0.25 %正实数指定扇束传感器的间距0.25
用函数ifanbeam根据扇[F3, sensor_pos3, fan_rot_angles3] = fanbeam(P,D,... 'FanSensorSpacing',dsensor3); %通过P，D等计算扇束的数据值，并得到扇束传感器的位
束投影数据重建图像，并置sensor_pos3和旋转角度fan_rot_angles3 figure, %创建窗口
imagesc(fan_rot_angles3, sensor_pos3, F3) %根据计算出的位置和角度展示F3的图片 colormap(hot); %设置色图为hot 计算重建图像和原图的colorbar; %显示色栏 xlabel('Fan Rotation Angle (degrees)') %定义x坐标轴 ylabel('Fan Sensor Position (degrees)') %定义y坐标轴
output_size = max(size(P)); %得到P维数的最大值，并赋值给output_size Ifan1 = ifanbeam(F1,D, ... 结束 'FanSensorSpacing',dsensor1,'OutputSize',output_size); %根据扇束投影数据F1及D等数据重建图像 figure, imshow(Ifan1) %创建窗口，并展示图片Ifan1 title('图一');
disp('图一和原图的性噪比为：'); result=psnr1(Ifan1,P); Ifan2 = ifanbeam(F2,D, ...
'FanSensorSpacing',dsensor2,'OutputSize',output_size);
%根据扇束投影数据F2及D等数据重建图像 figure, imshow(Ifan2) %创建窗口，并展示图片Ifan2 disp('图二和原图的性噪比为：'); result=psnr1(Ifan2,P); title('图二');
Ifan3 = ifanbeam(F3,D, ...
'FanSensorSpacing',dsensor3,'OutputSize',output_size);
%根据扇束投影数据F3及D等数据重建图像 figure, imshow(Ifan3) %创建窗口，并展示图片Ifan3 title('图三');
disp('图三和原图的性噪比为：'); result=psnr1(Ifan3,P);
function [p,ellipse]=phantom3d(varargin)
%PHANTOM3D Three-dimensional analogue of MATLAB Shepp-Logan phantom % P = PHANTOM3D(DEF,N) generates a 3D head phantom that can % be used to test 3-D reconstruction algorithms. %
% DEF is a string that specifies the type of head phantom to generate. % Valid values are: %
% 'Shepp-Logan' A test image used widely by researchers in % tomography
% 'Modified Shepp-Logan' (default) A variant of the Shepp-Logan phantom % in which the contrast is improved for better % visual perception. %
% N is a scalar that specifies the grid size of P. % If you omit the argument, N defaults to 64. %
% P = PHANTOM3D(E,N) generates a user-defined phantom, where each row % of the matrix E specifies an ellipsoid in the image. E has ten columns, % with each column containing a different parameter for the ellipsoids: %
% Column 1: A the additive intensity value of the ellipsoid % Column 2: a the length of the x semi-axis of the ellipsoid % Column 3: b the length of the y semi-axis of the ellipsoid % Column 4: c the length of the z semi-axis of the ellipsoid % Column 5: x0 the x-coordinate of the center of the ellipsoid % Column 6: y0 the y-coordinate of the center of the ellipsoid % Column 7: z0 the z-coordinate of the center of the ellipsoid % Column 8: phi phi Euler angle (in degrees) (rotation about z-axis) % Column 9: theta theta Euler angle (in degrees) (rotation about x-axis) % Column 10: psi psi Euler angle (in degrees) (rotation about z-axis) %
% For purposes of generating the phantom, the domains for the x-, y-, and % z-axes span [-1,1]. Columns 2 through 7 must be specified in terms % of this range. %
% [P,E] = PHANTOM3D(...) returns the matrix E used to generate the phantom. %
% Class Support % -------------
% All inputs must be of class double. All outputs are of class double. %
% Remarks % -------
% For any given voxel in the output image, the voxel's value is equal to the
% sum of the additive intensity values of all ellipsoids that the voxel is a % part of. If a voxel is not part of any ellipsoid, its value is 0. %
% The additive intensity value A for an ellipsoid can be positive or negative; % if it is negative, the ellipsoid will be darker than the surrounding pixels. % Note that, depending on the values of A, some voxels may have values outside % the range [0,1]. % % Example % -------
% ph = phantom3d(128);
% figure, imshow(squeeze(ph(64,:,:))) %
% Copyright 2005 Matthias Christian Schabel (matthias @ stanfordalumni . org) % University of Utah Department of Radiology % Utah Center for Advanced Imaging Research % 729 Arapeen Drive
% Salt Lake City, UT 84108-1218 %
% This code is released under the Gnu Public License (GPL). For more information, % see : http://www.gnu.org/copyleft/gpl.html %
% Portions of this code are based on phantom.m, copyrighted by the Mathworks %
[ellipse,n] = parse_inputs(varargin{:}); p = zeros([n n n]);
rng = ( (0:n-1)-(n-1)/2 ) / ((n-1)/2); [x,y,z] = meshgrid(rng,rng,rng);
coord = [flatten(x); flatten(y); flatten(z)]; p = flatten(p);
for k = 1:size(ellipse,1)
A = ellipse(k,1); % Amplitude change for this ellipsoid asq = ellipse(k,2)^2; % a^2 bsq = ellipse(k,3)^2; % b^2 csq = ellipse(k,4)^2; % c^2 x0 = ellipse(k,5); % x offset y0 = ellipse(k,6); % y offset z0 = ellipse(k,7); % z offset
phi = ellipse(k,8)*pi/180; % first Euler angle in radians theta = ellipse(k,9)*pi/180; % second Euler angle in radians psi = ellipse(k,10)*pi/180; % third Euler angle in radians cphi = cos(phi); sphi = sin(phi); ctheta = cos(theta); stheta = sin(theta); cpsi = cos(psi); spsi = sin(psi);
% Euler rotation matrix
CT图像三维重建（附源码）
