APPLIS
example_Thinning
example_BinaryDistanceMap
+ example_HausdorffDistance
example_ImageAlgorithmRegionGrow_00
example_ImageAlgorithmDijkstra_00
example_ImageAlgorithmFastMarching_00
#include <itkMinimumMaximumImageCalculator.h>
#include <itkDanielssonDistanceMapImageFilter.h>
-/*
-#include <cmath>
-#include <deque>
-#include <limits>
-#include <map>
-
-#include <itkConstNeighborhoodIterator.h>
-#include <itkNeighborhoodIterator.h>
-#include <itkImageToVTKImageFilter.h>
-
-
-#include <vtkPoints.h>
-#include <vtkCellArray.h>
-#include <vtkFloatArray.h>
-#include <vtkPolyData.h>
-#include <vtkSmartPointer.h>
-
-#include <fpa/Image/DijkstraWithSphereBacktracking.h>
-#include <fpa/VTK/ImageMPR.h>
-#include <fpa/VTK/Image3DObserver.h>
-*/
-
// -------------------------------------------------------------------------
const unsigned int Dim = 3;
typedef unsigned char TPixel;
--- /dev/null
+#include <ctime>
+#include <iostream>
+#include <string>
+
+#include <itkImage.h>
+#include <itkImageFileReader.h>
+#include <itkImageFileWriter.h>
+#include <itkHausdorffDistanceImageFilter.h>
+
+// -------------------------------------------------------------------------
+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$
#include <ctime>
#include <deque>
#include <iostream>
+#include <iomanip>
#include <limits>
#include <map>
#include <string>
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[] )
{
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( );
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<vtkLookupTable> lut =
- vtkSmartPointer<vtkLookupTable>::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( );
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( );
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;
// 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
// 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
+
}
// -------------------------------------------------------------------------
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 );
}
this->m_Data->GetPointData( )->
GetScalars( )->InsertNextTuple1( back_id );
this->m_Data->Modified( );
-
return;
} // fi
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff