5454321273

pending completion

Build # 5454321273

Build Type

push

github

Committed by

mgeplf

Commit Message

make readSamples standalone

Run Details

25 of 25 new or added lines in 1 file covered. (100.0%)

1982 of 2581 relevant lines covered (76.79%)

915.12 hits per line

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

93.36

/src/readers/morphologySWC.cpp

#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 {
// 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", 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;
}

}  // unnamed namespace

namespace morphio {
namespace readers {
namespace swc {
/**
  Parsing SWC according to this specification:
http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
 **/
class SWCBuilder
{
  public:
    explicit SWCBuilder(const std::string& path)
        : err_(path) {}

    void raiseIfBrokenSoma(const Sample& sample) {
        if (sample.type != SECTION_SOMA) {
            return;
        }

        if (sample.parentId != SWC_ROOT && !children_[sample.id].empty()) {
            std::vector<Sample> soma_bifurcations;
            for (unsigned int id : children_[sample.id]) {
                if (samples_[id].type == SECTION_SOMA) {
                    soma_bifurcations.push_back(samples_[id]);
                } else {
                    neurite_wrong_root_.push_back(samples_[id]);
                }
            }

            if (soma_bifurcations.size() > 1) {
                throw morphio::SomaError(err_.ERROR_SOMA_BIFURCATION(sample, soma_bifurcations));
            }
        }

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

    void checkSoma() {
        std::vector<Sample> somata;
        for (auto id : children_[SWC_ROOT]) {
            if (samples_.at(id).type == SECTION_SOMA) {
                somata.push_back(samples_[id]);
            }
        }

        if (somata.size() > 1) {
            throw morphio::SomaError(err_.ERROR_MULTIPLE_SOMATA(somata));
        }

        if (somata.empty()) {
            printError(Warning::NO_SOMA_FOUND, err_.WARNING_NO_SOMA_FOUND());
        } else {
            for (const auto& sample_pair : samples_) {
                const auto& sample = sample_pair.second;
                if (sample.parentId == SWC_ROOT && sample.type != SECTION_SOMA) {
                    printError(Warning::DISCONNECTED_NEURITE, err_.WARNING_DISCONNECTED_NEURITE(sample));
                }
            }
        }
    }

    bool isRootPoint(const Sample& sample) {
        bool isOrphanNeurite = sample.parentId == SWC_ROOT && sample.type != SECTION_SOMA;
        return isOrphanNeurite ||
               (sample.type != SECTION_SOMA &&
                samples_.at(sample.parentId).type == SECTION_SOMA);  // Exclude soma bifurcations
    }

    bool isSectionStart(const Sample& sample) {
        return (isRootPoint(sample) ||
                (sample.parentId != SWC_ROOT &&
                 isSectionEnd(samples_.at(sample.parentId))));  // Standard section
    }

    bool isSectionEnd(const Sample& sample) {
        return sample.id == lastSomaPoint_ ||        // End of soma
               children_[sample.id].empty() ||       // Reached leaf
               (children_[sample.id].size() >= 2 &&  // Reached neurite bifurcation
                sample.type != SECTION_SOMA);
    }

    template <typename T>
    void appendSample(const std::shared_ptr<T>& somaOrSection, const Sample& sample) {
        somaOrSection->points().push_back(sample.point);
        somaOrSection->diameters().push_back(sample.diameter);
    }

    std::vector<unsigned int> constructDepthFirstSamples() {
        std::vector<unsigned int> ret;
        ret.reserve(samples_.size());
        const auto pushChildren = [&](const auto& f, unsigned int id) -> void {
            for (unsigned int childId : children_[id]) {
                ret.push_back(childId);
                f(f, childId);
            }
        };

        pushChildren(pushChildren, SWC_ROOT);

        return ret;
    }

    void _checkNeuroMorphoSoma(const Sample& root, const std::vector<Sample>& _children) {
        floatType x = root.point[0];
        floatType y = root.point[1];
        floatType z = root.point[2];
        floatType d = root.diameter;
        floatType r = root.diameter / 2;
        const Sample& child1 = _children[0];
        const Sample& child2 = _children[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
        bool isSuited = 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;
        if (!isSuited) {
            return;
        }
        // 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) {
            printError(Warning::SOMA_NON_CONFORM,
                       err_.WARNING_NEUROMORPHO_SOMA_NON_CONFORM(root, child1, child2));
        }
    }

    SomaType somaType() {
        switch (morph_.soma()->points().size()) {
        case 0: {
            return SOMA_UNDEFINED;
        }
        case 1: {
            return SOMA_SINGLE_POINT;
        }
        case 2: {
            return SOMA_CYLINDERS;
        }
            // NeuroMorpho format is characterized by a 3 points soma
            // with a bifurcation at soma root
        case 3: {
            uint32_t somaRootId = children_[SWC_ROOT][0];
            const auto& somaChildren = children_[somaRootId];

            std::vector<Sample> children_soma_points;
            for (unsigned int child : somaChildren) {
                if (samples_.at(child).type == SECTION_SOMA) {
                    children_soma_points.push_back(samples_.at(child));
                }
            }

            if (children_soma_points.size() == 2) {
                //  NeuroMorpho is the main provider of morphologies, but they
                //  with SWC as their default file format: they convert all
                //  uploads to SWC.  In the process of conversion, they turn all
                //  somas into their custom 'Three-point soma representation':
                //   http://neuromorpho.org/SomaFormat.html

                if (!ErrorMessages::isIgnored(Warning::SOMA_NON_CONFORM)) {
                    _checkNeuroMorphoSoma(samples_.at(somaRootId), children_soma_points);
                }

                return SOMA_NEUROMORPHO_THREE_POINT_CYLINDERS;
            }
            return SOMA_CYLINDERS;
        }
        default:
            return SOMA_CYLINDERS;
        }
    }

    Property::Properties buildProperties(const std::string& contents, unsigned int options) {
        const auto samples = readSamples(contents, err_);

        for (const auto& sample: samples) {
            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 (!samples_.insert({sample.id, sample}).second) {
                throw RawDataError(err_.ERROR_REPEATED_ID(samples[sample.id], sample));
            }

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

            if (sample.type == SECTION_SOMA) {
                lastSomaPoint_ = sample.id;
            }
        }

        for (const auto& sample: samples) {
            raiseIfBrokenSoma(sample);

            if (sample.parentId != SWC_ROOT && samples_.count(sample.parentId) == 0) {
                throw morphio::MissingParentError(err_.ERROR_MISSING_PARENT(sample));
            }
        }

        checkSoma();

        // The process might occasionally creates empty section before
        // filling them so the warning is ignored
        bool originalIsIgnored = err_.isIgnored(morphio::Warning::APPENDING_EMPTY_SECTION);
        set_ignored_warning(morphio::Warning::APPENDING_EMPTY_SECTION, true);

        std::vector<unsigned int> depthFirstSamples = constructDepthFirstSamples();
        for (unsigned int id : depthFirstSamples) {
            const Sample& sample = samples_.at(id);

            // Bifurcation right at the start
            if (isRootPoint(sample) && isSectionEnd(sample)) {
                continue;
            }

            if (isSectionStart(sample)) {
                _processSectionStart(sample);
            } else if (sample.type != SECTION_SOMA) {
                swcIdToSectionId_[sample.id] = swcIdToSectionId_[sample.parentId];
            }

            if (sample.type == SECTION_SOMA) {
                appendSample(morph_.soma(), sample);
            } else {
                appendSample(morph_.section(swcIdToSectionId_.at(sample.id)), sample);
            }
        }

        if (morph_.soma()->points().size() == 3 && !neurite_wrong_root_.empty()) {
            printError(morphio::WRONG_ROOT_POINT, err_.WARNING_WRONG_ROOT_POINT(neurite_wrong_root_));
        }

        morph_.applyModifiers(options);

        Property::Properties properties = morph_.buildReadOnly();
        properties._cellLevel._somaType = somaType();

        set_ignored_warning(morphio::Warning::APPENDING_EMPTY_SECTION, originalIsIgnored);

        return properties;
    }

    /**
      - Append last point of previous section if current section is not a root section
      - Update the parent ID of the new section
     **/
    void _processSectionStart(const Sample& sample) {
        Property::PointLevel properties;

        uint32_t id = 0;

        if (isRootPoint(sample)) {
            id = morph_.appendRootSection(properties, sample.type)->id();
        } else {
            // Duplicating last point of previous section if there is not already a duplicate
            Sample& parent_sample = samples_.at(sample.parentId);
            if (sample.point != parent_sample.point) {
                properties._points.push_back(parent_sample.point);
                properties._diameters.push_back(parent_sample.diameter);
            }

            // Handle the case, bifurcatation at root point
            if (isRootPoint(parent_sample)) {
                id = morph_.appendRootSection(properties, sample.type)->id();
            } else {
                id = morph_.section(swcIdToSectionId_[sample.parentId])
                    ->appendSection(properties, sample.type)
                    ->id();
            }
        }

        swcIdToSectionId_[sample.id] = id;
    }

  private:
    // SWC Id of the last point of a section to morphio::mut::Section ID
    std::unordered_map<uint32_t, uint32_t> swcIdToSectionId_;

    // Neurite that do not have parent ID = 1, allowed for soma contour, not
    // 3-pts soma
    std::vector<Sample> neurite_wrong_root_;

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

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

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

}  // namespace swc
}  // namespace readers
}  // namespace morphio

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

BlueBrain / MorphIO / 5454321273

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