function [centro,titaXY] = ejes(B)
%EJES Obtiene los ejes de inercia de un objeto en una imagen binaria
%   [centro,titaXY] = ejes(B)
%   Devuelve el centro de gravedad centro(cx,cy) y el angulo
%   titaXY (en grados)del eje inercia mom. min con VERTICAL(eje X) 
%   positivo en sentido antihorario.
%   -------->Y
%   |
%   |
%   |
%   v  X
%   La imagen B debe contener solo un objeto y este debe estar en blanco.

% RECIBE: Imagen binaria con un solo objeto                            %
% DEVUELVE: Ang(en grados)del eje inercia mom. min con VERTICAL(eje X) %
%        en sentido antihorario y coordenadas del centro de gravedad   %
%        centroXY en coord (X,Y) filas y columnas                      %
% ALGORITMO EMPLEADO: Algoritmo página 86.                             %
% VARIABLES:   B = imagen introducida. Binaria                         %
%              filas = numero de filas de la imagen B                  %
%              columnas = numero de columnas de la imagen B            %
%              i = variable contador para recorrer las filas de la img.% 
%              j = variable contador para recorrer las columnas        %
%              ci = fila del centro de gravedad                        %
%              cj = columna del centro de gravedad                     %
%              Iii = sumatorio (i - ci)^2                              %
%              Ijj = sumatorio (j - cj)^2                              %
%              Iij = sumatorio (i - ci)*(j-cj)                         %
% BUGS: Solo un objeto en la imagen, sino da orientación del conjunto  %
% SUPOSICIONES : fondo negro y objeto en blanco. Hay 1 objeto en img   %
% LIBRO: VISION ARTIFICIAL INDUSTRIAL. UNIVERSIDAD DE VALLADOLID 2012  %   
% TEMA 4 IMAGENES BINARIAS                                             %
% AUTORES: Eusebio de la Fuente y Félix Miguel Trespaderne             %
%______________________________________________________________________%
%  Ejemplo de ejecucion:                                               %
%----------------------------------------------------------------------%
% >> Ient = imread('llavestr1.bmp');                                            
% >> umbralOpt=graythresh(Ient);  %busca umbral optimo
% >> BW = im2bw(Ient,umbralOpt);  %binariza
% >> BW = ~BW;
% >> [centro,tita] = ejes(BW);  
%______________________________________________________________________%
% >> Ient = imread('llavestr1.bmp');                                            
% >> umbralOpt=graythresh(Ient);  %busca umbral optimo
% >> BW = im2bw(Ient,umbralOpt);  %binariza
% >> BW = ~BW;
% >> s  = regionprops(BW, 'orientation'); %ojo! aquí ang resp horizontal
%______________________________________________________________________%

[filas,columnas] = size(B);   
Iii = 0;
Ijj = 0;
Iij = 0;
area=0;

%halla el centro de gravedad llamando a nuestra funcion centroide
centro= centroide(B);

for i = 1:filas
    for j = 1: columnas
        if B(i,j) ~= 0
            Iii = Iii + (i - centro(1))^2;
            Ijj = Ijj + (j - centro(2))^2;
            Iij = Iij + (i - centro(1))*(j - centro(2));
        end
    end
end

%t es el angulo del eje de inercia expresado en la forma polar (ang que 
%forma con eje x la perpendicular al eje de inercia)
t = 0.5 * atan2( 2*Iij, Iii-Ijj );
%Miramos si realmente es t o su perpendicular tp=t+pi/2 el eje de mom.minimo
tp=t+pi/2;
% Sumas S y Sp de las distancias de los puntos al eje y a su perpendicular
S = Iii*cos(t)*cos(t) +Ijj*sin(t)*sin(t) +  2*Iij*sin(t)*cos(t);
Sp = Iii*cos(tp)*cos(tp) +Ijj*sin(tp)*sin(tp) +  2*Iij*sin(tp)*cos(tp);
if S <= Sp  
    titaEnGrados = t*180/pi;
else
    titaEnGrados = tp*180/pi;
end
% titaEnGrados está referido a forma polar. El ang con eje x (vertical) 
% sera titaXY (positivo en sentido antihorario) 
titaXY = titaEnGrados - 90;  



