ANIMA  4.0
animaFibersCounter.cxx
Go to the documentation of this file.
1 #include <tclap/CmdLine.h>
2 
4 #include <animaShapesReader.h>
5 
6 #include <vtkSmartPointer.h>
7 #include <vtkPoints.h>
8 #include <vtkPointData.h>
9 #include <vtkPolyData.h>
10 #include <vtkCell.h>
11 
12 #include <itkCastImageFilter.h>
13 
14 int main(int argc, char **argv)
15 {
16  TCLAP::CmdLine cmd("Filters fibers from a vtp file using a label image and specifying with several -t and -f which labels should be touched or are forbidden for each fiber. INRIA / IRISA - VisAGeS/Empenn Team", ' ',ANIMA_VERSION);
17 
18  TCLAP::ValueArg<std::string> inArg("i","input","input tracks file",true,"","input tracks",cmd);
19  TCLAP::ValueArg<std::string> outArg("o","output","output mask image",true,"","output mask image",cmd);
20  TCLAP::ValueArg<std::string> geomArg("g","geometry","Geometry image",true,"","geometry image",cmd);
21 
22  TCLAP::SwitchArg proportionArg("P","proportion","Output proportion of fibers going through each pixel",cmd,false);
23 
24  try
25  {
26  cmd.parse(argc,argv);
27  }
28  catch (TCLAP::ArgException& e)
29  {
30  std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
31  return EXIT_FAILURE;
32  }
33 
34  typedef itk::Image <double, 3> OutputImageType;
35  OutputImageType::Pointer outputImage = anima::readImage <OutputImageType> (geomArg.getValue());
36  outputImage->FillBuffer(0.0);
37 
38  anima::ShapesReader trackReader;
39  trackReader.SetFileName(inArg.getValue());
40  trackReader.Update();
41 
42  vtkSmartPointer <vtkPolyData> tracks = trackReader.GetOutput();
43 
44  vtkIdType nbCells = tracks->GetNumberOfCells();
45  double incrementFactor = 1.0;
46  if (proportionArg.isSet())
47  incrementFactor /= nbCells;
48 
49  double ptVals[3];
50  OutputImageType::IndexType currentIndex;
51  OutputImageType::PointType currentPoint;
52  OutputImageType::RegionType region = outputImage->GetLargestPossibleRegion();
53 
54  // Explores individual fibers
55  for (int j = 0;j < nbCells;++j)
56  {
57  vtkCell *cell = tracks->GetCell(j);
58  vtkPoints *cellPts = cell->GetPoints();
59  vtkIdType nbPts = cellPts->GetNumberOfPoints();
60 
61  // Explores points in fibers
62  for (int i = 0;i < nbPts;++i)
63  {
64  cellPts->GetPoint(i,ptVals);
65 
66  for (unsigned int k = 0;k < 3;++k)
67  currentPoint[k] = ptVals[k];
68 
69  outputImage->TransformPhysicalPointToIndex(currentPoint,currentIndex);
70  bool insideBuffer = true;
71  for (unsigned int k = 0;k < 3;++k)
72  {
73  if ((currentIndex[k] < 0)||(currentIndex[k] >= region.GetSize(k)))
74  {
75  insideBuffer = false;
76  break;
77  }
78  }
79 
80  if (!insideBuffer)
81  continue;
82 
83  double countIndex = outputImage->GetPixel(currentIndex) + incrementFactor;
84  outputImage->SetPixel(currentIndex, countIndex);
85  }
86  }
87 
88  if (proportionArg.isSet())
89  anima::writeImage <OutputImageType> (outArg.getValue(),outputImage);
90  else
91  {
92  using MaskImageType = itk::Image <unsigned int, 3>;
93  using CastFilterType = itk::CastImageFilter <OutputImageType, MaskImageType>;
94 
95  CastFilterType::Pointer castFilter = CastFilterType::New();
96  castFilter->SetInput(outputImage);
97  castFilter->Update();
98 
99  anima::writeImage <MaskImageType> (outArg.getValue(), castFilter->GetOutput());
100  }
101 
102  return EXIT_SUCCESS;
103 }
int main(int argc, char **argv)
void SetFileName(std::string &name)
vtkPolyData * GetOutput()