diff --git a/arbor/include/arbor/morph/place_pwlin.hpp b/arbor/include/arbor/morph/place_pwlin.hpp index 73c0d5e0b5a532ee488d1e4fabc4de5c9c9d6ace..09cc179dc72537d0af1b77900595f767e45a84f4 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 b5492317b8d0944499068def431476438a34d574..6b55d26b81b511eddfd1b081a5806b4abd1779c3 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