From d4efd9f4b13f746951059be3f9a2730c0425d1a5 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Thu, 19 Mar 2015 18:32:09 -0500 Subject: [PATCH] Submission to GRETSI --- .../example_ImageAlgorithmDijkstra_03.cxx | 191 ++++++++++-------- .../Image/DijkstraWithSphereBacktracking.h | 6 +- .../Image/DijkstraWithSphereBacktracking.hxx | 70 ++++++- 3 files changed, 172 insertions(+), 95 deletions(-) diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx index 5fe6b24..dd3c07d 100644 --- a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx +++ b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx @@ -1,6 +1,7 @@ #include #include #include +#include #include #include #include @@ -48,24 +49,24 @@ typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs; struct TBranch { double Length; - TImage::PointType V1; - TImage::PointType V2; + TImage::PointType::VectorType V1; + TImage::PointType::VectorType V2; bool operator<( const TBranch& other ) const { return( other.Length < this->Length ); } }; -typedef std::set< TBranch > TBranches; +typedef std::multiset< TBranch > TBranches; // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) { - if( argc < 6 ) + if( argc < 8 ) { std::cerr << "Usage: " << argv[ 0 ] - << " input_image x y z output_image" + << " input_image x y z output_image output_imagej output_vtk" << " visual_debug" << std::endl; return( 1 ); @@ -77,9 +78,11 @@ int main( int argc, char* argv[] ) seed_pnt[ 1 ] = std::atof( argv[ 3 ] ); seed_pnt[ 2 ] = std::atof( argv[ 4 ] ); std::string output_image_fn = argv[ 5 ]; + std::string output_imagej_fn = argv[ 6 ]; + std::string output_vtk_fn = argv[ 7 ]; bool visual_debug = false; - if( argc > 6 ) - visual_debug = ( std::atoi( argv[ 6 ] ) == 1 ); + if( argc > 8 ) + visual_debug = ( std::atoi( argv[ 8 ] ) == 1 ); // Read image TImageReader::Pointer input_image_reader = TImageReader::New( ); @@ -106,25 +109,28 @@ int main( int argc, char* argv[] ) view.SetSize( 800, 800 ); view.SetImage( vtk_image->GetOutput( ) ); - /* - vtkSmartPointer< vtkImageMarchingCubes > mc = - vtkSmartPointer< vtkImageMarchingCubes >::New( ); - mc->SetInputData( vtk_image->GetOutput( ) ); - mc->SetValue( 0, 1e-2 ); - mc->Update( ); - view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 ); + /* TODO + vtkSmartPointer< vtkImageMarchingCubes > mc = + vtkSmartPointer< vtkImageMarchingCubes >::New( ); + mc->SetInputData( vtk_image->GetOutput( ) ); + mc->SetValue( 0, 1e-2 ); + mc->Update( ); + view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 ); */ // Allow some interaction view.Render( ); - view.Start( ); + if( visual_debug ) + view.Start( ); TImage::IndexType seed_idx; input_image->TransformPhysicalPointToIndex( seed_pnt, seed_idx ); - std::cout << seed_idx << " " << seed_pnt << std::endl; - seed_idx[ 0 ] = 256; - seed_idx[ 1 ] = 313; - seed_idx[ 2 ] = 381; + /* TODO + std::cout << seed_idx << " " << seed_pnt << std::endl; + seed_idx[ 0 ] = 256; + seed_idx[ 1 ] = 313; + seed_idx[ 2 ] = 381; + */ // Extract paths TDijkstra::Pointer paths = TDijkstra::New( ); @@ -156,44 +162,55 @@ int main( int argc, char* argv[] ) vtkSmartPointer< vtkFloatArray > scalars = vtkSmartPointer< vtkFloatArray >::New( ); + TDijkstra::TMarkImage::Pointer new_marks = TDijkstra::TMarkImage::New( ); + new_marks->SetLargestPossibleRegion( input_image->GetLargestPossibleRegion( ) ); + new_marks->SetRequestedRegion( input_image->GetRequestedRegion( ) ); + new_marks->SetBufferedRegion( input_image->GetBufferedRegion( ) ); + new_marks->SetOrigin( input_image->GetOrigin( ) ); + new_marks->SetSpacing( input_image->GetSpacing( ) ); + new_marks->SetDirection( input_image->GetDirection( ) ); + new_marks->Allocate( ); + new_marks->FillBuffer( 0 ); + const TDijkstra::TMarkImage* marks = paths->GetOutputMarkImage( ); TDijkstra::TMark max_mark = paths->GetNumberOfBranches( ); const TDijkstra::TVertices& endpoints = paths->GetEndPoints( ); const TDijkstra::TVertices& bifurcations = paths->GetBifurcationPoints( ); - const TDijkstra::TTree& tree = paths->GetFinalTree( ); + const TDijkstra::TTree& tree = paths->GetFullTree( ); + const TDijkstra::TTree& reduced_tree = paths->GetReducedTree( ); TDijkstra::TVertices::const_iterator epIt = endpoints.begin( ); - TDijkstra::TMarkImage::Pointer new_marks = TDijkstra::TMarkImage::New( ); - new_marks->SetLargestPossibleRegion( marks->GetLargestPossibleRegion( ) ); - new_marks->SetRequestedRegion( marks->GetRequestedRegion( ) ); - new_marks->SetBufferedRegion( marks->GetBufferedRegion( ) ); - new_marks->SetOrigin( marks->GetOrigin( ) ); - new_marks->SetSpacing( marks->GetSpacing( ) ); - new_marks->SetDirection( marks->GetDirection( ) ); - new_marks->Allocate( ); - new_marks->FillBuffer( 0 ); - - std::cout << max_mark << std::endl; - std::cout << endpoints.size( ) << std::endl; - TBranches branches; - TBranch branch; - for( ; epIt != endpoints.end( ); ++epIt ) + TImage::PointType ori = input_image->GetOrigin( ); + + TDijkstra::TTree::const_iterator rIt = reduced_tree.begin( ); + for( ; rIt != reduced_tree.end( ); ++rIt ) { - TDijkstra::TVertex idx = *epIt; - TDijkstra::TTree::const_iterator tIt = tree.find( idx ); - TImage::PointType pnt, prev; - input_image->TransformIndexToPhysicalPoint( idx, pnt ); + TDijkstra::TTree::const_iterator tIt = tree.find( rIt->first ); + + // Prepare branch "a la ImageJ" + TBranch branch; + TImage::PointType p1, p2; + input_image->TransformIndexToPhysicalPoint( rIt->second.first, p1 ); + input_image->TransformIndexToPhysicalPoint( rIt->first, p2 ); + branch.V1 = p1 - ori; + branch.V2 = p2 - ori; + branch.Length = double( 0 ); + + double pd = + ( double( tIt->second.second ) - double( 1 ) ) / + ( double( max_mark ) - double( 1 ) ); + bool start = true; + TImage::PointType prev; do { + TDijkstra::TVertex idx = tIt->first; + new_marks->SetPixel( idx, 255 ); + TImage::PointType pnt; input_image->TransformIndexToPhysicalPoint( idx, pnt ); - branch.V2 = pnt; - points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - scalars->InsertNextTuple1( double( tIt->second.second ) ); - new_marks->SetPixel( idx, tIt->second.second ); - - if( idx != *epIt ) + scalars->InsertNextTuple1( pd ); + if( !start ) { cells->InsertNextCell( 2 ); cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); @@ -201,47 +218,56 @@ int main( int argc, char* argv[] ) branch.Length += prev.EuclideanDistanceTo( pnt ); } // fi + start = false; prev = pnt; - idx = tIt->second.first; - - if( std::find( bifurcations.begin( ), bifurcations.end( ), idx ) != bifurcations.end( ) ) - { - branches.insert( branch ); - prev = pnt; - branch.V1 = pnt; - branch.Length = double( 0 ); - } + tIt = tree.find( tIt->second.first ); - tIt = tree.find( idx ); + } while( tIt->first != rIt->second.first ); - } while( idx != tIt->second.first ); + // Last point + TDijkstra::TVertex idx = tIt->first; + new_marks->SetPixel( idx, 255 ); + TImage::PointType pnt; + input_image->TransformIndexToPhysicalPoint( idx, pnt ); + points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + scalars->InsertNextTuple1( pd ); + if( !start ) + { + cells->InsertNextCell( 2 ); + cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); + cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); + branch.Length += prev.EuclideanDistanceTo( pnt ); + } // fi + branches.insert( branch ); } // rof - /* TODO - int i = 1; - TImage::PointType ori = input_image->GetOrigin( ); - std::cout << ori << std::endl; - for( TBranches::const_iterator bIt = branches.begin( ); bIt != branches.end( ); ++bIt, ++i ) - { - TImage::PointType::VectorType v1 = bIt->V1 - ori; - TImage::PointType::VectorType v2 = bIt->V2 - ori; - std::cout - << std::setiosflags( std::ios::fixed) << std::setprecision( 3 ) - << i << "\t1.000\t" - << bIt->Length << "\t" - << v1[ 0 ] << "\t" - << v1[ 1 ] << "\t" - << v1[ 2 ] << "\t" - << v2[ 0 ] << "\t" - << v2[ 1 ] << "\t" - << v2[ 2 ] << "\t" - << ( v2 - v1 ).GetNorm( ) - << std::endl; - - } // rof - */ + std::ofstream bout( output_imagej_fn.c_str( ) ); + if( !bout ) + { + std::cerr << "Error: could not open " << output_imagej_fn << std::endl; + return( 1 ); + + } // fi + int i = 1; + for( TBranches::const_iterator bIt = branches.begin( ); bIt != branches.end( ); ++bIt, ++i ) + { + bout + << std::setiosflags( std::ios::fixed) << std::setprecision( 3 ) + << i << "\t1.000\t" + << bIt->Length << "\t" + << bIt->V1[ 0 ] << "\t" + << bIt->V1[ 1 ] << "\t" + << bIt->V1[ 2 ] << "\t" + << bIt->V2[ 0 ] << "\t" + << bIt->V2[ 1 ] << "\t" + << bIt->V2[ 2 ] << "\t" + << ( bIt->V2 - bIt->V1 ).GetNorm( ) + << std::endl; + + } // rof + bout.close( ); vtkSmartPointer< vtkPolyData > vtk_tree = vtkSmartPointer< vtkPolyData >::New( ); @@ -251,18 +277,19 @@ int main( int argc, char* argv[] ) view.AddPolyData( vtk_tree ); view.Render( ); - view.Start( ); + if( visual_debug ) + view.Start( ); vtkSmartPointer< vtkPolyDataWriter > writer = vtkSmartPointer< vtkPolyDataWriter >::New( ); writer->SetInputData( vtk_tree ); - writer->SetFileName( output_image_fn.c_str( ) ); + writer->SetFileName( output_vtk_fn.c_str( ) ); writer->Update( ); itk::ImageFileWriter< TDijkstra::TMarkImage >::Pointer marks_writer = itk::ImageFileWriter< TDijkstra::TMarkImage >::New( ); marks_writer->SetInput( new_marks ); - marks_writer->SetFileName( "marks.mhd" ); + marks_writer->SetFileName( output_image_fn ); marks_writer->Update( ); /* TODO diff --git a/lib/fpa/Image/DijkstraWithSphereBacktracking.h b/lib/fpa/Image/DijkstraWithSphereBacktracking.h index 0dc630e..137017a 100644 --- a/lib/fpa/Image/DijkstraWithSphereBacktracking.h +++ b/lib/fpa/Image/DijkstraWithSphereBacktracking.h @@ -52,7 +52,8 @@ namespace fpa itkNewMacro( Self ); itkTypeMacro( DijkstraWithSphereBacktracking, Dijkstra ); - itkGetConstMacro( FinalTree, TTree ); + itkGetConstMacro( FullTree, TTree ); + itkGetConstMacro( ReducedTree, TTree ); itkGetConstMacro( EndPoints, TVertices ); itkGetConstMacro( BifurcationPoints, TVertices ); itkGetConstMacro( NumberOfBranches, TMark ); @@ -78,7 +79,8 @@ namespace fpa protected: _TCandidates m_Candidates; - TTree m_FinalTree; + TTree m_FullTree; + TTree m_ReducedTree; TVertices m_BifurcationPoints; TVertices m_EndPoints; TMark m_NumberOfBranches; diff --git a/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx b/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx index 7d50735..69836a9 100644 --- a/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx +++ b/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx @@ -61,7 +61,7 @@ _BeforeMainLoop( ) for( unsigned int d = 1; d < I::ImageDimension; ++d ) max_spac = ( max_spac < double( spac[ d ] ) )? double( spac[ d ] ): max_spac; - max_spac *= double( 3 ); + max_spac *= double( 30 ); // Correct seeds for( unsigned int s = 0; s < this->m_Seeds.size( ); ++s ) @@ -99,7 +99,8 @@ _AfterMainLoop( ) { // Finish base algorithm this->Superclass::_AfterMainLoop( ); - this->m_FinalTree.clear( ); + this->m_FullTree.clear( ); + this->m_ReducedTree.clear( ); this->m_EndPoints.clear( ); this->m_BifurcationPoints.clear( ); if( this->m_Candidates.size( ) == 0 ) @@ -161,7 +162,7 @@ _AfterMainLoop( ) { if( start ) { - if( this->m_FinalTree.find( max_vertex ) == this->m_FinalTree.end( ) ) + if( this->m_FullTree.find( max_vertex ) == this->m_FullTree.end( ) ) { // Mark a sphere around current point as visited double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) ); @@ -173,7 +174,7 @@ _AfterMainLoop( ) // Next vertex in current path this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) ); - this->m_FinalTree[ max_vertex ] = + this->m_FullTree[ max_vertex ] = TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); // std::cout << "New: " << this->m_NumberOfBranches << std::endl; } @@ -183,7 +184,7 @@ _AfterMainLoop( ) this->m_BifurcationPoints.push_back( max_vertex ); this->m_NumberOfBranches++; - this->m_FinalTree[ max_vertex ] = + this->m_FullTree[ max_vertex ] = TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); // std::cout << "Bifurcation: " << this->m_NumberOfBranches << std::endl; @@ -211,7 +212,7 @@ _AfterMainLoop( ) // Next vertex in current path this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) ); - this->m_FinalTree[ max_vertex ] = + this->m_FullTree[ max_vertex ] = TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); change = true; // std::cout << "Change: " << this->m_NumberOfBranches << std::endl; @@ -221,7 +222,7 @@ _AfterMainLoop( ) // A bifurcation point has been reached! // TODO: this->m_BifurcationPoints.push_back( max_vertex ); this->m_NumberOfBranches++; - this->m_FinalTree[ max_vertex ] = + this->m_FullTree[ max_vertex ] = TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); // std::cout << "Change bifurcation: " << this->m_NumberOfBranches << std::endl; @@ -233,6 +234,7 @@ _AfterMainLoop( ) } while( max_vertex != this->_Parent( max_vertex ) ); if( start || change ) this->m_NumberOfBranches++; + this->m_FullTree[ max_vertex ] = TTreeNode( max_vertex, this->m_NumberOfBranches ); this->InvokeEvent( TEndBacktrackingEvent( ) ); @@ -240,8 +242,8 @@ _AfterMainLoop( ) std::map< TMark, unsigned long > histo; for( - typename TTree::iterator treeIt = this->m_FinalTree.begin( ); - treeIt != this->m_FinalTree.end( ); + typename TTree::iterator treeIt = this->m_FullTree.begin( ); + treeIt != this->m_FullTree.end( ); ++treeIt ) histo[ treeIt->second.second ]++; @@ -257,8 +259,8 @@ _AfterMainLoop( ) this->m_NumberOfBranches = changes.size( ); for( - typename TTree::iterator treeIt = this->m_FinalTree.begin( ); - treeIt != this->m_FinalTree.end( ); + typename TTree::iterator treeIt = this->m_FullTree.begin( ); + treeIt != this->m_FullTree.end( ); ++treeIt ) { @@ -267,6 +269,52 @@ _AfterMainLoop( ) } // fi + // Construct reduced tree + for( + typename TVertices::const_iterator eIt = this->m_EndPoints.begin( ); + eIt != this->m_EndPoints.end( ); + ++eIt + ) + { + typename TTree::const_iterator tIt = this->m_FullTree.find( *eIt ); + if( tIt != this->m_FullTree.end( ) ) + { + TMark id = tIt->second.second; + do + { + tIt = this->m_FullTree.find( tIt->second.first ); + if( tIt == this->m_FullTree.end( ) ) + break; + + } while( tIt->second.second == id ); + this->m_ReducedTree[ *eIt ] = TTreeNode( tIt->first, id ); + + } // fi + + } // rof + + for( + typename TVertices::const_iterator bIt = this->m_BifurcationPoints.begin( ); + bIt != this->m_BifurcationPoints.end( ); + ++bIt + ) + { + typename TTree::const_iterator tIt = this->m_FullTree.find( *bIt ); + if( tIt != this->m_FullTree.end( ) ) + { + TMark id = tIt->second.second; + do + { + tIt = this->m_FullTree.find( tIt->second.first ); + if( tIt == this->m_FullTree.end( ) ) + break; + + } while( tIt->second.second == id ); + this->m_ReducedTree[ *bIt ] = TTreeNode( tIt->first, id ); + + } // fi + + } // rof } // ------------------------------------------------------------------------- -- 2.47.1