%--------------------------------------------------------------------------
% Program     : Rugosity by boxcounting (CalculaRugosidade.m)

% Written by  : Maysa Macedo

% email       : marcelo_at_cbpf.br ; maysagm_at_cbpf.br, dario_at_cbpf.br;
% albuquer_at_cbpf.br, luizlima_at_metalmat.ufrj.br

% Copyright(©): Centro Brasileiro de Pesquisas Fisicas - CBPF/MCT/Rio de Janeiro/Brazil

% Version     : 1.0

% Objective   : Rugosity Computing

% Compiler    : Matlab 7.1- R4-SP3

% Project     : -

% Authors     : Maysa Macedo, Dario oliveira, Marcelo Portes de Albquerque,
% Marcio Portes de Albuquerque, Luiz Lima

% Mode        : In Matlab 7.1- R4-SP3 prompt, at the current directory, you should execute the comand:
% main;

% Comments    : This program is distributed in the hope that it will be useful,
%                but WITHOUT ANY WARRANTY; without even the implied warranty of
%                MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. It is
%                FREEWARE.  
function [EstruturaRugosidade] = CalculaRugosidade(DadosImagem,tipoAlgoritmo)

warning off;
 EstruturaRugosidade = struct( 'inclinacao',	{0},...
                               'original',{0},...
                               'variancia',{0},...
                               'tamanhocaixa',{0},...
                               'dimensaoFractal',{0});
                    
if (tipoAlgoritmo == 1 )
       EstruturaRugosidade=BoxCounting2n(DadosImagem);
elseif (tipoAlgoritmo == 2)
       EstruturaRugosidade=BoxCounting2n4camadas(DadosImagem);
end  

%-------------------------------------------------------------------------

function [EstruturaRugosidade] = BoxCounting2n4camadas(DadosImagem)

warning off;
 EstruturaRugosidade = struct( 'inclinacao',	{0},...
                               'original',{0},...
                               'variancia',{0},...
                               'tamanhocaixa',{0},...
                               'dimensaoFractal',{0},...
                               'tamanhocaixaFit',{0},...
                               'varianciaFit',{0});

Imagem      = DadosImagem.imagem;
EstruturaRugosidade.original=Imagem;

resolucao   = str2num(DadosImagem.imageresolution); %Resolucao da imagem (numero de pixels por linha)
tamanho     = DadosImagem.imagesize; %Tamanho da imagem em microns ou nanometro
tipo        = strtrim(DadosImagem.zunit);

switch lower(tipo)
    
    case {'µm'}
         tamImagem   = str2num(tamanho)*10^-6;
    case {'nm'}
         tamImagem   = str2num(tamanho)*10^-9;
    case {'m'}
         tamImagem   = str2num(tamanho)*10^-12;
    otherwise
      disp('Erro! Tamanho de imagem desconhecido.')
      close
end

Soma=0;
SomaN=0;
MediaN=0;
n_interval=0;
MatNcaixa   = (1:floor(resolucao));
MatNcaixaFit(1:floor(resolucao)) = 0;
TotalVarianciaFit(1:floor(resolucao/2))=0 ;
TotalVariancia(1:floor(resolucao))=0 ;
maximo_divisao = floor(log2(resolucao));
%h = waitbar(0,'Calculando...');
for n=0:maximo_divisao
      %Define o lado do quadrado (2^n)
      lado = resolucao/2^n;
      if (lado < 2)
            n = n-1;
            break;
      end    
      intResoluLado = floor(resolucao/lado);
      MatVarianciaN = [];
      MatVarteste = [];
      DifLado = ((resolucao/2^n)+1)-(resolucao/(2^(n+1))-1);
      Incremento = round(DifLado/4);
      nCaixa = floor(DifLado/Incremento);
      MatCaixa = [];
      for i=0:(intResoluLado-1)
            for j=0:(intResoluLado-1)
                        minimo = floor((lado*j)+1);
                        centroTemp = floor(((lado*(j+1))-((lado*j)+1))/2);
                        
                        centro = minimo + centroTemp;
           
                        variancia = std2(Imagem((lado*i)+1:lado*(i+1),(lado*j)+1:lado*(j+1)));
                        MatVarianciaN(i*intResoluLado+j+1) = variancia;
                        if ~(maximo_divisao == n)
                            IncrementoTemp = Incremento;
                            
                            tamQuad = floor(IncrementoTemp/2);
                            
                            for(nNivel=1:nCaixa)
                                if(DifLado<=nNivel || (abs(centro-tamQuad)<2))
                                   break;
                                else
                                    %centro
                                    %tamQuad
                                   Caixa = std2(Imagem(centro-tamQuad:centro+tamQuad,centro-tamQuad:centro+tamQuad));
                                   MatCaixa((i*intResoluLado+j+1),nNivel) = Caixa;
                                   indice=i*intResoluLado+j+1;
                                                                   
                                   IncrementoTemp = IncrementoTemp + Incremento;
                                   
                                   tamQuad = floor(IncrementoTemp/2); 
                             
                                end
                            end
                        end 
            end 
            
      end
      %MatNcaixa(n)                = lado;
      TotalDesvioPadrao(floor(lado)) = std(MatVarianciaN); 
      TotalVariancia(floor(lado))    = mean(MatVarianciaN);
       
      TotalVarianciaFit(floor(lado)) = log(mean(MatVarianciaN));
      MatNcaixaFit(floor(lado))   = log((floor(lado)*tamImagem)/resolucao);
      
      for(nNivel=1:nCaixa)
           if ((~(maximo_divisao == n)) && (DifLado>nNivel)&& (mean(MatCaixa(:,nNivel))~=0))
             if floor((Incremento*nNivel))==1
                 d=9;
             end
             TotalVariancia(floor((Incremento*nNivel))) = mean(MatCaixa(:,nNivel));
             TotalVarianciaFit(floor((Incremento*nNivel))) = log(mean(MatCaixa(:,nNivel)));
               
             MatNcaixaFit(floor((Incremento*nNivel)))   = log((floor((Incremento*nNivel))*tamImagem)/resolucao);
           
           end
      end
     
      %  Soma= Soma + (log(TotalVariancia(floor(lado)))/log((lado*tamImagem)/resolucao));
       if (n==3)
         pontoMaximoY = log(TotalVariancia(floor(lado)));
         pontoMaximoX = log((lado*tamImagem)/resolucao);
       end
       if (n==7)
            pontoMinimoY = log(TotalVariancia(floor(lado)));
            pontoMinimoX = log((lado*tamImagem)/resolucao);
        end
      %waitbar(n/maximo_divisao);
  end 
  %close(h); 

MatLogVariancia  = log(TotalVariancia);
%converte o tamanho da caixa para microns
MatNcaixaConvert = (MatNcaixa.*tamImagem)./resolucao;
MatLogNcaixa     = log(MatNcaixaConvert);


Inclinacao = (pontoMaximoY-pontoMinimoY)/(pontoMaximoX-pontoMinimoX);
EstruturaRugosidade.varianciaFit = TotalVarianciaFit;
EstruturaRugosidade.tamanhocaixaFit = MatNcaixaFit;
EstruturaRugosidade.inclinacao = Inclinacao;
EstruturaRugosidade.variancia = MatLogVariancia;
EstruturaRugosidade.tamanhocaixa = MatLogNcaixa;
EstruturaRugosidade.dimensaoFractal = 3-Inclinacao;

%--------------------------------------------------------------------------

function [EstruturaRugosidade] = BoxCounting2n(DadosImagem)

warning off;
 EstruturaRugosidade = struct( 'inclinacao',{0},...
                               'original',{0},...
                               'variancia',{0},...
                               'tamanhocaixa',{0},...
                               'dimensaoFractal',{0},...
                               'tamanhocaixaFit',{0},...
                               'varianciaFit',{0});
                    
                             
%Le a imagem de entrada , considera-se que a imagem  sempre vai ser
%quadrada
%DadosImagem = LerDoc(strNomeArqImagem);
Imagem      = DadosImagem.imagem;
EstruturaRugosidade.original=Imagem;

resolucao   = str2num(DadosImagem.imageresolution); %Resolucao da imagem (numero de pixels por linha)
tamanho     = DadosImagem.imagesize; %Tamanho da imagem em microns ou nanometro
tipo        = strtrim(DadosImagem.zunit);

switch lower(tipo)
    
    case {'µm'}
         tamImagem   = str2num(tamanho)*10^-6;
    case {'nm'}
         tamImagem   = str2num(tamanho)*10^-9;
    case {'m'}
         tamImagem   = str2num(tamanho)*10^-12;
    otherwise
      disp('Erro! Tamanho de imagem desconhecido.')
      close
end

Soma=0;
SomaN=0;
MediaN=0;
n_interval=0;
MatNcaixa   = (1:floor(resolucao/2));
MatNcaixaFit(1:floor(resolucao/2)) = 0;
TotalVarianciaFit(1:floor(resolucao/2))=0 ;
TotalVariancia(1:floor(resolucao/2))=0 ;
maximo_divisao = floor(log2(resolucao));
%h = waitbar(0,'Calculando...');
for n=1:maximo_divisao
      %Define o lado do quadrado (2^n)
      lado = resolucao/2^n;
      if (lado < 2)
            n = n-1;
            break;
      end    
      intResoluLado = floor(resolucao/lado);
      MatVarianciaN = [];
      
      for i=0:(intResoluLado-1)
            for j=0:(intResoluLado-1)
                        variancia = std2(Imagem((lado*i)+1:lado*(i+1),(lado*j)+1:lado*(j+1)));
                        MatVarianciaN(i*intResoluLado+j+1) = variancia;
                      
                        
                 
            end 
            
      end
      %MatNcaixa(n)                = lado;
      TotalDesvioPadrao(floor(lado)) = std(MatVarianciaN); 
      TotalVariancia(floor(lado))    = mean(MatVarianciaN);
      TotalVarianciaFit(floor(lado)) = log(mean(MatVarianciaN));
      MatNcaixaFit(floor(lado))   = log((floor(lado)*tamImagem)/resolucao);
      
       if (n==3)
         pontoMaximoY = log(TotalVariancia(floor(lado)));
         pontoMaximoX = log((lado*tamImagem)/resolucao);
       end
       if (n==7)
            pontoMinimoY = log(TotalVariancia(floor(lado)));
            pontoMinimoX = log((lado*tamImagem)/resolucao);
        end
     
  end 
  

MatLogVariancia  = log(TotalVariancia);
%converte o tamanho da caixa para microns
MatNcaixaConvert = (MatNcaixa.*tamImagem)./resolucao;
MatLogNcaixa     = log(MatNcaixaConvert);
Inclinacao = (pontoMaximoY-pontoMinimoY)/(pontoMaximoX-pontoMinimoX);

EstruturaRugosidade.varianciaFit = TotalVarianciaFit;
EstruturaRugosidade.tamanhocaixaFit = MatNcaixaFit;
EstruturaRugosidade.inclinacao = Inclinacao;
EstruturaRugosidade.variancia = MatLogVariancia;
EstruturaRugosidade.tamanhocaixa = MatLogNcaixa;
EstruturaRugosidade.dimensaoFractal = 3-Inclinacao;

                        
%--------------------------------------------------------------------------

                           
