elvirrey elvirrey - 7 months ago 15
Python Question

Reordering image skeleton coordinates to make interp1d work better

For a school project I'm analyzing the centerlines of some C. elegans images. I've managed to generate a reasonable threshold and I'm using

skimage.morphology.skeletonize
to generate the centerline:

enter image description here

Then I use
np.nonzero
to get the coordinates of the centerline with the eventual goal of parametrizing these points to gain some insight into the centerline geometry.

However, when I use
scipy.interpolate.interp1d
, I get this mess:
enter image description here
I'm fairly sure this happens because when
np.nonzero
looks for the nonzero values, it goes bottom-up and right-left and orders the points as such, which is why I get the zig-zag effect during interpolation. Is there any way that I can reorder these points so that
interp1d
plays nice with them?

Here is my code:

import cv2
import numpy as np
from matplotlib import pyplot as plt
from scipy import ndimage
from scipy import stats
from scipy import misc
from skimage.morphology import skeletonize
from scipy.interpolate import interp1d

"""Curved Worm"""

img = misc.imread("model_image_crop_curved.tif")
plt.imshow(img)
plt.show()

imgThresh = img>200
plt.imshow(imgThresh)
plt.show()

misc.imsave('model_image_crop_curved_binary.tif',imgThresh)

imgSkel = skeletonize(imgThresh)
plt.imshow(imgSkel)
plt.show()
misc.imsave('model_image_crop_curved_skeleton.tif',imgSkel)

cv2Skel = cv2.imread('model_image_crop_curved_skeleton.tif',cv2.IMREAD_GRAYSCALE)

skelCoord = np.nonzero(imgSkel)

x = skelCoord[1]
y = skelCoord[0]

plt.plot(x,y,'.')
plt.show()

i = np.arange(len(x))
interp_i = np.linspace(0,i.max(),5*i.max())

interpKind = 'linear'

xi = interp1d(i,x,kind=interpKind)(interp_i)
yi = interp1d(i,y,kind=interpKind)(interp_i)

fig, ax = plt.subplots()
ax.plot(xi,yi,'b-')
ax.plot(x,y,'ko')
plt.show()


And here are the points I get from
np.nonzero
:

[51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 45 46 47 48 49 50 68 69
70 71 72 73 74 75 40 41 42 43 44 76 77 78 36 37 38 39 79 80 34 35 32 33 31
30 30 29 28 28 27 27 27 27 27 27 27 27 26 26 26 26 26 26 27 27 27 27 27 27
27 28 28 29 30 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
50 51 52 53 54 55 56 56]
[16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 17 17 17 17 17 17 17 17
17 17 17 17 17 17 18 18 18 18 18 18 18 18 19 19 19 19 19 19 20 20 21 21 22
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
48 49 50 51 52 53 54 55 56 57 58 58 59 59 60 61 61 61 62 62 62 63 64 64 65
65 66 67 68 69 70 71 72]


EDIT

To solve the ordering problem, I followed ev-br's suggestion and began with the end points. I used a series of possible end point orientations and the
ndimage.binary_hit_or_miss
function to isolate the end points. Then I wrote a function that walked down the skeleton by checking each pixel's neighbors for the next in the skeleton, moving to that one, and keeping a running list. This running list became the set of ordered points I was looking for.

However, once I accomplished this after many hours of frustration (including an entire half hour spent agonizing over a problem that was easily solved by replacing
and
with
or
in one of my
if
statements), I realized that when I interpolate these data points they really don't provide much additional information. Since the data points I collected by walking along the skeleton are parametric in themselves, any parametric analysis I do can use those points. So while it is good that I figured out how to order the points, interpolation didn't turn out to be the end goal anyways.

If anyone wants to see the code I wrote to accomplish this, send me a message and I'll gladly share.

Answer

As you noticed, you need to order the nonzero pixels. Since you've skeletonized the image, you can start from an endpoint (whixh has exactly one neighbor) and walk along the path until you reach the other one. This way, you get an ordered list of pixel coordinates, which you then interpolate. Notice though that the question of how to parameterize a curve is not trivial, and you might need to do things beyond just tossing in interp1d

One keyword for searching the internet is "analyze skeleton".

Comments