From 6a09951b1fc0a975f2cbf3263c1396fb153aedf2 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Tue, 24 Feb 2015 11:00:57 -0500 Subject: [PATCH] Test --- .../example_ImageAlgorithmDijkstra_03.cxx | 255 ++++++++++++------ ...example_ImageAlgorithm_Skeletonization.cxx | 1 + 2 files changed, 174 insertions(+), 82 deletions(-) diff --git a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx index 11fa961..f2d140b 100644 --- a/appli/examples/example_ImageAlgorithmDijkstra_03.cxx +++ b/appli/examples/example_ImageAlgorithmDijkstra_03.cxx @@ -112,7 +112,8 @@ protected: virtual void _AfterMainLoop( ) { - typedef itk::Image< unsigned char, TImage::ImageDimension > _TMarkImage; + typedef unsigned char _TMark; + typedef itk::Image< _TMark, TImage::ImageDimension > _TMarkImage; this->Superclass::_AfterMainLoop( ); if( this->m_Candidates.size( ) == 0 ) @@ -121,6 +122,7 @@ protected: const TImage* input = this->GetInput( ); TImage::SpacingType spacing = input->GetSpacing( ); + // Prepare an object to hold marks _TMarkImage::Pointer marks = _TMarkImage::New( ); marks->SetLargestPossibleRegion( input->GetLargestPossibleRegion( ) ); marks->SetRequestedRegion( input->GetRequestedRegion( ) ); @@ -129,93 +131,97 @@ protected: marks->SetSpacing( spacing ); marks->SetDirection( input->GetDirection( ) ); marks->Allocate( ); - marks->FillBuffer( 0 ); - - this->m_Axes = vtkSmartPointer< vtkPolyData >::New( ); - vtkSmartPointer< vtkPoints > points = - vtkSmartPointer< vtkPoints >::New( ); - vtkSmartPointer< vtkCellArray > lines = - vtkSmartPointer< vtkCellArray >::New( ); + marks->FillBuffer( _TMark( 0 ) ); + // Iterate over the candidates, starting fromt the candidates that + // are near thin branches _TCandidates::const_reverse_iterator cIt = this->m_Candidates.rbegin( ); - - for( unsigned long leo = 0; leo < 50 && cIt != this->m_Candidates.rend( ); ++cIt ) + for( int leo = 0; leo < 100 && cIt != this->m_Candidates.rend( ); ++cIt ) { - TVertex vIt = cIt->second; - if( marks->GetPixel( vIt ) != 0 ) + // If pixel hasn't been visited, continue + TVertex v = cIt->second; + if( marks->GetPixel( v ) != _TMark( 0 ) ) continue; - leo++; - std::cout << "Leaf: " << leo << " " << cIt->first << " " << vIt << std::endl; - bool start = true; - do - { - double dist = double( input->GetPixel( vIt ) ); + // Compute nearest start candidate TImage::SizeType radius; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - radius[ d ] = - ( unsigned long )( dist / spacing[ d ] ) + 1; - itk::NeighborhoodIterator< _TMarkImage > mIt( - radius, marks, marks->GetRequestedRegion( ) - ); - mIt.GoToBegin( ); - mIt.SetLocation( vIt ); - - TImage::SizeType size = mIt.GetSize( ); - unsigned int nN = 1; - for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) - nN *= size[ d ]; - for( unsigned int i = 0; i < nN; ++i ) - if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) ) - mIt.SetPixel( i, ( start )? 255: 100 ); + radius.Fill( 3 ); + itk::ConstNeighborhoodIterator< TImage > iIt( + radius, input, input->GetRequestedRegion( ) + ); + iIt.SetLocation( v ); + unsigned int nN = 1; + for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + nN *= iIt.GetSize( )[ d ]; + TImage::PixelType min_value = iIt.GetPixel( 0 ); + TVertex min_vertex = iIt.GetIndex( 0 ); + if( !( min_value > TImage::PixelType( 0 ) ) ) + min_value = std::numeric_limits< TImage::PixelType >::max( ); + for( unsigned int i = 1; i < nN; ++i ) + { + TImage::PixelType value = iIt.GetPixel( i ); + if( !( value > TImage::PixelType( 0 ) ) ) + value = std::numeric_limits< TImage::PixelType >::max( ); + if( value < min_value ) + { + min_value = value; + min_vertex = iIt.GetIndex( i ); - start = false; + } // fi - // Next vertex - vIt = this->_Parent( vIt ); + } // rof - } while( vIt != this->_Parent( vIt ) ); + if( min_value < std::numeric_limits< TImage::PixelType >::max( ) ) + { + if( marks->GetPixel( min_vertex ) != _TMark( 0 ) ) + continue; + leo++; + std::cout << "Leaf: " << leo << " " << min_value << " " << min_vertex << std::endl; - // Last vertex - // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl; + bool start = true; + do + { + double dist = double( 1.5 ) * std::sqrt( double( input->GetPixel( min_vertex ) ) ); + for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + radius[ d ] = + ( unsigned long )( dist / spacing[ d ] ) + 1; + itk::NeighborhoodIterator< _TMarkImage > mIt( + radius, marks, marks->GetRequestedRegion( ) + ); + mIt.SetLocation( min_vertex ); + nN = 1; + for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + nN *= mIt.GetSize( )[ d ]; + for( unsigned int i = 0; i < nN; ++i ) + if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) ) + { + mIt.SetPixel( i, ( start )? 255: 100 ); + start = false; + } + + /* + TImage::SizeType radius; + mIt.GoToBegin( ); + mIt.SetLocation( vIt ); + + TImage::SizeType size = mIt.GetSize( ); + unsigned int nN = 1; + for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + nN *= size[ d ]; + for( unsigned int i = 0; i < nN; ++i ) + if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) ) + mIt.SetPixel( i, ( start )? 255: 100 ); - /* TODO - TVertices pS; - TVertex parent = this->_Parent( vIt ); - while( parent != vIt ) - { - pS.push_front( vIt ); - vIt = parent; - parent = this->_Parent( vIt ); - - TImage::PointType pnt; - this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt ); - points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - if( points->GetNumberOfPoints( ) > 1 ) - { - lines->InsertNextCell( 2 ); - lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); - lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); - - } // fi - - } // elihw - pS.push_front( vIt ); - TImage::PointType pnt; - this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt ); - points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); - if( points->GetNumberOfPoints( ) > 1 ) - { - lines->InsertNextCell( 2 ); - lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); - lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); - - } // fi - - this->m_Axes->SetPoints( points ); - this->m_Axes->SetLines( lines ); - */ + start = false; + */ + // Next vertex in current path + min_vertex = this->_Parent( min_vertex ); + + } while( min_vertex != this->_Parent( min_vertex ) ); + } + else + marks->SetPixel( v, _TMark( 1 ) ); } // rof @@ -224,7 +230,92 @@ protected: w->SetInput( marks ); w->SetFileName( "marks.mhd" ); w->Update( ); + std::exit( 1 ); + /* + + this->m_Axes = vtkSmartPointer< vtkPolyData >::New( ); + vtkSmartPointer< vtkPoints > points = + vtkSmartPointer< vtkPoints >::New( ); + vtkSmartPointer< vtkCellArray > lines = + vtkSmartPointer< vtkCellArray >::New( ); + + { + + leo++; + std::cout << "Leaf: " << leo << " " << cIt->first << " " << vIt << std::endl; + bool start = true; + do + { + double dist = double( input->GetPixel( vIt ) ); + TImage::SizeType radius; + for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + radius[ d ] = + ( unsigned long )( dist / spacing[ d ] ) + 1; + itk::NeighborhoodIterator< _TMarkImage > mIt( + radius, marks, marks->GetRequestedRegion( ) + ); + mIt.GoToBegin( ); + mIt.SetLocation( vIt ); + + TImage::SizeType size = mIt.GetSize( ); + unsigned int nN = 1; + for( unsigned int d = 0; d < TImage::ImageDimension; ++d ) + nN *= size[ d ]; + for( unsigned int i = 0; i < nN; ++i ) + if( marks->GetRequestedRegion( ).IsInside( mIt.GetIndex( i ) ) ) + mIt.SetPixel( i, ( start )? 255: 100 ); + + start = false; + + // Next vertex + vIt = this->_Parent( vIt ); + + } while( vIt != this->_Parent( vIt ) ); + + // Last vertex + // std::cout << vIt << " " << this->GetSeed( 0 ) << std::endl; + */ + /* TODO + TVertices pS; + TVertex parent = this->_Parent( vIt ); + while( parent != vIt ) + { + pS.push_front( vIt ); + vIt = parent; + parent = this->_Parent( vIt ); + + TImage::PointType pnt; + this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt ); + points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + if( points->GetNumberOfPoints( ) > 1 ) + { + lines->InsertNextCell( 2 ); + lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); + lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); + + } // fi + + } // elihw + pS.push_front( vIt ); + TImage::PointType pnt; + this->GetInput( )->TransformIndexToPhysicalPoint( vIt, pnt ); + points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] ); + if( points->GetNumberOfPoints( ) > 1 ) + { + lines->InsertNextCell( 2 ); + lines->InsertCellPoint( points->GetNumberOfPoints( ) - 2 ); + lines->InsertCellPoint( points->GetNumberOfPoints( ) - 1 ); + + } // fi + + this->m_Axes->SetPoints( points ); + this->m_Axes->SetLines( lines ); + */ + + /* + } // rof + */ } virtual bool _UpdateNeigh( _TNode& nn, const _TNode& n ) @@ -246,13 +337,13 @@ protected: if( ret ) { TCost nc = this->_Cost( n.Vertex, n.Parent ); - TCost nc2 = nc * nc * nc * nc; - if( TCost( 0 ) < nc2 ) + if( TCost( 0 ) < nc ) { - n.Result = n.Cost / nc2; - this->GetOutput( )->SetPixel( n.Vertex, n.Result ); - this->m_Candidates.insert( _TCandidate( n.Result, n.Vertex ) ); - + TCost cc = n.Cost / nc; + this->m_Candidates.insert( _TCandidate( cc, n.Vertex ) ); + /* TODO: debug code + this->GetOutput( )->SetPixel( n.Vertex, cc ); + */ } // fi } // fi diff --git a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx index f2b8279..98023c5 100644 --- a/appli/examples/example_ImageAlgorithm_Skeletonization.cxx +++ b/appli/examples/example_ImageAlgorithm_Skeletonization.cxx @@ -124,6 +124,7 @@ int main( int argc, char* argv[] ) TDistanceFilter::Pointer distanceMap = TDistanceFilter::New( ); distanceMap->SetInput( segmentor->GetOutput( ) ); distanceMap->InputIsBinaryOn( ); + distanceMap->SquaredDistanceOn( ); start = std::clock( ); distanceMap->Update( ); end = std::clock( ); -- 2.45.1