#include <cmath>
#include <ctime>
#include <deque>
+#include <fstream>
#include <iostream>
#include <iomanip>
#include <limits>
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 );
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( );
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( );
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 );
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( );
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
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 )
{
// 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 )
{
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 ) ) );
// 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;
}
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;
// 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;
// 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;
} 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( ) );
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 ]++;
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
)
{
} // 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
}
// -------------------------------------------------------------------------
git://git.creatis.insa-lyon.fr / FrontAlgorithms.git / commitdiff