% Archivo principal y original con todos los algoritmos y opciones posibles

close all;
clear all;

% Cargamos las imágenes y las preparamos
oldR = imread('figdrecb.TIF'); % Se puede elegir la imagen que queramos
oldR = double (oldR); 
%oldR = rot90(oldR); % SOLO para las rw_i es necesaria una rotacion previa

% _____________________________________________________________________

% Podemos incluir un preprocesado de filtrado de ruido

oldR=wiener2(oldR); % En general suaviza la imagen
% oldR=abs(oldR-50); % Podemos quitar una cantidad fija para eliminar la luz
                     % ambiente en las imagenes experimentales.

% ______________________________________________________________________

% Metemos condiciones de paridad (Imponerlas es una buena eleccion, en general,
% puesto que las imagenes de las que partimos no estan absolutamente horizontales
% y al aplicar paridad arriba-abajo estamos falseando datos.
tam=size(oldR,1);
tam2=size(oldR',1);
temp=(oldR(1:floor(tam/2),:)+flipud(oldR((ceil(tam/2)+1):tam,:)))/2;
oldR(1:floor(tam/2),:)=temp;
oldR((ceil(tam/2)+1):tam,:)=flipud(temp);

% _______________________________________________________________________

% Algoritmo para igualar las intensidades de cada columna de la imagen
% puesto que en ocasiones tenemos unas franjas mas oscuras en la imagen
% experimental.
intensidadnorm=sum(oldR(:,1));
for(i=1:tam2),
   columna=oldR(:,i);
   intensidad=sum(columna);
   columna=columna*intensidadnorm/intensidad;
   oldR(:,i)=columna;
end

% ____________________________________________________________________

% Las imagenes experimentales tienen un cierto escalado en cada columna
% que podemos corregir (tampoco lo hemos utilizado a la hora de obtener
% los resultados finales)

% Definimos los parámetros experimentales
a=831;
z=426;
l=646;
f=-256.71;

xfactor=0:0.001:1;
yfactor=-a*(-1+f*tan(xfactor*pi/2)/z)./(l*(1-tan(xfactor*pi/2).*(f/z+f/l)));
%yfactor=1./yfactor; 
%figure
%plot(xfactor,yfactor);

% Es el momento ahora de cambiar la escala
for(i=1:tam2),
   p=(i-1)/tam2;
   factor=1;
   %factor=-a*(-1+f*tan(p*pi/2)/z)/(l*(1-tan(p*pi/2)*(f/z+f/l))); 
   %factor=1/factor; % Por si en vez de ampliar queremos disminuir
   columna=oldR(:,i);
   columna=columna';
   x=(1:tam);
   y=columna;
   xi=1:1/factor:tam;
   yi=interp1(x,y,xi,'spline');
   tam=size(y',1);
   tami=size(yi',1);
   if (factor<1), % Caso de que el factor sea menor que 1
      if (rem(tami,2)==0),
          % caso par
          dimension=ceil(tam/2)-ceil(tami/2)+1:ceil(tam/2)+ceil(tami/2);
          else
          % caso impar
          dimension=ceil(tam/2)-ceil(tami/2)+1:ceil(tam/2)+ceil(tami/2)-1;
      end
      columna=zeros(1,tam);
      columna(dimension)=yi;
   else % Si el caso es mayor que 1
      columna=yi(ceil(tami/2)-ceil(tam/2)+1:ceil(tami/2)+ceil(tam/2));
   end
   oldR(:,i)=columna';
end

% ____________________________________________________________________

% Reflejamos la imagen para obtener proyecciones de 90 a 180 grados
oldR(:,tam2:(tam2*2-1))= fliplr(oldR);
R = oldR;
%R = R.^2; % Se aplica si la imagen no esta tomada en intensidades
           % cosa que solo ocurre si estamos tratando con imagenes
           % obtenidas de manera teorica.

% __________________________________________________________________

% Una vez preparada la imagen hacemos la reconstuccion final
           
pasito=180/(tam2*2-2); % Define el paso de la transformada
cosita=0:pasito:180; % El vector de angulos

I = iradon(R,cosita); % Realiza una reconstruccion rapida
%I = iradon(R,cosita,'spline','Hamming'); % Realiza una reconstruccion precisa

% Mostramos los resultados con diferentes contrastes
figure
imagesc(R); colormap gray
figure;
imagesc(I); colormap gray
figure;
imagesc(abs(I)); colormap gray
figure;
imagesc(abs(I).^0.5); colormap gray
figure;
imagesc(abs(I).^0.1); colormap gray