ANIMA  4.0
animaBaseBlockMatcher.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <animaBobyqaOptimizer.h>
7 #include <itkPoolMultiThreader.h>
8 
9 namespace anima
10 {
11 
12 template <typename TInputImageType>
13 BaseBlockMatcher <TInputImageType>
15 {
16  m_ForceComputeBlocks = false;
17  m_NumberOfThreads = itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads();
18 
19  m_BlockPercentageKept = 0.8;
20  m_BlockSize = 5;
21  m_BlockSpacing = 2;
22 
23  m_BlockVarianceThreshold = 25;
24 
25  m_SearchRadius = 2;
26  m_FinalRadius = 0.001;
27  m_OptimizerMaximumIterations = 100;
28  m_StepSize = 1.0;
29 
30  m_OptimizerType = Bobyqa;
31  m_Verbose = true;
32 
33  m_HighestProcessedBlock = 0;
34 }
35 
36 template <typename TInputImageType>
37 void
39 ::InitializeBlocks()
40 {
41  // Init blocks on reference image
42  typedef typename TInputImageType::IOPixelType InputPixelType;
44  typedef typename InitializerType::Pointer InitializerPointer;
45 
46  InitializerPointer initPtr = InitializerType::New();
47  initPtr->AddReferenceImage(m_ReferenceImage);
48 
49  if (m_NumberOfThreads != 0)
50  initPtr->SetNumberOfThreads(m_NumberOfThreads);
51 
52  initPtr->SetPercentageKept(m_BlockPercentageKept);
53  initPtr->SetBlockSize(m_BlockSize);
54  initPtr->SetBlockSpacing(m_BlockSpacing);
55  initPtr->SetScalarVarianceThreshold(m_BlockVarianceThreshold);
56  initPtr->SetOrientedModelVarianceThreshold(m_BlockVarianceThreshold);
57  initPtr->AddGenerationMask(m_BlockGenerationMask);
58 
59  initPtr->SetRequestedRegion(m_ReferenceImage->GetLargestPossibleRegion());
60 
61  m_BlockRegions = initPtr->GetOutput();
62  m_BlockPositions = initPtr->GetOutputPositions();
63 
64  if (m_Verbose)
65  std::cout << "Generated " << m_BlockRegions.size() << " blocks..." << std::endl;
66 
67  m_BlockTransformPointers.resize(m_BlockRegions.size());
68  m_BlockWeights.resize(m_BlockRegions.size());
69  for (unsigned int i = 0;i < m_BlockRegions.size();++i)
70  m_BlockTransformPointers[i] = this->GetNewBlockTransform(m_BlockPositions[i]);
71 }
72 
73 template <typename TInputImageType>
76 ::SetupOptimizer()
77 {
78  OptimizerPointer optimizer;
79  switch (m_OptimizerType)
80  {
81  case Bobyqa:
82  {
83  typedef anima::BobyqaOptimizer LocalOptimizerType;
84  optimizer = LocalOptimizerType::New();
85  LocalOptimizerType *tmpOpt = (LocalOptimizerType *)optimizer.GetPointer();
86  tmpOpt->SetRhoBegin(m_SearchRadius);
87  tmpOpt->SetRhoEnd(m_FinalRadius);
88 
89  tmpOpt->SetNumberSamplingPoints(m_BlockTransformPointers[0]->GetNumberOfParameters() + 2);
90  tmpOpt->SetMaximumIteration(m_OptimizerMaximumIterations);
91  tmpOpt->SetMaximize(this->GetMaximizedMetric());
92 
93  break;
94  }
95 
96  case Exhaustive:
97  default:
98  {
99  typedef anima::VoxelExhaustiveOptimizer LocalOptimizerType;
100  optimizer = LocalOptimizerType::New();
101 
102  LocalOptimizerType *tmpOpt = (LocalOptimizerType *)optimizer.GetPointer();
103  LocalOptimizerType::StepsType steps(m_BlockTransformPointers[0]->GetNumberOfParameters());
104 
105  typename InputImageType::SpacingType fixedSpacing = m_ReferenceImage->GetSpacing();
106  typename InputImageType::DirectionType fixedDirection = m_ReferenceImage->GetDirection();
107 
108  vnl_matrix <double> geometry(InputImageType::ImageDimension,InputImageType::ImageDimension);
109  for (unsigned int i = 0;i < InputImageType::ImageDimension;++i)
110  {
111  double tmpVoxSize = fixedSpacing[i];
112 
113  for (unsigned int j = 0;j < InputImageType::ImageDimension;++j)
114  geometry(j,i) = tmpVoxSize*fixedDirection(j,i);
115  }
116 
117  tmpOpt->SetGeometry(geometry);
118 
119  LocalOptimizerType::ScalesType tmpScales(InputImageType::ImageDimension);
120 
121  for (unsigned i = 0; i < steps.Size(); ++i)
122  {
123  steps[i] = m_SearchRadius;
124  tmpScales[i] = m_StepSize;
125  }
126 
127  tmpOpt->SetNumberOfSteps(steps);
128  tmpOpt->SetScales(tmpScales);
129  tmpOpt->SetMaximize(this->GetMaximizedMetric());
130  break;
131  }
132  }
133 
134  this->TransformDependantOptimizerSetup(optimizer);
135 
136  return optimizer;
137 }
138 
139 template <typename TInputImageType>
140 void
142 ::Update()
143 {
144  // Generate blocks if needed on reference image
145  if ((m_ForceComputeBlocks) || (m_BlockTransformPointers.size() == 0))
146  this->InitializeBlocks();
147 
148  m_HighestProcessedBlock = 0;
149  itk::PoolMultiThreader::Pointer threadWorker = itk::PoolMultiThreader::New();
150  ThreadedMatchData *tmpStr = new ThreadedMatchData;
151  tmpStr->BlockMatch = this;
152 
153  threadWorker->SetNumberOfWorkUnits(m_NumberOfThreads);
154  threadWorker->SetSingleMethod(this->ThreadedMatching,tmpStr);
155  threadWorker->SingleMethodExecute();
156 
157  delete tmpStr;
158 }
159 
160 template <typename TInputImageType>
161 ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
163 ::ThreadedMatching(void *arg)
164 {
165  itk::MultiThreaderBase::WorkUnitInfo *threadArgs = (itk::MultiThreaderBase::WorkUnitInfo *)arg;
166  ThreadedMatchData* data = (ThreadedMatchData *)threadArgs->UserData;
167 
168  data->BlockMatch->ProcessBlockMatch();
169  return ITK_THREAD_RETURN_DEFAULT_VALUE;
170 }
171 
172 template <typename TInputImageType>
173 void
175 ::ProcessBlockMatch()
176 {
177  bool continueLoop = true;
178  unsigned int highestToleratedBlockIndex = m_BlockRegions.size();
179 
180  unsigned int stepData = std::ceil(m_BlockRegions.size() / (10 * this->GetNumberOfWorkUnits()));
181  unsigned int minNumBlocks = std::min(highestToleratedBlockIndex,static_cast <unsigned int> (10));
182  stepData = std::max(minNumBlocks,stepData);
183 
184  while (continueLoop)
185  {
186  m_LockHighestProcessedBlock.lock();
187 
188  if (m_HighestProcessedBlock >= highestToleratedBlockIndex)
189  {
190  m_LockHighestProcessedBlock.unlock();
191  continueLoop = false;
192  continue;
193  }
194 
195  unsigned int startPoint = m_HighestProcessedBlock;
196  unsigned int endPoint = m_HighestProcessedBlock + stepData;
197  if (endPoint > highestToleratedBlockIndex)
198  endPoint = highestToleratedBlockIndex;
199 
200  m_HighestProcessedBlock = endPoint;
201 
202  m_LockHighestProcessedBlock.unlock();
203 
204  this->BlockMatch(startPoint,endPoint);
205  }
206 }
207 
208 template <typename TInputImageType>
209 void
211 ::BlockMatch(unsigned int startIndex, unsigned int endIndex)
212 {
213  MetricPointer metric = this->SetupMetric();
214  OptimizerPointer optimizer = this->SetupOptimizer();
215 
216  // Loop over the desired blocks
217  for (unsigned int block = startIndex;block < endIndex;++block)
218  {
219  this->BlockMatchingSetup(metric, block);
220  optimizer->SetCostFunction(metric);
221  optimizer->SetInitialPosition(m_BlockTransformPointers[block]->GetParameters());
222 
223  try
224  {
225  optimizer->StartOptimization();
226  }
227  catch (itk::ExceptionObject & err)
228  {
229  m_BlockWeights[block] = 0;
230  continue;
231  }
232 
233  m_BlockTransformPointers[block]->SetParameters(optimizer->GetCurrentPosition());
234 
235  double val = optimizer->GetValue(optimizer->GetCurrentPosition());
236  m_BlockWeights[block] = this->ComputeBlockWeight(val,block);
237  }
238 }
239 
240 } // end namespace anima
MetricType::Pointer MetricPointer
void SetNumberOfWorkUnits(unsigned int val)
OptimizerType::Pointer OptimizerPointer