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

Add Fast Marching method #45

Closed masadcv closed 1 year ago

masadcv commented 1 year ago

Fixes #42 and issues raised in #40

This PR implements Fast Marching Method as indicated in #42

It has the following changes:

tvercaut commented 1 year ago

For clarity, it may be good to clean up a bit the pixel queue implementation. Maybe adding some comments to refer back to the algorithm described in https://doi.org/10.1016/j.patrec.2006.10.010

Also, I would replace the hand-crafted heaps by a standard use of a std::priority_queue: https://en.cppreference.com/w/cpp/container/priority_queue

For the fast marching, could you fetch the fast-marching-method repo or add it as a git submodule rather than copying the file here?

masadcv commented 1 year ago

Many thanks! For pixel queue related work, I have added an issue for general code cleanup and improvements as your proposed (issue #46 ) as I would like to keep this PR for adding complete FMM functionality and renaming existing method to pixel queue.

I will look into FMM submodule bit, I would ideally like to keep installations hassle-free (e.g. pip install git+repo or pip install FastGeodis), which is one of the reason why I copied the file over. I will try submodules, if not then copy the working directory of FMM git as dependency like done in SimpleCRF here

masadcv commented 1 year ago

It would be good if you could verify FMM methods' output (if your schedule permits).

Here is a quick overview of what I am getting: For 2d:

fig2d

For 3d:

fig3d

masadcv commented 1 year ago

This now uses fmm external dependency as it is (not a submodule though) copied within ./dependency folder.

masadcv commented 1 year ago

@tvercaut I have been testing the fmm method in this branch. It seems like there are certain scenarios where this fails with error: RuntimeError: invalid arrival time (distance) -nan at index [xxx, xxx]

This happens especially for signed distance calculation when image is bigger and only a single seed point is selected. Here is a failing example using the scaled image (if you set scale to 1, everything works as expected): https://colab.research.google.com/drive/1BXbJcoEf9Ib9aXIsWjcMwNWlV5HRgBaq?usp=sharing

I am unsure what is causing this, do you know?

tvercaut commented 1 year ago

Fast marching uses the inverse of the speed and even the inverse of the square of the speed. The most likely explanation is that your speed contains too low values.

Looking at https://github.com/masadcv/FastGeodis/blob/100a4e4dc363d063b08b5680b28a3cd1d0db67b1/FastGeodis/geodis_fastmarch.cpp#L183 I would be more conservative. At least use sqrt(machine_espilon) rather than machine_espilon.

There might be other issues but I'll start with this.

Also, unrelated to numerical stability but you should use the high-accuracy eikonal solver here: https://github.com/masadcv/FastGeodis/blob/100a4e4dc363d063b08b5680b28a3cd1d0db67b1/FastGeodis/geodis_fastmarch.cpp#L47-L51

tvercaut commented 1 year ago

I added a fallback in case the Eikonal equation is degenerate in this branch: https://github.com/tvercaut/fast-marching-method/tree/negativediscriminant

Can you give it a try?

masadcv commented 1 year ago

I added a fallback in case the Eikonal equation is degenerate in this branch: https://github.com/tvercaut/fast-marching-method/tree/negativediscriminant

Can you give it a try?

Many thanks! Let me give it a try now.

masadcv commented 1 year ago

I have included both sqrt(eps) and degenerate fallback within this PR. It is stable now and works fine for all cases. I have also switched to high accuracy eikonal solver but have left in commented code for lower accuracy method if needed. I will go ahead and merge this after some local testing. If any further issues then please raise a bug and I can look into it.

Thanks once again for your help with FMM, it wouldn't have been possible without it.

tvercaut commented 1 year ago

Nice. Thanks!