From 397a66b55c5ae7618e03eec7f9522c7510eb0d51 Mon Sep 17 00:00:00 2001
From: Sebastian Schmitt <sebastian.schmitt@kip.uni-heidelberg.de>
Date: Fri, 3 Dec 2021 13:16:44 +0100
Subject: [PATCH] lmorpho (#1746)
Add the lmorpho utility back to Arbor
The original utility was removed because it was a maintenance burden with no users.
Now users are requesting the feature, this PR adds the old code, and updates it
to work with the new morphology description API.
---
CMakeLists.txt | 1 +
lmorpho/CMakeLists.txt | 8 +
lmorpho/README.md | 115 ++++++++++++
lmorpho/lmorpho.cpp | 121 ++++++++++++
lmorpho/lsys_models.cpp | 252 +++++++++++++++++++++++++
lmorpho/lsys_models.hpp | 11 ++
lmorpho/lsystem.cpp | 403 ++++++++++++++++++++++++++++++++++++++++
lmorpho/lsystem.hpp | 134 +++++++++++++
8 files changed, 1045 insertions(+)
create mode 100644 lmorpho/CMakeLists.txt
create mode 100644 lmorpho/README.md
create mode 100644 lmorpho/lmorpho.cpp
create mode 100644 lmorpho/lsys_models.cpp
create mode 100644 lmorpho/lsys_models.hpp
create mode 100644 lmorpho/lsystem.cpp
create mode 100644 lmorpho/lsystem.hpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index ccdfe0db..ace085f3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -556,3 +556,4 @@ install(
cmake/FindUnwind.cmake
DESTINATION "${cmake_config_dir}")
+add_subdirectory(lmorpho)
diff --git a/lmorpho/CMakeLists.txt b/lmorpho/CMakeLists.txt
new file mode 100644
index 00000000..21262d16
--- /dev/null
+++ b/lmorpho/CMakeLists.txt
@@ -0,0 +1,8 @@
+add_executable(lmorpho lmorpho.cpp lsystem.cpp lsys_models.cpp)
+
+target_link_libraries(lmorpho PRIVATE arbor arbor-sup ext-tinyopt arborio)
+
+# TODO: resolve public headers
+target_link_libraries(lmorpho PRIVATE arbor-private-headers)
+
+install(TARGETS lmorpho RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
diff --git a/lmorpho/README.md b/lmorpho/README.md
new file mode 100644
index 00000000..561bbebd
--- /dev/null
+++ b/lmorpho/README.md
@@ -0,0 +1,115 @@
+# `lmporho` — generate random neuron morphologies
+
+`lmorpho` implements an L-system geometry generator for the purpose of
+generating a series of artificial neuron morphologies as described by a set of
+stasticial parameters.
+
+## Method
+
+The basic algorithm used is that of [Burke (1992)][burke1992], with extensions
+for embedding in 3-d taken from [Ascoli (2001)][ascoli2001].
+
+A dendritic tree is represented by a piece-wise linear embedding of a rooted
+tree into R³×Râº, describing the position in space and non-negative radius.
+Morphologies are represented as a sequence of these trees together with a point
+and radius describing a (spherical) soma.
+
+The generation process for a tree successively grows the terminal points by
+some length ΔL each step, changing its orientation each step according to a
+random distribution. After growing, there is a radius-dependent probability
+that it will bifurcate or terminate. Parameters describing these various
+probability distributions constitute the L-system model.
+
+The distributions used in [Ascoli (2001)][ascoli2001] can be constant (i.e. not
+random), uniform over some interval, truncated normal, or a mixture
+distribution of these. In the present version of the code, mixture
+distributions are not supported.
+
+## Output
+
+Generated morphologies can be output either in
+[SWC format](http://research.mssm.edu/cnic/swc.html), or as 'parent vectors',
+which describe the topology of the associated tree.
+
+Entry _i_ of the (zero-indexed) parent vector gives the index of the proximal
+unbranched section to which it connects. Somata have a parent index of -1.
+
+Morphology information can be written one each to separate files or
+concatenated. If concatenated, the parent vector indices will be shifted as
+required to maintain consistency.
+
+## Models
+
+Two L-system parameter sets are provided as 'built-in'; if neither model is
+selected, a simple, non-biological default model is used.
+
+As mixture distributions are not yet supported in `lmorpho`, these models are
+simplified approximations to those in the literature as described below.
+
+### Motor neuron model
+
+The 'motoneuron' model corresponds to the adult cat spinal alpha-motoneuron
+model described in [Ascoli (2001)][ascoli2001], based in turn on the model of
+[Burke (1992)][burke1992].
+
+The model parameters in Ascoli are not complete; missing parameters were taken
+from the corresponding ArborVitae model parameters in the same paper, and
+somatic diameters were taken from [Cullheim (1987)][cullheim1987].
+
+### Purkinke neuron model
+
+The 'purkinke' model corresponds to the guinea pig Purkinje cell model from
+[Ascoli (2001)][ascoli2001], which was fit to data from [Rapp
+(1994)][rapp1994]. Somatic diameters are a simple fit to the three measurements
+in the original source.
+
+Some required parameters are not recorded in [Ascoli (2001)][ascoli2001];
+the correlation component `diam_child_a` and the `length_step` ΔL are
+taken from the corresponding L-Neuron parameter description in the
+[L-Neuron database](l-neuron-database). These should be verified by
+fitting against the digitized morphologies as found in the database.
+
+Produced neurons from this model do not look especially realistic.
+
+### Other models
+
+There is not yet support for loading in arbitrary parameter sets, or for
+translating or approximating the parameters that can be found in the
+[L-Neuron database][l-neuron-database]. Support for parameter estimation
+from a population of discretized morphologies would be useful.
+
+## References
+
+1. Ascoli, G. A., Krichmar, J. L., Scorcioni, R., Nasuto, S. J., Senft, S. L.
+ and Krichmar, G. L. (2001). Computer generation and quantitative morphometric
+ analysis of virtual neurons. _Anatomy and Embryology_ _204_(4), 283–301.
+ [DOI: 10.1007/s004290100201][ascoli2001]
+
+2. Burke, R. E., Marks, W. B. and Ulfhake, B. (1992). A parsimonious description
+ of motoneuron dendritic morphology using computer simulation.
+ _The Journal of Neuroscience_ _12_(6), 2403–2416. [online][burke1992]
+
+3. Cullheim, S., Fleshman, J. W., Glenn, L. L., and Burke, R. E. (1987).
+ Membrane area and dendritic structure in type-identified triceps surae alpha
+ motoneurons. _The Journal of Comparitive Neurology_ _255_(1), 68–81.
+ [DOI: 10.1002/cne.902550106][cullheim1987]
+
+4. Rapp, M., Segev, I. and Yaom, Y. (1994). Physiology, morphology and detailed
+ passive models of guinea-pig cerebellar Purkinje cells.
+ _The Journal of Physiology_, _474_, 101–118.
+ [DOI: 10.1113/jphysiol.1994.sp020006][rapp1994].
+
+
+[ascoli2001]: http://dx.doi.org/10.1007/s004290100201
+ "Ascoli et al. (2001). Computer generation […] of virtual neurons."
+
+[burke1992]: http://www.jneurosci.org/content/12/6/2403
+ "Burke et al. (1992). A parsimonious description of motoneuron dendritic morphology […]"
+
+[cullheim1987]: http://dx.doi.org/10.1002/cne.902550106
+ "Cullheim et al. (1987). Membrane area and dendritic structure in […] alpha motoneurons."
+
+[rapp1994]: http://dx.doi.org/10.1113/jphysiol.1994.sp020006
+ "Rapp et al. (1994). Physiology, morphology […] of guinea-pig cerebellar Purkinje cells."
+
+[l-neuron-database]: http://krasnow1.gmu.edu/cn3/L-Neuron/database/index.html
diff --git a/lmorpho/lmorpho.cpp b/lmorpho/lmorpho.cpp
new file mode 100644
index 00000000..d6ed240f
--- /dev/null
+++ b/lmorpho/lmorpho.cpp
@@ -0,0 +1,121 @@
+#include <algorithm>
+#include <iostream>
+#include <fstream>
+#include <map>
+#include <sstream>
+#include <vector>
+
+#include <tinyopt/tinyopt.h>
+
+#include "lsystem.hpp"
+#include "lsys_models.hpp"
+
+#include <arbor/cable_cell.hpp>
+#include <arborio/cableio.hpp>
+#include <arborio/label_parse.hpp>
+
+const char* usage_str =
+"[OPTION]...\n"
+"\n"
+" -n, --count=N Number of morphologies to generate.\n"
+" -m, --model=MODEL Use L-system MODEL for generation (see below).\n"
+" --acc=FILE Output morphologies as ACC to FILE.\n"
+" -h, --help Emit this message and exit.\n"
+"\n"
+"Generate artificial neuron morphologies based on L-system descriptions.\n"
+"\n"
+"If a FILE argument contains a '%', then one file will be written for\n"
+"each generated morphology, with the '%' replaced by the index of the\n"
+"morphology, starting from zero.\n"
+"A FILE argument of '-' corresponds to standard output.\n"
+"\n"
+"Currently supported MODELs:\n"
+" motoneuron Adult cat spinal alpha-motoneurons, based on models\n"
+" and data in Burke 1992 and Ascoli 2001.\n"
+" purkinje Guinea pig Purkinje cells, basd on models and data\n"
+" from Rapp 1994 and Ascoli 2001.\n";
+
+int main(int argc, char** argv) {
+ // options
+ int n_morph = 1;
+ std::optional<unsigned> rng_seed = 0;
+ lsys_param P;
+ bool help = false;
+ std::optional<std::string> acc_file;
+
+ std::pair<const char*, const lsys_param*> models[] = {
+ {"motoneuron", &alpha_motoneuron_lsys},
+ {"purkinje", &purkinje_lsys}
+ };
+
+
+ try {
+ for (auto arg = argv+1; *arg; ) {
+ bool ok = false;
+ ok |= (n_morph << to::parse<int>(arg, "--n", "--count")).has_value();
+ ok |= (rng_seed << to::parse<int>(arg, "--s", "--seed")).has_value();
+ ok |= (acc_file << to::parse<std::string>(arg, "-f", "--acc")).has_value();
+ const lsys_param* o;
+ if((o << to::parse<const lsys_param*>(arg, to::keywords(models), "-m", "--model")).has_value()) {
+ P = *o;
+ }
+ ok |= (help << to::parse(arg, "-h", "--help")).has_value();
+ if (!ok) throw to::option_error("unrecognized argument", *arg);
+ }
+
+ if (help) {
+ to::usage(argv[0], usage_str);
+ return 0;
+ }
+
+ std::minstd_rand g;
+ if (rng_seed) {
+ g.seed(rng_seed.value());
+ }
+
+ using namespace arborio::literals;
+
+ auto apical = apical_lsys;
+ auto basal = basal_lsys;
+
+ for (int i = 0; i < n_morph; ++i) {
+ const arb::segment_tree morph = generate_morphology(P.diam_soma, std::vector<lsys_param>{apical, basal}, g);
+
+ arb::label_dict labels;
+ labels.set("soma", "(tag 1)"_reg);
+ labels.set("basal", "(tag 3)"_reg);
+ labels.set("apical", "(tag 4)"_reg);
+ arb::decor decor;
+
+ arb::cable_cell cell(morph, labels, decor);
+
+ if(acc_file) {
+ std::string filename = acc_file.value();
+ const std::string toReplace("%");
+ auto pos = filename.find(toReplace);
+ if(pos != std::string::npos) {
+ filename.replace(pos, toReplace.length(), std::to_string(i));
+ }
+
+ if(filename == "-") {
+ arborio::write_component(std::cout, cell);
+ } else {
+ std::ofstream of(filename);
+ arborio::write_component(of, cell);
+ }
+
+ }
+ }
+
+ } catch (to::option_error& e) {
+ std::cerr << argv[0] << ": " << e.what() << "\n";
+ std::cerr << "Try '" << argv[0] << " --help' for more information.\n";
+ std::exit(2);
+ }
+ catch (std::exception& e) {
+ std::cerr << "caught exception: " << e.what() << "\n";
+ std::exit(1);
+ }
+
+}
+
diff --git a/lmorpho/lsys_models.cpp b/lmorpho/lsys_models.cpp
new file mode 100644
index 00000000..6b560b90
--- /dev/null
+++ b/lmorpho/lsys_models.cpp
@@ -0,0 +1,252 @@
+#include <cmath>
+
+#include "lsystem.hpp"
+#include "lsys_models.hpp"
+
+static constexpr double inf = INFINITY;
+
+// Predefined parameters for two classes of neurons. Numbers taken primarily
+// from Ascoli et al. 2001, but some details (soma diameters for example)
+// taken from earlier literature.
+//
+// Without mixture distribution support, these aren't entirely faithful to
+// the original sources.
+//
+// Refer to `README.md` for references and more information.
+// https://krasnow1.gmu.edu/cn3/L-Neuron/database/moto/bu/motoburk.prm
+lsys_param make_alpha_motoneuron_lsys() {
+ lsys_param L;
+
+ // Soma diameter [µm].
+ L.diam_soma = { 47.0, 65.5, 57.6, 5.0 }; // somadiam
+
+ // Number of dendritic trees (rounded to nearest int).
+ L.n_tree = { 8.0, 16.0 }; // numtree
+
+ // Initial dendrite diameter [µm]. (Dstem)
+ L.diam_initial = { 3.0, 12.0 }; // initdiam
+
+ // Dendrite step length [µm]. (ΔL)
+ L.length_step = { 25 };
+
+ // Initial roll (intrinsic rotation about x-axis) [degrees].
+ L.roll_initial = { -180.0, 180.0 };
+
+ // Initial pitch (intrinsic rotation about y-axis) [degrees].
+ L.pitch_initial = { 0.0, 180.0 }; // biforient
+
+ // Tortuousness: roll within section over ΔL [degrees].
+ L.roll_section = { -180.0, 180.0 };
+
+ // Tortuousness: pitch within section over ΔL [degrees].
+ L.pitch_section = { -15.0, 15.0, 0.0, 5.0 };
+
+ // Taper rate: diameter decrease per unit length.
+ L.taper = { -1.25e-3 };
+
+ // Branching torsion: roll at branch point [degrees].
+ L.roll_at_branch = { 0.0, 180.0 };
+
+ // Branch angle between siblings [degrees].
+ L.branch_angle = { 1.0, 172.0, 45.0, 20.0 }; // bifamplitude
+
+ // Correlated child branch radius distribution parameters.
+ L.diam_child_a = -0.2087;
+ L.diam_child_r = { 0.2, 1.8, 0.8255, 0.2125 };
+
+ // Termination probability parameters [1/µm].
+ L.pterm_k1 = 2.62e-2;
+ L.pterm_k2 = -2.955;
+
+ // Bifurcation probability parameters [1/μm].
+ L.pbranch_ov_k1 = 2.58e-5;
+ L.pbranch_ov_k2 = 2.219;
+ L.pbranch_nov_k1 = 2.34e-3;
+ L.pbranch_nov_k2 = 0.194;
+
+ // Absolute maximum dendritic extent [µm].
+ L.max_extent = 2000;
+
+ return L;
+}
+
+lsys_param alpha_motoneuron_lsys = make_alpha_motoneuron_lsys();
+
+// Some parameters missing from the literature are taken from
+// the `purkburk.prm` L-Neuron parameter file:
+// http://krasnow1.gmu.edu/cn3/L-Neuron/database/purk/bu/purkburk.prm
+// Note 'Pov' and 'Pterm' numbers are swapped in Ascoli fig 3B.
+
+lsys_param make_purkinje_lsys() {
+ lsys_param L;
+
+ // Soma diameter [µm].
+ // (Gaussian fit to 3 samples rom Rapp 1994); truncate to 2σ.
+ L.diam_soma = { 22.12, 27.68, 24.9, 1.39 };
+
+ // Number of dendritic trees (rounded to nearest int).
+ L.n_tree = { 1 };
+
+ // Initial dendrite diameter [µm]. (Dstem)
+ L.diam_initial = { 4.8, 7.6, 6.167, 1.069 };
+
+ // Dendrite step length [µm]. (ΔL)
+ L.length_step = { 2.3 }; // from `purkburk.prm`
+
+ // Tortuousness: roll within section over ΔL [degrees].
+ L.roll_section = { 0 };
+
+ // Tortuousness: pitch within section over ΔL [degrees].
+ L.pitch_section = { -5, 5 }; // from `purkburk.prm`
+
+ // Taper rate: diameter decrease per unit length.
+ L.taper = { -inf, inf, -0.010, 0.022 };
+
+ // Branching torsion: roll at branch point [degrees].
+ L.roll_at_branch = { -inf, inf, 0.0, 6.5 }; // from `purkburj.prm`
+
+ // Branch angle between siblings [degrees].
+ L.branch_angle = { 1.0, 179.0, 71.0, 40.0 };
+
+ // Correlated child branch radius distribution parameters.
+ L.diam_child_a = 5.47078; // from `purkburk.prm`
+ L.diam_child_r = { 0, inf, 0.112, 0.035 };
+
+ // Termination probability parameters [1/µm].
+ L.pterm_k1 = 0.1973; // from `purkburj.prm`; 0.4539 in Ascoli.
+ L.pterm_k2 = -1.1533;
+
+ // Bifurcation probability parameters [1/μm].
+ L.pbranch_ov_k1 = 0.1021; // from `purkburk.prm`; 0.0355 in Ascoli.
+ L.pbranch_ov_k2 = 0.7237;
+ L.pbranch_nov_k1 = 0.1021; // from `purkburk.prm`; 0.2349 in Ascoli.
+ L.pbranch_nov_k2 = -0.0116;
+
+ // Absolute maximum dendritic extent [µm].
+ L.max_extent = 20000;
+
+ return L;
+}
+
+lsys_param purkinje_lsys = make_purkinje_lsys();
+
+
+lsys_param make_apical_lsys() {
+ lsys_param L;
+
+ // Soma diameter [µm].
+ L.diam_soma = { 47.0, 65.5, 57.6, 5.0 }; // somadiam
+
+ // Number of dendritic trees (rounded to nearest int).
+ L.n_tree = { 1 }; // numtree
+
+ // Initial dendrite diameter [µm]. (Dstem)
+ L.diam_initial = { 10, 20 }; // initdiam
+
+ // Dendrite step length [µm]. (ΔL)
+ L.length_step = { 25 };
+
+ // Initial roll (intrinsic rotation about x-axis) [degrees].
+ L.roll_initial = { 0 };
+
+ // Initial pitch (intrinsic rotation about y-axis) [degrees].
+ L.pitch_initial = { 0 }; // biforient
+
+ // Tortuousness: roll within section over ΔL [degrees].
+ L.roll_section = { -180.0, 180.0 };
+
+ // Tortuousness: pitch within section over ΔL [degrees].
+ L.pitch_section = { -15.0, 15.0, 0.0, 5.0 };
+
+ // Taper rate: diameter decrease per unit length.
+ L.taper = { -1.25e-3 };
+
+ // Branching torsion: roll at branch point [degrees].
+ L.roll_at_branch = { 0.0, 180.0 };
+
+ // Branch angle between siblings [degrees].
+ L.branch_angle = { 1.0, 172.0, 45.0, 20.0 }; // bifamplitude
+
+ // Correlated child branch radius distribution parameters.
+ L.diam_child_a = -0.2087;
+ L.diam_child_r = { 0.2, 1.8, 0.8255, 0.2125 };
+
+ // Termination probability parameters [1/µm].
+ L.pterm_k1 = 2.62e-2;
+ L.pterm_k2 = -2.955;
+
+ // Bifurcation probability parameters [1/μm].
+ L.pbranch_ov_k1 = 0.58e-5;
+ L.pbranch_ov_k2 = 2.219;
+ L.pbranch_nov_k1 = 2.34e-3;
+ L.pbranch_nov_k2 = 0.194;
+
+ // Absolute maximum dendritic extent [µm].
+ L.max_extent = 2000;
+
+ L.tag = 4;
+
+ return L;
+}
+
+lsys_param apical_lsys = make_apical_lsys();
+
+lsys_param make_basal_lsys() {
+ lsys_param L;
+
+ // Soma diameter [µm].
+ L.diam_soma = { 47.0, 65.5, 57.6, 5.0 }; // somadiam
+
+ // Number of dendritic trees (rounded to nearest int).
+ L.n_tree = { 1, 20, 6, 2}; // numtree
+
+ // Initial dendrite diameter [µm]. (Dstem)
+ L.diam_initial = { 3.0, 12.0 }; // initdiam
+
+ // Dendrite step length [µm]. (ΔL)
+ L.length_step = { 25 };
+
+ // Initial roll (intrinsic rotation about x-axis) [degrees].
+ L.roll_initial = { -180, 180 };
+
+ // Initial pitch (intrinsic rotation about y-axis) [degrees].
+ L.pitch_initial = { 140, 220 }; // biforient
+
+ // Tortuousness: roll within section over ΔL [degrees].
+ L.roll_section = { -180.0, 180.0 };
+
+ // Tortuousness: pitch within section over ΔL [degrees].
+ L.pitch_section = { -15.0, 15.0, 0.0, 5.0 };
+
+ // Taper rate: diameter decrease per unit length.
+ L.taper = { -1.25e-3 };
+
+ // Branching torsion: roll at branch point [degrees].
+ L.roll_at_branch = { 0.0, 180.0 };
+
+ // Branch angle between siblings [degrees].
+ L.branch_angle = { 1.0, 172.0, 45.0, 20.0 }; // bifamplitude
+
+ // Correlated child branch radius distribution parameters.
+ L.diam_child_a = -0.2087;
+ L.diam_child_r = { 0.2, 1.8, 0.8255, 0.2125 };
+
+ // Termination probability parameters [1/µm].
+ L.pterm_k1 = 2.62e-2;
+ L.pterm_k2 = -2.955;
+
+ // Bifurcation probability parameters [1/μm].
+ L.pbranch_ov_k1 = 1.58e-5;
+ L.pbranch_ov_k2 = 2.219;
+ L.pbranch_nov_k1 = 2.34e-3;
+ L.pbranch_nov_k2 = 0.194;
+
+ // Absolute maximum dendritic extent [µm].
+ L.max_extent = 1000;
+
+ L.tag = 3;
+
+ return L;
+}
+
+lsys_param basal_lsys = make_basal_lsys();
diff --git a/lmorpho/lsys_models.hpp b/lmorpho/lsys_models.hpp
new file mode 100644
index 00000000..35a02b8f
--- /dev/null
+++ b/lmorpho/lsys_models.hpp
@@ -0,0 +1,11 @@
+#pragma once
+
+#include "lsystem.hpp"
+
+// Predefined parameters for two classes of neurons.
+// See lsys_models.cc for details.
+
+extern lsys_param alpha_motoneuron_lsys;
+extern lsys_param apical_lsys;
+extern lsys_param basal_lsys;
+extern lsys_param purkinje_lsys;
diff --git a/lmorpho/lsystem.cpp b/lmorpho/lsystem.cpp
new file mode 100644
index 00000000..fc6ae025
--- /dev/null
+++ b/lmorpho/lsystem.cpp
@@ -0,0 +1,403 @@
+#include <cmath>
+#include <iostream>
+#include <random>
+#include <stack>
+#include <vector>
+#include <map>
+
+#include <arbor/math.hpp>
+
+#include "lsystem.hpp"
+
+using namespace arb::math;
+
+// L-system implementation.
+
+// Random distribution used in the specification of L-system parameters.
+// It can be a constant, uniform over an interval, or truncated normal.
+// Note that mixture distributions of the above should be, but isn't yet,
+// implemented.
+
+struct lsys_distribution {
+ using result_type = double;
+ using param_type = lsys_distribution_param;
+
+ lsys_distribution(const param_type& p): p_(p) { init(); }
+
+ lsys_distribution(): p_() { init(); }
+ lsys_distribution(result_type x): p_(x) { init(); }
+ lsys_distribution(result_type lb, result_type ub): p_(lb, ub) { init(); }
+ lsys_distribution(result_type lb, result_type ub, result_type mu, result_type sigma):
+ p_(lb, ub, mu, sigma)
+ {
+ init();
+ }
+
+ param_type param() const { return p_; }
+ void param(const param_type& pbis) { p_=pbis; }
+
+ result_type min() const { return p_.lb; }
+ result_type max() const {
+ return p_.kind==param_type::constant? min(): p_.ub;
+ }
+
+ bool operator==(const lsys_distribution& y) const {
+ return p_==y.p_;
+ }
+
+ bool operator!=(const lsys_distribution& y) const {
+ return p_!=y.p_;
+ }
+
+ // omitting out ostream, istream methods
+ // (...)
+
+ template <typename Gen>
+ result_type operator()(Gen& g, const param_type &p) const {
+ lsys_distribution dist(p);
+ return dist(g);
+ }
+
+ template <typename Gen>
+ result_type operator()(Gen& g) const {
+ switch (p_.kind) {
+ case param_type::constant:
+ return p_.lb;
+ case param_type::uniform:
+ return uniform_dist(g);
+ case param_type::normal:
+ // TODO: replace with a robust truncated normal generator!
+ for (;;) {
+ result_type r = normal_dist(g);
+ if (r>=p_.lb && r<=p_.ub) return r;
+ }
+ }
+ return 0;
+ }
+
+private:
+ param_type p_;
+ mutable std::uniform_real_distribution<double> uniform_dist;
+ mutable std::normal_distribution<double> normal_dist;
+
+ void init() {
+ switch (p_.kind) {
+ case param_type::uniform:
+ uniform_dist = std::uniform_real_distribution<double>(p_.lb, p_.ub);
+ break;
+ case param_type::normal:
+ normal_dist = std::normal_distribution<result_type>(p_.mu, p_.sigma);
+ break;
+ default: ;
+ }
+ }
+};
+
+class burke_correlated_r {
+ double a_;
+ lsys_distribution r_;
+
+public:
+ struct result_type { double rho1, rho2; };
+
+ burke_correlated_r(double a, lsys_distribution_param r): a_(a), r_(r) {}
+
+ template <typename Gen>
+ result_type operator()(Gen& g) const {
+ double r1 = r_(g);
+ double r2 = r_(g);
+ return { r1+a_*r2, r2+a_*r1 };
+ }
+};
+
+class burke_exp_test {
+ double k1_;
+ double k2_;
+
+public:
+ using result_type = int;
+
+ burke_exp_test(double k1, double k2): k1_(k1), k2_(k2) {}
+
+ template <typename Gen>
+ result_type operator()(double delta_l, double radius, Gen& g) const {
+ std::bernoulli_distribution B(delta_l*k1_*std::exp(k2_*radius*2.0));
+ return B(g);
+ }
+};
+
+class burke_exp2_test {
+ double ak1_;
+ double ak2_;
+ double bk1_;
+ double bk2_;
+
+public:
+ using result_type = int;
+
+ burke_exp2_test(double ak1, double ak2, double bk1, double bk2):
+ ak1_(ak1), ak2_(ak2), bk1_(bk1), bk2_(bk2) {}
+
+ template <typename Gen>
+ result_type operator()(double delta_l, double radius, Gen& g) const {
+ double diam = 2.*radius;
+ std::bernoulli_distribution B(
+ delta_l*std::min(ak1_*std::exp(ak2_*diam), bk1_*std::exp(bk2_*diam)));
+
+ return B(g);
+ }
+};
+
+// L-system parameter set with instantiated distributions.
+struct lsys_sampler {
+ // Create distributions once from supplied parameter set.
+ explicit lsys_sampler(const lsys_param& P):
+ diam_soma(P.diam_soma),
+ n_tree(P.n_tree),
+ diam_initial(P.diam_initial),
+ length_step(P.length_step),
+ roll_initial(P.roll_initial),
+ pitch_initial(P.pitch_initial),
+ roll_section(P.roll_section),
+ pitch_section(P.pitch_section),
+ taper(P.taper),
+ roll_at_branch(P.roll_at_branch),
+ branch_angle(P.branch_angle),
+ diam_child(P.diam_child_a, P.diam_child_r),
+ term_test(P.pterm_k1, P.pterm_k2),
+ branch_test(P.pbranch_ov_k1, P.pbranch_ov_k2,
+ P.pbranch_nov_k1, P.pbranch_nov_k2),
+ max_extent(P.max_extent)
+ {}
+
+ lsys_distribution diam_soma;
+ lsys_distribution n_tree;
+ lsys_distribution diam_initial;
+ lsys_distribution length_step;
+ lsys_distribution roll_initial;
+ lsys_distribution pitch_initial;
+ lsys_distribution roll_section;
+ lsys_distribution pitch_section;
+ lsys_distribution taper;
+ lsys_distribution roll_at_branch;
+ lsys_distribution branch_angle;
+ burke_correlated_r diam_child;
+ burke_exp_test term_test;
+ burke_exp2_test branch_test;
+ double max_extent;
+};
+
+struct section_point {
+ double x;
+ double y;
+ double z;
+ double r;
+ friend std::ostream& operator<<(std::ostream& os, const section_point& obj);
+};
+
+std::ostream& operator<<(std::ostream& os, const section_point& obj) {
+ os << obj.x << " " << obj.y << " " << obj.z << " " << obj.r << '\n';
+ return os;
+}
+
+struct section_tip {
+ section_point point = {0., 0., 0., 0.};
+ quaternion rotation = {1., 0., 0., 0.};
+ double somatic_distance = 0.;
+};
+
+
+struct section_geometry {
+ unsigned id = 0; // ids should be contigously numbered from 1 in the morphology.
+ unsigned parent_id = 0;
+ bool terminal = false;
+ std::vector<section_point> points;
+ double length = 0; // µm
+
+ section_geometry() = default;
+ section_geometry(unsigned id, unsigned parent_id, bool terminal, std::vector<section_point> points, double length) :
+ id(id), parent_id(parent_id), terminal(terminal), points(std::move(points)), length(length)
+ {}
+
+ friend std::ostream& operator<<(std::ostream& os, const section_geometry& obj);
+};
+
+std::ostream& operator<<(std::ostream& os, const section_geometry& obj) {
+ os << "id: " << obj.id << '\n';
+ os << "parent_id: " << obj.parent_id << '\n';
+ os << "terminal: " << obj.terminal << '\n';
+ os << "length: " << obj.length << '\n';
+ os << "points: " << '\n';
+ for (auto& p : obj.points) {
+ os << p << '\n';
+ }
+ return os;
+}
+
+struct grow_result {
+ std::vector<section_point> points;
+ std::vector<section_tip> children;
+ double length = 0.;
+};
+
+constexpr double deg_to_rad = 3.1415926535897932384626433832795/180.0;
+
+template <typename Gen>
+grow_result grow(section_tip tip, const lsys_sampler& S, Gen &g) {
+ constexpr quaternion xaxis = {0, 1, 0, 0};
+ static std::uniform_real_distribution<double> U;
+
+ grow_result result;
+ std::vector<section_point>& points = result.points;
+
+ points.push_back(tip.point);
+ for (;;) {
+ quaternion step = xaxis^tip.rotation;
+ double dl = S.length_step(g);
+ tip.point.x += dl*step.x;
+ tip.point.y += dl*step.y;
+ tip.point.z += dl*step.z;
+ tip.point.r += dl*0.5*S.taper(g);
+ tip.somatic_distance += dl;
+ result.length += dl;
+
+ if (tip.point.r<0) tip.point.r = 0;
+
+ double phi = S.roll_section(g);
+ double theta = S.pitch_section(g);
+ tip.rotation *= rotation_x(deg_to_rad*phi);
+ tip.rotation *= rotation_y(deg_to_rad*theta);
+
+ points.push_back(tip.point);
+
+ if (tip.point.r==0 || tip.somatic_distance>=S.max_extent) {
+ return result;
+ }
+
+ if (S.branch_test(dl, tip.point.r, g)) {
+ auto r = S.diam_child(g);
+ // ignore branch if we get non-positive radii
+ if (r.rho1>0 && r.rho2>0) {
+ double branch_phi = S.roll_at_branch(g);
+ double branch_angle = S.branch_angle(g);
+ double branch_theta1 = U(g)*branch_angle;
+ double branch_theta2 = branch_theta1-branch_angle;
+
+ tip.rotation *= rotation_x(deg_to_rad*branch_phi);
+
+ section_tip t1 = tip;
+ t1.point.r = t1.point.r * r.rho1;
+ t1.rotation *= rotation_y(deg_to_rad*branch_theta1);
+
+ section_tip t2 = tip;
+ t2.point.r = t2.point.r * r.rho2;
+ t2.rotation *= rotation_y(deg_to_rad*branch_theta2);
+
+ result.children = {t1, t2};
+ return result;
+ }
+ }
+
+ if (S.term_test(dl, tip.point.r, g)) {
+ return result;
+ }
+ }
+}
+
+std::vector<section_geometry> generate_sections(double soma_radius, lsys_param P, lsys_generator &g) {
+ constexpr quaternion xaxis = {0, 1, 0, 0};
+
+ lsys_sampler S(P);
+
+ struct section_start {
+ section_tip tip;
+ unsigned parent_id;
+ };
+ std::stack<section_start> starts;
+ unsigned next_id = 1u;
+
+ int n = (int)std::round(S.n_tree(g));
+ for (int i=0; i<n; ++i) {
+ double phi = S.roll_initial(g);
+ double theta = S.pitch_initial(g);
+ double radius = 0.5*S.diam_initial(g);
+
+ section_tip tip;
+ tip.rotation = rotation_x(deg_to_rad*phi)*rotation_y(deg_to_rad*theta);
+ tip.somatic_distance = 0.0;
+
+ auto p = (soma_radius*xaxis)^tip.rotation;
+ tip.point = {p.x, p.y, p.z, radius};
+
+ starts.push({tip, 0u});
+ }
+
+ std::vector<section_geometry> sections;
+
+ while (!starts.empty() && sections.size() < P.max_sections) {
+ auto start = starts.top();
+ starts.pop();
+
+ auto branch = grow(start.tip, S, g);
+ section_geometry section{next_id++, start.parent_id, branch.children.empty(), std::move(branch.points), branch.length};
+
+ for (auto child: branch.children) {
+ starts.push({child, section.id});
+ }
+
+ sections.push_back(std::move(section));
+ }
+
+ return sections;
+}
+
+arb::segment_tree generate_morphology(const lsys_distribution_param& soma, std::vector<lsys_param> Ps, lsys_generator &g) {
+
+ double const soma_diameter = lsys_distribution(soma)(g);
+ double const soma_radius = 0.5*soma_diameter;
+
+ arb::segment_tree morph;
+ morph.append(
+ arb::mnpos, arb::mpoint{-soma_diameter, 0, 0, soma_diameter},
+ arb::mpoint{soma_diameter, 0, 0, soma_diameter}, 1);
+
+ for(auto P : Ps) {
+ auto sections = generate_sections(soma_radius, P, g);
+
+ // map from internal representation to Arbor ids
+ std::map<size_t, arb::msize_t> parent_end_id;
+
+ // soma
+ parent_end_id[0] = 0;
+
+ // dendrites:
+ for (auto& sec: sections) {
+
+ // the first section must have the soma as parent
+ size_t parent = parent_end_id.at(sec.parent_id);
+
+ const auto& points = sec.points;
+ if (points.size() < 2) {
+ throw std::runtime_error("segment has only 1 point");
+ }
+
+ // Include first point only for dendrites segments attached to soma.
+ if (sec.parent_id==0) {
+ parent = morph.append(parent,
+ arb::mpoint{points[0].x, points[0].y, points[0].z, points[0].r},
+ arb::mpoint{points[1].x, points[1].y, points[1].z, points[1].r}, P.tag);
+ }
+
+ // Remaining points.
+ for (unsigned i = 1; i < points.size(); ++i) {
+ const auto& p = points[i];
+ parent = morph.append(parent, arb::mpoint{p.x, p.y, p.z, p.r}, P.tag);
+ }
+
+ parent_end_id[sec.id] = parent;
+ }
+ }
+
+ return morph;
+}
+
diff --git a/lmorpho/lsystem.hpp b/lmorpho/lsystem.hpp
new file mode 100644
index 00000000..3130ba64
--- /dev/null
+++ b/lmorpho/lsystem.hpp
@@ -0,0 +1,134 @@
+#pragma once
+
+#include <random>
+
+#include <arbor/morph/segment_tree.hpp>
+
+struct lsys_param;
+
+using lsys_generator = std::minstd_rand;
+
+class lsys_distribution_param;
+arb::segment_tree generate_morphology(const lsys_distribution_param& soma, std::vector<lsys_param> Ps, lsys_generator& g);
+
+// The distribution parameters used in the specification of the L-system parameters.
+// The distribution can be a constant, uniform over an interval, or truncated normal.
+// Note that mixture distributions of the above should be, but isn't yet,
+// implemented.
+
+struct lsys_distribution_param {
+ using result_type = double;
+ enum lsys_kind {
+ constant, uniform, normal
+ };
+
+ lsys_kind kind;
+ result_type lb; // lower bound
+ result_type ub; // upper bound (used for only non-constant)
+ result_type mu; // mean (used only for truncated Gaussian)
+ result_type sigma; // s.d. (used only for truncated Gaussian)
+
+ bool operator==(const lsys_distribution_param &p) const {
+ if (kind!=p.kind) return false;
+ switch (kind) {
+ case constant:
+ return lb == p.lb;
+ case uniform:
+ return lb == p.lb && ub == p.ub;
+ case normal:
+ return lb == p.lb && ub == p.ub && mu == p.mu && sigma == p.sigma;
+ }
+ return false;
+ }
+
+ bool operator!=(const lsys_distribution_param &p) const {
+ return !(*this==p);
+ }
+
+ // zero or one parameter => constant
+ lsys_distribution_param():
+ kind(constant), lb(0.0) {}
+
+ lsys_distribution_param(result_type x):
+ kind(constant), lb(x) {}
+
+ // two paramter => uniform
+ lsys_distribution_param(result_type lb, result_type ub):
+ kind(uniform), lb(lb), ub(ub) {}
+
+ // four paramter => truncated normal
+ lsys_distribution_param(result_type lb, result_type ub, result_type mu, result_type sigma):
+ kind(normal), lb(lb), ub(ub), mu(mu), sigma(sigma) {}
+};
+
+
+// Parameters as referenced in Ascoli 2001 are given in parentheses below.
+// Defaults below are chosen to be safe, boring and not biological.
+struct lsys_param {
+ // Soma diameter [µm].
+ lsys_distribution_param diam_soma = { 20 };
+
+ // Number of dendritic trees (rounded to nearest int). (Ntree)
+ lsys_distribution_param n_tree = { 1 };
+
+ // Initial dendrite diameter [µm]. (Dstem)
+ lsys_distribution_param diam_initial = { 2 };
+
+ // Dendrite step length [µm]. (ΔL)
+ lsys_distribution_param length_step = { 25 };
+
+ // Dendrites grow along x-axis after application of (instrinsic) rotations;
+ // roll (about x-axis) is applied first, followed by pitch (about y-axis).
+
+ // Initial roll (intrinsic rotation about x-axis) [degrees]. (Taz)
+ lsys_distribution_param roll_initial = { -45, 45 };
+
+ // Initial pitch (intrinsic rotation about y-axis) [degrees]. (Tel)
+ lsys_distribution_param pitch_initial = { -45, 45 };
+
+ // Tortuousness: roll within section over ΔL [degrees]. (Eaz)
+ lsys_distribution_param roll_section = { -45, 45 };
+
+ // Tortuousness: pitch within section over ΔL [degrees]. (Eel)
+ lsys_distribution_param pitch_section = { -45, 45 };
+
+ // Taper rate: diameter decrease per unit length. (TPRB)
+ lsys_distribution_param taper = { 0 };
+
+ // Branching torsion: roll at branch point [degrees]. (Btor)
+ lsys_distribution_param roll_at_branch = { 0 };
+
+ // Branch angle between siblings [degrees]. (Bamp)
+ lsys_distribution_param branch_angle = { 45 };
+
+ // Child branch diameter ratios are typically anticorrelated.
+ // The Burke 1992 parameterization describes the two child ratios
+ // as d1 = r1 + a*r2 and d2 = r2 + a*r1 where a is a constant
+ // determined by the correlation, and r1 and r2 are drawn from
+ // a common distribution (`d_child_r` below).
+ double diam_child_a = 0;
+ lsys_distribution_param diam_child_r = { 0.5 };
+
+ // Bifurcation and termination probabilities per unit-length below are given by
+ // P = k1 * exp(k2*diameter), described by the paremeters k1 and k2 [1/µm].
+
+ // Termination probability parameters [1/µm]. (k1trm, k2trm)
+ double pterm_k1 = 0.05;
+ double pterm_k2 = -2;
+
+ // Bifurcation probability is given by minimum of two probabilities, `pbranch_ov' and
+ // `pbranch_nov' below.
+ double pbranch_ov_k1 = 0.01;
+ double pbranch_ov_k2 = 0;
+ double pbranch_nov_k1 = 0.01;
+ double pbranch_nov_k2 = 0;
+
+ // Absolute maximum dendritic extent [µm]. (Forces termination of algorithm)
+ double max_extent = 2000;
+
+ // Absolute maximum number of unbranched sections. (Forces termination of algorithm)
+ unsigned max_sections = 10000;
+
+ size_t tag = 0;
+};
+
--
GitLab