CT图像三维重建(附源码)

程序流图： MATLAB 源码： 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,是从扇束顶点到旋转中心的像素距离。 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 [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 : % %  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 alpha = [cpsi*cphi-ctheta*sphi*spsi  cpsi*sphi+ctheta*cphi*spsi  spsi*stheta; -spsi*cphi-ctheta*sphi*cpsi  -spsi*sphi+ctheta*cphi*cpsi cpsi*stheta; stheta*sphi                  -stheta*cphi                ctheta];        % rotated ellipsoid coordinates coordp = alpha*coord; idx = find((coordp(1,:)-x0).^2./asq + (coordp(2,:)-y0).^2./bsq + (coordp(3,:)-z0).^2./csq <= 1); p(idx) = p(idx) + A; end p = reshape(p,[n n n]); return; function out = flatten(in) out = reshape(in,[1 prod(size(in))]); return; function [e,n] = parse_inputs(varargin) %  e is the m-by-10 array which defines ellipsoids %  n is the size of the phantom brain image n = 128;    % The default size e = []; defaults = {'shepp-logan', 'modified shepp-logan', 'yu-ye-wang'}; for i=1:nargin if ischar(varargin{i})        % Look for a default phantom def = lower(varargin{i}); idx = strmatch(def, defaults); if isempty(idx) eid = sprintf('Images:%s:unknownPhantom',mfilename); msg = 'Unknown default phantom selected.'; error(eid,'%s',msg); end switch defaults{idx} case 'shepp-logan' e = shepp_logan; case 'modified shepp-logan' e = modified_shepp_logan; case 'yu-ye-wang' e = yu_ye_wang; end elseif numel(varargin{i})==1 n = varargin{i};            % a scalar is the image size elseif ndims(varargin{i})==2 && size(varargin{i},2)==10 e = varargin{i};            % user specified phantom else eid = sprintf('Images:%s:invalidInputArgs',mfilename); msg = 'Invalid input arguments.'; error(eid,'%s',msg); end end % ellipse is not yet defined if isempty(e)                    e = modified_shepp_logan; end return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %  Default head phantoms:  % %%%%%%%%%%%%%%%%%%%%%%%%%%%%% function e = shepp_logan e = modified_shepp_logan;