I have read quite some research papers talking about extracting a skeleton from a 2D image using the Euclidian Distance Transform or other kinds of distance transforms. I'm trying to implement this on a 3D mesh, using the scipy.ndimage.distance_transform_edt
function, but I don't get the expected result.
import trimesh as tm
import tkinter as tk
from tkinter import filedialog
import pyvista as pv
import numpy as np
from scipy.ndimage import distance_transform_edt
root = tk.Tk()
root.withdraw()
root.title('Pick patient data')
np.set_printoptions(precision=5)
file_path_2 = filedialog.askopenfilename()
chosenMesh_2 = tm.load_mesh(f'{file_path_2}')
volume_2 = tm.voxel.creation.voxelize(mesh=chosenMesh_2,
pitch=chosenMesh_2.extents.max() / 200)
array_2 = volume_2.matrix
filled_2 = ndi.binary_fill_holes(array_2)
distance_map = distance_transform_edt(filled_2)
x_skel2, y_skel2, z_skel2= np.where(distance_map == 1) #According to distance_map
values
coordinates2 = []
for x2,y2,z2 in zip(x_skel2,y_skel2,z_skel2):
coordinates2.append((x2,y2,z2))
plotter2 = pv.Plotter()
plotter2.add_points(np.array(coordinates2),color='black')
plotter2.show()
Output of distance_map
is now:
[ 0. 1. 1.41421 1.73205 2. 2.23607 2.44949 2.82843
3. 3.16228 3.31662 3.4641 3.60555 3.74166 4. 4.12311
4.24264 4.3589 4.47214 4.58258 4.69042 4.89898 5. 5.09902
5.19615 5.38516 5.47723 5.65685 5.74456 5.83095 5.91608 6.
6.08276 6.16441 6.32456 6.40312 6.48074 6.55744 6.63325 6.7082
6.78233 6.9282 7. 7.07107 7.14143 7.2111 7.28011 7.34847
7.48331 7.54983 7.61577 7.68115 7.81025 7.87401 8. 8.06226
8.12404 8.18535 8.24621 8.30662 8.3666 8.48528 8.544 8.60233
8.66025 8.7178 8.77496 8.83176 8.94427 9. 9.05539 9.11043
9.16515 9.21954 9.27362 9.38083 9.43398 9.48683 9.53939 9.64365
9.69536 9.79796 9.84886 9.89949 9.94987 10. 10.04988 10.0995
10.19804 10.24695 10.29563 10.34408 10.3923 10.44031 10.48809 10.63015
10.67708 10.72381 10.77033 10.81665 10.86278 10.95445 11. 11.04536
11.09054 11.18034 11.22497 11.31371 11.35782 11.40175 11.44552 11.48913
11.53256 11.57584 11.6619 11.7047 11.74734 11.78983 11.83216 11.87434
11.91638 12. 12.04159 12.08305 12.20656 12.24745 12.32883 12.36932
12.40967 12.4499 12.52996 12.56981 12.64911 12.68858 12.72792 12.80625
12.84523 12.96148 13.11488 13.15295 13.37909]
As you can see, this basically plots the whole volume instead of the skeleton but there are not points with a distance smaller than 1. Is this a mistake in the scipy.ndimage.distance_transform_edt
function or am I doing something wrong?
If it's not possible to calculate the skeleton this way, are there other skeletonization algorithms using some kind of distance transforms that can be used in Python?
Thanks in advance!