diff --git a/arbor/include/arbor/morph/place_pwlin.hpp b/arbor/include/arbor/morph/place_pwlin.hpp
index da278fab7804ee74546da648dd34ff572d312fbc..223c287acbdf0b9f60bfc17bc0655aa0b6ed06e4 100644
--- a/arbor/include/arbor/morph/place_pwlin.hpp
+++ b/arbor/include/arbor/morph/place_pwlin.hpp
@@ -4,6 +4,7 @@
// sample points and interpolating linearly.
#include <cmath>
+#include <limits>
#include <utility>
#include <arbor/morph/morphology.hpp>
@@ -84,6 +85,9 @@ struct 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.
+ std::pair<mlocation, double> 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 9814aeb22c8f7fb01c0639a1b4fb9fdec3cadbe5..f42475bcfcdc38f6de857884a242ff734bc70228 100644
--- a/arbor/morph/place_pwlin.cpp
+++ b/arbor/morph/place_pwlin.cpp
@@ -1,10 +1,14 @@
+#include <iostream>
+
#include <cmath>
+#include <limits>
#include <memory>
#include <vector>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/place_pwlin.hpp>
#include <arbor/morph/primitives.hpp>
+#include <arbor/math.hpp>
#include "util/piecewise.hpp"
#include "util/rangeutil.hpp"
@@ -170,4 +174,82 @@ place_pwlin::place_pwlin(const arb::morphology& m, const isometry& iso) {
}
};
+struct p3d {
+ double x=0,y=0,z=0;
+ constexpr p3d() = default;
+ constexpr p3d(const mpoint& p): x(p.x), y(p.y), z(p.z) {}
+ constexpr p3d(double x, double y, double z): x(x), y(y), z(z) {}
+ friend constexpr p3d operator-(const p3d& l, const p3d& r) {
+ return {l.x-r.x, l.y-r.y, l.z-r.z};
+ }
+ friend constexpr p3d operator+(const p3d& l, const p3d& r) {
+ return {l.x+r.x, l.y+r.y, l.z+r.z};
+ }
+ friend constexpr p3d operator*(double l, const p3d& r) {
+ return {l*r.x, l*r.y, l*r.z};
+ }
+ friend constexpr double dot(const p3d& l, const p3d& r) {
+ return l.x*r.x + l.y*r.y + l.z*r.z;
+ }
+ friend double norm(const p3d p) {
+ return std::sqrt(dot(p, p));
+ }
+ friend std::ostream& operator<<(std::ostream& o, const p3d& p) {
+ return o << '[' << p.x << ' ' << p.y << ' ' << p.z << ']';
+ }
+};
+
+// 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 {
+ double mind = std::numeric_limits<double>::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 p3d vp = p-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<mind) {
+ mind = distance;
+ loc = {bid, s.lower_bound()};
+ }
+ }
+ else {
+ // Find the relative position of the orthogonal projection onto the line segment
+ // that along the axis of the segment:
+ // t=0 -> 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, vp) / wvs));
+ const double distance =
+ t<=0.? norm(p-v):
+ t>=1.? norm(p-w):
+ norm(p-(v + t*vw));
+ if (distance<mind) {
+ loc = {bid, math::lerp(s.lower_bound(), s.upper_bound(), t)};
+ mind = distance;
+ }
+ }
+ }
+ }
+ return {loc, mind};
+}
+
} // namespace arb
diff --git a/doc/cpp/morphology.rst b/doc/cpp/morphology.rst
index 581290e2e8b57ab6a7bc98c3ca349fff7115bef6..dc432b625b7c9b0d4998318aad2ba326dc9a72ff 100644
--- a/doc/cpp/morphology.rst
+++ b/doc/cpp/morphology.rst
@@ -183,6 +183,10 @@ zero-length components as a result of such discontinuities.
Return the maximal set of segments and partial segments whose
union is coterminous with the given :cpp:class:`mextent` in the placement.
+ .. cpp:function:: closest(double x, double y, double z) -> std::pair<mpoint, double>
+
+ Find the closest location to p. Returns the location and its distance from the input coordinates.
+
Isometries
^^^^^^^^^^
diff --git a/doc/python/morphology.rst b/doc/python/morphology.rst
index 972cf903b6a61e197ba2472f34ac1ac9c97d6338..6e233e183835029186b7255006789d8fd23112df 100644
--- a/doc/python/morphology.rst
+++ b/doc/python/morphology.rst
@@ -349,6 +349,10 @@ Cable cell morphology
union is coterminous with the sub-region of the morphology covered by
the given cables in the placement.
+ .. py:method:: closest(x: real, y: real, z: real) -> tuple[mpoint, real]
+
+ Find the closest location to p. Returns the location and its distance from the input coordinates.
+
.. py:class:: isometry
Isometries represent rotations and translations in space, and can be used with
diff --git a/python/morphology.cpp b/python/morphology.cpp
index 5a99baadd449c64c08214888ddd0049082321b53..a5d3fb0d15804616503c23b69c9109b67a5d19ca 100644
--- a/python/morphology.cpp
+++ b/python/morphology.cpp
@@ -199,7 +199,14 @@ void register_morphology(py::module& m) {
return self.all_segments(cables);
},
"Return maximal list of non-overlapping full or partial msegments whose union is coterminous "
- "with the extent of the given list of cables.");
+ "with the extent of the given list of cables.")
+ .def("closest",
+ [](const arb::place_pwlin& self, double x, double y, double z) {
+ auto [l, d] = self.closest(x, y, z);
+ return pybind11::make_tuple(l, d);
+ },
+ "Find the location on the morphology that is closest to a 3d point. "
+ "Returns the location and its distance from the point.");
//
// Higher-level data structures (segment_tree, morphology)
diff --git a/test/unit/test_morph_place.cpp b/test/unit/test_morph_place.cpp
index e7869fa705f5aca08871550b1e545e17ca2404b8..fff37b7461c7da3a956cdaba9d8409f903e6256e 100644
--- a/test/unit/test_morph_place.cpp
+++ b/test/unit/test_morph_place.cpp
@@ -570,3 +570,104 @@ TEST(place_pwlin, segments) {
EXPECT_TRUE(mpoint_almost_eq(p8p, x3all[2].prox));
EXPECT_TRUE(mpoint_almost_eq(p8p, x3all[2].dist));
}
+
+TEST(place_pwlin, nearest) {
+ segment_tree tree;
+
+ // the test morphology:
+ //
+ // x=-9 x=9
+ //
+ // _ _ y=40
+ // \ /
+ // seg4 \ / seg2
+ // branch 2 / / branch 1
+ // seg3 \ /
+ // | y=25
+ // |
+ // |
+ // | branch 0
+ // seg1 |
+ // |
+ // _ y=7
+ //
+ // - y=5
+ // seg0 |
+ // _ y=-5
+
+ // Root branch.
+ mpoint psoma_p{0, -5, 0, 5};
+ mpoint psoma_d{0, 5, 0, 5};
+
+ msize_t ssoma = tree.append(mnpos, psoma_p, psoma_d, 1);
+
+ // Main leg of y, of length 28 μm
+ // Note that there is a gap of 2 μm between the end of the soma segment
+ mpoint py1_p{0, 7, 0, 1};
+ mpoint py1_d{0, 25, 0, 1};
+
+ msize_t sy1 = tree.append(ssoma, py1_p, py1_d, 3);
+
+ // branch 1 of y: translation (9,15) in one segment
+ mpoint py2_d{ 9, 40, 0, 1};
+ tree.append(sy1, py2_d, 3);
+
+ // branch 2 of y: translation (-9,15) in 2 segments
+ mpoint py3_m{-6, 35, 0, 1};
+ mpoint py3_d{-9, 40, 0, 1};
+ tree.append(tree.append(sy1, py3_m, 3), py3_d, 3);
+
+ morphology m(tree);
+ place_pwlin place(m);
+
+ {
+ auto [l, d] = place.closest(0, -5, 0);
+ EXPECT_EQ((mlocation{0, 0.}), l);
+ EXPECT_EQ(0., d);
+ }
+ {
+ auto [l, d] = place.closest(10, -5, 0);
+ EXPECT_EQ((mlocation{0, 0.}), l);
+ EXPECT_EQ(10., d);
+ }
+ {
+ auto [l, d] = place.closest(0, 0, 0);
+ EXPECT_EQ((mlocation{0, 5./28.}), l);
+ EXPECT_EQ(0., d);
+ }
+ {
+ auto [l, d] = place.closest(10, 0, 0);
+ EXPECT_EQ((mlocation{0, 5./28.}), l);
+ EXPECT_EQ(10., d);
+ }
+ {
+ auto [l, d] = place.closest(0, 25, 0);
+ EXPECT_EQ((mlocation{0, 1.}), l);
+ EXPECT_EQ(0., d);
+ }
+ {
+ auto [l, d] = place.closest(0, 6, 0);
+ EXPECT_EQ((mlocation{0, 10./28.}), l);
+ EXPECT_EQ(1., d);
+ }
+ {
+ auto [l, d] = place.closest(3, 30, 0);
+ EXPECT_EQ((mlocation{1, 1./3.}), l);
+ EXPECT_EQ(0., d);
+ }
+ {
+ auto [l, d] = place.closest(-6, 35, 0);
+ EXPECT_EQ((mlocation{2, 2./3.}), l);
+ EXPECT_EQ(0., d);
+ }
+ {
+ auto [l, d] = place.closest(-9, 40, 0);
+ EXPECT_EQ((mlocation{2, 1.}), l);
+ EXPECT_EQ(0., d);
+ }
+ {
+ auto [l, d] = place.closest(-9, 41, 0);
+ EXPECT_EQ((mlocation{2, 1.}), l);
+ EXPECT_EQ(1., d);
+ }
+}