Expand the covered use cases by allowing to specify a edge weight function

[tvercaut]

Dijkstra's algorithm finds shortest paths based on weigths assigned to edges in a graph. As far as I understand it, in dijkstra3d, given a field image, a directed graph is implicitly constructed by using a voxel connectivity pattern and the edge weights are implicitly given by the value of the field image at the target of the edge (as per the readme):

For given input voxels A and B, the edge weight from A to B is B and from B to A is A.

While this covers some important use cases, it would also be interesting to allow for a more generic means of computing the edge weigths while retaining the implicit graph handling of dijkstra3d. In particular, allowing the user to specify a function taking as input the value of the field image at the edge source and edge target but also taking as input the type of edge (center->left, top-right->center, bottom->center, etc.) would be great.

I have created a simple python notebook to illustrate this concept in 2D but relying on networkx as the backend:
https://colab.research.google.com/drive/196EPtBO0BmIjYmi6RlBvLlD_JGyAE9zB?usp=sharing

I am not entirely sure how this should best be implemented in dijkstra3d but I guess one could create a functor in c++ (let's call it edgeweightfunc). edgeweightfunc would then be called instead of fetching only the edge target value (which seems to be done with delta = static_cast<float>(field[neighboridx]); in dijkstra3d:

dijkstra3d/dijkstra3d.hpp

Line 310 in 48dffd9

delta = static_cast<float>(field[neighboridx]); // high cache miss

For brainstorming purposes, a basic concept of what such an edgeweightfunc could look like (in 2D) and in c++ is below:

#include <iostream>
#include <cmath>
#include <array>

int main()
{
    // Encode edge type in 2D (xf: x forward - xb: x backward)
    // There are surely better ways of doing this including using
    // boolean operations to combine x, y and z modes
    // see e.g. https://en.cppreference.com/w/cpp/named_req/BitmaskType
    enum EdgeType : int { xf, xb, yf, yb, xfyf, xfyb, xbyf, xbyb };
    
    // This is the functor we would pass to dijkstra::distance_field3d
    auto edgeweightfunc = [](auto edgesourceval, auto edgetargetval, EdgeType edgetype ) {
        constexpr std::array<double,8> edgedists2 { 1., 1., 1., 1., 2., 2., 2., 2. };
        // we could also use std::hypot
        auto squarefunc = [](auto x) {return x*x;};
        return std::sqrt( edgedists2[edgetype] + squarefunc(edgesourceval - edgetargetval) );
    };
    
    std::cout << "Edge weight: " << edgeweightfunc(5.0,3.0,EdgeType::xbyb) << std::endl;
}

I haven't looked at how such c++ functors can be bound to python but I assume a solution can be found...

[william-silversmith]

Hi Tom, a C++ functor sounds pretty reasonable and a good way to extend the usefulness of the library. There could be a performance impact if the functor does not get inlined during compilation.

As for Python, it would be possible to pass parameters to precompiled functors, but it wouldn't be (practically) possible to pass a wholly new function to C++ (since it would have to be recompiled and hot linked (!)). A Julia reimplementation would probably give you the desired interactivity and performance (due to its integrated JIT) but I haven't written much of it. Is there a predefined functor structure you think would be particularly useful?

[tvercaut]

Thanks. For python, having only a predefined functor structure with only the parameters being exposed seems quite reasonable to me and would also make it easy for user to patch if they want a new structure.

For what it's worth it seems that pybind11 allows for wrapping functors:
https://pybind11.readthedocs.io/en/stable/advanced/cast/functional.html
But since dijkstra3d is not using pybind11, since it may incur a performance impact and since it may be enough to pass parameters and construct the functor in c++, there may not be a strong benefit in supporting this.

Personnally I would be interested in the following function (edited for the typo mentioned below) where I is the image field, x_s is the source voxel of the edge and x_t the target:

x_s - x_t would ideally take into account voxel spacing / anisotropy. alpha and beta are parameters.

[william-silversmith]

That pybind11 functionality is very cool. Thanks for the tip! I haven't used pybind11 much because I have felt more comfortable with using numpy (and other existing python functionality) in Cython, but I suspect that there's no real technical advantage.

That function looks like a good candidate. In 3D, I assume the equation would be (where X is a 3-vector):

sqrt( alpha^2 * (Xs - Xt) + beta^2( I(Xs) - I(Xt))^2 )

I'm a little confused by the missing squaring of (xs - xt) is that intentional? I think the vectorized version doesn't make sense without it.

[tvercaut]

Yes, X is a 3D vector in 3D, 2D in 2D.

I fixed the missing square typo. Thanks for spotting it

[william-silversmith]

Sorry was on vacation, still thinking about the right way to do this that won't clutter up the code or reduce performance.

[william-silversmith]

Thinking about this some more, it looks like if we don't want to blow away some careful optimization, we'll need to create a second function. The reason for this is neighborhood doesn't give easy access to the new x,y, and z indicies of the neighbor and so must be expensively recalculated. The metric is also more expensive to compute, so making sure that it gets inlined is probably important.

[william-silversmith]

I have a prototype of a C++ function that could be helpful. It only supports 26-connectivity though. If you find that helpful, I can add a python wrapper for it.

[tvercaut]

Sounds great. I'm on the go at the moment but @ReubenDo in my group has been looking into our specific use case and created a "hardcoded" version for it here if it is of interest:
https://github.com/ReubenDo/dijkstra3d

[william-silversmith]

That's really helpful to see! Thank you! I noticed that the equation is modified in the repo referenced. I doubt that the pseudo-probability that a voxel is background is the kind of detail that would make sense in a library implementation though.

I also saw that time trials were run in that repo of dijkstra's performance. Have you tried modifying the bidirectional search? In certain cases it can be much faster.