6533620342

Committed 16 Oct 2023 12:32PM UTC coverage: 76.051% (-0.05%) from 76.104%

Build # 6533620342

Build Type

Pull #476

github

Committed by

mgeplf

Commit Message

should apt-get update first

Pull Request Pull Request #476: Efficient swc build

Run Details

278 of 278 new or added lines in 4 files covered. (100.0%)

1972 of 2593 relevant lines covered (76.05%)

903.23 hits per line

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

66.54

/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>

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() && morphology.rootSections().empty()) {
        printError(Warning::WRITE_EMPTY_MORPHOLOGY,
                   readers::ErrorMessages().WARNING_WRITE_EMPTY_MORPHOLOGY());
        return;
    } 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);

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

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

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

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

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

    if (soma_points.empty()) {
        printError(Warning::WRITE_NO_SOMA, readers::ErrorMessages().WARNING_WRITE_NO_SOMA());
    }

    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 Morphology& morpho,
                               const std::shared_ptr<Section>& section,
                               const std::map<morphio::SectionType, std::string>& header,
                               size_t indentLevel) {
    // allowed types are only the ones available in the header
    if (header.count(section->type()) == 0) {
        throw WriterError(readers::ErrorMessages().ERROR_UNSUPPORTED_SECTION_TYPE(section->type()));
    }

    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, morpho, children[i], header, indentLevel + 2);
        }
        myfile << indent << ")\n";
    }
}

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

    for (const auto& root : morphology.rootSections()) {
        if (root->points().size() < 2) {
            throw morphio::SectionBuilderError("Root sections must have at least 2 points");
        }
    }


    if (soma->points().empty() && morphology.rootSections().empty()) {
        printError(Warning::WRITE_EMPTY_MORPHOLOGY,
                   readers::ErrorMessages().WARNING_WRITE_EMPTY_MORPHOLOGY());
        return;
    } 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());
    }

    std::ofstream myfile(filename);

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

    std::map<morphio::SectionType, std::string> header;
    header[SECTION_AXON] = "( (Color Cyan)\n  (Axon)\n";
    header[SECTION_DENDRITE] = "( (Color Red)\n  (Dendrite)\n";
    header[SECTION_APICAL_DENDRITE] = "( (Color Red)\n  (Apical)\n";

    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()) {
        // allowed types are only the ones available in the header
        if (header.count(section->type()) == 0) {
            throw WriterError(
                readers::ErrorMessages().ERROR_UNSUPPORTED_SECTION_TYPE(section->type()));
        }
        myfile << header.at(section->type());
        _write_asc_section(myfile, morphology, section, header, 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();

    for (const auto& root : morpho.rootSections()) {
        if (root->points().size() < 2) {
            throw morphio::SectionBuilderError("Root sections must have at least 2 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() != 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();

    bool hasPerimeterData_ = hasPerimeterData(morpho);

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

BlueBrain / MorphIO / 6533620342

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