图像重建迭代算法matlab

clear all;
clc;
I=phantom（256）;
imshow（I）;
theta1=0:179;
[R1xp1] = radon（I theta1）;
g=zeros（165536）;
for i=1:180
    if （mod（（i-1）2）==0）
        Q（i）=（i+1）/2;
        O（:i）=R1（:（（i+1）/2））;
    else 
        Q（i）=90+i/2;
        O（:i）=R1（:（90+i/2））;
    end
end
for y=1:180
    for x=1:367
        i=（y-1）*367+x;
        p（i）=O（xy）;
    end
end
for M=1:180*367
    a=fix（（M-1）/367）+1;
    theta=Q（a）;
    t=mod（（M-1）367）+1;
    A=cos（（theta-1）*pi/180）;
    B=sin（（theta-1）*pi/180）;
    C=-（t-184）;
    if （theta~=1）&&（theta~=91）
    i=1;    
    x=-128;
    y=-（A/B）*x-C/B;
    if（（y<=128）&&（y>=-128））
        tem（i）=x;
        i=i+1;
        tem（i）=y;
        i=i+1;
    end
    x=128;
    y=-（A/B）*x-C/B;
    if（（y<=128）&&（y>=-128））
        tem（i）=x;
        i=i+1;
        tem（i）=y;
        i=i+1;
    end;
    y=-128;
    x=-（B/A）*y-C/A;
    if（（x<=128）&&（x>=-128））
        tem（i）=x;
        i=i+1;
        tem（i）=y;
        i=i+1;
    end
    y=128;
     x=-（B/A）*y-C/A;
     if（（x<=128）&&（x>=-128））
        tem（i）=x;
        i=i+1;
        tem（i）=y;
        i=i+1;
     end
    end
    if（theta==91）
            tem（1）=-128;
            tem（2）=-C;
            tem（3）=128;
            tem（4）=-C;
            
    end
        if（theta==1）
            tem（1）=-C;
            tem（2）=-128;
            tem（3）=-C;
            tem（4）=128;
            
        end
         if（tem（1）>tem（3））
             tem（1）=ceil（tem（1））;
             tem（3）=floor（tem（3））;
             temp=tem（1）;
             tem（1）=tem（3）;
             tem（3）=temp;
        
         end
         if（tem（3）>tem（1））
             tem（3）=ceil（tem（3））;
             tem（1）=floor（tem（1））;
             
         end
         if（tem（2）>tem（4））
             tem（2）=ceil（tem（2））;
             tem（4）=floor（tem（4））;
             temp=tem（2）;
             tem（