Repository: BlueBrain/MorphIO (language: python)

Issue:
Vasculature section connectivity
Similarly to the neuronal morphologies, the section connectivity is also useful for vasculature datasets.

```python
from morphio.vasculature import Vasculature
v = Vasculature(filepath)
v.connectivity
```


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 055f5f31e..ba7d5cc2f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -47,3 +47,4 @@ build/
 ci/venv*
 cmake-build*
 env*
+venv-clang-format/
diff --git a/binds/python/bind_immutable.cpp b/binds/python/bind_immutable.cpp
index d30f3d0a0..e0f5f4dc5 100644
--- a/binds/python/bind_immutable.cpp
+++ b/binds/python/bind_immutable.cpp
@@ -108,14 +108,14 @@ void bind_immutable_module(py::module& m) {
         .def_property_readonly(
             "section_offsets",
             [](const morphio::Morphology& morpho) { return as_pyarray(morpho.sectionOffsets()); },
-            "Returns a list with offsets to access data of a specific section in the points\n"
-            "and diameters arrays.\n"
+            "Returns a list with N+1 offsets to access data of a specific section in the points\n"
+            "and diameters arrays (size N).\n"
             "\n"
-            "Example: accessing diameters of n'th section will be located in the DIAMETERS\n"
-            "array from DIAMETERS[sectionOffsets(n)] to DIAMETERS[sectionOffsets(n+1)-1]\n"
+            "Example: accessing diameters of n'th section will be located in the diameters\n"
+            "array: diameters[section_offsets[n]: section_offsets[n + 1]]\n"
             "\n"
-            "Note: for convenience, the last point of this array is the points array size\n"
-            "so that the above example works also for the last section.")
+            "Note: for convenience, the last element of this array is equal to the length of the\n"
+            "points/diameters arrays so that the above example works also for the last section.")
         .def_property_readonly(
             "section_types",
             [](const morphio::Morphology& morph) {
diff --git a/binds/python/bind_vasculature.cpp b/binds/python/bind_vasculature.cpp
index 601097865..365ec91b8 100644
--- a/binds/python/bind_vasculature.cpp
+++ b/binds/python/bind_vasculature.cpp
@@ -36,6 +36,21 @@ void bind_vasculature(py::module& m) {
              "throw RawDataError if the id is out of range",
              "section_id"_a)
 
+
+        .def_property_readonly(
+            "section_offsets",
+            [](morphio::vasculature::Vasculature& vasculature) {
+                return as_pyarray(vasculature.sectionOffsets());
+            },
+            "Returns a list with offsets to access data of a specific section in the points\n"
+            "and diameters arrays.\n"
+            "\n"
+            "Example: accessing diameters of n'th section will be located in the DIAMETERS\n"
+            "array from DIAMETERS[sectionOffsets(n)] to DIAMETERS[sectionOffsets(n+1)-1]\n"
+            "\n"
+            "Note: for convenience, the last point of this array is the points array size\n"
+            "so that the above example works also for the last section.")
+
         // Property accessors
         .def_property_readonly(
             "points",
@@ -59,6 +74,14 @@ void bind_vasculature(py::module& m) {
             },
             "Returns a vector with the section type of every section")
 
+        .def_property_readonly(
+            "section_connectivity",
+            [](morphio::vasculature::Vasculature* morpho) {
+                return py::array(static_cast<py::ssize_t>(morpho->sectionConnectivity().size()),
+                                 morpho->sectionConnectivity().data());
+            },
+            "Returns a 2D array of the section connectivity")
+
         // Iterators
         .def(
             "iter",
diff --git a/include/morphio/vasc/vasculature.h b/include/morphio/vasc/vasculature.h
index d38d36369..9061b5f32 100644
--- a/include/morphio/vasc/vasculature.h
+++ b/include/morphio/vasc/vasculature.h
@@ -47,6 +47,18 @@ class Vasculature
      */
     Section section(const uint32_t& id) const;
 
+    /**
+     * Returns a list with offsets to access data of a specific section in the points
+     * and diameters arrays.
+     p
+     * Example: accessing diameters of n'th section will be located in the Vasculature::diameters
+     * array from diameters[sectionOffsets(n)] to diameters[sectionOffsets(n+1)-1]
+     *
+     * Note: for convenience, the last point of this array is the points() array size
+     * so that the above example works also for the last section.
+     */
+    const std::vector<uint32_t> sectionOffsets() const noexcept;
+
     /**
      * Return a vector with all points from all sections
      **/
@@ -62,6 +74,13 @@ class Vasculature
      **/
     inline const std::vector<property::SectionType::Type>& sectionTypes() const noexcept;
 
+    /**
+     * Return a vector with all the connections between sections
+     **/
+    const std::vector<morphio::vasculature::property::Connection::Type>& sectionConnectivity() const
+        noexcept;
+
+
     /** graph iterator pointing to the begin */
     graph_iterator begin() const;
     /** graph iterator pointing to the end */
diff --git a/src/vasc/vasculature.cpp b/src/vasc/vasculature.cpp
index 213e0f52b..9c72124e4 100644
--- a/src/vasc/vasculature.cpp
+++ b/src/vasc/vasculature.cpp
@@ -55,6 +55,25 @@ std::vector<Section> Vasculature::sections() const {
     return sections_;
 }
 
+const std::vector<uint32_t> Vasculature::sectionOffsets() const noexcept {
+    // Vasculature section property is a single value representing the offset
+    const auto& offsets = _properties->get<property::VascSection>();
+
+    const auto size = offsets.size();
+    std::vector<uint32_t> indices(size + 1);
+
+    std::copy(offsets.begin(), offsets.end(), indices.begin());
+
+    indices[size] = static_cast<uint32_t>(this->points().size());
+
+    return indices;
+}
+
+const std::vector<morphio::vasculature::property::Connection::Type>&
+Vasculature::sectionConnectivity() const noexcept {
+    return _properties->get<property::Connection>();
+}
+
 graph_iterator Vasculature::begin() const {
     return graph_iterator(*this);
 }


Test patch:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index b3966905d..c74ed1c1c 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -3,6 +3,7 @@ set(TESTS_SRC
         test_morphology_readers.cpp
         test_immutable_morphology.cpp
         test_mutable_morphology.cpp
+        test_vasculature_morphology.cpp
         )
 set(TESTS_LINK_LIBRAIRIES morphio_static HighFive)
 
diff --git a/tests/test_10_vasculature.py b/tests/test_10_vasculature.py
index 8deb7a026..fabcb5737 100644
--- a/tests/test_10_vasculature.py
+++ b/tests/test_10_vasculature.py
@@ -1,6 +1,8 @@
 import os
+import h5py
 
 import numpy as np
+from numpy import testing as npt
 import pytest
 from numpy.testing import assert_array_almost_equal, assert_array_equal
 from pathlib import Path
@@ -63,6 +65,65 @@ def test_section_types():
         vasculature.Vasculature(os.path.join(_path, "h5/vasculature-broken-section-type.h5"))
 
 
+def test_section_offsets():
+    """
+    Example:
+
+    Array of points (listing ids instead of coordinates)
+    [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
+
+    And three sections:
+
+    S1: [0, 1, 2]
+    S2: [3, 4, 5]
+    S3: [6, 7, 8, 9]
+
+    The offsets will be:
+    [0, 3, 6, 10]
+
+    So the last entry is equal to the size of the points dataset. This allows for applying a diff
+    for instance to calculate the section length (in terms of points) for each section without
+    having to allocate an extra array to add that extra value.
+    """
+    morphology = vasculature.Vasculature(os.path.join(_path, "h5/vasculature1.h5"))
+
+    offset = 0
+    expected_offsets = [offset]
+    expected_lengths = []
+
+    for sec in morphology.sections:
+
+        ps = sec.points
+        offset += len(ps)
+        expected_lengths.append(len(ps))
+        expected_offsets.append(offset)
+
+    actual_section_offsets = morphology.section_offsets
+    npt.assert_allclose(actual_section_offsets, expected_offsets)
+
+    # last entry in the section offsets should be equal to the total
+    # number of points`
+    npt.assert_equal(actual_section_offsets[-1], len(morphology.points))
+
+    # and from the offsets we calculate the number of points per section
+    section_lengths = np.diff(actual_section_offsets)
+    npt.assert_allclose(section_lengths, expected_lengths)
+
+
+def test_section_connectivity():
+
+    path = os.path.join(_path, "h5/vasculature1.h5")
+    morphology = vasculature.Vasculature(path)
+
+    with h5py.File(path, 'r') as fd:
+        expected_connectivity = fd['connectivity'][:]
+
+    assert_array_almost_equal(
+        morphology.section_connectivity,
+        expected_connectivity
+    )
+
+
 def test_iterators_vasculature():
     morphology = vasculature.Vasculature(os.path.join(_path, "h5/vasculature1.h5"))
     assert_array_equal([sec.id for sec in morphology.sections], range(3080))
diff --git a/tests/test_vasculature_morphology.cpp b/tests/test_vasculature_morphology.cpp
new file mode 100644
index 000000000..deacbab80
--- /dev/null
+++ b/tests/test_vasculature_morphology.cpp
@@ -0,0 +1,63 @@
+#include "contrib/catch.hpp"
+
+#include <highfive/H5File.hpp>
+#include <morphio/vasc/section.h>
+#include <morphio/vasc/vasculature.h>
+
+
+using Vasculature = morphio::vasculature::Vasculature;
+
+namespace {
+
+struct Files {
+    std::string vasculature;
+
+    Files()
+        : vasculature("data/h5/vasculature1.h5") {}
+};
+
+}  // anonymous namespace
+
+
+TEST_CASE("vasculature_section_offsets", "[vasculature]") {
+    Files files;
+    morphio::vasculature::Vasculature morph(files.vasculature);
+
+    const auto& sections = morph.sections();
+    size_t n_offsets = sections.size() + 1;
+
+    std::vector<uint32_t> expected_section_offsets;
+    expected_section_offsets.reserve(n_offsets);
+
+    uint32_t offset = 0;
+    expected_section_offsets.push_back(offset);
+
+    for (const auto& section : sections) {
+        offset += static_cast<uint32_t>(section.points().size());
+        expected_section_offsets.push_back(offset);
+    }
+
+    REQUIRE(morph.sectionOffsets() == expected_section_offsets);
+}
+
+
+TEST_CASE("vasculature_section_connectivity", "[vasculature]") {
+    Files files;
+    morphio::vasculature::Vasculature morph(files.vasculature);
+
+    auto file = HighFive::File(files.vasculature, HighFive::File::ReadOnly);
+
+    const auto dset = file.getDataSet("/connectivity");
+    const auto dims = dset.getSpace().getDimensions();
+    const size_t n_connections = dims[0];
+
+    const auto section_connectivity = morph.sectionConnectivity();
+
+    REQUIRE(section_connectivity.size() == n_connections);
+
+    std::vector<std::array<unsigned int, 2>> expected_connectivity(n_connections);
+
+    dset.read(expected_connectivity.front().data());
+
+    REQUIRE(section_connectivity == expected_connectivity);
+}