function f=retroproj(R,tau,sigma,B); 

%              f=retroproj(R,tau,sigma,B); 
%
%Inverse la transformée de Radon A d'une image 2D  
%I=I(-R:tau:R-tau,-R:tau:R-tau) par rétroprojection ; le pas angulaire est
%fixé égal a sigma ; l'entrée B est la transformée de Radon du signal s
%calculee soit selon B=radon0(R,tau,sigma,I), soit selon 
%B=radon1(R,tau,sigma,I)(éventuellement bruitée), suivant que l'on utilise
%ou non la routine imrotate de l'environnement MATLAB.
%
%On se réfèrera pour plus de précisions à la documentation fournie dans le 
%fichier doc-retroproj.pdf. 

B=B'; 
[M,MM]=size(B);

% Calcule la dérivée du spectrogramme par rapport à s. 

C=zeros(M,MM); 
C(:,1)=(B(:,2)-B(:,1))/tau;
C(:,MM)=(B(:,MM)-B(:,MM-1))/tau;
B1=B(:,3:MM);
B2=B(:,1:MM-2);
CC=B1-B2;
C(:,2:MM-1)=CC/(2*tau); 

% Calcule la transformée de Hilbert (en s) après filtrage dans le domaine
% des fréquences (fenêtrage de Hamming). 

h=hamming(MM);
H=zeros(M,MM);
for k=1:M
   H(k,:)=hamming(MM)';
end
U=H.*spectremat(C,R);
MMM=fix(MM/2);
w1=- sqrt(-1)* ones(M,MMM);
w2=sqrt(-1)*ones(M,MM-MMM);
w=[w1 w2]; 
UU=w.*U;
Omega=pi/tau; 
BB=real(ispectremat(UU,Omega));

% Calcule l'image par R dièse (voir documentation).  

theta1=0:sigma:pi-sigma;
theta2=pi:sigma:2*pi-sigma;
theta=[theta1 theta2]; 
a=ones(1,MM);
x=-R:tau:R-tau;
x=x';
XX=x*a;
YY=a'*x';
L=zeros(MM,MM);
Laux=zeros(MM,MM);
for i=1:M;
   Q=(cos(theta(i))*XX+ sin(theta(i))*YY+R*ones(MM))/tau;
   Q=max(Q,zeros(MM));
   Q=min(Q,(MM-1)*ones(MM));
   QQ=round(Q+ones(MM));
   for k=1:MM
      for l=1:MM
         Laux(k,l)=BB(i,QQ(k,l));
      end
   end
   L=L+Laux;
end
f=(1/(4*pi))*L*sigma;
f=f.*disque(R,tau);