Bugfixes for `lmorpho` (#170)

* Ignore dendrite branches with negative radii arising from correlated child diameter distribution.
* Fix fencepost errors in morphology discretization.
* Rename `tip.p` to `tip.point`.

lmorpho/lsys_models.cpp

@@ -110,7 +110,7 @@ lsys_param make_purkinje_lsys() {
 
     // Correlated child branch radius distribution parameters.
     L.diam_child_a = 5.47078; // from `purkburk.prm`
-    L.diam_child_r = { -inf, inf, 0.112, 0.035 };
+    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.

lmorpho/lsystem.cpp

@@ -188,7 +188,7 @@ struct lsys_sampler {
 };
 
 struct section_tip {
-    section_point p = {0., 0., 0., 0.};
+    section_point point = {0., 0., 0., 0.};
     quaternion rotation = {1., 0., 0., 0.};
     double somatic_distance = 0.;
 };
@@ -209,52 +209,55 @@ grow_result grow(section_tip tip, const lsys_sampler& S, Gen &g) {
     grow_result result;
     std::vector<section_point>& points = result.points;
 
-    points.push_back(tip.p);
+    points.push_back(tip.point);
     for (;;) {
         quaternion step = xaxis^tip.rotation;
         double dl = S.length_step(g);
-        tip.p.x += dl*step.x;
-        tip.p.y += dl*step.y;
-        tip.p.z += dl*step.z;
-        tip.p.r += dl*0.5*S.taper(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.p.r<0) tip.p.r = 0;
+        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.p);
+        points.push_back(tip.point);
 
-        if (tip.p.r==0 || tip.somatic_distance>=S.max_extent) {
+        if (tip.point.r==0 || tip.somatic_distance>=S.max_extent) {
             return result;
         }
 
-        if (S.branch_test(dl, tip.p.r, g)) {
-            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;
+        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);
+                tip.rotation *= rotation_x(deg_to_rad*branch_phi);
 
-            section_tip t1 = tip;
-            t1.p.r = t1.p.r * r.rho1;
-            t1.rotation *= rotation_y(deg_to_rad*branch_theta1);
+                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.p.r = t2.p.r * r.rho2;
-            t2.rotation *= rotation_y(deg_to_rad*branch_theta2);
+                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;
+                result.children = {t1, t2};
+                return result;
+            }
         }
 
-        if (S.term_test(dl, tip.p.r, g)) {
+        if (S.term_test(dl, tip.point.r, g)) {
             return result;
         }
     }
@@ -286,7 +289,7 @@ morphology generate_morphology(const lsys_param& P, lsys_generator &g) {
         tip.somatic_distance = 0.0;
 
         auto p = (soma_radius*xaxis)^tip.rotation;
-        tip.p = {p.x, p.y, p.z, radius};
+        tip.point = {p.x, p.y, p.z, radius};
 
         starts.push({tip, 0u});
     }

lmorpho/morphology.cpp

@@ -39,16 +39,19 @@ void section_geometry::segment(double dx) {
     double x = length/nseg;
 
     // Scan segments for sample points.
-    for (unsigned i = 1; i<npoint;) {
+    unsigned i = 1;
+    for (;;) {
         if (right>x) {
             double u = (x-left)/(right-left);
             sampled.push_back(lerp(points[i-1], points[i], u));
             unsigned k = sampled.size();
-            sampled_length += distance(sampled[k-2], sampled[k]);
+            sampled_length += distance(sampled[k-2], sampled[k-1]);
             x = k*length/nseg;
         }
         else {
             ++i;
+            if (i>=npoint) break;
+
             left = right;
             right = left+distance(points[i-1], points[i]);
         }