%Esta funcao calcula a menor distancia entre pontos de de n curvas
function [EstruturaAjuste] =Fitting(Grafico,nretas)

warning off;
 EstruturaAjuste = struct( 'linhaFit',{0},...
                           'ptSaturacao',{0},...
                           'minimo',{0},...
                           'maximo',{0});
                           
                         
k=nretas; %define numero de retas que se deseja encontrar
Coordj=[];
tamgrafico=size(Grafico);
tamAnterior=0;
for i=1:tamgrafico(2);
   [CoordjTmpx CoordjTmpy CoordjTmpVal] = find(Grafico(i).tamanhocaixaFit); %eixo x
   [CoordiTmpx CoordiTmpy CoordiTmpVal] = find(Grafico(i).varianciaFit); %eixo x
   t=size(CoordjTmpy);
   tamAtual=i*t(2);
   Coordj(tamAnterior+1:tamAtual) = CoordjTmpVal;
   Coordi(tamAnterior+1:tamAtual) = CoordiTmpVal;
   tamAnterior = tamAtual;
  
end

numPontos = numel(Coordi) ;
conjunto =0;
tam = max(abs(max(Coordi)-(min(Coordi))));
tamcaixa = max(abs(max(Coordj)-(min(Coordj))));
lim_d = sqrt((tamcaixa^2)+(tam^2)); %calcula diagonal (valor maximo de rho)
EstruturaAjuste.minimo = min(Coordj);
EstruturaAjuste.maximo = max(Coordj);

%--------------------------
%gera retas iniciais 
%--------------------------
alocaPontos(find(Coordj>mean(Coordj)))=1;
alocaPontos(find(Coordj<mean(Coordj)))=2;



% figure;
% hold on;
% title('Primeiro Chute');
% 
% x=(0.1:0.1:tamcaixa);
% cores=['r' 'g' 'b' 'c' 'm' 'y' 'k' 'w']; %cria lista de cores pra plotar
% %plota pontos
% for (j=1:numPontos)
%      plot(Coordj(j),Coordi(j),strcat('+',cores(alocaPontos(j))));
%   
% end
% hold off

 
%---------------------%
%Line Fitting
%---------------------%
contador=0;
termina=0;

%iniciliza o comparador de iteracoes com o primeiro chute
for (m=1:k)
    comparador(1,m)=sum(find(alocaPontos==m));
end
alocaPontosFit=alocaPontos;

%faz refit das retas enquanto houver mudança de alocacao de pontos entre
%uma iteracao e outra
while (termina==0);

    contador=contador+1;
    
    for (i=1:k)
        %pega populacao de pontos alocados pra determinada reta
        x=Coordj(find(alocaPontosFit==i));
        y=Coordi(find(alocaPontosFit==i));
        
        %miniza a funcao de distancia quadratica utilizando o metodo de
        %total least squares, utilizando os autovetores como demonstrado em
        %sala de aula
        x2=x.^2;
        y2=y.^2;        
        xy=x.*y;
        mX2=mean(x2);
        mY2=mean(y2); 
        mXY=mean(xy);
        mX=mean(x);
        mY=mean(y);
        A=[mX2-mX*mX, mXY-mX*mY; mXY-mX*mY, mY2-mY*mY];
        [vec val]=eig(A);

        escala = 1/(vec(1,1)^2+vec(2,1)^2); %acha escala para (a^2+b^2=1)
        
        %verifica qual dos dois autovetores minimiza a funcao distancia
        %opcao 1
        opcao1A=escala*vec(1,1);
        opcao1B=escala*vec(2,1);
        opcao1C=-opcao1A*mX-opcao1B*mY;        
        %opcao 2
        opcao2A=escala*vec(1,1);
        opcao2B=escala*vec(2,1);
        opcao2C=-opcao2A*mX-opcao2B*mY;
        
        %calcula distancia dos pontos a cada opcao de reta
        for (j=1:numPontos)
            distancia1(j)=((opcao1A*Coordj(j)+opcao1B*Coordi(j)+opcao1C)^2);
            distancia2(j)=((opcao2A*Coordj(j)+opcao2B*Coordi(j)+opcao2C)^2);
        end  
        
        %escolhe a menor distancia media
        if (mean(distancia1)>mean(distancia2))
            valorA=opcao2A;
            valorB=opcao2B;
            valorC=opcao2C;
        else
            valorA=opcao1A;
            valorB=opcao1B;
            valorC=opcao1C;            
        end
        
        linhaFit(i,1)=valorA;
        linhaFit(i,2)=valorB;        
        linhaFit(i,3)=valorC;        
        %clear x,y,distancia1,distancia2;
    end
   
    %calcula distancia dos pontos às novas retas
    for (i=1:k)
        for (j=1:numPontos)
            distanciaQuadFit(i,j)=((linhaFit(i,1)*Coordj(j)+linhaFit(i,2)*Coordi(j)+linhaFit(i,3))^2);
        end
    end

    %realoca-se os pontos às retas mais próximas
    for (j=1:numPontos)
        vetor = distanciaQuadFit(:,j);
        [c ix] = min(vetor);
        alocaPontosFit(j)=ix;
    end
    alocaPtSatury = Coordi(find(alocaPontosFit==1));

% figure;    title(strcat('Iteracao-',int2str(contador)));
%     hold on;
%   
%     %x=(min(Coordj):0.1:max(Coordj)); %cria vetor pra ser o eixo x pra plotar as retas
%     %clear x
%     cores=['r' 'g' 'b' 'c' 'm' 'y' 'k' 'w']; %cria lista de cores pra plotar
% 
%     for (i=1:k)
%             fplot(@(x)(-linhaFit(i,1)*x-linhaFit(i,3))/linhaFit(i,2),[min(Coordj) max(Coordj)],cores(i));
%     end
% 
%     for (j=1:numPontos)
%         plot(Coordj(j),Coordi(j),strcat('+',cores(alocaPontosFit(j))));
%     end 
%     

    %verifica se houve mudanca na alocacao dos pontos. Se houve vai para
    %proxima iteracao, se nao houve termina.
    for (m=1:k)
        comparador(contador+1,m)=sum(find(alocaPontosFit==m));
    end
    
    
    %comp2=comparador(contador+1)
    if (comparador(contador,:)==comparador(contador+1,:))
        termina=1;
    end
    EstruturaAjuste.linhaFit = linhaFit;
    tamPtSatur= size(alocaPtSatury);
    if(tamPtSatur(2)<3)
        EstruturaAjuste.ptSaturacao = alocaPtSatury(1);
    else
        EstruturaAjuste.ptSaturacao = alocaPtSatury(end-2);
    end
end


