Add all_closest to place_pwlin.

arbor/include/arbor/morph/place_pwlin.hpp

@@ -86,9 +86,12 @@ struct ARB_ARBOR_API place_pwlin {
     // Maximal set of segments or part segments whose union is coterminous with extent.
     std::vector<msegment> all_segments(const mextent& extent) const;
 
-    // The closest location to p. Returns the location and its distance from the input coordinates.
+    // The closest location to p. Returns the location and its distance from the input coordinates. Ties are broken in favour of the most proximal point
     std::pair<mlocation, double> closest(double x, double y, double z) const;
 
+    // The closest location to p. Returns all possible locations and their shared distance from the input coordinates.
+    std::pair<std::vector<mlocation>, double> all_closest(double x, double y, double z) const;
+
 private:
     std::shared_ptr<place_pwlin_data> data_;
 };

arbor/include/arbor/morph/primitives.hpp

@@ -122,3 +122,4 @@ ARB_ARBOR_API bool test_invariants(const mcable_list&);
 
 ARB_DEFINE_HASH(arb::mcable, a.branch, a.prox_pos, a.dist_pos);
 ARB_DEFINE_HASH(arb::mlocation, a.branch, a.pos);
+ARB_DEFINE_HASH(arb::mpoint, a.x, a.y, a.z, a.radius);

arbor/morph/place_pwlin.cpp

@@ -199,16 +199,10 @@ struct p3d {
     }
 };
 
-// Policy:
-//  If two collated points are equidistant from the input point, take the
-//  proximal location.
-// Rationale:
-//  if the location is on a fork point, it makes sense to take the proximal
-//  location, which corresponds to the end of the parent branch.
-std::pair<mlocation, double> place_pwlin::closest(double x, double y, double z) const {
+std::pair<std::vector<mlocation>, double> place_pwlin::all_closest(double x, double y, double z) const {
     double mind = std::numeric_limits<double>::max();
     p3d p(x,y,z);
-    mlocation loc;
+    std::vector<mlocation> locs;
 
     // loop over each branch
     for (msize_t bid: util::count_along(data_->segment_index)) {
@@ -224,9 +218,13 @@ std::pair<mlocation, double> place_pwlin::closest(double x, double y, double z)
             const double wvs = dot(vw, vw);
             if (wvs==0.) { // zero length segment is a special case
                 const double distance = norm(p-v);
+                mlocation loc{bid, s.lower_bound()};
                 if (distance<mind) {
                     mind = distance;
-                    loc = {bid, s.lower_bound()};
+                    locs = {loc};
+                }
+                else if (distance == mind) {
+                    locs.push_back(loc);
                 }
             }
             else {
@@ -240,14 +238,29 @@ std::pair<mlocation, double> place_pwlin::closest(double x, double y, double z)
                     t<=0.? norm(p-v):
                     t>=1.? norm(p-w):
                            norm(p-(v + t*vw));
+                mlocation loc{bid, math::lerp(s.lower_bound(), s.upper_bound(), t)};
                 if (distance<mind) {
-                    loc = {bid, math::lerp(s.lower_bound(), s.upper_bound(), t)};
+                    locs = {loc};
                     mind = distance;
                 }
+                else if (distance == mind) {
+                    locs.push_back(loc);
+                }
             }
         }
     }
-    return {loc, mind};
+    return {locs, mind};
+}
+
+// Policy:
+//  If two collated points are equidistant from the input point, take the
+//  proximal location.
+// Rationale:
+//  if the location is on a fork point, it makes sense to take the proximal
+//  location, which corresponds to the end of the parent branch.
+std::pair<mlocation, double> place_pwlin::closest(double x, double y, double z) const {
+    const auto& [locs, delta] = all_closest(x, y, z);
+    return {locs.front(), delta};
 }
 
 } // namespace arb