ANIMA  4.0
animaIsosurface.cxx
Go to the documentation of this file.
1 #include <tclap/CmdLine.h>
2 
3 #include <vtkContourFilter.h>
4 #include <vtkDecimatePro.h>
5 #include <vtkSmoothPolyDataFilter.h>
6 #include <vtkXMLPolyDataWriter.h>
7 #include <vtkPolyDataWriter.h>
8 #include <vtkPolyData.h>
9 
10 #include <itkImageToVTKImageFilter.h>
12 
13 int main(int argc, char **argv)
14 {
15  TCLAP::CmdLine cmd("Performs isosurface extraction followed by mesh decimation and surface smoothing. INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
16 
17  TCLAP::ValueArg<std::string> inArg("i","input","Input binary image",true,"","input binary image",cmd);
18  TCLAP::ValueArg<std::string> outArg("o","output","Output surface (.vtk or .vtp)",true,"","output surface",cmd);
19 
20  TCLAP::SwitchArg preSmoothArg("P", "pre-smooth", "Smoothes also before decimation", cmd, false);
21  TCLAP::ValueArg<double> thrArg("t","thr","Threshold for isosurface generation (default: 0.5)",false,0.5,"generation threshold",cmd);
22  TCLAP::ValueArg<double> smoothIntensityArg("s","smooth","Smoothing intensity (default: 0.2)",false,0.2,"smoothing intensity",cmd);
23  TCLAP::ValueArg<unsigned int> smoothIterArg("I","iter","Number of smoothing iterations (default: 50)",false,50,"smoothing iterations",cmd);
24  TCLAP::ValueArg<double> decimateArg("d","dec","Fraction of surface point to decimate (default: 0.75)",false,0.75,"deicimation fraction",cmd);
25 
26  try
27  {
28  cmd.parse(argc,argv);
29  }
30  catch (TCLAP::ArgException& e)
31  {
32  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
33  return EXIT_FAILURE;
34  }
35 
36  typedef itk::Image <double,3> ImageType;
37  ImageType::Pointer inputImage = anima::readImage <ImageType> (inArg.getValue());
38 
39  typedef itk::ImageToVTKImageFilter <ImageType> GlueFilterType;
40  GlueFilterType::Pointer glueFilter = GlueFilterType::New();
41 
42  glueFilter->SetInput(inputImage);
43  glueFilter->Update();
44 
45  vtkSmartPointer <vtkContourFilter> contourExtractor = vtkContourFilter::New();
46  contourExtractor->SetInputData(glueFilter->GetOutput());
47  contourExtractor->SetValue(0,thrArg.getValue());
48  contourExtractor->ComputeGradientsOff();
49  contourExtractor->ComputeNormalsOff();
50  contourExtractor->ComputeScalarsOff();
51 
52  vtkSmartPointer <vtkDecimatePro> decimateFilter = vtkDecimatePro::New();
53  if (preSmoothArg.isSet())
54  {
55  vtkSmartPointer <vtkSmoothPolyDataFilter> smoothFilterPre = vtkSmoothPolyDataFilter::New();
56  smoothFilterPre->SetInputConnection(contourExtractor->GetOutputPort());
57  smoothFilterPre->SetFeatureAngle(90);
58  smoothFilterPre->SetNumberOfIterations(smoothIterArg.getValue());
59  smoothFilterPre->SetRelaxationFactor(smoothIntensityArg.getValue());
60 
61  decimateFilter->SetInputConnection(smoothFilterPre->GetOutputPort());
62  }
63  else
64  decimateFilter->SetInputConnection(contourExtractor->GetOutputPort());
65 
66  decimateFilter->SetTargetReduction(decimateArg.getValue());
67  decimateFilter->PreserveTopologyOn();
68 
69  vtkSmartPointer <vtkSmoothPolyDataFilter> smoothFilter = vtkSmoothPolyDataFilter::New();
70  smoothFilter->SetInputConnection(decimateFilter->GetOutputPort());
71  smoothFilter->SetFeatureAngle(90);
72  smoothFilter->SetNumberOfIterations(smoothIterArg.getValue());
73  smoothFilter->SetRelaxationFactor(smoothIntensityArg.getValue());
74 
75  smoothFilter->Update();
76 
77  vtkSmartPointer <vtkPolyData> outputSurface = smoothFilter->GetOutput();
78  itk::ContinuousIndex <double, 3> currentIndex;
79  double currentPoint[3];
80  typedef ImageType::PointType PointType;
81  typedef ImageType::SpacingType SpacingType;
82  PointType outputPoint;
83  PointType origin = inputImage->GetOrigin();
84  SpacingType spacing = inputImage->GetSpacing();
85 
86  for (unsigned int i = 0;i < outputSurface->GetNumberOfPoints();++i)
87  {
88  outputSurface->GetPoint(i,currentPoint);
89  for (unsigned int j = 0;j < 3;++j)
90  currentIndex[j] = (currentPoint[j] - origin[j]) / spacing[j];
91 
92  inputImage->TransformContinuousIndexToPhysicalPoint(currentIndex,outputPoint);
93 
94  for (unsigned int j = 0;j < 3;++j)
95  currentPoint[j] = outputPoint[j];
96 
97  outputSurface->GetPoints()->SetPoint(i,currentPoint);
98  }
99 
100  if (outArg.getValue().find(".vtk") != std::string::npos)
101  {
102  vtkSmartPointer <vtkPolyDataWriter> legacyWriter = vtkPolyDataWriter::New();
103  legacyWriter->SetInputData(outputSurface);
104  legacyWriter->SetFileName(outArg.getValue().c_str());
105 
106  legacyWriter->Update();
107  }
108  else
109  {
110  vtkSmartPointer <vtkXMLPolyDataWriter> vtkWriter = vtkXMLPolyDataWriter::New();
111  vtkWriter->SetInputData(outputSurface);
112  vtkWriter->SetFileName(outArg.getValue().c_str());
113  vtkWriter->SetDataModeToBinary();
114  vtkWriter->EncodeAppendedDataOff();
115  vtkWriter->SetCompressorTypeToZLib();
116  vtkWriter->Update();
117  }
118 
119  return EXIT_SUCCESS;
120 }
int main(int argc, char **argv)