ANIMA  4.0
animaPyramidalBlockMatchingBridge.hxx
Go to the documentation of this file.
1 #pragma once
2 #include "animaPyramidalBlockMatchingBridge.h"
3 
4 #include <animaReadWriteFunctions.h>
5 #include <itkTransformFileReader.h>
6 #include <itkTransformFileWriter.h>
7 #include <itkCenteredTransformInitializer.h>
8 #include <itkMinimumMaximumImageFilter.h>
9 
10 #include <animaAsymmetricBMRegistrationMethod.h>
11 #include <animaSymmetricBMRegistrationMethod.h>
12 #include <animaKissingSymmetricBMRegistrationMethod.h>
13 
14 #include <animaAnatomicalBlockMatcher.h>
15 
16 #include <animaLSWTransformAgregator.h>
17 #include <animaLTSWTransformAgregator.h>
18 #include <animaMEstTransformAgregator.h>
19 
20 namespace anima
21 {
22 
23 template <unsigned int ImageDimension>
24 PyramidalBlockMatchingBridge<ImageDimension>::PyramidalBlockMatchingBridge()
25 {
26  m_InitialTransform = nullptr;
27  m_DirectionTransform = nullptr;
28  m_ReferenceImage = nullptr;
29  m_FloatingImage = nullptr;
30 
31  m_OutputTransform = nullptr;
32  m_outputTransformFile = "";
33  m_outputNearestRigidTransformFile = "";
34  m_outputNearestSimilarityTransformFile = "";
35 
36  m_OutputImage = nullptr;
37 
38  m_ReferenceMinimalValue = 0.0;
39  m_FloatingMinimalValue = 0.0;
40 
41  m_BlockSize = 5;
42  m_BlockSpacing = 5;
43  m_StDevThreshold = 5;
44 
45  m_SymmetryType = Asymmetric;
46  m_Transform = Translation;
47  m_AffineDirection = 1;
48  m_Metric = SquaredCorrelation;
49  m_Optimizer = Bobyqa;
50 
51  m_MaximumIterations = 10;
52  m_MinimalTransformError = 0.01;
53  m_OptimizerMaximumIterations = 100;
54  m_SearchRadius = 2;
55  m_SearchAngleRadius = 5;
56  m_SearchScaleRadius = 0.1;
57  m_FinalRadius = 0.001;
58  m_StepSize = 1;
59  m_TranslateUpperBound = 50;
60  m_AngleUpperBound = 180;
61  m_ScaleUpperBound = 3;
62  m_Agregator = MEstimation;
63  m_OutputTransformType = outRigid;
64  m_AgregThreshold = 0.5;
65  m_SeStoppingThreshold = 0.01;
66  m_NumberOfPyramidLevels = 3;
67  m_LastPyramidLevel = 0;
68  m_PercentageKept = 0.8;
69  m_TransformInitializationType = ClosestTransform;
70 
71  this->SetNumberOfWorkUnits(itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads());
72 
73  m_Abort = false;
74  m_Verbose = true;
75 
76  m_callback = itk::CStyleCommand::New();
77  m_callback->SetClientData ((void *) this);
78  m_callback->SetCallback (ManageProgress);
79 }
80 
81 template <unsigned int ImageDimension>
82 PyramidalBlockMatchingBridge<ImageDimension>::~PyramidalBlockMatchingBridge()
83 {
84 }
85 
86 template <unsigned int ImageDimension>
87 void PyramidalBlockMatchingBridge<ImageDimension>::Abort()
88 {
89  m_Abort = true;
90 
91  if(m_bmreg)
92  m_bmreg->Abort();
93 }
94 
95 template <unsigned int ImageDimension>
96 void PyramidalBlockMatchingBridge<ImageDimension>::SetInitialTransform (std::string initialTransformFile)
97 {
98  if (initialTransformFile != "")
99  {
100  itk::TransformFileReader::Pointer tmpTrRead = itk::TransformFileReader::New();
101  tmpTrRead->SetFileName (initialTransformFile);
102 
103  try
104  {
105  tmpTrRead->Update();
106 
107  itk::TransformFileReader::TransformListType trsfList = * (tmpTrRead->GetTransformList());
108  itk::TransformFileReader::TransformListType::iterator tr_it = trsfList.begin();
109 
110  m_InitialTransform = dynamic_cast <AffineTransformType *> ((*tr_it).GetPointer());
111 
112  }
113  catch (itk::ExceptionObject &e)
114  {
115  std::cerr << "Unable to read initial transform... Exiting..." << std::endl;
116  }
117  }
118 }
119 
120 template <unsigned int ImageDimension>
121 void PyramidalBlockMatchingBridge<ImageDimension>::SetDirectionTransform(std::string directionTransformFile)
122 {
123  if (directionTransformFile != "")
124  {
125  itk::TransformFileReader::Pointer tmpTrRead = itk::TransformFileReader::New();
126  tmpTrRead->SetFileName(directionTransformFile);
127 
128  try
129  {
130  tmpTrRead->Update();
131 
132  itk::TransformFileReader::TransformListType trsfList = *(tmpTrRead->GetTransformList());
133  itk::TransformFileReader::TransformListType::iterator tr_it = trsfList.begin();
134 
135  m_DirectionTransform = dynamic_cast <AffineTransformType *> ((*tr_it).GetPointer());
136 
137  }
138  catch (itk::ExceptionObject &e)
139  {
140  std::cerr << "Unable to read direction transform... Exiting..." << std::endl;
141  }
142  }
143 }
144 
145 template <unsigned int ImageDimension>
146 void PyramidalBlockMatchingBridge<ImageDimension>::Update()
147 {
148  typedef BaseTransformAgregator<ImageDimension> BaseAgreg;
149 
150  m_Abort = false;
151 
152  // progress management
153  m_progressReporter = new itk::ProgressReporter(this, 0, GetNumberOfPyramidLevels()*this->m_MaximumIterations);
154  this->AddObserver(itk::ProgressEvent(), m_progressCallback);
155 
156  this->InvokeEvent(itk::StartEvent());
157 
158  // Compute minimal value of reference and Floating images
159  using MinMaxFilterType = itk::MinimumMaximumImageFilter <InputImageType>;
160  typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New();
161  minMaxFilter->SetInput(m_ReferenceImage);
162  if (this->GetNumberOfWorkUnits() != 0)
163  minMaxFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
164  minMaxFilter->Update();
165 
166  m_ReferenceMinimalValue = minMaxFilter->GetMinimum();
167 
168  minMaxFilter = MinMaxFilterType::New();
169  minMaxFilter->SetInput(m_FloatingImage);
170  if (this->GetNumberOfWorkUnits() != 0)
171  minMaxFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
172  minMaxFilter->Update();
173 
174  m_FloatingMinimalValue = minMaxFilter->GetMinimum();
175 
176  // Only CT images are below zero, little hack to set minimal values to either -1024 or 0
177  if (m_ReferenceMinimalValue < 0.0)
178  m_ReferenceMinimalValue = -1024;
179  else
180  m_ReferenceMinimalValue = 0.0;
181 
182  if (m_FloatingMinimalValue < 0.0)
183  m_FloatingMinimalValue = -1024;
184  else
185  m_FloatingMinimalValue = 0.0;
186 
187  // Set up pyramids of images and masks
188  this->SetupPyramids();
189 
190  typedef anima::AnatomicalBlockMatcher <InputImageType> BlockMatcherType;
191 
192  // Iterate over pyramid levels
193  for (unsigned int i = 0;i < GetNumberOfPyramidLevels() && !m_Abort; ++i)
194  {
195  if (i + GetLastPyramidLevel() >= m_ReferencePyramid->GetNumberOfLevels())
196  continue;
197 
198  typename InputImageType::Pointer refImage = m_ReferencePyramid->GetOutput(i);
199  refImage->DisconnectPipeline();
200 
201  typename InputImageType::Pointer floImage = m_FloatingPyramid->GetOutput(i);
202  floImage->DisconnectPipeline();
203 
204  typename MaskImageType::Pointer maskGenerationImage = ITK_NULLPTR;
205  if (m_BlockGenerationPyramid)
206  {
207  maskGenerationImage = m_BlockGenerationPyramid->GetOutput(i);
208  maskGenerationImage->DisconnectPipeline();
209  }
210 
211  BlockMatcherType *mainMatcher = new BlockMatcherType;
212  BlockMatcherType *reverseMatcher = 0;
213  mainMatcher->SetBlockPercentageKept(GetPercentageKept());
214  mainMatcher->SetBlockSize(GetBlockSize());
215  mainMatcher->SetBlockSpacing(GetBlockSpacing());
216  mainMatcher->SetBlockVarianceThreshold(GetStDevThreshold() * GetStDevThreshold());
217  mainMatcher->SetBlockGenerationMask(maskGenerationImage);
218  mainMatcher->SetDefaultBackgroundValue(m_FloatingMinimalValue);
219 
220  if (m_Verbose)
221  {
222  std::cout << "Processing pyramid level " << i << std::endl;
223  std::cout << "Image size: " << refImage->GetLargestPossibleRegion().GetSize() << std::endl;
224  }
225 
226  // Init bm registration method
227  switch (m_SymmetryType)
228  {
229  case Asymmetric:
230  {
231  typedef typename anima::AsymmetricBMRegistrationMethod <InputImageType> BlockMatchRegistrationType;
232  m_bmreg = BlockMatchRegistrationType::New();
233  break;
234  }
235 
236  case Symmetric:
237  {
238  typedef typename anima::SymmetricBMRegistrationMethod <InputImageType> BlockMatchRegistrationType;
239  typename BlockMatchRegistrationType::Pointer tmpReg = BlockMatchRegistrationType::New();
240 
241  reverseMatcher = new BlockMatcherType;
242  reverseMatcher->SetBlockPercentageKept(GetPercentageKept());
243  reverseMatcher->SetBlockSize(GetBlockSize());
244  reverseMatcher->SetBlockSpacing(GetBlockSpacing());
245  reverseMatcher->SetBlockVarianceThreshold(GetStDevThreshold() * GetStDevThreshold());
246  reverseMatcher->SetVerbose(m_Verbose);
247  reverseMatcher->SetBlockGenerationMask(maskGenerationImage);
248  reverseMatcher->SetDefaultBackgroundValue(m_ReferenceMinimalValue);
249 
250  tmpReg->SetReverseBlockMatcher(reverseMatcher);
251  m_bmreg = tmpReg;
252  break;
253  }
254 
255  case Kissing:
256  default:
257  {
258  typedef typename anima::KissingSymmetricBMRegistrationMethod <InputImageType> BlockMatchRegistrationType;
259  typename BlockMatchRegistrationType::Pointer tmpReg = BlockMatchRegistrationType::New();
260  tmpReg->SetReferenceBackgroundValue(m_ReferenceMinimalValue);
261  tmpReg->SetFloatingBackgroundValue(m_FloatingMinimalValue);
262 
263  m_bmreg = tmpReg;
264  break;
265  }
266  }
267 
268  mainMatcher->SetVerbose(m_Verbose);
269  m_bmreg->SetBlockMatcher(mainMatcher);
270 
271  if (m_progressCallback)
272  {
273  // we cannot connect directly bmreg to m_progressCallback
274  // we need to create a new progressReporter with more iterations (m_progressReporter),
275  // to listen to progress events from bmreg and to send new ones to m_progressCallback
276  m_bmreg->AddObserver(itk::ProgressEvent(), m_callback);
277  }
278 
279  if (this->GetNumberOfWorkUnits() != 0)
280  m_bmreg->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
281 
282  m_bmreg->SetFixedImage(refImage);
283  m_bmreg->SetMovingImage(floImage);
284 
285  typedef anima::ResampleImageFilter<InputImageType, InputImageType,
286  typename AgregatorType::ScalarType> ResampleFilterType;
287 
288  typename ResampleFilterType::Pointer refResampler = ResampleFilterType::New();
289  refResampler->SetSize(floImage->GetLargestPossibleRegion().GetSize());
290  refResampler->SetOutputOrigin(floImage->GetOrigin());
291  refResampler->SetOutputSpacing(floImage->GetSpacing());
292  refResampler->SetOutputDirection(floImage->GetDirection());
293  refResampler->SetDefaultPixelValue(m_ReferenceMinimalValue);
294  refResampler->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
295  m_bmreg->SetReferenceImageResampler(refResampler);
296 
297  typename ResampleFilterType::Pointer movingResampler = ResampleFilterType::New();
298  movingResampler->SetSize(refImage->GetLargestPossibleRegion().GetSize());
299  movingResampler->SetOutputOrigin(refImage->GetOrigin());
300  movingResampler->SetOutputSpacing(refImage->GetSpacing());
301  movingResampler->SetOutputDirection(refImage->GetDirection());
302  movingResampler->SetDefaultPixelValue(m_FloatingMinimalValue);
303  movingResampler->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
304  m_bmreg->SetMovingImageResampler(movingResampler);
305 
306  BaseAgreg *agreg = NULL;
307 
308  switch (GetAgregator())
309  {
310  case LeastSquares:
311  {
312  typedef LSWTransformAgregator<ImageDimension> Agreg;
313  agreg = new Agreg;
314  break;
315  }
316 
317  case LeastTrimmedSquares:
318  {
319  typedef LTSWTransformAgregator<ImageDimension> Agreg;
320  Agreg *tmpAg = new Agreg;
321 
322  tmpAg->SetLTSCut(GetAgregThreshold());
323  tmpAg->SeStoppingThreshold(GetSeStoppingThreshold());
324 
325  agreg = tmpAg;
326  break;
327  }
328 
329  case MEstimation:
330  {
331  typedef MEstTransformAgregator<ImageDimension> Agreg;
332  Agreg *tmpAg = new Agreg;
333 
334  tmpAg->SetMEstimateFactor(GetAgregThreshold());
335  tmpAg->SeStoppingThreshold(GetSeStoppingThreshold());
336 
337  agreg = tmpAg;
338  break;
339  }
340  }
341 
342  switch (GetOutputTransformType())
343  {
344  case outTranslation:
345  agreg->SetOutputTransformType(BaseAgreg::TRANSLATION);
346  break;
347  case outRigid:
348  agreg->SetOutputTransformType(BaseAgreg::RIGID);
349  break;
350  case outAnisotropic_Sim:
351  agreg->SetOutputTransformType(BaseAgreg::ANISOTROPIC_SIM);
352  if (m_DirectionTransform)
353  agreg->SetOrthogonalDirectionMatrix(m_DirectionTransform->GetMatrix());
354  break;
355  case outAffine:
356  default:
357  agreg->SetOutputTransformType(BaseAgreg::AFFINE);
358  break;
359  }
360 
361  agreg->SetVerboseAgregation(m_Verbose);
362  m_bmreg->SetAgregator(agreg);
363 
364  switch (GetTransform())
365  {
366  case Translation:
367  mainMatcher->SetBlockTransformType(BlockMatcherType::Superclass::Translation);
368  if (reverseMatcher)
369  reverseMatcher->SetBlockTransformType(BlockMatcherType::Superclass::Translation);
370  break;
371  case Rigid:
372  mainMatcher->SetBlockTransformType(BlockMatcherType::Superclass::Rigid);
373  if (reverseMatcher)
374  reverseMatcher->SetBlockTransformType(BlockMatcherType::Superclass::Rigid);
375  break;
376  case Directional_Affine:
377  mainMatcher->SetBlockTransformType(BlockMatcherType::Superclass::Directional_Affine);
378  mainMatcher->SetAffineDirection(m_AffineDirection);
379  if (reverseMatcher)
380  {
381  reverseMatcher->SetBlockTransformType(BlockMatcherType::Superclass::Directional_Affine);
382  reverseMatcher->SetAffineDirection(m_AffineDirection);
383  }
384  break;
385  case Affine:
386  default:
387  mainMatcher->SetBlockTransformType(BlockMatcherType::Superclass::Affine);
388  if (reverseMatcher)
389  reverseMatcher->SetBlockTransformType(BlockMatcherType::Superclass::Affine);
390  break;
391  }
392 
393 
394  switch(GetOptimizer())
395  {
396  case Exhaustive:
397  mainMatcher->SetOptimizerType(BlockMatcherType::Exhaustive);
398  if (reverseMatcher)
399  reverseMatcher->SetOptimizerType(BlockMatcherType::Exhaustive);
400  break;
401 
402  case Bobyqa:
403  default:
404  mainMatcher->SetOptimizerType(BlockMatcherType::Bobyqa);
405  if (reverseMatcher)
406  reverseMatcher->SetOptimizerType(BlockMatcherType::Bobyqa);
407  break;
408  }
409 
410  switch (GetMetric())
411  {
412  case Correlation:
413  mainMatcher->SetSimilarityType(BlockMatcherType::Correlation);
414  if (reverseMatcher)
415  reverseMatcher->SetSimilarityType(BlockMatcherType::Correlation);
416  break;
417  case SquaredCorrelation:
418  mainMatcher->SetSimilarityType(BlockMatcherType::SquaredCorrelation);
419  if (reverseMatcher)
420  reverseMatcher->SetSimilarityType(BlockMatcherType::SquaredCorrelation);
421  break;
422  case MeanSquares:
423  default:
424  mainMatcher->SetSimilarityType(BlockMatcherType::MeanSquares);
425  if (reverseMatcher)
426  reverseMatcher->SetSimilarityType(BlockMatcherType::MeanSquares);
427  break;
428  }
429 
430  m_bmreg->SetMaximumIterations(GetMaximumIterations());
431  m_bmreg->SetMinimalTransformError(GetMinimalTransformError());
432  m_bmreg->SetInitialTransform(m_OutputTransform);
433 
434  mainMatcher->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
435  mainMatcher->SetOptimizerMaximumIterations(GetOptimizerMaximumIterations());
436 
437  double sr = GetSearchRadius();
438  mainMatcher->SetSearchRadius(sr);
439 
440  double sar = GetSearchAngleRadius();
441  mainMatcher->SetSearchAngleRadius(sar);
442 
443  double scr = GetSearchScaleRadius();
444  mainMatcher->SetSearchScaleRadius(scr);
445 
446  double fr = GetFinalRadius();
447  mainMatcher->SetFinalRadius(fr);
448 
449  double ss = GetStepSize();
450  mainMatcher->SetStepSize(ss);
451 
452  double tub = GetTranslateUpperBound();
453  mainMatcher->SetTranslateMax(tub);
454 
455  double aub = GetAngleUpperBound();
456  mainMatcher->SetAngleMax(aub);
457 
458  double scub = GetScaleUpperBound();
459  mainMatcher->SetScaleMax(scub);
460 
461  if (reverseMatcher)
462  {
463  reverseMatcher->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
464  reverseMatcher->SetOptimizerMaximumIterations(GetOptimizerMaximumIterations());
465 
466  reverseMatcher->SetSearchRadius(sr);
467  reverseMatcher->SetSearchAngleRadius(sar);
468  reverseMatcher->SetSearchScaleRadius(scr);
469  reverseMatcher->SetFinalRadius(fr);
470  reverseMatcher->SetStepSize(ss);
471  reverseMatcher->SetTranslateMax(tub);
472  reverseMatcher->SetAngleMax(aub);
473  reverseMatcher->SetScaleMax(scub);
474  }
475 
476  m_bmreg->SetVerboseProgression(m_Verbose);
477 
478  try
479  {
480  m_bmreg->Update();
481  }
482  catch( itk::ExceptionObject & err )
483  {
484  std::cerr << "ExceptionObject caught in bmreg startregistration ! " << err << std::endl;
485  exit(-1);
486  }
487 
488  if ((GetOutputTransformType() == outAnisotropic_Sim)||(GetOutputTransformType() == outAffine))
489  m_EstimationBarycenter = agreg->GetEstimationBarycenter();
490 
491  // Polyrigid will have to be handled here
492  AffineTransformType *tmpTrsf = dynamic_cast<AffineTransformType *>(m_OutputTransform.GetPointer());
493  tmpTrsf->SetParameters(m_bmreg->GetOutput()->Get()->GetParameters());
494 
495  delete mainMatcher;
496  if (reverseMatcher)
497  delete reverseMatcher;
498  if (agreg)
499  delete agreg;
500  }
501 
502  if (m_Abort)
503  std::cout << "Process aborted" << std::endl;
504 
505  this->InvokeEvent(itk::EndEvent());
506  this->RemoveAllObservers();
507 
508  AffineTransformType *tmpTrsf = dynamic_cast<AffineTransformType *>(m_OutputTransform.GetPointer());
509 
510  if (m_SymmetryType == Kissing)
511  {
512  typename AffineTransformType::Pointer trsfCopy = AffineTransformType::New();
513  trsfCopy->SetMatrix(tmpTrsf->GetMatrix());
514  trsfCopy->SetOffset(tmpTrsf->GetOffset());
515 
516  tmpTrsf->Compose(trsfCopy);
517  }
518 
519  if (!m_InitialTransform.IsNull())
520  tmpTrsf->Compose(m_InitialTransform, false);
521 
522  typedef typename anima::ResampleImageFilter<InputImageType, InputImageType, typename AgregatorType::ScalarType> ResampleFilterType;
523  typename ResampleFilterType::Pointer tmpResample = ResampleFilterType::New();
524  tmpResample->SetTransform(m_OutputTransform);
525  tmpResample->SetInput(m_FloatingImage);
526 
527  tmpResample->SetSize(m_ReferenceImage->GetLargestPossibleRegion().GetSize());
528  tmpResample->SetOutputOrigin(m_ReferenceImage->GetOrigin());
529  tmpResample->SetOutputSpacing(m_ReferenceImage->GetSpacing());
530  tmpResample->SetOutputDirection(m_ReferenceImage->GetDirection());
531  tmpResample->SetDefaultPixelValue(m_FloatingMinimalValue);
532  tmpResample->Update();
533 
534  m_OutputImage = tmpResample->GetOutput();
535  m_OutputImage->DisconnectPipeline();
536 }
537 
538 template <unsigned int ImageDimension>
539 void PyramidalBlockMatchingBridge<ImageDimension>::EmitProgress(int prog)
540 {
541  if (m_progressReporter)
542  m_progressReporter->CompletedPixel();
543 }
544 
545 template <unsigned int ImageDimension>
546 void PyramidalBlockMatchingBridge<ImageDimension>::ManageProgress (itk::Object* caller, const itk::EventObject& event, void* clientData)
547 {
548  PyramidalBlockMatchingBridge * source = reinterpret_cast<PyramidalBlockMatchingBridge *> (clientData);
549  itk::ProcessObject *processObject = (itk::ProcessObject *) caller;
550 
551  if (source && processObject)
552  source->EmitProgress(processObject->GetProgress() * 100);
553 }
554 
555 template <unsigned int ImageDimension>
556 void PyramidalBlockMatchingBridge<ImageDimension>::WriteOutputs()
557 {
558  std::cout << "Writing output image to: " << GetResultFile() << std::endl;
559  anima::writeImage <InputImageType> (GetResultFile(),m_OutputImage);
560 
561  if (GetOutputTransformFile() != "")
562  {
563  std::cout << "Writing output transform to: " << GetOutputTransformFile() << std::endl;
564  itk::TransformFileWriter::Pointer writer = itk::TransformFileWriter::New();
565  writer->SetInput(m_OutputTransform);
566  writer->SetFileName(GetOutputTransformFile());
567  writer->Update();
568  }
569 
570  if ((m_outputNearestRigidTransformFile == "")&&(m_outputNearestSimilarityTransformFile == ""))
571  return;
572 
573  AffineTransformType *tmpTrsf = dynamic_cast<AffineTransformType *>(m_OutputTransform.GetPointer());
574 
575  typename AffineTransformType::MatrixType linearMatrix = tmpTrsf->GetMatrix();
576  typename AffineTransformType::OffsetType transformOffset = tmpTrsf->GetOffset();
577  vnl_svd<typename AffineTransformType::MatrixType::ValueType> UWVLinearMatrixSVD(linearMatrix.GetVnlMatrix().as_matrix());
578  vnl_matrix<typename AffineTransformType::MatrixType::ValueType> leftRot = UWVLinearMatrixSVD.U()*vnl_determinant(UWVLinearMatrixSVD.U());
579  vnl_matrix<typename AffineTransformType::MatrixType::ValueType> rightRot = UWVLinearMatrixSVD.V()*vnl_determinant(UWVLinearMatrixSVD.V());
580  vnl_diag_matrix<typename AffineTransformType::MatrixType::ValueType> scal = UWVLinearMatrixSVD.W();
581 
582  PointType ybar;
583  double isoScal = 0;
584  for (unsigned int i = 0;i < ImageDimension;++i)
585  {
586  isoScal += std::abs(scal(i, i)) / ImageDimension;
587  ybar[i] = transformOffset[i];
588 
589  for (unsigned int j = 0;j < ImageDimension;++j)
590  ybar[i] += linearMatrix(i, j) * m_EstimationBarycenter[j];
591  }
592 
593  typename AffineTransformType::OffsetType rigidOffset;
594  typename AffineTransformType::MatrixType linearPartRigid;
595  typename AffineTransformType::OffsetType similarityOffset;
596  typename AffineTransformType::MatrixType linearPartSimilarity;
597  typename AffineTransformType::OffsetType UOffset;
598  linearPartRigid.Fill(0);
599  linearPartSimilarity.Fill(0);
600 
601  for (unsigned int i = 0;i < ImageDimension;++i)
602  {
603  rigidOffset[i] = ybar[i];
604  similarityOffset[i] = ybar[i];
605  UOffset[i] = 0;
606  for (unsigned int j = 0;j < ImageDimension; ++j)
607  {
608  for (unsigned int k = 0;k < ImageDimension;++k)
609  {
610  linearPartRigid(i, j) += leftRot(i, k) * rightRot(j, k);
611  linearPartSimilarity(i, j) += isoScal * leftRot(i, k) * rightRot(j, k);
612  }
613 
614  rigidOffset[i] -= linearPartRigid(i, j) * m_EstimationBarycenter[j];
615  similarityOffset[i] -= linearPartSimilarity(i, j) * m_EstimationBarycenter[j];
616  }
617  }
618 
619  if (m_outputNearestRigidTransformFile != "")
620  {
621  AffineTransformPointer rigidTransform = AffineTransformType::New();
622  rigidTransform->SetMatrix(linearPartRigid);
623  rigidTransform->SetOffset(rigidOffset);
624  itk::TransformFileWriter::Pointer rigidWriter = itk::TransformFileWriter::New();
625  rigidWriter->SetInput(rigidTransform);
626  rigidWriter->SetFileName(GetOutputNearestRigidTransformFile());
627  std::cout << "Writing output nearest rigid transform to: " << GetOutputNearestRigidTransformFile() << std::endl;
628  rigidWriter->Update();
629  }
630 
631  if (m_outputNearestSimilarityTransformFile != "")
632  {
633  AffineTransformPointer similarityTransform = AffineTransformType::New();
634  similarityTransform->SetMatrix(linearPartSimilarity);
635  similarityTransform->SetOffset(similarityOffset);
636  itk::TransformFileWriter::Pointer similarityWriter = itk::TransformFileWriter::New();
637  similarityWriter->SetInput(similarityTransform);
638  similarityWriter->SetFileName(GetOutputNearestSimilarityTransformFile());
639  std::cout << "Writing output nearest similarity transform to: " << GetOutputNearestSimilarityTransformFile() << std::endl;
640  similarityWriter->Update();
641  }
642 }
643 
644 template <unsigned int ImageDimension>
645 void PyramidalBlockMatchingBridge<ImageDimension>::SetupPyramids()
646 {
647  // Create pyramid here, check images actually are of the same size.
648  typedef anima::ResampleImageFilter<InputImageType, InputImageType,
649  typename AgregatorType::ScalarType> ResampleFilterType;
650 
651  m_ReferencePyramid = PyramidType::New();
652 
653  m_ReferencePyramid->SetInput(m_ReferenceImage);
654  m_ReferencePyramid->SetNumberOfLevels(GetNumberOfPyramidLevels());
655  m_ReferencePyramid->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
656 
657  typename ResampleFilterType::Pointer refResampler = ResampleFilterType::New();
658  refResampler->SetDefaultPixelValue(m_ReferenceMinimalValue);
659  m_ReferencePyramid->SetImageResampler(refResampler);
660  m_ReferencePyramid->Update();
661 
662  InputImagePointer initialFloatingImage = const_cast <InputImageType *> (m_FloatingImage.GetPointer());
663 
664  // Compute initial transform if needed to get a decent initial Floating image
665  if (m_InitialTransform.IsNotNull())
666  {
667  typename ResampleFilterType::Pointer tmpResample = ResampleFilterType::New();
668  tmpResample->SetTransform(m_InitialTransform);
669  tmpResample->SetInput(m_FloatingImage);
670 
671  tmpResample->SetSize(m_ReferenceImage->GetLargestPossibleRegion().GetSize());
672  tmpResample->SetOutputOrigin(m_ReferenceImage->GetOrigin());
673  tmpResample->SetOutputSpacing(m_ReferenceImage->GetSpacing());
674  tmpResample->SetOutputDirection(m_ReferenceImage->GetDirection());
675  tmpResample->SetDefaultPixelValue(m_FloatingMinimalValue);
676  tmpResample->Update();
677 
678  initialFloatingImage = tmpResample->GetOutput();
679  initialFloatingImage->DisconnectPipeline();
680  }
681  else
682  {
683  m_InitialTransform = NULL;
684 
685  m_InitialTransform = AffineTransformType::New();
686  m_InitialTransform->SetIdentity();
687 
688  typedef itk::ImageMomentsCalculator <InputImageType> ImageCalculatorType;
689 
690  if (m_TransformInitializationType != Identity)
691  {
692  typename ImageCalculatorType::Pointer fixedCalculator = ImageCalculatorType::New();
693  fixedCalculator->SetImage(m_ReferenceImage);
694  fixedCalculator->Compute();
695  typename ImageCalculatorType::VectorType fixedBar = fixedCalculator->GetCenterOfGravity();
696 
697  typename ImageCalculatorType::Pointer movingCalculator = ImageCalculatorType::New();
698  movingCalculator->SetImage(m_FloatingImage);
699  movingCalculator->Compute();
700  typename ImageCalculatorType::VectorType movingBar = movingCalculator->GetCenterOfGravity();
701 
702  m_InitialTransform->SetOffset(movingBar - fixedBar);
703 
704  if ((m_TransformInitializationType == ClosestTransform)&&(m_OutputTransformType != outTranslation))
705  {
706  typename ImageCalculatorType::VectorType fixedPrincipalMom = fixedCalculator->GetPrincipalMoments();
707  typename ImageCalculatorType::ScalarType fixedMass = fixedCalculator->GetTotalMass();
708 
709  typename ImageCalculatorType::VectorType movingPrincipalMom = movingCalculator->GetPrincipalMoments();
710  typename ImageCalculatorType::ScalarType movingMass = movingCalculator->GetTotalMass();
711 
712  vnl_matrix<double> scalMatrix(ImageDimension, ImageDimension, 0);
713  itk::Vector<double, ImageDimension> scalOffset;
714  double fixedScalFactor = 0;
715  double movingScalFactor = 0;
716 
717  for (unsigned int i = 0; i < ImageDimension; ++i)
718  {
719  fixedScalFactor += pow(fixedPrincipalMom[i] / fixedMass, 0.5);
720  movingScalFactor += pow(movingPrincipalMom[i] / movingMass, 0.5);
721  }
722 
723  double scalingFactor = 1.0;
724  if ((m_OutputTransformType == outAnisotropic_Sim)||(m_OutputTransformType == outAffine))
725  scalingFactor = movingScalFactor / fixedScalFactor;
726 
727  scalMatrix.fill_diagonal(scalingFactor);
728 
729  for (unsigned int i = 0; i < ImageDimension; ++i)
730  scalOffset[i] = movingBar[i] - scalingFactor * fixedBar[i];
731 
732  m_InitialTransform->SetMatrix(scalMatrix);
733  m_InitialTransform->SetOffset(scalOffset);
734  }
735  }
736 
737  typename ResampleFilterType::Pointer tmpResample = ResampleFilterType::New();
738  tmpResample->SetTransform(m_InitialTransform);
739  tmpResample->SetInput(m_FloatingImage);
740 
741  tmpResample->SetSize(m_ReferenceImage->GetLargestPossibleRegion().GetSize());
742  tmpResample->SetOutputOrigin(m_ReferenceImage->GetOrigin());
743  tmpResample->SetOutputSpacing(m_ReferenceImage->GetSpacing());
744  tmpResample->SetOutputDirection(m_ReferenceImage->GetDirection());
745  tmpResample->SetDefaultPixelValue(m_FloatingMinimalValue);
746  tmpResample->Update();
747 
748  initialFloatingImage = tmpResample->GetOutput();
749  initialFloatingImage->DisconnectPipeline();
750 
751  using MinMaxImageFilterType = itk::MinimumMaximumImageFilter <InputImageType>;
752  typename MinMaxImageFilterType::Pointer minMaxFilter = MinMaxImageFilterType::New();
753  minMaxFilter->SetInput(initialFloatingImage);
754  minMaxFilter->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
755  minMaxFilter->Update();
756 
757  if ((minMaxFilter->GetMinimum() == minMaxFilter->GetMaximum())&&(m_TransformInitializationType == Identity))
758  {
759  std::cout << "Identity initialization outputs an empty image, initializing with centers of mass" << std::endl;
760  m_TransformInitializationType = GravityCenters;
761 
762  typename ImageCalculatorType::Pointer fixedCalculator = ImageCalculatorType::New();
763  fixedCalculator->SetImage(m_ReferenceImage);
764  fixedCalculator->Compute();
765  typename ImageCalculatorType::VectorType fixedBar = fixedCalculator->GetCenterOfGravity();
766 
767  typename ImageCalculatorType::Pointer movingCalculator = ImageCalculatorType::New();
768  movingCalculator->SetImage(m_FloatingImage);
769  movingCalculator->Compute();
770  typename ImageCalculatorType::VectorType movingBar = movingCalculator->GetCenterOfGravity();
771 
772  m_InitialTransform->SetOffset(movingBar - fixedBar);
773 
774  typename ResampleFilterType::Pointer tmpResample = ResampleFilterType::New();
775  tmpResample->SetTransform(m_InitialTransform);
776  tmpResample->SetInput(m_FloatingImage);
777 
778  tmpResample->SetSize(m_ReferenceImage->GetLargestPossibleRegion().GetSize());
779  tmpResample->SetOutputOrigin(m_ReferenceImage->GetOrigin());
780  tmpResample->SetOutputSpacing(m_ReferenceImage->GetSpacing());
781  tmpResample->SetOutputDirection(m_ReferenceImage->GetDirection());
782  tmpResample->SetDefaultPixelValue(m_FloatingMinimalValue);
783  tmpResample->Update();
784 
785  initialFloatingImage = tmpResample->GetOutput();
786  initialFloatingImage->DisconnectPipeline();
787  }
788  }
789 
790  // Create pyramid for Floating image
791  m_FloatingPyramid = PyramidType::New();
792 
793  m_FloatingPyramid->SetInput(initialFloatingImage);
794  m_FloatingPyramid->SetNumberOfLevels(GetNumberOfPyramidLevels());
795  m_FloatingPyramid->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
796 
797  typename ResampleFilterType::Pointer floResampler = ResampleFilterType::New();
798  floResampler->SetDefaultPixelValue(m_FloatingMinimalValue);
799  m_FloatingPyramid->SetImageResampler(floResampler);
800 
801  m_FloatingPyramid->Update();
802 
803  m_BlockGenerationPyramid = 0;
804  if (m_BlockGenerationMask)
805  {
806  typedef anima::ResampleImageFilter<MaskImageType, MaskImageType,
807  typename AgregatorType::ScalarType> MaskResampleFilterType;
808 
809  typename MaskResampleFilterType::Pointer maskResampler = MaskResampleFilterType::New();
810 
811  m_BlockGenerationPyramid = MaskPyramidType::New();
812  m_BlockGenerationPyramid->SetImageResampler(maskResampler);
813  m_BlockGenerationPyramid->SetInput(m_BlockGenerationMask);
814  m_BlockGenerationPyramid->SetNumberOfLevels(GetNumberOfPyramidLevels());
815  m_BlockGenerationPyramid->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
816  m_BlockGenerationPyramid->Update();
817  }
818 }
819 
820 } // end of namespace anima
anima::MEstTransformAgregator::SetMEstimateFactor
void SetMEstimateFactor(double mestFactor)
Definition: animaMEstTransformAgregator.h:26
anima::PyramidalBlockMatchingBridge::SetDirectionTransform
void SetDirectionTransform(AffineTransformPointer directionTransform)
Definition: animaPyramidalBlockMatchingBridge.h:121
anima::PyramidalBlockMatchingBridge::ManageProgress
static void ManageProgress(itk::Object *caller, const itk::EventObject &event, void *clientData)
Definition: animaPyramidalBlockMatchingBridge.hxx:546
anima::PyramidalBlockMatchingBridge::AffineTransformType
itk::AffineTransform< typename AgregatorType::ScalarType, ImageDimension > AffineTransformType
Definition: animaPyramidalBlockMatchingBridge.h:34
anima::BaseBlockMatcher::SetBlockPercentageKept
void SetBlockPercentageKept(double val)
Definition: animaBaseBlockMatcher.h:59
anima::PyramidalBlockMatchingBridge::MaskImageType
itk::Image< unsigned char, ImageDimension > MaskImageType
Definition: animaPyramidalBlockMatchingBridge.h:25
anima::PyramidalBlockMatchingBridge::InputImageType
itk::Image< double, ImageDimension > InputImageType
Definition: animaPyramidalBlockMatchingBridge.h:18
anima::PyramidalBlockMatchingBridge::SetInitialTransform
void SetInitialTransform(AffineTransformPointer initialTransform)
Definition: animaPyramidalBlockMatchingBridge.h:118
Generated by   doxygen 1.8.13