diff --git a/README.md b/README.md
index 470e581..908a096 100644
--- a/README.md
+++ b/README.md
@@ -45,9 +45,16 @@ print(path.shape)
 # the anisotropic euclidean chamfer distance from the source to all labeled vertices.
 # Source can be a single point or a list of points. Accepts bool, (u)int8 dtypes.
 dist_field = dijkstra3d.euclidean_distance_field(field, source=(0,0,0), anisotropy=(4,4,40))
+
+sources = [ (0,0,0), (10, 40, 232) ]
 dist_field = dijkstra3d.euclidean_distance_field(
-  field, source=[ (0,0,0), (10, 40, 232) ], anisotropy=(4,4,40)
+  field, source=sources, anisotropy=(4,4,40)
 )
+# You can return a map of source vertices to nearest voxels called
+# a feature map.
+dist_field, feature_map = dijkstra3d.euclidean_distance_field(
+  field, source=sources, return_feature_map=True,
+) 
 
 # 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
diff --git a/automated_test.py b/automated_test.py
index f107907..d65349d 100644
--- a/automated_test.py
+++ b/automated_test.py
@@ -903,6 +903,52 @@ def test_euclidean_distance_field_2d(free_space_radius):
   ])
   assert np.all(np.isclose(field, answer))
 
+def test_euclidean_distance_field_2d_feature_map():
+  values = np.ones((7, 7, 1), dtype=bool)
+  field, max_loc, feature_map = dijkstra3d.euclidean_distance_field(values, [0,0,0], return_max_location=True, return_feature_map=True)
+  
+  assert np.isclose(
+    np.max(field), 
+    np.sqrt((values.shape[0] - 1) ** 2 + (values.shape[1] - 1) ** 2)
+  )
+  assert max_loc == (6,6,0)
+
+  assert np.all(np.unique(feature_map) == 1)
+
+  sources = [
+    [0,0,0],
+    [6,6,0]
+  ]
+  field, max_loc, feature_map = dijkstra3d.euclidean_distance_field(values, sources, return_max_location=True, return_feature_map=True)
+  assert tuple(np.unique(feature_map)) == (1,2)
+
+  sources = [
+    [0,0,0],
+    [6,6,0],
+    [6,6,0]
+  ]
+  field, max_loc, feature_map = dijkstra3d.euclidean_distance_field(values, sources, return_max_location=True, return_feature_map=True)
+  assert tuple(np.unique(feature_map)) == (1,2)
+
+  sources = [
+    [0,0,0],
+    [6,6,0],
+    [3,3,0],
+  ]
+  field, max_loc, feature_map = dijkstra3d.euclidean_distance_field(values, sources, return_max_location=True, return_feature_map=True)
+  assert tuple(np.unique(feature_map)) == (1,2,3)
+  
+  result = np.array(
+    [[[1, 1, 1, 2, 2, 2, 2],
+      [1, 1, 2, 2, 2, 2, 2],
+      [1, 2, 2, 2, 2, 2, 2],
+      [2, 2, 2, 2, 2, 2, 3],
+      [2, 2, 2, 2, 2, 3, 3],
+      [2, 2, 2, 2, 3, 3, 3],
+      [2, 2, 2, 3, 3, 3, 3],]]).T
+
+  assert np.all(feature_map == result)
+
 @pytest.mark.parametrize('point', (np.random.randint(0,256, size=(3,)),))
 def test_euclidean_distance_field_3d_free_space_eqn(point):
   point = tuple(point)
diff --git a/dijkstra3d.hpp b/dijkstra3d.hpp
index a88534a..c3e5018 100644
--- a/dijkstra3d.hpp
+++ b/dijkstra3d.hpp
@@ -1835,6 +1835,184 @@ float* euclidean_distance_field3d(
   );
 }
 
+
+class HeapFeatureNode {
+public:
+  float key;
+  uint64_t value;
+  uint64_t source;
+
+  HeapFeatureNode() {
+    key = 0;
+    value = 0;
+  }
+
+  HeapFeatureNode (float k, uint64_t val, uint64_t src) {
+    key = k;
+    value = val;
+    source = src;
+  }
+
+  HeapFeatureNode (const HeapFeatureNode &h) {
+    key = h.key;
+    value = h.value;
+    source = h.source;
+  }
+};
+
+struct HeapFeatureNodeCompare {
+  bool operator()(const HeapFeatureNode &t1, const HeapFeatureNode &t2) const {
+    return t1.key >= t2.key;
+  }
+};
+
+// returns a map of the nearest source vertex
+template <typename OUT = uint64_t>
+OUT* edf_with_feature_map(
+  uint8_t* field, // really a boolean field
+  const uint64_t sx, const uint64_t sy, const uint64_t sz, 
+  const float wx, const float wy, const float wz, 
+  const std::vector<uint64_t> &sources,
+  float* dist = NULL, OUT* feature_map = NULL,
+  size_t &max_loc = _dummy_max_loc
+) {
+
+  const uint64_t voxels = sx * sy * sz;
+  const uint64_t sxy = sx * sy;
+
+  const libdivide::divider<uint64_t> fast_sx(sx); 
+  const libdivide::divider<uint64_t> fast_sxy(sxy); 
+
+  const bool power_of_two = !((sx & (sx - 1)) || (sy & (sy - 1))); 
+  const int xshift = std::log2(sx); // must use log2 here, not lg/lg2 to avoid fp errors
+  const int yshift = std::log2(sy);
+
+  bool clear_dist = false;
+  if (dist == NULL) {
+    dist = new float[voxels]();
+    clear_dist = true;
+  }
+  if (feature_map == NULL) {
+    feature_map = new OUT[voxels]();
+  }
+
+  fill(dist, +INFINITY, voxels);
+
+  int neighborhood[NHOOD_SIZE] = {};
+
+  float neighbor_multiplier[NHOOD_SIZE] = { 
+    wx, wx, wy, wy, wz, wz, // axial directions (6)
+    
+    // square diagonals (12)
+    _s(wx, wy), _s(wx, wy), _s(wx, wy), _s(wx, wy),  
+    _s(wy, wz), _s(wy, wz), _s(wy, wz), _s(wy, wz),
+    _s(wx, wz), _s(wx, wz), _s(wx, wz), _s(wx, wz),
+
+    // cube diagonals (8)
+    _c(wx, wy, wz), _c(wx, wy, wz), _c(wx, wy, wz), _c(wx, wy, wz), 
+    _c(wx, wy, wz), _c(wx, wy, wz), _c(wx, wy, wz), _c(wx, wy, wz)
+  };
+
+  std::priority_queue<
+    HeapFeatureNode, std::vector<HeapFeatureNode>, HeapFeatureNodeCompare
+  > queue;
+
+  uint64_t src = 1;
+  for (uint64_t source : sources) {
+    dist[source] = -0;
+    feature_map[source] = src;
+    queue.emplace(0.0, source, src);
+    src++;
+  }
+
+  uint64_t loc, next_loc;
+  float new_dist;
+  uint64_t neighboridx;
+
+  uint64_t x, y, z;
+
+  float max_dist = -1;
+  max_loc = voxels + 1;
+
+  while (!queue.empty()) {
+    src = queue.top().source;
+    loc = queue.top().value;
+    queue.pop();
+
+    if (max_dist < std::abs(dist[loc])) {
+      max_dist = std::abs(dist[loc]);
+      max_loc = loc;
+    }
+
+    if (std::signbit(dist[loc])) {
+      continue;
+    }
+
+    if (!queue.empty()) {
+      next_loc = queue.top().value;
+      if (!std::signbit(dist[next_loc])) {
+
+        // As early as possible, start fetching the
+        // data from RAM b/c the annotated lines below
+        // have 30-50% cache miss.
+        DIJKSTRA_3D_PREFETCH_26WAY(field, next_loc)
+        DIJKSTRA_3D_PREFETCH_26WAY(dist, next_loc)
+      }
+    }
+
+    if (power_of_two) {
+      z = loc >> (xshift + yshift);
+      y = (loc - (z << (xshift + yshift))) >> xshift;
+      x = loc - ((y + (z << yshift)) << xshift);
+    }
+    else {
+      z = loc / fast_sxy;
+      y = (loc - (z * sxy)) / fast_sx;
+      x = loc - sx * (y + z * sy);
+    }
+
+    compute_neighborhood(neighborhood, x, y, z, sx, sy, sz, NHOOD_SIZE, NULL);
+
+    for (int i = 0; i < NHOOD_SIZE; i++) {
+      if (neighborhood[i] == 0) {
+        continue;
+      }
+
+      neighboridx = loc + neighborhood[i];
+      if (field[neighboridx] == 0) {
+        continue;
+      }
+
+      new_dist = dist[loc] + neighbor_multiplier[i];
+      
+      // Visited nodes are negative and thus the current node
+      // will always be less than as field is filled with non-negative
+      // integers.
+      if (new_dist < dist[neighboridx]) { 
+        dist[neighboridx] = new_dist;
+        feature_map[neighboridx] = src;
+        queue.emplace(new_dist, neighboridx, src);
+      }
+      else if (new_dist == dist[neighboridx] && src > feature_map[neighboridx]) {
+        feature_map[neighboridx] = src;
+      }
+    }
+
+    dist[loc] = -dist[loc];
+  }
+
+  if (clear_dist) {
+    delete[] dist;
+  }
+  else {
+    for (size_t i = 0; i < voxels; i++) {
+      dist[i] = std::fabs(dist[i]);
+    }
+  }
+
+  return feature_map;
+}
+
 #undef sq
 #undef NHOOD_SIZE
 #undef DIJKSTRA_3D_PREFETCH_26WAY
diff --git a/dijkstra3d.pyx b/dijkstra3d.pyx
index f450576..2d6279d 100644
--- a/dijkstra3d.pyx
+++ b/dijkstra3d.pyx
@@ -124,6 +124,14 @@ cdef extern from "dijkstra3d.hpp" namespace "dijkstra":
     uint32_t* voxel_graph,
     size_t &max_loc
   )
+  cdef OUT* edf_with_feature_map[OUT](
+    uint8_t* field,
+    uint64_t sx, uint64_t sy, uint64_t sz, 
+    float wx, float wy, float wz, 
+    vector[uint64_t] &sources,
+    float* dist, OUT* feature_map,
+    size_t &max_loc
+  )
   cdef vector[T] query_shortest_path[T](
     T* parents, T target
   ) 
@@ -600,7 +608,8 @@ def parental_field(data, source, connectivity=26, voxel_graph=None):
 def euclidean_distance_field(
   data, source, anisotropy=(1,1,1), 
   free_space_radius=0, voxel_graph=None,
-  return_max_location=False
+  return_max_location=False,
+  return_feature_map=False,
 ):
   """
   Use dijkstra's shortest path algorithm
@@ -617,7 +626,7 @@ def euclidean_distance_field(
   negative zero).
   
   Parameters:
-   data: Input weights in a 2D or 3D numpy array. 
+   data: binary image represented as a 2D or 3D numpy array. 
    source: (x,y,z) coordinate or list of coordinates 
     of starting voxels
    anisotropy: (wx,wy,wz) weights for each axial direction.
@@ -625,15 +634,28 @@ def euclidean_distance_field(
     region surrounding the source is free space, we can use
     a much faster algorithm to fill in that volume. Value
     is physical radius (can get this from the EDT). 
+   voxel_graph: a bitwise representation of the premitted
+    directions of travel between voxels. Generated from
+    cc3d.voxel_connectivity_graph. 
+    (See https://github.com/seung-lab/connected-components-3d)
    return_max_location: returns the coordinates of one
      of the possibly multiple maxima.
+   return_feature_map: return labeling of the image such
+    that all the voxels nearest a given source point are 
+    a single label. At locations where two labels are 
+    equidistant, the numerically larger label wins (to
+    differences in behavior on different platforms).
   
   Returns: 
     let field = 2D or 3D numpy array with each index
       containing its distance from the source voxel.
 
-    if return_max_location:
+    if return_max_location and return_feature_map:
+      return (field, (x,y,z) of max distance, feature_map)
+    elif return_max_location:
       return (field, (x,y,z) of max distance)
+    elif return_feature_map:
+      return (field, feature_map)
     else:
       return field
   """
@@ -664,19 +686,38 @@ def euclidean_distance_field(
 
   data = np.asfortranarray(data)
 
-  field, max_loc = _execute_euclidean_distance_field(
-    data, source, anisotropy, 
-    free_space_radius, voxel_graph
-  )
+  if return_feature_map:
+    field, feature_map, max_loc = _execute_euclidean_distance_field_w_feature_map(
+      data, source, anisotropy,
+    )
+  else:
+    feature_map = None
+    field, max_loc = _execute_euclidean_distance_field(
+      data, source, anisotropy, 
+      free_space_radius, voxel_graph
+    )
+
   if dims < 3:
     field = np.squeeze(field, axis=2)
+    if feature_map:
+      feature_map = np.squeeze(feature_map, axis=2)
   if dims < 2:
     field = np.squeeze(field, axis=1)
+    if feature_map:
+      feature_map = np.squeeze(feature_map, axis=1)
+
+  ret = [ field ]
 
   if return_max_location:
-    return field, np.unravel_index(max_loc, data.shape, order="F")[:dims]
-  else:
-    return field
+    ret.append(
+      np.unravel_index(max_loc, data.shape, order="F")[:dims]
+    )
+  if return_feature_map:
+    ret.append(feature_map)
+  
+  if len(ret) == 1:
+    return ret[0]
+  return tuple(ret)
 
 def _validate_coord(data, coord):
   dims = len(data.shape)
@@ -1846,3 +1887,48 @@ def _execute_euclidean_distance_field(
     raise ValueError(f"Something went wrong during processing. max_loc: {max_loc}")
 
   return dist, max_loc
+
+def _execute_euclidean_distance_field_w_feature_map(
+  data, sources, anisotropy
+):
+  cdef uint8_t[:,:,:] arr_memview8
+
+  cdef uint64_t sx = data.shape[0]
+  cdef uint64_t sy = data.shape[1]
+  cdef uint64_t sz = data.shape[2]
+
+  cdef float wx = anisotropy[0]
+  cdef float wy = anisotropy[1]
+  cdef float wz = anisotropy[2]
+
+  sources = np.unique(sources, axis=0)
+
+  cdef vector[uint64_t] src
+  for source in sources:
+    src.push_back(source[0] + sx * (source[1] + sy * source[2]))
+
+  cdef cnp.ndarray[float, ndim=3] dist = np.zeros( (sx,sy,sz), dtype=np.float32, order='F' )
+  cdef cnp.ndarray[uint32_t, ndim=3] feature_map = np.zeros( (sx,sy,sz), dtype=np.uint32, order='F' )
+
+  dtype = data.dtype
+  cdef size_t max_loc = data.size + 1
+
+  if dtype in (np.int8, np.uint8, bool):
+    arr_memview8 = data.view(np.uint8)
+    edf_with_feature_map[uint32_t](
+      &arr_memview8[0,0,0],
+      sx, sy, sz,
+      wx, wy, wz,
+      src, 
+      &dist[0,0,0],
+      &feature_map[0,0,0],
+      max_loc
+    )
+  else:
+    raise TypeError(f"Type {dtype} not currently supported.")
+
+  if max_loc == data.size + 1:
+    raise ValueError(f"Something went wrong during processing. max_loc: {max_loc}")
+
+  return dist, feature_map, max_loc
+