Draw a center curve in a curved object using CV2

Question

I want to draw a line across the center of a curved object.

As an example: an image of a banana is given, the orientation might change with different images and also there might be more than one curve in the object, but the task is the same.

To determine the length of the object, it's required to calculate (interpolate) the center line of the contour of the object from start to end, then I can calculate the length of the interpolated line. This is what I thought so far.

But now there is the tricky part, determining the contour of the object using python and cv2 is no problem, this works well. But calculating the center line to determine its length is a struggle.

Edit

The goal of the script should be to measure the length and area of the worm, so I don't have to measure the values manually for hundreds of images.

Input image:

My computation so far for the contour (green line):

What I want (only drawn by hand as an example):

Code used so far (without the "center line" I want, because I don't have any idea how to start). My idea using the convex hull, building the skeleton and using it does not work as expected, because the convex hull is to large (caused by the concave part), combining it with polyDP also does not work because the polyDP often misses parts of the worm and the combined result is also bad.

import numpy as np
import cv2
import os

draw_windows = True  ## change fo False for no windows only calc


def calc_values(filename):
    img = cv2.imread(filename)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    gray = cv2.GaussianBlur(gray, (7, 7), 0)
    ret, thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)

    drawWindow('thresh', thresh)

    edged = cv2.Canny(gray, 50, 100)
    edged = cv2.dilate(edged, None, iterations=1)
    edged = cv2.erode(edged, None, iterations=1)

    drawWindow('edged', edged)


    contours, _ = cv2.findContours(edged, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    # Assume the largest contour corresponds to the worm
    if contours:
        largest_contour = max(contours, key=cv2.contourArea)

        # Draw the contour on the original image
        image_with_contour = cv2.cvtColor(edged, cv2.COLOR_GRAY2BGR)
        cv2.drawContours(image_with_contour, [largest_contour], -1, (0, 255, 0), 2)
        cv2.drawContours(image_with_contour, contours, -1, color=(255, 255, 255), thickness=cv2.FILLED)

        # Display the original image with the detected contour
        drawWindow('Worm with Contour', image_with_contour)


def drawWindow(window_name, image):
    if draw_windows:
        cv2.imshow(window_name, image)
        cv2.waitKey(0)
        cv2.destroyAllWindows()


def main():
    directory = "input"
    for filename in os.listdir(directory):
        file = os.path.join(directory, filename)
        calc_values(file)


if __name__ == "__main__":
    main()

(I know code quality is not the best so far, but it started as a quick and dirty "project" :D)

Answer 1

I googled a lot to figure out what is available. I ended up with using input as above input_77560561.jpg:

For the code, see the question and my previous answer and Python Image - Finding largest branch from image skeleton answer that uses fil_finder library FilFinder _ GitHub:

import numpy as np
import cv2

draw_windows = True  ## change fo False for no windows only calc

def ResizeWithAspectRatio(image, width=None, height=None, inter=cv2.INTER_AREA):
    dim = None
    (h, w) = image.shape[:2]

    if width is None and height is None:
        return image
    if width is None:
        r = height / float(h)
        dim = (int(w * r), height)
    else:
        r = width / float(w)
        dim = (width, int(h * r))

    return cv2.resize(image, dim, interpolation=inter)



def calc_values(filename):
    img = cv2.imread(filename, 1)
    
    print('FILENAME : ', type(img) , img.shape)
    
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    gray = cv2.GaussianBlur(gray, (7, 7), 0)
    ret, thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)

    drawWindow('thresh', thresh)
    
    cv2.imwrite("thresh_worm.png", thresh)

    edged = cv2.Canny(gray, 50, 100)
    edged = cv2.dilate(edged, None, iterations=1)
    edged = cv2.erode(edged, None, iterations=1)

    drawWindow('edged', edged)

    contours, hierarchy = cv2.findContours(edged, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    
    if contours:
        
        cntsSorted = sorted(contours, key=lambda x: cv2.contourArea(x) , reverse =  True)
        
        for index , i   in enumerate(cntsSorted) :
            
            print('i ************************************************************ ' ,index ,type(index))
            
            print('# \ncontour : ' ,index , 'Area ?? : ', cv2.contourArea(i) , "size : ", i.shape)
    
    # Assume the largest contour corresponds to the worm
    if contours:
        # largest_contour = max(contours, key=cv2.contourArea)
        
        largest_contour = cntsSorted[0]

        # Draw the contour on the original image
        image_with_contour = cv2.cvtColor(edged, cv2.COLOR_GRAY2BGR)
        cv2.drawContours(image_with_contour, [largest_contour], -1, (0, 255, 0), 2)
        cv2.drawContours(image_with_contour, contours, -1, color=(255, 255, 255), thickness=cv2.FILLED)
        
        print('CV2.filled : ' , cv2.FILLED)  # CV2.filled :  -1

        # Display the original image with the detected contour
        drawWindow('Worm with Contour', image_with_contour)
                
        print('largest_contour : ' ,largest_contour  , type(largest_contour))
        
        # drawing = np.zeros((largest_contour.shape[0], largest_contour.shape[1] , 3))
        drawing = np.zeros((edged.shape[0], edged.shape[1] , 3))
        
        print(drawing.shape)
                
        #### SAME THING before had largest_contour = cntsSorted[0]
        # cv2.drawContours(drawing, [largest_contour] , -1 , color = (0,255,0) , thickness = cv2.FILLED)
        cv2.drawContours(drawing, [cntsSorted[0]] , -1 , color = (0,255,0) , thickness = cv2.FILLED)
        
        cv2.imwrite("drawing.png", drawing)
        
        drawWindow('Worm Contour', drawing)
        
        drawing_gray = cv2.imread( 'drawing.png' , 0 )
        
        thinned = cv2.ximgproc.thinning(drawing_gray, thinningType = cv2.ximgproc.THINNING_ZHANGSUEN)
        
        cv2.imwrite("thinned_worm.png", thinned)
        
        drawWindow('thinned_worm', thinned)
        
        #### code from https://stackoverflow.com/questions/53481596/python-image-finding-largest-branch-from-image-skeleton
        from fil_finder import FilFinder2D
        import astropy.units as u
        
        skeleton = thinned
        
        fil = FilFinder2D(skeleton, distance=250 * u.pc, mask=skeleton)
        fil.preprocess_image(flatten_percent=85)
        fil.create_mask(border_masking=True, verbose=False,
        use_existing_mask=True)
        fil.medskel(verbose=False)
        fil.analyze_skeletons(branch_thresh=40* u.pix, skel_thresh=10 * u.pix, prune_criteria='length')
      
        drawWindow('skeleton', fil.skeleton_longpath)
        cv2.imwrite("skeleton.png", fil.skeleton_longpath*255)

        skel = fil.skeleton

        print('\nSkel : ',type(skel) , skel.shape, skel.size , skel.ndim , np.max(skel) , np.min(skel) , np.unique(skel))
            
        original = img
        
        mask = fil.skeleton_longpath
        
        print('\nmask : ',type(mask) , mask.shape, mask.size , mask.ndim , np.max(mask) , np.min(mask) , np.unique(mask))
        
        mask_dilated = cv2.dilate(mask, np.ones((4, 4)))
        
        result = original.copy()
        
        for i in range(original.shape[0]):
            for j in range(original.shape[1]):
                result[i, j] = [0,0,255] if mask_dilated[i, j] == 1.0 else result[i, j]

        print('RESULT : ', result.shape)
        cv2.imwrite('overlay_1.png', result)   # saves modified image to result.png
        drawWindow('overlay_1', result)
        
        original = drawing
        
        result = original.copy()
        
        for i in range(original.shape[0]):
            for j in range(original.shape[1]):
                result[i, j] = [0,0,255] if mask_dilated[i, j] == 1.0 else result[i, j]

        print('RESULT : ', result.shape)
        cv2.imwrite('overlay_2.png', result)   # saves modified image to result.png
        drawWindow('overlay_2', result)

def drawWindow(window_name, image):
    if draw_windows:
        
        resize = ResizeWithAspectRatio(image, width=600)
        
        cv2.imshow(window_name, resize)
        
        cv2.moveWindow(window_name, 600, 200)
        
        cv2.waitKey(0)
        cv2.destroyAllWindows()


def main():
    
    calc_values('input_77560561.jpg')


if __name__ == "__main__":
    main()

Output in order of images, some of the images shown by the script are missing:

Tresh image:

edged image, missing shown in script.

Worm with contour, missing shown in script.

Worm contour filled:

Thinned worm contour filled:

Skeleton, that is thinned worm contour filled with removed branch by FilFinder, because of how the contour is detected, this is not the best line, but as mentioned in the previous answer by me to this post I wasn't able to produce a better worm contour filled, by working on this we should be able to get desidered results:

Overlays to input and contour filled, line is dilated by OpenCV:

2 :

1 :

As in previous answer commenting out # gray = cv2.GaussianBlur(gray, (7, 7), 0) I can get a better single line:

with overlays:

As for the length of the line*** could it be approximated by the total number of pixels of the skeleton (i.e. print('\n\nskeleton_lenght_approx : ', np.sum(mask)) ? If not it will have to be calculated starting by one of the ends and adding the distance (+1 for x+-1 or y+-1 , or 1.4 for x,y+-1 for diagonal pixels) unless filFinder has the data we need hidden somewhere.

*** Note worth reading:

How different can the number of pixels in a straight line be to its real length

to calculate the skeleton length, here starting from saved file skeleton.png , I devised this code, probably not the best way but only one I was able to conceive:

import numpy as np
import cv2
import math
from PIL import Image

def neighbors_coords(matrix: np.ndarray, x: int, y: int):
    """
    
    stolen from https://stackoverflow.com/questions/73811239/query-the-value-of-the-four-neighbors-of-an-element-in-a-numpy-2d-array
    
    """

    x_len, y_len = np.array(matrix.shape) - 1
    nbr = []
    if x > x_len or y > y_len:
        return nbr
    if x != 0 : 
        if matrix[x-1][y] == 1:
            nbr.append((x-1,y))
        if y != 0:
            if matrix[x-1][y-1] == 1 :
                nbr.append((x-1,y-1))
        if y != y_len:
            if matrix[x-1][y+1] == 1:
                nbr.append((x-1,y+1))
    if y != 0:
        if matrix[x][y-1] == 1:
            nbr.append((x, y-1))
        if x != x_len:
            if matrix[x+1][y-1] == 1 :
                nbr.append((x+1 , y-1))
    if x != x_len:
        if matrix[x+1][y] == 1 :
            nbr.append((x+1 , y))
        if y != y_len:
            if matrix[x+1][y+1] == 1 :
                nbr.append((x+1 ,y+1))
    if y != y_len:
        if matrix[x][y+1] == 1:
            nbr.append((x, y+1))
            
    # print('nbr : ', nbr , x_len, y_len)
    
    nbr_dist = []
    
    for i in nbr :
        
        dist = math.dist([x,y], [i[0],i[1]])
        
        nbr_dist.append(((i[0],i[1]) , dist))
        

    nbr_dist.sort(key=lambda tup: tup[1] , reverse = False)  # sort points to get closest one first
    
    # print('nbr_dist : ', nbr_dist , x , y)
    
    return nbr_dist


img = cv2.imread("skeleton.png", cv2.IMREAD_UNCHANGED)

img[img==255] = 1


neighbor_kernel = np.uint8([
    [1, 1, 1],
    [1, 0, 1],
    [1, 1, 1]])

neighbors_count = cv2.filter2D(img.astype(np.uint8), cv2.CV_8U, neighbor_kernel)

endpoint_indices = [ (i, (y,x)) for (y,x) , i in np.ndenumerate(img) if img[y,x] == 1 and neighbors_count[y,x] == 1]

print('\n\nendpoint_indices on neighbors_count : ', endpoint_indices, type(endpoint_indices) , len(endpoint_indices))



start = endpoint_indices[0][1]

img[start[0]][start[1]] = 0



print('\nstart , : ', start)


cnt = 0

coords = []

lenght = 0



coords.append((start , 0))

while np.sum(img) > 0 :
    
    # print('first : ' , first)
    
    second_next = neighbors_coords(img, start[0] , start[1])[0]
    
    # print('second_next : ', second_next)
    
    img[second_next[0][0]][second_next[0][1]] = 0
    
    # print('second_next : ', second_next)
    
    start = second_next[0]
    
    # print('start : ', start)
    
    coords.append(second_next)
    
    lenght += second_next[1]
    
    cnt +=1
    


print('cnt : ', cnt  )

# print('\n\nCoordinates : ' , coords ,'len(coords) : ',  len(coords))

print('\n\nlen(coords) : ' , len(coords))
    
print('\n\nendpoint_indices: ' , endpoint_indices)

print('\n
"""

this bit just to save a file containing all identified coordinates

to check that my script is working right 

"""
\nlenght : ', lenght , '   len(coords) : ', len(coords))

points_img = np.zeros((img.shape[0], img.shape[1], 4)).astype(np.uint8)

for i in coords :
    
    # print(i)
    
    points_img[i[0][0]]  [i[0][1]] = (255,0,0,255) 
    

image3 = Image.fromarray(points_img.astype(np.uint8) , 'RGBA')

image3.save('check_Test_WALK.png')

One consideration about fil.FilFinder2D and its results fil.analyze_skeletons --> fil.skeleton_longpath, using below image , thinned_worm.png :

with following code:

import numpy as np
import cv2
import math
from PIL import Image

from fil_finder import FilFinder2D
import astropy.units as u

def neighbors_coords(matrix: np.ndarray, x: int, y: int):
    """
    
    stolen from https://stackoverflow.com/questions/73811239/query-the-value-of-the-four-neighbors-of-an-element-in-a-numpy-2d-array
    
    """

    x_len, y_len = np.array(matrix.shape) - 1
    nbr = []
    if x > x_len or y > y_len:
        return nbr
    if x != 0 : 
        if matrix[x-1][y] == 1:
            nbr.append((x-1,y))
        if y != 0:
            if matrix[x-1][y-1] == 1 :
                nbr.append((x-1,y-1))
        if y != y_len:
            if matrix[x-1][y+1] == 1:
                nbr.append((x-1,y+1))
    if y != 0:
        if matrix[x][y-1] == 1:
            nbr.append((x, y-1))
        if x != x_len:
            if matrix[x+1][y-1] == 1 :
                nbr.append((x+1 , y-1))
    if x != x_len:
        if matrix[x+1][y] == 1 :
            nbr.append((x+1 , y))
        if y != y_len:
            if matrix[x+1][y+1] == 1 :
                nbr.append((x+1 ,y+1))
    if y != y_len:
        if matrix[x][y+1] == 1:
            nbr.append((x, y+1))
            
    # print('nbr : ', nbr , x_len, y_len)
    
    nbr_dist = []
    
    for i in nbr :
        
        dist = math.dist([x,y], [i[0],i[1]])
        
        nbr_dist.append(((i[0],i[1]) , dist))
        

    nbr_dist.sort(key=lambda tup: tup[1] , reverse = False)  # sort points to get closest one first
    
    # print('nbr_dist : ', nbr_dist , x , y)
    
    return nbr_dist


def main(filename, thinned) :
    
    for i in (range(10))  :
        
        skeleton = cv2.imread(thinned , cv2.IMREAD_UNCHANGED)
        
        skeleton[skeleton == 255] = 1
        
        
        
        fil = FilFinder2D(skeleton, distance=250 * u.pc, mask=skeleton)
        fil.preprocess_image(flatten_percent=85)
        fil.create_mask(border_masking=True, verbose=False,
        use_existing_mask=True)
        fil.medskel(verbose=False)
        fil.analyze_skeletons(branch_thresh=40* u.pix, skel_thresh=10 * u.pix, prune_criteria='length')
      
        # drawWindow('skeleton', fil.skeleton_longpath)
        cv2.imwrite("skeleton"+str(i)+'.png', fil.skeleton_longpath*255)
    
    
    lenghts_calculated = []
    
    for i in (range(10))  :

        img = cv2.imread(filename.split('.')[0]+str(i)+'.png' , cv2.IMREAD_UNCHANGED)
        
        img[img==255] = 1
        
        
        neighbor_kernel = np.uint8([
            [1, 1, 1],
            [1, 0, 1],
            [1, 1, 1]])
        
        neighbors_count = cv2.filter2D(img.astype(np.uint8), cv2.CV_8U, neighbor_kernel)
        
        endpoint_indices = [ (i, (y,x)) for (y,x) , i in np.ndenumerate(img) if img[y,x] == 1 and neighbors_count[y,x] == 1]
        
        print('\n\nendpoint_indices on neighbors_count : ', endpoint_indices, type(endpoint_indices) , len(endpoint_indices))
        
        
        
        start = endpoint_indices[0][1]
        
        img[start[0]][start[1]] = 0
        
        
        
        print('\nstart , : ', start)
        
        
        cnt = 0
        
        coords = []
        
        lenght = 0
        
        
        
        coords.append((start , 0))
        
        while np.sum(img) > 0 :
            
            # print('first : ' , first)
            
            second_next = neighbors_coords(img, start[0] , start[1])[0]
            
            # print('second_next : ', second_next)
            
            img[second_next[0][0]][second_next[0][1]] = 0
            
            # print('second_next : ', second_next)
            
            start = second_next[0]
            
            # print('start : ', start)
            
            coords.append(second_next)
            
            lenght += second_next[1]
            
            cnt +=1
            
        
        
        print('cnt : ', cnt  )
        
        # print('\n\nCoordinates : ' , coords ,'len(coords) : ',  len(coords))
        
        print('\n\nlen(coords) : ' , len(coords))
            
        print('\n\nendpoint_indices: ' , endpoint_indices)
        
        print('\n\nlenght : ', lenght , '   len(coords) : ', len(coords))
        
        lenghts_calculated.append((lenght, len(coords)))
        
        
        # """
        
        # this bit just to save a file containing all identtified coordinates
        
        # to check that my script is working right 
        
        # """
        
        # points_img = np.zeros((img.shape[0], img.shape[1], 4)).astype(np.uint8)
        
        # for i in coords :
            
        #     # print(i)
            
        #     points_img[i[0][0]]  [i[0][1]] = (255,0,0,255) 
            
        
        # image3 = Image.fromarray(points_img.astype(np.uint8) , 'RGBA')
        
        # image3.save('check_Test_WALK.png')
        

        
        
    for i , value in enumerate(lenghts_calculated) :
        
        print(i+1, '____________' , value[0] ,' vs number of pixels ' , value[1])
        
        
    
    return lenghts_calculated
    

filename = "skeleton.png"

thinned = "thinned_worm.png"

val = main(filename, thinned)


val_measure = [i[0] for i in val]

val_numb = [i[1] for i in val]


def get_change(current, previous):
    if current == previous:
        return 100.0
    try:
        return (abs(current - previous) / previous) * 100.0
    except ZeroDivisionError:
        return 0


def mean(data):
    """Return the sample arithmetic mean of data."""
    n = len(data)
    if n < 1:
        raise ValueError('mean requires at least one data point')
    return sum(data)/n # in Python 2 use sum(data)/float(n)

def _ss(data):
    """Return sum of square deviations of sequence data."""
    c = mean(data)
    ss = sum((x-c)**2 for x in data)
    return ss

def stddev(data, ddof=0):
    """Calculates the population standard deviation
    by default; specify ddof=1 to compute the sample
    standard deviation."""
    n = len(data)
    if n < 2:
        raise ValueError('variance requires at least two data points')
    ss = _ss(data)
    pvar = ss/(n-ddof)
    return pvar**0.5


print('value measured , mean : ', mean(val_measure) ,'  SD : ' ,stddev(val_measure, ddof=1)) 

print('value number of pixel  , mean : ', mean(val_numb) ,'  SD : ' ,stddev(val_numb, ddof=1)) 

I get following output:

........
.......
lenght :  1316.467170879763    len(coords) :  1129
1 ____________ 1316.4671708797628  vs number of pixels  1129
2 ____________ 1315.295598004509  vs number of pixels  1127
3 ____________ 1317.0529573173897  vs number of pixels  1130
4 ____________ 1315.8813844421359  vs number of pixels  1128
5 ____________ 1318.8103166302703  vs number of pixels  1133
6 ____________ 1314.709811566882  vs number of pixels  1126
7 ____________ 1315.881384442136  vs number of pixels  1128
8 ____________ 1315.295598004509  vs number of pixels  1127
9 ____________ 1312.9524522540014  vs number of pixels  1123
10 ____________ 1316.467170879763  vs number of pixels  1129
value measured , mean :  1315.881384442136   SD :  1.5374956741211723
value number of pixel  , mean :  1128.0   SD :  2.6246692913372702
.....

Overlay of 5 of the 10 (0 to 9 results : skeleton#.png") longest path in files:

and magnification of a small part showing how the fil.skeleton_longpath changes with each run of the algorithm runs [do not know if I made any mistake somewhere so far, don't know how algorithm works]:

ADDENDUM

Had a chat with filFinder developer , was really kind to me, and pointed me to lengths that is :

def lengths(self, unit=u.pix):
        '''
        Return longest path lengths of the filaments.

        Parameters
        ----------
        unit : `~astropy.units.Unit`, optional
            Pixel, angular, or physical unit to convert to.
        '''
        pix_lengths = np.array([fil.length().value
                                for fil in self.filaments]) * u.pix
        return self.converter.from_pixel(pix_lengths, unit)

(it is an attribute of class fil_finder.FilFinder2D(...)),

as explained to me its values is different from the geometric lenght of the fil.skeleton_longpath image because :

There can be a moderate discrepancy between the geometric length and the one filfinder finds from the longest path. FilFinder generalizes the identification of the intersections between branches along the longest path to allow them to be made up of several pixels. The distance it uses if then a median of the intersection pixels, added to the geometric distances along each branch.

There's also an option when running FilFinder2D.find_widths where 2 * width is added to the filament length to account for the shortening from the skeletonization process: https://fil-finder.readthedocs.io/en/latest/api/fil_finder.FilFinder2D.html#fil_finder.FilFinder2D.find_widths. If you are running that step, you can disable this with add_width_to_length=False