fourier.m

几个图象处理程序

clear
sl=phantom;
the = 0:179;
[R,xp] = radon(sl,the);
[width,length]=size(R);      
width1=width;
length1=width;
rec=zeros(width1,length1);  
flou=zeros(width,length);
for i=1:length
flou(:,i) =fftshift(R(:,i)); 
flou(:,i)=fft(flou(:,i));  
end
o=round(width1/2);
 for i=1:width1
    for j=1:length1
        x=j-o;
        y=o-i;        
            if(x==0)       
                  if(y>=0)
                     theta=90;
                 else
                     theta=-90;
                  end
            else
                theta=(atan(y/x))*180/pi;
            end            
            if(x<=0)
                theta=180+theta;
            elseif(y<0)
                theta=2*180+theta;
            end            
           r=sqrt(x^2+y^2);
            if(theta<180)
            R=r+1;
            else
              R=width1-r;
              theta=theta-180;
            end
           theta=theta+1;
   %二维平面插值        
           theta1=floor(theta); theta2=ceil(theta);
           %当theta为整数的时候，theta1与theta2 相等
           R1=floor(R);R2=ceil(R);
           %当R为整数的时候，R1与R2 相等
           Re=abs(R-R1); 
           thet=abs(theta-theta1);
    %模方向的插补  
    if theta2<=180 && R2<=367
          Z1=(1-Re)*flou(R1,theta1)+Re*flou(R2,theta1);
          Z2=(1-Re)*flou(R1,theta2)+Re*flou(R2,theta2);
    end
    
    %角度方向的插补
          rec(i,j)=(1-thet)*Z1+thet*Z2;    %若theta与R都为整数，则rec(i,j)=flou(R,theta)
                                           %即取扇形区域左上角的一个样本点         
        
    end
end
rec=ifft2((rec));
rec2=abs(fftshift(rec));
figure;
imshow(rec2);