[Insight-users] Registration problem with MI metric and BSpline
transform.
Zien
superz9 at gmail.com
Sun Mar 12 12:05:24 EST 2006
Hi,all :
I want to register two image with viola MI metric,BSpline transform
and LBFGS Optimizer.
My code as below.
I don't know what happen,but the result is so strange.
The output image that is almost the same with moving image.
Have anybody give me some point ?
Best Regards.
The result :
Intial Parameters =
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Starting Registration
F = 0.302491, GNORM = 0.00259936
*************************************************
I NFN FUNC GNORM STEPLENGTH
1 3 0.224 0.003 1923.552
THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.
Last Transform Parameters
[3.48814e-006, -2.2697e-005, -0.000105602, 0.000187457, 4.73857e-005, -2.72309e-
005, -8.18973e-007, 0, 0.00152272, 0.0527738, 0.0645476, 0.0819562, 0.0512163, 0
.00437728, -2.74508e-005, 0, 0.0204939, 1.00828, 1.45459, 0.574319, 1.04031, 0.2
42405, 0.000572549, 0, 0.0660639, 1.48627, 1.77614, 0.724309, 2.13744, 0.443386,
-1.12093e-005, 0, 0.0912913, 1.45534, 1.33712, 0.00705034, -0.286602, -0.135778
, -0.000405499, 0, 0.023365, 0.390063, 0.379199, -0.14935, -0.133261, -0.0124785
, -6.45218e-005, 0, 0.000257838, 0.00251562, 0.00757467, -0.0225716, -0.0131375,
-0.000192109, -7.24246e-006, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4.34806e-006, -3.34649e
-005, -0.000235892, 2.62016e-005, -7.19323e-005, 3.263e-005, 5.19492e-007, 0, 0.
000590144, -0.000824211, 0.0544006, 0.379423, 0.0695502, -0.00139731, 1.35357e-0
05, 0, 0.00622242, 0.315731, 0.634949, 1.08377, -0.0456413, -0.0789634, -2.99736
e-005, 0, -0.0102174, 0.101388, 0.266624, -0.0263321, -0.576348, -0.156506, 2.41
444e-005, 0, -0.0356267, -0.624992, -0.761349, -0.545143, -0.208911, -0.00781882
, 0.000177452, 0, -0.0119721, -0.227898, -0.618297, -1.10105, -0.315192, -0.0007
10803, -7.29731e-006, 0, -8.15398e-005, -0.00427512, -0.0742942, -0.164705, -0.0
475115, -0.000819229, -6.3068e-006, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Probe Tag Starts Stops Time
Registration 1 1 2.30873
My code :
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkImageRegistrationMethod.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkImage.h"
#include "itkMutualInformationImageToImageMetric.h"//
#include "itkNormalizeImageFilter.h"//
#include "itkDiscreteGaussianImageFilter.h"//
#include "itkTimeProbesCollectorBase.h"
#include "itkBSplineDeformableTransform.h"
#include "itkLBFGSOptimizer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSquaredDifferenceImageFilter.h"
#include "itkBMPImageIO.h"
int main( int argc, char *argv[] )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile outputImagefile ";
std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";
std::cerr << " [deformationField] ";
return 1;
}
const unsigned int ImageDimension = 2;
typedef float PixelType;
typedef itk::Image< PixelType, ImageDimension > FixedImageType;
typedef itk::Image< PixelType, ImageDimension > MovingImageType;
typedef itk::Image< InternalPixelType, ImageDimension > InternalImageType;
typedef double CoordinateRepType;
const unsigned int SpaceDimension = ImageDimension;//Dimension=2
const unsigned int SplineOrder = 3;
typedef itk::BSplineDeformableTransform<
CoordinateRepType,
SpaceDimension,
SplineOrder > TransformType;
typedef itk::LBFGSOptimizer OptimizerType;
typedef itk:: LinearInterpolateImageFunction<
InternalImageType,
double > InterpolatorType;
//Inheritance diagram for itk::ImageRegistrationMethod< TFixedImage,
TMovingImage >:
typedef itk::ImageRegistrationMethod<
InternalImageType,
InternalImageType > RegistrationType;
//============================================================================
//Normalize the simple two image
typedef itk::MutualInformationImageToImageMetric<
InternalImageType,
InternalImageType > MetricType;
typedef itk::NormalizeImageFilter<
FixedImageType,
InternalImageType
> FixedNormalizeFilterType;
typedef itk::NormalizeImageFilter<
MovingImageType,
InternalImageType
> MovingNormalizeFilterType;
FixedNormalizeFilterType::Pointer fixedNormalizer =
FixedNormalizeFilterType::New();
MovingNormalizeFilterType::Pointer movingNormalizer =
MovingNormalizeFilterType::New();
//Gaussian Filter
typedef itk::DiscreteGaussianImageFilter<
InternalImageType,
InternalImageType
> GaussianFilterType;
GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
fixedSmoother->SetVariance( 2.0 );
movingSmoother->SetVariance( 2.0 );
//============================================================================
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
TransformType::Pointer transform = TransformType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
registration->SetTransform( transform );
metric->SetFixedImageStandardDeviation( 0.4 );
metric->SetMovingImageStandardDeviation( 0.4 );
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] );
fixedImageReader->Update();
movingImageReader->Update();
//FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
//===============================================================================
fixedNormalizer->SetInput( fixedImageReader->GetOutput() );
movingNormalizer->SetInput( movingImageReader->GetOutput() );
fixedSmoother->SetInput( fixedNormalizer->GetOutput() );
movingSmoother->SetInput( movingNormalizer->GetOutput() );
registration->SetFixedImage( fixedSmoother->GetOutput() );
registration->SetMovingImage( movingSmoother->GetOutput() );
fixedNormalizer->Update();
movingNormalizer->Update();
fixedSmoother->Update();
movingSmoother->Update();
FixedImageType::RegionType fixedImageRegion =
fixedNormalizer->GetOutput()->GetBufferedRegion();
registration->SetFixedImageRegion(fixedImageRegion);
//==========================================================================
//==========================================================================
const unsigned int numberOfPixels = fixedImageRegion.GetNumberOfPixels();
const unsigned int numberOfSamples =
static_cast< unsigned int >( numberOfPixels * 0.01 );
metric->SetNumberOfSpatialSamples( numberOfSamples );
//==========================================================================
typedef TransformType::RegionType RegionType;
RegionType bsplineRegion;
RegionType::SizeType gridSizeOnImage;
RegionType::SizeType gridBorderSize;
RegionType::SizeType totalGridSize;
gridSizeOnImage.Fill( 5 );
gridBorderSize.Fill( 3 );
totalGridSize = gridSizeOnImage + gridBorderSize;
bsplineRegion.SetSize( totalGridSize );
typedef TransformType::SpacingType SpacingType;
SpacingType spacing = fixedNormalizer->GetOutput()->GetSpacing();//
typedef TransformType::OriginType OriginType;
OriginType origin = fixedNormalizer->GetOutput()->GetOrigin();;//
FixedImageType::SizeType fixedImageSize = fixedImageRegion.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() );
std::cout << "Intial Parameters = " << std::endl;
std::cout << transform->GetParameters() << std::endl;
optimizer->SetGradientConvergenceTolerance( 0.05 );
optimizer->SetLineSearchAccuracy( 0.9 );
optimizer->SetDefaultStepLength( 1.5 );
optimizer->TraceOn();
optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
itk::TimeProbesCollectorBase collector;
std::cout << std::endl << "Starting Registration" << std::endl;
try
{
collector.Start( "Registration" );
registration->StartRegistration();
collector.Stop( "Registration" );
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
std::cout << "Last Transform Parameters" << std::endl;
std::cout << finalParameters << std::endl;
collector.Report();
// TransformType::Pointer finalTransform=TransformType::New();
// finalTransform->SetParameters(finalParameters);
transform->SetParameters(finalParameters);
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
ResampleFilterType::Pointer resample = ResampleFilterType::New();
resample->SetTransform( transform );
resample->SetInput( movingImageReader->GetOutput() );
FixedImageType::Pointer fixedImage=fixedImageReader->GetOutput();
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
resample->SetDefaultPixelValue( 100 );
typedef unsigned char OutputPixelType;
typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;
typedef itk::CastImageFilter<
FixedImageType,
OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName( argv[3] );
caster->SetInput( resample->GetOutput() );
writer->SetInput( caster->GetOutput() );
try
{
writer->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return -1;
}
More information about the Insight-users
mailing list