Как рассчитать взаимную информацию Mattes, используя мало памяти? - PullRequest
0 голосов
/ 14 октября 2019

Я получаю сообщение об ошибке при попытке увеличить количество контрольных точек в регистрации изображений с помощью B-сплайнов. Ошибка, которую я получил:

завершить вызов после выброса экземпляра 'itk :: ExceptionObject' what (): / home / felippe / SimpleITK-build / ITK / Modules / Core / Common / src/itkMultiThreader.cxx:399: itk :: ERROR: MultiThreader (0x4f58f40): во время SingleMethodExecute произошло исключение std :: bad_alloc

Код, который я использую:

// This one header will include all SimpleITK filters and external
// objects.
#include <SimpleITK.h>
#include <iostream>
#include <stdlib.h>
#include <iomanip>

namespace sitk = itk::simple;

// use sitk's output operator for std::vector etc..
using sitk::operator<<;


class IterationUpdate
  : public sitk::Command
{
public:
  IterationUpdate( const sitk::ImageRegistrationMethod &m, const sitk::BSplineTransform &tx)
    : m_Method(m),
      m_BSplineTransform(tx)
    {}

  // Override method from parent which is called at for the requested event
  virtual void Execute( )
    {
      if (m_Method.GetOptimizerIteration() == 0)
        {
        // The BSpline is resized before the first optimizer
        // iteration is completed per level. Print the transform object
        // to show the adapted BSpline transform.
        std::cout << m_BSplineTransform.ToString() << std::endl;
        }

      // stash the stream state
      std::ios  state(NULL);
      state.copyfmt(std::cout);
      std::cout << std::fixed << std::setfill(' ') << std::setprecision( 5 );
      std::cout << std::setw(3) << m_Method.GetOptimizerIteration();
      std::cout << " = " << std::setw(10) << m_Method.GetMetricValue() << std::endl;
      std::cout.copyfmt(state);
    }

private:
  const sitk::ImageRegistrationMethod &m_Method;
  const sitk::BSplineTransform &m_BSplineTransform;

};


class MultiResolutionIterationUpdate
  : public sitk::Command
{
public:
  MultiResolutionIterationUpdate( const sitk::ImageRegistrationMethod &m)
    : m_Method(m)
    {}

  // Override method from parent which is called at for the requested event
  virtual void Execute( )
    {
      // The sitkMultiResolutionIterationEvent occurs before the
      // resolution of the transform. This event is used here to print
      // the status of the optimizer from the previous registration level.
      if (m_Method.GetCurrentLevel() > 0)
        {
        std::cout << "Optimizer stop condition: " << m_Method.GetOptimizerStopConditionDescription() << std::endl;
        std::cout << " Iteration: " << m_Method.GetOptimizerIteration() << std::endl;
        std::cout << " Metric value: " << m_Method.GetMetricValue() << std::endl;
        }

      std::cout << "--------- Resolution Changing ---------" << std::endl;
    }

private:
  const sitk::ImageRegistrationMethod &m_Method;
};




int main()
{

    //******************************************************
    //Reading the Fixed Image
    sitk::Image fixed = sitk::ReadImage( "/home/felippe/Área de Trabalho/Felippe/Mestrado/REGISTRO/Deformation_Registering_MATLAB/Volumes/Teste1/MNI152_T1_0.5mm.nii", sitk::sitkFloat32 );
    //fixed = sitk::Normalize( fixed );

    //Printing the dimensions of the fixed image
    std::vector<unsigned int> fixed_dims = fixed.GetSize();
    std::cout << "Fixed Image Dimensions: ";
    for (auto i = fixed_dims.begin(); i != fixed_dims.end(); ++i)
        std::cout << *i << ' ';
    std::cout << std::endl << "Fixed Image Pixel Type: " << fixed.GetPixelIDTypeAsString() << std::endl;
    std::cout << std::endl;

    //Reading the Moving Image
    sitk::Image moving = sitk::ReadImage( "/home/felippe/Área de Trabalho/Felippe/Mestrado/C_plus_plus/Codigos/build-Registration_ITK_CMAKE-Desktop_Qt_5_12_3_GCC_64bit-Default/mri_vbm6_transformed_affine.mha", sitk::sitkFloat32 );
    //moving = sitk::Normalize( moving );

    //Printing the dimensions of the moving image
    std::vector<unsigned int> moving_dims = moving.GetSize();
    std::cout << "Moving Image Dimensions: ";
    for (auto i = moving_dims.begin(); i != moving_dims.end(); ++i)
        std::cout << *i << ' ';
    std::cout << std::endl << "Moving Image Pixel Type: " << moving.GetPixelIDTypeAsString() << std::endl;
    std::cout << std::endl;
    std::cout << std::endl;
    //****************************************************************************

    std::cout << "Spacing in the Fixed Image: ";
    std::vector<double> fixed_spacing = fixed.GetSpacing();
    for (auto i = fixed_spacing.begin(); i != fixed_spacing.end(); ++i)
        std::cout << *i << ' ';
    std::cout << std::endl << std::endl;

    std::vector<unsigned int> transformDomainMeshSize(fixed.GetDimension(), 4);
    sitk::BSplineTransform tx = sitk::BSplineTransformInitializer(fixed, transformDomainMeshSize);

    std::cout << "Transform Domain Mesh Size: ";
    for (auto i = transformDomainMeshSize.begin(); i != transformDomainMeshSize.end(); ++i)
        std::cout << *i << ' ';
    std::cout << std::endl << std::endl;


    std::cout << "Initial Number of Parameters:" << tx.GetNumberOfParameters() << std::endl;

    sitk::ImageRegistrationMethod R;
    unsigned int number_of_histogram_bins = 32;
    std::vector<double> samplingPercentage(3);
    samplingPercentage[0] = 0.01;
    samplingPercentage[1] = 0.01;
    samplingPercentage[2] = 0.01;
    R.SetMetricAsMattesMutualInformation(number_of_histogram_bins);
    R.SetMetricSamplingStrategy(R.RANDOM);
    R.SetMetricSamplingPercentagePerLevel(samplingPercentage);
    R.MetricUseFixedImageGradientFilterOn();

    const double learningRate = 5.0;
    const unsigned int numberOfIterations = 100u;
    const double convergenceMinimumValue = 1e-4;
    const unsigned int convergenceWindowSize = 5;

    R.SetOptimizerAsGradientDescentLineSearch( learningRate,
                                                numberOfIterations,
                                                convergenceMinimumValue,
                                                convergenceWindowSize);

    //R.SetOptimizerAsLBFGSB(convergenceMinimumValue, numberOfIterations);
    R.SetInterpolator(sitk::sitkLinear);

    const unsigned int numberOfLevels = 3;
    std::vector<unsigned int> scaleFactors(numberOfLevels);
    scaleFactors[0] = 1;
    scaleFactors[1] = 2;
    scaleFactors[2] = 5;
    const bool inPlace = true;
    R.SetInitialTransformAsBSpline(tx,
                                   inPlace,
                                   scaleFactors);

    std::vector<unsigned int> shrinkFactors( numberOfLevels );
    shrinkFactors[0] = 4;
    shrinkFactors[1] = 2;
    shrinkFactors[2] = 1;
      R.SetShrinkFactorsPerLevel( shrinkFactors );

      std::vector<double> smoothingSigmas( numberOfLevels );
      smoothingSigmas[0] = 4.0;
      smoothingSigmas[1] = 2.0;
      smoothingSigmas[2] = 1.0;
      R.SetSmoothingSigmasPerLevel( smoothingSigmas );

      IterationUpdate cmd1(R, tx);
      R.AddCommand( sitk::sitkIterationEvent, cmd1);

      MultiResolutionIterationUpdate cmd2(R);
      R.AddCommand( sitk::sitkMultiResolutionIterationEvent, cmd2);

      std::cout << "Initializing the Registration!!" << std::endl << std::endl;
      sitk::Transform outTx = R.Execute( fixed, moving );

      std::cout << "-------" << std::endl;
      std::cout << outTx.ToString() << std::endl;
      std::cout << "Optimizer stop condition: " << R.GetOptimizerStopConditionDescription() << std::endl;
      std::cout << " Iteration: " << R.GetOptimizerIteration() << std::endl;
      std::cout << " Metric value: " << R.GetMetricValue() << std::endl;

      //Saving the transformation
      sitk::WriteTransform(outTx, "simple_transform.tfm");

      //Applying the transformation
      sitk::Image transformed_image = sitk::Resample(moving, fixed, outTx, sitk::sitkLinear, 0, moving.GetPixelID());
      sitk::WriteImage(transformed_image, "mri_vbm6_transformed_Affine_BS_444.mha");

      return 0;transformDomainMeshSize
}

Когда я устанавливаю transformDomainMeshSize в 2, я могу получить успех при регистрации, но когда я устанавливаю 4, я получаю эту ошибку, которая, я думаю, связана с утечкой памяти. Так, как я могу решить эту проблему?

ОБНОВЛЕНИЕ: я исследую, что проблема с mattesMutualInformation. Есть ли способ рассчитать эту метрику, используя меньше памяти?

Спасибо.

...