From 72bbc30af686def919a97b45b259cf32f747856d Mon Sep 17 00:00:00 2001
From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com>
Date: Mon, 22 Aug 2022 12:49:22 +0200
Subject: [PATCH] Add all_closest to place_pwlin. (#1952)

---
 arbor/include/arbor/morph/place_pwlin.hpp |  5 +++-
 arbor/morph/place_pwlin.cpp               | 35 ++++++++++++++++-------
 2 files changed, 28 insertions(+), 12 deletions(-)

diff --git a/arbor/include/arbor/morph/place_pwlin.hpp b/arbor/include/arbor/morph/place_pwlin.hpp
index 73c0d5e0..09cc179d 100644
--- a/arbor/include/arbor/morph/place_pwlin.hpp
+++ b/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_;
 };
diff --git a/arbor/morph/place_pwlin.cpp b/arbor/morph/place_pwlin.cpp
index b5492317..6b55d26b 100644
--- a/arbor/morph/place_pwlin.cpp
+++ b/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
-- 
GitLab