滤波反投影matlab仿真程序

%%Accurate image reconstruction from few-views and limited-angle data indivergent-beam CT % % clear; clear all; M=256; N=256; I = phantom(M,N); figure(1); imshow(I); figure(2); I_st=zeros(M,N); for s=2:M for t=2:N I_st(s,t)= sqrt((I(s,t)-I(s-1,t))^2+((I(s,t)-I(s,t-1))^2)); end end I_st(1,:)=0.5*(I_st(2,:)+I_st(3,:)); I_st(:,1)=0.5*(I_st(:,2)+I_st(:,4)); imshow(I_st); xLength=20; yLength=20; xResolution=256; yResolution=256; xLengthPerUnit=xLength/xResolution; yLengthPerUnit=yLength/yResolution; xCenter=xLength/2; yCenter=yLength/2; xCoodinate=(([1:M]-1)*xLengthPerUnit-xCenter)'; yCoodinate=(([N:-1:1])*yLengthPerUnit-yCenter)'; xSource=-40; ySource= 0; coordinateSource=[xSource;ySource]; detecLength=41.3; detecResolution=512; unitDis=detecLength/(detecResolution-1); xDetector=40; yDetector=([1:detecResolution]-1)*unitDis-detecLength/2; coordinateDetector=zeros(2,detecResolution); coordinateDetector(1,:)=xDetector; coordinateDetector(2,:)=yDetector; disSourceToDetec=80; disFocal=40; theta=0; numSample=360; deltaTheta=2*pi/numSample; projMatrix=zeros(numSample,detecResolution); I_column=zeros(xResolution*yResolution,1); for i=1:yResolution for j=1:xResolution index=(i-1)*xResolution+j; I_column(index,1)=I(i,j); end end for sample=1:numSample theta=(sample-1)*deltaTheta; A=[cos(theta),sin(theta);-sin(theta),cos(theta)]; rotatedSource=A*coordinateSource; end for sample=1:numSample theta=(sample-1)*deltaTheta; A=[cos(theta),sin(theta);-sin(theta),cos(theta)]; rotatedSource=A*coordinateSource; rotatedDetector=A*coordinateDetector; R=zeros(detecResolution,xResolution*yResolution); for detecIndex=1:detecResolution; k=(rotatedDetector(2,detecIndex)-rotatedSource(2,1))/(rotatedDetector(1,detecIndex)-rotatedSource(1,1)); b=rotatedSource(2,1)-k*rotatedSource(1,1); for i=1:xResolution*yResolution xNum=floor((i-1)/xResolution)+1; yNum=rem(i,xResolution); if (yNum==0) yNum=xResolution; end rect=[xCoodinate(xNum,1), xCoodinate(xNum,1)+xLengthPerUnit,xCoodinate(xNum,1),xCoodinate(xNum,1)+xLengthPerUnit; yCoodinate(yNum,1), yCoodinate(yNum,1),yCoodinate(yNum,1)-yLengthPerUnit,yCoodinate(yNum,1)-yLengthPerUnit]; judge=sign(rect(2,:)-k*rect(1,:)-b); if(rem(i,65536)==0) mm=1; end if (abs(sum(judge))>2) continue; else if(sum((judge==1))==2) if(judge(1,1)==judge(1,2)) x1=rect(1,1); x2=rect(1,2); y1=k*x1+b; y2=k*x2+b; end if(judge(1,1)==judge(1,3)) y1=rect(2,1); y2=rect(2,3); x1=(y1-b)/k; x2=(y2-b)/k; end disIntersect=sqrt((y1-y2)^2+(x1-x2)^2); else if(sum((judge==1))==1) index=find(judge==1); end if(sum((judge==-1))==1) index =find(judge==-1); end x1=rect(1,index); y1=k*x1+b; y2=rect(2,index); x2=(y2-b)/k; disIntersect=sqrt((y1-y2)^2+(x1-x2)^2); end R(detecIndex,i)=disIntersect; end end % if(rem(detecIndex,10)==0) % mm=1; % end if(rem(detecIndex,50)==0) disp('Number of Detector ') detecIndex end end projMatrix(sample,:)=R*I_column; disp('Sample') sample end % for sample=1:numSample % theta=(sample-1)*deltaTheta; % A=[cos(theta),sin(theta);-sin(theta),cos(theta)]; % rotatedSource=A*coordinateSource; % rotatedDetector=A*coordinateDetector; % for detecIndex=1:detecResolution; % tempProj=0; % k=(rotatedDetector(2,detecIndex)-rotatedSource(2,1))/(rotatedDetector(1,detecIndex)-rotatedSource(1,1)); % b=rotatedSource(2,1)-k*rotatedSource(1,1); % for xNum=1: xResolution; % for yNum=1:yResolution % rect=[xCoodinate(xNum,1), xCoodinate(xNum,1)+xLengthPerUnit,xCoodinate(xNum,1),xCoodinate(xNum,1)+xLengthPerUnit; % yCoodinate(yNum,1), yCoodinate(yNum,1),yCoodinate(yNum,1)-yLengthPerUnit,yCoodinate(yNum,1)-yLengthPerUnit]; % judge=sign(rect(2,:)-k*rect(1,:)-b); % if (abs(sum(judge))>2) % continue; % else % if(sum((judge==1))==2) % if(judge(1,1)==judge(1,2)) % x1=rect(1,1); % x2=rect(1,2); % y1=k*x1+b; % y2=k*x2+b; % end % if(judge(1,1)==judge(1,3)) % y1=rect(2,1); % y2=rect(2,3); % x1=(y1-b)/k; % x2=(y2-b)/k; % end % disIntersect=sqrt((y1-y2)^2+(x1-x2)^2); % tempProj=tempProj+disIntersect*I(xNum,yNum); % % else % if(sum((judge==1))==1) % index=find(judge==1); % end % if(sum((judge==-1))==1) % index =find(judge==-1); % end % x1=rect(1,index); % y1=k*x1+b; % y2=rect(2,index); % x2=(y2-b)/k; % disIntersect=sqrt((y1-y2)^2+(x1-x2)^2); % % tempProj=tempProj+disIntersect*I(xNum,yNum); % end % % end % end % % end % % if ( rem(detecIndex,50)==0) % % m=1; % % end % % % % projMatrix(sample,detecIndex)=tempProj; % end % end % figure(3); imshow(projMatrix); % % %