function I = invradonatt(B,Theta,p,A,iter);

%             I = invradonatt(B,Theta,p,A,iter);
%Etant donnée une liste d'angles Theta (exprimés en degrés entiers entre 0
%et 359), I désigne l'inverse (par projection orthogonales successives dans
%les directions d'angles correspondant à la liste Theta, prise dans cet
%ordre de p en p, méthode A.R.T.), d'un sinogramme B correspondant à la
%transformée de Radon atténuée bidimensionnelle d'une image I, en présence
%d'une carte d'atténuation A. Source I et atténuation A sont supposées
%données dans [-1,1]^2 et nulles hors du disque unité du plan. Le nombre
%iter désigne le nombre d'itérations. L'algorithme utilise la routine
%imrotate (option 'bilinear') de l'environnement MATLAB. Les contraintes de
%positivité de l'image sont prises en compte à chaque étape.


[M,N]=size(B);
l=length(Theta);
tau=2/M;
x=-1+tau/2:tau:1-tau/2;
init=zeros(M,M);
I=init;

    for n=1:iter
    for k=1:p:l
        
        S=B(:,1+Theta(k))';
        II=radonatt(I,Theta(k),A);
        S=S-II'; 
       
        Arot=imrotate(A,-Theta(k),'bilinear','crop');
        A2=exp(-2*radonatt_aux(Arot)*2/M) ;
        A1=exp(-radonatt_aux(Arot)*2/M) ;
       
        S0=(radonatt(A2,Theta(k),zeros(M,M)))';
        U=S./S0;
        UU=ones(M,1)*U;
        V=imrotate(UU,Theta(k),'bilinear','crop');
        I=I+A1.*V; 
 %Prise en compte des contraintes de positivité 
        for i=1:M
         for j=1:M
          if (I(i,j)<0)
            I(i,j)=0;
         end
        end
    end  
    end
end
