Question on the meaning of parameters of generalised_geodesic3d

[anafyp]

Hi all,

I have recently discovered your implementation of FastGeodis. I am however unable to understand how I am supposed to make it work.

Say we have a 3D binary mask in tensor format (i.e [B, C, X, Y, Z]) if I wanted to simply do the distance transform of this mask as in scipy.ndimage.morphology.distance_transform_edt, how should I proceed?

To go further with my question I don't really understand the documentation on:

• image – input image, can be grayscale or multiple channels.
• softmask – softmask in range [0, 1] with seed information.
• spacing – spacing for 3D data
• v – weighting factor for establishing relationship between unary and spatial distances. To what should I set those if I just want to extract the euclidian distance transform of my mask ?

Thanks!

[masadcv]

Hi @anafyp ,

Thats a great question, many thanks for raising it.

This library is mainly aimed at calculating Geodesic distance transforms with a fast parallelized implementation. There is a trade-off between accuracy and speed, where most fast implementations are an approximation for the actual Geodesic distance transform. This also applies if you are doing Euclidean distance with the functions here, they may be fast but not exact computation but rather a close approximation.

Having said that, the library can be used to compute Euclidean distances. You need to provide an empty image (doesnt not matter what values you provide as for Euclidean distance this is not used but still required by the functions) and setting lamb = 0.0 will enable computing the Euclidean distance transform, which is similar to scipy edt function. Here is an example I implemented for your reference FastGeodisEuclideanDistanceTransformWithoutImage

Regarding your question about other parameters, please have a look at the simple demo colab samples provided here: Colab samples

I provide a brief overview about each parameter for your Euclidean distance use case:

• image: not relevant for Euclidean distance as it is needed for Geodesic distance only. You still need to provide an image here but it will not be used if you are only calculating Euclidean distance. This dummy image can be for example all 1s with same size as softmask below
• softmask: it could consist of points marked with 0s in a tensor of 1s for your Eucildean case. Points from which you wish to calculate distance (equivalent to (1-input) for scipy edt function)
• spacing: it is the spacing in your image for each dimension, e.g. [1.0, 1.0, 1.0] for isotropic spacing
• v: is weighting factor needed for Geodesic distance transforms that use soft masks. If you just need Euclidean distance and have discrete values in softmask i.e. either 0 or 1, then you can set this to a high value e.g. v=1e10
• lamb: lanb is 0 for Euclidean and 1.0 for Geodesic or in between if you want mix of both
• iter: for iterative methods these are iterations the method will run. Typically, it is set to 2 for 2d and 4 for 3d cases

I hope this helps.

[masadcv]

I forgot to add, if there is an interest in having specific Euclidean distance only functions, I can implement and include those within the library as well (so you wont need to set all the parameters that are specifically aimed at Geodesic distance)

[anafyp]

Thank you @masadcv for your response, that is also the behavior that I expected. The result however is very different from the one by scipy.

FastGeodis with the parameters you recommended produces the following output:

While scipy, with the expected result produces:

I think there is indeed an interest in having a specific Euclidian distance-only function, but for now, this distance looks more like a chessboard distance and not a euclidian one.

[masadcv]

Can you share this example so I can have a closer look at what you are trying?

[anafyp]

As i am working with medical data it might be hard to share data... I was mistaken in my previous comment I think FastGeodis computes the taxicab/manhattan distance and not the euclidian one. Here is a dummy example:

import torch
import numpy as np
import matplotlib.pyplot as plt
from FastGeodis import generalised_geodesic3d
from scipy.ndimage import distance_transform_edt, distance_transform_bf

inimage = torch.zeros((1, 1, 500, 500, 500))
inimage[0, 0, 250, 250, 250] = 1

dm_fastgeodis = generalised_geodesic3d(torch.zeros(inimage.shape),
                                       (1 - inimage),
                                       [1.0, 1.0, 1.0],
                                       1e10,
                                       0,
                                       4)
dm_fastgeodis = torch.squeeze(dm_fastgeodis).numpy()

inimage = np.zeros((500, 500, 500))
inimage[250, 250, 250] = 1

dm_scipy = distance_transform_bf(1 - inimage, 'euclidean')
dm_chess = distance_transform_bf(1 - inimage, "chessboard")
dm_taxi = distance_transform_bf(1 - inimage, "taxicab")

fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4)
ax1.imshow(dm_fastgeodis[250, :, :])
ax2.imshow(dm_scipy[250, :, :])
ax3.imshow(dm_chess[250, :, :])
ax4.imshow(dm_taxi[250, :, :])

ax1.set_title("FastGeodis")
ax2.set_title("Scipy Euclidian")
ax3.set_title("Scipy Chess")
ax4.set_title("Scipy taxicab")

And the output it produces:

[masadcv]

Many thanks for sharing the example.

Upon closer look, I found a minor bug in how local Euclidean distances were calculated. I have fixed this and released v1.0.1 on pypi with the fix.

Having said that, it still is an approximation and contains artifacts that are resulting from how these optimized algorithms compute distance transforms using raster scan methods. It is as good as these algorithms could be for Euclidean distance.

Here are the results for your 3d example using FastGeodis v1.0.1 (I have also added GeodisTK Euclidean case using their raster scan method for comparison):

[tvercaut]

As discussed earlier, on the long run, it would be good to have dedicated approaches for Euclidean distance maps.

Tensorflow has a GPU implementation of the Felzenszwalb & Huttenlocher Distance transform which could be ported here:

• https://www.tensorflow.org/addons/api_docs/python/tfa/image/euclidean_dist_transform
• https://github.com/tensorflow/addons/blob/master/tensorflow_addons/custom_ops/image/cc/kernels/euclidean_distance_transform_op_gpu.cu.cc

CUCIM has a GPU implementation of the Parallel Banding Algorithm (PBA) which may be even more appropriate for porting here:

• https://github.com/rapidsai/cucim/blob/branch-23.04/python/cucim/src/cucim/core/operations/morphology/_distance_transform.py

By the way, is the current output from Fast Marching any better?

[masadcv]

Many thanks @tvercaut for pointing out the reference implementations we can use for Euclidean distance maps. I will have a look and check how we can incorporate these within this package.

I have not checked Fast Marching for Euclidean distance maps yet. Will check and see if it is any better.

[tvercaut]

If it's of help for comparison purposes, the scikit-fmm package provide Euclidean distance map computation through fast marching: https://pythonhosted.org/scikit-fmm/#skfmm.distance

As far as their demo shows, their implementation doesn't seem to lead to star artefacts:

[masadcv]

Hi @anafyp , We have updated the Fast Marching Method to have accurate distance transform output (especially for Euclidean distance transform). Here is its output:

2D:

3D:

Could you try your examples with version 1.0.3 (make sure to refer to API docs to correct any function naming/arguments mismatch in new version) I believe this addresses your main concerns about Euclidean distance. Therefore I will close this issue, please reopen or create a new issue and I will look into helping you.