OpenCV detect partial circle with noise

Question

I have tried to use OpenCV HoughCircles and findContours to detect a circle however the circle isn't complete enough or there is too much noise in the image for these algorithms. Or perhaps we are just not familiar enough with the OpenCV. Attached is my image that I need to find the circle on. You should be able to see it clearly with your eyes however, none of the circle detection algorithms seem to work. I have found that applying a median filter cleans up most of the noise but even after median filtering the algorithms can't detect the circle.

Note I even looked at and tried the solution here so it is not a duplicate of that question: Detect semicircle in OpenCV

Any ideas? This is my source image that I need to use.

Also, the reason I want to detect the circle is I want to do a calculation using only the points that are a part of the circle.

Original Image: http://www.collegemobile.com/IMG_2021.JPG

Median Filtered Image: http://www.collegemobile.com/IMG_2022.JPG

Answer 1

Here you go:

I'm using my 2nd answer from Detect semi-circle in opencv and modify it a little. This version now detects the best found semi-circle (regarding completeness).

But first I want to tell you why the accepted answer of link to Detect semi-circle in opencv stack overflow question does not work here (beside noise): You have only edges of the circle! as stated in that question, HoughCircle function computes the gradient internally, which does not work well for edgy images.

But now how I do it:

using this as input (your own median filtered image (I've just cropped it):

First I "normalize" the image. I just stretch values, that smallest val is 0 and biggest val is 255, leading to this result: (maybe some real contrast enhancement is better)

after that I compute the threshold of that image with some fixed threshold (you might need to edit that and find a way to choose the threshold dynamically! a better contrast enhancement might help there)

from this image, I use some simple RANSAC circle detection(very similar to my answer in the linked semi-circle detection question), giving you this result as a best semi-sircle:

and here's the code:

int main()
{
    //cv::Mat color = cv::imread("../inputData/semi_circle_contrast.png");
    cv::Mat color = cv::imread("../inputData/semi_circle_median.png");
    cv::Mat gray;

    // convert to grayscale
    cv::cvtColor(color, gray, CV_BGR2GRAY);

    // now map brightest pixel to 255 and smalles pixel val to 0. this is for easier finding of threshold
    double min, max;
    cv::minMaxLoc(gray,&min,&max);
    float sub = min;
    float mult = 255.0f/(float)(max-sub);
    cv::Mat normalized = gray - sub;
    normalized = mult * normalized;
    cv::imshow("normalized" , normalized);
    //--------------------------------


    // now compute threshold
    // TODO: this might ne a tricky task if noise differs...
    cv::Mat mask;
    //cv::threshold(input, mask, 0, 255, CV_THRESH_BINARY | CV_THRESH_OTSU);
    cv::threshold(normalized, mask, 100, 255, CV_THRESH_BINARY);



    std::vector<cv::Point2f> edgePositions;
    edgePositions = getPointPositions(mask);

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

    //TODO: maybe seed random variable for real random numbers.

    unsigned int nIterations = 0;

    cv::Point2f bestCircleCenter;
    float bestCircleRadius;
    float bestCirclePercentage = 0;
    float minRadius = 50;   // TODO: ADJUST THIS PARAMETER TO YOUR NEEDS, otherwise smaller circles wont be detected or "small noise circles" will have a high percentage of completion

    //float minCirclePercentage = 0.2f;
    float minCirclePercentage = 0.05f;  // at least 5% of a circle must be present? maybe more...

    int maxNrOfIterations = edgePositions.size();   // TODO: adjust this parameter or include some real ransac criteria with inlier/outlier percentages to decide when to stop

    for(unsigned int its=0; its< maxNrOfIterations; ++its)
    {
        //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);

        // inlier set unused at the moment but could be used to approximate a (more robust) circle from alle inlier
        std::vector<cv::Point2f> inlierSet;

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

        // update best circle information if necessary
        if(cPerc >= bestCirclePercentage)
            if(radius >= minRadius)
        {
            bestCirclePercentage = cPerc;
            bestCircleRadius = radius;
            bestCircleCenter = center;
        }

    }

    // draw if good circle was found
    if(bestCirclePercentage >= minCirclePercentage)
        if(bestCircleRadius >= minRadius);
        cv::circle(color, bestCircleCenter,bestCircleRadius, cv::Scalar(255,255,0),1);


        cv::imshow("output",color);
        cv::imshow("mask",mask);
        cv::waitKey(0);

        return 0;
    }

with these helper functions:

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);
}


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;
}

edit : one more thing: speed performance highly depends on maxNrOfIterations. If that matters you really should read about RANSAC an when to stop it. So you might be able to decide early that a found circle is the right one and dont need to test any other ones ;)

Answer 2

I implemented a least square algorithm to fit a circle to 2d points. I applied the algorithm to the thresholded image from Micka but used a morphological opening to remove the outliers before.

Mat img;
img = imread("iokqh.png");

if (img.empty())
{
    cout << "Could not open image..." << endl;
    return -1;
}

cvtColor(img, img, COLOR_BGR2GRAY);

int dilation_type = 0;
int dilation_elem = 0;

if (dilation_elem == 0) { dilation_type = MORPH_RECT; }
else if (dilation_elem == 1) { dilation_type = MORPH_CROSS; }
else if (dilation_elem == 2) { dilation_type = MORPH_ELLIPSE; }

int size = 1;

Mat element = getStructuringElement(dilation_type, Size(2 * size + 1, 2 * size + 1), Point(size, size));
morphologyEx(img, img, MORPH_OPEN, element);

vector<Point2f> points;
for (int x = 0; x < img.cols; x++)
{
    for (int y = 0; y < img.rows; y++)
    {
        if (img.at<uchar>(y, x) > 0)
        {
            points.push_back(cv::Point2f(x, y));
        }
    }
}

//// Least Square Algorithm 

float xn = 0, xsum = 0;
float yn = 0, ysum = 0;
float n = points.size();

for (int i = 0; i < n; i++)
{
    xsum = xsum + points[i].x;
    ysum = ysum + points[i].y;
}

xn = xsum / n;
yn = ysum / n;

float ui = 0;
float vi = 0;
float suu = 0, suuu = 0;
float svv = 0, svvv = 0;
float suv = 0;
float suvv = 0, svuu = 0;

for (int i = 0; i < n; i++)
{
    ui = points[i].x - xn;
    vi = points[i].y - yn;

    suu = suu + (ui * ui);
    suuu = suuu + (ui * ui * ui);

    svv = svv + (vi * vi);
    svvv = svvv + (vi * vi * vi);

    suv = suv + (ui * vi);

    suvv = suvv + (ui * vi * vi);
    svuu = svuu + (vi * ui * ui);
}

cv::Mat A = (cv::Mat_<float>(2, 2) <<
    suu, suv,
    suv, svv);

cv::Mat B = (cv::Mat_<float>(2, 1) <<
    0.5*(suuu + suvv),
    0.5*(svvv + svuu));

cv::Mat abc;
cv::solve(A, B, abc);

float u = abc.at<float>(0);
float v = abc.at<float>(1);

float x = u + xn;
float y = v + yn;

float alpha = u * u + v * v + ((suu + svv) / n);
float r = sqrt(alpha);

////

cvtColor(img, img, COLOR_GRAY2BGR);

// Draw circle
circle(img, Point(x, y), r, Scalar(255, 0, 0), 1, 8, 0);
imshow("window", img);
waitKey(0);

Here is the result: