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:
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:
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.
