[Insight-users] NonRigidRegistration using MutualInformation

[=?ISO-2022-JP?B?GyRCPj5IeDcwGyhC?=]

Dear Insight user,

I am trying to 3D NonRigidRegistration(CT-CT).
The components i'm using are
BSplineDeformableTransform,
LBFGSBOptimizer,
NormalizedMutualInformationHistogramImageToImageMetric
(don't use MultiResolution)

I have a 3D images 280*280*201 and spacing 1*1*1.
The time it takes to calculate each itertion is veeeeery slow.(Processing
doesn't end......)

How can I solve that????

Regards.
            Matsuo

Code----------------------------------------------------------------------------------
    const unsigned int ImageDimension = 3;
    typedef signed short PixelType;

    typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
    typedef itk::Image< PixelType, ImageDimension >  MovingImageType;

    const unsigned int SpaceDimension = ImageDimension;
    const unsigned int SplineOrder = 3;
    typedef double CoordinateRepType;
    typedef itk::BSplineDeformableTransform< CoordinateRepType,
SpaceDimension, SplineOrder > TransformType;
    typedef itk::LBFGSBOptimizer OptimizerType;
    typedef itk::NormalizedMutualInformationHistogramImageToImageMetric<
FixedImageType, MovingImageType > MetricType;
    typedef itk::ImageRegistrationMethod< FixedImageType, MovingImageType >
RegistrationType;
    typedef itk::LinearInterpolateImageFunction< MovingImageType, double >
InterpolatorType;

    MetricType::Pointer metric = MetricType::New();
    OptimizerType::Pointer optimizer = OptimizerType::New();
    InterpolatorType::Pointer interpolator = InterpolatorType::New();
    RegistrationType::Pointer registration = RegistrationType::New();
    TransformType::Pointer transform = TransformType::New();

    registration->SetMetric( metric );
    registration->SetOptimizer( optimizer );
    registration->SetInterpolator( interpolator );
    registration->SetTransform( transform );

    typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
    typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;

    FixedImageReaderType::Pointer fixedImageReader =
FixedImageReaderType::New();
    MovingImageReaderType::Pointer movingImageReader =
MovingImageReaderType::New();

    fixedImageReader->SetFileName( argv[1] );
    movingImageReader->SetFileName( argv[2] );
    FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();

    registration->SetFixedImage( fixedImage );
    registration->SetMovingImage( movingImageReader->GetOutput() );

    fixedImageReader->Update();

    FixedImageType::RegionType fixedRegion =
fixedImage->GetBufferedRegion();
    registration->SetFixedImageRegion( fixedRegion );

    typedef TransformType::RegionType RegionType;
    RegionType bsplineRegion;
    RegionType::SizeType gridSizeOnImage;
    RegionType::SizeType gridBorderSize;
    RegionType::SizeType totalGridSize;

    gridSizeOnImage.Fill( 8 );
    gridBorderSize.Fill( 3 );
    totalGridSize = gridSizeOnImage + gridBorderSize;

    bsplineRegion.SetSize( totalGridSize );

    typedef TransformType::SpacingType SpacingType;
    SpacingType spacing = fixedImage->GetSpacing();

    typedef TransformType::OriginType OriginType;
    OriginType origin = fixedImage->GetOrigin();;

    FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();

    for(unsigned int r=0; r<ImageDimension; r++)
    {
        spacing[r] *= floor( static_cast<double>(fixedImageSize[r] - 1)  /
                    static_cast<double>(gridSizeOnImage[r] - 1) );
        origin[r]  -=  spacing[r];
    }

    transform->SetGridSpacing( spacing );
    transform->SetGridOrigin( origin );
    transform->SetGridRegion( bsplineRegion );


    typedef TransformType::ParametersType     ParametersType;
    const unsigned int numberOfParameters =
transform->GetNumberOfParameters();
    ParametersType parameters( numberOfParameters );
    parameters.Fill( 0.0 );
    transform->SetParameters( parameters );
    registration->SetInitialTransformParameters( transform->GetParameters()
);

    unsigned int numberOfHistogramBins = 256;
    MetricType::HistogramType::SizeType histogramSize;
    histogramSize[0] = numberOfHistogramBins;
    histogramSize[1] = numberOfHistogramBins;
    metric->SetHistogramSize( histogramSize );

    typedef MetricType::ScalesType ScalesType;
    ScalesType scales( numberOfParameters );

    scales.Fill( 1.0 );

    metric->SetDerivativeStepLengthScales(scales);

    OptimizerType::BoundSelectionType boundSelect(
transform->GetNumberOfParameters() );
    OptimizerType::BoundValueType upperBound(
transform->GetNumberOfParameters() );
    OptimizerType::BoundValueType lowerBound(
transform->GetNumberOfParameters() );

    boundSelect.Fill( 0 );
    upperBound.Fill( 0.0 );
    lowerBound.Fill( 0.0 );

    optimizer->SetBoundSelection( boundSelect );
    optimizer->SetUpperBound( upperBound );
    optimizer->SetLowerBound( lowerBound );

    optimizer->SetCostFunctionConvergenceFactor( 1e+7);
    optimizer->SetProjectedGradientTolerance(1e-4 );
    optimizer->SetMaximumNumberOfIterations( 500 );
    optimizer->SetMaximumNumberOfEvaluations( 500 );
    optimizer->SetMaximumNumberOfCorrections( 12 );

    CommandIterationUpdate::Pointer observer =
CommandIterationUpdate::New();
    optimizer->AddObserver( itk::IterationEvent(), observer );

    try
    {
        registration->StartRegistration();
    }
    catch( itk::ExceptionObject & err )
    {
        std::cerr << "ExceptionObject caught !" << std::endl;
        std::cerr << err << std::endl;
        return -1;
    }
    ::
    ::
--------------------------------------------------------------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20071101/8d91894e/attachment.html