1

I'm working on finding the inflection points of an image contour, which for a typical function are the set of points where the double derivative vanishes. Although I seem to have been able to tackle the problem for synthetic data but for actual image the code messes up. Here I present the working of my code with the example of a sine wave.

First I generate a list of coordinates for the sine wave in the format of list of list of a list (This peculiar format is required for the functioning of my code).

import numpy as np
all_unique =[[ [i,10*np.sin(np.pi*i/5)] for i in  range(210)]]

The sine wave is generated for x in range (0-210). Since the wave is horizontal, the x coordinate of the inflection points must be multiples of 5,while the y coordinate should be zero.

To find the inflection points I try to find the curvature at each point. Points with zero curvature shall be labelled points of inflection since curvature is directly proportional to the second derivative of a function. Although accounting for the fact that a contour may not always be a function I try to parameterize the X and Y Coordinates with respect to parameter time via Cubic Spline interpolation. Then the formula for curvature is:Formula

# A function to return separate list of x and y coordinates:
from scipy.interpolate import CubicSpline as CS
def getXY(segment):
    x = list()
    y = list()
    for p in segment:
        x.append(p[0])
        y.append(p[1])
    x = np.array(x)
    y = np.array(y)
    return x,y


functionsX = list()
functionsY = list()

for c in all_unique:
    ux,uy = getXY(c) 
    # Parameter time ranges from 0 to the length of the contour ie the number of pixels
    time = np.array(list(range(len(ux))))
    # Fitting the Cubic Spline
    fx = CS(time,ux)
    fy = CS(time,uy)
    functionsX.append(fx)
    functionsY.append(fy)
criticalX = list()
criticalY = list()
curvature_all = list()
for i in range(len(all_unique)):
    ux,uy = getXY(all_unique[i])
    #generating more number of time steps
    time = np.linspace(0,len(ux),num=300)
    # First derivative of x with respect to time
    xp = functionsX[i](time,1)
    # Second derivative of x with respect to time
    xpp = functionsX[i](time,2)
    # First derivative of y with respect to time
    yp = functionsY[i](time,1)
    # Second derivative of x with respect to time
    ypp = functionsY[i](time,2)
    # velocity
    v = np.sqrt(np.square(xp)+np.square(yp))
    # Calculating curvature for 300 distinct timesteps
    curvature = np.array([(xp[i]*ypp[i]-yp[i]*xpp[i])/np.power(v[i],3) for i in range(len(v))])
    curvature_all.extend(curvature)
    abs_curvature = np.absolute(curvature)
    min_curv = abs_curvature.min()
    threshold = factor*min_curv
    timeSteps = list()
    fact = len(ux)/300
    # Filtering the time steps where curvature changes sign and is below a threshold
    for j in range(len(abs_curvature)-1):
        if(curvature[j]*curvature[j+1]<0 and abs_curvature[j]<threshold):
            timeSteps.append(fact * (j+1))
    timeSteps = np.floor(timeSteps)
    timeSteps = np.unique(timeSteps)
    cx = functionsX[i](timeSteps)
    cy = functionsY[i](timeSteps)
    criticalX.append(cx)
    criticalY.append(cy)

This is implemented as a function called findCritical() in future as:

def findCritical(all_unique,factor):
    functionsX = list()
    functionsY = list()
    for c in all_unique:
        ux,uy = getXY(c)
        time = np.array(list(range(len(ux))))
        fx = CS(time,ux)
        fy = CS(time,uy)
        functionsX.append(fx)
        functionsY.append(fy)
    criticalX = list()
    criticalY = list()
    curvature_all = list()
    for i in range(len(all_unique)):
        ux,uy = getXY(all_unique[i])
        time = np.linspace(0,len(ux),num=300)
        xp = functionsX[i](time,1)
        xpp = functionsX[i](time,2)
        yp = functionsY[i](time,1)
        ypp = functionsY[i](time,2)
        v = np.sqrt(np.square(xp)+np.square(yp))
        curvature = np.array([(xp[i]*ypp[i]-yp[i]*xpp[i])/np.power(v[i],3) for i in range(len(v))])
        curvature_all.extend(curvature)
        min_curv = abs_curvature.min()
        threshold = factor*min_curv
        timeSteps = list()
        fact = len(ux)/300
        for j in range(len(abs_curvature)-1):
            if(curvature[j]*curvature[j+1]<0 and abs_curvature[j]<threshold):
                timeSteps.append(fact * (j+1))
        timeSteps = np.floor(timeSteps)
        timeSteps = np.unique(timeSteps)
        cx = functionsX[i](timeSteps)
        cy = functionsY[i](timeSteps)
        criticalX.append(cx)
        criticalY.append(cy)
    return criticalX,criticalY

This returns a pretty decent list of inflection points.

([array([  5.,  10.,  15.,  20.,  25.,  30.,  35.,  39.,  45.,  50.,  55.,
          60.,  65.,  70.,  74.,  79.,  85.,  90.,  95., 100., 105., 109.,
         114., 119., 124., 130., 135., 140., 144., 149., 154., 159., 164.,
         170., 175., 179., 184., 189., 194., 199., 204.])],
 [array([ 1.22464680e-15, -2.44929360e-15,  3.67394040e-15, -4.89858720e-15,
          6.12323400e-15, -7.34788079e-15,  8.57252759e-15, -5.87785252e+00,
          1.10218212e-14, -1.22464680e-14, -2.20560220e-14, -1.46957616e-14,
          5.14475452e-14, -1.71450552e-14,  5.87785252e+00, -5.87785252e+00,
         -1.47081412e-14, -2.20436424e-14, -1.22588476e-14, -2.44929360e-14,
         -9.80955401e-15, -5.87785252e+00,  5.87785252e+00, -5.87785252e+00,
          5.87785252e+00, -1.02895090e-13,  6.85926004e-14, -3.42901104e-14,
          5.87785252e+00, -5.87785252e+00,  5.87785252e+00, -5.87785252e+00,
          5.87785252e+00,  2.94162824e-14,  7.83897748e-14, -5.87785252e+00,
          5.87785252e+00, -5.87785252e+00,  5.87785252e+00, -5.87785252e+00,
          5.87785252e+00])])

However if I have an image of a sine wave this code snippet misbehaves. Consider a sine wave generated as follows. I make a matrix of zeros and make a skeleton sine wave of amplitude 4 and time period 10. The I extract Contours using OpenCV. Since the contour may repeat points I extract unique points using get_unique() function. The aforementioned Code snippet is used for defining a function findCritical() which returns the inflection points.

sine = np.zeros((210,210))
y0 = 105
for i in range(210):
    sine[y0 + int(4*np.sin(np.pi*i/10))][i] = 1

This gives a sine image as:Sine image generated

Now I use functions of OpenCV namely get_points() to extract contours. Since the function to extract contours tends to repeat points I define extra functions like num_to_str() and get_unique() to get a unique list of contours. The order of points is preserved via storing the points and index in a dictionary as key value pairs, reversing the keys and values so that the key is the index and value is the point and then traversing through the dictionary to extract the points. Here is the code:

import cv2
def get_points(skeleton_image):
    cnts = cv2.findContours(skeleton_image.copy().astype('uint8'), cv2.RETR_LIST, cv2.CHAIN_APPROX_NONE)
    cnts = imutils.grab_contours(cnts)
    return cnts
def num_to_str(contour):
    str_list = list()
    for c in contour:
        str_list.append(str(c[0][0])+','+str(c[0][1]))
    return str_list

def get_unique(contour):
    cont = contour.copy()
    str_list = num_to_str(cont)
    index_dict = dict()
    for i in range(len(str_list)):
        if str_list[i] not in index_dict.keys():
            index_dict[str_list[i]] = i
    new_dict = dict((val,key) for key,val in index_dict.items())
    sub_cont = list()
    keys = list(new_dict.keys())
    keys.sort()
    sub_cont.append(list())
    last_ind = 0
    sub_cont[last_ind].append(new_dict[keys[0]])
    for i in range(len(keys)-1):
        if keys[i+1]-keys[i]!=1:
            last_ind+=1
            sub_cont.append(list())
            sub_cont[last_ind].append(new_dict[keys[i+1]])
        else:
            sub_cont[last_ind].append(new_dict[keys[i+1]])
    num_cont = sub_cont.copy()
    for l in range(len(num_cont)):
        for s in range(len(num_cont[l])):
            nums = num_cont[l][s].split(',')
            x = int(nums[0])
            y = int(nums[1])
            num_cont[l][s] = [x,y]
    return num_cont

Finally I compress the function to calculate Inflection points into a function named findCritical() which shall return the inflection points as:

masked_data = sine

all_contours = get_points(masked_data)

all_unique = list()

for cont in all_contours:
    u = get_unique(cont)
    for s in u:
        if len(s)>3:
            all_unique.append(s)


findCritical(all_unique,60)

The sine wave is detected as a combination of two contours one from range (0-15) and the other (16-210):Contour detected Contour detected

However this is the set of inflections points detected:

# The first two lists refer to one coordinate axis and the later two to the other
([array([14., 13., 12., 10.,  9.,  8.,  4.,  3.,  2.]),
  array([ 20.,  21.,  30.,  50.,  56.,  70.,  71.,  90., 100., 109., 116.,
         120., 140., 149., 158., 169., 189., 191., 199.])],
 [array([102., 102., 103., 105., 106., 107., 108., 108., 107.]),
  array([105., 106., 105., 105., 102., 105., 104., 105., 105., 106., 102.,
         105., 105., 106., 103., 106., 106., 104., 104.])])

As you can see that the coordinates are consecutive points. The same code snippet would produce desirable results for synthetic data of the sine wave with floating point coordinates but for the image with integer coordinate(which is the problem of interest), it fails. I really need inputs and suggestions on this.

  • There’s a lot of stuff going on in your code, it’s hard to read it all. But I think that you are sorting points in the function that removes duplicate points? If so, you are likely changing the order of the points in the contour, thereby screwing up any possibility to compute derivatives. Did you try plotting your sequence of points as a polygon? – Cris Luengo Jun 07 '20 at 15:39
  • Also, if the extracted contour duplicates points it is because it goes all the way around your shape, it doesn’t just follow a line, it follows the outline of the line. The second half of the contour would be redundant. So deleting that would lead to a unique set of points *in the right order*. – Cris Luengo Jun 07 '20 at 15:39
  • No, the relative arrangement of the points is not changed.. I've traced all the points and whenever the points tend to repeat I make a new contour instead – gowreesh mago Jun 07 '20 at 17:29

0 Answers0