FDK算法程序

functiondR=FDK(x,y,z,u_off,v_off,du,dv,Beta,P,A,SDD,SAD,varargin);%%Inputs:%%Beta:theprojectionanglewhichis%%fixedforthewholefunction.%%P:2-DMatrixofsizenu*nv%%u_off,v_off:giveusthefarthestpointfrom%%theoriginofthedetector%%inhorizontalandvertical%%axesrespectively.%%{Inotherwords,ourdetectorisofsize%%[2*u_off]*[2*v_off]}(?)%%du,dv:stepsize%%x,y,z:The3-Dgridforreconstructingtheimage.%%A:ProjectionMatrix%%SDD:distancefromsourcetodetector%%radius:radiusofthefield-of-view%%Xs:???%%Phihat:1-DRampFilter,Differentoptionsoffilters%%-Willbeaddedlater.%%-willbeavailable(inFourierspace)%%-Ithastobeofthesamelengthasu%%-length(ghat)=nu%%Outputs:%%Afterfixingtheprojectionangle,theroutinereturnsthe%%valuecalculatedbyFDKequationforeachvoxel.%%Description:%%ReconstructsingleCone-Beamprojectionusing3DFeldkamp%Dimensionsofreconstructiongridnx=length(x);ny=length(y);nz=length(z);%Remark:TheremaybesomesavingsinfftcomputationifPistransposed[nv,nu]=size(P);u=(-u_off+[0:nu-1]')*du;v=(-v_off+[0:nv-1]')*dv;%Remark:ifPistransposed,'meshgrid'shouldbechangedto'ndgrid'[uu,vv]=meshgrid(u,v);weight=SDD./sqrt(SDD^2+vv.^2+uu.^2);%Remark:Pisover-writtentoconservememoryinfollowingcomputations.%TheinterpretationateachstageshouldbeclearfromcontextP=P.*weight;%Filtering:Phihat=varargin{1};if~isempty(Phihat),%Skiptobackprojectionisfilterisempty%Remark:Thiscalltofftcanlikelybeoptimisedusingthetranspose%ofPinsteadtocapitaliseonMatlab'scolumn-orientedstorage.P=fft(P,length(Phihat),2);%Rowszero-paddedautomatically%Remark:Thisloopcanlikelybevectorisedusingrepmatforj=1:nvP(j,:)=P(j,:).*Phihat;end;%Remark:Thiscalltofftcanlikelybeoptimisedusingthetranspose%ofPinsteadtocapitaliseonMatlab'scolumn-orientedstorage.P=real(ifft(P,[],2));P=P(:,1:nu);%Trimzero-paddingonrowsend%ofif-block%IncrementtoaddtoreconstructioninbackprojectionstagedR=zeros(ny,nx,nz);%Vectorisedcomputationofbackprojection[yy,xx,zz]=ndgrid(y,x,z);%Useprojectionmatrixtoprojectreconstruction(x,y,z)grid%intodetector(u,v)gridUV=A*[xx(:)';yy(:)';zz(:)';ones(1,nx*ny*nz)];？%Nearestneighbourinterpolationtofinddetectorcoordinates(u,v)%thatareclosesttoprojectionsofvoxels(x,y,z)U=round(UV(1,:)./UV(3,:))+1;？V=round(UV(2,:)./UV(3,:))+1;？%Identifyindicesofvoxelswithprojectionsstrictlywithindetectorgridind=find((U=1)&(V=1)&(U=nu)&(V=nv));%Remark:ThiswillhavetobechangedifPistransposedP_ind=(U(ind)-1)*nv+V(ind);？co=cos(Beta*pi/180);si=sin(Beta*pi/180);%Grid-dependentweightingfactors(onlycomputedforvoxelsneeded)W=(SAD./(SAD-co*xx(ind)-si*yy(ind))).^2;dR(ind)=W.*P(P_ind);returnfunctionA=Orbit2ProjectionMatrix(theta_G,du,dv,u_off,v_off,SDD,SAD)%%%%Description:%%OpenaniView3DIMGfileintoMatlabmatrixobject%%References:%disp('Enteringopen_IMG');%Assumption:SAD,SDD,IADallconstantsforidealisedcirculartrajectorytheta_G=theta_G(:);u_off=u_off(:);v_off=v_off(:);N_proj=length(theta_G);%SDD=SAD+IAD;A=zeros(3,4,N_proj);alpha=0;%Skewangle(fornon-rectangularpixels)fork=1:N_proj%focallengthf=SDD;C=[1/dutan(alpha)u_off(k)01/dvv_off(k)001];P=[-f0000-f000010];%transformationfromworldcoordinatestoimagingcoordinates%Ri=pmhRotationMatrix(90,0,theta_G(k)+180);%forMVwhichis90degreesoutofphaseRi=pmhRotationMatrix(90,0,theta_G(k)+90);Xs=[SAD*cosd(theta_G(k));SAD*sind(theta_G(k));0];Tw=[Ri-Ri*Xszeros(1,3)1];A(:,:,k)=C*P*Tw;%NormalizeProjectiondistancetoimageplaneA(:,:,k)=A(:,:,k)./A(3,4,k);endreturnfunctionR=pmhRotationMatrix(phi,theta,psi);%%Description:%%Computerotationmatrix%%%%[cos(theta)*cos(psi)cos(phi)*cos(psi)Theta=theta*pi/180;Phi=phi*pi/180;Psi=psi*pi/180;R=[1000cos(Phi)sin(Phi)0-sin(Phi)cos(Phi)]...*...[cos(Theta)0-sin(Theta)010sin(Theta)0cos(Theta)]...*...[cos(Psi)sin(Psi)0-sin(Psi)cos(Psi)0001];return