% MOMENTS
%
% Function calculates the moments of a binary image and returns 
% the centroid, the angle of axis of minimum inertia, and a measure 
% of 'roundness'.  The function assumes that there is only one object
% in the binary image.
%
% function [centroid, thetamin, roundness] = moments(im)
%
% Argument:  im   - a binary image containing values of 0 or 1
%
% Returns:   centroid  - a 2 element vector
%            thetamin  - the angle of axis of minimum inertia (radians)
%            roundness - ratio of minimum inertia/maximum inertia.
%
% Note that positive x is to the right and positive y is downwards
% thus angles are positive clockwise.
%
% The function also displays the image and overlays the position of
% the centroid and the axis of minimum inertia.

 function [centroid, thetamin, roundness] = moments(im)
 
    [rows,cols] = size(im);
    x = ones(rows,1)*[1:cols];    % Matrix with each pixel set to its x coordinate
    y = [1:rows]'*ones(1,cols);   % Matrix with each pixel set to its y coordinate

    area = sum(sum(im));
    meanx = sum(sum(double(im).*x))/area;
    meany = sum(sum(double(im).*y))/area;
    centroid = [meanx, meany];
    
    % coordinates changed with respect to centre of mass
    x = x - meanx;
    y = y - meany;
    
    a = sum(sum(double(im).*x.^2));
    b = sum(sum(double(im).*x.*y))*2;
    c = sum(sum(double(im).*y.^2));
    
    denom = b^2 + (a-c)^2;
    
    if denom == 0
        % let thetas equal arbitrary angles
        thetamin = 2*pi*rand;
        thetamax = 2*pi*rand;
        roundness = 1;
    else
        sin2thetamin = b/sqrt(denom);       %positive solution
        sin2thetamax = -sin2thetamin;
        cos2thetamin = (a-c)/sqrt(denom);   %positive solution
        cos2thetamax = -cos2thetamin;
    
        thetamin = atan2(sin2thetamin, cos2thetamin)/2;
        thetamax = atan2(sin2thetamax, cos2thetamax)/2;
        Imin = 0.5*(c+a) - 0.5*(a-c)*cos2thetamin - 0.5*b*sin2thetamin;
        Imax = 0.5*(c+a) - 0.5*(a-c)*cos2thetamax - 0.5*b*sin2thetamax;
        roundness = Imin/Imax;
    end
    
    % draw an axis proportional to object size
    % 0.5 takes into acount lines with roundness = 0
    % 5   takes into acount small objects, so axis is still visible.
    rho = sqrt(area)/(roundness + 0.5) + 5 ; 
    [X1,Y1] = pol2cart(thetamin,      rho);
    [X2,Y2] = pol2cart(thetamin + pi, rho);

    imshow(im);
    hold;
    line([X1 + meanx, X2 + meanx],[Y1 + meany, Y2 + meany])
    plot(meanx, meany,'r+')
    hold;