Distance transform to calculate centerline/skeleton

27 views Asked by At

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]

I get the following plot: This should plot the centerline but it shows the whole volume

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!

0

There are 0 answers