《FDK算法程序.docx》由会员分享，可在线阅读，更多相关《FDK算法程序.docx（6页珍藏版）》请在三一办公上搜索。
1、FDK算法程序function dR = FDK(x,y,z,u_off,v_off,du,dv,Beta,P,A,SDD,SAD,varargin); % Inputs: % Beta: the projection angle which is % fixed for the whole function. % P: 2-D Matrix of size nu*nv % u_off, v_off: give us the farthest point from % the origin of the detector % in horizontal and vertical % axes
2、respectively. % In other words, our detector is of size % 2*u_off*2*v_off(?) % du, dv: step size % x,y,z: The 3-D grid for reconstructing the image. % A: Projection Matrix % SDD: distance from source to detector % radius: radius of the field-of-view % Xs: ? % Phihat: 1-D Ramp Filter, Different optio
3、ns of filters % -Will be added later. % -will be available (in Fourier space) % -It has to be of the same length as u % -length(ghat)=nu % Outputs: % After fixing the projection angle, the routine returns the % value calculated by FDK equation for each voxel. % Description: % Reconstruct single Cone
4、-Beam projection using 3D Feldkamp % Dimensions of reconstruction grid nx = length(x); ny = length(y); nz = length(z); % Remark: There may be some savings in fft computation if P is transposed nv,nu = size(P); u = (-u_off + 0:nu-1)*du; v = (-v_off + 0:nv-1)*dv; % Remark: if P is transposed, meshgrid
5、 should be changed to ndgrid uu,vv = meshgrid(u,v); weight = SDD./sqrt(SDD2 + vv.2 + uu.2); % Remark: P is over-written to conserve memory in following computations. % The interpretation at each stage should be clear from context P = P .* weight; % Filtering: Phihat = varargin1; if isempty(Phihat),
6、% Skip to backprojection is filter is empty % Remark: This call to fft can likely be optimised using the transpose % of P instead to capitalise on Matlabs column-oriented storage. P = fft( P, length(Phihat), 2 ); % Rows zero-padded automatically % Remark: This loop can likely be vectorised using rep
7、mat for j=1:nv P(j,:) = P(j,:) .* Phihat ; end; % Remark: This call to fft can likely be optimised using the transpose % of P instead to capitalise on Matlabs column-oriented storage. P = real( ifft( P, , 2 ) ); P = P(:,1:nu); % Trim zero-padding on rows end % of if-block % Increment to add to recon
8、struction in backprojection stage dR = zeros(ny,nx,nz); % Vectorised computation of backprojection yy, xx, zz = ndgrid( y, x, z); % Use projection matrix to project reconstruction (x,y,z) grid % into detector (u,v) grid UV = A * xx(:); yy(:); zz(:); ones(1,nx*ny*nz) ;？ % Nearest neighbour interpolat
9、ion to find detector coordinates (u,v) % that are closest to projections of voxels (x,y,z) U = round( UV(1,:)./UV(3,:) ) + 1 ;？ V = round( UV(2,:)./UV(3,:) ) + 1 ;？ % Identify indices of voxels with projections strictly within detector grid ind = find( (U=1) & (V=1) & (U=nu) & (V=nv) ); % Remark: Th
10、is will have to be changed if P is transposed P_ind = ( U(ind) - 1 )*nv + V(ind) ;？ co = cos(Beta*pi/180); si = sin(Beta*pi/180); % Grid-dependent weighting factors (only computed for voxels needed) W = ( SAD ./ ( SAD - co*xx(ind) - si*yy(ind) ) ).2; dR( ind ) = W .* P( P_ind ); return function A =
11、Orbit2ProjectionMatrix(theta_G,du,dv,u_off,v_off,SDD,SAD) % % Description: % Open an iView3D IMG file into Matlab matrix object % References: %disp(Entering open_IMG); % Assumption: SAD, SDD, IAD all constants for idealised circular trajectory theta_G = theta_G(:); u_off = u_off(:); v_off = v_off(:)
12、; N_proj = length(theta_G); %SDD = SAD+IAD; A = zeros(3,4,N_proj); alpha = 0; % Skew angle (for non-rectangular pixels) for k=1:N_proj % focal length f = SDD; C = 1/du tan(alpha) u_off(k) 0 1/dv v_off(k) 0 0 1 ; P = -f 0 0 0 0 -f 0 0 0 0 1 0 ; % transformation from world coordinates to imaging coord
13、inates % Ri = pmhRotationMatrix(90,0,theta_G(k)+180); % for MV which is 90 degrees out of phase Ri = pmhRotationMatrix(90,0,theta_G(k)+90); Xs = SAD*cosd(theta_G(k); SAD*sind(theta_G(k); 0 ; Tw = Ri -Ri*Xs zeros(1,3) 1 ; A(:,:,k) = C*P*Tw; % Normalize Projection distance to image plane A(:,:,k) = A(
14、:,:,k)./A(3,4,k); end return function R = pmhRotationMatrix(phi,theta,psi); % Description: % Compute rotation matrix % % cos(theta)*cos(psi) cos(phi)*cos(psi) Theta = theta*pi/180; Phi = phi*pi/180; Psi = psi*pi/180; R = 1 0 0 0 cos(Phi) sin(Phi) 0 -sin(Phi) cos(Phi) . * . cos(Theta) 0 -sin(Theta) 0 1 0 sin(Theta) 0 cos(Theta) . * . cos(Psi) sin(Psi) 0 -sin(Psi) cos(Psi) 0 0 0 1 ; return
