Repository: BlueBrain/MorphIO (language: python)

Issue:
Non-recursive deletion of a child changes the children order
```python
import tempfile

from numpy import testing as npt

from morphio.mut import Morphology

with tempfile.NamedTemporaryFile(suffix=".asc") as tfile:
    filepath = tfile.name
    with open(filepath, "w") as f:
        f.write(
            """
            ("CellBody"
              (Color Red)
              (CellBody)
              ( 0.1  0.1 0.0 0.1)
              (-0.1  0.1 0.0 0.1)
              (-0.1 -0.1 0.0 0.1)
              ( 0.1 -0.1 0.0 0.1)
            )

            ( (Color Cyan)
              (Axon)
              (0.0  0.0 0.0 2.0)
              (0.0 -4.0 0.0 2.0)
              (
                (0.0 -4.0 0.0 4.0)
                (0.0 -4.0 0.0 4.0)
                (0.0 -4.0 0.0 4.0)
                (
                    (6.0 -4.0 0.0 4.0)
                    (7.0 -5.0 0.0 4.0)
                |
                    (6.0 -4.0 0.0 4.0)
                    (8.0 -4.0 0.0 4.0)
                )
              |
                ( 0.0 -4.0 0.0 4.0)
                (-5.0 -4.0 0.0 4.0)
              )
            )
            """
        )

    morph = Morphology(filepath)
 
    axon = morph.root_sections[0]

    # child 1 points should remain the same after deletion
    child2_points = axon.children[1].points

    print(axon.children[0].points)
    print(axon.children[1].points)

    morph.delete_section(axon.children[0], recursive=False)

    print()
    print(axon.children[0].points)
    print(axon.children[1].points)

    npt.assert_array_equal(child2_points, axon.children[1].points)

```
The second child which is not deleted, becomes the first after the deletion.
```
Before:
[[ 0. -4.  0.]
 [ 0. -4.  0.]
 [ 0. -4.  0.]]

[[ 0. -4.  0.]
 [-5. -4.  0.]]

After:
[[ 0. -4.  0.]
 [-5. -4.  0.]]

[[ 0. -4.  0.]
 [ 6. -4.  0.]
 [ 7. -5.  0.]]

```


Discussion / hints:

Fix patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 86141481d..153e95cf6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,9 @@
+v3.3.4
+======
+
+Fixes:
+* Fix wrong section order when section is deleted (#412)
+
 v3.3.3
 ======
 
diff --git a/include/morphio/mut/morphology.h b/include/morphio/mut/morphology.h
index b93492b3c..138cb440e 100644
--- a/include/morphio/mut/morphology.h
+++ b/include/morphio/mut/morphology.h
@@ -180,7 +180,7 @@ class Morphology
     }
 
     /// Write file to H5, SWC, ASC format depending on filename extension
-    void write(const std::string& filename);
+    void write(const std::string& filename) const;
 
     void addAnnotation(const morphio::Property::Annotation& annotation) {
         _cellProperties->_annotations.push_back(annotation);
diff --git a/src/mut/morphology.cpp b/src/mut/morphology.cpp
index c1a8f2e5e..bb1326c04 100644
--- a/src/mut/morphology.cpp
+++ b/src/mut/morphology.cpp
@@ -1,5 +1,6 @@
+#include <algorithm>  // std::find_if
 #include <cassert>
-#include <cctype>  // std::tolower
+#include <cctype>    // std::tolower
 #include <iterator>  // std::back_inserter
 #include <sstream>
 #include <string>
@@ -15,9 +16,24 @@
 #include "../shared_utils.hpp"
 
 namespace {
-void _appendProperties(morphio::Property::PointLevel& to,
-                       const morphio::Property::PointLevel& from,
-                       int offset = 0) {
+using SectionP = std::shared_ptr<morphio::mut::Section>;
+
+void insertSectionsBeforeSection(std::vector<SectionP>& sections_to_update,
+                                 const std::vector<SectionP>& sections,
+                                 const SectionP& target_section) {
+    auto is_equal = [&target_section](const SectionP& section) {
+        return (section->id() == target_section->id()) && section->hasSameShape(*target_section);
+    };
+
+    const auto it = std::find_if(sections_to_update.begin(), sections_to_update.end(), is_equal);
+
+    // Note that the section with section_id is not removed at this step
+    sections_to_update.insert(it, sections.begin(), sections.end());
+}
+
+void appendProperties(morphio::Property::PointLevel& to,
+                      const morphio::Property::PointLevel& from,
+                      int offset = 0) {
     morphio::_appendVector(to._points, from._points, offset);
     morphio::_appendVector(to._diameters, from._diameters, offset);
 
@@ -71,8 +87,7 @@ Morphology::Morphology(const morphio::Morphology& morphology, unsigned int optio
 /**
    Return false if there is no duplicate point
  **/
-bool _checkDuplicatePoint(const std::shared_ptr<Section>& parent,
-                          const std::shared_ptr<Section>& current) {
+bool _checkDuplicatePoint(const SectionP& parent, const SectionP& current) {
     // Weird edge case where parent is empty: skipping it
     if (parent->points().empty()) {
         return true;
@@ -151,7 +166,7 @@ std::shared_ptr<Section> Morphology::appendRootSection(const Property::PointLeve
 }
 
 uint32_t Morphology::_register(const std::shared_ptr<Section>& section_) {
-    if (_sections.count(section_->id())) {
+    if (_sections.count(section_->id()) > 0) {
         throw SectionBuilderError("Section already exists");
     }
     _counter = std::max(_counter, section_->id()) + 1;
@@ -175,13 +190,12 @@ void Morphology::eraseByValue(std::vector<std::shared_ptr<Section>>& vec,
     vec.erase(std::remove(vec.begin(), vec.end(), section), vec.end());
 }
 
-void Morphology::deleteSection(std::shared_ptr<Section> section_, bool recursive) {
+
+void Morphology::deleteSection(const std::shared_ptr<Section> section_, bool recursive) {
     if (!section_) {
         return;
     }
 
-    unsigned int id = section_->id();
-
     if (recursive) {
         // The deletion must start by the furthest leaves, otherwise you may cut
         // the topology and forget to remove some sections
@@ -192,22 +206,33 @@ void Morphology::deleteSection(std::shared_ptr<Section> section_, bool recursive
             deleteSection(*it, false);
         }
     } else {
+        const auto section_id = section_->id();
+
         // Careful not to use a reference here or you will face reference invalidation problem
         // with vector resize
-        for (auto child : section_->children()) {
-            if (section_->isRoot()) {
-                _rootSections.push_back(child);
-            } else {
-                // Re-link children to their "grand-parent"
-                _children[_parent[id]].push_back(child);
-                _parent[child->id()] = _parent[id];
+        const auto children = section_->children();
+
+        if (section_->isRoot()) {
+            // put root section's children in its place
+            insertSectionsBeforeSection(_rootSections, children, section_);
+            eraseByValue(_rootSections, section_);
+        } else {
+            // set grandparent as section children's parent
+            for (const auto& child : children) {
+                _parent[child->id()] = _parent[section_id];
             }
+
+            auto& children_of_parent = _children[_parent[section_id]];
+
+            // put grandchildren at the position of the deleted section
+            insertSectionsBeforeSection(children_of_parent, children, section_);
+            eraseByValue(children_of_parent, section_);
         }
-        eraseByValue(_rootSections, section_);
-        eraseByValue(_children[_parent[id]], section_);
-        _children.erase(id);
-        _parent.erase(id);
-        _sections.erase(id);
+
+        // remove section id from connectivity and section lists
+        _children.erase(section_id);
+        _parent.erase(section_id);
+        _sections.erase(section_id);
     }
 }
 
@@ -275,7 +300,7 @@ Property::Properties Morphology::buildReadOnly() const {
 
     properties._cellLevel = *_cellProperties;
     properties._cellLevel._somaType = _soma->type();
-    _appendProperties(properties._somaLevel, _soma->point_properties_);
+    appendProperties(properties._somaLevel, _soma->point_properties_);
 
     for (auto it = depth_begin(); it != depth_end(); ++it) {
         const std::shared_ptr<Section>& section = *it;
@@ -286,7 +311,7 @@ Property::Properties Morphology::buildReadOnly() const {
         properties._sectionLevel._sections.push_back({start, parentOnDisk});
         properties._sectionLevel._sectionTypes.push_back(section->type());
         newIds[sectionId] = sectionIdOnDisk++;
-        _appendProperties(properties._pointLevel, section->point_properties_);
+        appendProperties(properties._pointLevel, section->point_properties_);
     }
 
     mitochondria()._buildMitochondria(properties);
@@ -357,7 +382,7 @@ std::unordered_map<int, std::vector<unsigned int>> Morphology::connectivity() {
     return connectivity;
 }
 
-void Morphology::write(const std::string& filename) {
+void Morphology::write(const std::string& filename) const {
     for (const auto& root : rootSections()) {
         if (root->points().size() < 2) {
             throw morphio::SectionBuilderError("Root sections must have at least 2 points");


Test patch:
diff --git a/tests/test_5_mut.py b/tests/test_5_mut.py
index 131bddf3a..82e047a41 100644
--- a/tests/test_5_mut.py
+++ b/tests/test_5_mut.py
@@ -1,8 +1,9 @@
+import tempfile
 from collections import OrderedDict
 
 import numpy as np
+from numpy import testing as npt
 from pathlib import Path
-from tempfile import TemporaryDirectory
 
 import morphio
 from morphio import CellFamily, IterType, MitochondriaPointLevel, MorphioError
@@ -494,7 +495,7 @@ def test_endoplasmic_reticulum():
     assert_array_equal(reticulum.surface_areas, [3, 3])
     assert_array_equal(reticulum.filament_counts, [4, 4])
 
-    with TemporaryDirectory('test-endoplasmic-reticulum') as folder:
+    with tempfile.TemporaryDirectory('test-endoplasmic-reticulum') as folder:
         path = Path(folder, 'with-reticulum.h5')
         neuron.write(path)
 
@@ -577,7 +578,7 @@ def test_glia():
 
 
 def test_glia_round_trip():
-    with TemporaryDirectory() as folder:
+    with tempfile.TemporaryDirectory() as folder:
         g = GlialCell(Path(DATA_DIR, 'astrocyte.h5'))
         filename = Path(folder, 'glial-cell.h5')
         g.write(filename)
@@ -606,7 +607,7 @@ def test_dendritic_spine():
 
 
 def test_dendritic_spine_round_trip():
-    with TemporaryDirectory() as folder:
+    with tempfile.TemporaryDirectory() as folder:
         h5v1 = DATA_DIR / "h5/v1/"
         d = DendriticSpine(h5v1 / 'simple-dendritric-spine.h5')
         filename = Path(folder, 'test-dendritic-spine.h5')
@@ -635,7 +636,7 @@ def test_dendritic_spine_round_trip():
 
 def test_dendritic_spine_round_trip_empty_postsynaptic_density():
     h5v1 = DATA_DIR / "h5/v1/"
-    with TemporaryDirectory() as folder:
+    with tempfile.TemporaryDirectory() as folder:
         d = DendriticSpine(h5v1 / 'simple-dendritric-spine.h5')
         assert d.cell_family == CellFamily.SPINE
 
@@ -702,3 +703,356 @@ def test_lifetime_iteration_fails_with_orphan_section():
     for iter_type in IterType.depth_first, IterType.breadth_first, IterType.upstream:
         with pytest.raises(RuntimeError):
             section.iter(iter_type)
+
+
+def _assert_sections_equal(morph1, morph2, section_ids):
+
+    for section_id in section_ids:
+
+        s1 = morph1.section(section_id)
+        s2 = morph2.section(section_id)
+
+        assert s1.type == s2.type
+        npt.assert_array_equal(s1.points, s2.points)
+        npt.assert_array_equal(s1.diameters, s2.diameters)
+
+
+def test_delete_section__maintain_children_order_non_root_section():
+    string = """
+    ("CellBody"
+      (Color Red)
+      (CellBody)
+      ( 0.1  0.1 0.0 0.1)
+      ( 0.1 -0.1 0.0 0.1)
+    )
+
+    ( (Color Cyan)
+      (Axon)
+      (0.0  0.0 0.0 2.0) ; section 0
+      (0.0 -4.0 1.0 2.0)
+      (
+        (0.0 -4.0 1.0 4.0) ; section 1
+        (0.0 -4.0 2.0 4.0)
+        (0.0 -4.0 3.0 4.0)
+        (
+            (6.0 -4.0 0.0 4.0) ; section 2
+            (7.0 -5.0 0.0 4.0)
+        |
+            (6.0 -4.0 0.0 4.0) ; section 3
+            (8.0 -4.0 0.0 4.0)
+        )
+      |
+        ( 0.0 -4.0 1.0 4.0) ; section 4
+        (-5.0 -4.0 0.0 4.0)
+      )
+    )
+    """
+    morph = ImmutableMorphology(string, "asc").as_mutable()
+
+    axon = morph.root_sections[0]
+    to_be_removed_section = axon.children[0]
+    morph.delete_section(to_be_removed_section, recursive=False)
+
+    assert [s.id for s in morph.iter()] == [0, 2, 3, 4]
+
+    # a trifurcation is made
+    assert [s.id for s in axon.children] == [2, 3, 4]
+
+    # ensure nothing has been added to the leaves
+    for child in axon.children:
+        assert child.children == []
+
+    # check that data has not been messed up
+    original = ImmutableMorphology(string, "asc").as_mutable()
+    _assert_sections_equal(morph, original, [0, 2, 3, 4])
+
+
+def test_delete_section__maintain_children_order_non_root_section_multiple():
+    string = """
+    ("CellBody"
+      (Color Red)
+      (CellBody)
+      ( 0.1  0.1 0.0 0.1)
+      ( 0.1 -0.1 0.0 0.1)
+    )
+
+    ( (Color Cyan)
+      (Axon)
+      (0.0  0.0 0.0 2.0) ; section 0
+      (0.0 -4.0 1.0 2.0)
+      (
+        (0.0 -4.0 1.0 4.0) ; section 1
+        (0.0 -4.0 2.0 4.0)
+        (0.0 -4.0 3.0 4.0)
+        (
+            (6.0 -4.0 0.0 4.0) ; section 2
+            (7.0 -5.0 0.0 4.0)
+            (
+                (8.0 -4.0 0.0 4.0) ; section 3
+                (8.0 -5.0 0.0 4.0)
+            |
+                (8.0 -4.0 0.0 4.0) ; section 4
+                (8.0 -5.0 0.0 4.0)
+            )
+        |
+            (6.0 -4.0 0.0 4.0) ; section 5
+            (8.0 -4.0 0.0 4.0)
+        )
+      |
+        ( 0.0 -4.0 1.0 4.0) ; section 6
+        (-5.0 -4.0 0.0 4.0)
+        (
+            (3.0 2.0 1.0 4.0) ; section 7
+            (3.0 2.0 0.0 4.0)
+        |
+            (2.0 3.0 1.0 4.0) ; section 8
+            (2.0 3.0 0.0 4.0)
+        )
+      )
+    )
+    """
+    morph = ImmutableMorphology(string, "asc").as_mutable()
+
+    # for reference and sanity
+    assert [s.id for s in morph.iter(iter_type=morphio.IterType.depth_first)] == [0, 1, 2, 3, 4, 5, 6, 7, 8]
+    assert [s.id for s in morph.iter(iter_type=morphio.IterType.breadth_first)] == [0, 1, 6, 2, 5, 7, 8, 3, 4]
+
+    morph.delete_section(morph.section(0), recursive=False)
+    morph.delete_section(morph.section(1), recursive=False)
+    morph.delete_section(morph.section(4), recursive=False)
+    morph.delete_section(morph.section(7), recursive=False)
+
+    assert [s.id for s in morph.iter()] == [2, 3, 5, 6, 8]
+    #assert [s.id for s in morph.iter(iter_type=morphio.IterType.breadth_first)] == [6, 2, 5, 8, 3]
+    assert [s.id for s in morph.root_sections] == [2, 5, 6]
+
+    # made unifurcations
+    assert [s.id for s in morph.section(2).children] == [3]
+    assert [s.id for s in morph.section(6).children] == [8]
+
+    # ensure nothing has been added to the leaves
+    for leaf_id in [3, 5, 8]:
+        assert morph.section(leaf_id).children == []
+
+    # check that data has not been messed up
+    original = ImmutableMorphology(string, "asc").as_mutable()
+    _assert_sections_equal(morph, original, [2, 3, 5, 6, 8])
+
+
+def test_delete_section__maintain_children_order_root_section():
+    string = """
+    ("CellBody"
+      (Color Red)
+      (CellBody)
+      ( 0.1  0.1 0.0 0.1)
+      ( 0.1 -0.1 0.0 0.1)
+    )
+
+    ( (Color Cyan)
+      (Axon)
+      (0.0  0.0 0.0 2.0) ; section 0
+      (0.0 -4.0 1.0 2.0)
+      (
+        (0.0 -4.0 1.0 4.0) ; section 1
+        (0.0 -4.0 2.0 4.0)
+        (0.0 -4.0 3.0 4.0)
+        (
+            (6.0 -4.0 0.0 4.0) ; section 2
+            (7.0 -5.0 0.0 4.0)
+        |
+            (6.0 -4.0 0.0 4.0) ; section 3
+            (8.0 -4.0 0.0 4.0)
+        )
+      |
+        ( 0.0 -4.0 1.0 4.0) ; section 4
+        (-5.0 -4.0 0.0 4.0)
+      )
+    )
+
+    ( (Color Cyan)
+      (Dendrite)
+      (1.0  1.0 1.0 2.0) ; section 5
+      (1.0  1.0 2.0 2.0)
+    )
+    """
+    morph = ImmutableMorphology(string, "asc").as_mutable()
+
+    root = morph.root_sections[0]
+    morph.delete_section(root, recursive=False)
+
+    assert [s.id for s in morph.iter()] == [1, 2, 3, 4, 5]
+
+    # the children of root 0 are put in its place
+    assert [s.id for s in morph.root_sections] == [1, 4, 5]
+    assert [s.id for s in morph.root_sections[0].children] == [2, 3]
+    assert [s.id for s in morph.root_sections[1].children] == []
+    assert [s.id for s in morph.root_sections[2].children] == []
+
+    # check that nothing has been added to the leaves
+    for section_id in [2, 3, 4, 5]:
+        assert morph.section(section_id).children == []
+
+    # check that data has not been messed up
+    original = ImmutableMorphology(string, "asc").as_mutable()
+    _assert_sections_equal(morph, original, [1, 2, 3, 4, 5])
+
+
+def test_delete_section__unifurcation_bridge():
+    pl = morphio.PointLevel(
+        [[1., 2., 3.], [1., 2., 3.]],
+        [0.1, 0.1],
+    )
+
+    basal = morphio.SectionType.basal_dendrite
+    axon =  morphio.SectionType.axon
+
+    morph = morphio.mut.Morphology()
+
+    s0 = morph.append_root_section(pl, basal) # 0
+    s1 = morph.append_root_section(pl, axon) # 1
+    s2 = morph.append_root_section(pl, basal) # 2
+
+    # a unifurcation, only one child
+    s3 = s1.append_section(pl, axon) # 3
+    s4 = s3.append_section(pl, axon) # 4
+
+    s4.append_section(pl, axon) # 5
+    s4.append_section(pl, axon) # 6
+
+    assert [s.id for s in morph.iter()] == [0, 1, 3, 4, 5, 6, 2]
+
+    morph.delete_section(s3, recursive=False)
+
+    assert [s.id for s in morph.iter()] == [0, 1, 4, 5, 6, 2]
+    assert [s.parent.id for s in s1.children] == [1]
+    assert [s.id for s in s1.children] == [4]
+
+    morph.delete_section(s4, recursive=False)
+
+    assert [s.id for s in morph.iter()] == [0, 1, 5, 6, 2]
+    assert [s.parent.id for s in s1.children] == [1, 1]
+    assert [s.id for s in s1.children] == [5, 6]
+
+
+def test_delete_section__remove_all_children():
+
+    pl = morphio.PointLevel(
+        [[1., 2., 3.], [1., 2., 3.]],
+        [0.1, 0.1],
+    )
+
+    basal = morphio.SectionType.basal_dendrite
+
+    morph = morphio.mut.Morphology()
+
+    s0 = morph.append_root_section(pl, basal) # 0
+    s1 = s0.append_section(pl, basal) # 1
+    s2 = s0.append_section(pl, basal) # 2
+
+    morph.delete_section(s1, recursive=False)
+    morph.delete_section(s2, recursive=False)
+
+    assert [s.id for s in morph.iter()] == [0]
+    assert not s0.children
+
+    s3 = s0.append_section(pl, basal) # 3
+    s4 = s0.append_section(pl, basal) # 4
+
+    assert [s.id for s in morph.iter()] == [0, 3, 4]
+    assert [s.id for s in s0.children] == [3, 4]
+
+
+def test_delete_section__append_delete():
+
+    pl = morphio.PointLevel(
+        [[1., 2., 3.], [1., 2., 3.]],
+        [0.1, 0.1],
+    )
+
+    basal = morphio.SectionType.basal_dendrite
+    axon =  morphio.SectionType.axon
+
+    morph = morphio.mut.Morphology()
+
+    s0 = morph.append_root_section(pl, basal) # 0
+    s1 = morph.append_root_section(pl, axon) # 1
+    s2 = s1.append_section(pl, axon) # 2
+
+    s3 = s1.append_section(pl, axon) # 3
+
+    s4 = s3.append_section(pl, axon) # 4
+    s5 = s3.append_section(pl, axon) # 5
+    s6 = s3.append_section(pl, axon) # 6
+
+    morph.delete_section(s4, recursive=False)
+    s3.append_section(pl, axon) # 7
+    morph.delete_section(s3, recursive=False)
+
+    assert [s.id for s in morph.iter()] == [0, 1, 2, 5, 6, 7]
+
+    assert [s.parent.id for s in s1.children] == [1, 1, 1, 1]
+    assert [s.id for s in s1.children] == [2, 5, 6, 7]
+
+
+def test_delete_section__traversals():
+    string = """
+    ("CellBody"
+      (Color Red)
+      (CellBody)
+      ( 0.1  0.1 0.0 0.1)
+      ( 0.1 -0.1 0.0 0.1)
+    )
+
+    ( (Color Cyan)
+      (Axon)
+      (0.0  0.0 0.0 2.0) ; section 0
+      (0.0 -4.0 1.0 2.0)
+      (
+        (0.0 -4.0 1.0 4.0) ; section 1
+        (0.0 -4.0 2.0 4.0)
+        (0.0 -4.0 3.0 4.0)
+        (
+            (6.0 -4.0 0.0 4.0) ; section 2
+            (7.0 -5.0 0.0 4.0)
+        |
+            (6.0 -4.0 0.0 4.0) ; section 3
+            (8.0 -4.0 0.0 4.0)
+        )
+      |
+        ( 0.0 -4.0 1.0 4.0) ; section 4
+        (-5.0 -4.0 0.0 4.0)
+      )
+    )
+
+    ( (Color Cyan)
+      (Dendrite)
+      (1.0  1.0 1.0 2.0) ; section 5
+      (1.0  1.0 2.0 2.0)
+    )
+    """
+    morph = ImmutableMorphology(string, "asc").as_mutable()
+
+    assert [s.id for s in morph.iter(morphio.IterType.depth_first)] == [0, 1, 2, 3, 4, 5]
+    assert [s.id for s in morph.iter(morphio.IterType.breadth_first)] == [0, 5, 1, 4, 2, 3]
+    assert [s.id for s in morph.section(3).iter(morphio.IterType.upstream)] == [3, 1, 0]
+
+    morph.delete_section(morph.section(1), recursive=False)
+
+    assert [s.id for s in morph.iter(morphio.IterType.depth_first)] == [0, 2, 3, 4, 5]
+    assert [s.id for s in morph.iter(morphio.IterType.breadth_first)] == [0, 5, 2, 3, 4]
+    assert [s.id for s in morph.section(3).iter(morphio.IterType.upstream)] == [3, 0]
+
+    morph.delete_section(morph.section(3), recursive=False)
+
+    assert [s.id for s in morph.iter(morphio.IterType.depth_first)] == [0, 2, 4, 5]
+    assert [s.id for s in morph.iter(morphio.IterType.breadth_first)] == [0, 5, 2, 4]
+
+    morph.delete_section(morph.section(0), recursive=False)
+
+    assert [s.id for s in morph.iter(morphio.IterType.depth_first)] == [2, 4, 5]
+    assert [s.id for s in morph.iter(morphio.IterType.breadth_first)] == [2, 4, 5]
+
+    morph.delete_section(morph.section(4), recursive=False)
+
+    assert [s.id for s in morph.iter(morphio.IterType.depth_first)] == [2, 5]
+    assert [s.id for s in morph.iter(morphio.IterType.breadth_first)] == [2, 5]
diff --git a/tests/test_morphology_readers.cpp b/tests/test_morphology_readers.cpp
index eab572030..74b702a02 100644
--- a/tests/test_morphology_readers.cpp
+++ b/tests/test_morphology_readers.cpp
@@ -3,6 +3,7 @@
 
 #include <highfive/H5File.hpp>
 #include <morphio/dendritic_spine.h>
+#include <morphio/enums.h>
 #include <morphio/morphology.h>
 #include <morphio/mut/morphology.h>
 #include <morphio/soma.h>
@@ -42,6 +43,15 @@ TEST_CASE("LoadH5Morphology", "[morphology]") {
         REQUIRE(m.somaType() == morphio::SomaType::SOMA_UNDEFINED);
     }
 
+    {
+        // This is to cover the appendProperties perimeters line in mut/morphology.cpp,
+        // which is triggered if modifiers are used in a morphology that has perimeters
+        const morphio::Morphology m("data/h5/v1/glia.h5", morphio::enums::Option::NRN_ORDER);
+        REQUIRE(m.soma().points().size() == 2);
+        REQUIRE(m.points().size() == 2);
+        REQUIRE(m.perimeters().size() == 2);
+    }
+
     {  // file is an not a valid h5 file
         CHECK_THROWS_AS(morphio::Morphology("data/h5/non-valid.h5"), morphio::RawDataError);
     }
@@ -89,7 +99,6 @@ TEST_CASE("LoadH5Morphology", "[morphology]") {
     }
 }
 
-
 TEST_CASE("LoadH5Glia", "[morphology]") {
     {
         const morphio::Morphology m("data/h5/v1/glia.h5");
@@ -97,7 +106,6 @@ TEST_CASE("LoadH5Glia", "[morphology]") {
         REQUIRE(m.points().size() == 2);
         REQUIRE(m.perimeters().size() == 2);
     }
-
     {
         const morphio::Morphology m("data/h5/v1/glia_soma_only.h5");
         REQUIRE(m.soma().points().size() == 1);