ANIMA  4.0
animaMCMAverageImagesImageFilter.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 MCMAverageImagesImageFilter <TPixelType>
13 ::BeforeThreadedGenerateData()
14 {
15  this->Superclass::BeforeThreadedGenerateData();
16 
17  m_ReferenceInputModels.resize(this->GetNumberOfIndexedInputs());
18 
19  InputImageType *input = const_cast <InputImageType *> (this->GetInput(0));
20  MCModelPointer tmpMCM = input->GetDescriptionModel();
21  unsigned int numIsoCompartments = tmpMCM->GetNumberOfIsotropicCompartments();
22 
23  std::vector <anima::DiffusionModelCompartmentType> refIsoVec(numIsoCompartments);
24  for (unsigned int j = 0;j < numIsoCompartments;++j)
25  refIsoVec[j] = m_ReferenceOutputModel->GetCompartment(j)->GetCompartmentType();
26 
27  for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i)
28  {
29  InputImageType *input = const_cast <InputImageType *> (this->GetInput(i));
30  MCModelPointer model = input->GetDescriptionModel();
31 
32  if (model->GetNumberOfIsotropicCompartments() != numIsoCompartments)
33  itkExceptionMacro("All input images should have the same isotropic compartment structure as the output model");
34 
35  for (unsigned int j = 0;j < numIsoCompartments;++j)
36  {
37  if (refIsoVec[j] != model->GetCompartment(j)->GetCompartmentType())
38  itkExceptionMacro("All input images should have the same isotropic compartment structure as the output model");
39  }
40 
41  m_ReferenceInputModels[i] = model;
42  }
43 
44  this->GetOutput()->SetDescriptionModel(m_ReferenceOutputModel);
45 }
46 
47 template <class TPixelType>
48 void
50 ::SetReferenceOutputModel(MCModelPointer &model)
51 {
52  m_ReferenceOutputModel = model;
53 }
54 
55 template <class TPixelType>
58 ::CreateAverager()
59 {
61  mcmAverager->SetOutputModel(m_ReferenceOutputModel);
62 
63  return mcmAverager;
64 }
65 
66 template <class TPixelType>
67 void
69 ::DynamicThreadedGenerateData(const InputRegionType &region)
70 {
71  unsigned int numInputs = this->GetNumberOfIndexedInputs();
72 
73  typedef itk::ImageRegionConstIterator <InputImageType> InputIteratorType;
74  typedef itk::ImageRegionIterator <OutputImageType> OutputIteratorType;
75 
76  std::vector<InputIteratorType> inputIterators (numInputs);
77 
78  for (unsigned int i = 0; i < numInputs; ++i)
79  inputIterators[i] = InputIteratorType (this->GetInput(i), region);
80 
81  OutputIteratorType outputIterator (this->GetOutput(), region);
82 
83  MCMAveragerPointer mcmAverager = this->CreateAverager();
84 
85  std::vector <MCModelPointer> workInputModels(numInputs);
86  for (unsigned int i = 0;i < numInputs;++i)
87  workInputModels[i] = m_ReferenceInputModels[i]->Clone();
88 
89  std::vector <double> workInputWeights(numInputs,0.0);
90  PixelType voxelOutputValue(m_ReferenceOutputModel->GetSize());
91 
92  using MaskIteratorType = itk::ImageRegionConstIterator <MaskImageType>;
93  std::vector <MaskIteratorType> maskIterators;
94  if (m_MaskImages.size() == numInputs)
95  {
96  for (unsigned int i = 0; i < numInputs; ++i)
97  maskIterators.push_back(MaskIteratorType (m_MaskImages[i], region));
98  }
99 
100  unsigned int numMasks = maskIterators.size();
101  while (!inputIterators[0].IsAtEnd())
102  {
103  voxelOutputValue.Fill(0.0);
104  std::fill(workInputWeights.begin(),workInputWeights.end(),0.0);
105 
106  unsigned int trueInputNumber = 0;
107  for (unsigned int i = 0;i < numInputs;++i)
108  {
109  if (isZero(inputIterators[i].Get()))
110  continue;
111 
112  if (numMasks == numInputs)
113  {
114  if (maskIterators[i].Get() == 0)
115  continue;
116  }
117 
118  workInputModels[i]->SetModelVector(inputIterators[i].Get());
119  workInputWeights[i] = 1.0;
120 
121  ++trueInputNumber;
122  }
123 
124  if (trueInputNumber == 0)
125  {
126  outputIterator.Set(voxelOutputValue);
127  for (unsigned int i = 0;i < numInputs;++i)
128  ++inputIterators[i];
129  ++outputIterator;
130 
131  for (unsigned int i = 0;i < numMasks;++i)
132  ++maskIterators[i];
133 
134  continue;
135  }
136 
137  mcmAverager->SetInputModels(workInputModels);
138  mcmAverager->SetInputWeights(workInputWeights);
139 
140  // weights are normalized inside averager
141  mcmAverager->Update();
142 
143  voxelOutputValue = mcmAverager->GetOutputModel()->GetModelVector();
144 
145  //Set output value then upgrade iterators
146  outputIterator.Set(voxelOutputValue);
147  for (unsigned int i = 0;i < numInputs;++i)
148  ++inputIterators[i];
149  ++outputIterator;
150 
151  for (unsigned int i = 0;i < numMasks;++i)
152  ++maskIterators[i];
153  }
154 }
155 
156 } // end namespace anima
MCMPointer & GetDescriptionModel()