ANIMA  4.0
animaMCMScalarMapsImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkImageRegionConstIterator.h>
5 #include <itkImageRegionIterator.h>
6 
7 namespace anima
8 {
9 
10 template <class TPixelType>
11 void
12 MCMScalarMapsImageFilter <TPixelType>
13 ::DynamicThreadedGenerateData(const InputRegionType &outputRegionForThread)
14 {
15  typedef itk::ImageRegionConstIterator <InputImageType> InputIteratorType;
16  typedef itk::ImageRegionIterator <OutputImageType> OutputIteratorType;
17 
18  InputIteratorType inItr(this->GetInput(), outputRegionForThread);
19  unsigned int numOutputs = this->GetNumberOfIndexedOutputs();
20  std::vector <OutputIteratorType> outItrs(numOutputs);
21 
22  for (unsigned int i = 0;i < numOutputs;++i)
23  outItrs[i] = OutputIteratorType(this->GetOutput(i), outputRegionForThread);
24 
25  InputImageType *input = const_cast <InputImageType *> (this->GetInput());
26  MCModelPointer mcmPtr = input->GetDescriptionModel()->Clone();
27  PixelType voxelValue;
28  unsigned int numCompartments = mcmPtr->GetNumberOfCompartments();
29  unsigned int numIsoCompartments = mcmPtr->GetNumberOfIsotropicCompartments();
30 
31  while (!inItr.IsAtEnd())
32  {
33  voxelValue = inItr.Get();
34  if (isZero(voxelValue))
35  {
36  ++inItr;
37  for (unsigned int i = 0;i < numOutputs;++i)
38  {
39  outItrs[i].Set(0.0);
40  ++outItrs[i];
41  }
42 
43  continue;
44  }
45 
46  mcmPtr->SetModelVector(voxelValue);
47  double anisoWeight = 0.0;
48  double faCompartments = 0.0;
49  double mdCompartments = 0.0;
50  double parDiffCompartments = 0.0;
51  double perpDiffCompartments = 0.0;
52  double isoRWeight = 0.0;
53  double fwWeight = 0.0;
54 
55  for (unsigned int i = numIsoCompartments;i < numCompartments;++i)
56  {
57  double weight = mcmPtr->GetCompartmentWeight(i);
58  if (weight <= 0.0)
59  continue;
60  anima::BaseCompartment *currentComp = mcmPtr->GetCompartment(i);
61  anisoWeight += weight;
62  faCompartments += weight * currentComp->GetApparentFractionalAnisotropy();
63  mdCompartments += weight * currentComp->GetApparentMeanDiffusivity();
64  parDiffCompartments += weight * currentComp->GetApparentParallelDiffusivity();
65  perpDiffCompartments += weight * currentComp->GetApparentPerpendicularDiffusivity();
66  }
67 
68  for (unsigned int i = 0;i < numIsoCompartments;++i)
69  {
70  double weight = mcmPtr->GetCompartmentWeight(i);
71  if (weight <= 0.0)
72  continue;
73 
74  anima::BaseCompartment *currentComp = mcmPtr->GetCompartment(i);
75 
76  if (m_IncludeIsotropicWeights)
77  {
78  mdCompartments += weight * currentComp->GetApparentMeanDiffusivity();
79  parDiffCompartments += weight * currentComp->GetApparentParallelDiffusivity();
80  perpDiffCompartments += weight * currentComp->GetApparentPerpendicularDiffusivity();
81  }
82 
83  if (currentComp->GetCompartmentType() == anima::FreeWater)
84  fwWeight += weight;
85  else if (currentComp->GetCompartmentType() == anima::IsotropicRestrictedWater)
86  isoRWeight += weight;
87  else if (currentComp->GetCompartmentType() == anima::Stanisz)
88  isoRWeight += weight;
89  }
90 
91  if ((!m_IncludeIsotropicWeights)&&(anisoWeight > 0.0))
92  {
93  mdCompartments /= anisoWeight;
94  faCompartments /= anisoWeight;
95  parDiffCompartments /= anisoWeight;
96  perpDiffCompartments /= anisoWeight;
97  }
98 
99  outItrs[0].Set(fwWeight);
100  outItrs[1].Set(isoRWeight);
101  outItrs[2].Set(anisoWeight);
102  outItrs[3].Set(faCompartments);
103  outItrs[4].Set(mdCompartments);
104  outItrs[5].Set(parDiffCompartments);
105  outItrs[6].Set(perpDiffCompartments);
106 
107  ++inItr;
108  for (unsigned int i = 0;i < numOutputs;++i)
109  ++outItrs[i];
110  }
111 }
112 
113 } // end namespace anima
MCMPointer & GetDescriptionModel()
virtual double GetApparentFractionalAnisotropy()
virtual double GetApparentMeanDiffusivity()
virtual double GetApparentPerpendicularDiffusivity()
Get approximation to perpendicular diffusivity of the compartment (may be different from radial diffu...
virtual DiffusionModelCompartmentType GetCompartmentType()=0
Utility function to return compartment type without needing dynamic casts.
virtual double GetApparentParallelDiffusivity()
Get approximation to parallel diffusivity of the compartment (may be different from axial diffusivity...