[ITK] Bspline registration with Mutual Information v4 : memory concern

[Timothee Evain]

Hello everyone,

I'm trying a Bspline registration based on the Mutual Information metric (v4) and I'm facing a memory problem.
When running the registration, the memory usage skyrockets in a second until freezing the PC by lack of memory. 
I don't know why : this doesn't happen when I'm switching to another metric (e.g. Correlation), the images aren't that big (~400k voxels) and a previous rigid registration phase behaved correctly.
Here is my code snippet, what am I doing wrong ? (Image overlap had been ensured)

   typedef itk::Image<float,3> ImageTypef;

   itk::ImageRegistrationMethodv4<ImageTypef,ImageTypef>::Pointer WarpHelper=
      itk::ImageRegistrationMethodv4<ImageTypef,ImageTypef>::New();
   itk::MattesMutualInformationImageToImageMetricv4<ImageTypef,ImageTypef>::Pointer MIMetrics=
      itk::MattesMutualInformationImageToImageMetricv4<ImageTypef,ImageTypef>::New();
   itk::BSplineTransform<double,3,3>::Pointer BSplineTransform=
      itk::BSplineTransform<double,3,3>::New();
   itk::LBFGSOptimizerv4::Pointer Optimizer=
      itk::LBFGSOptimizerv4::New();
   
   MIMetrics->SetNumberOfHistogramBins(50);
   MIMetrics->SetUseMovingImageGradientFilter(false);
   MIMetrics->SetUseFixedImageGradientFilter(false);
   MIMetrics->SetUseFixedSampledPointSet(false);
   MIMetrics->SetVirtualDomainFromImage(TestIm);
   
   
   WarpHelper->SetFixedImage(TestIm);
   WarpHelper->SetMetricSamplingPercentage(0.002);
   WarpHelper->SetMovingImage(MovIm);
   WarpHelper->SetMetric(MIMetrics);
   //WarpHelper->SetMetricSamplingStrategy(WarpHelper->RANDOM);
   WarpHelper->InPlaceOn();
   WarpHelper->SetInitialTransform(BSplineTransform);
   WarpHelper->SetOptimizer(Optimizer);

   itk::Array<double> ShrinkFactorPerLevel(1);
   ShrinkFactorPerLevel[0]=1;
   itk::Array<double> SmoothSigmaPerLevel(1);
   SmoothSigmaPerLevel[0]=0;
   WarpHelper->SetNumberOfLevels(1);
   WarpHelper->SetSmoothingSigmasPerLevel(SmoothSigmaPerLevel);
   WarpHelper->SetShrinkFactorsPerLevel(ShrinkFactorPerLevel);

   /* ----- Initialize the Bspline Transform ----- */
   BSplineTransform->SetTransformDomainOrigin(TestIm->GetOrigin());
   BSplineTransform->SetTransformDomainDirection(TestIm->GetDirection());

   itk::Vector<double,3> PhysicalDim;
   for (int dim = 0; dim < 3; dim++)
   {
      PhysicalDim=(TestIm->GetSpacing()[dim])*float((TestIm->GetLargestPossibleRegion().GetSize()[dim])-1);
   }
   BSplineTransform->SetTransformDomainPhysicalDimensions(PhysicalDim);

   itk::Size<3> meshSize;
   meshSize.Fill(10-3);// NbGridNodesEachDim-SplineOrder
   BSplineTransform->SetTransformDomainMeshSize(meshSize);

   BSplineTransform->SetIdentity();

   /* ----- Optimizer setup ----- */
   itk::Array<double> ScaleBspline(BSplineTransform->GetNumberOfParameters());
   ScaleBspline.fill(1);

   Optimizer->SetLineSearchAccuracy(0.9);
   Optimizer->SetDefaultStepLength(1.5);
   Optimizer->TraceOn();
   Optimizer->SetGradientConvergenceTolerance(1e-3);
   Optimizer->SetMaximumNumberOfFunctionEvaluations(20);
   Optimizer->SetDoEstimateScales(false);
   Optimizer->SetScales(ScaleBspline);


   try
   {
      WarpHelper->Update();
   }
   catch (itk::ExceptionObject &err)
   {
      std::cerr<<err;
      return 0;
   }

Thanks,

Tim