From d19e958e2f7f4bd4bc3be76a2ab38b891428d14c Mon Sep 17 00:00:00 2001 From: Juan Hernando Vieites Date: Mon, 19 Sep 2016 14:39:50 +0200 Subject: [PATCH] Fix for #90: SWC parser creates multiple soma sections. This happened when neurites were connected to arbitrary soma samples. --- brion/plugin/morphologySWC.cpp | 78 ++++++++++++++++++++++------ doc/Changelog.md | 4 ++ tests/data/swc/first_order_sections_ring.swc | 16 ++++++ tests/morphology.cpp | 44 ++++++++++++---- 4 files changed, 116 insertions(+), 26 deletions(-) create mode 100644 tests/data/swc/first_order_sections_ring.swc diff --git a/brion/plugin/morphologySWC.cpp b/brion/plugin/morphologySWC.cpp index 14724db..204e121 100644 --- a/brion/plugin/morphologySWC.cpp +++ b/brion/plugin/morphologySWC.cpp @@ -348,8 +348,14 @@ void MorphologySWC::_buildSampleTree( RawSWCInfo& info ) if( !root ) { - Sample& parent = samples[sample.parent]; - if( !parent.valid ) + if( sample.parent == int( currentSample )) + { + LBTHROW( std::runtime_error( + "Reading swc morphology file: " + info.filename + + ", found a sample point to itself" )); + } + Sample* parent = &samples[sample.parent]; + if( !parent->valid ) { std::stringstream msg; msg << "Reading swc morphology file: " << info.filename @@ -358,17 +364,51 @@ void MorphologySWC::_buildSampleTree( RawSWCInfo& info ) LBTHROW( std::runtime_error( msg.str( ))); } - if( parent.nextID != -1 ) + if( parent->type == SWC_SECTION_SOMA ) { - // The parent was already connected. Linking this sample - // to its sibling. - sample.siblingID = parent.nextID; - // This also means that a sequence of samples is now split - // in three different sections (a parent and two children). - info.numSections += 2; + // When the parent is the soma we have to handle it differently + // as we don't want to split a soma ring where neurites + // connect to arbitrary soma points in multiple sections. + if( sample.type == SWC_SECTION_SOMA ) + { + if( parent->nextID != -1 ) + { + LBWARN << "Reading swc morphology file: " + << info.filename + << ", found bifurcation in soma section"; + sample.siblingID = parent->nextID; + } + // Linking the parent to this sample. + parent->nextID = int(currentSample); + } + else + { + info.roots.push_back( currentSample ); + // Sections whose parent is the soma need their parent + // section to be assigned at this point. + sample.parentSection = 0; + } + } + else + { + if( sample.type == SWC_SECTION_SOMA ) + { + LBTHROW( std::runtime_error( + "Reading swc morphology file: " + info.filename + + ", found soma sample with neurite parent" )); + } + if( parent->nextID != -1 ) + { + // The parent was already connected. Linking this sample + // to its sibling. + sample.siblingID = parent->nextID; + // This also means that a sequence of samples is now split + // in three different sections (a parent and two children). + info.numSections += 2; + } + // Linking the parent to this sample. + parent->nextID = int(currentSample); } - // Linking the parent to this sample - parent.nextID = int(currentSample); } else { @@ -456,10 +496,18 @@ void MorphologySWC::_buildStructure( RawSWCInfo& info ) // Pushing first point of the section using the parent sample // if necessary - const Sample* parent = - sample->parent == -1 ? 0 : &samples[sample->parent]; - if( parent && parent->type != SWC_SECTION_SOMA ) - _points->push_back( parent->point ); + if( sample->parent != SWC_UNDEFINED_PARENT ) + { + const Sample* parent = &samples[sample->parent]; + // If the parent sections is the soma, we connect this section + // to the soma only if the soma is described with more than + // one sample (that is, sections are not connected to point somas). + if( parent->type != SWC_SECTION_SOMA || + parent->nextID != -1 || parent->parent != -1 ) + { + _points->push_back( parent->point ); + } + } // Iterate while we stay on the same section and push points // to the point vector. diff --git a/doc/Changelog.md b/doc/Changelog.md index 7736961..fca506c 100644 --- a/doc/Changelog.md +++ b/doc/Changelog.md @@ -3,6 +3,10 @@ Changelog {#Changelog} # Release 1.9.0 (git master) +* [94](https://github.com/BlueBrain/Brion/pull/94): + Fixed SWC morphology parser for morphologies with soma contour. The parser was + creating invalid soma sections when the first order sections where connected + to arbitrary soma sample points. * [89](https://github.com/BlueBrain/Brion/pull/89): Python wrapping of brain classes. * [88](https://github.com/BlueBrain/Brion/pull/88): diff --git a/tests/data/swc/first_order_sections_ring.swc b/tests/data/swc/first_order_sections_ring.swc new file mode 100644 index 0000000..0c55d34 --- /dev/null +++ b/tests/data/swc/first_order_sections_ring.swc @@ -0,0 +1,16 @@ +# Soma with neurites branching out of the middle +# This is the layout of the "contour soma". +1 1 0 0 1 0.0 -1 +2 1 0 0 2 0.0 1 +3 1 0 0 3 0.0 2 +4 1 0 0 4 0.0 3 +5 1 0 0 5 0.0 4 +# first neurite branch, from point 2 +6 2 1 1 6 0.5 2 +7 2 1 2 7 0.5 6 +# second neurite branch, from point 4 +8 3 2 1 8 0.5 3 +9 3 2 2 9 0.5 8 +# third neurite branch, from point 6 +10 4 3 1 10 0.5 4 +11 4 3 2 11 0.5 10 diff --git a/tests/morphology.cpp b/tests/morphology.cpp index e2a9f74..4cc3322 100644 --- a/tests/morphology.cpp +++ b/tests/morphology.cpp @@ -557,16 +557,38 @@ BOOST_AUTO_TEST_CASE( swc_first_order_sections ) const brion::Morphology source( path.string( )); checkEqualArrays( *source.readSections( stage ), 4, - V2i( 0, -1 ), V2i( 1, 0 ), V2i( 2, 0 ), V2i( 3, 0 )); - // The tree construction algorithm reserves the order of two sections + V2i( 0, -1 ), V2i( 1, 0 ), V2i( 2, 0 ), V2i( 3, 0 )); + // The tree construction algorithm reverses the order of the sections // compared to how they appear in the file checkEqualArrays( *source.readPoints( stage ), 4, V4f( 0, 0, 0, 20 ), - V4f( 0, 0, 1, 4 ), V4f( 0, 0, 3, 4 ), V4f( 0, 0, 2, 4 )); + V4f( 0, 0, 3, 4 ), V4f( 0, 0, 2, 4 ), V4f( 0, 0, 1, 4 )); checkEqualArrays( *source.readSectionTypes(), 4, - SOMA, AXON, APICAL_DENDRITE, DENDRITE ); + SOMA, APICAL_DENDRITE, DENDRITE, AXON ); } +BOOST_AUTO_TEST_CASE( swc_first_order_sections_from_arbitrary_points ) +{ + boost::filesystem::path path( BRION_TESTDATA ); + const brion::MorphologyRepairStage stage = brion::MORPHOLOGY_REPAIRED; + path /= "swc/first_order_sections_ring.swc"; + + const brion::Morphology source( path.string( )); + + checkEqualArrays( *source.readSections( stage ), 4, + V2i( 0, -1 ), V2i( 5, 0 ), V2i( 8, 0 ), V2i( 11, 0 )); + // The tree construction algorithm reverses the order of the sections + // compared to how they appear in the file + checkEqualArrays( *source.readPoints( stage ), 14, V4f( 0, 0, 1, 0 ), + V4f( 0, 0, 2, 0 ), V4f( 0, 0, 3, 0 ), V4f( 0, 0, 4, 0 ), + V4f( 0, 0, 5, 0 ), V4f( 0, 0, 4, 0 ), V4f( 3, 1, 10, 1 ), + V4f( 3, 2, 11, 1 ), V4f( 0, 0, 3, 0 ), V4f( 2, 1, 8, 1 ), + V4f( 2, 2, 9, 1 ), V4f( 0, 0, 2, 0 ), V4f( 1, 1, 6, 1 ), + V4f( 1, 2, 7, 1 )); + checkEqualArrays( *source.readSectionTypes(), 4, + SOMA, APICAL_DENDRITE, DENDRITE, AXON ); +} + BOOST_AUTO_TEST_CASE( swc_bifurcation ) { boost::filesystem::path path( BRION_TESTDATA ); @@ -594,11 +616,11 @@ BOOST_AUTO_TEST_CASE( swc_end_points ) const brion::Morphology source( path.string( )); checkEqualArrays( *source.readSections( stage ), 6, - V2i( 0, -1 ), V2i( 1, 0 ), V2i( 2, 1 ), V2i( 4, 1 ), - V2i( 7, 0 ), V2i( 8, 0 )); + V2i( 0, -1 ), V2i( 1, 0 ), V2i( 2, 0 ), V2i( 3, 0 ), + V2i( 4, 3 ), V2i( 6, 3 )); checkEqualArrays( *source.readSectionTypes(), 6, - SOMA, AXON, AXON, AXON, UNDEFINED, UNDEFINED ); + SOMA, UNDEFINED, UNDEFINED, AXON, AXON, AXON ); } BOOST_AUTO_TEST_CASE( swc_fork_points ) @@ -610,11 +632,11 @@ BOOST_AUTO_TEST_CASE( swc_fork_points ) const brion::Morphology source( path.string( )); checkEqualArrays( *source.readSections( stage ), 6, - V2i( 0, -1 ), V2i( 1, 0 ), V2i( 2, 1 ), V2i( 4, 1 ), - V2i( 7, 0 ), V2i( 8, 0 )); + V2i( 0, -1 ), V2i( 1, 0 ), V2i( 2, 0 ), V2i( 3, 0 ), + V2i( 4, 3 ), V2i( 6, 3 )); checkEqualArrays( *source.readSectionTypes(), 6, - SOMA, AXON, AXON, AXON, UNDEFINED, UNDEFINED ); + SOMA, UNDEFINED, UNDEFINED, AXON, AXON, AXON ); } BOOST_AUTO_TEST_CASE( swc_neuron ) @@ -624,7 +646,7 @@ BOOST_AUTO_TEST_CASE( swc_neuron ) brion::Morphology neuron( path.string( )); BOOST_CHECK_EQUAL( neuron.readPoints( brion::MORPHOLOGY_REPAIRED )->size(), - 927 ); + 933 ); BOOST_CHECK_EQUAL( neuron.getCellFamily(), brion::FAMILY_NEURON ); BOOST_CHECK( neuron.readPerimeters()->empty( )); }