diff --git a/brain/CMakeLists.txt b/brain/CMakeLists.txt index 8677add..848a701 100644 --- a/brain/CMakeLists.txt +++ b/brain/CMakeLists.txt @@ -8,16 +8,22 @@ source_group(\\ FILES CMakeLists.txt) set(BRAIN_PUBLIC_HEADERS circuit.h - morphology.h types.h + neuron/morphology.h + neuron/section.h + neuron/soma.h ) set(BRAIN_HEADERS + neuron/morphologyImpl.h ) set(BRAIN_SOURCES circuit.cpp - morphology.cpp + neuron/morphology.cpp + neuron/morphologyImpl.cpp + neuron/section.cpp + neuron/soma.cpp ) set(BRAIN_PUBLIC_INCLUDE_DIRECTORIES ${Boost_INCLUDE_DIRS}) diff --git a/brain/circuit.cpp b/brain/circuit.cpp index 81a8745..fa75ab0 100644 --- a/brain/circuit.cpp +++ b/brain/circuit.cpp @@ -18,7 +18,7 @@ */ #include "circuit.h" -#include "morphology.h" +#include "neuron/morphology.h" #include #include @@ -102,11 +102,11 @@ URIs Circuit::getMorphologyURIs( const GIDSet& gids ) const return uris; } -Morphologies Circuit::loadMorphologies( const GIDSet& gids, +neuron::Morphologies Circuit::loadMorphologies( const GIDSet& gids, const Coordinates coords ) const { const URIs& uris = getMorphologyURIs( gids ); - Morphologies result; + neuron::Morphologies result; result.reserve( uris.size( )); if( coords == COORDINATES_GLOBAL ) @@ -116,22 +116,22 @@ Morphologies Circuit::loadMorphologies( const GIDSet& gids, { const URI& uri = uris[i]; const brion::Morphology raw( uri.getPath( )); - result.push_back( MorphologyPtr( new Morphology( raw, - transforms[i] ))); + result.push_back( neuron::MorphologyPtr( + new neuron::Morphology( raw, transforms[i] ))); } return result; } - std::map< std::string, MorphologyPtr > loaded; + std::map< std::string, neuron::MorphologyPtr > loaded; for( size_t i = 0; i < uris.size(); ++i ) { const URI& uri = uris[i]; - MorphologyPtr& morphology = loaded[uri.getPath()]; + neuron::MorphologyPtr& morphology = loaded[uri.getPath()]; if( !morphology ) { const brion::Morphology raw( uri.getPath( )); - morphology.reset( new Morphology( raw )); + morphology.reset( new neuron::Morphology( raw )); } result.push_back( morphology ); diff --git a/brain/circuit.h b/brain/circuit.h index 9aef83b..97f682f 100644 --- a/brain/circuit.h +++ b/brain/circuit.h @@ -75,8 +75,8 @@ class Circuit : public boost::noncopyable * will shared the same Morphology object in the list. If global * coordinates are requested, all Morphology objects are unique. */ - BRAIN_API Morphologies loadMorphologies( const GIDSet& gids, - Coordinates coords ) const; + BRAIN_API neuron::Morphologies loadMorphologies( const GIDSet& gids, + Coordinates coords ) const; /** @return The positions of the given cells. */ BRAIN_API Vector3fs getPositions( const GIDSet& gids ) const; diff --git a/brain/morphology.cpp b/brain/morphology.cpp deleted file mode 100644 index 4feb440..0000000 --- a/brain/morphology.cpp +++ /dev/null @@ -1,215 +0,0 @@ - -/* Copyright (c) 2013-2015, EPFL/Blue Brain Project - * Juan Hernando - * - * This file is part of Brion - * - * This library is free software; you can redistribute it and/or modify it under - * the terms of the GNU Lesser General Public License version 3.0 as published - * by the Free Software Foundation. - * - * This library is distributed in the hope that it will be useful, but WITHOUT - * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with this library; if not, write to the Free Software Foundation, Inc., - * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - */ - -#include "morphology.h" - -#include -#include - -#include - -#include - -#include - -namespace brain -{ - -class Morphology::Impl -{ -public: - const brion::Vector4fsPtr points; - const brion::Vector2isPtr sections; - const brion::SectionTypesPtr types; - const brion::Vector2isPtr apicals; - Vector3f somaPosition; - - Impl( const brion::Morphology& morphology ) - : points( morphology.readPoints( brion::MORPHOLOGY_UNDEFINED )) - , sections( morphology.readSections( brion::MORPHOLOGY_UNDEFINED )) - , types( morphology.readSectionTypes( )) - , apicals( morphology.readApicals( )) - , somaPosition( Vector3f::ZERO ) // the soma is assumed to be - // centered at the 0, 0, 0 - { - } - - uint32_ts getSectionIDs( const SectionTypes& requestedTypes ) const - { - std::bitset< brion::SECTION_APICAL_DENDRITE > bits; - BOOST_FOREACH( const SectionType type, requestedTypes ) - bits[size_t( type )] = true; - - uint32_ts result; - for( size_t i = 0; i != types->size(); ++i ) - { - const SectionType type = ( *types )[i]; - if( type == brion::SECTION_ALL ) - LBWARN << "Unknown section type " << int(type) << std::endl; - else - if( bits[size_t( type )] ) - result.push_back( i ); - } - return result; - } - - Vector4fs getSectionSamples( - const size_t sectionID, const floats& samplePoints ) const - { - // If the section is the soma return directly the soma position. - if(( *types )[sectionID] == brion::SECTION_SOMA ) - { - Vector4fs result; - for( size_t i = 0; i != samplePoints.size(); ++i ) - result.push_back( somaPosition ); - return result; - } - - const size_t start = ( *sections )[sectionID][0]; - const size_t end = sectionID == sections->size() - 1 ? - points->size() : ( *sections )[sectionID + 1][0]; - - Vector4fs result; - - if( end <= start ) - { - LBWARN << "Trying to sample broken morphology or empty section " - << sectionID << std::endl; - return result; - } - - result.reserve( samplePoints.size( )); - - // Dealing with the degenerate case of single point sections. - if( start + 1 == end ) - { - for( size_t i = 0; i != samplePoints.size(); ++i) - result.push_back(( *points )[start] ); - return result; - } - - // Computing the accumulated length at each segment of the section. - floats accumLengths; - accumLengths.reserve( end - start ); - accumLengths.push_back(0); - for( size_t i = start; i != end - 1; ++i ) - { - const Vector4f& segmentStart = ( *points )[i]; - const Vector4f& segmentEnd = ( *points )[i + 1]; - accumLengths.push_back( accumLengths.back() + - ( segmentEnd - segmentStart ).length( )); - } - const float totalLength = accumLengths.back(); - - BOOST_FOREACH( const float point, samplePoints ) - { - // Finding the segment index for the requested sampling position. - const float length = - std::max( 0.f, std::min( 1.f, point )) * totalLength; - size_t index = 0; - for ( ; accumLengths[index + 1] < length && - index < accumLengths.size() - 1; ++index ) - ; - - // Interpolating the cross section at point. - const float alpha = ( length - accumLengths[index] ) / - ( accumLengths[index + 1] - accumLengths[index] ); - const Vector4f sample = ( *points )[start + index + 1] * alpha + - ( *points )[start + index] * (1 - alpha ); - result.push_back( sample ); - } - - return result; - } - - void transform( const Matrix4f& transformation ) - { - #pragma omp parallel for - for( size_t i = 0; i < points->size(); ++i) - { - Vector4f& p = ( *points )[i]; - const Vector3f pp = transformation * p.get_sub_vector< 3 >(); - p.get_sub_vector< 3 >() = pp; - } - somaPosition = transformation * somaPosition; - } -}; - -Morphology::Morphology( const URI& source, const Matrix4f& transform ) - : _impl( new Impl( brion::Morphology( source.getPath( )))) -{ - _impl->transform( transform ); -} - -Morphology::Morphology( const brion::Morphology& morphology, - const Matrix4f& transform ) - : _impl( new Impl( morphology )) -{ - _impl->transform( transform ); -} - -Morphology::Morphology( const URI& source ) - : _impl( new Impl( brion::Morphology( source.getPath( )))) -{ -} - -Morphology::Morphology( const brion::Morphology& morphology ) - : _impl( new Impl( morphology )) -{ -} - -Morphology::~Morphology() -{ - delete _impl; -} - -const Vector4fs& Morphology::getPoints() const -{ - return *_impl->points; -} - -const Vector2is& Morphology::getSections() const -{ - return *_impl->sections; -} - -const SectionTypes& Morphology::getSectionTypes() const -{ - return *_impl->types; -} - -const Vector2is& Morphology::getApicals() const -{ - return *_impl->apicals; -} - -uint32_ts Morphology::getSectionIDs( const SectionTypes& types ) const -{ - return _impl->getSectionIDs( types ); -} - -Vector4fs Morphology::getSectionSamples( - const size_t sectionID, const floats& points ) const -{ - return _impl->getSectionSamples( sectionID, points ); -} - -} - diff --git a/brain/neuron/morphology.cpp b/brain/neuron/morphology.cpp new file mode 100644 index 0000000..24839e5 --- /dev/null +++ b/brain/neuron/morphology.cpp @@ -0,0 +1,142 @@ +/* Copyright (c) 2013-2015, EPFL/Blue Brain Project + * Juan Hernando + * + * This file is part of Brion + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License version 3.0 as published + * by the Free Software Foundation. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + */ + +#include "morphology.h" +#include "section.h" +#include "soma.h" + +#include "morphologyImpl.h" + +#include + +#include + +#include + +namespace brain +{ +namespace neuron +{ + +Morphology::Morphology( const URI& source, const Matrix4f& transform ) + : _impl( new Impl( brion::Morphology( source.getPath( )))) +{ + _impl->ref(); + _impl->transform( transform ); + _impl->transformation = transform; +} + +Morphology::Morphology( const brion::Morphology& morphology, + const Matrix4f& transform ) + : _impl( new Impl( morphology )) +{ + _impl->ref(); + _impl->transform( transform ); + _impl->transformation = transform; +} + +Morphology::Morphology( const URI& source ) + : _impl( new Impl( brion::Morphology( source.getPath( )))) +{ + _impl->ref(); +} + +Morphology::Morphology( const brion::Morphology& morphology ) + : _impl( new Impl( morphology )) +{ + _impl->ref(); +} + +Morphology::~Morphology() +{ + _impl->unref(); +} + +const Vector4fs& Morphology::getPoints() const +{ + return *_impl->points; +} + +const Vector2is& Morphology::getSections() const +{ + return *_impl->sections; +} + +const SectionTypes& Morphology::getSectionTypes() const +{ + return *_impl->types; +} + +const Vector2is& Morphology::getApicals() const +{ + return *_impl->apicals; +} + +uint32_ts Morphology::getSectionIDs( const SectionTypes& types ) const +{ + return _impl->getSectionIDs( types, false ); +} + +Sections Morphology::getSections( const SectionType type ) const +{ + const SectionTypes types( 1, type ); + const uint32_ts ids = _impl->getSectionIDs( types, true ); + Sections result; + BOOST_FOREACH( const uint32_t id, ids ) + result.push_back( Section( id, _impl )); + return result; +} + +Sections Morphology::getSections( const SectionTypes& types ) const +{ + const uint32_ts ids = _impl->getSectionIDs( types, true ); + Sections result; + BOOST_FOREACH( const uint32_t id, ids ) + result.push_back( Section( id, _impl )); + return result; +} + +Section Morphology::getSection( const uint32_t& id ) const +{ + if(( *_impl->types )[id] == SECTION_SOMA ) + LBTHROW( + std::runtime_error( "The soma cannot be accessed as a Section" )); + + if( _impl->sections->size() <= id ) + { + std::stringstream msg; + msg << "Section ID out of range: " << id; + LBTHROW( std::runtime_error( msg.str( ))); + } + + return Section( id, _impl ); +} + +Soma Morphology::getSoma() const +{ + return Soma( _impl ); +} + +const Matrix4f& Morphology::getTransformation() const +{ + return _impl->transformation; +} + +} +} diff --git a/brain/morphology.h b/brain/neuron/morphology.h similarity index 70% rename from brain/morphology.h rename to brain/neuron/morphology.h index 66c10d8..5478a32 100644 --- a/brain/morphology.h +++ b/brain/neuron/morphology.h @@ -1,4 +1,3 @@ - /* Copyright (c) 2013-2015, EPFL/Blue Brain Project * Juan Hernando * @@ -18,8 +17,8 @@ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */ -#ifndef BRAIN_MORPHOLOGY -#define BRAIN_MORPHOLOGY +#ifndef BRAIN_NEURON_MORPHOLOGY +#define BRAIN_NEURON_MORPHOLOGY #include #include @@ -28,6 +27,8 @@ namespace brain { +namespace neuron +{ /** * Wrapper around brion::Morphology with higher level functions. @@ -46,12 +47,16 @@ namespace brain class Morphology : public boost::noncopyable { public: + /** @internal */ + class Impl; + /** * Create a morphology from a URI, load all the data and transform * the points. * @param source URI of the morphology data source. * @param transform the transformation matrix to apply to the points. * Radii will not be affected by this transformation. + * @throw runtime_error if an inconsistency is detected in the input file. */ BRAIN_API Morphology( const URI& source, const Matrix4f& transform ); @@ -61,6 +66,7 @@ class Morphology : public boost::noncopyable * @param morphology the brion::Morphology to load from. * @param transform the transformation matrix to apply to the points. * Radii will not be affected by this transformation. + * @throw runtime_error if an inconsistency is detected in the input file. */ BRAIN_API Morphology( const brion::Morphology& morphology, const Matrix4f& transform ); @@ -68,12 +74,14 @@ class Morphology : public boost::noncopyable /** * Create a morphology from a URI and load all the data. * @param source URI of the morphology data source. + * @throw runtime_error if an inconsistency is detected in the input file. */ BRAIN_API Morphology( const URI& source ); /** * Create a morphology from a brion::Morphology and load all the data. * @param morphology the brion::Morphology to load from. + * @throw runtime_error if an inconsistency is detected in the input file. */ BRAIN_API Morphology( const brion::Morphology& morphology ); @@ -91,31 +99,39 @@ class Morphology : public boost::noncopyable /** @copydoc brion::Morphology::getApicals */ BRAIN_API const Vector2is& getApicals() const; + /** Return the list of ids for the given section types. */ + BRAIN_API uint32_ts getSectionIDs( const SectionTypes& types ) const; + /** - * Return the list of ids for the given section types. + * Return the sections which have the given section type. + * If type is SECTION_SOMA an empty list is returned. */ - BRAIN_API uint32_ts getSectionIDs( const SectionTypes& types ) const; + BRAIN_API Sections getSections( SectionType type ) const; /** - * Return a list of points sampling a sections at discrete locations or all - * the available points of no sample locations are given. - * - * @param sectionID Section to sample. If the section is the soma this - * function will return the soma position for all sampling positions. - * The soma position is assumed to be (0, 0, 0) unless the morphology - * has been transformed. - * @param points Normalized positions of the sample points along the - * section. Values will be clampled to [0, 1] before sampling. - * @return The vector with all the points of the section if points is empty. - * The section sampled at the given points otherwise. + * Return the sections which have any of the given section types. + * No sections are returned for the type SECTION_SOMA. */ - BRAIN_API Vector4fs getSectionSamples( - size_t sectionID, const floats& points = brion::floats( )) const; + BRAIN_API Sections getSections( const SectionTypes& types ) const; + + /** + * Return the Section with the given id. + * @throw runtime_error if the id is out of range or the given id refers to + * a soma section. + */ + BRAIN_API Section getSection( const uint32_t& id ) const; + + /** Return the object with the information about the neuron soma */ + BRAIN_API Soma getSoma() const; + + /** Return the transformation that was passed to the constructor or the + identity matrix is no transformation was given. */ + BRAIN_API const Matrix4f& getTransformation() const; private: - class Impl; Impl* const _impl; }; } +} #endif diff --git a/brain/neuron/morphologyImpl.cpp b/brain/neuron/morphologyImpl.cpp new file mode 100644 index 0000000..98b980f --- /dev/null +++ b/brain/neuron/morphologyImpl.cpp @@ -0,0 +1,253 @@ + +/* Copyright (c) 2013-2015, EPFL/Blue Brain Project + * Juan Hernando + * + * This file is part of Brion + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License version 3.0 as published + * by the Free Software Foundation. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + */ + +#include "morphologyImpl.h" +#include "section.h" + +#include + +#include + +#include + +#include + +namespace brain +{ +namespace neuron +{ + +Morphology::Impl::Impl( const brion::Morphology& morphology ) + : points( morphology.readPoints( MORPHOLOGY_UNDEFINED )) + , sections( morphology.readSections( MORPHOLOGY_UNDEFINED )) + , types( morphology.readSectionTypes( )) + , apicals( morphology.readApicals( )) + , transformation( Matrix4f::IDENTITY ) +{ + _extractChildrenLists(); + + const uint32_ts ids = + getSectionIDs( SectionTypes( 1, SECTION_SOMA ), false ); + + if( ids.size() != 1 ) + LBTHROW( std::runtime_error( + "Bad input morphology. None or more than one soma found" )); + somaSection = ids[0]; +} + +SectionRange Morphology::Impl::getSectionRange( const uint32_t sectionID ) const +{ + const size_t start = ( *sections )[sectionID][0]; + const size_t end = sectionID == sections->size() - 1 ? + points->size() : ( *sections )[sectionID + 1][0]; + return std::make_pair( start, end ); +} + +uint32_ts Morphology::Impl::getSectionIDs( + const SectionTypes& requestedTypes, const bool excludeSoma ) const +{ + std::bitset< SECTION_APICAL_DENDRITE > bits; + BOOST_FOREACH( const SectionType type, requestedTypes ) + { + if( type != SECTION_SOMA || !excludeSoma ) + bits[size_t( type )] = true; + } + + uint32_ts result; + for( size_t i = 0; i != types->size(); ++i ) + { + const SectionType type = ( *types )[i]; + if( type == SECTION_ALL ) + LBWARN << "Unknown section type " << int(type) << std::endl; + else + if( bits[size_t( type )] ) + result.push_back( i ); + } + return result; +} + +float Morphology::Impl::getSectionLength( const uint32_t sectionID ) const +{ + if( _sectionLengths.size() <= sectionID ) + _sectionLengths.resize( sectionID + 1 ); + + float& length = _sectionLengths[sectionID]; + + if( length == 0 && ( *types )[sectionID] != SECTION_SOMA ) + length = _computeSectionLength( sectionID ); + return length; +} + +Vector4fs Morphology::Impl::getSectionSamples( const uint32_t sectionID ) const +{ + const SectionRange range = getSectionRange( sectionID ); + Vector4fs result; + result.reserve( range.second - range.first ); + result.insert( result.end(), points->begin() + range.first, + points->begin() + range.second); + return result; +} + +Vector4fs Morphology::Impl::getSectionSamples( const uint32_t sectionID, + const floats& samplePoints ) const +{ + const SectionRange range = getSectionRange( sectionID ); + + // If the section is the soma return directly the soma position. + if(( *types )[sectionID] == SECTION_SOMA ) + // This code shouldn't be reached. + LBTHROW( std::runtime_error( "Invalid method called on soma section" )); + + // Dealing with the degenerate case of single point sections. + if( range.first + 1 == range.second ) + return Vector4fs(samplePoints.size( ), ( *points )[range.first] ); + + Vector4fs result; + result.reserve( samplePoints.size( )); + + const floats accumLengths = _computeAccumulatedLengths( range ); + const float totalLength = accumLengths.back(); + + BOOST_FOREACH( const float point, samplePoints ) + { + // Finding the segment index for the requested sampling position. + const float length = + std::max( 0.f, std::min( 1.f, point )) * totalLength; + size_t index = 0; + for ( ; accumLengths[index + 1] < length && + index < accumLengths.size() - 1; ++index ) + ; + + // Interpolating the cross section at point. + const float alpha = ( length - accumLengths[index] ) / + ( accumLengths[index + 1] - accumLengths[index] ); + const size_t start = range.first + index; + const Vector4f sample = ( *points )[start + 1] * alpha + + ( *points )[start] * (1 - alpha ); + result.push_back( sample ); + } + + return result; +} + +float Morphology::Impl::getDistanceToSoma( const uint32_t sectionID ) const +{ + if( _distancesToSoma.size() <= sectionID ) + _distancesToSoma.resize( sectionID + 1 ); + + float& distance = _distancesToSoma[sectionID]; + if( distance == 0) + { + // This is the soma, a first order section or the distance hasn't + // been computed yet. Soma and first order sections are cheap + // to detect and compute. + const int32_t parent = ( *sections )[sectionID][1]; + if( parent == -1 || ( *types )[parent] == SECTION_SOMA ) + return 0; + // For the other cases it doesn't matter to have concurrent updates + // because they will yield the same result (and it's probably + // cheaper to go ahead with the computation than to contend for a + // mutex). + distance = getSectionLength( parent ) + getDistanceToSoma( parent ); + } + return distance; +} + +floats Morphology::Impl::getSampleDistancesToSoma( + const uint32_t sectionID ) const +{ + const SectionRange range = getSectionRange( sectionID ); + const floats accumLengths = _computeAccumulatedLengths( range ); + floats result; + result.reserve( accumLengths.size( )); + const float distance = getDistanceToSoma( sectionID ); + BOOST_FOREACH( const float length, accumLengths ) + result.push_back( distance + length ); + + return result; +} + +const uint32_ts& Morphology::Impl::getChildren( const uint32_t sectionID ) const +{ + return _sectionChildren[sectionID]; +} + +void Morphology::Impl::transform( const Matrix4f& matrix ) +{ + #pragma omp parallel for + for( size_t i = 0; i < points->size(); ++i) + { + Vector4f& p = ( *points )[i]; + const Vector3f pp = matrix * p.get_sub_vector< 3 >(); + p.get_sub_vector< 3 >() = pp; + } +} + +void Morphology::Impl::_extractChildrenLists() +{ + typedef std::map< uint32_t, uint32_ts > ChildrenMap; + ChildrenMap children; + for( size_t i = 0; i < sections->size(); ++i ) + { + const int32_t parent = ( *sections )[i][1]; + if( parent != -1 ) + children[parent].push_back( i ); + } + _sectionChildren.resize( sections->size( )); + BOOST_FOREACH( ChildrenMap::value_type& sectionAndChildren, + children ) + { + _sectionChildren[sectionAndChildren.first].swap( + sectionAndChildren.second ); + } +} + +float Morphology::Impl::_computeSectionLength( const uint32_t sectionID ) const +{ + const SectionRange range = getSectionRange( sectionID ); + float length = 0; + for( size_t i = range.first; i != range.second - 1; ++i ) + { + const Vector4f& start = ( *points )[i]; + const Vector4f& end = ( *points )[i + 1]; + const Vector3f diff = ( end - start ).get_sub_vector< 3 >(); + length += diff.length( ); + } + return length; +} + +floats Morphology::Impl::_computeAccumulatedLengths( + const SectionRange& range ) const +{ + floats result; + result.reserve( range.second - range.first ); + result.push_back(0); + for( size_t i = range.first; i != range.second - 1; ++i ) + { + const Vector4f& start = ( *points )[i]; + const Vector4f& end = ( *points )[i + 1]; + const Vector3f diff = ( end - start ).get_sub_vector< 3 >(); + result.push_back( result.back() + diff.length( )); + } + return result; +} + +} +} diff --git a/brain/neuron/morphologyImpl.h b/brain/neuron/morphologyImpl.h new file mode 100644 index 0000000..3a7fd38 --- /dev/null +++ b/brain/neuron/morphologyImpl.h @@ -0,0 +1,88 @@ + +/* Copyright (c) 2013-2015, EPFL/Blue Brain Project + * Juan Hernando + * + * This file is part of Brion + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License version 3.0 as published + * by the Free Software Foundation. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + */ + +#ifndef BRAIN_NEURON_MORPHOLOGYIMPL +#define BRAIN_NEURON_MORPHOLOGYIMPL + +#include "morphology.h" + +#include +#include + +namespace brain +{ +namespace neuron +{ + +typedef std::pair< size_t, size_t > SectionRange; + +class Morphology::Impl : public lunchbox::Referenced +{ +public: + const brion::Vector4fsPtr points; + const brion::Vector2isPtr sections; + const brion::SectionTypesPtr types; + const brion::Vector2isPtr apicals; + + Matrix4f transformation; + + uint32_t somaSection; + + Impl( const brion::Morphology& morphology ); + + SectionRange getSectionRange( const uint32_t sectionID ) const; + + uint32_ts getSectionIDs( const SectionTypes& requestedTypes, + bool excludeSoma ) const; + + float getSectionLength( const uint32_t sectionID ) const; + + Vector4fs getSectionSamples( const uint32_t sectionID ) const; + + Vector4fs getSectionSamples( const uint32_t sectionID, + const floats& samplePoints ) const; + + float getDistanceToSoma( const uint32_t sectionID ) const; + + floats getSampleDistancesToSoma( const uint32_t sectionID ) const; + + const uint32_ts& getChildren( const uint32_t sectionID ) const; + + void transform( const Matrix4f& matrix ); + +private: + + // Distances caches. These caches need to be thread-safe to follow the + // recommendations for C++11 about mutable and const correctness. + // (http://herbsutter.com/2013/05/24/gotw-6a-const-correctness-part-1-3/) + typedef lunchbox::LFVector< float > LFFloats; + mutable LFFloats _distancesToSoma; + mutable LFFloats _sectionLengths; + + std::vector< uint32_ts > _sectionChildren; + + void _extractChildrenLists(); + float _computeSectionLength( const uint32_t sectionID ) const; + floats _computeAccumulatedLengths( const SectionRange& range ) const; +}; + +} +} +#endif diff --git a/brain/neuron/section.cpp b/brain/neuron/section.cpp new file mode 100644 index 0000000..f3b866f --- /dev/null +++ b/brain/neuron/section.cpp @@ -0,0 +1,137 @@ + +/* Copyright (c) 2013-2015, EPFL/Blue Brain Project + * Juan Hernando + * + * This file is part of Brion + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License version 3.0 as published + * by the Free Software Foundation. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + */ + +#include "section.h" + +#include "morphology.h" +#include "morphologyImpl.h" + +#include + +namespace brain +{ +namespace neuron +{ + +Section::Section( const uint32_t id, Morphology::Impl* morphology ) + : _id(id) + , _morphology(morphology) +{ + _morphology->ref(); + + SectionRange range = morphology->getSectionRange( id ); + if( range.second <= range.first ) + LBWARN << "Dereferencing broken morphology section " << _id + << std::endl; +} + +Section::~Section() +{ + _morphology->unref(); +} + +Section::Section( const Section& section) + : _id(section._id) + , _morphology(section._morphology) +{ + _morphology->ref(); +} + +Section& Section::operator=( const Section& section) +{ + if( §ion == this ) + return *this; + if( _morphology ) + _morphology->unref(); + _id = section._id; + _morphology = section._morphology; + _morphology->ref(); + return *this; +} + +bool Section::operator==( const Section& other ) const +{ + return other._id == _id && other._morphology == _morphology; +} + +bool Section::operator!=( const Section& other ) const +{ + return !( *this == other ); +} +uint32_t Section::getID() const +{ + return _id; +} + +SectionType Section::getType() const +{ + return ( *_morphology->types )[_id]; +} + +float Section::getLength() const +{ + return _morphology->getSectionLength( _id ); +} + +Vector4fs Section::getSamples() const +{ + return _morphology->getSectionSamples( _id ); +} + +Vector4fs Section::getSamples( const floats& points ) const +{ + return _morphology->getSectionSamples( _id, points ); +} + +float Section::getDistanceToSoma() const +{ + return _morphology->getDistanceToSoma( _id ); +} + +floats Section::getSampleDistancesToSoma() const +{ + return _morphology->getSampleDistancesToSoma( _id ); +} + +bool Section::hasParent() const +{ + const int32_t parent = ( *_morphology->sections )[_id][1]; + return parent != -1 && uint32_t(parent) != _morphology->somaSection; +} +Section Section::getParent() const +{ + const int32_t parent = ( *_morphology->sections )[_id][1]; + if( parent == -1 || uint32_t(parent) == _morphology->somaSection ) + LBTHROW( std::runtime_error( "Cannot access parent section" )); + return Section( parent, _morphology ); +} + +Sections Section::getChildren() const +{ + const uint32_ts& children = _morphology->getChildren( _id ); + Sections result; + result.reserve( children.size( )); + BOOST_FOREACH( uint32_t id, children ) + result.push_back( Section( id, _morphology )); + return result; +} + +} +} diff --git a/brain/neuron/section.h b/brain/neuron/section.h new file mode 100644 index 0000000..e3b76ea --- /dev/null +++ b/brain/neuron/section.h @@ -0,0 +1,146 @@ +/* Copyright (c) 2013-2015, EPFL/Blue Brain Project + * Juan Hernando + * + * This file is part of Brion + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License version 3.0 as published + * by the Free Software Foundation. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + */ + +#ifndef BRAIN_NEURON_SECTION +#define BRAIN_NEURON_SECTION + +#include +#include + +#include + +namespace brain +{ +namespace neuron +{ + +/** + * A class to represent a morphological section. + * + * A Section is an unbranched piece of a morphological skeleton. + * This class provides functions to query information about the sample points + * that compose the section and functions to obtain the parent and children + * sections. + * + * The cell soma is also considered a section, but some functions have + * special meaning for it. + * + * Sections cannot be directly created, but are returned by several + * brain::Morphology and brain::Section methods. + * + * This is a lightweight object with STL container style thread safety. + * It is also safe to use a section after the morphology from where it comes + * has been deallocated. The morphological data will be kept as long as there + * is a Section referring to it. + */ +class Section +{ +public: + friend class Morphology; + friend class Soma; + + BRAIN_API Section( const Section& section ); + + BRAIN_API ~Section(); + + BRAIN_API Section& operator=( const Section& section ); + + BRAIN_API bool operator==( const Section& section ) const; + BRAIN_API bool operator!=( const Section& section ) const; + + /** Return the ID of this section. */ + BRAIN_API uint32_t getID() const; + + /** Return the morphological type of this section (dedrite, axon, ...). */ + BRAIN_API SectionType getType() const; + + /** + * Return the total length of this section in microns. + * + * If this section is a soma section the length is ill-defined and this + * function will return 0. + */ + BRAIN_API float getLength() const; + + /** + * Return the list of all point samples that define this section. + * + * If this sections is a soma section return the list of points of the + * soma profile poly-line. + * + * @return A list of point positions with radius. For a section consisting + * of n segments, this list will have n + 1 points. + */ + BRAIN_API Vector4fs getSamples() const; + + /** + * Return a list of points sampling this section at discrete locations. + * + * If the section is a soma section this function will return the soma + * position for all sampling positions. The soma position is assumed to + * be (0, 0, 0) unless the origin morphology has been transformed. + * + * @param points Normalized positions of the sample points along the + * section. Values will be clampled to [0, 1] before sampling. + * @return The section sampled at the given relative positions. + */ + BRAIN_API Vector4fs getSamples( const floats& points ) const; + + /** + * Return the absolute distance from the start of the section to the soma. + */ + BRAIN_API float getDistanceToSoma() const; + + /** + * Return the absolute distances to the soma in microns for all sample + * positions. + * + * @return A list of distances. For a section consisting + * of n segments, this list will have n + 1 values. The section + * length is equal to the difference between the first and last + * values of the list. + */ + BRAIN_API floats getSampleDistancesToSoma() const; + + /** Return true if this section has a parent section, false otherwise. */ + BRAIN_API bool hasParent() const; + + /** + * Return the parent section of this section. + * @throw runtime_error is the section doesn't have a parent. + */ + BRAIN_API Section getParent() const; + + /** + * Return a vector with all the direct children of this section. + * The container will be empty for terminal sections. + */ + BRAIN_API Sections getChildren() const; + +protected: + BRAIN_API Section( uint32_t id, Morphology::Impl* morphology ); + +private: + uint32_t _id; + Morphology::Impl* _morphology; +}; + +} +} +#endif diff --git a/brain/neuron/soma.cpp b/brain/neuron/soma.cpp new file mode 100644 index 0000000..28f0a65 --- /dev/null +++ b/brain/neuron/soma.cpp @@ -0,0 +1,104 @@ +/* Copyright (c) 2013-2015, EPFL/Blue Brain Project + * Juan Hernando + * + * This file is part of Brion + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License version 3.0 as published + * by the Free Software Foundation. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + */ + +#include "soma.h" +#include "section.h" +#include "morphologyImpl.h" + +#include + +namespace brain +{ +namespace neuron +{ + +namespace +{ + +Vector3f _computeCentroid( const Vector4fs& points ) +{ + Vector3f centroid( Vector3f::ZERO ); + BOOST_FOREACH( const Vector4f& point, points ) + centroid += point.get_sub_vector< 3 >(); + centroid /= float( points.size( )); + return centroid; +} + +} + +Soma::Soma( Morphology::Impl* morphology ) + : _morphology(morphology) +{ + _morphology->ref(); +} + +Soma::~Soma() +{ + _morphology->unref(); +} + +Soma::Soma( const Soma& soma ) + : _morphology(soma._morphology) +{ + _morphology->ref(); +} + +Soma& Soma::operator=( const Soma& soma ) +{ + if( &soma == this ) + return *this; + if( _morphology ) + _morphology->unref(); + _morphology = soma._morphology; + _morphology->ref(); + return *this; +} + +Vector4fs Soma::getProfilePoints() const +{ + return _morphology->getSectionSamples( _morphology->somaSection ); +} + +float Soma::getMeanRadius() const +{ + const Vector4fs points = getProfilePoints(); + const Vector3f centroid = _computeCentroid( points ); + float radius = 0; + BOOST_FOREACH( const Vector4f point, points ) + radius += ( point.get_sub_vector< 3 >() - centroid ).length(); + return radius /= float( points.size( )); +} + +Vector3f Soma::getCentroid() const +{ + return _computeCentroid( getProfilePoints( )); +} + +Sections Soma::getChildren() const +{ + const uint32_ts& children = + _morphology->getChildren( _morphology->somaSection ); + Sections result; + BOOST_FOREACH( uint32_t id, children ) + result.push_back( Section( id, _morphology )); + return result; +} + +} +} diff --git a/brain/neuron/soma.h b/brain/neuron/soma.h new file mode 100644 index 0000000..d3d5176 --- /dev/null +++ b/brain/neuron/soma.h @@ -0,0 +1,87 @@ + +/* Copyright (c) 2013-2015, EPFL/Blue Brain Project + * Juan Hernando + * + * This file is part of Brion + * + * This library is free software; you can redistribute it and/or modify it under + * the terms of the GNU Lesser General Public License version 3.0 as published + * by the Free Software Foundation. + * + * This library is distributed in the hope that it will be useful, but WITHOUT + * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this library; if not, write to the Free Software Foundation, Inc., + * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + */ + +#ifndef BRAIN_NEURON_SOMA +#define BRAIN_NEURON_SOMA + +#include +#include + +#include + +namespace brain +{ +namespace neuron +{ + +/** + * A class to represent a neuron soma. + * + * This class provides functions to query information about the soma of a + * neuron. + * + * Typically the soma is described as the poly-line of the projection + * of the soma onto a plane, where the plane normal points in the vertical + * direction in the local coordinate system of the morphology. In other cases + * the poly-line is not projected onto a plane, but is an approximation of + * the countour of the soma as seen in an orhogonal projection down the + * vertical axis (this is basically the same as before, but the vertical + * coordinate is not 0 for all the points). + * This class can also be used for both descriptions as well as somas simply + * approximated as spheres. + * + * The coordinates system used by a soma will be in the same as the + * brain::Morphology from where it comes. + * + * @version unstable + */ +class Soma +{ +public: + friend class Morphology; + + ~Soma(); + + BRAIN_API Soma( const Soma& soma ); + + BRAIN_API Soma& operator=( const Soma& soma ); + + /** Return the points of the soma profile. */ + BRAIN_API Vector4fs getProfilePoints() const; + + /** Return the mean distance between the profile points and the centroid. */ + BRAIN_API float getMeanRadius() const; + + /** Return the average of the profile points. */ + BRAIN_API Vector3f getCentroid() const; + + /** Return the first order sections starting from the soma. */ + BRAIN_API Sections getChildren() const; + +protected: + BRAIN_API Soma( Morphology::Impl* morhology ); + +private: + Morphology::Impl* _morphology; +}; + +} +} +#endif diff --git a/brain/types.h b/brain/types.h index ee34ef8..63b2740 100644 --- a/brain/types.h +++ b/brain/types.h @@ -26,10 +26,9 @@ namespace brain { +using namespace brion::enums; + class Circuit; -class Morphology; -typedef boost::shared_ptr< Morphology > MorphologyPtr; -typedef std::vector< MorphologyPtr > Morphologies; using vmml::Vector2i; using vmml::Vector3f; @@ -43,13 +42,22 @@ typedef std::vector< Matrix4f > Matrix4fs; using brion::uint32_ts; using brion::floats; -using brion::MorphologyRepairStage; -using brion::SectionType; using brion::SectionTypes; using brion::GIDSet; using brion::URI; typedef std::vector< URI > URIs; +namespace neuron +{ + class Morphology; + typedef boost::shared_ptr< Morphology > MorphologyPtr; + typedef std::vector< MorphologyPtr > Morphologies; + class Section; + typedef std::vector< Section > Sections; + class Soma; +} + + } #endif diff --git a/doc/Changelog.md b/doc/Changelog.md index a13781a..920a292 100644 --- a/doc/Changelog.md +++ b/doc/Changelog.md @@ -5,16 +5,20 @@ Changelog {#Changelog} * [39](https://github.com/BlueBrain/Brion/pull/39): Add compartment report converter tool +* [30](https://github.com/BlueBrain/Brion/pull/30), + [35](https://github.com/BlueBrain/Brion/pull/35): + Added a new library, Brain, to provide higher level functions and classes. + The library provides: + - A Circuit class to get basic information from cells and targets and load + morphologies at circuit level. + - A brain::cell namspace with a Morphology class plus other related classes + to access morphological data about neurons. * [38](https://github.com/BlueBrain/Brion/pull/38): Fix crash while reading more than `ulimit -Sn` (1024 default) morphologies * [37](https://github.com/BlueBrain/Brion/pull/37): Added support for synapse nrn_extra.h5 files. * [31](https://github.com/BlueBrain/Brion/pull/31): Fix crash while reading more than `ulimit -Sn` (1024 default) NEST gdf files -* [30](https://github.com/BlueBrain/Brion/pull/30): - Added a new library, Brain, to provide higher level functions. The library - provides the Circuit and Morphology classes to deal with morphologies at - circuit level. * [29](https://github.com/BlueBrain/Brion/pull/29): New member functions in brion::BlueConfig to provide a semantic API. * [28](https://github.com/BlueBrain/Brion/pull/28): diff --git a/tests/circuit.cpp b/tests/circuit.cpp index b6479af..584ebfc 100644 --- a/tests/circuit.cpp +++ b/tests/circuit.cpp @@ -201,7 +201,7 @@ BOOST_AUTO_TEST_CASE( brain_circuit_positions ) namespace { -void _checkMorphology( const brain::Morphology& morphology, +void _checkMorphology( const brain::neuron::Morphology& morphology, const std::string& other ) { const brion::Morphology reference( @@ -209,11 +209,11 @@ void _checkMorphology( const brain::Morphology& morphology, BOOST_CHECK( morphology.getPoints() == *reference.readPoints( brion::MORPHOLOGY_UNDEFINED )); } -void _checkMorphology( const brain::Morphology& morphology, +void _checkMorphology( const brain::neuron::Morphology& morphology, const std::string& other, const brain::Matrix4f& transform ) { - const brain::Morphology reference( + const brain::neuron::Morphology reference( brion::URI( BBP_TESTDATA + ( "/local/morphologies/01.07.08/h5/" + other)), transform ); @@ -244,7 +244,7 @@ BOOST_AUTO_TEST_CASE( load_local_morphologies ) for( uint32_t gid = 1; gid < 500; gid += 75) gids.insert(gid); // This call also tests brain::Circuit::getMorphologyURIs - const brain::Morphologies morphologies = + const brain::neuron::Morphologies morphologies = circuit.loadMorphologies( gids, brain::Circuit::COORDINATES_LOCAL ); BOOST_CHECK_EQUAL( morphologies.size(), gids.size( )); @@ -256,7 +256,7 @@ BOOST_AUTO_TEST_CASE( load_local_morphologies ) gids.insert(2); gids.insert(4); gids.insert(6); - const brain::Morphologies repeated = + const brain::neuron::Morphologies repeated = circuit.loadMorphologies( gids, brain::Circuit::COORDINATES_LOCAL ); BOOST_CHECK_EQUAL( repeated.size(), gids.size( )); @@ -271,7 +271,7 @@ BOOST_AUTO_TEST_CASE( load_global_morphologies ) brion::GIDSet gids; for( uint32_t gid = 1; gid < 500; gid += 75) gids.insert(gid); - const brain::Morphologies morphologies = + const brain::neuron::Morphologies morphologies = circuit.loadMorphologies( gids, brain::Circuit::COORDINATES_GLOBAL ); BOOST_CHECK_EQUAL( morphologies.size(), gids.size( )); diff --git a/tests/morphology.cpp b/tests/morphology.cpp index 25715e4..e50f75b 100644 --- a/tests/morphology.cpp +++ b/tests/morphology.cpp @@ -40,6 +40,7 @@ #include #include +#include // typdef for brevity typedef brion::Vector4f V4f; @@ -313,8 +314,8 @@ void checkEqualArrays( const std::vector< T >& array, const size_t length, ... ) } template< typename T > -void _checkCloseVectorArrays( const std::vector< T >& array, - const size_t length, va_list args ) +void _checkCloseArrays( const std::vector< T >& array, + const size_t length, va_list args ) { for( size_t i = 0; i != length; ++i ) { @@ -324,27 +325,53 @@ void _checkCloseVectorArrays( const std::vector< T >& array, } template< typename T > -void checkCloseVectorArrays( const std::vector< T >& array, - const size_t length, ... ) +void checkCloseArrays( const std::vector< T >& array, + const size_t length, ... ) { BOOST_CHECK_EQUAL( array.size(), length ); va_list args; va_start( args, length ); - _checkCloseVectorArrays( array, length, args ); + _checkCloseArrays( array, length, args ); va_end( args ); } template< typename T > -void checkCloseVectorArraysUptoN( const std::vector< T >& array, - const size_t length, ... ) +void checkCloseArrays( const std::vector< T >& array1, + const std::vector< T >& array2 ) +{ + BOOST_CHECK_EQUAL( array1.size(), array2.size() ); + for( size_t i = 0; i != std::min( array1.size(), array2.size( )); ++i ) + BOOST_CHECK_CLOSE( array1[i], array2[i], 2e-5f); +} + +template< typename T, long unsigned int M > +void checkCloseArrays( const std::vector< vmml::vector< M, T > >& array1, + const std::vector< vmml::vector< M, T > >& array2 ) +{ + BOOST_CHECK_EQUAL( array1.size(), array2.size() ); + for( size_t i = 0; i != std::min( array1.size(), array2.size( )); ++i ) + BOOST_CHECK_SMALL(( array1[i] - array2[i] ).length( ), 0.00001f); +} + +template< typename T > +void checkCloseArraysUptoN( const std::vector< T >& array, + const size_t length, ... ) { BOOST_CHECK( array.size() >= length ); va_list args; va_start( args, length ); - _checkCloseVectorArrays( array, length, args ); + _checkCloseArrays( array, length, args ); va_end( args ); } +brion::uint32_ts getSectionIDs( const brain::neuron::Sections& sections ) +{ + brion::uint32_ts result; + BOOST_FOREACH( const brain::neuron::Section& section, sections ) + result.push_back( section.getID( )); + return result; +} + BOOST_AUTO_TEST_CASE( swc_soma ) { boost::filesystem::path path( BRION_TESTDATA ); @@ -528,7 +555,7 @@ BOOST_AUTO_TEST_CASE( swc_neuron ) namespace { -void checkEqualMorphologies( const brain::Morphology& first, +void checkEqualMorphologies( const brain::neuron::Morphology& first, const brion::Morphology& second ) { BOOST_CHECK( *second.readPoints( brion::MORPHOLOGY_UNDEFINED ) == @@ -545,20 +572,24 @@ BOOST_AUTO_TEST_CASE( v2_morphology_constructors ) boost::shared_ptr< brion::Morphology > raw( new brion::Morphology( TEST_MORPHOLOGY_FILENAME )); - checkEqualMorphologies( brain::Morphology( TEST_MORPHOLOGY_URI ), *raw ); - checkEqualMorphologies( brain::Morphology( *raw ), *raw ); + brain::neuron::Morphology morphology( TEST_MORPHOLOGY_URI ); + BOOST_CHECK_EQUAL( morphology.getTransformation(), + brain::Matrix4f::IDENTITY ); + checkEqualMorphologies( morphology, *raw ); + checkEqualMorphologies( brain::neuron::Morphology( *raw ), *raw ); - BOOST_CHECK_THROW( brain::Morphology( brion::URI( "/mars" )), + BOOST_CHECK_THROW( brain::neuron::Morphology( brion::URI( "/mars" )), std::runtime_error); } BOOST_AUTO_TEST_CASE( get_section_ids ) { - brain::Morphology morphology( TEST_MORPHOLOGY_URI ); + brain::neuron::Morphology morphology( TEST_MORPHOLOGY_URI ); brion::SectionTypes types; types.push_back( brion::SECTION_SOMA ); checkEqualArrays( morphology.getSectionIDs( types ), 1, 0 ); + types.push_back( brion::SECTION_DENDRITE ); checkEqualArrays( morphology.getSectionIDs( types ), 7, 0, 4, 5, 6, 7, 8, 9 ); @@ -570,62 +601,187 @@ BOOST_AUTO_TEST_CASE( get_section_ids ) types.push_back( brion::SECTION_DENDRITE ); checkEqualArrays( morphology.getSectionIDs( types ), 9, 1, 2, 3, 4, 5, 6, 7, 8, 9 ); +} + +BOOST_AUTO_TEST_CASE( get_sections ) +{ + brain::neuron::Morphology morphology( TEST_MORPHOLOGY_URI ); + + BOOST_CHECK_THROW( morphology.getSection( 0 ), std::runtime_error ); + + for( size_t i = 1; i < 13; ++i ) + BOOST_CHECK_EQUAL( morphology.getSection( i ).getID(), i ); + + brain::neuron::Section section = morphology.getSection( 1 ); + BOOST_CHECK( section == morphology.getSection( 1 )); + section = morphology.getSection( 2 ); + BOOST_CHECK( section != morphology.getSection( 1 )); + BOOST_CHECK( section == morphology.getSection( 2 )); + for( size_t i = 1; i < 4; ++i ) + BOOST_CHECK_EQUAL( morphology.getSection( i ).getType(), + brion::SECTION_AXON ); + for( size_t i = 4; i < 10; ++i ) + BOOST_CHECK_EQUAL( morphology.getSection( i ).getType(), + brion::SECTION_DENDRITE ); + for( size_t i = 10; i < 13; ++i ) + BOOST_CHECK_EQUAL( morphology.getSection( i ).getType(), + brion::SECTION_APICAL_DENDRITE ); } BOOST_AUTO_TEST_CASE( get_section_samples ) { - brain::Morphology morphology( TEST_MORPHOLOGY_URI ); + brain::neuron::Morphology morphology( TEST_MORPHOLOGY_URI ); + + brion::Vector4fs points; + for( size_t i = 0; i != 11; ++i) + { + float i2 = i * i; + points.push_back( + brion::Vector4f(0, -i2 / 20.0, i2 / 20.0, 0.5 + i2 /1000.0)); + } + checkCloseArrays( morphology.getSection( 1 ).getSamples(), points ); + + points.clear(); + for( size_t i = 0; i != 11; ++i) + { + float i2 = i * i; + points.push_back( + brion::Vector4f(i2 / 20.0, 0, i2 / 20.0, 0.5 + i2 /1000.0)); + } + checkCloseArrays( morphology.getSection( 4 ).getSamples(), points ); + + points.clear(); + for( size_t i = 0; i != 11; ++i) + { + float i2 = i * i; + points.push_back( + brion::Vector4f(-i2 / 20.0, 0, i2 / 20.0, 0.5 + i2 /1000.0)); + } + checkCloseArrays( morphology.getSection( 7 ).getSamples(), points ); + + points.clear(); + for( size_t i = 0; i != 11; ++i) + { + float i2 = i * i; + points.push_back( + brion::Vector4f(0, i2 / 20.0, i2 / 20.0, 0.5 + i2 /1000.0)); + } + checkCloseArrays( morphology.getSection( 10 ).getSamples(), points ); +} + +BOOST_AUTO_TEST_CASE( get_section_distances_to_soma ) +{ + brain::neuron::Morphology morphology( TEST_MORPHOLOGY_URI ); + + uint32_t sections[] = {1, 4, 7, 10}; + + for( size_t i = 0; i != 4; ++i) + { + uint32_t section = sections[i]; + BOOST_CHECK_EQUAL( + morphology.getSection( section ).getDistanceToSoma(), 0 ); + const float length = std::sqrt( 5 * 5 * 2 ); + BOOST_CHECK_CLOSE( + morphology.getSection( section ).getLength(), length, 1e-5 ); + + // The distance to the soma of the next section is equal to the length + // of its parent + BOOST_CHECK_CLOSE( + morphology.getSection( section + 1 ).getDistanceToSoma(), + length, 1e-5 ); + + brion::floats reference; + for( size_t j = 0; j != 11; ++j) + { + const float p = j*j / 20.0; + reference.push_back( std::sqrt( p * p * 2 )); + } + checkCloseArrays( + morphology.getSection( section ).getSampleDistancesToSoma( ), + reference ); + } +} + +BOOST_AUTO_TEST_CASE( get_soma_geomery ) +{ + brain::neuron::Morphology morphology( TEST_MORPHOLOGY_URI ); + + const brain::neuron::Soma soma = morphology.getSoma(); + checkEqualArrays( soma.getProfilePoints(), 4, + V4f( .1, 0, 0, .1 ), V4f( 0, .1, 0, .1 ), + V4f( -.1, 0, 0, .1 ), V4f( 0, -.1, 0, .1 )); + + BOOST_CHECK_CLOSE( soma.getMeanRadius(), 0.1, 1e-5 ); + BOOST_CHECK_EQUAL( soma.getCentroid(), V3f::ZERO ); + + brain::Matrix4f matrix( brain::Matrix4f::IDENTITY ); + matrix.set_translation( V3f( 2, 0, 0 )); + brain::neuron::Morphology transformed( TEST_MORPHOLOGY_URI, matrix ); + BOOST_CHECK_EQUAL( transformed.getSoma().getCentroid(), V3f( 2, 0, 0 )); + +} + +BOOST_AUTO_TEST_CASE( get_section_samples_by_positions ) +{ + brain::neuron::Morphology morphology( TEST_MORPHOLOGY_URI ); brion::floats points; for( float p = 0.0; p <= 1.0; p += 0.2 ) points.push_back( p ); - checkCloseVectorArrays( morphology.getSectionSamples( 1, points ), 6, + checkCloseArrays( morphology.getSection( 1 ).getSamples( points ), 6, V4f( 0, 0, 0, .5 ), V4f( 0, -1, 1, .52 ), V4f( 0, -2, 2, .54 ), V4f( 0, -3, 3, .56 ), V4f( 0, -4, 4, .58 ), V4f( 0, -5, 5, .6 )); - checkCloseVectorArrays( morphology.getSectionSamples( 4, points ), 6, + checkCloseArrays( morphology.getSection( 4 ).getSamples( points ), 6, V4f( 0, 0, 0, .5 ), V4f( 1, 0, 1, .52 ), V4f( 2, 0, 2, .54 ), V4f( 3, 0, 3, .56 ), V4f( 4, 0, 4, .58 ), V4f( 5, 0, 5, .6 )); - checkCloseVectorArrays( morphology.getSectionSamples( 7, points ), 6, + checkCloseArrays( morphology.getSection( 7 ).getSamples( points ), 6, V4f( 0, 0, 0, .5 ), V4f( -1, 0, 1, .52 ), V4f( -2, 0, 2, .54 ), V4f( -3, 0, 3, .56 ), V4f( -4, 0, 4, .58 ), V4f( -5, 0, 5, .6 )); - checkCloseVectorArrays( morphology.getSectionSamples( 10, points ), 6, + checkCloseArrays( morphology.getSection( 10 ).getSamples( points ), 6, V4f( 0, 0, 0, .5 ), V4f( 0, 1, 1, .52 ), V4f( 0, 2, 2, .54 ), V4f( 0, 3, 3, .56 ), V4f( 0, 4, 4, .58 ), V4f( 0, 5, 5, .6 )); } -BOOST_AUTO_TEST_CASE( soma_position ) +BOOST_AUTO_TEST_CASE( morphology_hierarchy ) { - brain::Morphology local( TEST_MORPHOLOGY_URI ); - brion::floats points; - points.push_back( 0.232 ); // arbitrary - BOOST_CHECK_EQUAL( local.getSectionSamples( 0, points )[0], V3f::ZERO ); + brain::neuron::Morphology morphology( TEST_MORPHOLOGY_URI ); - brain::Matrix4f matrix( brain::Matrix4f::IDENTITY ); - matrix.set_translation( V3f( 2, 0, 0 )); - brain::Morphology transformed( TEST_MORPHOLOGY_URI, matrix ); - BOOST_CHECK_EQUAL( transformed.getSectionSamples( 0, points )[0], - V3f( 2, 0, 0 )); + BOOST_CHECK( !morphology.getSection( 1 ).hasParent( )); + BOOST_CHECK( !morphology.getSection( 4 ).hasParent( )); + BOOST_CHECK_EQUAL( morphology.getSection( 2 ).getParent().getID(), 1 ); + BOOST_CHECK_EQUAL( morphology.getSection( 3 ).getParent().getID(), 1 ); + BOOST_CHECK_EQUAL( morphology.getSection( 5 ).getParent().getID(), 4 ); + BOOST_CHECK_EQUAL( morphology.getSection( 6 ).getParent().getID(), 4 ); + + checkEqualArrays( getSectionIDs( morphology.getSoma().getChildren( )), + 4, 1, 4, 7, 10 ); + checkEqualArrays( getSectionIDs( morphology.getSection( 1 ).getChildren( )), + 2, 2, 3 ); + checkEqualArrays( getSectionIDs( morphology.getSection( 4 ).getChildren( )), + 2, 5, 6 ); + BOOST_CHECK( morphology.getSection( 5 ).getChildren().empty( )); } BOOST_AUTO_TEST_CASE( transform_with_matrix ) { brain::Matrix4f matrix( brain::Matrix4f::IDENTITY ); matrix.rotate_z( M_PI * 0.5 ); - brain::Morphology rotated( TEST_MORPHOLOGY_URI, matrix ); - checkCloseVectorArraysUptoN( rotated.getPoints(), 4, + brain::neuron::Morphology rotated( TEST_MORPHOLOGY_URI, matrix ); + checkCloseArraysUptoN( rotated.getPoints(), 4, V4f( .0, .1, .0, .1 ), V4f( -.1, .0, .0, .1 ), V4f( .0, -.1, .0, .1 ), V4f( .1, .0, .0, .1 )); matrix = brain::Matrix4f::IDENTITY; matrix.rotate_z( M_PI * 0.5 ); matrix.set_translation( V3f( 2, 0, 0 )); - brain::Morphology transformed( TEST_MORPHOLOGY_URI, matrix ); - checkCloseVectorArraysUptoN( transformed.getPoints(), 4, + brain::neuron::Morphology transformed( TEST_MORPHOLOGY_URI, matrix ); + BOOST_CHECK_EQUAL( transformed.getTransformation(), matrix ); + checkCloseArraysUptoN( transformed.getPoints(), 4, V4f( 2., .1, .0, .1 ), V4f( 1.9, .0, .0, .1 ), V4f( 2., -.1, .0, .1 ), V4f( 2.1, .0, .0, .1 )); }