From a64ae6db9b5a17890bd524f8d5959061c9acf446 Mon Sep 17 00:00:00 2001 From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com> Date: Thu, 18 Aug 2022 10:37:13 +0200 Subject: [PATCH 1/2] Add all_closest to place_pwlin. --- arbor/include/arbor/morph/place_pwlin.hpp | 4 +- arbor/include/arbor/morph/primitives.hpp | 1 + arbor/morph/place_pwlin.cpp | 52 +++++++++++++++++++++++ 3 files changed, 56 insertions(+), 1 deletion(-) diff --git a/arbor/include/arbor/morph/place_pwlin.hpp b/arbor/include/arbor/morph/place_pwlin.hpp index 73c0d5e0b5..2e8eb19367 100644 --- a/arbor/include/arbor/morph/place_pwlin.hpp +++ b/arbor/include/arbor/morph/place_pwlin.hpp @@ -86,9 +86,11 @@ struct ARB_ARBOR_API place_pwlin { // Maximal set of segments or part segments whose union is coterminous with extent. std::vector 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 brokwn in favour of the most proximal point std::pair closest(double x, double y, double z) const; + std::pair, double> all_closest(double x, double y, double z) const; + private: std::shared_ptr data_; }; diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp index 6c5d64bfc5..cdb6472acb 100644 --- a/arbor/include/arbor/morph/primitives.hpp +++ b/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); diff --git a/arbor/morph/place_pwlin.cpp b/arbor/morph/place_pwlin.cpp index b5492317b8..f86fb66c89 100644 --- a/arbor/morph/place_pwlin.cpp +++ b/arbor/morph/place_pwlin.cpp @@ -250,4 +250,56 @@ std::pair place_pwlin::closest(double x, double y, double z) return {loc, mind}; } +std::pair, double> place_pwlin::all_closest(double x, double y, double z) const { + double mind = std::numeric_limits::max(); + p3d p(x,y,z); + std::vector locs; + + // loop over each branch + for (msize_t bid: util::count_along(data_->segment_index)) { + const auto b = data_->segment_index[bid]; + // loop over the segments in the branch + for (auto s: b) { + const auto& seg = data_->segments[s.value]; + + // v and w are the proximal and distal ends of the segment. + const p3d v = seg.prox; + const p3d w = seg.dist; + const p3d vw = w-v; + 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 proximal end of the segment + // t=1 -> distal end of the segment + // values are clamped to the range [0, 1] + const double t = std::max(0., std::min(1., dot(vw, p-v) / wvs)); + const double distance = + 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 Date: Mon, 22 Aug 2022 11:26:08 +0200 Subject: [PATCH 2/2] Review. --- arbor/include/arbor/morph/place_pwlin.hpp | 3 +- arbor/morph/place_pwlin.cpp | 63 +++++------------------ 2 files changed, 14 insertions(+), 52 deletions(-) diff --git a/arbor/include/arbor/morph/place_pwlin.hpp b/arbor/include/arbor/morph/place_pwlin.hpp index 2e8eb19367..09cc179dc7 100644 --- a/arbor/include/arbor/morph/place_pwlin.hpp +++ b/arbor/include/arbor/morph/place_pwlin.hpp @@ -86,9 +86,10 @@ struct ARB_ARBOR_API place_pwlin { // Maximal set of segments or part segments whose union is coterminous with extent. std::vector all_segments(const mextent& extent) const; - // The closest location to p. Returns the location and its distance from the input coordinates. Ties are brokwn in favour of the most proximal point + // 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 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, double> all_closest(double x, double y, double z) const; private: diff --git a/arbor/morph/place_pwlin.cpp b/arbor/morph/place_pwlin.cpp index f86fb66c89..6b55d26b81 100644 --- a/arbor/morph/place_pwlin.cpp +++ b/arbor/morph/place_pwlin.cpp @@ -199,57 +199,6 @@ 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 place_pwlin::closest(double x, double y, double z) const { - double mind = std::numeric_limits::max(); - p3d p(x,y,z); - mlocation loc; - - // loop over each branch - for (msize_t bid: util::count_along(data_->segment_index)) { - const auto b = data_->segment_index[bid]; - // loop over the segments in the branch - for (auto s: b) { - const auto& seg = data_->segments[s.value]; - - // v and w are the proximal and distal ends of the segment. - const p3d v = seg.prox; - const p3d w = seg.dist; - const p3d vw = w-v; - const double wvs = dot(vw, vw); - if (wvs==0.) { // zero length segment is a special case - const double distance = norm(p-v); - if (distance proximal end of the segment - // t=1 -> distal end of the segment - // values are clamped to the range [0, 1] - const double t = std::max(0., std::min(1., dot(vw, p-v) / wvs)); - const double distance = - t<=0.? norm(p-v): - t>=1.? norm(p-w): - norm(p-(v + t*vw)); - if (distance, double> place_pwlin::all_closest(double x, double y, double z) const { double mind = std::numeric_limits::max(); p3d p(x,y,z); @@ -302,4 +251,16 @@ std::pair, double> place_pwlin::all_closest(double x, dou } 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 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