close all
%-------------------------------------------------------------------------%
%                   Traitement de l'image                                 %
%                                                                         %
%   Ce premier fichier est le pogramme de traitement. Il lit une image,   %
%   puis appelle differentes fonctions afin de la traiter.                %
clear figure
Im = imread('faisc_31_03_a.jpg'); %Im : variable contenant 3 matrices qui 
                                %contiennent les intensites de rouge, vert
                                %pour chaque pixel de l'image de base

mint = [135 29 69]; %Cette variable correspond a un pixel du rayon ayant la 
                    %plus intensite (prise de maniere completement
                    %empirique

image(Im) %Affichage de l'image de base

imf = faisceau(Im, mint); %appelle de la fonction faisceau, resultat stocke
                          %dans imf (plus d'explication dans le fichier
                          %faisceau.m)

figure %nouvelle fenetre pour une figure
image(imf) %affichage de l'image imf ne selectionnant que le faisceau

[imf2, zx] = affinage(imf); %appelle de la fonction affinage. imf2 est une 
                            %image et zx une matrice 2xp, p etant la 
                            %largeur de l'image (plus d'explication dans le 
                            %fichier affinage.m)
                          

figure %nouvelle fenetre pour une figure
image(imf2) %affichage de l'image imf2

ime = imread('cb_31_03_a.jpg');
[imech, ech] = echelle(ime); %appelle de la fonction echelle. imech est une 
                            %image et ech un flottant (plus d'explication  
                            %dans le fichier affinage.m)

figure %nouvelle fenetre pour une figure
image(imech) %affichage de l'image imech

res = pos_deriv(ech, zx, 0.018, 0.066); %appelle de la fonction pos_deriv. res  
                                  %est une matrice nx5 (plus d'explication   
                                  %dans le fichier affinage.m)

figure  %nouvelle fenetre pour une figure
plot(res(:,1), res(:,2)) %graphe du faisceau

title('Height of the beam as a function of x')
xlabel('x (m)')
ylabel('Height (m)')


figure
plot(res(:,2),res(:,4), '.')
title('Refractive index as a function of height')
xlabel('Height (m)')
ylabel('Refractive index')
%errorbar(res(:,4),res(:,5))

% figure
% ft = fit(res(:,1), res(:,2), 'poly9')
% plot(ft)
% save('fit_31_03_b.txt', 'ft');
save('31_03_a_fit.txt', 'res', '-ascii'); %ecrit les donnes de res dans un
                                          %fichier .txt

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

function [ image_fais ] = faisceau( im, min )
%FAISCEAU Mat(3*Mat(int)) * Mat(int) -> Mat(3*Mat(int))
%   Retourne une image ne selectionnant "theoriquement" que le faisceau du
% laser grace a la coueur rouge et prend en argument l'image de depart et
% le pixel minimal


n= size(im,1); %nombre de ligne
p= size(im,2); %nombre de colonne

im = uint8(im); %convertit l'image pour que l'intensite
                %varie de 0 a 255 au lieu de 0 a 1
                
image_fais= zeros(n,p,3); %initialise l'image en noir complet
image_fais = uint8(image_fais); %convertit l'image pour que l'intensite
                                %varie de 0 a 255 au lieu de 0 a 1

I = 1:n; 
J = 1:p;

for i = I
    
    for j = J
        
         if im(i,j,1)>min(1)-1 %selection de la couleur rouge
            image_fais(i,p+1-j,1)=im(i,j,1); %mise a l'endroit
         end   
         
    end
            
end
end

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

function [ imf, xm ] = affinage( im )
%AFFINAGE Mat(3*Mat(int)) -> Mat(3*Mat(int)) * Mat(int)
% Retourne l'image en moyennant les positions verticales des pixels rouges 
%pour avoir un trait fin et le tableau contenant ces valeurs et
%l'ecart entre la valeur max et la valeur min en valeur absolue et prend en
%argument l'image avec seulement le faisceau

n= size(im,1); %nombre de ligne
p= size(im,2); %nombre de colonne

im = uint8(im); %convertit l'image pour que l'intensite
                %varie de 0 a 255 au lieu de 0 a 1
                
xm=zeros(2,p); %initizalisation d'une matrice 2xp

I = 1:n;
J = 1:p;

%-------------------------------------------------------------------------%
%Remplissage du tableau xm


for j = J
    
    s=0;        %
    x=0;        %initialisation de 3 variables
    sx = 0;     %
    
    for i = I
        
        if im(i,j, 1) > 0 %Un pixel non noir
            
            if x == 0 %verification pour la premiere valeur
                xmin = i; %Recuperation de la position du coefficent 
                          %minimale
            end
            
            x = x+i;        %
            xmax = x-sx;    % Pour le calcul de la posisition maximale qui
            sx=sx+i;        % doit se faire a chaque boucle
            
            s=s+1;          %Somme pour le calcul de la moyenne
            
        end
        
    end
    
    ecart = abs(xmin-xmax)/2; %calcul de l'ecart entre valeur max et min
    
    xm(1,j)=floor(x/s);       %un pixel doit être un entier
    xm(2,j)=ecart;

end

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

imf = zeros(n,p,3); %initialisation de l'image contenant le trait fin du 
                    %faisceau

%-------------------------------------------------------------------------%
% Elimination des valeurs tres aberantes

x=xm(1,p-1);
xm(1,p)=x;  
test = xm(1,p);

for j=2:p
    if abs(xm(1,p+1-j)-test)>8
        xm(1,p+1-j)=NaN; %La valeur NaN permet de faire en sorte que la 
                         %valeur ne sera pas prise en compte lors des
                         %prochains calculs
                         
    else
        test = xm(1,p+1-j);
    end
    
end
%-------------------------------------------------------------------------%

imf = uint8(imf); %convertit l'image pour que l'intensite
                  %varie de 0 a 255 au lieu de 0 a 1
                  
%-------------------------------------------------------------------------%
%Creation du faisceau affine dans une nouvelle image

for j = J
    if xm(j)>0
        imf(xm(1,j), j,1)=255; 
    end
end
end

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

function [ image_ech, ech] = echelle( im )
%ECHELLE Mat(3*Mat(int)) -> Mat(3*Mat(int)) * float 
% Renvoie une image avec seulement le carre temoin en vert ainsi que
% l'echelle correspondant et prend en argument l'image de base

im = uint8(im); %convertit l'image pour que l'intensite
                %varie de 0 a 255 au lieu de 0 a 1

n= size(im,1); %nombre de ligne
p= size(im,2); %nombre de colonne

cbv = 60; %Intensite minimale du vert sur le carre blanc prise 
          %empiriquement

image_ech= zeros(n,p,3); %initialise l'image en noir complet
image_ech = uint8(image_ech); %convertit l'image pour que l'intensite
                              %varie de 0 a 255 au lieu de 0 a 1

tcb = 0.01; %m :taille du carre blanc 

I = 1:n;
J = 1:p;

ech = 0;    % 
s= 0;       %intialisation des variables

%-------------------------------------------------------------------------%
% Meme principe que pour la selecton du faisceau sauf qu'ici le test se
% fait sur le vert de ne pas selectionner le laser

for i = I
    for j = J
         if im(i,j,2)>cbv %selection du vert
            image_ech(i,j,2)=im(i,j,2); 
         end        
    end
            
end
%-------------------------------------------------------------------------%

%-------------------------------------------------------------------------%
% Calcul de l'echelle

cpt = 0;    %
cpti = 1;   %Initialisation des variable

for i = I
    for j = J
        if image_ech(i,j,2)>0 %Pixel qui n'est pas noir
            cpt = cpt + 1; %Compteur augment
        else
            if and(cpt>0, cpt<10) %test pour eviter les pixels aui ont pu 
                                  %rester apres la selction mais qui
                                  %forment une chaine pas assez longue pour
                                  %faire partie du carre
            cpt = 0; %remise a zero du compteur
            end
        end
            
    end
    if cpt>0                %test pour savoir si on a eu une tranche du 
                            %carre
        s=s + cpt;          %somme de tout les cotes pour minimiser les
                            %incertitudes
        cpti = cpti + 1;    %nombre de cote +1
        cpt=0;              %remise du compteur a zero
    end

end

ech = (cpti-1)*tcb/s;   %regle de trois pour calculer l'echelle
end
%-------------------------------------------------------------------------%

function [ tab_donnees ] = pos_deriv( ech, zx, alpha, z0)
%POS_DERIV float * Mat(int) * float * float-> Mat(float)
% Renvoie une matrice de 5 colonnes aves en premiere colonne la position en
% x (prise tout les dix pixels), la deuxieme contient l'altitude, la
% troisieme la derive spatiale, la quatrieme l'indice et en dernier
% l'incertitude sur l'indice. Elle prend en argument l'echelle, la matrice
% contenant l'altitude en pixel et l'ecart valeur max-min, l'angle
% initiale et l'altitude initiale

neau = 1.33; %on suppose qu'initialement vu que le faisceau est en haut de 
            %la cuve est donc que la concentration en sucre est faible,
            %l'indice optique sera celui de l'eau.

q = size(zx,2); %nombre de colonne de zx

z = zx(1,1);

%-------------------------------------------------------------------------%
%Mise a l'echelle de zx en prenant en compte que la premiere valeur devra
%etre z0

for i = 1:q
    zx(1,i) = z0 - ech * (zx(1,i)-z);
    zx(2,i) = zx(2,i)*ech;
end


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

C0 = neau^2/(sin(alpha)^2+1);    %calcul de la constante C0 de la formule 
                                %en tenant compte des valeurs initiales
                                
%-------------------------------------------------------------------------%
%Remplissage du tableau de donnees

j=1;

for i = 1:10:q %10 est le pas afin de ne pas avoir trop de valeur
    if i+10<q %test pour eviter le probleme en fin de matrice
        tab_donnees(j,1)=(i-1)*ech; %Remplissage des x
        tab_donnees(j,2)=zx(1,i); %Remplissage des z
        tab_donnees(j,3)=devfite((i-1)*ech); %derive spatiale du fit
        %Ce fit a ete obtenu a l'aide des donnees contenu dans zx.
        %La ligne pour l'obtenir est : ft = fit(zx(2,:),zx(1,:), 'poly9')
        tab_donnees(j,4)=sqrt(C0*(devfite((i-1)*ech)^2+1)); %indice
        %tab_donnees(j,5)=1/2*(zx2(2,i)/zx2(1,i)); %incertidue de l'indice
        j = j+1;
    end
end

end
%-------------------------------------------------------------------------%