3

I have following shape.

enter image description here

It may be rotated by unknown angle. I want to determine its rotation in reference to horizontal axis (so shape above would have rotation equal to 0). Best idea I have come up so far is to determine contours of the shape, find minimum area rectangle and then take its rotation as rotation of shape itself.

Mat mask = imread("path_to_image");
vector<vector<Point>> contours;
vector<Vec4i> hierarchy;
vector<RotatedRect> rotatedRects;

findContours(mask, contours, hierarchy, RetrievalModes::RETR_TREE, ContourApproximationModes::CHAIN_APPROX_SIMPLE);

const auto& largestContour = max_element(contours.begin(), contours.end(),
    [](const auto& e1, const auto& e2) { return e1.size() < e2.size(); });
RotatedRect rotatedRect = minAreaRect(*largestContour);

The problem is that rectangle doesn't border the shape in expected way.

enter image description here

I'm not sure if I can go with that and simply calculate rotation from it anyway, because shape comes from other image processing and I don't know if rectangle would not laid on a different diagonal.

Is there more reliable and better way of finding rotation of this shape?

Edit: Image with shape can be in different scale.

serwus
  • 155
  • 1
  • 14
  • 2
    No experience with with OpenCV but if you have a monochromatic image it should be straightforward to perform singular vector decomposition. Then the angle of rotation of your object corresponds to the first eigenvector. – Konrad Rudolph Dec 20 '19 at 16:42
  • 1
    @KonradRudolph I assume you meant [principal component analysis](https://en.wikipedia.org/wiki/Principal_component_analysis)!? This would probably be the simplest approach here indeed… – Michael Kenzel Dec 20 '19 at 16:46
  • 1
    Thanks for fast reply. I'm new in image processing world and I need to become familiar with these concepts. If I find satisfying solution, I will post it here. – serwus Dec 20 '19 at 16:53
  • Does the speed of the algorithm matters and can the shape be in different scale? I can think of a really stupid solution but it would be slow I guess... – Hawky Dec 20 '19 at 17:22
  • @MichaelKenzel The two are mathematically identical (or rather, PCA is solvable with SVD) but you don't actually need the rotated data, which PCA gives you. – Konrad Rudolph Dec 20 '19 at 18:32
  • how do you define the rotation of that shape? What is the reference? If you compare the rotatedRect of the reference to the rotatedRect your current image, do you get the desired rotation angle? – Micka Dec 20 '19 at 19:31
  • WHY is your sample 0 degrees? For that special shape you could detect the half-circle/semi-circle and compute the contour point furthest away from the circle center. – Micka Dec 20 '19 at 19:36
  • @Hawky For now I'm trying to find working solution, but I believe it will matter. I'm afraid It can be in different scale. – serwus Dec 20 '19 at 19:40
  • or you could try to detect the two straight line segments and derive the middle line between them – Micka Dec 20 '19 at 19:43
  • @Micka This shape is part of bigger image. Rotation I showed in the first picture is rotation correct for that bigger image, so I need to determine its value. I don't think this solution with two straight segments could be easly done for rotated images. But I'm grateful for any ideas, I will try them later on! – serwus Dec 20 '19 at 20:01
  • Is the image size constant? Before and after rotation? – nayab Dec 20 '19 at 20:44

2 Answers2

6

I adapted my answer from here: https://stackoverflow.com/a/23993030/2393191 It gives quite good results:

inline void getCircle(cv::Point2f& p1, cv::Point2f& p2, cv::Point2f& p3, cv::Point2f& center, float& radius)
{
    float x1 = p1.x;
    float x2 = p2.x;
    float x3 = p3.x;

    float y1 = p1.y;
    float y2 = p2.y;
    float y3 = p3.y;

    // PLEASE CHECK FOR TYPOS IN THE FORMULA :)
    center.x = (x1*x1 + y1*y1)*(y2 - y3) + (x2*x2 + y2*y2)*(y3 - y1) + (x3*x3 + y3*y3)*(y1 - y2);
    center.x /= (2 * (x1*(y2 - y3) - y1*(x2 - x3) + x2*y3 - x3*y2));

    center.y = (x1*x1 + y1*y1)*(x3 - x2) + (x2*x2 + y2*y2)*(x1 - x3) + (x3*x3 + y3*y3)*(x2 - x1);
    center.y /= (2 * (x1*(y2 - y3) - y1*(x2 - x3) + x2*y3 - x3*y2));

    radius = sqrt((center.x - x1)*(center.x - x1) + (center.y - y1)*(center.y - y1));
}



std::vector<cv::Point2f> getPointPositions(cv::Mat binaryImage)
{
    std::vector<cv::Point2f> pointPositions;

    for (unsigned int y = 0; y<binaryImage.rows; ++y)
    {
        //unsigned char* rowPtr = binaryImage.ptr<unsigned char>(y);
        for (unsigned int x = 0; x<binaryImage.cols; ++x)
        {
            //if(rowPtr[x] > 0) pointPositions.push_back(cv::Point2i(x,y));
            if (binaryImage.at<unsigned char>(y, x) > 0) pointPositions.push_back(cv::Point2f(x, y));
        }
    }

    return pointPositions;
}


float verifyCircle(cv::Mat dt, cv::Point2f center, float radius, std::vector<cv::Point2f> & inlierSet)
{
    unsigned int counter = 0;
    unsigned int inlier = 0;
    float minInlierDist = 2.0f;
    float maxInlierDistMax = 100.0f;
    float maxInlierDist = radius / 25.0f;
    if (maxInlierDist<minInlierDist) maxInlierDist = minInlierDist;
    if (maxInlierDist>maxInlierDistMax) maxInlierDist = maxInlierDistMax;

    // choose samples along the circle and count inlier percentage
    for (float t = 0; t<2 * 3.14159265359f; t += 0.05f)
    {
        counter++;
        float cX = radius*cos(t) + center.x;
        float cY = radius*sin(t) + center.y;

        if (cX < dt.cols)
            if (cX >= 0)
                if (cY < dt.rows)
                    if (cY >= 0)
                        if (dt.at<float>(cY, cX) < maxInlierDist)
                        {
                            inlier++;
                            inlierSet.push_back(cv::Point2f(cX, cY));
                        }
    }

    return (float)inlier / float(counter);
}

float evaluateCircle(cv::Mat dt, cv::Point2f center, float radius)
{

    float completeDistance = 0.0f;
    int counter = 0;

    float maxDist = 1.0f;   //TODO: this might depend on the size of the circle!

    float minStep = 0.001f;
    // choose samples along the circle and count inlier percentage

    //HERE IS THE TRICK that no minimum/maximum circle is used, the number of generated points along the circle depends on the radius.
    // if this is too slow for you (e.g. too many points created for each circle), increase the step parameter, but only by factor so that it still depends on the radius

    // the parameter step depends on the circle size, otherwise small circles will create more inlier on the circle
    float step = 2 * 3.14159265359f / (6.0f * radius);
    if (step < minStep) step = minStep; // TODO: find a good value here.

    //for(float t =0; t<2*3.14159265359f; t+= 0.05f) // this one which doesnt depend on the radius, is much worse!
    for (float t = 0; t<2 * 3.14159265359f; t += step)
    {
        float cX = radius*cos(t) + center.x;
        float cY = radius*sin(t) + center.y;

        if (cX < dt.cols)
            if (cX >= 0)
                if (cY < dt.rows)
                    if (cY >= 0)
                        if (dt.at<float>(cY, cX) <= maxDist)
                        {
                            completeDistance += dt.at<float>(cY, cX);
                            counter++;
                        }

    }

    return counter;
}




int main(int argc, char* argv[])
{

    cv::Mat input = cv::imread("C:/StackOverflow/Input/rotatedShape1.png", cv::IMREAD_GRAYSCALE);
    std::string outString = "C:/StackOverflow/Output/rotatedShape1.png";

    cv::Mat output;
    cv::cvtColor(input, output, cv::COLOR_GRAY2BGR);

    std::vector<std::vector<cv::Point> > contours;
    cv::findContours(input, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE);

    std::vector<cv::Point> biggestContour;
    double biggestArea = 0;
    for (int i = 0; i < contours.size(); ++i)
    {
        double cArea = cv::contourArea(contours[i]);
        if (cArea > biggestArea)
        {
            biggestArea = cArea;
            biggestContour = contours[i];
        }
    }

    if (biggestContour.size() == 0)
    {
        std::cout << "error: no contour found. Press enter to quit." << std::endl;
        std::cin.get();
        return 0;
    }



    cv::Mat mask = cv::Mat::zeros(input.size(), input.type());
    std::vector < std::vector<cv::Point> > tmp;
    tmp.push_back(biggestContour);
    cv::drawContours(mask, tmp, 0, cv::Scalar::all(255), 1); // contour points in the image

    std::vector<cv::Point2f> circlesList;

    unsigned int numberOfCirclesToDetect = 2;   // TODO: if unknown, you'll have to find some nice criteria to stop finding more (semi-) circles

    for (unsigned int j = 0; j<numberOfCirclesToDetect; ++j)
    {
        std::vector<cv::Point2f> edgePositions;
        //for (int i = 0; i < biggestContour.size(); ++i) edgePositions.push_back(biggestContour[i]);
        edgePositions = getPointPositions(mask);



        std::cout << "number of edge positions: " << edgePositions.size() << std::endl;

        // create distance transform to efficiently evaluate distance to nearest edge
        cv::Mat dt;
        cv::distanceTransform(255 - mask, dt, CV_DIST_L1, 3);



        unsigned int nIterations = 0;

        cv::Point2f bestCircleCenter;
        float bestCircleRadius;
        //float bestCVal = FLT_MAX;
        float bestCVal = -1;

        //float minCircleRadius = 20.0f; // TODO: if you have some knowledge about your image you might be able to adjust the minimum circle radius parameter.
        float minCircleRadius = 0.0f;

        //TODO: implement some more intelligent ransac without fixed number of iterations
        for (unsigned int i = 0; i<2000; ++i)
        {
            //RANSAC: randomly choose 3 point and create a circle:
            //TODO: choose randomly but more intelligent,
            //so that it is more likely to choose three points of a circle.
            //For example if there are many small circles, it is unlikely to randomly choose 3 points of the same circle.
            unsigned int idx1 = rand() % edgePositions.size();
            unsigned int idx2 = rand() % edgePositions.size();
            unsigned int idx3 = rand() % edgePositions.size();

            // we need 3 different samples:
            if (idx1 == idx2) continue;
            if (idx1 == idx3) continue;
            if (idx3 == idx2) continue;

            // create circle from 3 points:
            cv::Point2f center; float radius;
            getCircle(edgePositions[idx1], edgePositions[idx2], edgePositions[idx3], center, radius);

            if (radius < minCircleRadius)continue;


            //verify or falsify the circle by inlier counting:
            //float cPerc = verifyCircle(dt,center,radius, inlierSet);
            float cVal = evaluateCircle(dt, center, radius);

            if (cVal > bestCVal)
            {
                bestCVal = cVal;
                bestCircleRadius = radius;
                bestCircleCenter = center;
            }

            ++nIterations;
        }
        std::cout << "current best circle: " << bestCircleCenter << " with radius: " << bestCircleRadius << " and nInlier " << bestCVal << std::endl;
        cv::circle(output, bestCircleCenter, bestCircleRadius, cv::Scalar(0, 0, 255));

        //TODO: hold and save the detected circle.

        //TODO: instead of overwriting the mask with a drawn circle it might be better to hold and ignore detected circles and dont count new circles which are too close to the old one.
        // in this current version the chosen radius to overwrite the mask is fixed and might remove parts of other circles too!

        // update mask: remove the detected circle!
        cv::circle(mask, bestCircleCenter, bestCircleRadius, 0, 10); // here the thickness is fixed which isnt so nice.

        circlesList.push_back(bestCircleCenter);
    }



    if (circlesList.size() < 2)
    {
        std::cout << "error: not enough circles found. Press enter." << std::endl;
        std::cin.get();
        return 0;
    }

    cv::Point2f centerOfMass = circlesList[0];
    cv::Point2f cogFP = circlesList[1];
    std::cout << cogFP - centerOfMass << std::endl;
    float angle = acos((cogFP - centerOfMass).x / cv::norm(cogFP - centerOfMass)); // scalar product of [1,0] and point
    std::cout << angle * 180 / CV_PI << std::endl;

    cv::line(output, centerOfMass, cogFP, cv::Scalar(0, 255, 0), 1);
    cv::circle(output, centerOfMass, 5, cv::Scalar(0, 0, 255), 1);
    cv::circle(output, cogFP, 3, cv::Scalar(255, 0, 0), 1);


    cv::imwrite(outString, output);

    cv::imshow("input", input);
    cv::imshow("output", output);
    cv::waitKey(0);
    return 0;
}

results:

enter image description here enter image description here enter image description here enter image description here enter image description here

Micka
  • 19,585
  • 4
  • 56
  • 74
  • 1
    I like this solution, because I had to find center of that larger circle as next step. :-) It's getting lost a little though if shape is much smaller than background, but it is not a problem, because I just simply use boundingRect function and then crop area with a shape. Thanks! – serwus Jan 02 '20 at 10:55
5

here's the simple logic of finding the center of gravity and the furthest contour point from it. It has an offset of 6 degrees for that contour, either because of the actual contour shape, or because of a slightly wrong center of gravity.

int main(int argc, char* argv[])
{

    //cv::Mat input = cv::imread("C:/StackOverflow/Input/rotatedShape1.png", cv::IMREAD_GRAYSCALE);
    cv::Mat input = cv::imread("C:/StackOverflow/Input/rotatedShape5.png", cv::IMREAD_GRAYSCALE);
    std::string outString = "C:/StackOverflow/Output/rotatedShape5.png";

    cv::Mat output;
    cv::cvtColor(input, output, cv::COLOR_GRAY2BGR);

    std::vector<std::vector<cv::Point> > contours;
    cv::findContours(input, contours, cv::RETR_EXTERNAL, cv::CHAIN_APPROX_NONE);

    std::vector<cv::Point> biggestContour;
    double biggestArea = 0;
    for (int i = 0; i < contours.size(); ++i)
    {
        double cArea = cv::contourArea(contours[i]);
        if (cArea > biggestArea)
        {
            biggestArea = cArea;
            biggestContour = contours[i];
        }
    }

    if (biggestContour.size() == 0)
    {
        std::cout << "error: no contour found. Press enter to quit." << std::endl;
        std::cin.get();
        return 0;
    }

    cv::Point2f centerOfMass(0,0);
    float invContourSize = 1.0f / biggestContour.size();
    for (int i = 0; i < biggestContour.size(); ++i)
    {
        centerOfMass = centerOfMass + (invContourSize * cv::Point2f(biggestContour[i]));
    }

    float furthestDist = 0;
    cv::Point2f furthestPoint = centerOfMass;
    for (int i = 0; i < biggestContour.size(); ++i)
    {
        float cDist = cv::norm(cv::Point2f(biggestContour[i]) - centerOfMass);
        if (cDist > furthestDist)
        {
            furthestDist = cDist;
            furthestPoint = biggestContour[i];
        }
    }

    // find points with very similar distance
    float maxDifference = 20; // magic number
    std::vector<cv::Point2f> listOfFurthestPoints;
    for (int i = 0; i < biggestContour.size(); ++i)
    {
        float cDist = cv::norm(cv::Point2f(biggestContour[i]) - furthestPoint);
        if (cDist < maxDifference)
        {
            listOfFurthestPoints.push_back( biggestContour[i] );
            // render:
            cv::circle(output, biggestContour[i], 0, cv::Scalar(255, 0, 255), 0);
        }
    }

    cv::Point2f cogFP(0, 0);
    float invListSize = 1.0f / listOfFurthestPoints.size();
    for (int i = 0; i < listOfFurthestPoints.size(); ++i)
    {
        cogFP = cogFP + (invListSize * cv::Point2f(listOfFurthestPoints[i]));
    }

    std::cout << cogFP - centerOfMass << std::endl;
    float angle = acos((cogFP - centerOfMass).x / cv::norm(cogFP - centerOfMass)); // scalar product of [1,0] and point
    std::cout << angle * 180 / CV_PI << std::endl;

    cv::line(output, centerOfMass, cogFP, cv::Scalar(0, 255, 0), 1);
    cv::circle(output, centerOfMass, 5, cv::Scalar(0, 0, 255), 1);
    cv::circle(output, cogFP, 3, cv::Scalar(255, 0, 0), 1);


    cv::imwrite(outString, output);
    cv::imshow("input", input);
    cv::imshow("output", output);
    cv::waitKey(0);
    return 0;
}

this is the ouput for several rotations:

enter image description here

enter image description here

enter image description here

enter image description here

enter image description here

I would love to try the circle method, using RANSAC to find the best 2 circles, but maybe won't have the time...

Another way could be to find the turning points of the smoothed contour.

Micka
  • 19,585
  • 4
  • 56
  • 74