[Insight-users] Demons 3D deformable registration problem

[姜华]

Hi,luis,
   I have several questions about ITK deformable registration.
  First, I use the Demons for 3D deformable registration, when I set the
parameter "SetStandardDeviations" smaller(for example 1.0),
  the image whose size is 256*176*16 will have black regions in each slice,
but when I set this paremeter bigger (for example 10.0),
  only the first and the last slice have  black regions, so I do not know
why this phenomena will happen?
  second, since the phenomena happened, so I use the "MirrorPadImageFilter"
to pad the image in the slice direction in order to
  make its size 256*176*18, and set the "SetStandardDeviations"  parameter
10.0, then after registration,
  I use the WarpImageFilter to warp the movingImage according to the
deformable field, at last I use the "RegionOfInterestImageFilter" extract
the image from second slice to the seventeenth slice. Thus I can get better
result, but I don't know this method is whether correct.
    In order to explain these problem, I copy the part programs in the
email, and the result images in the inclosure.
 Thank you for your answers as soon as possible.

 With best regards!

                 HuaJiang,


The Program:(this is the part of the my program, the read and write parts
arenot attached)
FixedImageCasterType::Pointer fixedImageCaster   =
FixedImageCasterType::New();
    MovingImageCasterType::Pointer movingImageCaster =
MovingImageCasterType::New();
 fixedImageCaster->SetInput( fixedImage );
    movingImageCaster->SetInput( movingImage );
/*
 HistogramMatchFilterType::Pointer matcher =
HistogramMatchFilterType::New();
 matcher->SetInput( movingImageCaster->GetOutput() );
    matcher->SetReferenceImage( fixedImageCaster->GetOutput() );

 matcher->SetNumberOfHistogramLevels( 1024 );
    matcher->SetNumberOfMatchPoints( 7 );
    matcher->ThresholdAtMeanIntensityOn();
    */


 Demons3RegistrationFilterType::Pointer filter =
Demons3RegistrationFilterType::New();

 filter->SetStandardDeviations( 10.0 );
 MultiResRegistrationFilterType:: Pointer multires =
MultiResRegistrationFilterType::New();

 multires->SetRegistrationFilter( filter );
 multires->SetNumberOfLevels( 3 );

 unsigned int nIterations[4] = {32, 16, 4 };
    multires->SetNumberOfIterations( nIterations );
 DemonsCommandIterationUpdate::Pointer observer =
DemonsCommandIterationUpdate::New();
 filter->AddObserver(itk::IterationEvent(), observer);

 CommandResolutionLevelUpdate::Pointer levelobserver =
CommandResolutionLevelUpdate::New();
    multires->AddObserver( itk::IterationEvent(), levelobserver );

 // Pad the image.
    unsigned long upFactor[3] = {0,0,1};
 unsigned long lowFactor[3] = {0,0,1};
    InternalPixelType defaultPixel = 0;
 // define warp filter.
 WarperType::Pointer warper = WarperType::New();
 typedef itk::LinearInterpolateImageFunction< InternalImageType, double >
LinearInterpolatorType;
 LinearInterpolatorType::Pointer linearInterpolator =
LinearInterpolatorType::New();

 warper->SetInterpolator( linearInterpolator );
    ConstantPadImageFilterType::Pointer constantPadFixedFilter;
 ConstantPadImageFilterType::Pointer constantPadMovingFilter;
 MirrorPadImageFilterType::Pointer mirrorPadFixedFilter;
 MirrorPadImageFilterType::Pointer mirrorPadMovingFilter;
 if (m_padMode == ConstantPad)
 {
     constantPadFixedFilter = ConstantPadImageFilterType::New();
     constantPadMovingFilter = ConstantPadImageFilterType::New();

  constantPadFixedFilter->SetInput(fixedImageCaster->GetOutput());
  constantPadMovingFilter->SetInput(movingImageCaster->GetOutput());

  constantPadFixedFilter->SetPadLowerBound(lowFactor);
  constantPadFixedFilter->SetPadUpperBound(upFactor);
        constantPadFixedFilter->SetConstant(defaultPixel);
        constantPadMovingFilter->SetPadLowerBound(lowFactor);
  constantPadMovingFilter->SetPadUpperBound(upFactor);
        constantPadMovingFilter->SetConstant(defaultPixel);
  multires->SetFixedImage( constantPadFixedFilter->GetOutput() );
     multires->SetMovingImage( constantPadMovingFilter->GetOutput() );

  try
  {
   multires->Update();
  }
  catch( itk::ExceptionObject & excp )
  {
   std::cerr << excp << std::endl;
   return ;
    }
        warper->SetInput(constantPadMovingFilter->GetOutput());
        warper->SetOutputSpacing(
constantPadFixedFilter->GetOutput()->GetSpacing() );
  warper->SetOutputOrigin( constantPadFixedFilter->GetOutput()->GetOrigin()
);

  warper->SetDeformationField( multires->GetOutput() );
        warper->Update();
 }
 else if (m_padMode == MirrorPad)
 {
     mirrorPadFixedFilter = MirrorPadImageFilterType::New();
  mirrorPadMovingFilter = MirrorPadImageFilterType::New();

  mirrorPadFixedFilter->SetInput(fixedImageCaster->GetOutput());
  mirrorPadMovingFilter->SetInput(movingImageCaster->GetOutput());

  mirrorPadFixedFilter->SetPadLowerBound(lowFactor);
  mirrorPadFixedFilter->SetPadUpperBound(upFactor);

        mirrorPadMovingFilter->SetPadLowerBound(lowFactor);
  mirrorPadMovingFilter->SetPadUpperBound(upFactor);

  multires->SetFixedImage( mirrorPadFixedFilter->GetOutput() );
  multires->SetMovingImage( mirrorPadMovingFilter->GetOutput() );

  try
  {
   multires->Update();
  }
  catch( itk::ExceptionObject & excp )
  {
   std::cerr << excp << std::endl;
   return ;
  }

  warper->SetInput(mirrorPadMovingFilter->GetOutput());
        warper->SetOutputSpacing(
mirrorPadFixedFilter->GetOutput()->GetSpacing() );
  warper->SetOutputOrigin( mirrorPadFixedFilter->GetOutput()->GetOrigin() );

  warper->SetDeformationField( multires->GetOutput() );
        warper->Update();
 }

    // extract the image.
    RegionOfInterestImageFilterType::Pointer regionExtractFitler
  = RegionOfInterestImageFilterType::New();

    regionExtractFitler->SetInput(warper->GetOutput());
 InternalImageType::RegionType desiredRegion ;
 InternalImageType::SizeType size =
fixedImage->GetLargestPossibleRegion().GetSize();
 InternalImageType::IndexType index;
 index[0] = 0;
 index[1] = 0;
 index[2] = 0;
    desiredRegion.SetSize(size);
 desiredRegion.SetIndex(index);
 regionExtractFitler->SetRegionOfInterest(desiredRegion);
 InternalImageCasterType::Pointer internalCaster =
InternalImageCasterType::New();
 internalCaster->SetInput(regionExtractFitler->GetOutput());
    internalCaster->Update();
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091230/bcc66d39/attachment.htm>