algorithmcurvepolynomials

How to detect in real time a "knee/elbow" (maximal curvature) in a curve


In the following curve (blue line) I'm trying to detect the "knee/elbow" which should be located around x = 2.5

enter image description here

This is the set of values I'm using:

x = {-10, -9, -8, -7, -6, -5, -4, -3, -2 , -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}

y = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 107, 122, 145, 176, 215, 262, 317, 380, 451, 530, 617}

I have tried the Kneedle algorithm and the formal definition of the curvature of a graph (signed curvature). The problem I have with the Kneedle algorithm is that in the real time application (embedded system) I don't know which will be the maximal value of the y axis, so I can't properly normalize the points nor find a slope value which works for all the cases. When using the formal definition of the curvature of a graph I try to fit the curve with an order 5 polynomial (green line) and then get the values of the derivatives to calculate the curvature. Nevertheless, this method finds the curvature around x = -2 as there is a curvature around that point because of the polynomial.

enter image description here

Could somebody suggest me a way to detect the knee/elbow?


Solution

  • This is not a continuous curve, nor an approximation of a continuous curve, so we don't have to waste time here: you have a simple polygon. The polygon equivalent to "radius of curvature" is the difference in angle of incidence and angle of refraction: the smaller that difference is, the larger the radius of "curvature" is.

    Provided you sample your data properly, we can compute the angular difference at each data point:

    for (i=1, i<p.length -1):
      vector1 = p[i] - p[i-1] // assuming your language of choice has points as primitives
      vector2 = p[1+1] - p[i] // if not, you'll have to extract x/y separately.
      p[i].angle = findAngle(vector1,vector2)
    

    The findAngle function should be easy enough to implement, there are a million tutorials on how to implement this in your favourite language (it's basic arithemtic in many, and even built-in for some, languages). And that's it, we've turned our 2D data into 3D data with (x,y,z) coordinates in which z represents the local angle in travel across the point.

    To then find any "knee", we can look at all sets of three data points (x-1), (x), and (x+1) and consider those x where the local angular difference is larger than the angular differences of its neighbours. The x with the biggest difference is the "winner": you've found your knee. (or really, a knee, since point data makes zero promises not to go up and down many many times, resulting in a polygon with lots of inflections)

    knee = undefined
    max_diff = 0;
    for (i=1, i<p.length -1):
      a = p[i].angle
    
      a1 = p[i-1].angle
      d1 = a1-a
    
      a3 = p[i+1].angle
      d2 = a3-a
    
      diff = ... // there's a few ways to compute this
      p[i].diff = diff // always useful if we need to do more work later
    
      if (diff > max_diff):
        max_diff = diff
        knee = p[i]
    

    I've left the diff computation up to you, because maybe you just want d1+d2, or (d1+d2)/2, or maybe you want to switch based on whether d1 or d2 (but not both) are 0, etc.

    The important thing to note here, of course, is that we can do all these things "in one go", as we gather data points, since new points don't affect the position of old points, so at some point n, we already know the angles up to n-1 and we already know the knees up to n-2, so we can compute all these values in a single linear pass. Nice and efficient.

    This approach is analogous to finding the root(s) of the derivative of a data-fitting curve, but when all we have is numerical data, sometimes the most correct approach is to work with that data, not with "a presumed reconstruction of the original thing you derived the data from". In this case, you have no idea what your data source is doing between sample points. Maybe it's well behaved, maybe it's not, but you don't know, so it's well worth not wasting cycles on complicating the problem and instead work directly with the properties you know are true (and of course, what they are true for: these are sensor readings and your sensor can absolutely be noisy or even faulty).