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

88.8

/src/readers/morphologySWC.cpp

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

#include "morphologySWC.h"
#include "utils.h"

#include <cctype>         // isdigit
#include <cstdint>        // uint32_t
#include <memory>         // std::shared_ptr
#include <string>         // std::string
#include <unordered_map>  // std::unordered_map
#include <vector>         // std::vector

#include <morphio/errorMessages.h>
#include <morphio/mut/morphology.h>
#include <morphio/mut/section.h>
#include <morphio/mut/soma.h>
#include <morphio/properties.h>

namespace details {
using namespace morphio;

// It's not clear if -1 is the only way of identifying a root section.
const int SWC_UNDEFINED_PARENT = -1;
const unsigned int SWC_ROOT = 0xFFFFFFFD;

/* simple stream parser for SWC file format which is a line oriented format
 *
 * This parser advances across comments and blank lines, and allows the caller
 * to get integers and floats
 */
class SWCTokenizer
{
public:
  explicit SWCTokenizer(std::string contents, const morphio::readers::ErrorMessages& err)
      : contents_(std::move(contents))
      , err_(err) {
      // ensure null termination
      (void) contents_.c_str();
  }

  bool done() const noexcept {
      return pos_ >= contents_.size();
  }

  size_t lineNumber() const noexcept {
      return line_;
  }

  void skip_to(char value) {
      std::size_t pos = contents_.find_first_of(value, pos_);
      if (pos == std::string::npos) {
          pos_ = contents_.size();
      }
      pos_ = pos;
  }

  void advance_to_non_whitespace() {
      std::size_t pos = contents_.find_first_not_of(" \t\r", pos_);
      if (pos == std::string::npos) {
          pos_ = contents_.size();
      }
      pos_ = pos;
  }

  void advance_to_number() {
      while (consume_line_and_trailing_comments()) {
      }

      if (done()) {
          throw morphio::RawDataError(err_.EARLY_END_OF_FILE(line_));
      }

      auto c = contents_.at(pos_);
      if (std::isdigit(c) != 0 || c == '-' || c == '+' || c == '.') {
          return;
      }

      throw morphio::RawDataError(err_.ERROR_LINE_NON_PARSABLE(line_));
  }

  int64_t read_int() {
      advance_to_number();
      auto parsed = stn_.toInt(contents_, pos_);
      pos_ = std::get<1>(parsed);
      return std::get<0>(parsed);
  }

  morphio::floatType read_float() {
      advance_to_number();
      auto parsed = stn_.toFloat(contents_, pos_);
      pos_ = std::get<1>(parsed);
      return std::get<0>(parsed);
  }

  bool consume_line_and_trailing_comments() {
      bool found_newline = false;

      advance_to_non_whitespace();
      while (!done() && (contents_.at(pos_) == '#' || contents_.at(pos_) == '\n')) {
          switch (contents_.at(pos_)) {
          case '#':
              skip_to('\n');
              break;
          case '\n':
              ++line_;
              ++pos_;
              found_newline = true;
              break;
          }
          advance_to_non_whitespace();
      }
      return found_newline || done();
  }

private:
  size_t pos_ = 0;
  size_t line_ = 1;
  std::string contents_;
  morphio::StringToNumber stn_{};
  morphio::readers::ErrorMessages err_;
};

std::vector<morphio::readers::Sample> readSamples(const std::string& contents,
                                                  const morphio::readers::ErrorMessages& err) {
    std::vector<morphio::readers::Sample> samples;
    morphio::readers::Sample sample;

    SWCTokenizer tokenizer{contents, err};
    tokenizer.consume_line_and_trailing_comments();

    while (!tokenizer.done()) {
        sample.lineNumber = static_cast<unsigned int>(tokenizer.lineNumber());

        int64_t id = tokenizer.read_int();
        if (id < 0) {
            throw morphio::RawDataError(err.ERROR_NEGATIVE_ID(sample.lineNumber));
        }

        sample.id = static_cast<unsigned int>(id);

        sample.type = static_cast<morphio::SectionType>(tokenizer.read_int());

        for (auto& point : sample.point) {
            point = tokenizer.read_float();
        }

        sample.diameter = 2 * tokenizer.read_float();

        int64_t parentId = tokenizer.read_int();
        if (parentId < -1) {
            throw morphio::RawDataError(err.ERROR_NEGATIVE_ID(sample.lineNumber));
        } else if (parentId == SWC_UNDEFINED_PARENT) {
            sample.parentId = SWC_ROOT;
        } else {
            sample.parentId = static_cast<unsigned int>(parentId);
        }

        if (!tokenizer.consume_line_and_trailing_comments()) {
            throw morphio::RawDataError(err.ERROR_LINE_NON_PARSABLE(sample.lineNumber));
        }
        samples.push_back(sample);
    }
    return samples;
}

/**
  Parsing SWC according to this specification:
http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
 **/
using morphio::readers::ErrorMessages;
using morphio::readers::Sample;

enum class DeclaredID : unsigned int {};
unsigned int operator+(DeclaredID type) {
    return static_cast<unsigned int>(type);
}

class SWCBuilder
{
    using MorphioID = uint32_t;
    using PhysicalID = size_t;
    using Samples = std::vector<Sample>;

  public:
    explicit SWCBuilder(const std::string& path)
        : err_(path) {}

    Property::Properties buildProperties(const std::string& contents, unsigned int options) {
        const Samples samples = readSamples(contents, err_);
        buildSWC(samples);
        morph_.applyModifiers(options);
        return morph_.buildReadOnly();
    }

  private:
    void _checkNeuromorph3PointSoma(const Samples& soma_samples){
        // First point is the 'center'; has 2 children
        const Sample& center = soma_samples[0];
        const Sample& child1 = soma_samples[1];
        const Sample& child2 = soma_samples[2];

        floatType x = center.point[0];
        floatType z = center.point[2];
        floatType d = center.diameter;
        floatType r = center.diameter / 2;
        floatType y = center.point[1];

        // whether the soma should be checked for the special case of 3 point soma
        // for details see https://github.com/BlueBrain/MorphIO/issues/273
        // If the 2nd and the 3rd point have the same x, z, d values then the only valid soma
        // is: 1 1 x   y   z r -1 2 1 x (y-r) z r  1 3 1 x (y+r) z r  1
        if ((child1.point[0] != x || child2.point[0] != x || child1.point[1] != y - r ||
             child2.point[1] != y + r || child1.point[2] != z || child2.point[2] != z ||
             child1.diameter != d || child2.diameter != d) &&
            std::fabs(child1.diameter - d) < morphio::epsilon &&
            std::fabs(child2.diameter - d) < morphio::epsilon &&
            std::fabs(child1.point[0] - x) < morphio::epsilon &&
            std::fabs(child2.point[0] - x) < morphio::epsilon &&
            std::fabs(child1.point[2] - z) < morphio::epsilon &&
            std::fabs(child2.point[2] - z) < morphio::epsilon
            ) {
            printError(Warning::SOMA_NON_CONFORM,
                       err_.WARNING_NEUROMORPHO_SOMA_NON_CONFORM(center, child1, child2));
        }
    }

    void build_soma(const Samples& soma_samples) {
        auto& soma = morph_.soma();

        if (soma_samples.empty()) {
            soma->type() = SOMA_UNDEFINED;
            return;
        } else if (soma_samples.size() == 1) {
            Sample sample = soma_samples[0];

            if (sample.parentId != SWC_ROOT && samples_.at(sample.parentId).type != SECTION_SOMA) {
                throw morphio::SomaError(err_.ERROR_SOMA_WITH_NEURITE_PARENT(sample));
            }

            soma->type() = SOMA_SINGLE_POINT;
            soma->points() = {sample.point};
            soma->diameters() = {sample.diameter};
            return;
        } else if (soma_samples.size() == 3) {
            const Sample& center = soma_samples[0];
            const Sample& child1 = soma_samples[1];
            const Sample& child2 = soma_samples[2];
            // All soma that bifurcate with the first parent having two children are considered
            // SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS
            if(center.id == child1.parentId && center.id == child2.parentId){
                soma->type() = SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS;
                soma->points() = {center.point, child1.point, child2.point};
                soma->diameters() = {center.diameter, child1.diameter, child2.diameter};
                _checkNeuromorph3PointSoma(soma_samples);
               return;
           }
        }
        // might also have 3 points at this point, as well

        // a "normal" SWC soma
        soma->type() = SOMA_CYLINDERS;
        auto& points = soma->points();
        auto& diameters = soma->diameters();
        points.reserve(soma_samples.size());
        diameters.reserve(soma_samples.size());

        size_t parent_count = 0;
        for (const auto& s : soma_samples) {
            if (s.parentId == SWC_ROOT) {
                parent_count++;
            } else if (samples_.count(s.parentId) == 0) {
                throw morphio::MissingParentError(err_.ERROR_MISSING_PARENT(s));
            } else if (samples_.at(s.parentId).type != SECTION_SOMA) {
                throw morphio::SomaError(err_.ERROR_SOMA_WITH_NEURITE_PARENT(s));
            }

            if (children_.count(s.id) > 0 && children_.at(s.id).size() > 1) {
                std::vector<Sample> soma_bifurcations;
                for (auto id : children_.at(s.id)) {
                    if (samples_[id].type == SECTION_SOMA) {
                        soma_bifurcations.push_back(samples_[id]);
                    }
                }
                if (soma_bifurcations.size() > 1) {
                    throw morphio::SomaError(err_.ERROR_SOMA_BIFURCATION(s, soma_bifurcations));
                }
            }
            points.push_back(s.point);
            diameters.push_back(s.diameter);
        }

        if (parent_count > 1) {
            throw morphio::SomaError(err_.ERROR_MULTIPLE_SOMATA(soma_samples));
        }
    }

    void buildSWC(const Samples& samples) {
        Samples soma_samples;
        Samples root_samples;

        for (const auto& sample: samples) {
            // { checks
            if (sample.diameter < morphio::epsilon) {
                printError(Warning::ZERO_DIAMETER, err_.WARNING_ZERO_DIAMETER(sample));
            }

            if (sample.parentId == sample.id) {
                throw morphio::RawDataError(err_.ERROR_SELF_PARENT(sample));
            }

            if (sample.type >= morphio::SECTION_OUT_OF_RANGE_START || sample.type <= 0) {
                throw morphio::RawDataError(
                    err_.ERROR_UNSUPPORTED_SECTION_TYPE(sample.lineNumber, sample.type));
            }
            if (sample.parentId == SWC_ROOT && sample.type != SECTION_SOMA) {
                printError(Warning::DISCONNECTED_NEURITE,
                           err_.WARNING_DISCONNECTED_NEURITE(sample));
            }
            // } checks

            if (sample.type == SECTION_SOMA) {
                soma_samples.push_back(sample);
            }

            if (sample.parentId == SWC_ROOT || sample.type == SECTION_SOMA) {
                root_samples.push_back(sample);
            }

            if (!samples_.insert({sample.id, sample}).second) {
                throw RawDataError(err_.ERROR_REPEATED_ID(samples[sample.id], sample));
            }

            children_[sample.parentId].push_back(sample.id);
        }

        // can only check for missing parents once all samples are loaded
        // since it's possible there may be forward references
        for (const auto& sample: samples) {
            if(sample.parentId != SWC_ROOT && samples_.count(sample.parentId) == 0){
                throw morphio::MissingParentError(err_.ERROR_MISSING_PARENT(sample));
            }
        }

        build_soma(soma_samples);

        std::unordered_map<DeclaredID, std::shared_ptr<morphio::mut::Section>> declared_to_swc;
        declared_to_swc.reserve(samples.size());

        for (const Sample& root_sample : root_samples) {
            if (children_.count(root_sample.id) == 0) {
                continue;
            }

            // https://neuromorpho.org/SomaFormat.html
            // "The second and third soma points, as well as all starting points
            // (roots) of dendritic and axonal arbors have this first point as
            // the parent (parent ID 1)."
            if (morph_.soma()->type() == SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS &&
                root_sample.type == SECTION_SOMA && root_sample.id != 1) {
                printError(Warning::WRONG_ROOT_POINT, err_.WARNING_WRONG_ROOT_POINT(root_sample));
            }

            for (unsigned int child_id : children_.at(root_sample.id)) {
                if (samples_.at(child_id).type == SECTION_SOMA) {
                    continue;
                }
                if (root_sample.type == SECTION_SOMA) {
                    assembleSections(child_id,
                                     DeclaredID(root_sample.id),
                                     declared_to_swc,
                                     morph_.soma()->points()[0],
                                     morph_.soma()->diameters()[0],
                                     true);
                } else {
                    // this is neurite as the start
                    assembleSections(root_sample.id,
                                     DeclaredID(SWC_ROOT),
                                     declared_to_swc,
                                     root_sample.point,
                                     root_sample.diameter,
                                     true);
                    break;
                }
            }
        }
    }

    void assembleSections(
        unsigned int id,
        DeclaredID parent_id,
        std::unordered_map<DeclaredID, std::shared_ptr<morphio::mut::Section>>& declared_to_swc,
        const Point& start_point,
        floatType start_diameter,
        bool is_root) {

        morphio::Property::PointLevel properties;
        auto& points = properties._points;
        auto& diameters = properties._diameters;

        auto appendSection = [&](DeclaredID section_id_, DeclaredID parent_id_, SectionType starting_section_type) {
            std::shared_ptr<morphio::mut::Section> new_section;
            if (is_root) {
                new_section = morph_.appendRootSection(properties, starting_section_type);
            } else {
                new_section = declared_to_swc.at(parent_id_)
                                  ->appendSection(properties, starting_section_type);
            }
            declared_to_swc[section_id_] = new_section;
        };

        auto get_child_count = [&](unsigned int child_id) {
            return children_.count(child_id) == 0 ? 0 : children_.at(child_id).size();
        };

        const Sample* sample = &samples_.at(id);

        // create duplicate point if needed
        if (!is_root && sample->point != start_point /*|| sample->diameter != start_diameter */) {
        //if (!(is_root && sample->type == SECTION_SOMA) && sample->point != start_point /*|| sample->diameter != start_diameter */) {
            points.push_back(start_point);
            diameters.push_back(start_diameter);
        }

        // try and combine as many single samples into a single section as possible
        size_t children_count = get_child_count(id);
        while (children_count == 1) {
            sample = &samples_.at(id);
            if(sample->type != samples_.at(children_.at(id)[0]).type){
                break;
            }
            points.push_back(sample->point);
            diameters.push_back(sample->diameter);
            id = children_.at(id)[0];
            children_count = get_child_count(id);
        }
        sample = &samples_.at(id);
        points.push_back(sample->point);
        diameters.push_back(sample->diameter);
        appendSection(DeclaredID(id), parent_id, sample->type);

        if (children_count == 0) {
            // section was already appended above, nothing to do
        } else if (children_count == 1) {
            // section_type changed
            size_t offset = properties._points.size() - 1;
            const Point& new_start_point = properties._points[offset];
            floatType new_start_diameter = properties._diameters[offset];
            assembleSections(children_.at(id)[0], DeclaredID(id), declared_to_swc, new_start_point, new_start_diameter, false);
        } else {
            size_t offset = properties._points.size() - 1;
            const Point& new_start_point = properties._points[offset];
            floatType new_start_diameter = properties._diameters[offset];
            for (unsigned int child_id : children_.at(id)) {
                assembleSections(child_id,
                                 DeclaredID(id),
                                 declared_to_swc,
                                 new_start_point,
                                 new_start_diameter,
                                 false);
            }
        }
    }

    std::unordered_map<unsigned int, std::vector<unsigned int>> children_;
    std::unordered_map<unsigned int, Sample> samples_;
    mut::Morphology morph_;
    ErrorMessages err_;
};


}  // namespace details

namespace morphio {
namespace readers {
namespace swc {
Property::Properties load(const std::string& path,
                          const std::string& contents,
                          unsigned int options) {
    auto properties = details::SWCBuilder(path).buildProperties(contents, options);

    properties._cellLevel._cellFamily = NEURON;
    properties._cellLevel._version = {"swc", 1, 0};
    return properties;
}

}  // namespace swc
}  // 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 "morphologySWC.h"
7	#include "utils.h"
8
9	#include <cctype> // isdigit
10	#include <cstdint> // uint32_t
11	#include <memory> // std::shared_ptr
12	#include <string> // std::string
13	#include <unordered_map> // std::unordered_map
14	#include <vector> // std::vector
15
16	#include <morphio/errorMessages.h>
17	#include <morphio/mut/morphology.h>
18	#include <morphio/mut/section.h>
19	#include <morphio/mut/soma.h>
20	#include <morphio/properties.h>
21
22	namespace details {
23	using namespace morphio;
24
25	// It's not clear if -1 is the only way of identifying a root section.
26	const int SWC_UNDEFINED_PARENT = -1;
27	const unsigned int SWC_ROOT = 0xFFFFFFFD;
28
29	/* simple stream parser for SWC file format which is a line oriented format
30	*
31	* This parser advances across comments and blank lines, and allows the caller
32	* to get integers and floats
33	*/
34	class SWCTokenizer
35	{
36	public:
37	explicit SWCTokenizer(std::string contents, const morphio::readers::ErrorMessages& err)	66✔
38	: contents_(std::move(contents))	66✔
39	, err_(err) {	66✔
40	// ensure null termination
41	(void) contents_.c_str();	66✔
42	}	66✔
43
44	bool done() const noexcept {	13,728✔
45	return pos_ >= contents_.size();	13,728✔
46	}
47
48	size_t lineNumber() const noexcept {	506✔
49	return line_;	506✔
50	}
51
52	void skip_to(char value) {	684✔
53	std::size_t pos = contents_.find_first_of(value, pos_);	684✔
54	if (pos == std::string::npos) {	684✔
55	pos_ = contents_.size();	×
56	}
57	pos_ = pos;	684✔
58	}	684✔
59
60	void advance_to_non_whitespace() {	6,072✔
61	std::size_t pos = contents_.find_first_not_of(" \t\r", pos_);	6,072✔
62	if (pos == std::string::npos) {	6,072✔
63	pos_ = contents_.size();	64✔
64	}
65	pos_ = pos;	6,072✔
66	}	6,072✔
67
68	void advance_to_number() {	3,542✔
69	while (consume_line_and_trailing_comments()) {	3,542✔
70	}
71
72	if (done()) {	3,542✔
73	throw morphio::RawDataError(err_.EARLY_END_OF_FILE(line_));	×
74	}
75
76	auto c = contents_.at(pos_);	3,542✔
77	if (std::isdigit(c) != 0 \|\| c == '-' \|\| c == '+' \|\| c == '.') {	3,542✔
78	return;	3,540✔
79	}
80
81	throw morphio::RawDataError(err_.ERROR_LINE_NON_PARSABLE(line_));	2✔
82	}
83
84	int64_t read_int() {	1,518✔
85	advance_to_number();	1,518✔
86	auto parsed = stn_.toInt(contents_, pos_);	1,516✔
87	pos_ = std::get<1>(parsed);	1,516✔
88	return std::get<0>(parsed);	1,516✔
89	}
90
91	morphio::floatType read_float() {	2,024✔
92	advance_to_number();	2,024✔
93	auto parsed = stn_.toFloat(contents_, pos_);	2,024✔
94	pos_ = std::get<1>(parsed);	2,024✔
95	return std::get<0>(parsed);	2,024✔
96	}
97
98	bool consume_line_and_trailing_comments() {	4,112✔
99	bool found_newline = false;	4,112✔
100
101	advance_to_non_whitespace();	4,112✔
102	while (!done() && (contents_.at(pos_) == '#' \|\| contents_.at(pos_) == '\n')) {	6,072✔
103	switch (contents_.at(pos_)) {	1,960✔
104	case '#':	684✔
105	skip_to('\n');	684✔
106	break;	684✔
107	case '\n':	1,276✔
108	++line_;	1,276✔
109	++pos_;	1,276✔
110	found_newline = true;	1,276✔
111	break;	1,276✔
112	}
113	advance_to_non_whitespace();	1,960✔
114	}
115	return found_newline \|\| done();	4,112✔
116	}
117
118	private:
119	size_t pos_ = 0;
120	size_t line_ = 1;
121	std::string contents_;
122	morphio::StringToNumber stn_{};
123	morphio::readers::ErrorMessages err_;
124	};
125
126	std::vector<morphio::readers::Sample> readSamples(const std::string& contents,	66✔
127	const morphio::readers::ErrorMessages& err) {
128	std::vector<morphio::readers::Sample> samples;	66✔
129	morphio::readers::Sample sample;	66✔
130
131	SWCTokenizer tokenizer{contents, err};	132✔
132	tokenizer.consume_line_and_trailing_comments();	66✔
133
134	while (!tokenizer.done()) {	570✔
135	sample.lineNumber = static_cast<unsigned int>(tokenizer.lineNumber());	506✔
136
137	int64_t id = tokenizer.read_int();	506✔
138	if (id < 0) {	506✔
139	throw morphio::RawDataError(err.ERROR_NEGATIVE_ID(sample.lineNumber));	×
140	}
141
142	sample.id = static_cast<unsigned int>(id);	506✔
143
144	sample.type = static_cast<morphio::SectionType>(tokenizer.read_int());	506✔
145
146	for (auto& point : sample.point) {	2,024✔
147	point = tokenizer.read_float();	1,518✔
148	}
149
150	sample.diameter = 2 * tokenizer.read_float();	506✔
151
152	int64_t parentId = tokenizer.read_int();	506✔
153	if (parentId < -1) {	504✔
154	throw morphio::RawDataError(err.ERROR_NEGATIVE_ID(sample.lineNumber));	×
155	} else if (parentId == SWC_UNDEFINED_PARENT) {	504✔
156	sample.parentId = SWC_ROOT;	66✔
157	} else {
158	sample.parentId = static_cast<unsigned int>(parentId);	438✔
159	}
160
161	if (!tokenizer.consume_line_and_trailing_comments()) {	504✔
162	throw morphio::RawDataError(err.ERROR_LINE_NON_PARSABLE(sample.lineNumber));	×
163	}
164	samples.push_back(sample);	504✔
165	}
166	return samples;	128✔
167	}
168
169	/**
170	Parsing SWC according to this specification:
171	http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
172	**/
173	using morphio::readers::ErrorMessages;
174	using morphio::readers::Sample;
175
176	enum class DeclaredID : unsigned int {};
177	unsigned int operator+(DeclaredID type) {	×
178	return static_cast<unsigned int>(type);	×
179	}
180
181	class SWCBuilder
182	{
183	using MorphioID = uint32_t;
184	using PhysicalID = size_t;
185	using Samples = std::vector<Sample>;
186
187	public:
188	explicit SWCBuilder(const std::string& path)	66✔
189	: err_(path) {}	66✔
190
191	Property::Properties buildProperties(const std::string& contents, unsigned int options) {	66✔
192	const Samples samples = readSamples(contents, err_);	130✔
193	buildSWC(samples);	64✔
194	morph_.applyModifiers(options);	50✔
195	return morph_.buildReadOnly();	100✔
196	}
197
198	private:
199	void _checkNeuromorph3PointSoma(const Samples& soma_samples){	4✔
200	// First point is the 'center'; has 2 children
201	const Sample& center = soma_samples[0];	4✔
202	const Sample& child1 = soma_samples[1];	4✔
203	const Sample& child2 = soma_samples[2];	4✔
204
205	floatType x = center.point[0];	4✔
206	floatType z = center.point[2];	4✔
207	floatType d = center.diameter;	4✔
208	floatType r = center.diameter / 2;	4✔
209	floatType y = center.point[1];	4✔
210
211	// whether the soma should be checked for the special case of 3 point soma
212	// for details see https://github.com/BlueBrain/MorphIO/issues/273
213	// If the 2nd and the 3rd point have the same x, z, d values then the only valid soma
214	// is: 1 1 x y z r -1 2 1 x (y-r) z r 1 3 1 x (y+r) z r 1
215	if ((child1.point[0] != x \|\| child2.point[0] != x \|\| child1.point[1] != y - r \|\|	8✔
216	child2.point[1] != y + r \|\| child1.point[2] != z \|\| child2.point[2] != z \|\|	4✔
217	child1.diameter != d \|\| child2.diameter != d) &&	4✔
218	std::fabs(child1.diameter - d) < morphio::epsilon &&	×
219	std::fabs(child2.diameter - d) < morphio::epsilon &&	×
220	std::fabs(child1.point[0] - x) < morphio::epsilon &&	×
221	std::fabs(child2.point[0] - x) < morphio::epsilon &&	×
222	std::fabs(child1.point[2] - z) < morphio::epsilon &&	8✔
223	std::fabs(child2.point[2] - z) < morphio::epsilon	×
224	) {
225	printError(Warning::SOMA_NON_CONFORM,	×
226	err_.WARNING_NEUROMORPHO_SOMA_NON_CONFORM(center, child1, child2));	×
227	}
228	}	4✔
229
230	void build_soma(const Samples& soma_samples) {	56✔
231	auto& soma = morph_.soma();	56✔
232
233	if (soma_samples.empty()) {	56✔
234	soma->type() = SOMA_UNDEFINED;	×
235	return;	×
236	} else if (soma_samples.size() == 1) {	56✔
237	Sample sample = soma_samples[0];	38✔
238
239	if (sample.parentId != SWC_ROOT && samples_.at(sample.parentId).type != SECTION_SOMA) {	38✔
240	throw morphio::SomaError(err_.ERROR_SOMA_WITH_NEURITE_PARENT(sample));	2✔
241	}
242
243	soma->type() = SOMA_SINGLE_POINT;	36✔
244	soma->points() = {sample.point};	36✔
245	soma->diameters() = {sample.diameter};	36✔
246	return;	36✔
247	} else if (soma_samples.size() == 3) {	18✔
248	const Sample& center = soma_samples[0];	4✔
249	const Sample& child1 = soma_samples[1];	4✔
250	const Sample& child2 = soma_samples[2];	4✔
251	// All soma that bifurcate with the first parent having two children are considered
252	// SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS
253	if(center.id == child1.parentId && center.id == child2.parentId){	4✔
254	soma->type() = SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS;	4✔
255	soma->points() = {center.point, child1.point, child2.point};	4✔
256	soma->diameters() = {center.diameter, child1.diameter, child2.diameter};	4✔
257	_checkNeuromorph3PointSoma(soma_samples);	4✔
258	return;	4✔
259	}
260	}
261	// might also have 3 points at this point, as well
262
263	// a "normal" SWC soma
264	soma->type() = SOMA_CYLINDERS;	14✔
265	auto& points = soma->points();	14✔
266	auto& diameters = soma->diameters();	14✔
267	points.reserve(soma_samples.size());	14✔
268	diameters.reserve(soma_samples.size());	14✔
269
270	size_t parent_count = 0;	14✔
271	for (const auto& s : soma_samples) {	56✔
272	if (s.parentId == SWC_ROOT) {	44✔
273	parent_count++;	16✔
274	} else if (samples_.count(s.parentId) == 0) {	28✔
275	throw morphio::MissingParentError(err_.ERROR_MISSING_PARENT(s));	×
276	} else if (samples_.at(s.parentId).type != SECTION_SOMA) {	28✔
277	throw morphio::SomaError(err_.ERROR_SOMA_WITH_NEURITE_PARENT(s));	×
278	}
279
280	if (children_.count(s.id) > 0 && children_.at(s.id).size() > 1) {	44✔
281	std::vector<Sample> soma_bifurcations;	20✔
282	for (auto id : children_.at(s.id)) {	30✔
283	if (samples_[id].type == SECTION_SOMA) {	20✔
284	soma_bifurcations.push_back(samples_[id]);	12✔
285	}
286	}
287	if (soma_bifurcations.size() > 1) {	10✔
288	throw morphio::SomaError(err_.ERROR_SOMA_BIFURCATION(s, soma_bifurcations));	2✔
289	}
290	}
291	points.push_back(s.point);	42✔
292	diameters.push_back(s.diameter);	42✔
293	}
294
295	if (parent_count > 1) {	12✔
296	throw morphio::SomaError(err_.ERROR_MULTIPLE_SOMATA(soma_samples));	2✔
297	}
298	}
299
300	void buildSWC(const Samples& samples) {	64✔
301	Samples soma_samples;	128✔
302	Samples root_samples;	128✔
303
304	for (const auto& sample: samples) {	562✔
305	// { checks
306	if (sample.diameter < morphio::epsilon) {	504✔
307	printError(Warning::ZERO_DIAMETER, err_.WARNING_ZERO_DIAMETER(sample));	10✔
308	}
309
310	if (sample.parentId == sample.id) {	504✔
311	throw morphio::RawDataError(err_.ERROR_SELF_PARENT(sample));	2✔
312	}
313
314	if (sample.type >= morphio::SECTION_OUT_OF_RANGE_START \|\| sample.type <= 0) {	502✔
315	throw morphio::RawDataError(
316	err_.ERROR_UNSUPPORTED_SECTION_TYPE(sample.lineNumber, sample.type));	2✔
317	}
318	if (sample.parentId == SWC_ROOT && sample.type != SECTION_SOMA) {	500✔
319	printError(Warning::DISCONNECTED_NEURITE,	4✔
320	err_.WARNING_DISCONNECTED_NEURITE(sample));	8✔
321	}
322	// } checks
323
324	if (sample.type == SECTION_SOMA) {	500✔
325	soma_samples.push_back(sample);	104✔
326	}
327
328	if (sample.parentId == SWC_ROOT \|\| sample.type == SECTION_SOMA) {	500✔
329	root_samples.push_back(sample);	108✔
330	}
331
332	if (!samples_.insert({sample.id, sample}).second) {	500✔
333	throw RawDataError(err_.ERROR_REPEATED_ID(samples[sample.id], sample));	2✔
334	}
335
336	children_[sample.parentId].push_back(sample.id);	498✔
337	}
338
339	// can only check for missing parents once all samples are loaded
340	// since it's possible there may be forward references
341	for (const auto& sample: samples) {	548✔
342	if(sample.parentId != SWC_ROOT && samples_.count(sample.parentId) == 0){	492✔
343	throw morphio::MissingParentError(err_.ERROR_MISSING_PARENT(sample));	2✔
344	}
345	}
346
347	build_soma(soma_samples);	56✔
348
349	std::unordered_map<DeclaredID, std::shared_ptr<morphio::mut::Section>> declared_to_swc;	100✔
350	declared_to_swc.reserve(samples.size());	50✔
351
352	for (const Sample& root_sample : root_samples) {	134✔
353	if (children_.count(root_sample.id) == 0) {	84✔
354	continue;	18✔
355	}
356
357	// https://neuromorpho.org/SomaFormat.html
358	// "The second and third soma points, as well as all starting points
359	// (roots) of dendritic and axonal arbors have this first point as
360	// the parent (parent ID 1)."
361	if (morph_.soma()->type() == SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS &&	66✔
362	root_sample.type == SECTION_SOMA && root_sample.id != 1) {	66✔
363	printError(Warning::WRONG_ROOT_POINT, err_.WARNING_WRONG_ROOT_POINT(root_sample));	×
364	}
365
366	for (unsigned int child_id : children_.at(root_sample.id)) {	218✔
367	if (samples_.at(child_id).type == SECTION_SOMA) {	152✔
368	continue;	34✔
369	}
370	if (root_sample.type == SECTION_SOMA) {	118✔
371	assembleSections(child_id,	236✔
372	DeclaredID(root_sample.id),	118✔
373	declared_to_swc,
374	morph_.soma()->points()[0],	118✔
375	morph_.soma()->diameters()[0],	118✔
376	true);
377	} else {
378	// this is neurite as the start
379	assembleSections(root_sample.id,	×
380	DeclaredID(SWC_ROOT),
381	declared_to_swc,
382	root_sample.point,	×
383	root_sample.diameter,	×
384	true);
385	break;	×
386	}
387	}
388	}
389	}	50✔
390
391	void assembleSections(	238✔
392	unsigned int id,
393	DeclaredID parent_id,
394	std::unordered_map<DeclaredID, std::shared_ptr<morphio::mut::Section>>& declared_to_swc,
395	const Point& start_point,
396	floatType start_diameter,
397	bool is_root) {
398
399	morphio::Property::PointLevel properties;	476✔
400	auto& points = properties._points;	238✔
401	auto& diameters = properties._diameters;	238✔
402
403	auto appendSection = [&](DeclaredID section_id_, DeclaredID parent_id_, SectionType starting_section_type) {	238✔
404	std::shared_ptr<morphio::mut::Section> new_section;	238✔
405	if (is_root) {	238✔
406	new_section = morph_.appendRootSection(properties, starting_section_type);	118✔
407	} else {
408	new_section = declared_to_swc.at(parent_id_)	120✔
409	->appendSection(properties, starting_section_type);	120✔
410	}
411	declared_to_swc[section_id_] = new_section;	238✔
412	};	238✔
413
414	auto get_child_count = [&](unsigned int child_id) {	388✔
415	return children_.count(child_id) == 0 ? 0 : children_.at(child_id).size();	388✔
416	};	238✔
417
418	const Sample* sample = &samples_.at(id);	238✔
419
420	// create duplicate point if needed
421	if (!is_root && sample->point != start_point /\|\| sample->diameter != start_diameter /) {	238✔
422	//if (!(is_root && sample->type == SECTION_SOMA) && sample->point != start_point /\|\| sample->diameter != start_diameter /) {
423	points.push_back(start_point);	112✔
424	diameters.push_back(start_diameter);	112✔
425	}
426
427	// try and combine as many single samples into a single section as possible
428	size_t children_count = get_child_count(id);	238✔
429	while (children_count == 1) {	388✔
430	sample = &samples_.at(id);	150✔
431	if(sample->type != samples_.at(children_.at(id)[0]).type){	150✔
432	break;	×
433	}
434	points.push_back(sample->point);	150✔
435	diameters.push_back(sample->diameter);	150✔
436	id = children_.at(id)[0];	150✔
437	children_count = get_child_count(id);	150✔
438	}
439	sample = &samples_.at(id);	238✔
440	points.push_back(sample->point);	238✔
441	diameters.push_back(sample->diameter);	238✔
442	appendSection(DeclaredID(id), parent_id, sample->type);	238✔
443
444	if (children_count == 0) {	238✔
445	// section was already appended above, nothing to do
446	} else if (children_count == 1) {	60✔
447	// section_type changed
448	size_t offset = properties._points.size() - 1;	×
449	const Point& new_start_point = properties._points[offset];	×
450	floatType new_start_diameter = properties._diameters[offset];	×
451	assembleSections(children_.at(id)[0], DeclaredID(id), declared_to_swc, new_start_point, new_start_diameter, false);	×
452	} else {
453	size_t offset = properties._points.size() - 1;	60✔
454	const Point& new_start_point = properties._points[offset];	60✔
455	floatType new_start_diameter = properties._diameters[offset];	60✔
456	for (unsigned int child_id : children_.at(id)) {	180✔
457	assembleSections(child_id,	120✔
458	DeclaredID(id),
459	declared_to_swc,
460	new_start_point,
461	new_start_diameter,
462	false);
463	}
464	}
465	}	238✔
466
467	std::unordered_map<unsigned int, std::vector<unsigned int>> children_;
468	std::unordered_map<unsigned int, Sample> samples_;
469	mut::Morphology morph_;
470	ErrorMessages err_;
471	};
472
473
474	} // namespace details
475
476	namespace morphio {
477	namespace readers {
478	namespace swc {
479	Property::Properties load(const std::string& path,	66✔
480	const std::string& contents,
481	unsigned int options) {
482	auto properties = details::SWCBuilder(path).buildProperties(contents, options);	82✔
483
484	properties._cellLevel._cellFamily = NEURON;	50✔
485	properties._cellLevel._version = {"swc", 1, 0};	50✔
486	return properties;	50✔
487	}
488
489	} // namespace swc
490	} // namespace readers
491	} // 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