FDK重建算法的代码（matlab）

function [vol_tmp,H] = iconebeamdemo2007(dir,flprefix,theta,interp,filter,d,N,step,Ro,range_slices,cut) % author: Gianni Schena Sept. 2004 (schena@univ.trieste.it - I welcome suggestions !) % Reconstructs image/s for cone beam (CB) geometry according to % the classical Feldkamp et al. algorithm. Individual slices are stored in a volume % % FDK10_CB reconstructs the 3D volume image vol from planar projection % data. The routine assumes that the center of rotation % is the center point of the projections. % % THETA describes the angles (in degrees) at which the projections % were taken. It it is a vector containing angle between % projections. THETA is a vector and must contain angles with % equal spacing between them. % % The routine uses the filtered backprojection.The filter is designed directly % in the frequency domain and then multiplied by the FFT of the % projections. The projections are zero-padded to a power of 2 % before filtering to prevent spatial domain aliasing and to % speed up the FFT. % % INTERP specifies the type of interpolation to use in the % backprojection. The available options are listed in order % of increasing accuracy and computational complexity: % % 'nearest' - nearest neighbor interpolation % 'linear' - linear interpolation (suggested) % % FILTER specifies the filter to use for frequency domain filtering. % FILTER is a string that specifies any of the following standard % filters: % % 'Ram-Lak' The cropped Ram-Lak or ramp filter (default). The % frequency response of this filter is |f|. Because % this filter is sensitive to noise in the projections, % one of the filters listed below may be preferable. % 'Shepp-Logan' The Shepp-Logan filter multiplies the Ram-Lak % filter by a sinc function. % 'Cosine' The cosine filter multiplies the Ram-Lak filter % by a cosine function. % 'Hamming' The Hamming filter multiplies the Ram-Lak filter % by a Hamming window. % 'Hann' The Hann filter multiplies the Ram-Lak filter by % a Hann window. % % D is a scalar in the range (0,1] that modifies the filter by % rescaling its frequency axis. The default is 1. If D is less % than 1, the filter is compressed to fit into the frequency range % [0,D], in normalized frequencies; all frequencies above D are set % to 0. % % N is a scalar that specifies the number of rows and columns in the % reconstructed 2D slice images. If N is not specified, the size is determined % % [vol,H] = FDKX... (...) returns also the frequency response of the filter % in the vector H. % % Ro is the orbital radius i.e. the distance source to object % Ds2dct is the distance between source and planar detector % Ds2dct = Ro + Do2dct = Ds2obj + Do2dct % pxlsz is the size of the (square) pixel of the detector dct % Class Support : All input arguments must be of class double. % The output arguments are of class double. % References: % A. C. Kak, Malcolm Slaney, "Principles of Computerized Tomographic % Imaging", IEEE Press 1988. Chapter 3 (available by downloading) %[p,theta,filter,d,interp,N] = parse_inputs(varargin{:}); thetag=theta; theta = pi*theta/180; % from degrees to radiants A=read_1rad(1.,dir,flprefix); % read in image frame figure , imshow(A,[]); title(' first projection') % disp first image file %global flat , global dark % get size of the data files [nr_p , nc_p] = size(A) ; % get the number of rows and colums of the data files % range_slices=range_slices(range_slices>-nr_p/2 & range_slices<nr_p/2); % remove slices out of range % N is a scalar that specifies the number of rows and columns in the % reconstructed 2D slice images. If N is not specified, the size is determined % If the user didn't specify the size of the reconstruction, so % deduce it from the length of projections % nr_p; % sinogram size i.e. one size of the projection if N==0 || isempty(N) % if N is not given then ... N = 2*floor( nc_p /(2*sqrt(2)) ); % take the default value end imgDiag = 2*ceil(N/sqrt(2))+1; % largest distance through image. %np = length(theta); W=pix_weight(nr_p,nc_p,Ro,0,0); % weights to account for a planar detector %size(W) % Design the filter H len=nc_p; H = designFilter(filter, len, d); LH=length(H); % length filter in the line vector H % Define the x & y axes for the reconstructed image so that the origin xax = (1:N)-ceil(N/2); x = repmat(xax, N, 1); y= repmat(xax', 1, N); %y = rot90(x); % x coordinates, the y coordinates are rot90(x) mask= ((x-0).^2+(y-0).^2) <= (N/2).^2 ; % mask the inner circle of diam N mask=true(N); % no mask : all the pixels are to be reconstructed maskV=mask(:); % mask in a vector costheta = cos(theta); sintheta = sin(theta); ctrIdx = ceil(len/2); % index of the center of the projections (on the hor axis) % range slices to be recons if exist('range_slices','var') & ~isempty(range_slices) slices=range_slices; % slices to be recons if isempty( find(slices==0) ) % se non c'e' la slice centrale aggiungila slices=[slices ,0]; slices=sort(slices) ; end else % if no slices then recns the slice z=0 only slices=[0];% [-floor(nc_p/2):step:+floor(nc_p/2)]; end cnt_slc=find(slices==0); slices %vol = repmat(0,[N,N,length(slices)]) ;% allocate 3D volume matrix vol_tmp=repmat(single(0),[(N.^2),length(slices)]); % allocation szproj=[nr_p,nc_p;]; % size(proj); % % Filtered Back-projection - i.e. from proj pixel (u,v) to voxel (x,y,z) if strcmp(interp, 'nearest neighbor') x=x(maskV); y=y(maskV); % mask ! nw_slices = real(sort(complex(slices))) ; % sort and reordering for faster reconstruction % eg. +5 -5 + 6 -6 ...it allows to change the sign of v rather than recalculating v ZFlag= abs(nw_slices(2:end)) == nw_slices(1:(end-1)) ; % logical variable ZFlag=real([0, ZFlag] ); % true if +z is followed by -z and one can exploit symmetry for ipos=1:length(slices) % right positioning the re-orded slices in the volume iposz_v(ipos)=find( nw_slices(ipos)== slices) ; % slice correspondence end % data type and array allocation statements proj=repmat(0,nr_p,nc_p); tmp=zeros(nnz(mask),1); tmp2=zeros(N.^2,1); % allocate memory to temporary arrays % tmp1 v0=zeros([size(x)]); v1=v0; u1=v1; % u=u1; v=u; good_u=logical(v0); idx=zeros(size(v0)); % allocation for i=1:length(theta) % loop over the projections % i dgri=thetag(i); % theta in degrees proj=read_1rad(i,dir,flprefix); % read in the i-th radiograph proj = radio2proj(proj,W,H,len); % tranforms radiog. into projection and multiply by W %%% proj(:,length(H))=0; % Zero pad projections for for ir = 1:nr_p; % filtering row by row with row vector H proj(ir,:)= real (ifft ( fft(proj(ir,:)).*H ) ); % fast frequency domain filtering end proj(:,len+1:end) = [] ; % Truncate the filtered projections %%% s = x.*costheta(i) + y.*sintheta(i); t = -x.*sintheta(i) + y.*costheta(i); % Goddard & Trepanier % s = x.*costheta(i) + y.*sintheta(i); % also in matlab iradon R_s = Ro-s; u=round( Ro*t./(R_s) + ctrIdx) ; % u is the 1st projection pixel coordinate good_u = (u>0 & u<(nc_p+1)) ; % only valid pixels -i.e. within the projection frame lengoodu=length(good_u); u1=u(good_u);% inv_R_s_quad=1./(R_s(good_u).^2); Ro_dev_R_s=Ro./(R_s(good_u)); % for generally .* is faster than ./ iz=0; id_1_term=(u1-1)*szproj(1); % does not depend up on v