feat: add dijkstra_anisotropy to support anisotropy in 3d dijkstra field

[turtleizzy]

Test:

import matplotlib.pyplot as plt
import numpy as np
import cv2
import dijkstra3d

def arclen(pnts):
    return np.sum(np.linalg.norm(np.diff(pnts.astype('float'), axis=0), axis=1))

weight = np.ones([1, 50, 50])
stPnt = [0, 0, 0]
edPnt = [0, 30, 49]
path = dijkstra3d.dijkstra(weight, stPnt, edPnt, anisotropy=(1, 1, 1), connectivity=26)

canvas = np.zeros([1, 50, 50], dtype=np.uint8)
cv2.line(canvas[0], stPnt[1:][::-1], edPnt[1:][::-1], color=2)
canvas[path[:, 0], path[:, 1], path[:, 2]] += 1

plt.imshow(canvas[0])

print(f"Euclidean length of dijkstra path: {arclen(path)}")
print(f"Euclidean length of bresenham line: {arclen(np.stack(np.where(canvas & 2 > 0)).T)}")

path = dijkstra3d.dijkstra(weight, stPnt, edPnt, connectivity=26)

canvas = np.zeros([1, 50, 50], dtype=np.uint8)
cv2.line(canvas[0], stPnt[1:][::-1], edPnt[1:][::-1], color=2)
canvas[path[:, 0], path[:, 1], path[:, 2]] += 1

plt.imshow(canvas[0])

print(f"Euclidean length of dijkstra path: {arclen(path)}")
print(f"Euclidean length of bresenham line: {arclen(np.stack(np.where(canvas & 2 > 0)).T)}")

Output:

Euclidean length of dijkstra path: 61.42640687119285
Euclidean length of bresenham line: 61.426406871192846
Euclidean length of dijkstra path: 66.39696961966999
Euclidean length of bresenham line: 61.426406871192846

[william-silversmith]

Hi turtleizzy! Thanks for the contribution! I think this could make a lot of sense, but need to think about the implications a bit more and maybe run some tests for non-uniform fields. I can tell you that a change like this makes a lot of sense for binary_dijkstra. I'm just a bit more sensitive about the regular dijkstra function because it is used with (and was designed for) https://github.com/seung-lab/kimimaro . I am wondering if the centerlines generated with this metric will be nicer, which is quite possible as they also ultimately exist in a euclidean space (but dijkstra is applied to a penalty field derived from the euclidean distance transform).

[turtleizzy]

I am wondering if the centerlines generated with this metric will be nicer, which is quite possible as they also ultimately exist in a euclidean space (but dijkstra is applied to a penalty field derived from the euclidean distance transform).

I was also using dijkstra3d for this use case (dijkstra on 3d edt) and I agree that the implication should be evaluated more carefully.

Currently in my implementation, I kept dijkstra3d unchanged and added a distinct dijkstra3d_anisotropy function for this feature to keep backward compatibility in C++. Also at the Cython level, I made anisotropy parameter optional, and that the behavior will be identical to the previous dijkstra3d.dijkstra if not provided, which should be backward compatible since this parameter was not supported previously.

I expect that this feature should make a difference in the 'flat' regions in the weight field, but should not influence much in more 'fluctuating' regions. For a weight field derived from the distance map, most neighboring voxels have different values, especially after some kind of augmentation (for example, pdrf_scale and pdrf_exponent in TEASER). So I expect that the result might not differ much.

Also on the other hand, even after adding the anisotropy feature in dijkstra3d, the output is still significantly different from the more apparent straight line generated by the Bresenham algorithm, even if two path have identical euclidean arclen. I personally don't know an algorithm to generate Bresenham-like line in shortest path problem on grids. More research may be necessary for a more beautiful solution.

[william-silversmith]

Will take another look soon! So sorry, just been a little busy at the moment.

[william-silversmith]

Okay, one last thing! Would you mind adding a few basic tests to make sure it gives a reasonable result? For example, that the path contacts start/end and/or that its a shorter path than dijkstra would give you?

[turtleizzy]

Okay, one last thing! Would you mind adding a few basic tests to make sure it gives a reasonable result? For example, that the path contacts start/end and/or that its a shorter path than dijkstra would give you?

Sure, will do it soon!

[turtleizzy]

Okay, one last thing! Would you mind adding a few basic tests to make sure it gives a reasonable result? For example, that the path contacts start/end and/or that its a shorter path than dijkstra would give you?

Sure, will do it soon!

Ahhh sooorryy I’ve been quite busy that days and I completely forgot this issue >_< Thank you so much @william-silversmith