perf: euclidean distance field acceleration

README.md

@@ -28,6 +28,13 @@ print(path.shape)
 # the anisotropic euclidean distance from the source to all labeled vertices.
 dist_field = dijkstra3d.euclidean_distance_field(field, source=(0,0,0), 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
+# is unobstructed space. For that region, we use an equation instead
+# of dijkstra's algorithm. Hybrid algorithm! free_space_radius is a physical
+# distance, meaning you must account for anisotropy in setting it.
+dist_field = dijkstra3d.euclidean_distance_field(field, source=(0,0,0), anisotropy=(4,4,40), free_space_radius=300) 
+
 # 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.
@@ -82,7 +89,7 @@ float* field = dijkstra::euclidean_distance_field3d<float>(
   labels, 
   /*sx=*/512, /*sy=*/512, /*sz=*/512, 
   /*wx=*/4, /*wy=*/4, /*wz=*/40, 
-  source
+  source, /*free_space_radius=*/0 // set to > 0 to switch on
 );
 
 float* field = dijkstra::distance_field3d<float>(labels, /*sx=*/512, /*sy=*/512, /*sz=*/512, source);

automated_test.py

@@ -563,7 +563,8 @@ def test_distance_field_2d_asymmetric(dtype):
   field = dijkstra3d.distance_field(values, (0,1))
   assert np.all(field == answer)
 
-def test_euclidean_distance_field_2d():
+@pytest.mark.parametrize('free_space_radius', (0,1,2,3,4,5,10))
+def test_euclidean_distance_field_2d(free_space_radius):
   values = np.ones((2, 2), dtype=bool)
 
   sq2 = sqrt(2)
@@ -574,7 +575,7 @@ def test_euclidean_distance_field_2d():
     [1, sq2],
   ], dtype=np.float32)
 
-  field = dijkstra3d.euclidean_distance_field(values, (0,0))
+  field = dijkstra3d.euclidean_distance_field(values, (0,0), free_space_radius=free_space_radius)
   assert np.all(np.abs(field - answer) < 0.00001)
 
   values = np.ones((5, 5), dtype=bool)
@@ -587,7 +588,7 @@ def test_euclidean_distance_field_2d():
      [4,        4.4142137 , 4.8284273, 5.2426405 , 5.656854 ]], 
   dtype=np.float32)
 
-  field = dijkstra3d.euclidean_distance_field(values, (0,0,0), (1,1,1))
+  field = dijkstra3d.euclidean_distance_field(values, (0,0,0), (1,1,1), free_space_radius=free_space_radius)
   assert np.all(np.abs(field - answer) < 0.00001) 
 
   answer = np.array([
@@ -602,11 +603,11 @@ def test_euclidean_distance_field_2d():
   ], dtype=np.float32)
 
   values = np.ones((2, 2, 2), dtype=bool)
-  field = dijkstra3d.euclidean_distance_field(values, (0,0,0), (1,1,1))
+  field = dijkstra3d.euclidean_distance_field(values, (0,0,0), (1,1,1), free_space_radius=free_space_radius)
   assert np.all(np.abs(field - answer) < 0.00001)   
 
   values = np.ones((2, 2, 2), dtype=bool)
-  field = dijkstra3d.euclidean_distance_field(values, (1,1,1), (1,1,1))
+  field = dijkstra3d.euclidean_distance_field(values, (1,1,1), (1,1,1), free_space_radius=free_space_radius)
 
   answer = np.array([
     [
@@ -621,6 +622,15 @@ def test_euclidean_distance_field_2d():
 
   assert np.all(np.abs(field - answer) < 0.00001)   
 
+@pytest.mark.parametrize('point', (np.random.randint(0,256, size=(3,)),))
+def test_euclidean_distance_field_3d_free_space_eqn(point):
+  point = tuple(point)
+  print(point)
+  values = np.ones((256, 256, 256), dtype=bool)
+  field_dijk = dijkstra3d.euclidean_distance_field(values, point, free_space_radius=0) # free space off
+  field_free = dijkstra3d.euclidean_distance_field(values, point, free_space_radius=10000) # free space 100% on
+
+  assert np.all(np.abs(field_free - field_dijk) < 0.001) # there's some difference below this
 
 def test_compass():
   field = np.array([