
function [angle, nim] = findang(rows, cols, cr, cc, image)
% Usage: [angle, newimage] = findangle(r, c, crow, ccol, image)
%    r: row vector
%    c: column vector
%    crow: row coordinate of center
%    ccol: column coordinate of center
%    image: object's image
%    angle: angle indicating orientation of object (in degrees)
%    newimag: image with line marked

pa = dot((cols - cc),(cols - cc));
pc = dot((rows - cr),(rows - cr));
pb = 2*dot((cols - cc),(rows - cr));

theta = 0.5*(pi + atan(pb/(pa - pc)));
psi_2 = -2*((pa - pc)*cos(2*theta) + pb*sin(2*theta));

if (psi_2 < 0) 
  theta = -pi/2 + theta
end

angle = theta*180/pi;
m = -cos(theta)/sin(theta);
rho = cc*cos(theta) + cr*sin(theta);	
b = rho/sin(theta);

nim = image;

ys = round(m*[1:255] + b);

i = 1; 
while (i < 256)
  if ((ys(i) > 0) & (ys(i) < 256))
    nim(ys(i), i)=255;
  end
  i = i+1;
end


