ANIMA  4.0
animaGraph3DFilter.hxx
Go to the documentation of this file.
1 #pragma once
2 
3 #include "animaGraph3DFilter.h"
4 
5 namespace anima
6 {
7 
8 template <typename TInput, typename TOutput>
9 void Graph3DFilter<TInput, TOutput>::SetInputImage(unsigned int i, const TInput* image)
10 {
11  if (i == 0)
12  this->SetNthInput(0, const_cast<TInput*>(image));
13  else
14  {
15  this->SetNthInput(m_NbInputs, const_cast<TInput*>(image));
16 
17  switch(i)
18  {
19  case 1:
20  m_IndexImage2 = m_NbInputs;
21  break;
22 
23  case 2:
24  m_IndexImage3 = m_NbInputs;
25  break;
26 
27  case 3:
28  m_IndexImage4 = m_NbInputs;
29  break;
30 
31  case 4:
32  m_IndexImage5 = m_NbInputs;
33  break;
34 
35  default:
36  itkExceptionMacro("Only a maximum of 5 input images are accepted");
37  break;
38  }
39 
40  m_NbInputs++;
41  }
42 }
43 
44 template <typename TInput, typename TOutput>
46 {
47  this->SetNthInput(1, const_cast<TSeedProba*>(proba));
48 }
49 
50 template <typename TInput, typename TOutput>
52 {
53  this->SetNthInput(2, const_cast<TSeedProba*>(proba));
54 }
55 
56 template <typename TInput, typename TOutput>
58 {
59  this->SetNthInput(3, const_cast<TMask*>(mask));
60 }
61 
62 template <typename TInput, typename TOutput>
63 typename TInput::ConstPointer Graph3DFilter<TInput, TOutput>::GetInputImage1()
64 {
65  return static_cast< const TInput * >
66  ( this->itk::ProcessObject::GetInput(0) );
67 }
68 
69 template <typename TInput, typename TOutput>
70 itk::Image <double,3>::ConstPointer Graph3DFilter<TInput, TOutput>::GetInputSeedProbaSources()
71 {
72  return static_cast< const TSeedProba * >
73  ( this->itk::ProcessObject::GetInput(1) );
74 }
75 
76 template <typename TInput, typename TOutput>
77 itk::Image <double,3>::ConstPointer Graph3DFilter<TInput, TOutput>::GetInputSeedProbaSinks()
78 {
79  return static_cast< const TSeedProba * >
80  ( this->itk::ProcessObject::GetInput(2) );
81 }
82 
83 template <typename TInput, typename TOutput>
84 itk::Image <unsigned char,3>::ConstPointer Graph3DFilter<TInput, TOutput>::GetMask()
85 {
86  return static_cast< const TMask * >
87  ( this->itk::ProcessObject::GetInput(3) );
88 }
89 
90 template <typename TInput, typename TOutput>
91 typename TInput::ConstPointer Graph3DFilter<TInput, TOutput>::GetInputImage2()
92 {
93  return static_cast< const TInput * >
94  ( this->itk::ProcessObject::GetInput(m_IndexImage2) );
95 }
96 
97 template <typename TInput, typename TOutput>
98 typename TInput::ConstPointer Graph3DFilter<TInput, TOutput>::GetInputImage3()
99 {
100  return static_cast< const TInput * >
101  ( this->itk::ProcessObject::GetInput(m_IndexImage3) );
102 }
103 
104 template <typename TInput, typename TOutput>
105 typename TInput::ConstPointer Graph3DFilter<TInput, TOutput>::GetInputImage4()
106 {
107  return static_cast< const TInput * >
108  ( this->itk::ProcessObject::GetInput(m_IndexImage4) );
109 }
110 
111 template <typename TInput, typename TOutput>
112 typename TInput::ConstPointer Graph3DFilter<TInput, TOutput>::GetInputImage5()
113 {
114  return static_cast< const TInput * >
115  ( this->itk::ProcessObject::GetInput(m_IndexImage5) );
116 }
117 
118 template <typename TInputImage, typename TOutput>
119 itk::DataObject::Pointer Graph3DFilter <TInputImage, TOutput>::MakeOutput(unsigned int idx)
120 {
121  itk::DataObject::Pointer output;
122 
123  switch ( idx )
124  {
125  case 0:
126  output = ( TOutput::New() ).GetPointer();
127  break;
128  case 1:
129  output = ( TOutput::New() ).GetPointer();
130  break;
131  default:
132  std::cerr << "No output " << idx << std::endl;
133  output = NULL;
134  break;
135  }
136  return output.GetPointer();
137 }
138 
139 template <typename TInputImage, typename TOutput>
141 {
142  return dynamic_cast< TOutput* >( this->itk::ProcessObject::GetOutput(0) );
143 }
144 template <typename TInputImage, typename TOutput>
146 {
147  return dynamic_cast< TOutput* >( this->itk::ProcessObject::GetOutput(1) );
148 }
149 
150 
151 template <typename TInput, typename TOutput>
152 void
155 {
156  // find out if graph fit into memory
157  bool mem = this->CheckMemory();
158 
159  if(mem)
160  {
161  // process graph with without downsampling
162  this->ProcessGraphCut();
163  }
164  else
165  {
166  // if graph does not fit into memory, decimate images until we can process the graph cut.
167  this->FindDownsampleFactor(); // find out how much we need to downsample
168 
169  int current_count = m_Count;
170 
171  std::cout << "-- Downsampling images for the graph... " << std::endl;
172  this->InitResampleFilters();
173 
174  typename TInput::SizeType outputSize;
175  typename TInput::SpacingType outputSpacing;
176 
177  while(current_count != -1)
178  {
179  // Resize
180  outputSize[0] = m_InputSize[0] / m_DownsamplingFactor;
181  outputSize[1] = m_InputSize[1] / m_DownsamplingFactor;
182  outputSize[2] = m_InputSize[2] / m_DownsamplingFactor;
183 
184  outputSpacing[0] = this->GetMask()->GetSpacing()[0] * (static_cast<double>(m_InputSize[0]) / static_cast<double>(outputSize[0]));
185  outputSpacing[1] = this->GetMask()->GetSpacing()[1] * (static_cast<double>(m_InputSize[1]) / static_cast<double>(outputSize[1]));
186  outputSpacing[2] = this->GetMask()->GetSpacing()[2] * (static_cast<double>(m_InputSize[2]) / static_cast<double>(outputSize[2]));
187 
188  m_ResampleMask->SetSize(outputSize);
189  m_ResampleMask->SetOutputSpacing(outputSpacing);
190  m_ResampleMask->Update();
191 
192  m_ResampleSources->SetSize(outputSize);
193  m_ResampleSources->SetOutputSpacing(outputSpacing);
194 
195  m_ResampleSinks->SetSize(outputSize);
196  m_ResampleSinks->SetOutputSpacing(outputSpacing);
197 
198  m_Resample1->SetSize(outputSize);
199  m_Resample1->SetOutputSpacing(outputSpacing);
200 
201  m_Resample2->SetSize(outputSize);
202  m_Resample2->SetOutputSpacing(outputSpacing);
203 
204  m_Resample3->SetSize(outputSize);
205  m_Resample3->SetOutputSpacing(outputSpacing);
206 
207  m_Resample4->SetSize(outputSize);
208  m_Resample4->SetOutputSpacing(outputSpacing);
209 
210  m_Resample5->SetSize(outputSize);
211  m_Resample5->SetOutputSpacing(outputSpacing);
212 
213  if(current_count == m_Count)
214  {
215  m_CurrentMask = m_ResampleMask->GetOutput();
216  }
217 
218  // Process graph cut with downsampled images
219  this->ProcessDownsampledGraphCut(current_count);
220 
221  // Process the graph cut with the following downsampling factor
222  m_DownsamplingFactor/=2;
223 
224  typename TInput::SizeType upSize;
225  upSize[0] = m_InputSize[0] / m_DownsamplingFactor;
226  upSize[1] = m_InputSize[1] / m_DownsamplingFactor;
227  upSize[2] = m_InputSize[2] / m_DownsamplingFactor;
228 
229  if(current_count > 0)
230  {
231  TMask::Pointer upImage = Upsample(m_NLinksFilterDecim->GetOutput(), upSize, m_OutputDirection, m_OutputOrigin);
232 
233  typename TMask::SpacingType maskSpacing;
234  maskSpacing[0] = upImage->GetSpacing()[0];
235  maskSpacing[1] = upImage->GetSpacing()[1];
236  maskSpacing[2] = upImage->GetSpacing()[2];
237 
238  m_ResampleMask->SetSize(upSize);
239  m_ResampleMask->SetOutputSpacing(maskSpacing);
240  m_ResampleMask->Update();
241 
242  // Dilate image inside the mask
243  m_CurrentMask = Dilate(upImage, m_ResampleMask->GetOutput());
244  }
245 
246  current_count--;
247  }
248  }
249 }
250 
251 template <typename TInput, typename TOutput>
252 bool
255 {
256  unsigned int nb_vox = 0;
257  bool mem = true;
258  MaskRegionConstIteratorType maskIt (this->GetMask(),this->GetMask()->GetLargestPossibleRegion() );
259  while (!maskIt.IsAtEnd())
260  {
261  if (maskIt.Get() != 0)
262  {
263  nb_vox++;
264  }
265  ++maskIt;
266  }
267 
268  unsigned int nb_edges = 7*nb_vox;
269 
270  try
271  {
272  m_graph = new GraphType(nb_vox, nb_edges);
273  }
274  catch (std::bad_alloc& ba)
275  {
276  std::cerr << "-- In Graph3DFilter: insufficient memory to create the graph: " << ba.what() << '\n';
277  mem = false;
278  }
279 
280  if (m_graph)
281  {
282  delete m_graph;
283  m_graph = NULL;
284  }
285  return mem;
286 }
287 
288 template <typename TInput, typename TOutput>
289 void
292 {
293  m_NLinksFilter->SetSigma( this->GetSigma() );
294  m_NLinksFilter->SetUseSpectralGradient( this->GetUseSpectralGradient() );
295  m_NLinksFilter->SetMatrix( this->GetMatrix() );
296  m_NLinksFilter->SetMatFilename( this->GetMatFilename() );
297  m_NLinksFilter->SetVerbose( this->GetVerbose() );
298  m_NLinksFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
299  m_NLinksFilter->SetTol( m_Tol );
300 
301  m_NLinksFilter->SetInputSeedProbaSources( this->GetInputSeedProbaSources() );
302  m_NLinksFilter->SetInputSeedProbaSinks( this->GetInputSeedProbaSinks() );
303  m_NLinksFilter->SetMask( this->GetMask() );
304 
305  m_NLinksFilter->SetInputImage1( this->GetInputImage1() );
306  if(m_IndexImage2!=m_NbMaxImages){m_NLinksFilter->SetInputImage2( this->GetInputImage2() );}
307  if(m_IndexImage3!=m_NbMaxImages){m_NLinksFilter->SetInputImage3( this->GetInputImage3() );}
308  if(m_IndexImage4!=m_NbMaxImages){m_NLinksFilter->SetInputImage4( this->GetInputImage4() );}
309  if(m_IndexImage5!=m_NbMaxImages){m_NLinksFilter->SetInputImage5( this->GetInputImage5() );}
310 
311  m_NLinksFilter->GraftNthOutput( 0, this->GetOutput() );
312  m_NLinksFilter->GraftNthOutput( 1, this->GetOutputBackground() );
313 
314  try
315  {
316  m_NLinksFilter->Update();
317  }
318  catch (itk::ExceptionObject &e)
319  {
320  std::cerr << e << std::endl;
321  exit(-1);
322  }
323 
324  this->GraftNthOutput( 0 , m_NLinksFilter->GetOutput() );
325  this->GraftNthOutput( 1 , m_NLinksFilter->GetOutputBackground() );
326 }
327 
328 template <typename TInput, typename TOutput>
329 void
331 {
332  m_DownsamplingFactor = 2.0;
333  m_Count = 1;
334 
335  bool mem2 = false;
336  while(!mem2)
337  {
338  typename TInput::SizeType inputSize = this->GetMask()->GetLargestPossibleRegion().GetSize();
339  typename TInput::DirectionType outputDirection = this->GetMask()->GetDirection();
340  typename TInput::PointType outputOrigin = this->GetMask()->GetOrigin();
341 
342  // Resize images
343  typename TInput::SizeType outputSize;
344  outputSize[0] = inputSize[0] / m_DownsamplingFactor;
345  outputSize[1] = inputSize[1] / m_DownsamplingFactor;
346  outputSize[2] = inputSize[2] / m_DownsamplingFactor;
347  typename TInput::SpacingType outputSpacing;
348  outputSpacing[0] = this->GetMask()->GetSpacing()[0] * (static_cast<double>(inputSize[0]) / static_cast<double>(outputSize[0]));
349  outputSpacing[1] = this->GetMask()->GetSpacing()[1] * (static_cast<double>(inputSize[1]) / static_cast<double>(outputSize[1]));
350  outputSpacing[2] = this->GetMask()->GetSpacing()[2] * (static_cast<double>(inputSize[2]) / static_cast<double>(outputSize[2]));
351 
352  ResampleImageFilterMaskType::Pointer resampleMask = ResampleImageFilterMaskType::New();
353  resampleMask->SetOutputDirection( outputDirection );
354  resampleMask->SetOutputOrigin( outputOrigin );
355  resampleMask->SetInput( this->GetMask() );
356  resampleMask->SetSize( outputSize );
357  resampleMask->SetOutputSpacing( outputSpacing );
358  resampleMask->SetTransform( TransformType::New() );
359  resampleMask->SetNumberOfWorkUnits( this->GetNumberOfWorkUnits() );
360  resampleMask->SetCoordinateTolerance( m_Tol );
361  resampleMask->SetDirectionTolerance( m_Tol );
362  resampleMask->Update();
363 
364  unsigned int nb_vox = 0;
365  MaskRegionIteratorType maskIt (resampleMask->GetOutput(),resampleMask->GetOutput()->GetLargestPossibleRegion() );
366  while (!maskIt.IsAtEnd())
367  {
368  if (maskIt.Get() != 0)
369  {
370  nb_vox++;
371  }
372  ++maskIt;
373  }
374 
375  unsigned int nb_edges = 7*nb_vox;
376  mem2 = true;
377  try
378  {
379  m_graph = new GraphType(nb_vox, nb_edges);
380  }
381  catch (std::bad_alloc& ba)
382  {
383  std::cerr << "-- In Graph3DFilter: insufficient memory to create the graph ICI: " << ba.what() << '\n';
384  mem2 = false;
385  m_Count++;
386  m_DownsamplingFactor*=2.0;
387  }
388 
389  if (m_graph)
390  {
391  delete m_graph;
392  m_graph = NULL;
393  }
394  }
395 }
396 
397 template <typename TInput, typename TOutput>
398 void
401 {
402  m_InputSize = this->GetMask()->GetLargestPossibleRegion().GetSize();
403  m_OutputDirection = this->GetMask()->GetDirection();
404  m_OutputOrigin = this->GetMask()->GetOrigin();
405 
406  m_ResampleMask = ResampleImageFilterMaskType::New();
407  m_ResampleMask->SetOutputDirection(m_OutputDirection);
408  m_ResampleMask->SetOutputOrigin(m_OutputOrigin);
409  m_ResampleMask->SetInput(this->GetMask());
410  m_ResampleMask->SetTransform(TransformType::New());
411  m_ResampleMask->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
412  m_ResampleMask->SetCoordinateTolerance(m_Tol);
413  m_ResampleMask->SetDirectionTolerance(m_Tol);
414 
415  m_ResampleSources = ResampleImageFilterdoubleType::New();
416  m_ResampleSources->SetOutputDirection(m_OutputDirection);
417  m_ResampleSources->SetOutputOrigin(m_OutputOrigin);
418  m_ResampleSources->SetInput(this->GetInputSeedProbaSources());
419  m_ResampleSources->SetTransform(TransformType::New());
420  m_ResampleSources->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
421  m_ResampleSources->SetCoordinateTolerance(m_Tol);
422  m_ResampleSources->SetDirectionTolerance(m_Tol);
423 
424  m_ResampleSinks = ResampleImageFilterdoubleType::New();
425  m_ResampleSinks->SetOutputDirection(m_OutputDirection);
426  m_ResampleSinks->SetOutputOrigin(m_OutputOrigin);
427  m_ResampleSinks->SetInput(this->GetInputSeedProbaSinks());
428  m_ResampleSinks->SetTransform(TransformType::New());
429  m_ResampleSinks->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
430  m_ResampleSinks->SetCoordinateTolerance(m_Tol);
431  m_ResampleSinks->SetDirectionTolerance(m_Tol);
432 
433  m_Resample1 = ResampleImageFilterType::New();
434  m_Resample1->SetOutputDirection(m_OutputDirection);
435  m_Resample1->SetOutputOrigin(m_OutputOrigin);
436  m_Resample1->SetInput(this->GetInputImage1());
437  m_Resample1->SetTransform(TransformType::New());
438  m_Resample1->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
439  m_Resample1->SetCoordinateTolerance(m_Tol);
440  m_Resample1->SetDirectionTolerance(m_Tol);
441 
442  m_Resample2 = ResampleImageFilterType::New();
443  m_Resample2->SetOutputDirection(m_OutputDirection);
444  m_Resample2->SetOutputOrigin(m_OutputOrigin);
445  m_Resample2->SetInput(this->GetInputImage2());
446  m_Resample2->SetTransform(TransformType::New());
447  m_Resample2->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
448  m_Resample2->SetCoordinateTolerance(m_Tol);
449  m_Resample2->SetDirectionTolerance(m_Tol);
450 
451  m_Resample3 = ResampleImageFilterType::New();
452  m_Resample3->SetOutputDirection(m_OutputDirection);
453  m_Resample3->SetOutputOrigin(m_OutputOrigin);
454  m_Resample3->SetInput(this->GetInputImage4());
455  m_Resample3->SetTransform(TransformType::New());
456  m_Resample3->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
457  m_Resample3->SetCoordinateTolerance(m_Tol);
458  m_Resample3->SetDirectionTolerance(m_Tol);
459 
460  m_Resample4 = ResampleImageFilterType::New();
461  m_Resample4->SetOutputDirection(m_OutputDirection);
462  m_Resample4->SetOutputOrigin(m_OutputOrigin);
463  m_Resample4->SetInput(this->GetInputImage4());
464  m_Resample4->SetTransform(TransformType::New());
465  m_Resample4->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
466  m_Resample4->SetCoordinateTolerance(m_Tol);
467  m_Resample4->SetDirectionTolerance(m_Tol);
468 
469  m_Resample5 = ResampleImageFilterType::New();
470  m_Resample5->SetOutputDirection(m_OutputDirection);
471  m_Resample5->SetOutputOrigin(m_OutputOrigin);
472  m_Resample5->SetInput(this->GetInputImage5());
473  m_Resample5->SetTransform(TransformType::New());
474  m_Resample5->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
475  m_Resample5->SetCoordinateTolerance(m_Tol);
476  m_Resample5->SetDirectionTolerance(m_Tol);
477 }
478 
479 template <typename TInput, typename TOutput>
480 void
483 {
484  m_NLinksFilterDecim = NLinksFilterType::New();
485  m_NLinksFilterDecim->SetSigma( this->GetSigma() );
486  m_NLinksFilterDecim->SetUseSpectralGradient( this->GetUseSpectralGradient() );
487  m_NLinksFilterDecim->SetMatrix( this->GetMatrix() );
488  m_NLinksFilterDecim->SetMatFilename( this->GetMatFilename() );
489  m_NLinksFilterDecim->SetVerbose( this->GetVerbose() );
490  m_NLinksFilterDecim->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
491  m_NLinksFilterDecim->SetTol( m_Tol );
492  m_NLinksFilterDecim->SetMask( m_CurrentMask ); // mandatory brain mask
493 
494  m_NLinksFilterDecim->SetInputSeedProbaSources( m_ResampleSources->GetOutput() );
495  m_NLinksFilterDecim->SetInputSeedProbaSinks( m_ResampleSinks->GetOutput() );
496 
497  m_NLinksFilterDecim->SetInputImage1( m_Resample1->GetOutput() );
498  if(m_IndexImage2!=m_NbMaxImages){m_NLinksFilterDecim->SetInputImage2( m_Resample2->GetOutput() );}
499  if(m_IndexImage3!=m_NbMaxImages){m_NLinksFilterDecim->SetInputImage3( m_Resample3->GetOutput() );}
500  if(m_IndexImage4!=m_NbMaxImages){m_NLinksFilterDecim->SetInputImage4( m_Resample4->GetOutput() );}
501  if(m_IndexImage5!=m_NbMaxImages){m_NLinksFilterDecim->SetInputImage5( m_Resample5->GetOutput() );}
502 
503  if(current_count==0)
504  {
505  m_NLinksFilterDecim->GraftNthOutput( 0, this->GetOutput() );
506  m_NLinksFilterDecim->GraftNthOutput( 1, this->GetOutputBackground() );
507  }
508 
509  try
510  {
511  m_NLinksFilterDecim->Update();
512  }
513  catch (itk::ExceptionObject &e)
514  {
515  std::cerr << e << std::endl;
516  exit(-1);
517  }
518 
519  if(current_count==0)
520  {
521  this->GraftNthOutput( 0 , m_NLinksFilterDecim->GetOutput() );
522  this->GraftNthOutput( 1 , m_NLinksFilterDecim->GetOutputBackground() );
523  }
524 }
525 
526 template <typename TInput, typename TOutput>
527 itk::Image <unsigned char,3>::Pointer
529 ::Dilate(const TMask* input, const TMask* mask)
530 {
531  unsigned int radius = 4;
532 
533  typedef itk::BinaryBallStructuringElement<TMask::PixelType, TMask::ImageDimension> StructuringElementType;
534  StructuringElementType structuringElement;
535  structuringElement.SetRadius(radius);
536  structuringElement.CreateStructuringElement();
537 
538  typedef itk::GrayscaleDilateImageFilter <TMask,TMask,StructuringElementType> DilateFilterType;
539  DilateFilterType::Pointer dilateFilter = DilateFilterType::New();
540  dilateFilter->SetInput(input);
541  dilateFilter->SetKernel(structuringElement);
542  dilateFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
543  dilateFilter->SetCoordinateTolerance(m_Tol);
544  dilateFilter->SetDirectionTolerance(m_Tol);
545 
546  typedef itk::MaskImageFilter< TMask, TMask > MaskFilterType;
547  MaskFilterType::Pointer maskFilter = MaskFilterType::New();
548  maskFilter->SetInput(dilateFilter->GetOutput());
549  maskFilter->SetMaskImage(mask);
550  maskFilter->Update();
551 
552  return maskFilter->GetOutput();
553 }
554 
555 template <typename TInput, typename TOutput>
556 itk::Image <unsigned char,3>::Pointer
558 ::Upsample(const TMask* input, TMask::SizeType upSize, TMask::DirectionType outputDirection, TMask::PointType outputOrigin)
559 {
560  // We don't want any transform on our image except rescaling which is not
561  // specified by a transform but by the input/output spacing.
562  // So no transform will be specified.
563  typedef itk::IdentityTransform<double, 3> T_Transform;
564  // Instantiate the transform and specify it should be the id transform.
565  T_Transform::Pointer Transform = T_Transform::New();
566  Transform->SetIdentity();
567 
568  // If ITK resampler determines there is something to interpolate which is
569  // usually the case when upscaling (!) then we must specify the interpolation
570  // algorithm. In our case, we want bicubic interpolation. One way to implement
571  // it is with a third order b-spline. So the type is specified here and the
572  // order will be specified with a method call later on.
573 
574  typedef itk::BSplineInterpolateImageFunction<TMask, double, double> T_Interpolator;
575  // Instantiate the b-spline interpolator and set it as the third orderfor bicubic.
576  T_Interpolator::Pointer Interpolator = T_Interpolator::New();
577  Interpolator->SetSplineOrder(3);
578 
579  // The resampler type itself.
580  typedef itk::ResampleImageFilter<TMask, TMask> T_ResampleFilter;
581  // Instantiate the resampler. Wire in the transform and the interpolator.
582  T_ResampleFilter::Pointer ResizeFilter = T_ResampleFilter::New();
583  ResizeFilter->SetTransform(Transform);
584  ResizeFilter->SetInterpolator(Interpolator);
585 
586  // Set the output origin.
587  ResizeFilter->SetOutputOrigin(m_OutputOrigin);
588  ResizeFilter->SetOutputDirection(m_OutputDirection);
589 
590  // Compute and set the output spacing
591  // Compute the output spacing from input spacing and old and new sizes.
592  //
593  // The computation must be so that the following holds:
594  //
595  // new width old x spacing
596  // ---------- = ---------------
597  // old width new x spacing
598  //
599  // we specify new height and width and compute new spacings
600 
601  unsigned int nNewSize_0 = upSize[0];
602  unsigned int nNewSize_1 = upSize[1];
603  unsigned int nNewSize_2 = upSize[2];
604 
605  // Fetch original image size.
606  const TMask::RegionType& inputRegion = input->GetLargestPossibleRegion();
607  const TMask::SizeType& vnInputSize = inputRegion.GetSize();
608  unsigned int nOldSize_0 = vnInputSize[0];
609  unsigned int nOldSize_1 = vnInputSize[1];
610  unsigned int nOldSize_2 = vnInputSize[2];
611 
612  // Fetch original image spacing.
613  const TMask::SpacingType& vfInputSpacing = input->GetSpacing();
614 
615  double vfOutputSpacing[3];
616  vfOutputSpacing[0] = vfInputSpacing[0] * (double) nOldSize_0 / (double) nNewSize_0;
617  vfOutputSpacing[1] = vfInputSpacing[1] * (double) nOldSize_1 / (double) nNewSize_1;
618  vfOutputSpacing[2] = vfInputSpacing[2] * (double) nOldSize_2 / (double) nNewSize_2;
619 
620  // Set the output spacing.
621  ResizeFilter->SetOutputSpacing(vfOutputSpacing);
622 
623  // Set the output size
624  itk::Size<3> vnOutputSize = { {nNewSize_0, nNewSize_1, nNewSize_2} };
625  ResizeFilter->SetSize(vnOutputSize);
626  ResizeFilter->SetInput(input);
627  ResizeFilter->SetCoordinateTolerance(m_Tol);
628  ResizeFilter->SetDirectionTolerance(m_Tol);
629  ResizeFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
630  ResizeFilter->Update();
631 
632  return ResizeFilter->GetOutput();
633 }
634 
635 } //end of namespace anima
TSeedProba::ConstPointer GetInputSeedProbaSinks()
TSeedProba::ConstPointer GetInputSeedProbaSources()
void SetInputSeedProbaSources(const TSeedProba *mask)
Class allowing the decimation of the images if necessary (if 3D graph size causes memory problems)...
itk::Image< PixelTypeUC, 3 > TMask
void SetMask(const TMask *mask)
void SetInputSeedProbaSinks(const TSeedProba *mask)
TMask::Pointer Dilate(const TMask *input, const TMask *mask)
TMask::ConstPointer GetMask()
itk::ImageRegionConstIterator< TMask > MaskRegionConstIteratorType
TInput::ConstPointer GetInputImage4()
void ProcessDownsampledGraphCut(int current_count)
void GenerateData() ITK_OVERRIDE
TInput::ConstPointer GetInputImage2()
TInput::ConstPointer GetInputImage3()
TMask::Pointer Upsample(const TMask *input, TMask::SizeType upSize, TMask::DirectionType outputDirection, TMask::PointType outputOrigin)
void SetInputImage(unsigned int i, const TInput *image)
itk::Image< PixelTypeD, 3 > TSeedProba
TInput::ConstPointer GetInputImage1()
TInput::ConstPointer GetInputImage5()
itk::ImageRegionIterator< TMask > MaskRegionIteratorType