% *** ---------------------------------------------------------------------------------------- *** %
% ***| Funcion que crea la subimagen correspondiente y devuelve su Mapa de Coste Local,      | *** %
% ***| asi como las posiciones relativas a la subimagen de los pixeles semilla y posicion    | *** %
% *** ---------------------------------------------------------------------------------------- *** %

function [L,sem_rel,pos_rel,margenes] = mcl(imag,semilla,posicion,dist_margen)

wz = 7/10;                              % Pesos de las tres funciones de coste
wg = 2/10;
wd = 1/10;

im = double(imag);
im_R = double(imag(:,:,1));             % Lectura de las componentes RGB
im_G = double(imag(:,:,2));
im_B = double(imag(:,:,3));

% Creacion de la Sub-Imagen de la imagen y de cada componente RGB

[imag,sem_rel,pos_rel,margenes,ini_sub,fin_sub] = subimagen(im,semilla,posicion,dist_margen);
imag_R = im_R(ini_sub(1):fin_sub(1) , ini_sub(2):fin_sub(2));
imag_G = im_G(ini_sub(1):fin_sub(1) , ini_sub(2):fin_sub(2));
imag_B = im_B(ini_sub(1):fin_sub(1) , ini_sub(2):fin_sub(2));

% __________________________________________________________________________________________________
%                               Matriz de la funcion de 'Zero-crossing' 

fz_R = edge(imag_R,'canny');                % Deteccion de bordes en las 3 componentes
fz_G = edge(imag_G,'canny');
fz_B = edge(imag_B,'canny');

fz_inv = fz_R | fz_G | fz_B;                % Con un '1' en una componente ya tenemos frontera

fz = 1 - double(fz_inv);                    % Invirtiendo la salida de 'edge'

% __________________________________________________________________________________________________
%                           Matriz de la funcion Magnitud del Gradiente 

fg_R = grad_mag(imag_R);                    % Calculo de las 3 Matrices de la Magnitud del Gradiente
fg_G = grad_mag(imag_G);                    % en las 3 componentes
fg_B = grad_mag(imag_B);

fg = min(min(fg_R,fg_G),fg_B);              % Solo se tendra en cuenta el valor minimo entre R,G,B

% __________________________________________________________________________________________________
%                           Matriz de la funcion Direccion del Gradiente 

fd = grad_dir(imag);

% __________________________________________________________________________________________________
%                                      MATRIZ DE COSTE LOCAL

[filas,columnas,dim] = size(imag);     

for i = 1:filas                     % Bucle para la formacion del Mapa de Coste Local (MCL)
    for j = 1:columnas 
        
        L{i,j} = inf*ones(3,3);             % Inicializacion de las celdas de L (MCL)
        
        ini_fila = 1;                       % Indices necesarios en la creacion la matriz del MCL
        fin_fila = 3; 
        ini_columna = 1; 
        fin_columna = 3;  
         
        arriba = i - 1;
        abajo = i + 1;
        izquierda = j - 1;
        derecha = j + 1;

        if i == 1                           % Si nos encontramos en los bordes 1ª y ultima fila
            ini_fila = 3;
            arriba = i + 1;
        elseif i == filas
            fin_fila = 1;
            abajo = i - 1;
        end
     
        if j == 1                           % Si nos encontramos en los bordes de 1ª y ultima columna
            ini_columna = 3;
            izquierda = j + 1;
        elseif j == columnas
            fin_columna = 1;
            derecha = j - 1;
        end
        
        %Formacion de los costes del pixel (i,j) hacia sus D-vecinos
        L{i,j}(ini_fila:2:fin_fila,ini_columna:2:fin_columna) = ...
            wd*fd{i,j}(ini_fila:2:fin_fila,ini_columna:2:fin_columna) + wz*fz(arriba:2:abajo,izquierda:2:derecha) + wg*fg(arriba:2:abajo,izquierda:2:derecha);
        
        %Formacion de los costes del pixel (i,j) hacia sus 4-vecinos
        if i ~= 1              
            L{i,j}(1,2) = wd*fd{i,j}(1,2) + wz*fz(i-1,j) + wg*fg(i-1,j);        % vecino arriba
        end
        if i ~= filas
            L{i,j}(3,2) = wd*fd{i,j}(3,2) + wz*fz(i+1,j) + wg*fg(i+1,j);        % vecino abajo
        end
        if j ~= 1
            L{i,j}(2,1) = wd*fd{i,j}(2,1) + wz*fz(i,j-1) + wg*fg(i,j-1);        % vecino izquierda
        end
        if j ~= columnas
            L{i,j}(2,3) = wd*fd{i,j}(2,3) + wz*fz(i,j+1) + wg*fg(i,j+1);        % vecino derecha
        end
               
    end
end
