itk.BinaryThinningImageFilter3D output an empty image

[qiang-zhang-neu]

Thank you very much for sharing the code. However, I met an error when using BinaryThinningImageFilter3D, which output an empty image.

My environment is:
win 10
python 3.5.6

My code is:

import itk, os
import SimpleITK as sitk
import numpy as np
import scipy.io as sio
npImg = np.zeros(shape=[256,256,256],dtype=np.int16)
npImg[50:100, :, :] = 1
sitkimg = sitk.GetImageFromArray(npImg, isVector=False)
print('-------------sitk image --------------------')
print(sitkimg)

itkImage = itk.GetImageFromArray(sitk.GetArrayFromImage(sitkimg), is_vector =False)
print('-------------itk image --------------------')
print(itkImage)
skeleton = itk.BinaryThinningImageFilter3D.New(itkImage)
skeletonImg = skeleton.GetOutput()
print('-------------skeleton image --------------------')
print(skeletonImg)

From the result of print, we can find that:
the size of sitk image is 256256256
the size of itk image is: 256256256
the size of skeleton image is: 256256256

Those, I don't know what's the wrong for my code.
Thanks for any suggestion!

[T4mmi]

Hi @qiang-zhang-neu
thanks for your feedback.
It's right to say that your example returns an empty image/volume (here is the minimal code to reproduce, needing itk itk-thickness3d numpy imageio):

import itk
import numpy as np
import imageio

# create the original test array
array = np.zeros((256,256,256),dtype=np.int16)
array[50:100, :, :] = 1
# cast to itk image
image = itk.GetImageViewFromArray(array)
# extract the skeleton
skeleton = itk.BinaryThinningImageFilter3D(image)
# cast back the skeleton from itk to numpy
output = itk.GetArrayViewFromImage(skeleton)
# simple way to see if output is empty:
print("original object size : {} vx".format(np.sum(array)))
print("skeleton size : {} vx".format(np.sum(output)))
# export as tiff objects fot visualization
imageio.volwrite("in.tif", array)
imageio.volwrite("out.tif", output)

It is right that with several "perfect" geometrical objects this algorithm has a strange behavior (see comments in the original post).
I personnaly observed these with spheres … but did not impact the computations with actual non-geomatricaly-perfect objects …

I Don't have any suggestions or answer for you … try testing using your real data and see how it works …
Alternatively you can try using the test data provided along with the module to see how it works with branch-like objects

[fusentasticus]

@T4mmi I managed to overcome Swig problems, see #33 . Now I'm hit by this problem.

For a real medical boolean 3D image, the result was actually sort of the expected 3D surface (for both methods in the module) except not really: almost all points in the skeleton are not there.

That is, the skeleton has withered into thousands of tiny disconnected islands, most of size 1. The erosion clearly is too aggressive.

Moreover, for the little tiff test data you mentioned above (part of the package), I get an all black result image.

So, yes, I see empty or almost empty outputs, too.

That's a shame, because the existing ITK thinning algorithm seems to be 2D based and it produces very bad artefacts when applied to 3D, so it's not really useable.

Do you know of any alternatives in other morphological software? The algorithmic problem of skeletonization appears to be solved in a number of ways, but the frustrating thing is that I can't find an implementation that works...

[phcerdan]

@fusentasticus You can try sgext, as commented in this ITK discourse post.
It was born in my PhD thesis, but recently has matured quite a lot, including 100% python bindings of al C++ functionality.

The documentation is good, but the "intros" and a documentation with a bit more narrative is lacking. However for a 3D thinning, the README might suffice. In python, you won't need more than:

pip install sgext

import sgext
sgext.read_binary_image?
sgext.create_distance_map? # To ensure you get a centered skeleton.
sgext.thin? 

Hope it helps, the algorithm is more reliable than the 3D thinning used in this module (it's just more recent research from digital topology groups), however it might use more memory and CPU.

If you have any doubt/question, open an issue there.
Hope it helps.

[fusentasticus]

@phcerdan Much thanks! I'll try it out. Also, the paper https://hal.archives-ouvertes.fr/hal-01217974/document that you cite on your site is a nice reference on trade-offs in how you even define the skeleton. It also has plenty of comparisons to at least some of the other work. And I gather that's the starting point for your own work.

I don't know how it relates to Hanno Homann's work Implementation of a 3D thinning algorithm. That's what this ITKThickness3D module is based on apparently.

[thewtex]

skeleton = itk.BinaryThinningImageFilter3D.New(itkImage)
skeletonImg = skeleton.GetOutput()

Before calling .GetOutput(), call .Update() on the produce the output.

Or, use the functional interface, i.e.

skeletonImg = itk.binary_thinning_image_filter3_d(itkImage)

[fusentasticus]

Before calling .GetOutput(), call .Update() on the produce the output.

Yes, it works now and it produces a wire skeleton. Thanks @thewtex.

However, I still think the claim in the README:

The skeletonization insure a minimal set of measurements to fully describe the object.

is incorrect (and not only for grammar!). For a thin sheet, a filamentous tree is calculated, but he union of the corresponding calculated spheres is not the object. That's because voxels are deleted as long as they don't destroy connectivity with no regard to geometry.

(@phcerdan's sgext algorithm is able to deal with surfaces that are to be reflected in the skeleton)

[charlotte10170]

@phcerdan Hi! thank you for your contribution. I would like to calculate the medial axis from the skeleton, so i need an extra information related to the distance information. Do you have any suggested function in sgext compatible with python?

[phcerdan]

@charlotte10170 you just need to pass a distance_map image (with the distance information) to the thin algorithm and the result will be the medial axis. Read sgext.thin?.
Please note that your input image needs to be isotropic (the spacing must be equal in all the axes) for the result to be centered in physical space.

Open an issue in https://github.com/phcerdan/sgext if you have further questions.

[charlotte10170]

@phcerdan Thanks a lot for your help, i solved my problem

[thewtex]

is incorrect (and not only for grammar!). For a thin sheet, a filamentous tree is calculated, but he union of the corresponding calculated spheres is not the object. That's because voxels are deleted as long as they don't destroy connectivity with no regard to geometry.

@fusentasticus makes sense. Could you please create a PR to update the grammar, description, and add a pointer to sgext in the README?

As it looks like the original issue is now resolved, closing this. If there are additional follow ups or other issues, let's please create new separate issues to track them.