Skip to content
Snippets Groups Projects
Unverified Commit 0e8875a9 authored by Benjamin Cumming's avatar Benjamin Cumming Committed by GitHub
Browse files

Find morphology location from a coordinate (#1751)

Extend the `pwlin_place` interface to find the location of a morphology that is closest to a 3d coordinate.
* extend `arb::pwlin_place` interface
* python wrapper
* unit tests
* documentation

Addresses a feature request by @Helveg that I don't think ever had a ticket assigned to it.

Fixes #1661 and #1108.
parent e0c13e16
No related branches found
No related tags found
No related merge requests found
......@@ -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_;
};
......
#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
......@@ -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
^^^^^^^^^^
......
......@@ -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
......
......@@ -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)
......
......@@ -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);
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment