appli/examples/example_CircleOfWillis.cxx

   1 #include <iostream>
   2 #include <fstream>
   3 #include <limits>
   4 #include <string>
   5
   6 #include <itkImage.h>
   7 #include <itkImageToVTKImageFilter.h>
   8 #include <itkImageFileReader.h>
   9
  10 #include <vtkCellArray.h>
  11 #include <vtkFloatArray.h>
  12 #include <vtkPointData.h>
  13 #include <vtkPoints.h>
  14 #include <vtkPolyData.h>
  15 #include <vtkSmartPointer.h>
  16
  17 #include <fpa/VTK/ImageMPR.h>
  18
  19 // -------------------------------------------------------------------------
  20 const   unsigned int Dimension = 3;
  21 typedef short        TPixel;
  22
  23 typedef itk::Image< TPixel, Dimension > TImage;
  24
  25 // -------------------------------------------------------------------------
  26 int main( int argc, char* argv[] )
  27 {
  28   if( argc < 3 )
  29   {
  30     std::cerr
  31       << "Usage: " << argv[ 0 ]
  32       << " input_image input_seeds"
  33       << std::endl;
  34       return( 1 );
  35
  36   } // fi
  37   std::string input_image_fn = argv[ 1 ];
  38   std::string input_seeds_fn = argv[ 2 ];
  39
  40   // Read image
  41   itk::ImageFileReader< TImage >::Pointer input_image_reader =
  42     itk::ImageFileReader< TImage >::New( );
  43   input_image_reader->SetFileName( input_image_fn );
  44   try
  45   {
  46     input_image_reader->Update( );
  47   }
  48   catch( itk::ExceptionObject& err )
  49   {
  50     std::cerr << "Error caugth: " << err << std::endl;
  51     return( 1 );
  52
  53   } // yrt
  54   TImage::ConstPointer input_image = input_image_reader->GetOutput( );
  55   itk::ImageToVTKImageFilter< TImage >::Pointer input_image_vtk =
  56     itk::ImageToVTKImageFilter< TImage >::New( );
  57   input_image_vtk->SetInput( input_image );
  58   input_image_vtk->Update( );
  59
  60   // Read seeds
  61   std::ifstream input_seeds_fs( input_seeds_fn.c_str( ) );
  62   if( !input_seeds_fs )
  63   {
  64     std::cerr
  65       << "Error opening file \"" << input_seeds_fn << "\"" << std::endl;
  66     return( 1 );
  67
  68   } // fi
  69
  70   // Read seed points
  71   std::string str_val;
  72   unsigned int number_of_points;
  73   input_seeds_fs >> str_val >> number_of_points;
  74   input_seeds_fs >> str_val; // X
  75   input_seeds_fs >> str_val; // Y
  76   input_seeds_fs >> str_val; // Z
  77   input_seeds_fs >> str_val; // value
  78   input_seeds_fs >> str_val; // Label
  79
  80   vtkSmartPointer< vtkPoints > input_seeds_points =
  81     vtkSmartPointer< vtkPoints >::New( );
  82   vtkSmartPointer< vtkCellArray > input_seeds_cells =
  83     vtkSmartPointer< vtkCellArray >::New( );
  84   vtkSmartPointer< vtkFloatArray > input_seeds_ids =
  85     vtkSmartPointer< vtkFloatArray >::New( );
  86
  87   input_seeds_ids->SetNumberOfComponents( 1 );
  88   input_seeds_ids->SetNumberOfTuples( number_of_points );
  89
  90   double min_value = std::numeric_limits< double >::max( );
  91   double max_value = double( 0 );
  92   for( unsigned int i = 0; i < number_of_points; ++i )
  93   {
  94     unsigned int x, y, z;
  95     double value;
  96     input_seeds_fs >> x >> y >> z >> value >> str_val;
  97
  98     min_value = ( value < min_value )? value: min_value;
  99     max_value = ( value > max_value )? value: max_value;
 100
 101     TImage::IndexType idx;
 102     idx[ 0 ] = x;
 103     idx[ 1 ] = y;
 104     idx[ 2 ] = z;
 105     TImage::PointType pnt;
 106     input_image->TransformIndexToPhysicalPoint( idx, pnt );
 107
 108     input_seeds_points->InsertNextPoint( pnt[ 0 ], pnt[ 1 ], pnt[ 2 ] );
 109     input_seeds_cells->InsertNextCell( 1 );
 110     input_seeds_cells->InsertCellPoint( i );
 111     input_seeds_ids->SetTuple1(
 112       i, double( i ) / double( number_of_points - 1 )
 113       );
 114
 115
 116   } // rof
 117   input_seeds_fs.close( );
 118
 119   vtkSmartPointer< vtkPolyData > input_seeds =
 120     vtkSmartPointer< vtkPolyData >::New( );
 121   input_seeds->SetPoints( input_seeds_points );
 122   input_seeds->SetVerts( input_seeds_cells );
 123   input_seeds->GetPointData( )->SetScalars( input_seeds_ids );
 124
 125   // Show input image and let some interaction
 126   fpa::VTK::ImageMPR view;
 127   view.SetBackground( 0.3, 0.2, 0.8 );
 128   view.SetSize( 600, 600 );
 129   view.SetImage( input_image_vtk->GetOutput( ) );
 130   view.SetWindowLevel(
 131     max_value - min_value,
 132     ( min_value + max_value ) / double( 2 )
 133     );
 134   view.AddPolyData( input_seeds );
 135   view.Render( );
 136   view.Start( );
 137
 138   return( 0 );
 139 }
 140
 141 // eof - $RCSfile$