[Insight-users] White Output after Rigid Registration

[Andy Eow]

Hi,

I'm trying to perform a straightforward rigid registration in 3D using sample code from the MultiResMIRegistration application. The transform parameters after registration are as follows for my MRI (T1 & T2) dataset: [0.61162, 0.0627385, 0.0105601, 0.78859, 1.63678, -0.662944, 0.223458]. However, I suspect my post registration resampling is not working because the resulting volume from that is completely white. Attached is the exact code that I'm running. Any advice is appreciated. Thanks!!!

VolumeType::Pointer MakeSymmetrical( VolumeType::Pointer fixedVol, VolumeType::Pointer movingVol ) {
  typedef itk::ChangeInformationImageFilter< VolumeType > CentererType;
  typedef itk::NormalizeImageFilter< VolumeType, VolumeType > NormalizerType;
  typedef itk::QuaternionRigidTransform< double > TransformType;
  typedef itk::QuaternionRigidTransformGradientDescentOptimizer OptimizerType;
  typedef itk::MutualInformationImageToImageMetric< VolumeType, VolumeType > MetricType;
  typedef itk::LinearInterpolateImageFunction< VolumeType, double > InterpolatorType;
  typedef itk::RecursiveMultiResolutionPyramidImageFilter< VolumeType, VolumeType > FixedImagePyramidType;
  typedef itk::RecursiveMultiResolutionPyramidImageFilter< VolumeType, VolumeType > MovingImagePyramidType;
  typedef itk::MultiResolutionImageRegistrationMethod< VolumeType, VolumeType > RegistrationType;
  typedef itk::FixedArray< unsigned int, 3 > ShrinkFactorsArray;
  typedef itk::Array< unsigned int > UnsignedIntArray;
  typedef itk::Array< double > DoubleArray;
  typedef itk::AffineTransform< double, 3 > AffineTransformType;
  typedef AffineTransformType::ScalarType CoordRepType;
  typedef itk::LinearInterpolateImageFunction< VolumeType, CoordRepType > ResampleInterpolatorType;
  typedef itk::ResampleImageFilter< VolumeType, VolumeType > ResamplerType;

  // Parameters
  unsigned short noOfLevels = 1;
  double translationScale = 1.0;
  double movingVolStdDev = 0.4;
  double fixedVolStdDev = 0.4;
  unsigned short noOfSpatialSamples = 50;
  DoubleArray learningRates = DoubleArray( noOfLevels );
  learningRates.Fill( 1e-4 );
  UnsignedIntArray noOfIterations = UnsignedIntArray( noOfLevels );
  noOfIterations.Fill( 2500 );
  ShrinkFactorsArray fixedVolShrinkFactors;
  fixedVolShrinkFactors.Fill( 1 );
  ShrinkFactorsArray movingVolShrinkFactors;
  movingVolShrinkFactors.Fill( 1 );

  // Fixed Volume
  CentererType::Pointer centerer =  CentererType::New();
  centerer->CenterImageOn();
  centerer->SetInput( fixedVol );

  NormalizerType::Pointer normalizer = NormalizerType::New();  
  normalizer->SetInput( centerer->GetOutput() );
  ITK_TRY_MACRO( normalizer->Update() );
  VolumeType::Pointer tFixedVol = normalizer->GetOutput();

  // Moving Volume
  CentererType::Pointer centerer2 =  CentererType::New();
  centerer2->CenterImageOn();
  centerer2->SetInput( movingVol );

  NormalizerType::Pointer normalizer2 = NormalizerType::New();  
  normalizer2->SetInput( centerer2->GetOutput() );
  ITK_TRY_MACRO( normalizer2->Update() );
  VolumeType::Pointer tMovingVol = normalizer2->GetOutput();

  // Setup the registrator
  TransformType::Pointer transform = TransformType::New();
  RegistrationType::ParametersType transformInitialParameters = RegistrationType::ParametersType( transform->GetNumberOfParameters() );
  transformInitialParameters.Fill( 0.0 );
  transformInitialParameters[3] = 1.0;

  OptimizerType::Pointer optimizer = OptimizerType::New();
  OptimizerType::ScalesType scales( transform->GetNumberOfParameters() );
  scales.Fill( 1.0 );
  for ( int i = 4; i < 7; i++ )
    scales[i] = translationScale;
  optimizer->SetScales( scales );
  optimizer->MaximizeOn();

  MetricType::Pointer metric = MetricType::New();
  metric->SetMovingImageStandardDeviation( movingVolStdDev );
  metric->SetFixedImageStandardDeviation( fixedVolStdDev );
  metric->SetNumberOfSpatialSamples( noOfSpatialSamples );

  FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
  fixedImagePyramid->SetNumberOfLevels( noOfLevels );
  fixedImagePyramid->SetStartingShrinkFactors( fixedVolShrinkFactors.GetDataPointer() );

  MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();  
  movingImagePyramid->SetNumberOfLevels( noOfLevels );
  movingImagePyramid->SetStartingShrinkFactors( movingVolShrinkFactors.GetDataPointer() );

  InterpolatorType::Pointer interpolator = InterpolatorType::New();

  RegistrationType::Pointer registration  = RegistrationType::New();
  registration->SetTransform( transform );
  registration->SetOptimizer( optimizer );
  registration->SetMetric( metric );
  registration->SetInterpolator( interpolator );
  registration->SetFixedImagePyramid( fixedImagePyramid );
  registration->SetMovingImagePyramid( movingImagePyramid );
  registration->SetFixedImage( tFixedVol );
  registration->SetMovingImage( tMovingVol );
  registration->SetNumberOfLevels( noOfLevels );
  registration->SetInitialTransformParameters( transformInitialParameters );
  registration->SetFixedImageRegion( tFixedVol->GetBufferedRegion() );
  ITK_TRY_MACRO( registration->StartRegistration() );
  MESSAGE( "Transform Parameters (post-reg): " << transform->GetParameters() );
  MESSAGE( "Rotation (post-reg): " << transform->GetRotation() );
  MESSAGE( "Translation (post-reg): " << transform->GetTranslation() );

  // Setup the resampler
  AffineTransformType::Pointer resamTransform = AffineTransformType::New();
  resamTransform->SetIdentity();
  resamTransform->SetParameters( transform->GetParameters() );
  resamTransform->SetMatrix( transform->GetRotationMatrix() );
  resamTransform->SetOffset( transform->GetOffset() );
  
  ResampleInterpolatorType::Pointer resamInterpolator = ResampleInterpolatorType::New();

  ResamplerType::Pointer resampler = ResamplerType::New();
  resampler->SetInput( tMovingVol );
  resampler->SetTransform( resamTransform );
  resampler->SetInterpolator( resamInterpolator );
  resampler->SetSize( tFixedVol->GetLargestPossibleRegion().GetSize() );
  resampler->SetOutputOrigin( tFixedVol->GetOrigin() );
  resampler->SetOutputSpacing( tFixedVol->GetSpacing() );
  resampler->SetDefaultPixelValue( 0 );
  ITK_TRY_MACRO( resampler->Update() );
  
  return resampler->GetOutput();
}


Cheers,
Andy
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20040927/8349cf7e/attachment-0001.html