masadcv / FastGeodis

Fast Implementation of Generalised Geodesic Distance Transform for CPU (OpenMP) and GPU (CUDA)
https://fastgeodis.readthedocs.io
BSD 3-Clause "New" or "Revised" License
90 stars 14 forks source link

how to get equivalent result #53

Closed lfz closed 4 months ago

lfz commented 1 year ago

Thanks for the fast fix for the padding issue!

I found the definition of "geodesic distance transform (GDT)" seems different from what I understand.

I need a function like "bwdistgeodesic" in (matlab), in this senario, it is very close to the traditional distance transform, which is also used for binary image. The only difference is that, the distance calculation cannot transmit through the "0" values in GDT.

but when I try FastGeodis, I found the basic logic differs from "bwdistgeodesic". the distance value only grows when face image gradient

import numpy as np
import torch
import FastGeodis
BW = np.array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 11],
                            [1, 1, 1, 1, 1, 1, 0, 0, 1, 1],
                            [1, 1, 1, 1, 1, 1, 0, 0, 1, 1],
                            [1, 1, 1, 1, 1, 1, 0, 0, 1, 1],
                            [0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
                            [0, 0, 0, 0, 1, 1, 0, 1, 1, 0],
                            [0, 1, 0, 0, 1, 1, 0, 0, 0, 0],
                            [0, 1, 1, 1, 1, 1, 1, 0, 1, 0],
                            [0, 1, 1, 0, 0, 0, 1, 1, 1, 0],
                            [0, 0, 0, 0, 1, 0, 0, 0, 0, 0]])

C = [1, 2, 3, 3, 3]
R = [3, 3, 3, 1, 2]
seed_img = np.zeros_like(BW)
seed_img[C,R] = 1
v = 1e10
# lamb = 0.0 (Euclidean) or 1.0 (Geodesic) (or (0.0, 1.0) (mixture)
lamb = 1.0
iterations = 2
geodesic_dist = FastGeodis.generalised_geodesic2d(
                        torch.from_numpy(BW).cuda()[None,None ].float(),
                        1-torch.from_numpy(seed_img).cuda()[None,None],
                        v, lamb, iterations)

but i get the following result: image

what I expect is: image

as far as I understand, this lib is designed to use image gradient as energy barrier, is there a simple way for me to reproduce the expected result?

lfz commented 1 year ago

I figured out a walkaround:

  1. do traditional 3d edt firstly (edt)
  2. set the masked area of edt to a large value, e.g. edt[mask] = 1e4
  3. then we use generalised_geodesic, set the result mask to nan or others

but I think this is not a good way to do so, is it possible to use more elegant method and incorporate this function in the lib

besides, I also want to get the nearest seed point for each point in the image, is it possible

lfz commented 1 year ago

In my scenario, it is close to a watershed algorithm, I got several seed points, and what to segment the binary mask to different basins that bound to those seed points. So I want the function can return the nearest seed point index