How can I calculate the curvature of an extracted contour by opencv?

Question

I did use the findcontours() method to extract contour from the image, but I have no idea how to calculate the curvature from a set of contour points. Can somebody help me? Thank you very much!

Answer 1

While the theory behind Gombat's answer is correct, there are some errors in the code as well as in the formulae (the denominator t+n-x should be t+n-t). I have made several changes:

• use symmetric derivatives to get more precise locations of curvature maxima
• allow to use a step size for derivative calculation (can be used to reduce noise from noisy contours)
• works with closed contours

Fixes: * return infinity as curvature if denominator is 0 (not 0) * added square calculation in denominator * correct checking for 0 divisor

std::vector<double> getCurvature(std::vector<cv::Point> const& vecContourPoints, int step)
{
  std::vector< double > vecCurvature( vecContourPoints.size() );

  if (vecContourPoints.size() < step)
    return vecCurvature;

  auto frontToBack = vecContourPoints.front() - vecContourPoints.back();
  std::cout << CONTENT_OF(frontToBack) << std::endl;
  bool isClosed = ((int)std::max(std::abs(frontToBack.x), std::abs(frontToBack.y))) <= 1;

  cv::Point2f pplus, pminus;
  cv::Point2f f1stDerivative, f2ndDerivative;
  for (int i = 0; i < vecContourPoints.size(); i++ )
  {
      const cv::Point2f& pos = vecContourPoints[i];

      int maxStep = step;
      if (!isClosed)
        {
          maxStep = std::min(std::min(step, i), (int)vecContourPoints.size()-1-i);
          if (maxStep == 0)
            {
              vecCurvature[i] = std::numeric_limits<double>::infinity();
              continue;
            }
        }


      int iminus = i-maxStep;
      int iplus = i+maxStep;
      pminus = vecContourPoints[iminus < 0 ? iminus + vecContourPoints.size() : iminus];
      pplus = vecContourPoints[iplus > vecContourPoints.size() ? iplus - vecContourPoints.size() : iplus];


      f1stDerivative.x =   (pplus.x -        pminus.x) / (iplus-iminus);
      f1stDerivative.y =   (pplus.y -        pminus.y) / (iplus-iminus);
      f2ndDerivative.x = (pplus.x - 2*pos.x + pminus.x) / ((iplus-iminus)/2*(iplus-iminus)/2);
      f2ndDerivative.y = (pplus.y - 2*pos.y + pminus.y) / ((iplus-iminus)/2*(iplus-iminus)/2);

      double curvature2D;
      double divisor = f1stDerivative.x*f1stDerivative.x + f1stDerivative.y*f1stDerivative.y;
      if ( std::abs(divisor) > 10e-8 )
        {
          curvature2D =  std::abs(f2ndDerivative.y*f1stDerivative.x - f2ndDerivative.x*f1stDerivative.y) /
                pow(divisor, 3.0/2.0 )  ;
        }
      else
        {
          curvature2D = std::numeric_limits<double>::infinity();
        }

      vecCurvature[i] = curvature2D;


  }
  return vecCurvature;
}

Answer 2

For me curvature is:

where t is the position inside the contour and x(t) resp. y(t) return the related x resp. y value. See here.

So, according to my definition of curvature, one can implement it this way:

std::vector< float > vecCurvature( vecContourPoints.size() );

cv::Point2f posOld, posOlder;
cv::Point2f f1stDerivative, f2ndDerivative;   
for (size_t i = 0; i < vecContourPoints.size(); i++ )
{
    const cv::Point2f& pos = vecContourPoints[i];

    if ( i == 0 ){ posOld = posOlder = pos; }

    f1stDerivative.x =   pos.x -        posOld.x;
    f1stDerivative.y =   pos.y -        posOld.y;
    f2ndDerivative.x = - pos.x + 2.0f * posOld.x - posOlder.x;
    f2ndDerivative.y = - pos.y + 2.0f * posOld.y - posOlder.y;

    float curvature2D = 0.0f;
    if ( std::abs(f2ndDerivative.x) > 10e-4 && std::abs(f2ndDerivative.y) > 10e-4 )
    {
        curvature2D = sqrt( std::abs( 
            pow( f2ndDerivative.y*f1stDerivative.x - f2ndDerivative.x*f1stDerivative.y, 2.0f ) / 
            pow( f2ndDerivative.x + f2ndDerivative.y, 3.0 ) ) );
    }
    
    vecCurvature[i] = curvature2D;
    
    posOlder = posOld;
    posOld = pos;
}

It works on non-closed pointlists as well. For closed contours, you may would like to change the boundary behavior (for the first iterations).

UPDATE:

Explanation for the derivatives:

A derivative for a continuous 1 dimensional function f(t) is:

But we are in a discrete space and have two discrete functions f_x(t) and f_y(t) where the smallest step for t is one.

The second derivative is the derivative of the first derivative:

Using the approximation of the first derivative, it yields to:

There are other approximations for the derivatives, if you google it, you will find a lot.