6545228678

Committed 17 Oct 2023 08:31AM UTC coverage: 76.668% (+0.6%) from 76.104%

Build # 6545228678

Build Type

Pull #476

github

Committed by

mgeplf

Commit Message

fix missing data

Pull Request Pull Request #476: Efficient swc build

Run Details

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

1988 of 2593 relevant lines covered (76.67%)

914.88 hits per line

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

95.2

/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() {
      if (done()) {
          return;
      }
      std::size_t pos = contents_.find_first_not_of(" \t\r", pos_);
      if (pos == std::string::npos) {
          pos_ = contents_.size();
          return;
      }
      pos_ = pos;
  }

  void advance_to_number() {
      while (!done() && 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 {};

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

BlueBrain / MorphIO / 6545228678

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