7111656197

Committed 06 Dec 2023 07:50AM UTC coverage: 76.555% (+0.5%) from 76.104%

Build # 7111656197

Build Type

push

github

Committed by

web-flow

Commit Message

Better read/write for soma types (#411)

* Be strict about soma types: ASC/H5 are contours, SWC is single point, cylinders or neuromorpho 3 point
* when writing, an incorrect soma type will throw ERROR_UNSUPPORTED_SOMA_TYPE
* removed unused ERROR_OPENING_FILE
* changed single point h5/asc files used in tests to be 4 point contours
* For compat, we will for now print a warning if the soma_type isn't set

Run Details

103 of 138 new or added lines in 6 files covered. (74.64%)

6 existing lines in 3 files now uncovered.

1982 of 2589 relevant lines covered (76.55%)

882.46 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.93

/src/readers/morphologyHDF5.cpp

/* Copyright (c) 2013-2023, EPFL/Blue Brain Project
 *
 * SPDX-License-Identifier: Apache-2.0
 */

#include <cassert>

#include "morphologyHDF5.h"

#include <highfive/H5Utility.hpp>  // HighFive::SilenceHDF5

#include <morphio/errorMessages.h>

namespace {

constexpr size_t SECTION_START_OFFSET = 0;
constexpr size_t SECTION_TYPE = 1;
constexpr size_t SECTION_PARENT_OFFSET = 2;

constexpr int SOMA_ONLY = -1;

//{v1
const std::string _d_structure("structure");
const std::string _d_points("points");

//{ v1.1
const std::string _a_version("version");
const std::string _g_metadata("metadata");
const std::string _a_family("cell_family");
const std::string _d_perimeters("perimeters");
//} v1.1

//{ v1.2
const std::string _g_mitochondria("organelles/mitochondria");

// endoplasmic reticulum
const std::string _g_endoplasmic_reticulum("organelles/endoplasmic_reticulum");
const std::string _d_section_index("section_index");
const std::string _d_volume("volume");
const std::string _d_surface_area("surface_area");
const std::string _d_filament_count("filament_count");
// } v1.2

//{ v1.3

// Dendritic Spine
const std::string _g_postsynaptic_density("organelles/postsynaptic_density");
const std::string _d_dendritic_spine_section_id("section_id");
const std::string _d_dendritic_spine_segment_id("segment_id");
const std::string _d_dendritic_spine_offset("offset");

//}

// } v1

//{ v2
const std::string _g_v2root("neuron1");
//} v2

}  // namespace

namespace morphio {
namespace readers {
namespace h5 {

MorphologyHDF5::MorphologyHDF5(const HighFive::Group& group, const std::string& uri)
    : _group(group)
    , _uri(uri) {}

Property::Properties load(const std::string& uri) {
    try {
        std::lock_guard<std::recursive_mutex> lock(morphio::readers::h5::global_hdf5_mutex());
        HighFive::SilenceHDF5 silence;
        auto file = HighFive::File(uri, HighFive::File::ReadOnly);
        return MorphologyHDF5(file.getGroup("/"), uri).load();

    } catch (const HighFive::FileException& exc) {
        throw RawDataError("Could not open morphology file " + uri + ": " + exc.what());
    }
}

Property::Properties load(const HighFive::Group& group) {
    std::lock_guard<std::recursive_mutex> lock(morphio::readers::h5::global_hdf5_mutex());
    return MorphologyHDF5(group).load();
}

Property::Properties MorphologyHDF5::load() {
    _readMetadata();

    int firstSectionOffset = _readSections();

    _readPoints(firstSectionOffset);

    if (_properties._cellLevel.minorVersion() >= 1) {
        _readPerimeters(firstSectionOffset);

        if (_properties._cellLevel.minorVersion() >= 2) {
            _readMitochondria();
            _readEndoplasmicReticulum();
        }

        if (_properties._cellLevel.minorVersion() >= 3 &&
            _properties._cellLevel._cellFamily == CellFamily::SPINE) {
            _readDendriticSpinePostSynapticDensity();
        }
    }

    switch (_properties._somaLevel._points.size()) {
    case 0:
        _properties._cellLevel._somaType = enums::SOMA_UNDEFINED;
        break;
    case 1:
        throw RawDataError("Morphology contour with only a single point is not valid: " + _uri);
    case 2:
        _properties._cellLevel._somaType = enums::SOMA_UNDEFINED;
        break;
    default:
        _properties._cellLevel._somaType = enums::SOMA_SIMPLE_CONTOUR;
        break;
    }

    return _properties;
}

void MorphologyHDF5::_readMetadata() {
    // default to h5v1.0
    uint32_t majorVersion = 1;
    uint32_t minorVersion = 0;
    _properties._cellLevel._cellFamily = CellFamily::NEURON;

    if (!_group.exist(_d_points) || !_group.exist(_d_structure)) {
        // h5v2 is deprecated, but it can be detected, throw a custom error messages if it is
        if (_group.exist(_g_v2root)) {
            throw RawDataError(
                "Error in " + _uri +
                "\nh5v2 is no longer supported, see: https://github.com/BlueBrain/MorphIO#H5v2");
        }
        throw RawDataError("Missing " + _d_points + " or " + _d_structure +
                           " datasets, cannot load morphology without them");
    }

    // if there is metadata, perhaps it's h5v1, minor version 1, 2
    if (_group.exist(_g_metadata)) {
        const auto metadata = _group.getGroup(_g_metadata);
        if (metadata.hasAttribute(_a_version)) {
            const auto attr = metadata.getAttribute(_a_version);

            std::array<uint32_t, 2> versions = {0, 0};
            attr.read(versions);
            majorVersion = versions[0];
            minorVersion = versions[1];

            if (majorVersion == 1 &&
                (minorVersion == 1 || minorVersion == 2 || minorVersion == 3)) {
                uint32_t family;
                metadata.getAttribute(_a_family).read(family);
                _properties._cellLevel._cellFamily = static_cast<CellFamily>(family);
            } else {
                throw RawDataError("Error in " + _uri +
                                   "\nUnsupported h5 version: " + std::to_string(majorVersion) +
                                   "." + std::to_string(minorVersion) +
                                   " See "
                                   "https://bbpteam.epfl.ch/documentation/projects/"
                                   "Morphology%20Documentation/latest/"
                                   "index.html for the list of supported versions.");
            }
        } else {
            throw RawDataError("Missing " + _a_version +
                               " attribute, cannot load morphology without them");
        }
    }

    _properties._cellLevel._version = {"h5", majorVersion, minorVersion};
}

void MorphologyHDF5::_readPoints(int firstSectionOffset) {
    constexpr size_t pointColumns = 4;

    const auto pointsDataSet = _group.getDataSet(_d_points);
    const auto pointsDims = pointsDataSet.getSpace().getDimensions();
    const size_t numberPoints = pointsDims[0];

    if (pointsDims.size() != 2) {
        throw RawDataError("Opening morphology '" + _uri +
                           "': incorrect number of dimensions in 'points'.");
    } else if (pointsDims[1] != pointColumns) {
        throw RawDataError("Opening morphology '" + _uri +
                           "': incorrect number of columns for points");
    }

    std::vector<std::array<floatType, pointColumns>> hdf5Data(numberPoints);

    if (!hdf5Data.empty()) {
        pointsDataSet.read(hdf5Data.front().data());
    }

    const bool hasSoma = firstSectionOffset != 0;
    const bool hasNeurites = static_cast<size_t>(firstSectionOffset) < numberPoints;
    const size_t somaPointCount = hasNeurites ? static_cast<size_t>(firstSectionOffset)
                                              : hdf5Data.size();

    auto& somaPoints = _properties._somaLevel._points;
    auto& somaDiameters = _properties._somaLevel._diameters;

    if (hasSoma) {
        somaPoints.resize(somaPointCount);
        somaDiameters.resize(somaPointCount);

        for (size_t i = 0; i < somaPointCount; ++i) {
            const auto& p = hdf5Data[i];
            somaPoints[i] = {p[0], p[1], p[2]};
            somaDiameters[i] = p[3];
        }
    }

    auto& points = _properties.get_mut<Property::Point>();
    auto& diameters = _properties.get_mut<Property::Diameter>();

    if (hasNeurites) {
        const size_t size = (hdf5Data.size() - somaPointCount);
        points.resize(size);
        diameters.resize(size);
        for (size_t i = somaPointCount; i < hdf5Data.size(); ++i) {
            const auto& p = hdf5Data[i];
            const size_t section = i - somaPointCount;
            points[section] = {p[0], p[1], p[2]};
            diameters[section] = p[3];
        }
    }
}

int MorphologyHDF5::_readSections() {
    // Important: The code used to split the reading of the sections and types
    //            into two separate fine-grained H5 selections. This does not
    //            reduce the number of I/O operations, but increases them by
    //            forcing HDF5 + MPI-IO to read in 4-byte groups. Thus, we now
    //            read the whole dataset at once, and split it in memory.

    constexpr size_t structureV1Columns = 3;

    const auto structure = _group.getDataSet(_d_structure);
    const auto dims = structure.getSpace().getDimensions();

    if (dims.size() != 2 || dims[1] != structureV1Columns) {
        throw(RawDataError("Error reading morphologies " + _uri +
                           " bad number of dimensions in 'structure' dataspace"));
    }

    std::vector<std::array<int, structureV1Columns>> vec(dims[0]);
    if (dims[0] > 0) {
        structure.read(vec.front().data());
    }

    assert(!vec.empty());

    bool hasSoma = true;
    if (static_cast<SectionType>(vec[0][SECTION_TYPE]) != SECTION_SOMA) {
        hasSoma = false;
    } else if (vec.size() == 1) {
        return SOMA_ONLY;
    }

    const size_t firstSection = hasSoma ? 1 : 0;
    const int firstSectionOffset = vec[firstSection][SECTION_START_OFFSET];

    auto& sections = _properties.get_mut<Property::Section>();
    sections.reserve(vec.size() - firstSection);

    auto& types = _properties.get_mut<Property::SectionType>();
    types.reserve(vec.size() - firstSection);

    // The first section is skipped if it corresponds to a soma
    for (size_t i = firstSection; i < vec.size(); ++i) {
        const auto& section = vec[i];
        SectionType type = static_cast<SectionType>(section[SECTION_TYPE]);

        if (section[SECTION_TYPE] >= SECTION_OUT_OF_RANGE_START || section[SECTION_TYPE] <= 0) {
            ErrorMessages err;
            throw RawDataError(err.ERROR_UNSUPPORTED_SECTION_TYPE(0, type));
        } else if (!hasSoma && type == SECTION_SOMA) {
            throw(RawDataError("Error reading morphology " + _uri +
                               ": it has soma section that doesn't come first"));
        } else if (hasSoma && type == SECTION_SOMA) {
            throw(RawDataError("Error reading morphology " + _uri +
                               ": it has multiple soma sections"));
        }

        sections.emplace_back(
            Property::Section::Type{section[SECTION_START_OFFSET] - firstSectionOffset,
                                    section[SECTION_PARENT_OFFSET] - (hasSoma ? 1 : 0)});
        types.emplace_back(type);
    }

    return firstSectionOffset;
}

void MorphologyHDF5::_readPerimeters(int firstSectionOffset) {
    if (!(_properties._cellLevel.majorVersion() == 1 &&
          _properties._cellLevel.minorVersion() > 0)) {
        throw RawDataError("Perimeter information is available starting at v1.1");
    }

    // soma only, won't have perimeters
    if (firstSectionOffset == SOMA_ONLY) {
        return;
    }

    if (!_group.exist(_d_perimeters)) {
        if (_properties._cellLevel._cellFamily == GLIA) {
            throw RawDataError("No empty perimeters allowed for glia morphology");
        }
        return;
    }

    auto& perimeters = _properties.get_mut<Property::Perimeter>();
    _read("", _d_perimeters, 1, perimeters);
    perimeters.erase(perimeters.begin(), perimeters.begin() + firstSectionOffset);
}


template <typename T>
void MorphologyHDF5::_read(const std::string& groupName,
                           const std::string& datasetName,
                           unsigned int expectedDimension,
                           T& data) {
    if (groupName != "" && !_group.exist(groupName)) {
        throw(
            RawDataError("Reading morphology '" + _uri + "': Missing required group " + groupName));
    }
    const auto group = groupName == "" ? _group : _group.getGroup(groupName);

    if (!group.exist(datasetName)) {
        throw(RawDataError("Reading morphology '" + _uri + "': Missing required dataset " +
                           datasetName));
    }
    const HighFive::DataSet dataset = group.getDataSet(datasetName);

    const auto dims = dataset.getSpace().getDimensions();
    if (dims.size() != expectedDimension) {
        throw(RawDataError("Reading morphology '" + _uri + "': bad number of dimensions in " +
                           datasetName));
    }

    data.resize(dims[0]);
    dataset.read(data);
}

void MorphologyHDF5::_readDendriticSpinePostSynapticDensity() {
    std::vector<morphio::Property::DendriticSpine::SectionId_t> sectionIds;
    _read(_g_postsynaptic_density, _d_dendritic_spine_section_id, 1, sectionIds);

    std::vector<morphio::Property::DendriticSpine::SegmentId_t> segmentIds;
    _read(_g_postsynaptic_density, _d_dendritic_spine_segment_id, 1, segmentIds);

    std::vector<morphio::Property::DendriticSpine::Offset_t> offsets;
    _read(_g_postsynaptic_density, _d_dendritic_spine_offset, 1, offsets);

    if (sectionIds.size() != segmentIds.size() || offsets.size() != segmentIds.size()) {
        throw(RawDataError(
            "Dendritic datasets must match in size:"
            " sectionIds: " +
            std::to_string(sectionIds.size()) + " segmentIds: " +
            std::to_string(segmentIds.size()) + " offsets: " + std::to_string(offsets.size())));
    }

    auto& properties = _properties._dendriticSpineLevel._post_synaptic_density;

    properties.reserve(sectionIds.size());
    for (size_t i = 0; i < sectionIds.size(); ++i) {
        properties.push_back({sectionIds[i], segmentIds[i], offsets[i]});
    }
}

void MorphologyHDF5::_readEndoplasmicReticulum() {
    if (!_group.exist(_g_endoplasmic_reticulum)) {
        return;
    }

    _read(_g_endoplasmic_reticulum,
          _d_section_index,
          1,
          _properties._endoplasmicReticulumLevel._sectionIndices);
    _read(_g_endoplasmic_reticulum,
          _d_volume,
          1,
          _properties._endoplasmicReticulumLevel._volumes);
    _read(_g_endoplasmic_reticulum,
          _d_surface_area,
          1,
          _properties._endoplasmicReticulumLevel._surfaceAreas);
    _read(_g_endoplasmic_reticulum,
          _d_filament_count,
          1,
          _properties._endoplasmicReticulumLevel._filamentCounts);
}

void MorphologyHDF5::_readMitochondria() {
    if (!_group.exist(_g_mitochondria)) {
        return;
    }

    std::vector<std::vector<floatType>> points;
    _read(_g_mitochondria, _d_points, 2, points);

    auto& mitoSectionId = _properties.get_mut<Property::MitoNeuriteSectionId>();
    auto& pathlength = _properties.get_mut<Property::MitoPathLength>();
    auto& diameters = _properties.get_mut<Property::MitoDiameter>();
    mitoSectionId.reserve(mitoSectionId.size() + points.size());
    pathlength.reserve(pathlength.size() + points.size());
    diameters.reserve(diameters.size() + points.size());
    for (const auto& p : points) {
        mitoSectionId.push_back(static_cast<Property::MitoNeuriteSectionId::Type>(p[0]));
        pathlength.push_back(p[1]);
        diameters.push_back(p[2]);
    }

    std::vector<std::vector<int32_t>> structure;
    _read(_g_mitochondria, _d_structure, 2, structure);

    auto& mitoSection = _properties.get_mut<Property::MitoSection>();
    mitoSection.reserve(mitoSection.size() + structure.size());
    for (auto& s : structure)
        mitoSection.emplace_back(Property::MitoSection::Type{s[0], s[1]});
}

}  // namespace h5
}  // namespace readers
}  // namespace morphio

1	/* Copyright (c) 2013-2023, EPFL/Blue Brain Project
2	*
3	* SPDX-License-Identifier: Apache-2.0
4	*/
5
6	#include <cassert>
7
8	#include "morphologyHDF5.h"
9
10	#include <highfive/H5Utility.hpp> // HighFive::SilenceHDF5
11
12	#include <morphio/errorMessages.h>
13
14	namespace {
15
16	constexpr size_t SECTION_START_OFFSET = 0;
17	constexpr size_t SECTION_TYPE = 1;
18	constexpr size_t SECTION_PARENT_OFFSET = 2;
19
20	constexpr int SOMA_ONLY = -1;
21
22	//{v1
23	const std::string _d_structure("structure");
24	const std::string _d_points("points");
25
26	//{ v1.1
27	const std::string _a_version("version");
28	const std::string _g_metadata("metadata");
29	const std::string _a_family("cell_family");
30	const std::string _d_perimeters("perimeters");
31	//} v1.1
32
33	//{ v1.2
34	const std::string _g_mitochondria("organelles/mitochondria");
35
36	// endoplasmic reticulum
37	const std::string _g_endoplasmic_reticulum("organelles/endoplasmic_reticulum");
38	const std::string _d_section_index("section_index");
39	const std::string _d_volume("volume");
40	const std::string _d_surface_area("surface_area");
41	const std::string _d_filament_count("filament_count");
42	// } v1.2
43
44	//{ v1.3
45
46	// Dendritic Spine
47	const std::string _g_postsynaptic_density("organelles/postsynaptic_density");
48	const std::string _d_dendritic_spine_section_id("section_id");
49	const std::string _d_dendritic_spine_segment_id("segment_id");
50	const std::string _d_dendritic_spine_offset("offset");
51
52	//}
53
54	// } v1
55
56	//{ v2
57	const std::string _g_v2root("neuron1");
58	//} v2
59
60	} // namespace
61
62	namespace morphio {
63	namespace readers {
64	namespace h5 {
65
66	MorphologyHDF5::MorphologyHDF5(const HighFive::Group& group, const std::string& uri)	250✔
67	: _group(group)
68	, _uri(uri) {}	250✔
69
70	Property::Properties load(const std::string& uri) {	182✔
71	try {
72	std::lock_guard<std::recursive_mutex> lock(morphio::readers::h5::global_hdf5_mutex());	364✔
73	HighFive::SilenceHDF5 silence;	364✔
74	auto file = HighFive::File(uri, HighFive::File::ReadOnly);	210✔
75	return MorphologyHDF5(file.getGroup("/"), uri).load();	590✔
76
77	} catch (const HighFive::FileException& exc) {	8✔
78	throw RawDataError("Could not open morphology file " + uri + ": " + exc.what());	4✔
79	}
80	}
81
82	Property::Properties load(const HighFive::Group& group) {	70✔
83	std::lock_guard<std::recursive_mutex> lock(morphio::readers::h5::global_hdf5_mutex());	70✔
84	return MorphologyHDF5(group).load();	210✔
85	}
86
87	Property::Properties MorphologyHDF5::load() {	248✔
88	_readMetadata();	248✔
89
90	int firstSectionOffset = _readSections();	240✔
91
92	_readPoints(firstSectionOffset);	232✔
93
94	if (_properties._cellLevel.minorVersion() >= 1) {	228✔
95	_readPerimeters(firstSectionOffset);	220✔
96
97	if (_properties._cellLevel.minorVersion() >= 2) {	212✔
98	_readMitochondria();	206✔
99	_readEndoplasmicReticulum();	206✔
100	}
101
102	if (_properties._cellLevel.minorVersion() >= 3 &&	294✔
103	_properties._cellLevel._cellFamily == CellFamily::SPINE) {	82✔
104	_readDendriticSpinePostSynapticDensity();	26✔
105	}
106	}
107
108	switch (_properties._somaLevel._points.size()) {	220✔
109	case 0:	30✔
110	_properties._cellLevel._somaType = enums::SOMA_UNDEFINED;	30✔
111	break;	30✔
NEW 112	case 1:	×
NEW 113	throw RawDataError("Morphology contour with only a single point is not valid: " + _uri);	×
114	case 2:	96✔
115	_properties._cellLevel._somaType = enums::SOMA_UNDEFINED;	96✔
116	break;	96✔
117	default:	94✔
118	_properties._cellLevel._somaType = enums::SOMA_SIMPLE_CONTOUR;	94✔
119	break;	94✔
120	}
121
122	return _properties;	220✔
123	}
124
125	void MorphologyHDF5::_readMetadata() {	248✔
126	// default to h5v1.0
127	uint32_t majorVersion = 1;	248✔
128	uint32_t minorVersion = 0;	248✔
129	_properties._cellLevel._cellFamily = CellFamily::NEURON;	248✔
130
131	if (!_group.exist(_d_points) \|\| !_group.exist(_d_structure)) {	248✔
132	// h5v2 is deprecated, but it can be detected, throw a custom error messages if it is
133	if (_group.exist(_g_v2root)) {	4✔
134	throw RawDataError(
135	"Error in " + _uri +	4✔
136	"\nh5v2 is no longer supported, see: https://github.com/BlueBrain/MorphIO#H5v2");	6✔
137	}
138	throw RawDataError("Missing " + _d_points + " or " + _d_structure +	4✔
139	" datasets, cannot load morphology without them");	6✔
140	}
141
142	// if there is metadata, perhaps it's h5v1, minor version 1, 2
143	if (_group.exist(_g_metadata)) {	244✔
144	const auto metadata = _group.getGroup(_g_metadata);	452✔
145	if (metadata.hasAttribute(_a_version)) {	226✔
146	const auto attr = metadata.getAttribute(_a_version);	448✔
147
148	std::array<uint32_t, 2> versions = {0, 0};	224✔
149	attr.read(versions);	224✔
150	majorVersion = versions[0];	224✔
151	minorVersion = versions[1];	224✔
152
153	if (majorVersion == 1 &&	224✔
154	(minorVersion == 1 \|\| minorVersion == 2 \|\| minorVersion == 3)) {	224✔
155	uint32_t family;
156	metadata.getAttribute(_a_family).read(family);	222✔
157	_properties._cellLevel._cellFamily = static_cast<CellFamily>(family);	222✔
158	} else {
159	throw RawDataError("Error in " + _uri +	4✔
160	"\nUnsupported h5 version: " + std::to_string(majorVersion) +	8✔
161	"." + std::to_string(minorVersion) +	8✔
162	" See "
163	"https://bbpteam.epfl.ch/documentation/projects/"
164	"Morphology%20Documentation/latest/"
165	"index.html for the list of supported versions.");	6✔
166	}
167	} else {
168	throw RawDataError("Missing " + _a_version +	4✔
169	" attribute, cannot load morphology without them");	6✔
170	}
171	}
172
173	_properties._cellLevel._version = {"h5", majorVersion, minorVersion};	240✔
174	}	240✔
175
176	void MorphologyHDF5::_readPoints(int firstSectionOffset) {	232✔
177	constexpr size_t pointColumns = 4;	232✔
178
179	const auto pointsDataSet = _group.getDataSet(_d_points);	464✔
180	const auto pointsDims = pointsDataSet.getSpace().getDimensions();	464✔
181	const size_t numberPoints = pointsDims[0];	232✔
182
183	if (pointsDims.size() != 2) {	232✔
184	throw RawDataError("Opening morphology '" + _uri +	4✔
185	"': incorrect number of dimensions in 'points'.");	6✔
186	} else if (pointsDims[1] != pointColumns) {	230✔
187	throw RawDataError("Opening morphology '" + _uri +	4✔
188	"': incorrect number of columns for points");	6✔
189	}
190
191	std::vector<std::array<floatType, pointColumns>> hdf5Data(numberPoints);	456✔
192
193	if (!hdf5Data.empty()) {	228✔
194	pointsDataSet.read(hdf5Data.front().data());	228✔
195	}
196
197	const bool hasSoma = firstSectionOffset != 0;	228✔
198	const bool hasNeurites = static_cast<size_t>(firstSectionOffset) < numberPoints;	228✔
199	const size_t somaPointCount = hasNeurites ? static_cast<size_t>(firstSectionOffset)	228✔
200	: hdf5Data.size();	4✔
201
202	auto& somaPoints = _properties._somaLevel._points;	228✔
203	auto& somaDiameters = _properties._somaLevel._diameters;	228✔
204
205	if (hasSoma) {	228✔
206	somaPoints.resize(somaPointCount);	198✔
207	somaDiameters.resize(somaPointCount);	198✔
208
209	for (size_t i = 0; i < somaPointCount; ++i) {	88,872✔
210	const auto& p = hdf5Data[i];	88,674✔
211	somaPoints[i] = {p[0], p[1], p[2]};	88,674✔
212	somaDiameters[i] = p[3];	88,674✔
213	}
214	}
215
216	auto& points = _properties.get_mut<Property::Point>();	228✔
217	auto& diameters = _properties.get_mut<Property::Diameter>();	228✔
218
219	if (hasNeurites) {	228✔
220	const size_t size = (hdf5Data.size() - somaPointCount);	224✔
221	points.resize(size);	224✔
222	diameters.resize(size);	224✔
223	for (size_t i = somaPointCount; i < hdf5Data.size(); ++i) {	48,778✔
224	const auto& p = hdf5Data[i];	48,554✔
225	const size_t section = i - somaPointCount;	48,554✔
226	points[section] = {p[0], p[1], p[2]};	48,554✔
227	diameters[section] = p[3];	48,554✔
228	}
229	}
230	}	228✔
231
232	int MorphologyHDF5::_readSections() {	240✔
233	// Important: The code used to split the reading of the sections and types
234	// into two separate fine-grained H5 selections. This does not
235	// reduce the number of I/O operations, but increases them by
236	// forcing HDF5 + MPI-IO to read in 4-byte groups. Thus, we now
237	// read the whole dataset at once, and split it in memory.
238
239	constexpr size_t structureV1Columns = 3;	240✔
240
241	const auto structure = _group.getDataSet(_d_structure);	480✔
242	const auto dims = structure.getSpace().getDimensions();	480✔
243
244	if (dims.size() != 2 \|\| dims[1] != structureV1Columns) {	240✔
245	throw(RawDataError("Error reading morphologies " + _uri +	4✔
246	" bad number of dimensions in 'structure' dataspace"));	6✔
247	}
248
249	std::vector<std::array<int, structureV1Columns>> vec(dims[0]);	476✔
250	if (dims[0] > 0) {	238✔
251	structure.read(vec.front().data());	238✔
252	}
253
254	assert(!vec.empty());	238✔
255
256	bool hasSoma = true;	238✔
257	if (static_cast<SectionType>(vec[0][SECTION_TYPE]) != SECTION_SOMA) {	238✔
258	hasSoma = false;	34✔
259	} else if (vec.size() == 1) {	204✔
260	return SOMA_ONLY;	6✔
261	}
262
263	const size_t firstSection = hasSoma ? 1 : 0;	232✔
264	const int firstSectionOffset = vec[firstSection][SECTION_START_OFFSET];	232✔
265
266	auto& sections = _properties.get_mut<Property::Section>();	232✔
267	sections.reserve(vec.size() - firstSection);	232✔
268
269	auto& types = _properties.get_mut<Property::SectionType>();	232✔
270	types.reserve(vec.size() - firstSection);	232✔
271
272	// The first section is skipped if it corresponds to a soma
273	for (size_t i = firstSection; i < vec.size(); ++i) {	6,968✔
274	const auto& section = vec[i];	6,742✔
275	SectionType type = static_cast<SectionType>(section[SECTION_TYPE]);	6,742✔
276
277	if (section[SECTION_TYPE] >= SECTION_OUT_OF_RANGE_START \|\| section[SECTION_TYPE] <= 0) {	6,742✔
278	ErrorMessages err;	4✔
279	throw RawDataError(err.ERROR_UNSUPPORTED_SECTION_TYPE(0, type));	2✔
280	} else if (!hasSoma && type == SECTION_SOMA) {	6,740✔
281	throw(RawDataError("Error reading morphology " + _uri +	4✔
282	": it has soma section that doesn't come first"));	6✔
283	} else if (hasSoma && type == SECTION_SOMA) {	6,738✔
284	throw(RawDataError("Error reading morphology " + _uri +	4✔
285	": it has multiple soma sections"));	6✔
286	}
287
288	sections.emplace_back(	6,736✔
289	Property::Section::Type{section[SECTION_START_OFFSET] - firstSectionOffset,	6,736✔
290	section[SECTION_PARENT_OFFSET] - (hasSoma ? 1 : 0)});	6,736✔
291	types.emplace_back(type);	6,736✔
292	}
293
294	return firstSectionOffset;	226✔
295	}
296
297	void MorphologyHDF5::_readPerimeters(int firstSectionOffset) {	220✔
298	if (!(_properties._cellLevel.majorVersion() == 1 &&	440✔
299	_properties._cellLevel.minorVersion() > 0)) {	220✔
300	throw RawDataError("Perimeter information is available starting at v1.1");	×
301	}
302
303	// soma only, won't have perimeters
304	if (firstSectionOffset == SOMA_ONLY) {	220✔
305	return;	4✔
306	}
307
308	if (!_group.exist(_d_perimeters)) {	216✔
309	if (_properties._cellLevel._cellFamily == GLIA) {	114✔
310	throw RawDataError("No empty perimeters allowed for glia morphology");	4✔
311	}
312	return;	110✔
313	}
314
315	auto& perimeters = _properties.get_mut<Property::Perimeter>();	102✔
316	_read("", _d_perimeters, 1, perimeters);	110✔
317	perimeters.erase(perimeters.begin(), perimeters.begin() + firstSectionOffset);	98✔
318	}
319
320
321	template <typename T>
322	void MorphologyHDF5::_read(const std::string& groupName,	480✔
323	const std::string& datasetName,
324	unsigned int expectedDimension,
325	T& data) {
326	if (groupName != "" && !_group.exist(groupName)) {	480✔
327	throw(
328	RawDataError("Reading morphology '" + _uri + "': Missing required group " + groupName));	×
329	}
330	const auto group = groupName == "" ? _group : _group.getGroup(groupName);	960✔
331
332	if (!group.exist(datasetName)) {	480✔
333	throw(RawDataError("Reading morphology '" + _uri + "': Missing required dataset " +	×
334	datasetName));	×
335	}
336	const HighFive::DataSet dataset = group.getDataSet(datasetName);	960✔
337
338	const auto dims = dataset.getSpace().getDimensions();	484✔
339	if (dims.size() != expectedDimension) {	480✔
340	throw(RawDataError("Reading morphology '" + _uri + "': bad number of dimensions in " +	4✔
341	datasetName));	4✔
342	}
343
344	data.resize(dims[0]);	476✔
345	dataset.read(data);	476✔
346	}	476✔
347
348	void MorphologyHDF5::_readDendriticSpinePostSynapticDensity() {	26✔
349	std::vector<morphio::Property::DendriticSpine::SectionId_t> sectionIds;	52✔
350	_read(_g_postsynaptic_density, _d_dendritic_spine_section_id, 1, sectionIds);	26✔
351
352	std::vector<morphio::Property::DendriticSpine::SegmentId_t> segmentIds;	52✔
353	_read(_g_postsynaptic_density, _d_dendritic_spine_segment_id, 1, segmentIds);	26✔
354
355	std::vector<morphio::Property::DendriticSpine::Offset_t> offsets;	52✔
356	_read(_g_postsynaptic_density, _d_dendritic_spine_offset, 1, offsets);	26✔
357
358	if (sectionIds.size() != segmentIds.size() \|\| offsets.size() != segmentIds.size()) {	26✔
359	throw(RawDataError(
360	"Dendritic datasets must match in size:"
361	" sectionIds: " +	×
362	std::to_string(sectionIds.size()) + " segmentIds: " +	×
363	std::to_string(segmentIds.size()) + " offsets: " + std::to_string(offsets.size())));	×
364	}
365
366	auto& properties = _properties._dendriticSpineLevel._post_synaptic_density;	26✔
367
368	properties.reserve(sectionIds.size());	26✔
369	for (size_t i = 0; i < sectionIds.size(); ++i) {	78✔
370	properties.push_back({sectionIds[i], segmentIds[i], offsets[i]});	52✔
371	}
372	}	26✔
373
374	void MorphologyHDF5::_readEndoplasmicReticulum() {	206✔
375	if (!_group.exist(_g_endoplasmic_reticulum)) {	206✔
376	return;	180✔
377	}
378
379	_read(_g_endoplasmic_reticulum,	26✔
380	_d_section_index,
381	1,
382	_properties._endoplasmicReticulumLevel._sectionIndices);	26✔
383	_read(_g_endoplasmic_reticulum,	26✔
384	_d_volume,
385	1,
386	_properties._endoplasmicReticulumLevel._volumes);	26✔
387	_read(_g_endoplasmic_reticulum,	26✔
388	_d_surface_area,
389	1,
390	_properties._endoplasmicReticulumLevel._surfaceAreas);	26✔
391	_read(_g_endoplasmic_reticulum,	26✔
392	_d_filament_count,
393	1,
394	_properties._endoplasmicReticulumLevel._filamentCounts);	26✔
395	}
396
397	void MorphologyHDF5::_readMitochondria() {	206✔
398	if (!_group.exist(_g_mitochondria)) {	206✔
399	return;	108✔
400	}
401
402	std::vector<std::vector<floatType>> points;	196✔
403	_read(_g_mitochondria, _d_points, 2, points);	98✔
404
405	auto& mitoSectionId = _properties.get_mut<Property::MitoNeuriteSectionId>();	98✔
406	auto& pathlength = _properties.get_mut<Property::MitoPathLength>();	98✔
407	auto& diameters = _properties.get_mut<Property::MitoDiameter>();	98✔
408	mitoSectionId.reserve(mitoSectionId.size() + points.size());	98✔
409	pathlength.reserve(pathlength.size() + points.size());	98✔
410	diameters.reserve(diameters.size() + points.size());	98✔
411	for (const auto& p : points) {	1,078✔
412	mitoSectionId.push_back(static_cast<Property::MitoNeuriteSectionId::Type>(p[0]));	980✔
413	pathlength.push_back(p[1]);	980✔
414	diameters.push_back(p[2]);	980✔
415	}
416
417	std::vector<std::vector<int32_t>> structure;	196✔
418	_read(_g_mitochondria, _d_structure, 2, structure);	98✔
419
420	auto& mitoSection = _properties.get_mut<Property::MitoSection>();	98✔
421	mitoSection.reserve(mitoSection.size() + structure.size());	98✔
422	for (auto& s : structure)	392✔
423	mitoSection.emplace_back(Property::MitoSection::Type{s[0], s[1]});	294✔
424	}
425
426	} // namespace h5
427	} // namespace readers
428	} // namespace morphio

BlueBrain / MorphIO / 7111656197

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous