% 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, theta, roundness] = moments(im)
   %
   % Argument:  im   - a binary image containing values of 0 or 1
   %
   % Returns:   centroid  - a 2 element vector
   %            theta     - 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, theta, roundness,figNo] = moments(im,figNo)
   [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);   %   "     "     "    "    "  "   "  y    "

   area = sum(sum(im));
if area == 0
   meanx =0; meany =0;
else
   meanx = sum(sum(double(im) .* x))/area;
   meany = sum(sum(double(im) .* y))/area;
end
   centroid = [meanx meany];
   
   xp=x-(ones(size(im))*meanx);
   yp=y-(ones(size(im))*meany);
 
   a = sum(sum(double(im) .* xp .* xp));
   b = 2*sum(sum(double(im) .* xp .* yp));
   c = sum(sum(double(im) .* yp .* yp));
   
   theta = 0.5* atan2(b,(a-c));
   
if a-c == 0
   roundness = 1;
else
   cos2theta = (a-c)/sqrt(b^2+(a-c)^2);
   sin2theta = b/sqrt(b^2+(a-c)^2);
   Imin = 0.5*(c+a)-0.5*(a-c)*(cos2theta)-0.5*b*(sin2theta); % Positive solutions minimize I
   Imax = 0.5*(c+a)-0.5*(a-c)*(-cos2theta)-0.5*b*(-sin2theta); % Negative solutions maximize I
   
   roundness = Imin/Imax;
end

   t = double(im) .* x *cos(theta) + double(im) .* y *sin(theta);
   rho = meany*cos(theta)-meanx*sin(theta);
   x0 = -rho*sin(theta) + t*cos(theta);
   y0 =  rho*cos(theta) + t*sin(theta);

if nargin == 2
	if figNo               % We have a valid figure number
	    figure(figNo);     % Reuse or create a figure window with this number
            imshow(im);
            line(x0,y0);
	else
	    newWindow=0;       % figNo == 0
    end
else
end