diff --git a/Modules/DiffusionImaging/Connectomics/Algorithms/mitkConnectomicsNetworkCreator.cpp b/Modules/DiffusionImaging/Connectomics/Algorithms/mitkConnectomicsNetworkCreator.cpp index 17b866066d..e71c975016 100644 --- a/Modules/DiffusionImaging/Connectomics/Algorithms/mitkConnectomicsNetworkCreator.cpp +++ b/Modules/DiffusionImaging/Connectomics/Algorithms/mitkConnectomicsNetworkCreator.cpp @@ -1,847 +1,864 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkConnectomicsNetworkCreator.h" #include #include #include "mitkConnectomicsConstantsManager.h" #include "mitkImageAccessByItk.h" #include "mitkImageStatisticsHolder.h" #include "mitkImageCast.h" #include "itkImageRegionIteratorWithIndex.h" // VTK #include #include #include mitk::ConnectomicsNetworkCreator::ConnectomicsNetworkCreator() : m_FiberBundle() , m_Segmentation() , m_ConNetwork( mitk::ConnectomicsNetwork::New() ) , idCounter(0) , m_LabelToVertexMap() , m_LabelToNodePropertyMap() , allowLoops( false ) , m_UseCoMCoordinates( false ) , m_LabelsToCoordinatesMap() , m_MappingStrategy( EndElementPositionAvoidingWhiteMatter ) , m_EndPointSearchRadius( 10.0 ) +, m_ZeroLabelInvalid( true ) +, m_AbortConnection( false ) { } mitk::ConnectomicsNetworkCreator::ConnectomicsNetworkCreator( mitk::Image::Pointer segmentation, mitk::FiberBundleX::Pointer fiberBundle ) : m_FiberBundle(fiberBundle) , m_Segmentation(segmentation) , m_ConNetwork( mitk::ConnectomicsNetwork::New() ) , idCounter(0) , m_LabelToVertexMap() , m_LabelToNodePropertyMap() , allowLoops( false ) , m_LabelsToCoordinatesMap() , m_MappingStrategy( EndElementPositionAvoidingWhiteMatter ) , m_EndPointSearchRadius( 10.0 ) +, m_ZeroLabelInvalid( true ) +, m_AbortConnection( false ) { mitk::CastToItkImage( segmentation, m_SegmentationItk ); } mitk::ConnectomicsNetworkCreator::~ConnectomicsNetworkCreator() { } void mitk::ConnectomicsNetworkCreator::SetFiberBundle(mitk::FiberBundleX::Pointer fiberBundle) { m_FiberBundle = fiberBundle; } void mitk::ConnectomicsNetworkCreator::SetSegmentation(mitk::Image::Pointer segmentation) { m_Segmentation = segmentation; mitk::CastToItkImage( segmentation, m_SegmentationItk ); } itk::Point mitk::ConnectomicsNetworkCreator::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } void mitk::ConnectomicsNetworkCreator::CreateNetworkFromFibersAndSegmentation() { //empty graph m_ConNetwork->clear(); m_LabelToVertexMap.clear(); m_LabelToNodePropertyMap.clear(); idCounter = 0; vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData(); vtkSmartPointer vLines = fiberPolyData->GetLines(); vLines->InitTraversal(); int numFibers = m_FiberBundle->GetNumFibers(); for( int fiberID( 0 ); fiberID < numFibers; fiberID++ ) { vtkIdType numPointsInCell(0); vtkIdType* pointsInCell(NULL); vLines->GetNextCell ( numPointsInCell, pointsInCell ); TractType::Pointer singleTract = TractType::New(); for( int pointInCellID( 0 ); pointInCellID < numPointsInCell ; pointInCellID++) { // push back point PointType point = GetItkPoint( fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ) ); singleTract->InsertElement( singleTract->Size(), point ); } if ( singleTract && ( singleTract->Size() > 0 ) ) { AddConnectionToNetwork( ReturnAssociatedVertexPairForLabelPair( ReturnLabelForFiberTract( singleTract, m_MappingStrategy ) ) ); + m_AbortConnection = false; } } // Prune unconnected nodes //m_ConNetwork->PruneUnconnectedSingleNodes(); // provide network with geometry m_ConNetwork->SetGeometry( dynamic_cast(m_Segmentation->GetGeometry()->Clone().GetPointer()) ); m_ConNetwork->UpdateBounds(); m_ConNetwork->SetIsModified( true ); MBI_INFO << mitk::ConnectomicsConstantsManager::CONNECTOMICS_WARNING_INFO_NETWORK_CREATED; } void mitk::ConnectomicsNetworkCreator::AddConnectionToNetwork(ConnectionType newConnection) { + if( m_AbortConnection ) + { + MITK_DEBUG << "Connection aborted"; + return; + } + VertexType vertexA = newConnection.first; VertexType vertexB = newConnection.second; // check for loops (if they are not allowed) if( allowLoops || !( vertexA == vertexB ) ) { // If the connection already exists, increment weight, else create connection if ( m_ConNetwork->EdgeExists( vertexA, vertexB ) ) { m_ConNetwork->IncreaseEdgeWeight( vertexA, vertexB ); } else { m_ConNetwork->AddEdge( vertexA, vertexB ); } } } mitk::ConnectomicsNetworkCreator::VertexType mitk::ConnectomicsNetworkCreator::ReturnAssociatedVertexForLabel( ImageLabelType label ) { + if( m_ZeroLabelInvalid && ( label == 0 ) ) + { + m_AbortConnection = true; + return NULL; + } + // if label is not known, create entry if( ! ( m_LabelToVertexMap.count( label ) > 0 ) ) { VertexType newVertex = m_ConNetwork->AddVertex( idCounter ); idCounter++; SupplyVertexWithInformation(label, newVertex); m_LabelToVertexMap.insert( std::pair< ImageLabelType, VertexType >( label, newVertex ) ); } //return associated vertex return m_LabelToVertexMap.find( label )->second; } mitk::ConnectomicsNetworkCreator::ConnectionType mitk::ConnectomicsNetworkCreator::ReturnAssociatedVertexPairForLabelPair( ImageLabelPairType labelpair ) { //hand both labels through to the single label function ConnectionType connection( ReturnAssociatedVertexForLabel(labelpair.first), ReturnAssociatedVertexForLabel(labelpair.second) ); return connection; } mitk::ConnectomicsNetworkCreator::ImageLabelPairType mitk::ConnectomicsNetworkCreator::ReturnLabelForFiberTract( TractType::Pointer singleTract, mitk::ConnectomicsNetworkCreator::MappingStrategy strategy) { switch( strategy ) { case EndElementPosition: { return EndElementPositionLabel( singleTract ); } case PrecomputeAndDistance: { return PrecomputeVertexLocationsBySegmentation( singleTract ); } case JustEndPointVerticesNoLabel: { return JustEndPointVerticesNoLabelTest( singleTract ); } case EndElementPositionAvoidingWhiteMatter: { return EndElementPositionLabelAvoidingWhiteMatter( singleTract ); } } // To remove warnings, this code should never be reached MBI_ERROR << mitk::ConnectomicsConstantsManager::CONNECTOMICS_ERROR_INVALID_MAPPING; ImageLabelPairType nullPair( NULL, NULL ); return nullPair; } mitk::ConnectomicsNetworkCreator::ImageLabelPairType mitk::ConnectomicsNetworkCreator::EndElementPositionLabel( TractType::Pointer singleTract ) { ImageLabelPairType labelpair; {// Note: .fib image tracts are safed using index coordinates mitk::Point3D firstElementFiberCoord, lastElementFiberCoord; mitk::Point3D firstElementSegCoord, lastElementSegCoord; mitk::Index3D firstElementSegIndex, lastElementSegIndex; if( singleTract->front().Size() != 3 ) { MBI_ERROR << mitk::ConnectomicsConstantsManager::CONNECTOMICS_ERROR_INVALID_DIMENSION_NEED_3; } for( int index = 0; index < singleTract->front().Size(); index++ ) { firstElementFiberCoord.SetElement( index, singleTract->front().GetElement( index ) ); lastElementFiberCoord.SetElement( index, singleTract->back().GetElement( index ) ); } // convert from fiber index coordinates to segmentation index coordinates FiberToSegmentationCoords( firstElementFiberCoord, firstElementSegCoord ); FiberToSegmentationCoords( lastElementFiberCoord, lastElementSegCoord ); for( int index = 0; index < 3; index++ ) { firstElementSegIndex.SetElement( index, firstElementSegCoord.GetElement( index ) ); lastElementSegIndex.SetElement( index, lastElementSegCoord.GetElement( index ) ); } int firstLabel = m_SegmentationItk->GetPixel(firstElementSegIndex); int lastLabel = m_SegmentationItk->GetPixel(lastElementSegIndex ); labelpair.first = firstLabel; labelpair.second = lastLabel; // Add property to property map CreateNewNode( firstLabel, firstElementSegIndex, m_UseCoMCoordinates ); CreateNewNode( lastLabel, lastElementSegIndex, m_UseCoMCoordinates ); } return labelpair; } mitk::ConnectomicsNetworkCreator::ImageLabelPairType mitk::ConnectomicsNetworkCreator::PrecomputeVertexLocationsBySegmentation( TractType::Pointer singleTract ) { ImageLabelPairType labelpair; return labelpair; } mitk::ConnectomicsNetworkCreator::ImageLabelPairType mitk::ConnectomicsNetworkCreator::EndElementPositionLabelAvoidingWhiteMatter( TractType::Pointer singleTract ) { ImageLabelPairType labelpair; {// Note: .fib image tracts are safed using index coordinates mitk::Point3D firstElementFiberCoord, lastElementFiberCoord; mitk::Point3D firstElementSegCoord, lastElementSegCoord; mitk::Index3D firstElementSegIndex, lastElementSegIndex; if( singleTract->front().Size() != 3 ) { MBI_ERROR << mitk::ConnectomicsConstantsManager::CONNECTOMICS_ERROR_INVALID_DIMENSION_NEED_3; } for( int index = 0; index < singleTract->front().Size(); index++ ) { firstElementFiberCoord.SetElement( index, singleTract->front().GetElement( index ) ); lastElementFiberCoord.SetElement( index, singleTract->back().GetElement( index ) ); } // convert from fiber index coordinates to segmentation index coordinates FiberToSegmentationCoords( firstElementFiberCoord, firstElementSegCoord ); FiberToSegmentationCoords( lastElementFiberCoord, lastElementSegCoord ); for( int index = 0; index < 3; index++ ) { firstElementSegIndex.SetElement( index, firstElementSegCoord.GetElement( index ) ); lastElementSegIndex.SetElement( index, lastElementSegCoord.GetElement( index ) ); } int firstLabel = m_SegmentationItk->GetPixel(firstElementSegIndex); int lastLabel = m_SegmentationItk->GetPixel(lastElementSegIndex ); // Check whether the labels belong to the white matter (which means, that the fibers ended early) bool extendFront(false), extendEnd(false), retractFront(false), retractEnd(false); extendFront = !IsNonWhiteMatterLabel( firstLabel ); extendEnd = !IsNonWhiteMatterLabel( lastLabel ); retractFront = IsBackgroundLabel( firstLabel ); retractEnd = IsBackgroundLabel( lastLabel ); //if( extendFront || extendEnd ) //{ //MBI_INFO << "Before Start: " << firstLabel << " at " << firstElementSegIndex[ 0 ] << " " << firstElementSegIndex[ 1 ] << " " << firstElementSegIndex[ 2 ] << " End: " << lastLabel << " at " << lastElementSegIndex[ 0 ] << " " << lastElementSegIndex[ 1 ] << " " << lastElementSegIndex[ 2 ]; //} if ( extendFront ) { std::vector< int > indexVectorOfPointsToUse; //Use first two points for direction indexVectorOfPointsToUse.push_back( 1 ); indexVectorOfPointsToUse.push_back( 0 ); // label and coordinate temp storage int tempLabel( firstLabel ); mitk::Index3D tempIndex = firstElementSegIndex; LinearExtensionUntilGreyMatter( indexVectorOfPointsToUse, singleTract, tempLabel, tempIndex ); firstLabel = tempLabel; firstElementSegIndex = tempIndex; } if ( extendEnd ) { std::vector< int > indexVectorOfPointsToUse; //Use last two points for direction indexVectorOfPointsToUse.push_back( singleTract->Size() - 2 ); indexVectorOfPointsToUse.push_back( singleTract->Size() - 1 ); // label and coordinate temp storage int tempLabel( lastLabel ); mitk::Index3D tempIndex = lastElementSegIndex; LinearExtensionUntilGreyMatter( indexVectorOfPointsToUse, singleTract, tempLabel, tempIndex ); lastLabel = tempLabel; lastElementSegIndex = tempIndex; } if ( retractFront ) { // label and coordinate temp storage int tempLabel( firstLabel ); mitk::Index3D tempIndex = firstElementSegIndex; RetractionUntilBrainMatter( true, singleTract, tempLabel, tempIndex ); firstLabel = tempLabel; firstElementSegIndex = tempIndex; } if ( retractEnd ) { // label and coordinate temp storage int tempLabel( lastLabel ); mitk::Index3D tempIndex = lastElementSegIndex; RetractionUntilBrainMatter( false, singleTract, tempLabel, tempIndex ); lastLabel = tempLabel; lastElementSegIndex = tempIndex; } //if( extendFront || extendEnd ) //{ // MBI_INFO << "After Start: " << firstLabel << " at " << firstElementSegIndex[ 0 ] << " " << firstElementSegIndex[ 1 ] << " " << firstElementSegIndex[ 2 ] << " End: " << lastLabel << " at " << lastElementSegIndex[ 0 ] << " " << lastElementSegIndex[ 1 ] << " " << lastElementSegIndex[ 2 ]; //} labelpair.first = firstLabel; labelpair.second = lastLabel; // Add property to property map CreateNewNode( firstLabel, firstElementSegIndex, m_UseCoMCoordinates ); CreateNewNode( lastLabel, lastElementSegIndex, m_UseCoMCoordinates ); } return labelpair; } mitk::ConnectomicsNetworkCreator::ImageLabelPairType mitk::ConnectomicsNetworkCreator::JustEndPointVerticesNoLabelTest( TractType::Pointer singleTract ) { ImageLabelPairType labelpair; {// Note: .fib image tracts are safed using index coordinates mitk::Point3D firstElementFiberCoord, lastElementFiberCoord; mitk::Point3D firstElementSegCoord, lastElementSegCoord; mitk::Index3D firstElementSegIndex, lastElementSegIndex; if( singleTract->front().Size() != 3 ) { MBI_ERROR << mitk::ConnectomicsConstantsManager::CONNECTOMICS_ERROR_INVALID_DIMENSION_NEED_3; } for( int index = 0; index < singleTract->front().Size(); index++ ) { firstElementFiberCoord.SetElement( index, singleTract->front().GetElement( index ) ); lastElementFiberCoord.SetElement( index, singleTract->back().GetElement( index ) ); } // convert from fiber index coordinates to segmentation index coordinates FiberToSegmentationCoords( firstElementFiberCoord, firstElementSegCoord ); FiberToSegmentationCoords( lastElementFiberCoord, lastElementSegCoord ); for( int index = 0; index < 3; index++ ) { firstElementSegIndex.SetElement( index, firstElementSegCoord.GetElement( index ) ); lastElementSegIndex.SetElement( index, lastElementSegCoord.GetElement( index ) ); } int firstLabel = 1 * firstElementSegIndex[ 0 ] + 1000 * firstElementSegIndex[ 1 ] + 1000000 * firstElementSegIndex[ 2 ]; int lastLabel = 1 * firstElementSegIndex[ 0 ] + 1000 * firstElementSegIndex[ 1 ] + 1000000 * firstElementSegIndex[ 2 ]; labelpair.first = firstLabel; labelpair.second = lastLabel; // Add property to property map CreateNewNode( firstLabel, firstElementSegIndex, m_UseCoMCoordinates ); CreateNewNode( lastLabel, lastElementSegIndex, m_UseCoMCoordinates ); } return labelpair; } void mitk::ConnectomicsNetworkCreator::SupplyVertexWithInformation( ImageLabelType& label, VertexType& vertex ) { // supply a vertex with the additional information belonging to the label // TODO: Implement additional information acquisition m_ConNetwork->SetLabel( vertex, m_LabelToNodePropertyMap.find( label )->second.label ); m_ConNetwork->SetCoordinates( vertex, m_LabelToNodePropertyMap.find( label )->second.coordinates ); } std::string mitk::ConnectomicsNetworkCreator::LabelToString( ImageLabelType& label ) { int tempInt = (int) label; std::stringstream ss;//create a stringstream std::string tempString; ss << tempInt;//add number to the stream tempString = ss.str(); return tempString;//return a string with the contents of the stream } mitk::ConnectomicsNetwork::Pointer mitk::ConnectomicsNetworkCreator::GetNetwork() { return m_ConNetwork; } void mitk::ConnectomicsNetworkCreator::FiberToSegmentationCoords( mitk::Point3D& fiberCoord, mitk::Point3D& segCoord ) { mitk::Point3D tempPoint; // convert from fiber index coordinates to segmentation index coordinates m_FiberBundle->GetGeometry()->IndexToWorld( fiberCoord, tempPoint ); m_Segmentation->GetGeometry()->WorldToIndex( tempPoint, segCoord ); } void mitk::ConnectomicsNetworkCreator::SegmentationToFiberCoords( mitk::Point3D& segCoord, mitk::Point3D& fiberCoord ) { mitk::Point3D tempPoint; // convert from fiber index coordinates to segmentation index coordinates m_Segmentation->GetGeometry()->IndexToWorld( segCoord, tempPoint ); m_FiberBundle->GetGeometry()->WorldToIndex( tempPoint, fiberCoord ); } bool mitk::ConnectomicsNetworkCreator::IsNonWhiteMatterLabel( int labelInQuestion ) { bool isWhite( false ); isWhite = ( ( labelInQuestion == freesurfer_Left_Cerebral_White_Matter ) || ( labelInQuestion == freesurfer_Left_Cerebellum_White_Matter ) || ( labelInQuestion == freesurfer_Right_Cerebral_White_Matter ) || ( labelInQuestion == freesurfer_Right_Cerebellum_White_Matter )|| ( labelInQuestion == freesurfer_Left_Cerebellum_Cortex ) || ( labelInQuestion == freesurfer_Right_Cerebellum_Cortex ) || ( labelInQuestion == freesurfer_Brain_Stem ) ); return !isWhite; } bool mitk::ConnectomicsNetworkCreator::IsBackgroundLabel( int labelInQuestion ) { bool isBackground( false ); isBackground = ( labelInQuestion == 0 ); return isBackground; } void mitk::ConnectomicsNetworkCreator::LinearExtensionUntilGreyMatter( std::vector & indexVectorOfPointsToUse, TractType::Pointer singleTract, int & label, mitk::Index3D & mitkIndex ) { if( indexVectorOfPointsToUse.size() > singleTract->Size() ) { MBI_WARN << mitk::ConnectomicsConstantsManager::CONNECTOMICS_WARNING_MORE_POINTS_THAN_PRESENT; return; } if( indexVectorOfPointsToUse.size() < 2 ) { MBI_WARN << mitk::ConnectomicsConstantsManager::CONNECTOMICS_WARNING_ESTIMATING_LESS_THAN_2; return; } for( int index( 0 ); index < indexVectorOfPointsToUse.size(); index++ ) { if( indexVectorOfPointsToUse[ index ] > singleTract->Size() ) { MBI_WARN << mitk::ConnectomicsConstantsManager::CONNECTOMICS_WARNING_ESTIMATING_BEYOND_END; return; } if( indexVectorOfPointsToUse[ index ] < 0 ) { MBI_WARN << mitk::ConnectomicsConstantsManager::CONNECTOMICS_WARNING_ESTIMATING_BEYOND_START; return; } } mitk::Point3D startPoint, endPoint; std::vector< double > differenceVector; differenceVector.resize( singleTract->front().Size() ); { // which points to use, currently only last two //TODO correct using all points int endPointIndex = indexVectorOfPointsToUse.size() - 1; int startPointIndex = indexVectorOfPointsToUse.size() - 2; // convert to segmentation coords mitk::Point3D startFiber, endFiber; for( int index = 0; index < singleTract->front().Size(); index++ ) { endFiber.SetElement( index, singleTract->GetElement( indexVectorOfPointsToUse[ endPointIndex ] ).GetElement( index ) ); startFiber.SetElement( index, singleTract->GetElement( indexVectorOfPointsToUse[ startPointIndex ] ).GetElement( index ) ); } FiberToSegmentationCoords( endFiber, endPoint ); FiberToSegmentationCoords( startFiber, startPoint ); // calculate straight line for( int index = 0; index < singleTract->front().Size(); index++ ) { differenceVector[ index ] = endPoint.GetElement( index ) - startPoint.GetElement( index ); } // normalizing direction vector double length( 0.0 ); double sum( 0.0 ); for( int index = 0; index < differenceVector.size() ; index++ ) { sum = sum + differenceVector[ index ] * differenceVector[ index ]; } length = std::sqrt( sum ); for( int index = 0; index < differenceVector.size() ; index++ ) { differenceVector[ index ] = differenceVector[ index ] / length; } // follow line until first non white matter label mitk::Index3D tempIndex; int tempLabel( label ); bool keepOn( true ); itk::ImageRegion<3> itkRegion = m_SegmentationItk->GetLargestPossibleRegion(); for( int parameter( 0 ) ; keepOn ; parameter++ ) { if( parameter > 1000 ) { MBI_WARN << mitk::ConnectomicsConstantsManager::CONNECTOMICS_WARNING_DID_NOT_FIND_WHITE; break; } for( int index( 0 ); index < 3; index++ ) { tempIndex.SetElement( index, endPoint.GetElement( index ) + parameter * differenceVector[ index ] ); } if( itkRegion.IsInside( tempIndex ) ) { tempLabel = m_SegmentationItk->GetPixel( tempIndex ); } else { tempLabel = -1; } if( IsNonWhiteMatterLabel( tempLabel ) ) { if( tempLabel < 1 ) { keepOn = false; //MBI_WARN << mitk::ConnectomicsConstantsManager::CONNECTOMICS_WARNING_NOT_EXTEND_TO_WHITE; } else { label = tempLabel; mitkIndex = tempIndex; keepOn = false; } } } } } void mitk::ConnectomicsNetworkCreator::RetractionUntilBrainMatter( bool retractFront, TractType::Pointer singleTract, int & label, mitk::Index3D & mitkIndex ) { int retractionStartIndex( singleTract->Size() - 1 ); int retractionStepIndexSize( -1 ); int retractionTerminationIndex( 0 ); if( retractFront ) { retractionStartIndex = 0; retractionStepIndexSize = 1; retractionTerminationIndex = singleTract->Size() - 1; } int currentRetractionIndex = retractionStartIndex; bool keepRetracting( true ); mitk::Point3D currentPoint, nextPoint; std::vector< double > differenceVector; differenceVector.resize( singleTract->front().Size() ); while( keepRetracting && ( currentRetractionIndex != retractionTerminationIndex ) ) { // convert to segmentation coords mitk::Point3D currentPointFiberCoord, nextPointFiberCoord; for( int index = 0; index < singleTract->front().Size(); index++ ) { currentPointFiberCoord.SetElement( index, singleTract->GetElement( currentRetractionIndex ).GetElement( index ) ); nextPointFiberCoord.SetElement( index, singleTract->GetElement( currentRetractionIndex + retractionStepIndexSize ).GetElement( index ) ); } FiberToSegmentationCoords( currentPointFiberCoord, currentPoint ); FiberToSegmentationCoords( nextPointFiberCoord, nextPoint ); // calculate straight line for( int index = 0; index < singleTract->front().Size(); index++ ) { differenceVector[ index ] = nextPoint.GetElement( index ) - currentPoint.GetElement( index ); } // calculate length of direction vector double length( 0.0 ); double sum( 0.0 ); for( int index = 0; index < differenceVector.size() ; index++ ) { sum = sum + differenceVector[ index ] * differenceVector[ index ]; } length = std::sqrt( sum ); // retract mitk::Index3D tempIndex; int tempLabel( label ); for( int parameter( 0 ) ; parameter < length ; parameter++ ) { for( int index( 0 ); index < 3; index++ ) { tempIndex.SetElement( index, currentPoint.GetElement( index ) + ( 1.0 + parameter ) / ( 1.0 + length ) * differenceVector[ index ] ); } tempLabel = m_SegmentationItk->GetPixel( tempIndex ); if( !IsBackgroundLabel( tempLabel ) ) { // check whether result is within the search space { mitk::Point3D endPoint, foundPointSegmentation, foundPointFiber; for( int index = 0; index < singleTract->front().Size(); index++ ) { // this is in fiber (world) coordinates endPoint.SetElement( index, singleTract->GetElement( retractionStartIndex ).GetElement( index ) ); } for( int index( 0 ); index < 3; index++ ) { foundPointSegmentation.SetElement( index, currentPoint.GetElement( index ) + ( 1.0 + parameter ) / ( 1.0 + length ) * differenceVector[ index ] ); } SegmentationToFiberCoords( foundPointSegmentation, foundPointFiber ); std::vector< double > finalDistance; finalDistance.resize( singleTract->front().Size() ); for( int index = 0; index < singleTract->front().Size(); index++ ) { finalDistance[ index ] = foundPointFiber.GetElement( index ) - endPoint.GetElement( index ); } // calculate length of direction vector double finalLength( 0.0 ); double finalSum( 0.0 ); for( int index = 0; index < finalDistance.size() ; index++ ) { finalSum = finalSum + finalDistance[ index ] * finalDistance[ index ]; } finalLength = std::sqrt( finalSum ); if( finalLength > m_EndPointSearchRadius ) { // the found point was not within the search space return; } } label = tempLabel; mitkIndex = tempIndex; return; } // hit next point without finding brain matter currentRetractionIndex = currentRetractionIndex + retractionStepIndexSize; if( ( currentRetractionIndex < 1 ) || ( currentRetractionIndex > ( singleTract->Size() - 2 ) ) ) { keepRetracting = false; } } } } void mitk::ConnectomicsNetworkCreator::CalculateCenterOfMass() { const int dimensions = 3; int max = m_Segmentation->GetStatistics()->GetScalarValueMax(); int min = m_Segmentation->GetStatistics()->GetScalarValueMin(); int range = max - min +1; // each label owns a vector of coordinates std::vector< std::vector< std::vector< double> > > coordinatesPerLabelVector; coordinatesPerLabelVector.resize( range ); itk::ImageRegionIteratorWithIndex it_itkImage( m_SegmentationItk, m_SegmentationItk->GetLargestPossibleRegion() ); for( it_itkImage.GoToBegin(); !it_itkImage.IsAtEnd(); ++it_itkImage ) { std::vector< double > coordinates; coordinates.resize(dimensions); itk::Index< dimensions > index = it_itkImage.GetIndex(); for( int loop(0); loop < dimensions; loop++) { coordinates.at( loop ) = index.GetElement( loop ); } // add the coordinates to the corresponding label vector coordinatesPerLabelVector.at( it_itkImage.Value() - min ).push_back( coordinates ); } for(int currentIndex(0), currentLabel( min ); currentIndex < range; currentIndex++, currentLabel++ ) { std::vector< double > currentCoordinates; currentCoordinates.resize(dimensions); int numberOfPoints = coordinatesPerLabelVector.at( currentIndex ).size(); std::vector< double > sumCoords; sumCoords.resize( dimensions ); for( int loop(0); loop < numberOfPoints; loop++ ) { for( int loopDimension( 0 ); loopDimension < dimensions; loopDimension++ ) { sumCoords.at( loopDimension ) += coordinatesPerLabelVector.at( currentIndex ).at( loop ).at( loopDimension ); } } for( int loopDimension( 0 ); loopDimension < dimensions; loopDimension++ ) { currentCoordinates.at( loopDimension ) = sumCoords.at( loopDimension ) / numberOfPoints; } m_LabelsToCoordinatesMap.insert( std::pair< int, std::vector >( currentLabel, currentCoordinates ) ); } //can now use center of mass coordinates m_UseCoMCoordinates = true; } std::vector< double > mitk::ConnectomicsNetworkCreator::GetCenterOfMass( int label ) { // if label is not known, warn if( ! ( m_LabelsToCoordinatesMap.count( label ) > 0 ) ) { MITK_ERROR << "Label " << label << " not found. Could not return coordinates."; std::vector< double > nullVector; nullVector.resize(3); return nullVector; } //return associated coordinates return m_LabelsToCoordinatesMap.find( label )->second; } void mitk::ConnectomicsNetworkCreator::CreateNewNode( int label, mitk::Index3D index, bool useCoM ) { // Only create node if it does not exist yet if( ! ( m_LabelToNodePropertyMap.count( label ) > 0 ) ) { NetworkNode newNode; newNode.coordinates.resize( 3 ); if( !useCoM ) { for( unsigned int loop = 0; loop < newNode.coordinates.size() ; loop++ ) { newNode.coordinates[ loop ] = index[ loop ] ; } } else { std::vector labelCoords = GetCenterOfMass( label ); for( unsigned int loop = 0; loop < newNode.coordinates.size() ; loop++ ) { newNode.coordinates[ loop ] = labelCoords.at( loop ) ; } } newNode.label = LabelToString( label ); m_LabelToNodePropertyMap.insert( std::pair< ImageLabelType, NetworkNode >( label, newNode ) ); } } diff --git a/Modules/DiffusionImaging/Connectomics/Algorithms/mitkConnectomicsNetworkCreator.h b/Modules/DiffusionImaging/Connectomics/Algorithms/mitkConnectomicsNetworkCreator.h index 954725a234..0223e49bac 100644 --- a/Modules/DiffusionImaging/Connectomics/Algorithms/mitkConnectomicsNetworkCreator.h +++ b/Modules/DiffusionImaging/Connectomics/Algorithms/mitkConnectomicsNetworkCreator.h @@ -1,240 +1,248 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkConnectomicsNetworkCreator_h #define mitkConnectomicsNetworkCreator_h #include #include #include #include "mitkCommon.h" #include "mitkImage.h" #include "mitkFiberBundleX.h" #include "mitkConnectomicsNetwork.h" #include "ConnectomicsExports.h" namespace mitk { /** * \brief Creates connectomics networks from fibers and parcellation * * This class needs a parcellation image and a fiber image to be set. Then you can create * a connectomics network from the two, using different strategies. */ class Connectomics_EXPORT ConnectomicsNetworkCreator : public itk::Object { public: /** Enum for different ways to create the mapping from fibers to network */ enum MappingStrategy { EndElementPosition, PrecomputeAndDistance, JustEndPointVerticesNoLabel, EndElementPositionAvoidingWhiteMatter }; /** Standard class typedefs. */ /** Method for creation through the object factory. */ mitkClassMacro(ConnectomicsNetworkCreator, itk::Object); itkNewMacro(Self); /** Type for Images **/ typedef itk::Image ITKImageType; /** Types for the standardized Tract **/ typedef itk::Point PointType; typedef itk::VectorContainer TractType; typedef itk::VectorContainer< unsigned int, TractType::Pointer > TractContainerType; //init via smartpointer /** Types for Network **/ typedef mitk::ConnectomicsNetwork::VertexDescriptorType VertexType; typedef mitk::ConnectomicsNetwork::EdgeDescriptorType EdgeType; typedef mitk::ConnectomicsNetwork::NetworkNode NetworkNode; typedef std::pair< VertexType, VertexType > ConnectionType; /** Types for labels **/ typedef int ImageLabelType; typedef std::pair< ImageLabelType, ImageLabelType > ImageLabelPairType; /** Given a fiber bundle and a parcellation are set, this will create a network from both */ void CreateNetworkFromFibersAndSegmentation(); void SetFiberBundle(mitk::FiberBundleX::Pointer fiberBundle); void SetSegmentation(mitk::Image::Pointer segmentation); mitk::ConnectomicsNetwork::Pointer GetNetwork(); itkSetMacro(MappingStrategy, MappingStrategy); itkSetMacro(EndPointSearchRadius, double); + itkSetMacro(ZeroLabelInvalid, bool); /** \brief Calculate the locations of vertices * * Calculate the center of mass for each label and store the information. This will need a set parcellation image. * Unless this function is called the first location where a label is encountered will be used. After calling this function * the center of mass will be used instead. */ void CalculateCenterOfMass(); protected: //////////////////// Functions /////////////////////// ConnectomicsNetworkCreator(); ConnectomicsNetworkCreator( mitk::Image::Pointer segmentation, mitk::FiberBundleX::Pointer fiberBundle ); ~ConnectomicsNetworkCreator(); /** Add a connection to the network */ void AddConnectionToNetwork(ConnectionType newConnection); /** Determine if a label is already identified with a vertex, otherwise create a new one */ VertexType ReturnAssociatedVertexForLabel( ImageLabelType label ); /** Return the vertexes associated with a pair of labels */ ConnectionType ReturnAssociatedVertexPairForLabelPair( ImageLabelPairType labelpair ); /** Return the pair of labels which identify the areas connected by a single fiber */ ImageLabelPairType ReturnLabelForFiberTract( TractType::Pointer singleTract, MappingStrategy strategy ); /** Assign the additional information which should be part of the vertex */ void SupplyVertexWithInformation( ImageLabelType& label, VertexType& vertex ); /** Create a string from the label */ std::string LabelToString( ImageLabelType& label ); /** Check whether the label in question belongs to white matter according to the freesurfer table */ bool IsNonWhiteMatterLabel( int labelInQuestion ); /** Check whether the label in question belongs to background according to the freesurfer table */ bool IsBackgroundLabel( int labelInQuestion ); /** Extend a straight line through the given points and look for the first non white matter label It will try extend in the direction of the points in the vector so a vector {B,C} will result in extending from C in the direction C-B */ void LinearExtensionUntilGreyMatter( std::vector & indexVectorOfPointsToUse, TractType::Pointer singleTract, int & label, mitk::Index3D & mitkIndex ); /** Retract fiber until the first brain matter label is hit The bool parameter controls whether the front or the end is retracted */ void RetractionUntilBrainMatter( bool retractFront, TractType::Pointer singleTract, int & label, mitk::Index3D & mitkIndex ); /** \brief Get the location of the center of mass for a specific label * This can throw an exception if the label is not found. */ std::vector< double > GetCenterOfMass( int label ); /** Convert point to itk point */ itk::Point GetItkPoint(double point[3]); /** \brief Creates a new node * * A new node will be created, using the label and either the supplied coordinates * or the center of mass coordinates, depending on the supplied bool. */ void CreateNewNode( int label, mitk::Index3D, bool useIndex ); ///////// Mapping strategies ////////// /** Use the position of the end and starting element only to map to labels Map a fiber to a vertex by taking the value of the parcellation image at the same world coordinates as the last and first element of the tract.*/ ImageLabelPairType EndElementPositionLabel( TractType::Pointer singleTract ); /** Map by distance between elements and vertices depending on their volume First go through the parcellation and compute the coordinates of the future vertices. Assign a radius according on their volume. Then map an edge to a label by considering the nearest vertices and comparing the distance to them to their radii. */ ImageLabelPairType PrecomputeVertexLocationsBySegmentation( TractType::Pointer singleTract ); /** Use the position of the end and starting element only to map to labels Just take first and last position, no labelling, nothing */ ImageLabelPairType JustEndPointVerticesNoLabelTest( TractType::Pointer singleTract ); /** Use the position of the end and starting element unless it is in white matter, then search for nearby parcellation to map to labels Map a fiber to a vertex by taking the value of the parcellation image at the same world coordinates as the last and first element of the tract. If this happens to be white matter, then try to extend the fiber in a line and take the first non-white matter parcel, that is intersected. */ ImageLabelPairType EndElementPositionLabelAvoidingWhiteMatter( TractType::Pointer singleTract ); ///////// Conversions ////////// /** Convert fiber index to segmentation index coordinates */ void FiberToSegmentationCoords( mitk::Point3D& fiberCoord, mitk::Point3D& segCoord ); /** Convert segmentation index to fiber index coordinates */ void SegmentationToFiberCoords( mitk::Point3D& segCoord, mitk::Point3D& fiberCoord ); /////////////////////// Variables //////////////////////// mitk::FiberBundleX::Pointer m_FiberBundle; mitk::Image::Pointer m_Segmentation; ITKImageType::Pointer m_SegmentationItk; // the graph itself mitk::ConnectomicsNetwork::Pointer m_ConNetwork; // the id counter int idCounter; // the map mapping labels to vertices std::map< ImageLabelType, VertexType > m_LabelToVertexMap; // mapping labels to additional information std::map< ImageLabelType, NetworkNode > m_LabelToNodePropertyMap; // toggles whether edges between a node and itself can exist bool allowLoops; // toggles whether to use the center of mass coordinates bool m_UseCoMCoordinates; // stores the coordinates of labels std::map< int, std::vector< double> > m_LabelsToCoordinatesMap; // the straty to use for mapping MappingStrategy m_MappingStrategy; // search radius for finding a non white matter/background area. Should be in mm double m_EndPointSearchRadius; + // toggles whether a node with the label 0 may be present + bool m_ZeroLabelInvalid; + + // used internally to communicate a connection should not be added if the a problem + // is encountered while adding it + bool m_AbortConnection; + //////////////////////// IDs //////////////////////////// // These IDs are the freesurfer ids used in parcellation static const int freesurfer_Left_Cerebral_White_Matter = 2; static const int freesurfer_Left_Cerebellum_White_Matter = 7; static const int freesurfer_Left_Cerebellum_Cortex = 8; static const int freesurfer_Brain_Stem = 16; static const int freesurfer_Right_Cerebral_White_Matter = 41; static const int freesurfer_Right_Cerebellum_White_Matter = 46; static const int freesurfer_Right_Cerebellum_Cortex = 47; }; }// end namespace mitk #endif // _mitkConnectomicsNetworkCreator_H_INCLUDED