#include <fpa/VTK/ImageMPR.h>
#include <fpa/VTK/Image3DObserver.h>
-/*
- #include <itkMinimumMaximumImageCalculator.h>
- #include <itkInvertIntensityImageFilter.h>
- #include <itkDanielssonDistanceMapImageFilter.h>
- #include <itkImageFileWriter.h>
-
- #include <vtkCamera.h>
- #include <vtkImageActor.h>
- #include <vtkInteractorStyleImage.h>
- #include <vtkPointHandleRepresentation3D.h>
- #include <vtkProperty.h>
- #include <vtkRenderer.h>
- #include <vtkRenderWindow.h>
- #include <vtkRenderWindowInteractor.h>
- #include <vtkSeedRepresentation.h>
- #include <vtkSeedWidget.h>
- #include <vtkSmartPointer.h>
-
- #include <fpa/Image/Dijkstra.h>
-*/
+#include <vtkCellArray.h>
+#include <vtkImageMarchingCubes.h>
+#include <vtkPoints.h>
+#include <vtkPolyData.h>
+#include <vtkSmartPointer.h>
// -------------------------------------------------------------------------
const unsigned int Dim = 3;
template< class I >
void ReadImage( typename I::Pointer& image, const std::string& filename );
+template< class I >
+void SaveImage( const I* image, const std::string& filename );
+
template< class I, class O >
void DistanceMap(
const typename I::Pointer& input, typename O::Pointer& output
// -------------------------------------------------------------------------
int main( int argc, char* argv[] )
{
- if( argc < 2 )
+ if( argc < 5 )
{
std::cerr
<< "Usage: " << argv[ 0 ]
- << " input_image"
+ << " input_image distancemap output_costmap output_labels"
<< std::endl;
return( 1 );
} // fi
std::string input_image_fn = argv[ 1 ];
+ std::string distancemap_fn = argv[ 2 ];
+ std::string output_costmap_fn = argv[ 3 ];
+ std::string output_labels_fn = argv[ 4 ];
// Read image
TImage::Pointer input_image;
return( 1 );
} // yrt
+
TVTKInputImage::Pointer vtk_input_image = TVTKInputImage::New( );
vtk_input_image->SetInput( input_image );
vtk_input_image->Update( );
view.SetSize( 800, 800 );
view.SetImage( vtk_input_image->GetOutput( ) );
+ vtkSmartPointer< vtkImageMarchingCubes > mc =
+ vtkSmartPointer< vtkImageMarchingCubes >::New( );
+ mc->SetInputData( vtk_input_image->GetOutput( ) );
+ mc->SetValue( 0, 1e-1 );
+ mc->Update( );
+ view.AddPolyData( mc->GetOutput( ), 1, 1, 1, 0.4 );
+
// Allow some interaction and wait for, at least, one seed
view.Render( );
while( view.GetNumberOfSeeds( ) == 0 )
<< std::difftime( end, start )
<< " s." << std::endl;
- /* TODO
- // Save final total cost map
- itk::ImageFileWriter< TOutputImage >::Pointer output_image_writer =
- itk::ImageFileWriter< TOutputImage >::New( );
- output_image_writer->SetFileName( output_image_fn );
- output_image_writer->SetInput( filter->GetOutput( ) );
- try
+ // Minimum spanning tree
+ const TFilter::TMinimumSpanningTree* mst =
+ filter->GetMinimumSpanningTree( );
+
+ // Build branches from endpoints to seed
+ vtkSmartPointer< vtkPoints > endpoints_points =
+ vtkSmartPointer< vtkPoints >::New( );
+ vtkSmartPointer< vtkCellArray > endpoints_cells =
+ vtkSmartPointer< vtkCellArray >::New( );
+ vtkSmartPointer< vtkCellArray > endpoints_vertices =
+ vtkSmartPointer< vtkCellArray >::New( );
+
+ const TFilter::TUniqueVertices& endpoints = filter->GetEndPoints( );
+ TFilter::TUniqueVertices::const_iterator eIt = endpoints.begin( );
+ for( ; eIt != endpoints.end( ); ++eIt )
{
- output_image_writer->Update( );
- }
- catch( itk::ExceptionObject& err )
+ TFilter::TVertices path;
+ mst->GetPath( path, *eIt, filter->GetSeed( 0 ) );
+
+ std::cout << *eIt << std::endl;
+ TFilter::TVertices::const_iterator pIt = path.begin( );
+ for( ; pIt != path.end( ); ++pIt )
+ {
+ TImage::PointType pnt;
+ input_image->TransformIndexToPhysicalPoint( *pIt, pnt );
+ endpoints_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
+
+ if( pIt != path.begin( ) )
+ {
+ endpoints_cells->InsertNextCell( 2 );
+ endpoints_cells->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 2 );
+ endpoints_cells->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
+ }
+ else
+ {
+ std::cout << "ok" << std::endl;
+ endpoints_vertices->InsertNextCell( 1 );
+ endpoints_vertices->InsertCellPoint( endpoints_points->GetNumberOfPoints( ) - 1 );
+
+ } // fi
+
+
+
+
+ } // rof
+
+ } // rof
+ std::cout << endpoints.size( ) << std::endl;
+
+ vtkSmartPointer< vtkPolyData > endpoints_polydata =
+ vtkSmartPointer< vtkPolyData >::New( );
+ endpoints_polydata->SetPoints( endpoints_points );
+ endpoints_polydata->SetLines( endpoints_cells );
+ endpoints_polydata->SetVerts( endpoints_vertices );
+
+
+
+
+
+
+
+
+
+
+
+ // Bifurcations
+ vtkSmartPointer< vtkPoints > bifurcations_points =
+ vtkSmartPointer< vtkPoints >::New( );
+ vtkSmartPointer< vtkCellArray > bifurcations_vertices =
+ vtkSmartPointer< vtkCellArray >::New( );
+
+ const TFilter::TUniqueVertices& bifurcations = filter->GetBifurcationPoints( );
+ TFilter::TUniqueVertices::const_iterator bIt = bifurcations.begin( );
+ for( ; bIt != bifurcations.end( ); ++bIt )
{
- std::cerr
- << "Error while writing image to " << output_image_fn << ": "
- << err << std::endl;
- return( 1 );
+ std::cout << *bIt << std::endl;
- } // yrt
- */
+ TImage::PointType pnt;
+ input_image->TransformIndexToPhysicalPoint( *bIt, pnt );
+ bifurcations_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
+ bifurcations_vertices->InsertNextCell( 1 );
+ bifurcations_vertices->InsertCellPoint( bifurcations_points->GetNumberOfPoints( ) - 1 );
+
+ } // rof
+
+ vtkSmartPointer< vtkPolyData > bifurcations_polydata =
+ vtkSmartPointer< vtkPolyData >::New( );
+ bifurcations_polydata->SetPoints( bifurcations_points );
+ bifurcations_polydata->SetVerts( bifurcations_vertices );
+
+
+
+
+
+ // Build branches from endpoints to seed
+ vtkSmartPointer< vtkPoints > branches_points =
+ vtkSmartPointer< vtkPoints >::New( );
+ vtkSmartPointer< vtkCellArray > branches_cells =
+ vtkSmartPointer< vtkCellArray >::New( );
+
+ const TFilter::TBranches& branches = filter->GetBranches( );
+ TFilter::TBranches::const_iterator brIt = branches.begin( );
+ for( ; brIt != branches.end( ); ++brIt )
+ {
+ TFilter::TBranch::const_iterator brsIt = brIt->second.begin( );
+ for( ; brsIt != brIt->second.end( ); ++brsIt )
+ {
+ TImage::PointType pnt0, pnt1;
+ input_image->TransformIndexToPhysicalPoint( brIt->first, pnt0 );
+ input_image->TransformIndexToPhysicalPoint( brsIt->first, pnt1 );
+ branches_points->InsertNextPoint( pnt0[ 0 ], pnt0[ 1 ], pnt0[ 2 ] );
+ branches_points->InsertNextPoint( pnt1[ 0 ], pnt1[ 1 ], pnt1[ 2 ] );
+
+ branches_cells->InsertNextCell( 2 );
+ branches_cells->InsertCellPoint( branches_points->GetNumberOfPoints( ) - 2 );
+ branches_cells->InsertCellPoint( branches_points->GetNumberOfPoints( ) - 1 );
+
+ } // rof
+
+ } // rof
+
+ vtkSmartPointer< vtkPolyData > branches_polydata =
+ vtkSmartPointer< vtkPolyData >::New( );
+ branches_polydata->SetPoints( branches_points );
+ branches_polydata->SetLines( branches_cells );
+
+ view.AddPolyData( endpoints_polydata, 1, 0, 0, 1 );
+ view.AddPolyData( bifurcations_polydata, 0, 0, 1, 1 );
+ view.AddPolyData( branches_polydata, 0, 1, 0, 1 );
+ view.Render( );
+ view.Start( );
+
+ // Save final total cost map
+ SaveImage( filter->GetOutput( ), output_costmap_fn );
+ SaveImage( dmap.GetPointer( ), distancemap_fn );
+ SaveImage( filter->GetLabelImage( ), output_labels_fn );
return( 0 );
}
image->DisconnectPipeline( );
}
+// -------------------------------------------------------------------------
+template< class I >
+void SaveImage( const I* image, const std::string& filename )
+{
+ typename itk::ImageFileWriter< I >::Pointer writer =
+ itk::ImageFileWriter< I >::New( );
+ writer->SetInput( image );
+ writer->SetFileName( filename );
+ try
+ {
+ writer->Update( );
+ }
+ catch( itk::ExceptionObject& err )
+ {
+ std::cerr
+ << "Error saving \"" << filename << "\": " << err
+ << std::endl;
+
+ } // yrt
+}
+
// -------------------------------------------------------------------------
template< class I, class O >
void DistanceMap(
}
// eof - $RCSfile$
+
+
+
+
+
+
+
+
+
+