From 52c34a2aded47a2a3ce068887d9e4eeaf671568b Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Wed, 18 Mar 2015 18:22:57 -0500 Subject: [PATCH] Big bug stamped on wall --- appli/examples/CMakeLists.txt | 1 + appli/examples/example_BinaryDistanceMap.cxx | 22 --- appli/examples/example_HausdorffDistance.cxx | 59 +++++++ .../example_ImageAlgorithmDijkstra_03.cxx | 98 +++++++++-- .../Image/DijkstraWithSphereBacktracking.h | 6 +- .../Image/DijkstraWithSphereBacktracking.hxx | 165 +++++++++++------- lib/fpa/VTK/Image3DObserver.hxx | 1 - 7 files changed, 250 insertions(+), 102 deletions(-) create mode 100644 appli/examples/example_HausdorffDistance.cxx diff --git a/appli/examples/CMakeLists.txt b/appli/examples/CMakeLists.txt index 93d3c1d..254ef0b 100644 --- a/appli/examples/CMakeLists.txt +++ b/appli/examples/CMakeLists.txt @@ -3,6 +3,7 @@ IF(BUILD_EXAMPLES) APPLIS example_Thinning example_BinaryDistanceMap + example_HausdorffDistance example_ImageAlgorithmRegionGrow_00 example_ImageAlgorithmDijkstra_00 example_ImageAlgorithmFastMarching_00 diff --git a/appli/examples/example_BinaryDistanceMap.cxx b/appli/examples/example_BinaryDistanceMap.cxx index 5fac440..7628ebf 100644 --- a/appli/examples/example_BinaryDistanceMap.cxx +++ b/appli/examples/example_BinaryDistanceMap.cxx @@ -9,28 +9,6 @@ #include #include -/* -#include -#include -#include -#include - -#include -#include -#include - - -#include -#include -#include -#include -#include - -#include -#include -#include -*/ - // ------------------------------------------------------------------------- const unsigned int Dim = 3; typedef unsigned char TPixel; diff --git a/appli/examples/example_HausdorffDistance.cxx b/appli/examples/example_HausdorffDistance.cxx new file mode 100644 index 0000000..ebcede0 --- /dev/null +++ b/appli/examples/example_HausdorffDistance.cxx @@ -0,0 +1,59 @@ +#include +#include +#include + +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +const unsigned int Dim = 3; +typedef unsigned short TPixel; +typedef float TScalar; +typedef itk::Image< TPixel, Dim > TImage; +typedef itk::ImageFileReader< TImage > TImageReader; +typedef itk::HausdorffDistanceImageFilter< TImage, TImage > TDistance; + +// ------------------------------------------------------------------------- +int main( int argc, char* argv[] ) +{ + if( argc < 3 ) + { + std::cerr + << "Usage: " << argv[ 0 ] + << " first_image second_image" + << std::endl; + return( 1 ); + + } // fi + std::string first_image_fn = argv[ 1 ]; + std::string second_image_fn = argv[ 2 ]; + + // Read image + TImageReader::Pointer first_image_reader = TImageReader::New( ); + TImageReader::Pointer second_image_reader = TImageReader::New( ); + first_image_reader->SetFileName( first_image_fn ); + second_image_reader->SetFileName( second_image_fn ); + try + { + first_image_reader->Update( ); + second_image_reader->Update( ); + } + catch( itk::ExceptionObject& err ) + { + std::cerr << "Error caught: " << err << std::endl; + return( 1 ); + + } // yrt + + TDistance::Pointer distance = TDistance::New( ); + distance->SetInput1( first_image_reader->GetOutput( ) ); + distance->SetInput2( second_image_reader->GetOutput( ) ); + distance->Update( ); + std::cout << "Hausdorff distance = " << distance->GetHausdorffDistance( ) << std::endl; + + return( 0 ); +} + +// eof - $RCSfile$ diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx index 613378b..5fe6b24 100644 --- a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx +++ b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx @@ -2,6 +2,7 @@ #include #include #include +#include #include #include #include @@ -44,6 +45,19 @@ typedef fpa::Image::DijkstraWithSphereBacktracking< TImage, TScalar > TDijkstra; typedef fpa::VTK::ImageMPR TMPR; typedef fpa::VTK::Image3DObserver< TDijkstra, vtkRenderWindow > TDijkstraObs; +struct TBranch +{ + double Length; + TImage::PointType V1; + TImage::PointType V2; + + bool operator<( const TBranch& other ) const + { + return( other.Length < this->Length ); + } +}; +typedef std::set< TBranch > TBranches; + // ------------------------------------------------------------------------- int main( int argc, char* argv[] ) { @@ -92,12 +106,14 @@ 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 ); + */ // Allow some interaction view.Render( ); @@ -143,47 +159,97 @@ int main( int argc, char* argv[] ) 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( ); TDijkstra::TVertices::const_iterator epIt = endpoints.begin( ); - for( unsigned int epId = 0; epIt != endpoints.end( ); ++epIt, ++epId ) + + 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 ) { TDijkstra::TVertex idx = *epIt; + TDijkstra::TTree::const_iterator tIt = tree.find( idx ); + TImage::PointType pnt, prev; + input_image->TransformIndexToPhysicalPoint( idx, pnt ); do { - TImage::PointType pnt; input_image->TransformIndexToPhysicalPoint( idx, pnt ); - - TDijkstra::TMark mark = marks->GetPixel( idx ); - double pd = double( mark ); + branch.V2 = pnt; points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - scalars->InsertNextTuple1( pd ); + scalars->InsertNextTuple1( double( tIt->second.second ) ); + new_marks->SetPixel( idx, tIt->second.second ); + if( idx != *epIt ) { cells->InsertNextCell( 2 ); cells->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); cells->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); + branch.Length += prev.EuclideanDistanceTo( pnt ); } // fi - idx = tree.find( idx )->second; + 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( idx ); + + } while( idx != tIt->second.first ); - } while( idx != tree.find( idx )->second ); } // 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 + */ + vtkSmartPointer< vtkPolyData > vtk_tree = vtkSmartPointer< vtkPolyData >::New( ); vtk_tree->SetPoints( points ); vtk_tree->SetLines( cells ); vtk_tree->GetPointData( )->SetScalars( scalars ); - vtkSmartPointer lut = - vtkSmartPointer::New( ); - lut->SetNumberOfTableValues( max_mark + 1 ); - lut->SetTableRange( 0, max_mark ); - lut->Build( ); - - view.AddPolyData( vtk_tree, lut ); + view.AddPolyData( vtk_tree ); view.Render( ); view.Start( ); @@ -195,7 +261,7 @@ int main( int argc, char* argv[] ) itk::ImageFileWriter< TDijkstra::TMarkImage >::Pointer marks_writer = itk::ImageFileWriter< TDijkstra::TMarkImage >::New( ); - marks_writer->SetInput( marks ); + marks_writer->SetInput( new_marks ); marks_writer->SetFileName( "marks.mhd" ); marks_writer->Update( ); diff --git a/lib/fpa/Image/DijkstraWithSphereBacktracking.h b/lib/fpa/Image/DijkstraWithSphereBacktracking.h index febfcc9..0dc630e 100644 --- a/lib/fpa/Image/DijkstraWithSphereBacktracking.h +++ b/lib/fpa/Image/DijkstraWithSphereBacktracking.h @@ -28,12 +28,12 @@ namespace fpa typedef typename Superclass::InputImageType TImage; typedef std::deque< TVertex > TVertices; - typedef typename Superclass::TTraits::TVertexCmp TVertexCmp; - typedef std::map< TVertex, TVertex, TVertexCmp > TTree; - typedef unsigned short TMark; typedef itk::Image< TMark, I::ImageDimension > TMarkImage; + typedef typename Superclass::TTraits::TVertexCmp TVertexCmp; + typedef std::pair< TVertex, TMark > TTreeNode; + typedef std::map< TVertex, TTreeNode, TVertexCmp > TTree; typedef typename Superclass::TEndEvent TEndEvent; typedef typename Superclass::TBacktrackingEvent TBacktrackingEvent; diff --git a/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx b/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx index 346103b..7d50735 100644 --- a/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx +++ b/lib/fpa/Image/DijkstraWithSphereBacktracking.hxx @@ -123,7 +123,7 @@ _AfterMainLoop( ) // are near thin branches typename _TCandidates::const_reverse_iterator cIt = this->m_Candidates.rbegin( ); - this->m_NumberOfBranches = 0; + this->m_NumberOfBranches = 1; for( ; cIt != this->m_Candidates.rend( ); ++cIt ) { // If pixel hasn't been visited, continue @@ -152,81 +152,121 @@ _AfterMainLoop( ) // Keep endpoint if( marks->GetPixel( max_vertex ) != 0 ) continue; - this->m_NumberOfBranches++; this->m_EndPoints.push_back( max_vertex ); // Construct path (at least the part that hasn't been iterated) + bool start = true; + bool change = false; do { - if( - std::find( - this->m_BifurcationPoints.begin( ), - this->m_BifurcationPoints.end( ), - max_vertex - ) == this->m_BifurcationPoints.end( ) - ) + if( start ) { - // Mark a sphere around current point as visited - double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) ); - region = this->_Region( max_vertex, dist * double( 1.1 ) ); - itk::ImageRegionIteratorWithIndex< TMarkImage > - mIt( marks, region ); - for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) - mIt.Set( this->m_NumberOfBranches ); - - // Next vertex in current path - this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) ); - this->m_FinalTree[ max_vertex ] = this->_Parent( max_vertex ); + if( this->m_FinalTree.find( max_vertex ) == this->m_FinalTree.end( ) ) + { + // Mark a sphere around current point as visited + double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) ); + region = this->_Region( max_vertex, dist * double( 1.5 ) ); + itk::ImageRegionIteratorWithIndex< TMarkImage > + mIt( marks, region ); + for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) + mIt.Set( this->m_NumberOfBranches ); + + // Next vertex in current path + this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) ); + this->m_FinalTree[ max_vertex ] = + TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); + // std::cout << "New: " << this->m_NumberOfBranches << std::endl; + } + else + { + // A bifurcation point has been reached! + this->m_BifurcationPoints.push_back( max_vertex ); + this->m_NumberOfBranches++; + + this->m_FinalTree[ max_vertex ] = + TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); + // std::cout << "Bifurcation: " << this->m_NumberOfBranches << std::endl; + + start = false; + + } // fi } else { - // A bifurcation point has been reached! - this->m_BifurcationPoints.push_back( max_vertex ); - this->m_NumberOfBranches++; + if( + std::find( + this->m_BifurcationPoints.begin( ), + this->m_BifurcationPoints.end( ), + max_vertex + ) == this->m_BifurcationPoints.end( ) + ) + { + // Mark a sphere around current point as visited + double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) ); + region = this->_Region( max_vertex, dist * double( 1.5 ) ); + itk::ImageRegionIteratorWithIndex< TMarkImage > + mIt( marks, region ); + for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) + mIt.Set( this->m_NumberOfBranches ); + + // Next vertex in current path + this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) ); + this->m_FinalTree[ max_vertex ] = + TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); + change = true; + // std::cout << "Change: " << this->m_NumberOfBranches << std::endl; + } + else + { + // A bifurcation point has been reached! + // TODO: this->m_BifurcationPoints.push_back( max_vertex ); + this->m_NumberOfBranches++; + this->m_FinalTree[ max_vertex ] = + TTreeNode( this->_Parent( max_vertex ), this->m_NumberOfBranches ); + // std::cout << "Change bifurcation: " << this->m_NumberOfBranches << std::endl; + + } // fi } // fi max_vertex = this->_Parent( max_vertex ); } while( max_vertex != this->_Parent( max_vertex ) ); + if( start || change ) + this->m_NumberOfBranches++; - /* TODO - bool terminate = false; - do - { - if( this->m_FinalTree.find( max_vertex ) == this->m_FinalTree.end( ) ) - { - // Mark a sphere around current point as visited - double dist = std::sqrt( double( input->GetPixel( max_vertex ) ) ); - region = this->_Region( max_vertex, dist * double( 1.25 ) ); - itk::ImageRegionIteratorWithIndex< TMarkImage > - mIt( marks, region ); - for( mIt.GoToBegin( ); !mIt.IsAtEnd( ); ++mIt ) - mIt.Set( true ); - - // Next vertex in current path - this->InvokeEvent( TBacktrackingEvent( max_vertex, this->m_NumberOfBranches ) ); - this->m_FinalTree[ max_vertex ] = this->_Parent( max_vertex ); - } - else - { - // A bifurcation point has been reached! - this->m_BifurcationPoints.push_back( max_vertex ); - terminate = true; - - } // fi - max_vertex = this->_Parent( max_vertex ); - - } while( max_vertex != this->_Parent( max_vertex ) && !terminate ); - - if( !terminate ) - { - this->m_FinalTree[ max_vertex ] = max_vertex; - this->InvokeEvent( TEndBacktrackingEvent( this->m_NumberOfBranches ) ); - - } // fi - */ + this->InvokeEvent( TEndBacktrackingEvent( ) ); } // rof + + std::map< TMark, unsigned long > histo; + for( + typename TTree::iterator treeIt = this->m_FinalTree.begin( ); + treeIt != this->m_FinalTree.end( ); + ++treeIt + ) + histo[ treeIt->second.second ]++; + + std::map< TMark, TMark > changes; + TMark last_change = 1; + for( TMark i = 1; i <= this->m_NumberOfBranches; ++i ) + { + if( histo[ i ] != 0 ) + changes[ i ] = last_change++; + + } // rof + this->m_NumberOfBranches = changes.size( ); + + for( + typename TTree::iterator treeIt = this->m_FinalTree.begin( ); + treeIt != this->m_FinalTree.end( ); + ++treeIt + ) + { + TMark old = treeIt->second.second; + treeIt->second.second = changes[ old ]; + + } // fi + } // ------------------------------------------------------------------------- @@ -237,8 +277,13 @@ _UpdateNeigh( _TNode& nn, const _TNode& n ) C nc = this->_Cost( nn.Vertex, n.Vertex ); if( TCost( 0 ) < nc ) { + typename I::PointType pnn, pn; + this->GetInput( )->TransformIndexToPhysicalPoint( nn.Vertex, pnn ); + this->GetInput( )->TransformIndexToPhysicalPoint( n.Vertex, pn ); + + nc += TCost( 1 ); - nn.Cost = n.Cost + ( TCost( 1 ) / std::pow( nc, 4 ) ); + nn.Cost = n.Cost + ( TCost( pnn.EuclideanDistanceTo( pn ) ) / std::pow( nc, 4 ) ); nn.Result = nn.Cost; return( true ); } diff --git a/lib/fpa/VTK/Image3DObserver.hxx b/lib/fpa/VTK/Image3DObserver.hxx index e9f9770..46ef743 100644 --- a/lib/fpa/VTK/Image3DObserver.hxx +++ b/lib/fpa/VTK/Image3DObserver.hxx @@ -176,7 +176,6 @@ Execute( const itk::Object* c, const itk::EventObject& e ) this->m_Data->GetPointData( )-> GetScalars( )->InsertNextTuple1( back_id ); this->m_Data->Modified( ); - return; } // fi -- 2.47.1