feat: multi-source distance field

README.md

@@ -33,8 +33,12 @@ path = dijkstra3d.path_from_parents(parents, target=(511, 511, 511))
 print(path.shape)
 
 # Given a boolean label "field" and a source vertex, compute 
-# the anisotropic euclidean distance from the source to all labeled vertices.
+# the anisotropic euclidean chamfer distance from the source to all labeled vertices.
+# Source can be a single point or a list of points.
 dist_field = dijkstra3d.euclidean_distance_field(field, source=(0,0,0), anisotropy=(4,4,40))
+dist_field = dijkstra3d.euclidean_distance_field(
+  field, source=[ (0,0,0), (10, 40, 232) ], anisotropy=(4,4,40)
+)
 
 # To make the EDF go faster add the free_space_radius parameter. It's only
 # safe to use if you know that some distance around the source point
@@ -45,8 +49,9 @@ dist_field = dijkstra3d.euclidean_distance_field(field, source=(0,0,0), anisotro
 
 # Given a numerical field, for each directed edge from adjacent voxels A and B, 
 # use B as the edge weight. In this fashion, compute the distance from a source 
-# point for all finite voxels.
-dist_field = dijkstra3d.distance_field(field, source=(0,0,0))
+# point for all finite voxels. 
+dist_field = dijkstra3d.distance_field(field, source=(0,0,0)) # single source
+dist_field = dijkstra3d.distance_field(field, source=[ (0,0,0), (52, 55, 23) ]) # multi-source
 
 # You can also provide a voxel connectivity graph to provide customized
 # constraints on the permissible directions of travel. The graph is a
@@ -109,13 +114,16 @@ uint32_t* parents = dijkstra::parental_field3d<float>(
 vector<unsigned int> path = dijkstra::query_shortest_path(parents, target);
 
 
+// Really a chamfer distance.
+// source can be a size_t (single source) or a std::vector<size_t> (multi-source)
 float* field = dijkstra::euclidean_distance_field3d<float>(
   labels, 
   /*sx=*/512, /*sy=*/512, /*sz=*/512, 
   /*wx=*/4, /*wy=*/4, /*wz=*/40, 
   source, /*free_space_radius=*/0 // set to > 0 to switch on
 );
 
+// source can be a size_t (single source) or a std::vector<size_t> (multi-source)
 float* field = dijkstra::distance_field3d<float>(labels, /*sx=*/512, /*sy=*/512, /*sz=*/512, source);
 ```
 

automated_test.py

@@ -463,6 +463,30 @@ def test_distance_field_2d(dtype):
     ]
   ]))
 
+  field = dijkstra3d.distance_field(values, [ (0,0), (4,4) ])
+
+  assert np.all(field == np.array([
+    [
+      [0, 1, 2, 3, 4],
+      [1, 1, 2, 3, 3],
+      [2, 2, 2, 2, 2],
+      [3, 3, 2, 1, 1],
+      [4, 3, 2, 1, 0],
+    ]
+  ]))
+
+  field = dijkstra3d.distance_field(values, [ (0,0), (2,2), (4,4) ])
+
+  assert np.all(field == np.array([
+    [
+      [0, 1, 2, 2, 2],
+      [1, 1, 1, 1, 2],
+      [2, 1, 0, 1, 2],
+      [2, 1, 1, 1, 1],
+      [2, 2, 2, 1, 0],
+    ]
+  ]))
+
   field = dijkstra3d.distance_field(values, (4,4))
 
   assert np.all(field == np.array([
@@ -500,6 +524,8 @@ def test_distance_field_2d(dtype):
     ]
   ]))
 
+
+
 @pytest.mark.parametrize("dtype", TEST_TYPES)
 def test_distance_field_2d_symmetric_26(dtype):
   values = np.ones((5, 5), dtype=dtype)
@@ -680,6 +706,18 @@ def test_euclidean_distance_field_2d(free_space_radius):
 
   assert np.all(np.abs(field - answer) < 0.00001)   
 
+  # Multi-source
+  values = np.ones((4,4), dtype=bool)
+  field = dijkstra3d.euclidean_distance_field(values, [ (0,0), (3,3) ], free_space_radius=free_space_radius)
+
+  answer = np.array([
+    [0, 1,        2,        3],
+    [1, sq2,      (1+sq2),  2],
+    [2, (1+sq2),  sq2,      1],
+    [3, 2,        1,        0],
+  ])
+  assert np.all(np.isclose(field, answer))
+
 @pytest.mark.parametrize('point', (np.random.randint(0,256, size=(3,)),))
 def test_euclidean_distance_field_3d_free_space_eqn(point):
   point = tuple(point)