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

71.53

/src/mut/writers.cpp

/* Copyright (c) 2013-2023, EPFL/Blue Brain Project
 *
 * SPDX-License-Identifier: Apache-2.0
 */
#include <cassert>
#include <fstream>
#include <iomanip>  // std::fixed, std::setw, std::setprecision

#include <morphio/errorMessages.h>
#include <morphio/mut/mitochondria.h>
#include <morphio/mut/morphology.h>
#include <morphio/mut/section.h>
#include <morphio/mut/writers.h>
#include <morphio/version.h>

#include <highfive/H5DataSet.hpp>
#include <highfive/H5File.hpp>
#include <highfive/H5Object.hpp>

#include "../shared_utils.hpp"

namespace {

/**
   A structure to get the base type of nested vectors
   https://stackoverflow.com/a/30960730
 **/
template <typename T>
struct base_type {
    using type = T;
};

template <typename T>
struct base_type<std::vector<T>>: base_type<T> {};

constexpr int FLOAT_PRECISION_PRINT = 9;

bool hasPerimeterData(const morphio::mut::Morphology& morpho) {
    return !morpho.rootSections().empty() && !morpho.rootSections().front()->perimeters().empty();
}

void writeLine(std::ofstream& myfile,
               int id,
               int parentId,
               morphio::SectionType type,
               const morphio::Point& point,
               morphio::floatType diameter) {
    using std::setw;

    myfile << std::to_string(id) << setw(12) << std::to_string(type) << ' ' << setw(12);
    myfile << std::fixed
#if !defined(MORPHIO_USE_DOUBLE)
           << std::setprecision(FLOAT_PRECISION_PRINT)
#endif
           << point[0] << ' ' << setw(12) << point[1] << ' ' << setw(12) << point[2] << ' '
           << setw(12) << diameter / 2 << setw(12);
    myfile << std::to_string(parentId) << '\n';
}

std::string version_string() {
    return std::string("Created by MorphIO v") + morphio::getVersionString();
}

/**
   Only skip duplicate if it has the same diameter
 **/
bool _skipDuplicate(const std::shared_ptr<morphio::mut::Section>& section) {
    return section->diameters().front() == section->parent()->diameters().back();
}

void checkSomaHasSameNumberPointsDiameters(const morphio::mut::Soma& soma) {
    const size_t n_points = soma.points().size();
    const size_t n_diameters = soma.diameters().size();

    if (n_points != n_diameters) {
        throw morphio::WriterError(morphio::readers::ErrorMessages().ERROR_VECTOR_LENGTH_MISMATCH(
            "soma points", n_points, "soma diameters", n_diameters));
    }
}

void raiseIfUnifurcations(const morphio::mut::Morphology& morph) {
    for (auto it = morph.depth_begin(); it != morph.depth_end(); ++it) {
        std::shared_ptr<morphio::mut::Section> section_ = *it;
        if (section_->isRoot()) {
            continue;
        }

        unsigned int parentId = section_->parent()->id();

        auto parent = section_->parent();
        bool isUnifurcation = parent->children().size() == 1;

        if (isUnifurcation) {
            throw morphio::WriterError(
                morphio::readers::ErrorMessages().ERROR_ONLY_CHILD_SWC_WRITER(parentId));
        }
    }
}
}  // anonymous namespace

namespace morphio {
namespace mut {
namespace writer {

void swc(const Morphology& morphology, const std::string& filename) {
    const auto& soma = morphology.soma();
    const auto& soma_points = soma->points();
    if (soma_points.empty()) {
        if (morphology.rootSections().empty()) {
            printError(Warning::WRITE_EMPTY_MORPHOLOGY,
                       readers::ErrorMessages().WARNING_WRITE_EMPTY_MORPHOLOGY());
            return;
        }
        printError(Warning::WRITE_NO_SOMA, readers::ErrorMessages().WARNING_WRITE_NO_SOMA());
    } else if (soma->type() == SOMA_UNDEFINED) {
        printError(Warning::WRITE_UNDEFINED_SOMA,
                   readers::ErrorMessages().WARNING_UNDEFINED_SOMA());
    } else if (!(soma->type() == SomaType::SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS ||
                 soma->type() == SomaType::SOMA_CYLINDERS ||
                 soma->type() == SomaType::SOMA_SINGLE_POINT)) {
        printError(Warning::SOMA_NON_CYLINDER_OR_POINT,
                   readers::ErrorMessages().WARNING_SOMA_NON_CYLINDER_OR_POINT());
    } else if (soma->type() == SomaType::SOMA_SINGLE_POINT && soma_points.size() != 1) {
        throw WriterError(readers::ErrorMessages().ERROR_SOMA_INVALID_SINGLE_POINT());
    } else if (soma->type() == SomaType::SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS &&
               soma_points.size() != 3) {
        throw WriterError(readers::ErrorMessages().ERROR_SOMA_INVALID_THREE_POINT_CYLINDER());
    }

    checkSomaHasSameNumberPointsDiameters(*soma);

    if (hasPerimeterData(morphology)) {
        throw WriterError(readers::ErrorMessages().ERROR_PERIMETER_DATA_NOT_WRITABLE());
    }

    raiseIfUnifurcations(morphology);

    if (!morphology.mitochondria().rootSections().empty()) {
        printError(Warning::MITOCHONDRIA_WRITE_NOT_SUPPORTED,
                   readers::ErrorMessages().WARNING_MITOCHONDRIA_WRITE_NOT_SUPPORTED());
    }

    std::ofstream myfile(filename);
    using std::setw;

    myfile << "# " << version_string() << '\n';
    myfile << "# index" << setw(9) << "type" << setw(10) << 'X' << setw(13) << 'Y' << setw(13)
           << 'Z' << setw(13) << "radius" << setw(13) << "parent" << '\n';

    int segmentIdOnDisk = 1;
    std::map<uint32_t, int32_t> newIds;

    const auto& soma_diameters = soma->diameters();

    if (soma->type() == SomaType::SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS) {
        const std::array<Point, 3> points = {
            soma_points[0],
            soma_points[1],
            soma_points[2],
        };

        const details::ThreePointSomaStatus status =
            details::checkNeuroMorphoSoma(points, soma_diameters[0] / 2);

        if (status != details::ThreePointSomaStatus::Conforms) {
            std::stringstream stream;
            stream << status;
            printError(Warning::SOMA_NON_CONFORM,
                       readers::ErrorMessages().WARNING_NEUROMORPHO_SOMA_NON_CONFORM(stream.str()));
        }

        writeLine(myfile, 1, -1, SECTION_SOMA, soma_points[0], soma_diameters[0]);
        writeLine(myfile, 2, 1, SECTION_SOMA, soma_points[1], soma_diameters[1]);
        writeLine(myfile, 3, 1, SECTION_SOMA, soma_points[2], soma_diameters[2]);
        segmentIdOnDisk += 3;
    } else {
        for (unsigned int i = 0; i < soma_points.size(); ++i) {
            writeLine(myfile,
                      segmentIdOnDisk,
                      i == 0 ? -1 : segmentIdOnDisk - 1,
                      SECTION_SOMA,
                      soma_points[i],
                      soma_diameters[i]);
            ++segmentIdOnDisk;
        }
    }

    for (auto it = morphology.depth_begin(); it != morphology.depth_end(); ++it) {
        const std::shared_ptr<Section>& section = *it;
        const auto& points = section->points();
        const auto& diameters = section->diameters();

        assert(!points.empty() && "Empty section");
        bool isRootSection = section->isRoot();

        // skips duplicate point for non-root sections
        unsigned int firstPoint = ((isRootSection || !_skipDuplicate(section)) ? 0 : 1);
        for (unsigned int i = firstPoint; i < points.size(); ++i) {
            int parentIdOnDisk;
            if (i > firstPoint) {
                parentIdOnDisk = segmentIdOnDisk - 1;
            } else {
                parentIdOnDisk = (isRootSection ? (soma->points().empty() ? -1 : 1)
                                                : newIds[section->parent()->id()]);
            }

            writeLine(
                myfile, segmentIdOnDisk, parentIdOnDisk, section->type(), points[i], diameters[i]);

            ++segmentIdOnDisk;
        }
        newIds[section->id()] = segmentIdOnDisk - 1;
    }
}

static void _write_asc_points(std::ofstream& myfile,
                              const Points& points,
                              const std::vector<morphio::floatType>& diameters,
                              size_t indentLevel) {
    for (unsigned int i = 0; i < points.size(); ++i) {
        myfile << std::fixed << std::setprecision(FLOAT_PRECISION_PRINT)
               << std::string(indentLevel, ' ') << '(' << points[i][0] << ' ' << points[i][1] << ' '
               << points[i][2] << ' ' << diameters[i] << ")\n";
    }
}

static void _write_asc_section(std::ofstream& myfile,
                               const std::shared_ptr<Section>& section,
                               size_t indentLevel) {
    std::string indent(indentLevel, ' ');
    _write_asc_points(myfile, section->points(), section->diameters(), indentLevel);

    if (!section->children().empty()) {
        auto children = section->children();
        size_t nChildren = children.size();
        for (unsigned int i = 0; i < nChildren; ++i) {
            myfile << indent << (i == 0 ? "(\n" : "|\n");
            _write_asc_section(myfile, children[i], indentLevel + 2);
        }
        myfile << indent << ")\n";
    }
}

void asc(const Morphology& morphology, const std::string& filename) {
    const auto& soma = morphology.soma();
    const auto& somaPoints = soma->points();

    if (soma->points().empty() && morphology.rootSections().empty()) {
        printError(Warning::WRITE_EMPTY_MORPHOLOGY,
                   readers::ErrorMessages().WARNING_WRITE_EMPTY_MORPHOLOGY());
        return;
    } else if (soma->type() == SOMA_UNDEFINED) {
        printError(Warning::WRITE_UNDEFINED_SOMA,
                   readers::ErrorMessages().WARNING_UNDEFINED_SOMA());
    } else if (soma->type() != SomaType::SOMA_SIMPLE_CONTOUR) {
        printError(Warning::SOMA_NON_CONTOUR, readers::ErrorMessages().WARNING_SOMA_NON_CONTOUR());
    } else if (somaPoints.empty()) {
        printError(Warning::WRITE_NO_SOMA, readers::ErrorMessages().WARNING_WRITE_NO_SOMA());
    } else if (somaPoints.size() < 3) {
        throw WriterError(readers::ErrorMessages().ERROR_SOMA_INVALID_CONTOUR());
    }

    checkSomaHasSameNumberPointsDiameters(*soma);

    if (hasPerimeterData(morphology)) {
        throw WriterError(readers::ErrorMessages().ERROR_PERIMETER_DATA_NOT_WRITABLE());
    }

    if (!morphology.mitochondria().rootSections().empty()) {
        printError(Warning::MITOCHONDRIA_WRITE_NOT_SUPPORTED,
                   readers::ErrorMessages().WARNING_MITOCHONDRIA_WRITE_NOT_SUPPORTED());
    }

    std::ofstream myfile(filename);

    if (!soma->points().empty()) {
        myfile << "(\"CellBody\"\n  (Color Red)\n  (CellBody)\n";
        _write_asc_points(myfile, soma->points(), soma->diameters(), 2);
        myfile << ")\n\n";
    }

    for (const std::shared_ptr<Section>& section : morphology.rootSections()) {
        const auto type = section->type();
        if (type == SECTION_AXON) {
            myfile << "( (Color Cyan)\n  (Axon)\n";
        } else if (type == SECTION_DENDRITE) {
            myfile << "( (Color Red)\n  (Dendrite)\n";
        } else if (type == SECTION_APICAL_DENDRITE) {
            myfile << "( (Color Red)\n  (Apical)\n";
        } else {
            throw WriterError(readers::ErrorMessages().ERROR_UNSUPPORTED_SECTION_TYPE(type));
        }
        _write_asc_section(myfile, section, 2);
        myfile << ")\n\n";
    }

    myfile << "; " << version_string() << '\n';
}

template <typename T>
HighFive::Attribute write_attribute(HighFive::File& file,
                                    const std::string& name,
                                    const T& version) {
    HighFive::Attribute a_version =
        file.createAttribute<typename T::value_type>(name, HighFive::DataSpace::From(version));
    a_version.write(version);
    return a_version;
}

template <typename T>
HighFive::Attribute write_attribute(HighFive::Group& group,
                                    const std::string& name,
                                    const T& version) {
    HighFive::Attribute a_version =
        group.createAttribute<typename T::value_type>(name, HighFive::DataSpace::From(version));
    a_version.write(version);
    return a_version;
}

template <typename T>
void write_dataset(HighFive::File& file, const std::string& name, const T& raw) {
    HighFive::DataSet dpoints =
        file.createDataSet<typename base_type<T>::type>(name, HighFive::DataSpace::From(raw));

    dpoints.write(raw);
}

template <typename T>
void write_dataset(HighFive::Group& file, const std::string& name, const T& raw) {
    HighFive::DataSet dpoints =
        file.createDataSet<typename base_type<T>::type>(name, HighFive::DataSpace::From(raw));

    dpoints.write(raw);
}

static void mitochondriaH5(HighFive::File& h5_file, const Mitochondria& mitochondria) {
    if (mitochondria.rootSections().empty()) {
        return;
    }

    Property::Properties properties;
    mitochondria._buildMitochondria(properties);
    auto& p = properties._mitochondriaPointLevel;
    size_t size = p._diameters.size();

    std::vector<std::vector<morphio::floatType>> points;
    std::vector<std::vector<int32_t>> structure;
    points.reserve(size);
    for (unsigned int i = 0; i < size; ++i) {
        points.push_back({static_cast<morphio::floatType>(p._sectionIds[i]),
                          p._relativePathLengths[i],
                          p._diameters[i]});
    }

    auto& s = properties._mitochondriaSectionLevel;
    structure.reserve(s._sections.size());
    for (const auto& section : s._sections) {
        structure.push_back({section[0], section[1]});
    }

    HighFive::Group g_organelles = h5_file.createGroup("organelles");
    HighFive::Group g_mitochondria = g_organelles.createGroup("mitochondria");

    write_dataset(g_mitochondria, "points", points);
    write_dataset(g_mitochondria, "structure", structure);
}


static void endoplasmicReticulumH5(HighFive::File& h5_file, const EndoplasmicReticulum& reticulum) {
    if (reticulum.sectionIndices().empty()) {
        return;
    }

    HighFive::Group g_organelles = h5_file.createGroup("organelles");
    HighFive::Group g_reticulum = g_organelles.createGroup("endoplasmic_reticulum");

    write_dataset(g_reticulum, "section_index", reticulum.sectionIndices());
    write_dataset(g_reticulum, "volume", reticulum.volumes());
    write_dataset(g_reticulum, "filament_count", reticulum.filamentCounts());
    write_dataset(g_reticulum, "surface_area", reticulum.surfaceAreas());
}

static void dendriticSpinePostSynapticDensityH5(HighFive::File& h5_file,
                                                const Property::DendriticSpine::Level& l) {
    const auto& psd = l._post_synaptic_density;

    HighFive::Group g_organelles = h5_file.createGroup("organelles");
    HighFive::Group g_postsynaptic_density = g_organelles.createGroup("postsynaptic_density");

    std::vector<morphio::Property::DendriticSpine::SectionId_t> sectionIds;
    sectionIds.reserve(psd.size());
    std::vector<morphio::Property::DendriticSpine::SegmentId_t> segmentIds;
    segmentIds.reserve(psd.size());
    std::vector<morphio::Property::DendriticSpine::Offset_t> offsets;
    offsets.reserve(psd.size());

    for (const auto& v : psd) {
        sectionIds.push_back(v.sectionId);
        segmentIds.push_back(v.segmentId);
        offsets.push_back(v.offset);
    }
    write_dataset(g_postsynaptic_density, "section_id", sectionIds);
    write_dataset(g_postsynaptic_density, "segment_id", segmentIds);
    write_dataset(g_postsynaptic_density, "offset", offsets);
}


void h5(const Morphology& morpho, const std::string& filename) {
    const auto& soma = morpho.soma();
    const auto& somaPoints = soma->points();

    if (somaPoints.empty()) {
        if (morpho.rootSections().empty()) {
            printError(Warning::WRITE_EMPTY_MORPHOLOGY,
                       readers::ErrorMessages().WARNING_WRITE_EMPTY_MORPHOLOGY());
            return;
        }
        printError(Warning::WRITE_NO_SOMA, readers::ErrorMessages().WARNING_WRITE_NO_SOMA());
    } else if (soma->type() == SOMA_UNDEFINED) {
        printError(Warning::WRITE_UNDEFINED_SOMA,
                   readers::ErrorMessages().WARNING_UNDEFINED_SOMA());
    } else if (soma->type() != SomaType::SOMA_SIMPLE_CONTOUR) {
        printError(Warning::SOMA_NON_CONTOUR, readers::ErrorMessages().WARNING_SOMA_NON_CONTOUR());
    } else if (somaPoints.size() < 3) {
        throw WriterError(readers::ErrorMessages().ERROR_SOMA_INVALID_CONTOUR());
    }

    checkSomaHasSameNumberPointsDiameters(*soma);

    HighFive::File h5_file(filename,
                           HighFive::File::ReadWrite | HighFive::File::Create |
                               HighFive::File::Truncate);

    int sectionIdOnDisk = 1;
    std::map<uint32_t, int32_t> newIds;

    std::vector<std::vector<morphio::floatType>> raw_points;
    std::vector<std::vector<int32_t>> raw_structure;
    std::vector<morphio::floatType> raw_perimeters;

    const auto& somaDiameters = soma->diameters();

    for (unsigned int i = 0; i < somaPoints.size(); ++i) {
        raw_points.push_back(
            {somaPoints[i][0], somaPoints[i][1], somaPoints[i][2], somaDiameters[i]});

        // If the morphology has some perimeter data, we need to fill some
        // perimeter dummy value in the soma range of the data structure to keep
        // the length matching
        if (hasPerimeterData(morpho)) {
            raw_perimeters.push_back(0);
        }
    }

    raw_structure.push_back({0, SECTION_SOMA, -1});
    size_t offset = 0;
    offset += morpho.soma()->points().size();

    for (auto it = morpho.depth_begin(); it != morpho.depth_end(); ++it) {
        const std::shared_ptr<Section>& section = *it;
        int parentOnDisk = (section->isRoot() ? 0 : newIds[section->parent()->id()]);

        const auto& points = section->points();
        const auto& diameters = section->diameters();
        const auto& perimeters = section->perimeters();

        const auto numberOfPoints = points.size();
        const auto numberOfPerimeters = perimeters.size();
        raw_structure.push_back({static_cast<int>(offset), section->type(), parentOnDisk});

        for (unsigned int i = 0; i < numberOfPoints; ++i) {
            raw_points.push_back({points[i][0], points[i][1], points[i][2], diameters[i]});
        }

        if (numberOfPerimeters > 0) {
            if (numberOfPerimeters != numberOfPoints) {
                throw WriterError(readers::ErrorMessages().ERROR_VECTOR_LENGTH_MISMATCH(
                    "points", numberOfPoints, "perimeters", numberOfPerimeters));
            }
            for (unsigned int i = 0; i < numberOfPerimeters; ++i) {
                raw_perimeters.push_back(perimeters[i]);
            }
        }

        newIds[section->id()] = sectionIdOnDisk++;
        offset += numberOfPoints;
    }

    write_dataset(h5_file, "/points", raw_points);
    write_dataset(h5_file, "/structure", raw_structure);

    HighFive::Group g_metadata = h5_file.createGroup("metadata");

    write_attribute(g_metadata, "version", std::array<uint32_t, 2>{1, 3});
    write_attribute(g_metadata,
                    "cell_family",
                    std::vector<uint32_t>{static_cast<uint32_t>(morpho.cellFamily())});
    write_attribute(h5_file, "comment", std::vector<std::string>{version_string()});

    if (hasPerimeterData(morpho)) {
        write_dataset(h5_file, "/perimeters", raw_perimeters);
    }

    mitochondriaH5(h5_file, morpho.mitochondria());
    endoplasmicReticulumH5(h5_file, morpho.endoplasmicReticulum());
    if (morpho.cellFamily() == SPINE) {
        dendriticSpinePostSynapticDensityH5(h5_file, morpho._dendriticSpineLevel);
    }
}

}  // end namespace writer
}  // end namespace mut
}  // end namespace morphio

1	/* Copyright (c) 2013-2023, EPFL/Blue Brain Project
2	*
3	* SPDX-License-Identifier: Apache-2.0
4	*/
5	#include <cassert>
6	#include <fstream>
7	#include <iomanip> // std::fixed, std::setw, std::setprecision
8
9	#include <morphio/errorMessages.h>
10	#include <morphio/mut/mitochondria.h>
11	#include <morphio/mut/morphology.h>
12	#include <morphio/mut/section.h>
13	#include <morphio/mut/writers.h>
14	#include <morphio/version.h>
15
16	#include <highfive/H5DataSet.hpp>
17	#include <highfive/H5File.hpp>
18	#include <highfive/H5Object.hpp>
19
20	#include "../shared_utils.hpp"
21
22	namespace {
23
24	/**
25	A structure to get the base type of nested vectors
26	https://stackoverflow.com/a/30960730
27	**/
28	template <typename T>
29	struct base_type {
30	using type = T;
31	};
32
33	template <typename T>
34	struct base_type<std::vector<T>>: base_type<T> {};
35
36	constexpr int FLOAT_PRECISION_PRINT = 9;
37
38	bool hasPerimeterData(const morphio::mut::Morphology& morpho) {	24✔
39	return !morpho.rootSections().empty() && !morpho.rootSections().front()->perimeters().empty();	24✔
40	}
41
42	void writeLine(std::ofstream& myfile,	118✔
43	int id,
44	int parentId,
45	morphio::SectionType type,
46	const morphio::Point& point,
47	morphio::floatType diameter) {
48	using std::setw;
49
50	myfile << std::to_string(id) << setw(12) << std::to_string(type) << ' ' << setw(12);	118✔
51	myfile << std::fixed	118✔
52	#if !defined(MORPHIO_USE_DOUBLE)
53	<< std::setprecision(FLOAT_PRECISION_PRINT)	118✔
54	#endif
55	<< point[0] << ' ' << setw(12) << point[1] << ' ' << setw(12) << point[2] << ' '	118✔
56	<< setw(12) << diameter / 2 << setw(12);	118✔
57	myfile << std::to_string(parentId) << '\n';	118✔
58	}	118✔
59
60	std::string version_string() {	16✔
61	return std::string("Created by MorphIO v") + morphio::getVersionString();	32✔
62	}
63
64	/**
65	Only skip duplicate if it has the same diameter
66	**/
67	bool _skipDuplicate(const std::shared_ptr<morphio::mut::Section>& section) {	8✔
68	return section->diameters().front() == section->parent()->diameters().back();	8✔
69	}
70
71	void checkSomaHasSameNumberPointsDiameters(const morphio::mut::Soma& soma) {	18✔
72	const size_t n_points = soma.points().size();	18✔
73	const size_t n_diameters = soma.diameters().size();	18✔
74
75	if (n_points != n_diameters) {	18✔
76	throw morphio::WriterError(morphio::readers::ErrorMessages().ERROR_VECTOR_LENGTH_MISMATCH(	4✔
77	"soma points", n_points, "soma diameters", n_diameters));	6✔
78	}
79	}	16✔
80
81	void raiseIfUnifurcations(const morphio::mut::Morphology& morph) {	12✔
82	for (auto it = morph.depth_begin(); it != morph.depth_end(); ++it) {	42✔
83	std::shared_ptr<morphio::mut::Section> section_ = *it;	30✔
84	if (section_->isRoot()) {	30✔
85	continue;	22✔
86	}
87
88	unsigned int parentId = section_->parent()->id();	8✔
89
90	auto parent = section_->parent();	16✔
91	bool isUnifurcation = parent->children().size() == 1;	8✔
92
93	if (isUnifurcation) {	8✔
94	throw morphio::WriterError(
95	morphio::readers::ErrorMessages().ERROR_ONLY_CHILD_SWC_WRITER(parentId));	×
96	}
97	}
98	}	12✔
99	} // anonymous namespace
100
101	namespace morphio {
102	namespace mut {
103	namespace writer {
104
105	void swc(const Morphology& morphology, const std::string& filename) {	20✔
106	const auto& soma = morphology.soma();	20✔
107	const auto& soma_points = soma->points();	20✔
108	if (soma_points.empty()) {	20✔
109	if (morphology.rootSections().empty()) {	4✔
110	printError(Warning::WRITE_EMPTY_MORPHOLOGY,	2✔
111	readers::ErrorMessages().WARNING_WRITE_EMPTY_MORPHOLOGY());	4✔
112	return;	2✔
113	}
114	printError(Warning::WRITE_NO_SOMA, readers::ErrorMessages().WARNING_WRITE_NO_SOMA());	2✔
115	} else if (soma->type() == SOMA_UNDEFINED) {	16✔
116	printError(Warning::WRITE_UNDEFINED_SOMA,	2✔
117	readers::ErrorMessages().WARNING_UNDEFINED_SOMA());	4✔
118	} else if (!(soma->type() == SomaType::SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS \|\|	22✔
119	soma->type() == SomaType::SOMA_CYLINDERS \|\|	8✔
120	soma->type() == SomaType::SOMA_SINGLE_POINT)) {	8✔
121	printError(Warning::SOMA_NON_CYLINDER_OR_POINT,	2✔
122	readers::ErrorMessages().WARNING_SOMA_NON_CYLINDER_OR_POINT());	4✔
123	} else if (soma->type() == SomaType::SOMA_SINGLE_POINT && soma_points.size() != 1) {	12✔
124	throw WriterError(readers::ErrorMessages().ERROR_SOMA_INVALID_SINGLE_POINT());	2✔
125	} else if (soma->type() == SomaType::SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS &&	16✔
126	soma_points.size() != 3) {	6✔
127	throw WriterError(readers::ErrorMessages().ERROR_SOMA_INVALID_THREE_POINT_CYLINDER());	2✔
128	}
129
130	checkSomaHasSameNumberPointsDiameters(*soma);	14✔
131
132	if (hasPerimeterData(morphology)) {	12✔
133	throw WriterError(readers::ErrorMessages().ERROR_PERIMETER_DATA_NOT_WRITABLE());	×
134	}
135
136	raiseIfUnifurcations(morphology);	12✔
137
138	if (!morphology.mitochondria().rootSections().empty()) {	12✔
NEW 139	printError(Warning::MITOCHONDRIA_WRITE_NOT_SUPPORTED,	×
NEW 140	readers::ErrorMessages().WARNING_MITOCHONDRIA_WRITE_NOT_SUPPORTED());	×
141	}
142
143	std::ofstream myfile(filename);	24✔
144	using std::setw;
145
146	myfile << "# " << version_string() << '\n';	12✔
147	myfile << "# index" << setw(9) << "type" << setw(10) << 'X' << setw(13) << 'Y' << setw(13)
148	<< 'Z' << setw(13) << "radius" << setw(13) << "parent" << '\n';	12✔
149
150	int segmentIdOnDisk = 1;	12✔
151	std::map<uint32_t, int32_t> newIds;	24✔
152
153	const auto& soma_diameters = soma->diameters();	12✔
154
155	if (soma->type() == SomaType::SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS) {	12✔
156	const std::array<Point, 3> points = {
157	soma_points[0],	4✔
158	soma_points[1],	4✔
159	soma_points[2],	4✔
160	};	4✔
161
162	const details::ThreePointSomaStatus status =
163	details::checkNeuroMorphoSoma(points, soma_diameters[0] / 2);	4✔
164
165	if (status != details::ThreePointSomaStatus::Conforms) {	4✔
166	std::stringstream stream;	2✔
167	stream << status;	2✔
168	printError(Warning::SOMA_NON_CONFORM,	2✔
169	readers::ErrorMessages().WARNING_NEUROMORPHO_SOMA_NON_CONFORM(stream.str()));	4✔
170	}
171
172	writeLine(myfile, 1, -1, SECTION_SOMA, soma_points[0], soma_diameters[0]);	4✔
173	writeLine(myfile, 2, 1, SECTION_SOMA, soma_points[1], soma_diameters[1]);	4✔
174	writeLine(myfile, 3, 1, SECTION_SOMA, soma_points[2], soma_diameters[2]);	4✔
175	segmentIdOnDisk += 3;	4✔
176	} else {
177	for (unsigned int i = 0; i < soma_points.size(); ++i) {	14✔
178	writeLine(myfile,	6✔
179	segmentIdOnDisk,
180	i == 0 ? -1 : segmentIdOnDisk - 1,
181	SECTION_SOMA,
182	soma_points[i],	6✔
183	soma_diameters[i]);	6✔
184	++segmentIdOnDisk;	6✔
185	}
186	}
187
188	for (auto it = morphology.depth_begin(); it != morphology.depth_end(); ++it) {	42✔
189	const std::shared_ptr<Section>& section = *it;	30✔
190	const auto& points = section->points();	30✔
191	const auto& diameters = section->diameters();	30✔
192
193	assert(!points.empty() && "Empty section");	30✔
194	bool isRootSection = section->isRoot();	30✔
195
196	// skips duplicate point for non-root sections
197	unsigned int firstPoint = ((isRootSection \|\| !_skipDuplicate(section)) ? 0 : 1);	30✔
198	for (unsigned int i = firstPoint; i < points.size(); ++i) {	130✔
199	int parentIdOnDisk;
200	if (i > firstPoint) {	100✔
201	parentIdOnDisk = segmentIdOnDisk - 1;	70✔
202	} else {
203	parentIdOnDisk = (isRootSection ? (soma->points().empty() ? -1 : 1)	30✔
204	: newIds[section->parent()->id()]);	30✔
205	}
206
207	writeLine(	100✔
208	myfile, segmentIdOnDisk, parentIdOnDisk, section->type(), points[i], diameters[i]);	100✔
209
210	++segmentIdOnDisk;	100✔
211	}
212	newIds[section->id()] = segmentIdOnDisk - 1;	30✔
213	}
214	}
215
216	static void _write_asc_points(std::ofstream& myfile,	14✔
217	const Points& points,
218	const std::vector<morphio::floatType>& diameters,
219	size_t indentLevel) {
220	for (unsigned int i = 0; i < points.size(); ++i) {	46✔
221	myfile << std::fixed << std::setprecision(FLOAT_PRECISION_PRINT)	32✔
222	<< std::string(indentLevel, ' ') << '(' << points[i][0] << ' ' << points[i][1] << ' '	64✔
223	<< points[i][2] << ' ' << diameters[i] << ")\n";	32✔
224	}
225	}	14✔
226
227	static void _write_asc_section(std::ofstream& myfile,	12✔
228	const std::shared_ptr<Section>& section,
229	size_t indentLevel) {
230	std::string indent(indentLevel, ' ');	24✔
231	_write_asc_points(myfile, section->points(), section->diameters(), indentLevel);	12✔
232
233	if (!section->children().empty()) {	12✔
234	auto children = section->children();	8✔
235	size_t nChildren = children.size();	4✔
236	for (unsigned int i = 0; i < nChildren; ++i) {	12✔
237	myfile << indent << (i == 0 ? "(\n" : "\|\n");	8✔
238	_write_asc_section(myfile, children[i], indentLevel + 2);	8✔
239	}
240	myfile << indent << ")\n";	4✔
241	}
242	}	12✔
243
244	void asc(const Morphology& morphology, const std::string& filename) {	2✔
245	const auto& soma = morphology.soma();	2✔
246	const auto& somaPoints = soma->points();	2✔
247
248	if (soma->points().empty() && morphology.rootSections().empty()) {	2✔
249	printError(Warning::WRITE_EMPTY_MORPHOLOGY,	×
250	readers::ErrorMessages().WARNING_WRITE_EMPTY_MORPHOLOGY());	×
251	return;	×
252	} else if (soma->type() == SOMA_UNDEFINED) {	2✔
NEW 253	printError(Warning::WRITE_UNDEFINED_SOMA,	×
NEW 254	readers::ErrorMessages().WARNING_UNDEFINED_SOMA());	×
255	} else if (soma->type() != SomaType::SOMA_SIMPLE_CONTOUR) {	2✔
UNCOV 256	printError(Warning::SOMA_NON_CONTOUR, readers::ErrorMessages().WARNING_SOMA_NON_CONTOUR());	×
257	} else if (somaPoints.empty()) {	2✔
258	printError(Warning::WRITE_NO_SOMA, readers::ErrorMessages().WARNING_WRITE_NO_SOMA());	×
259	} else if (somaPoints.size() < 3) {	2✔
260	throw WriterError(readers::ErrorMessages().ERROR_SOMA_INVALID_CONTOUR());	×
261	}
262
263	checkSomaHasSameNumberPointsDiameters(*soma);	2✔
264
265	if (hasPerimeterData(morphology)) {	2✔
266	throw WriterError(readers::ErrorMessages().ERROR_PERIMETER_DATA_NOT_WRITABLE());	×
267	}
268
269	if (!morphology.mitochondria().rootSections().empty()) {	2✔
UNCOV 270	printError(Warning::MITOCHONDRIA_WRITE_NOT_SUPPORTED,	×
271	readers::ErrorMessages().WARNING_MITOCHONDRIA_WRITE_NOT_SUPPORTED());	×
272	}
273
274	std::ofstream myfile(filename);	2✔
275
276	if (!soma->points().empty()) {	2✔
277	myfile << "(\"CellBody\"\n (Color Red)\n (CellBody)\n";	2✔
278	_write_asc_points(myfile, soma->points(), soma->diameters(), 2);	2✔
279	myfile << ")\n\n";	2✔
280	}
281
282	for (const std::shared_ptr<Section>& section : morphology.rootSections()) {	6✔
283	const auto type = section->type();	4✔
284	if (type == SECTION_AXON) {	4✔
285	myfile << "( (Color Cyan)\n (Axon)\n";	2✔
286	} else if (type == SECTION_DENDRITE) {	2✔
287	myfile << "( (Color Red)\n (Dendrite)\n";	2✔
NEW 288	} else if (type == SECTION_APICAL_DENDRITE) {	×
NEW 289	myfile << "( (Color Red)\n (Apical)\n";	×
290	} else {
NEW 291	throw WriterError(readers::ErrorMessages().ERROR_UNSUPPORTED_SECTION_TYPE(type));	×
292	}
293	_write_asc_section(myfile, section, 2);	4✔
294	myfile << ")\n\n";	4✔
295	}
296
297	myfile << "; " << version_string() << '\n';	2✔
298	}
299
300	template <typename T>
301	HighFive::Attribute write_attribute(HighFive::File& file,	2✔
302	const std::string& name,
303	const T& version) {
304	HighFive::Attribute a_version =	2✔
305	file.createAttribute<typename T::value_type>(name, HighFive::DataSpace::From(version));
306	a_version.write(version);	2✔
307	return a_version;	2✔
308	}
309
310	template <typename T>
311	HighFive::Attribute write_attribute(HighFive::Group& group,	4✔
312	const std::string& name,
313	const T& version) {
314	HighFive::Attribute a_version =	4✔
315	group.createAttribute<typename T::value_type>(name, HighFive::DataSpace::From(version));
316	a_version.write(version);	4✔
317	return a_version;	4✔
318	}
319
320	template <typename T>
321	void write_dataset(HighFive::File& file, const std::string& name, const T& raw) {	4✔
322	HighFive::DataSet dpoints =	4✔
323	file.createDataSet<typename base_type<T>::type>(name, HighFive::DataSpace::From(raw));
324
325	dpoints.write(raw);	4✔
326	}	4✔
327
328	template <typename T>
329	void write_dataset(HighFive::Group& file, const std::string& name, const T& raw) {	×
330	HighFive::DataSet dpoints =	×
331	file.createDataSet<typename base_type<T>::type>(name, HighFive::DataSpace::From(raw));
332
333	dpoints.write(raw);	×
334	}	×
335
336	static void mitochondriaH5(HighFive::File& h5_file, const Mitochondria& mitochondria) {	2✔
337	if (mitochondria.rootSections().empty()) {	2✔
338	return;	2✔
339	}
340
341	Property::Properties properties;	×
342	mitochondria._buildMitochondria(properties);	×
343	auto& p = properties._mitochondriaPointLevel;	×
344	size_t size = p._diameters.size();	×
345
346	std::vector<std::vector<morphio::floatType>> points;	×
347	std::vector<std::vector<int32_t>> structure;	×
348	points.reserve(size);	×
349	for (unsigned int i = 0; i < size; ++i) {	×
350	points.push_back({static_cast<morphio::floatType>(p._sectionIds[i]),	×
351	p._relativePathLengths[i],	×
352	p._diameters[i]});	×
353	}
354
355	auto& s = properties._mitochondriaSectionLevel;	×
356	structure.reserve(s._sections.size());	×
357	for (const auto& section : s._sections) {	×
358	structure.push_back({section[0], section[1]});	×
359	}
360
361	HighFive::Group g_organelles = h5_file.createGroup("organelles");	×
362	HighFive::Group g_mitochondria = g_organelles.createGroup("mitochondria");	×
363
364	write_dataset(g_mitochondria, "points", points);	×
365	write_dataset(g_mitochondria, "structure", structure);	×
366	}
367
368
369	static void endoplasmicReticulumH5(HighFive::File& h5_file, const EndoplasmicReticulum& reticulum) {	2✔
370	if (reticulum.sectionIndices().empty()) {	2✔
371	return;	2✔
372	}
373
374	HighFive::Group g_organelles = h5_file.createGroup("organelles");	×
375	HighFive::Group g_reticulum = g_organelles.createGroup("endoplasmic_reticulum");	×
376
377	write_dataset(g_reticulum, "section_index", reticulum.sectionIndices());	×
378	write_dataset(g_reticulum, "volume", reticulum.volumes());	×
379	write_dataset(g_reticulum, "filament_count", reticulum.filamentCounts());	×
380	write_dataset(g_reticulum, "surface_area", reticulum.surfaceAreas());	×
381	}
382
383	static void dendriticSpinePostSynapticDensityH5(HighFive::File& h5_file,	×
384	const Property::DendriticSpine::Level& l) {
385	const auto& psd = l._post_synaptic_density;	×
386
387	HighFive::Group g_organelles = h5_file.createGroup("organelles");	×
388	HighFive::Group g_postsynaptic_density = g_organelles.createGroup("postsynaptic_density");	×
389
390	std::vector<morphio::Property::DendriticSpine::SectionId_t> sectionIds;	×
391	sectionIds.reserve(psd.size());	×
392	std::vector<morphio::Property::DendriticSpine::SegmentId_t> segmentIds;	×
393	segmentIds.reserve(psd.size());	×
394	std::vector<morphio::Property::DendriticSpine::Offset_t> offsets;	×
395	offsets.reserve(psd.size());	×
396
397	for (const auto& v : psd) {	×
398	sectionIds.push_back(v.sectionId);	×
399	segmentIds.push_back(v.segmentId);	×
400	offsets.push_back(v.offset);	×
401	}
402	write_dataset(g_postsynaptic_density, "section_id", sectionIds);	×
403	write_dataset(g_postsynaptic_density, "segment_id", segmentIds);	×
404	write_dataset(g_postsynaptic_density, "offset", offsets);	×
405	}	×
406
407
408	void h5(const Morphology& morpho, const std::string& filename) {	2✔
409	const auto& soma = morpho.soma();	2✔
410	const auto& somaPoints = soma->points();	2✔
411
412	if (somaPoints.empty()) {	2✔
413	if (morpho.rootSections().empty()) {	×
414	printError(Warning::WRITE_EMPTY_MORPHOLOGY,	×
415	readers::ErrorMessages().WARNING_WRITE_EMPTY_MORPHOLOGY());	×
416	return;	×
417	}
418	printError(Warning::WRITE_NO_SOMA, readers::ErrorMessages().WARNING_WRITE_NO_SOMA());	×
419	} else if (soma->type() == SOMA_UNDEFINED) {	2✔
NEW 420	printError(Warning::WRITE_UNDEFINED_SOMA,	×
NEW 421	readers::ErrorMessages().WARNING_UNDEFINED_SOMA());	×
422	} else if (soma->type() != SomaType::SOMA_SIMPLE_CONTOUR) {	2✔
UNCOV 423	printError(Warning::SOMA_NON_CONTOUR, readers::ErrorMessages().WARNING_SOMA_NON_CONTOUR());	×
424	} else if (somaPoints.size() < 3) {	2✔
425	throw WriterError(readers::ErrorMessages().ERROR_SOMA_INVALID_CONTOUR());	×
426	}
427
428	checkSomaHasSameNumberPointsDiameters(*soma);	2✔
429
430	HighFive::File h5_file(filename,
431	HighFive::File::ReadWrite \| HighFive::File::Create \|
432	HighFive::File::Truncate);	4✔
433
434	int sectionIdOnDisk = 1;	2✔
435	std::map<uint32_t, int32_t> newIds;	4✔
436
437	std::vector<std::vector<morphio::floatType>> raw_points;	4✔
438	std::vector<std::vector<int32_t>> raw_structure;	4✔
439	std::vector<morphio::floatType> raw_perimeters;	4✔
440
441	const auto& somaDiameters = soma->diameters();	2✔
442
443	for (unsigned int i = 0; i < somaPoints.size(); ++i) {	10✔
444	raw_points.push_back(	8✔
445	{somaPoints[i][0], somaPoints[i][1], somaPoints[i][2], somaDiameters[i]});	8✔
446
447	// If the morphology has some perimeter data, we need to fill some
448	// perimeter dummy value in the soma range of the data structure to keep
449	// the length matching
450	if (hasPerimeterData(morpho)) {	8✔
451	raw_perimeters.push_back(0);	×
452	}
453	}
454
455	raw_structure.push_back({0, SECTION_SOMA, -1});	2✔
456	size_t offset = 0;	2✔
457	offset += morpho.soma()->points().size();	2✔
458
459	for (auto it = morpho.depth_begin(); it != morpho.depth_end(); ++it) {	14✔
460	const std::shared_ptr<Section>& section = *it;	12✔
461	int parentOnDisk = (section->isRoot() ? 0 : newIds[section->parent()->id()]);	12✔
462
463	const auto& points = section->points();	12✔
464	const auto& diameters = section->diameters();	12✔
465	const auto& perimeters = section->perimeters();	12✔
466
467	const auto numberOfPoints = points.size();	12✔
468	const auto numberOfPerimeters = perimeters.size();	12✔
469	raw_structure.push_back({static_cast<int>(offset), section->type(), parentOnDisk});	12✔
470
471	for (unsigned int i = 0; i < numberOfPoints; ++i) {	36✔
472	raw_points.push_back({points[i][0], points[i][1], points[i][2], diameters[i]});	24✔
473	}
474
475	if (numberOfPerimeters > 0) {	12✔
476	if (numberOfPerimeters != numberOfPoints) {	×
477	throw WriterError(readers::ErrorMessages().ERROR_VECTOR_LENGTH_MISMATCH(	×
478	"points", numberOfPoints, "perimeters", numberOfPerimeters));	×
479	}
480	for (unsigned int i = 0; i < numberOfPerimeters; ++i) {	×
481	raw_perimeters.push_back(perimeters[i]);	×
482	}
483	}
484
485	newIds[section->id()] = sectionIdOnDisk++;	12✔
486	offset += numberOfPoints;	12✔
487	}
488
489	write_dataset(h5_file, "/points", raw_points);	2✔
490	write_dataset(h5_file, "/structure", raw_structure);	2✔
491
492	HighFive::Group g_metadata = h5_file.createGroup("metadata");	6✔
493
494	write_attribute(g_metadata, "version", std::array<uint32_t, 2>{1, 3});	2✔
495	write_attribute(g_metadata,	2✔
496	"cell_family",
497	std::vector<uint32_t>{static_cast<uint32_t>(morpho.cellFamily())});	4✔
498	write_attribute(h5_file, "comment", std::vector<std::string>{version_string()});	4✔
499
500	if (hasPerimeterData(morpho)) {	2✔
501	write_dataset(h5_file, "/perimeters", raw_perimeters);	×
502	}
503
504	mitochondriaH5(h5_file, morpho.mitochondria());	2✔
505	endoplasmicReticulumH5(h5_file, morpho.endoplasmicReticulum());	2✔
506	if (morpho.cellFamily() == SPINE) {	2✔
507	dendriticSpinePostSynapticDensityH5(h5_file, morpho._dendriticSpineLevel);	×
508	}
509	}
510
511	} // end namespace writer
512	} // end namespace mut
513	} // end 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