9 #include <itkResampleImageFilter.h> 21 template <
unsigned int ImageDimension>
24 m_ReferenceImage = NULL;
25 m_FloatingImage = NULL;
27 m_OutputTransform = BaseTransformType::New();
28 m_OutputTransform->SetIdentity();
30 m_outputTransformFile =
"";
40 m_FiniteStrainImageReorientation =
true;
45 m_MaximumIterations = 10;
46 m_MinimalTransformError = 0.01;
47 m_OptimizerMaximumIterations = 100;
49 m_SearchAngleRadius = 5;
50 m_SearchScaleRadius = 0.1;
51 m_FinalRadius = 0.001;
53 m_TranslateUpperBound = 50;
54 m_AngleUpperBound = 180;
55 m_ScaleUpperBound = 3;
57 m_ExtrapolationSigma = 3;
60 m_MEstimateConvergenceThreshold = 0.01;
61 m_NeighborhoodApproximation = 2.5;
62 m_BCHCompositionOrder = 1;
63 m_ExponentiationOrder = 1;
64 m_NumberOfPyramidLevels = 3;
65 m_LastPyramidLevel = 0;
66 m_PercentageKept = 0.8;
67 this->SetNumberOfWorkUnits(itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads());
70 template <
unsigned int ImageDimension>
75 template <
unsigned int ImageDimension>
79 this->SetupPyramids();
82 for (
unsigned int i = 0;i < m_ReferencePyramid->GetNumberOfLevels();++i)
84 if (i + m_LastPyramidLevel >= m_ReferencePyramid->GetNumberOfLevels())
87 typename InputImageType::Pointer refImage = m_ReferencePyramid->GetOutput(i);
88 refImage->DisconnectPipeline();
90 typename InputImageType::Pointer floImage = m_FloatingPyramid->GetOutput(i);
91 floImage->DisconnectPipeline();
93 typename MaskImageType::Pointer maskGenerationImage = ITK_NULLPTR;
94 if (m_BlockGenerationPyramid)
96 maskGenerationImage = m_BlockGenerationPyramid->GetOutput(i);
97 maskGenerationImage->DisconnectPipeline();
101 if (m_OutputTransform->GetParametersAsVectorField() != NULL)
103 typedef itk::ResampleImageFilter<VelocityFieldType,VelocityFieldType> VectorResampleFilterType;
104 typedef typename VectorResampleFilterType::Pointer VectorResampleFilterPointer;
107 tmpIdentity->SetIdentity();
109 VectorResampleFilterPointer tmpResample = VectorResampleFilterType::New();
110 tmpResample->SetTransform(tmpIdentity);
111 tmpResample->SetInput(m_OutputTransform->GetParametersAsVectorField());
113 tmpResample->SetSize(refImage->GetLargestPossibleRegion().GetSize());
114 tmpResample->SetOutputOrigin(refImage->GetOrigin());
115 tmpResample->SetOutputSpacing(refImage->GetSpacing());
116 tmpResample->SetOutputDirection(refImage->GetDirection());
118 tmpResample->Update();
121 m_OutputTransform->SetParametersAsVectorField(tmpOut);
122 tmpOut->DisconnectPipeline();
125 std::cout <<
"Processing pyramid level " << i << std::endl;
126 std::cout <<
"Image size: " << refImage->GetLargestPossibleRegion().GetSize() << std::endl;
128 double meanSpacing = 0;
129 for (
unsigned int j = 0;j < ImageDimension;++j)
130 meanSpacing += refImage->GetSpacing()[j];
131 meanSpacing /= ImageDimension;
143 if (this->GetNumberOfWorkUnits() != 0)
149 agreg->
SetDistanceBoundary(m_ExtrapolationSigma * meanSpacing * m_NeighborhoodApproximation);
161 if (this->GetNumberOfWorkUnits() != 0)
172 BlockMatcherType *mainMatcher =
new BlockMatcherType;
173 BlockMatcherType *reverseMatcher = 0;
175 mainMatcher->SetBlockSize(GetBlockSize());
176 mainMatcher->SetBlockSpacing(GetBlockSpacing());
177 mainMatcher->SetBlockVarianceThreshold(GetStDevThreshold() * GetStDevThreshold());
178 mainMatcher->SetBlockGenerationMask(maskGenerationImage);
180 switch (m_SymmetryType)
185 m_bmreg = BlockMatchRegistrationType::New();
192 typename BlockMatchRegistrationType::Pointer tmpReg = BlockMatchRegistrationType::New();
193 reverseMatcher =
new BlockMatcherType;
194 reverseMatcher->SetBlockPercentageKept(GetPercentageKept());
195 reverseMatcher->SetBlockSize(GetBlockSize());
196 reverseMatcher->SetBlockSpacing(GetBlockSpacing());
197 reverseMatcher->SetBlockVarianceThreshold(GetStDevThreshold() * GetStDevThreshold());
198 reverseMatcher->SetBlockGenerationMask(maskGenerationImage);
200 tmpReg->SetReverseBlockMatcher(reverseMatcher);
208 m_bmreg = BlockMatchRegistrationType::New();
213 m_bmreg->SetBlockMatcher(mainMatcher);
214 m_bmreg->SetAgregator(agregPtr);
215 m_bmreg->SetBCHCompositionOrder(m_BCHCompositionOrder);
216 m_bmreg->SetExponentiationOrder(m_ExponentiationOrder);
218 if (this->GetNumberOfWorkUnits() != 0)
219 m_bmreg->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
221 m_bmreg->SetFixedImage(refImage);
222 m_bmreg->SetMovingImage(floImage);
224 m_bmreg->SetSVFElasticRegSigma(m_ElasticSigma * meanSpacing);
228 typename ResampleFilterType::Pointer refResampler = ResampleFilterType::New();
229 refResampler->SetOutputLargestPossibleRegion(floImage->GetLargestPossibleRegion());
230 refResampler->SetOutputOrigin(floImage->GetOrigin());
231 refResampler->SetOutputSpacing(floImage->GetSpacing());
232 refResampler->SetOutputDirection(floImage->GetDirection());
233 refResampler->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
234 refResampler->SetFiniteStrainReorientation(this->GetFiniteStrainImageReorientation());
235 m_bmreg->SetReferenceImageResampler(refResampler);
237 typename ResampleFilterType::Pointer movingResampler = ResampleFilterType::New();
238 movingResampler->SetOutputLargestPossibleRegion(refImage->GetLargestPossibleRegion());
239 movingResampler->SetOutputOrigin(refImage->GetOrigin());
240 movingResampler->SetOutputSpacing(refImage->GetSpacing());
241 movingResampler->SetOutputDirection(refImage->GetDirection());
242 movingResampler->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
243 movingResampler->SetFiniteStrainReorientation(this->GetFiniteStrainImageReorientation());
244 m_bmreg->SetMovingImageResampler(movingResampler);
246 switch (GetTransform())
266 switch (GetOptimizer())
307 switch (m_MetricOrientation)
327 m_bmreg->SetMaximumIterations(m_MaximumIterations);
328 m_bmreg->SetMinimalTransformError(m_MinimalTransformError);
329 m_bmreg->SetInitialTransform(m_OutputTransform.GetPointer());
331 mainMatcher->SetOptimizerMaximumIterations(m_OptimizerMaximumIterations);
333 double sr = m_SearchRadius;
334 mainMatcher->SetSearchRadius(sr);
336 double sar = m_SearchAngleRadius;
337 mainMatcher->SetSearchAngleRadius(sar);
339 double scr = m_SearchScaleRadius;
340 mainMatcher->SetSearchScaleRadius(scr);
342 double fr = m_FinalRadius;
343 mainMatcher->SetFinalRadius(fr);
345 double ss = m_StepSize;
346 mainMatcher->SetStepSize(ss);
348 double tub = m_TranslateUpperBound;
349 mainMatcher->SetTranslateMax(tub);
351 double aub = m_AngleUpperBound;
352 mainMatcher->SetAngleMax(aub);
354 double scub = m_ScaleUpperBound;
355 mainMatcher->SetScaleMax(scub);
359 reverseMatcher->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
360 reverseMatcher->SetOptimizerMaximumIterations(GetOptimizerMaximumIterations());
362 reverseMatcher->SetSearchRadius(sr);
363 reverseMatcher->SetSearchAngleRadius(sar);
364 reverseMatcher->SetSearchScaleRadius(scr);
365 reverseMatcher->SetFinalRadius(fr);
366 reverseMatcher->SetStepSize(ss);
367 reverseMatcher->SetTranslateMax(tub);
368 reverseMatcher->SetAngleMax(aub);
369 reverseMatcher->SetScaleMax(scub);
376 catch( itk::ExceptionObject & err )
378 std::cout <<
"ExceptionObject caught !" << err << std::endl;
383 m_OutputTransform->SetParametersAsVectorField(resTrsf->GetParametersAsVectorField());
387 delete reverseMatcher;
395 typedef itk::MultiplyImageFilter <VelocityFieldType,itk::Image <double,ImageDimension>,
VelocityFieldType> MultiplyFilterType;
397 typename MultiplyFilterType::Pointer fieldMultiplier = MultiplyFilterType::New();
398 fieldMultiplier->SetInput(finalTrsfField);
399 fieldMultiplier->SetConstant(2.0);
400 fieldMultiplier->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
401 fieldMultiplier->InPlaceOn();
403 fieldMultiplier->Update();
406 m_OutputTransform->SetParametersAsVectorField(fieldMultiplier->GetOutput());
407 outputField->DisconnectPipeline();
411 anima::GetSVFExponential(m_OutputTransform.GetPointer(), outputDispTrsf.GetPointer(), m_ExponentiationOrder, GetNumberOfWorkUnits(),
false);
414 typename ResampleFilterType::Pointer tmpResample = ResampleFilterType::New();
416 typedef itk::Transform<typename BaseAgregatorType::ScalarType,ImageDimension,ImageDimension>
BaseTransformType;
417 typename BaseTransformType::Pointer baseTrsf = outputDispTrsf.GetPointer();
418 tmpResample->SetTransform(baseTrsf);
419 tmpResample->SetFiniteStrainReorientation(this->GetFiniteStrainImageReorientation());
420 tmpResample->SetInput(m_FloatingImage);
422 if (this->GetNumberOfWorkUnits() != 0)
423 tmpResample->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
425 tmpResample->SetOutputLargestPossibleRegion(m_ReferenceImage->GetLargestPossibleRegion());
426 tmpResample->SetOutputOrigin(m_ReferenceImage->GetOrigin());
427 tmpResample->SetOutputSpacing(m_ReferenceImage->GetSpacing());
428 tmpResample->SetOutputDirection(m_ReferenceImage->GetDirection());
429 tmpResample->Update();
432 typename ExpTensorFilterType::Pointer tensorExper = ExpTensorFilterType::New();
434 typename InputImageType::Pointer tmpImage = tmpResample->GetOutput();
435 tmpImage->DisconnectPipeline();
437 tensorExper->SetInput(tmpImage);
438 tensorExper->SetScaleNonDiagonal(
true);
440 if (this->GetNumberOfWorkUnits() != 0)
441 tensorExper->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
443 tensorExper->Update();
445 m_OutputImage = tensorExper->GetOutput();
446 m_OutputImage->DisconnectPipeline();
449 template <
unsigned int ImageDimension>
453 std::cout <<
"Writing output image to: " << m_resultFile << std::endl;
454 anima::writeImage <InputImageType> (m_resultFile,m_OutputImage);
456 if (m_outputTransformFile !=
"")
458 std::cout <<
"Writing output SVF to: " << m_outputTransformFile << std::endl;
459 anima::writeImage <VelocityFieldType> (m_outputTransformFile,
460 const_cast <
VelocityFieldType *> (m_OutputTransform->GetParametersAsVectorField()));
464 template <
unsigned int ImageDimension>
469 m_ReferencePyramid = PyramidType::New();
471 m_ReferencePyramid->SetInput(m_ReferenceImage);
472 m_ReferencePyramid->SetNumberOfLevels(m_NumberOfPyramidLevels);
474 if (this->GetNumberOfWorkUnits() != 0)
475 m_ReferencePyramid->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
479 typename ResampleFilterType::Pointer refResampler = ResampleFilterType::New();
480 refResampler->SetFiniteStrainReorientation(this->GetFiniteStrainImageReorientation());
481 m_ReferencePyramid->SetImageResampler(refResampler);
483 m_ReferencePyramid->Update();
486 m_FloatingPyramid = PyramidType::New();
488 m_FloatingPyramid->SetInput(m_FloatingImage);
489 m_FloatingPyramid->SetNumberOfLevels(m_NumberOfPyramidLevels);
491 if (this->GetNumberOfWorkUnits() != 0)
492 m_FloatingPyramid->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
494 typename ResampleFilterType::Pointer floResampler = ResampleFilterType::New();
495 floResampler->SetFiniteStrainReorientation(this->GetFiniteStrainImageReorientation());
496 m_FloatingPyramid->SetImageResampler(floResampler);
498 m_FloatingPyramid->Update();
500 m_BlockGenerationPyramid = 0;
501 if (m_BlockGenerationMask)
506 typename MaskResampleFilterType::Pointer maskResampler = MaskResampleFilterType::New();
508 m_BlockGenerationPyramid = MaskPyramidType::New();
509 m_BlockGenerationPyramid->SetImageResampler(maskResampler);
510 m_BlockGenerationPyramid->SetInput(m_BlockGenerationMask);
511 m_BlockGenerationPyramid->SetNumberOfLevels(GetNumberOfPyramidLevels());
512 m_BlockGenerationPyramid->SetNumberOfWorkUnits(GetNumberOfWorkUnits());
513 m_BlockGenerationPyramid->Update();
PyramidalDenseTensorSVFMatchingBridge()
itk::Image< unsigned char, ImageDimension > MaskImageType
virtual ~PyramidalDenseTensorSVFMatchingBridge()
DisplacementFieldTransformType::Pointer DisplacementFieldTransformPointer
MEstimateAgregatorType::BaseOutputTransformType BaseTransformType
void Update() ITK_OVERRIDE
void SetBlockPercentageKept(double val)
void GetSVFExponential(itk::StationaryVelocityFieldTransform< ScalarType, NDimensions > *baseTrsf, rpi::DisplacementFieldTransform< ScalarType, NDimensions > *resultTransform, unsigned int exponentiationOrder, unsigned int numThreads, bool invert)
AffineTransformType::Pointer AffineTransformPointer
BaseTransformType::VectorFieldType VelocityFieldType