% *** ---------------------------------------------------------------------------------------- *** %
% ***| Funcion que devuelve la Matriz de la funcion de la Direccion del gradiente            | *** %
% *** ---------------------------------------------------------------------------------------- *** %

function matriz_fd = grad_dir(imag)

[filas,columnas] = size(imag);

[Gx,Gy] = gradient(imag);

G = sqrt(Gx.^2 + Gy.^2) + eps;
Gx_n = Gx./G;                               % Gradientes normalizados
Gy_n = Gy./G;

for i = 1:filas                             
	for j = 1:columnas
       	       
        Dp = [Gy_n(i,j),-Gx_n(i,j)];                    % Vector direccion del pixel p(i,j), perpendicular al gradiente
        
        for dy = -1:1
            if (i+dy) >= 1 & (i+dy) <= filas
                for dx = -1:1
                    if (j+dx) >= 1 & (j+dx) <= columnas
                       if dy ~= 0 | dx ~= 0                             % 8-vecinos del pixel en cuestion
                           
                            Dq = [Gy_n(i+dy,j+dx),-Gx_n(i+dy,j+dx)];    % Vector direccion del pixel q(i+dx,j+dy), perpendicular al gradiente
                            dr = [dy,dx];                               % Vector incremento
                            
                            Lpq = dot(Dp,dr);                           % Vector direccion entre p y q
                            if Lpq >= 0
                                Lpq = dr;
                            else
                                Lpq = -dr;
                            end
                            Lpq_n = Lpq/(sqrt(Lpq(1)^2 + Lpq(2)^2));    % Vector direccion entre p y q normalizado
                            
                            dp = dot(Dp,Lpq_n);
                            dq = dot(Lpq_n,Dq);
                            
                            fd = (1/pi)*(acos(dp) + acos(dq));          % Valor de la direccion del gradiente entre p y q
                            bloque(dy+2,dx+2) = fd; 
                            
                        end
                    end
                end
            end
        end
        
        matriz_fd{i,j} = bloque;
            
    end
end
