[BUG] "Fast marching" implementation doe not use Eikonal update formula

[tvercaut]

Describe the bug Having looked through the code, it seems that the implementation of geodis_fastmarch.cpp does not use any Eikonal update formula: https://github.com/masadcv/FastGeodis/blob/7101eb32a0541e0d5657c1cf7aba92811a640c57/FastGeodis/geodis_fastmarch.cpp#L120-L138 As such, it probably isn't really implementing fast marching.

To Reproduce Look at the output of the current fastmarching for Euclidean distance: https://github.com/masadcv/FastGeodis/blob/7101eb32a0541e0d5657c1cf7aba92811a640c57/samples/simpledemo2d.py#L106 Star artefacts are as pronounced as those in Toivanen approach. See also #40

Expected behavior A much better approximation of the distance map should be achievable through fast marching.

Additional context A good description of fast marching can be found in the following technical report: Bærentzen, J. A. (2001). On the implementation of fast marching methods for 3D lattices. http://www2.imm.dtu.dk/pubdb/pubs/841-full.html

A python wrapper for a basic c++ implementation of fast marching can be found in the scikit-fmm package.

A standalone c++, header-only implementation of fast marcing can be found here: https://github.com/tvercaut/fast-marching-method It would make sense to replace teh current implementation in the FastGeodis repo by that provided in fast-marching-method.

[tvercaut]

I have created some proof of concept pybind11 based python bindings for the above fast marching code: https://github.com/tvercaut/fast-marching-method/blob/master/bindings/python/py-fast-marching-method.cpp https://github.com/tvercaut/fast-marching-method/blob/master/tests/py-bindings-test.py

It would make sense to try and replace the current fast marching implementation by something along the same lines but for pytorch tensors. As such https://github.com/masadcv/FastGeodis/blob/master/FastGeodis/geodis_fastmarch.cpp would just provide some thin wrappers.

[masadcv]

Many thanks @tvercaut for providing proof of concept python bindings to the fast marching code!

I will have a closer look in the next few days to propose changes that utilize fast marching method. I may get back here if I have any questions on the method...

[masadcv]

Hi @tvercaut I had a look through the code. Many thanks!

I can understand to some extent how this is working - as far as I understand, this will only be implementing Euclidean distance transform (if not, please correct me with some pointers to what I may be missing)...

[tvercaut]

Fast marching takes a "speed" scalar field as input as well. So while you can't define a speed between 2 pixels directly, geodesic distances can still be approximated by assigning the inverse of the gradient magnitude at each pixel as the speed.

I am not sure how exactly "speed" is handled in the current ad-hoc version but a quick look seems to indicate L1 is used. https://github.com/masadcv/FastGeodis/blob/c8c51a9e122b40e86ee4bd908047d38a82f5247f/FastGeodis/geodis_fastmarch.cpp#L35-L38

Roughly replicating this with a proper fast marching can be seen in this notebook: https://colab.research.google.com/drive/1PNhl0jr-v5iJw0gEmLLJK4oeUlGqzDb7?usp=sharing

For a better approximation, an anisotropic fast marching methods would be needed but I think we can leave that for future work.

[masadcv]

Many thanks @tvercaut for clearing this. It is clear to me now what is needed here, I should have a fix in next few days...

[tvercaut]

Interestingly it looks like the algorithm currently referred to as "fastmarch" in this repo is essentially the pixel queue algorithm of Ikonen and Toivanen (rather than fast marching): Ikonen, L., & Toivanen, P. (2007). Distance and nearest neighbor transforms on gray-level surfaces. Pattern Recognition Letters, 28(5), 604-612. https://doi.org/10.1016/j.patrec.2006.10.010

[masadcv]

Ohhh great, we need a major renaming of methods then!

By the way, I have fast marching initial prototype working for 2D, for 3D I am currently debugging... Should have a PR this weekend! Thanks for the help, I basically added a light wrapper + speed for geodesic using varying speed method.

[masadcv]

I have added a PR with fix for this (still WIP but has initial prototype for FMM based method). Can you have a look and try when possible?

I have another question about generalised geodesic distance implementation using Fast Marching. I could not figure out how we could do it, for example a good description of generalised geodesic distance is in (Equation 9): Criminisi, Antonio, Toby Sharp, and Andrew Blake. "Geos: Geodesic image segmentation." European Conference on Computer Vision, Berlin, Heidelberg, 2008.

In particular, generalised variant utilises soft masks and parameter v that "establishes mapping between the unary beliefs and spatial distances."

For parameter v, i believe this could be incorporated in speed? I am unsure about soft masks as FMM uses boundary points instead of a mask... what do you think?

[tvercaut]

Good question re generalised geodesic distances. Actually how is it done with the pixel queue? I quickly went back to the journal paper

Criminisi, A., Sharp, T., Rother, C., & Pérez, P. (2010). Geodesic image and video editing. ACM Trans. Graph., 29(5), 134-1. https://doi.org/10.1145/1857907.1857910

Noting that wave-front methods encompass both pixel queue and fast-marching, it states

Application of wave-front methods to this new definition would be computationally prohibitive. Raster-scan methods, on the contrary, extends naturally to this case by simply adding a “no-move” option in each minimization step of the scan.

[tvercaut]

Also, it might be worth flagging that the pixel queue method is essentially a specific implementation of Dijkstra's algorithm. It should thus be relatively easy to unit test against say the scipy.sparse.csgraph.dijkstra implementation.

A quick example:

import skimage as skim
import scipy

def edge_distance_func(values0, values1, eucl_distances, lambdaw):
    """A parametrised edge function for complete image graphs.
    Parameters
    ----------
    values0 : array
        The pixel values for each node.
    values1 : array
        The pixel values for each neighbor.
    eucl_distances : array
        The Euclidean distance between each node and its neighbor.
    Returns
    -------
    edge_values : array of float
        The computed values: weighted pixel difference and euclidean distance
    """
    return lambdaw * np.abs(values0 - values1) + (1 - lambdaw) * eucl_distances

lamb = 0.8  # <-- Geodesic distance weight (1.0 purely geodesic, 0 purelay Euclidean)
edge_func = lambda v0, v1, d: edge_distance_func(v0, v1, d, lamb)
g, v = skim.graph.pixel_graph(image, 
                              mask=np.full_like(image, True, dtype=bool),
                              edge_function=edge_func,
                              connectivity=2)
geodesic_dist_dijkstra = scipy.sparse.csgraph.dijkstra(g,
                                                       directed=False,
                                                       indices=seed_indices,
                                                       min_only=True)
geodesic_dist_dijkstra = geodesic_dist_dijkstra.reshape(image.shape)

These parallels would eventually allow using additional optimisation tricks such as using an untidy queue instead of a proper priority queue: Yatziv, L., Bartesaghi, A., & Sapiro, G. (2006). O (N) implementation of the fast marching algorithm. Journal of computational physics, 212(2), 393-399. https://doi.org/10.1016/j.jcp.2005.08.005

[masadcv]

Good question re generalised geodesic distances. Actually how is it done with the pixel queue? I quickly went back to the journal paper

Criminisi, A., Sharp, T., Rother, C., & Pérez, P. (2010). Geodesic image and video editing. ACM Trans. Graph., 29(5), 134-1. https://doi.org/10.1145/1857907.1857910

Noting that wave-front methods encompass both pixel queue and fast-marching, it states

Application of wave-front methods to this new definition would be computationally prohibitive. Raster-scan methods, on the contrary, extends naturally to this case by simply adding a “no-move” option in each minimization step of the scan.

I need to read the paper again I think to see why it wont be possible.

For the 'pixel queue' method, it simply is an initialization of distance to softmask * v as outlined in the same paper for raster scan (section 2.2.2): "Raster scan methods extends naturally to this case, as follows. After initializing the distance map with scaled seed mask (D = νM), a first scan is executed from top-left to bottom-right, updating the distance map at the current pixel x according to...."

Here is where it is updated in comparison to GeodisTK (from which the code for pixel queue is based):

Initialize with v*softmask: https://github.com/masadcv/FastGeodis/blob/c8c51a9e122b40e86ee4bd908047d38a82f5247f/FastGeodis/geodis_fastmarch.cpp#L583

Update initial points with this: https://github.com/masadcv/FastGeodis/blob/c8c51a9e122b40e86ee4bd908047d38a82f5247f/FastGeodis/geodis_fastmarch.cpp#L495

In the non-generalised variants (e.g. in GeodisTK) this update was done for hard mask with a very large distance value.

Now coming back to the point that this may not work in wave-front methods, I believe the pixel queue implementation is not purely wave-front prop, but I need to double-check this with the code and linked papers...

[tvercaut]

I don't have a good intuition as to why a Dijkstra-type method (which includes pixel queue) would work for the generalised geodesic distance. It it does then the paper missed something. Do we have tests showing that it works?

As per the generalised geodesic distance formula

We could use a somewhat brute force approach for it: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csgraph.shortest_path.html

[tvercaut]

A simple test confirms that the current pixel queue method doesn't work for generalised geodesic distances: https://colab.research.google.com/drive/1PNhl0jr-v5iJw0gEmLLJK4oeUlGqzDb7?usp=sharing

[masadcv]

Thanks @tvercaut for reporting this. I have added an additional task to remove 'generalised' implementation attempt from wave-front methods. I will get back to this after added fmm functionality.

[masadcv]

I have added the updated FMM method with Eikonal solver as mentioned here. some of the smaller tasks mentioned here relevant to FMM have been address with this pr #45 along with adding FMM method.

For other tasks, I have opened issues and put them in ToDo - which I will look into later. In case I missed anything, please open an issue with it, so I can track and address them.

Thanks once again for your help with this.