Binary Image Orientation

Question

I'm trying to find the orientation of a binary image (where orientation is defined to be the axis of least moment of inertia, i.e. least second moment of area). I'm using Dr. Horn's book (MIT) on Robot Vision which can be found here as reference.

Using OpenCV, here is my function, where a, b, and c are the second moments of area as found on page 15 of the pdf above (page 60 of the text):

Point3d findCenterAndOrientation(const Mat& src)
{
    Moments m = cv::moments(src, true);
    double cen_x = m.m10/m.m00; //Centers are right
    double cen_y = m.m01/m.m00;

    double a = m.m20-m.m00*cen_x*cen_x;
    double b = 2*m.m11-m.m00*(cen_x*cen_x+cen_y*cen_y);
    double c = m.m02-m.m00*cen_y*cen_y;

    double theta = a==c?0:atan2(b, a-c)/2.0;

    return Point3d(cen_x, cen_y, theta);
}

OpenCV calculates the second moments around the origin (0,0) so I have to use the Parallel Axis Theorem to move the axis to the center of the shape, mr^2.

The center looks right when I call

Point3d p = findCenterAndOrientation(src);
rectangle(src, Point(p.x-1,p.y-1), Point(p.x+1, p.y+1), Scalar(0.25), 1);

But when I try to draw the axis with lowest moment of inertia, using this function, it looks completely wrong:

line(src, (Point(p.x,p.y)-Point(100*cos(p.z), 100*sin(p.z))), (Point(p.x, p.y)+Point(100*cos(p.z), 100*sin(p.z))), Scalar(0.5), 1); 

Here are some examples of input and output:

(I'd expect it to be vertical)

(I'd expect it to be horizontal)

Answer 1

I worked with the orientation sometimes back and coded the following. It returns me the exact orientation of the object. largest_contour is the shape that is detected.

CvMoments moments1,cenmoments1;
           double M00, M01, M10;

           cvMoments(largest_contour,&moments1);
           M00 = cvGetSpatialMoment(&moments1,0,0);
           M10 = cvGetSpatialMoment(&moments1,1,0);
           M01 = cvGetSpatialMoment(&moments1,0,1);
           posX_Yellow = (int)(M10/M00);
           posY_Yellow = (int)(M01/M00);

          double theta = 0.5 * atan(
                    (2 * cvGetCentralMoment(&moments1, 1, 1)) /
                    (cvGetCentralMoment(&moments1, 2, 0) -  cvGetCentralMoment(&moments1, 0, 2)));
                theta = (theta / PI) * 180;

                // fit an ellipse (and draw it)

                if (largest_contour->total >= 6) // can only do an ellipse fit
                                                 // if we have > 6 points
                {
                    CvBox2D box = cvFitEllipse2(largest_contour);
                    if ((box.size.width < imgYellowThresh->width) &&  (box.size.height < imgYellowThresh->height))
                    {

                        cvEllipseBox(imgYellowThresh, box, CV_RGB(255, 255 ,255), 3, 8, 0 );
                    }
                }