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

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

For me curvature is:

$|K| = \sqrt{ \frac{ \left(\frac{d^2y(t)}{dt^2}\frac{dx(t))}{dt}-\frac{d^2x(t)}{dt^2}\frac{dy(t)}{dt} \right )^2 }{ \left( \frac{d^2x(t)}{dt^2}+\frac{d^2y(t)}{dt^2} \right )^{3} } }$

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:

$f'(t)=\lim\limits_{n \rightarrow 0}{\frac{f(t+n)-f(t))}{t+n-x}} = \lim\limits_{n \rightarrow 0}{\frac{f(t+n)-f(t))}{n}}$

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.

$f_x'(t)=\lim\limits_{n \rightarrow 0}{\frac{f_x(t+n)-f_x(t))}{t+n-x}} = \lim\limits_{n \rightarrow 0}{\frac{f_x(t+n)-f_x(t))}{n}} \approx \frac{f_x(t+1)-f_x(t))}{1}$

The second derivative is the derivative of the first derivative:

$f_x''(t)=\lim\limits_{n \rightarrow 0}{\frac{f_x'(t+n)-f_x'(t))}{t+n-x}} = \lim\limits_{n \rightarrow 0}{\frac{f_x'(t+n)-f_x'(t))}{n}} \approx \frac{f_x'(t+1)-f_x'(t))}{1}$

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

$f_x''(t)=(f_x(t+2)-f_x(t+1)) - (f_x(t+1)-f_x(t)) = f_x(t+2) - 2f_x(t+1) + f_x(t)$

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

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

Tags:

C++

Opencv

Related

Recent Posts