% Exponente de Hurst en anillos de Saturno
% Programa para estudiar exponente de Hurst y dimensión fractal de la estructura de anillos

% Trabajo de Tesis - UNLP
% Horacio Daniel Salomone

% Start

%Image reading
AS1=imread('PIA22418.jpg'); 

%Correct to a circle
image(AS1), axis image;

%Shows the figure
subplot(2,2,1);
imshow(AS1),title('Imágen Original-Cassini - Anillos de Saturno 5');

%Cuts the figure and shows the cut
AS2=imcrop;
subplot(2,2,2);
imshow(AS2),title('Imagen Recortada-Anillos de Saturno 5');

%Loop to get the grayscale profiles

    %image size AS2
    sz=size(AS2);
    n=sz(1,1);
            
    %Intensity profile
    Perfil='Vectores Perfil Intensidad.txt';
    fid=fopen(Perfil,'a+');
    formatSpec='%9.5f\r\n';
      
    %Int vectors
        for i=1:n
        Int=improfile(AS2,[1 1018],[1 (0+i)]);
        %Recorta ceros del vector Int
        Int(isnan(Int))=0;
        %Int(Int==0)=[];
        l=1;
            while Int(l)==0
                l=l+1;
            end
        L=length(Int);
        r=1;
            while Int(L-r)==0
                r=r+1;
            end
        Int(L-r+1:L)=[];    
        Int(1:l)=[];
        end
        
%Hurst exponent calculation
    for i=1:120
    m=floor(L/(2*i));
        for j=1:i
        r=Int(1+(j-1)*m:j*m);
        M=mean(r);
        x=(r-M);
        V=cumsum(x);
        R(j)=max(V)-min(V);
        S(j)=std(r);
        end
    tau(i)=m;
    RS(i)=mean(R./S);
    end
    
%Linear regression to obtain Hurst
    subplot (2,2,3);
    plot(log10(tau),log10(RS),'+')
    xlabel('log(\tau)','FontSize',12)
    ylabel('log(R/S)','FontSize',12)
    hold on
    q=polyfit(log10(tau),log10(RS),1);
    t=0.5:0.1:3;
    y=q(1)*t+(q(2));
    plot(t,y,'r','LineWidth',2)
    %text(1,2,['y =' num2str(q(1)),' x ' num2str(q(2))],'FontSize',12)
    hold off