[Insight-developers] [slicer-devel] bspline registration final cost value using lbfgs wrong
Gao, Yi
gaoyi.cn at gmail.com
Thu Nov 11 17:02:13 EST 2010
Hi Steve,
I compiled the ITK 3.18 using the same config as it is used in Slicer
trunk, namely:
ITK_USE_REVIEW:BOOL=ON
ITK_USE_REVIEW_STATISTICS:BOOL=ON
ITK_USE_OPTIMIZED_REGISTRATION_METHODS:BOOL=ON
ITK_USE_PORTABLE_ROUND:BOOL=ON *** didn't find this
ITK_USE_CENTERED_PIXEL_COORDINATES_CONSISTENTLY:BOOL=ON
ITK_USE_TRANSFORM_IO_FACTORIES:BOOL=ON
BUILD_SHARED_LIBS:BOOL=ON
CMAKE_SKIP_RPATH:BOOL=ON
BUILD_EXAMPLES:BOOL=OFF
BUILD_TESTING:BOOL=OFF
CMAKE_BUILD_TYPE:STRING=RELEASE $::VTK_BUILD_TYPE
ITK_LEGACY_REMOVE:BOOL=ON
And run the bspline registration using MI, the returned final cost
function value is always 0, which is wrong.
With the same config for itk 3.20, the returned value is correct.
So I guess it might be a itk 3.18 problem....
Please let me know if you want me to do further tests.
Thanks!
Best,
yi
On Mon, Nov 8, 2010 at 12:49 PM, Steve Pieper <pieper at bwh.harvard.edu> wrote:
> Hi Yi -
>
> Wow, that's an odd issue - did you compare the versions and configurations
> of ITK? Slicer3 trunk uses ITK 3.18 with the following options:
>
>
> -DITK_USE_REVIEW:BOOL=ON \
> -DITK_USE_REVIEW_STATISTICS:BOOL=ON \
> -DITK_USE_OPTIMIZED_REGISTRATION_METHODS:BOOL=ON \
> -DITK_USE_PORTABLE_ROUND:BOOL=ON \
> -DITK_USE_CENTERED_PIXEL_COORDINATES_CONSISTENTLY:BOOL=ON \
> -DITK_USE_TRANSFORM_IO_FACTORIES:BOOL=ON \
> -DBUILD_SHARED_LIBS:BOOL=ON \
> -DCMAKE_SKIP_RPATH:BOOL=ON \
> -DBUILD_EXAMPLES:BOOL=OFF \
> -DBUILD_TESTING:BOOL=OFF \
> -DCMAKE_BUILD_TYPE:STRING=$::VTK_BUILD_TYPE \
> -DITK_LEGACY_REMOVE:BOOL=ON \
>
> Maybe the optimized registration or one of the other configurations is
> different and might explain the difference?
>
> -Steve
>
> On 11/08/2010 12:39 PM, Gao, Yi wrote:
>>
>> Hi all,
>>
>> I think I got wrong final cost value using lbfgs to do bspline
>> registration over MI. What I did was:
>>
>> 1. locally build Slicer3 using the getbuildtest.tcl, in Slicer3-build
>>
>> 2.
>> Using lbfgs optimization for bspline MI registration, in compiling,
>> point to Slicer3-build for SLICER_DIR, then link to the itk libs in
>> that slicer build. In doing that, I got a 0 when calling
>> optimizer->GetValue(); no matter what images are used.
>>
>> 3. However, the same code, when linked to a separate stand alone ITK
>> build, the final cost is non-zero values.
>>
>> 4. Then, I used the affine MI reg with the
>> regularSteepestDescentOptimizer, linked to either Slicer3-build, or
>> separate ITK, both times the GetValue gives the same reasonable
>> values.
>>
>> 5. So I'm guessing the lbfgs optimizer got some problem in the itk in
>> the Slicer build.
>>
>> The code is just the example MI bpsline reg code:
>>
>> #include "itkImageRegistrationMethod.h"
>> #include "itkMattesMutualInformationImageToImageMetric.h"
>> #include "itkLinearInterpolateImageFunction.h"
>> #include "itkRegularStepGradientDescentOptimizer.h"
>> #include "itkImage.h"
>>
>> #include "itkCenteredTransformInitializer.h"
>>
>> #include "itkAffineTransform.h"
>>
>> #include "itkImageFileReader.h"
>> #include "itkImageFileWriter.h"
>>
>> #include "itkResampleImageFilter.h"
>> #include "itkCastImageFilter.h"
>> #include "itkSubtractImageFilter.h"
>> #include "itkRescaleIntensityImageFilter.h"
>>
>> #include "itkCommand.h"
>>
>> #include "itkImageRegistrationMethod.h"
>> #include "itkMattesMutualInformationImageToImageMetric.h"
>> #include "itkLinearInterpolateImageFunction.h"
>> #include "itkImage.h"
>>
>> #include "itkBSplineDeformableTransform.h"
>> #include "itkLBFGSBOptimizer.h"
>>
>> #include "itkResampleImageFilter.h"
>> #include "itkCastImageFilter.h"
>>
>> //std
>> #include<string>
>> #include<csignal>
>>
>>
>> int main(int argc, char** argv)
>> {
>> std::string fixImage(argv[1]);
>> std::string movingImage(argv[2]);
>>
>>
>> typedef float pixel_t;
>> typedef itk::Image<pixel_t, 3> image_t;
>>
>> typedef itk::ImageFileReader< image_t> ImageReaderType;
>>
>>
>> // Read in fixed image
>> ImageReaderType::Pointer reader1 = ImageReaderType::New();
>> reader1->SetFileName(fixImage.c_str());
>>
>> reader1->Update();
>> image_t::Pointer itkFixImage = reader1->GetOutput();
>>
>>
>> // Read in moving image
>> ImageReaderType::Pointer reader2 = ImageReaderType::New();
>> reader2->SetFileName(movingImage.c_str());
>>
>> reader2->Update();
>> image_t::Pointer itkMovingImage = reader2->GetOutput();
>>
>>
>> double finalCostValue = 12.3; // just a random init value
>>
>>
>>
>> {
>> const unsigned int ImageDimension = 3;
>> const unsigned int SpaceDimension = ImageDimension;
>> const unsigned int SplineOrder = 3;
>> typedef double CoordinateRepType;
>>
>> typedef itk::BSplineDeformableTransform< CoordinateRepType,
>> SpaceDimension, SplineOrder> TransformType;
>>
>> /**
>> * set up the optimizer
>> */
>> typedef itk::LBFGSBOptimizer OptimizerType;
>>
>>
>> /**
>> * MI metric
>> */
>> typedef itk::MattesMutualInformationImageToImageMetric< image_t,
>> image_t> MetricType;
>>
>> /**
>> * Linear interpolator
>> */
>> typedef itk::LinearInterpolateImageFunction< image_t, double>
>> InterpolatorType;
>>
>>
>> /**
>> * Registratio type
>> */
>> typedef itk::ImageRegistrationMethod< image_t, image_t>
>> RegistrationType;
>>
>> MetricType::Pointer metric = MetricType::New();
>> OptimizerType::Pointer optimizer = OptimizerType::New();
>> InterpolatorType::Pointer interpolator = InterpolatorType::New();
>> RegistrationType::Pointer registration = RegistrationType::New();
>>
>>
>> registration->SetMetric( metric );
>> registration->SetOptimizer( optimizer );
>> registration->SetInterpolator( interpolator );
>>
>> TransformType::Pointer transform = TransformType::New();
>> registration->SetTransform( transform );
>>
>> registration->SetFixedImage( itkFixImage );
>> registration->SetMovingImage( itkMovingImage );
>>
>> image_t::RegionType fixedRegion = itkFixImage->GetBufferedRegion();
>> registration->SetFixedImageRegion( fixedRegion );
>>
>>
>> typedef TransformType::RegionType RegionType;
>> RegionType bsplineRegion;
>> RegionType::SizeType gridSizeOnImage;
>> RegionType::SizeType gridBorderSize;
>> RegionType::SizeType totalGridSize;
>>
>> int gridNum = 10;
>> gridSizeOnImage.Fill( gridNum );
>> gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1
>> lower, 2 upper )
>> totalGridSize = gridSizeOnImage + gridBorderSize;
>>
>> bsplineRegion.SetSize( totalGridSize );
>>
>> typedef TransformType::SpacingType SpacingType;
>> SpacingType spacing = itkFixImage->GetSpacing();
>>
>> typedef TransformType::OriginType OriginType;
>> OriginType origin = itkFixImage->GetOrigin();
>>
>> image_t::SizeType fixedImageSize = fixedRegion.GetSize();
>>
>> for(unsigned int r=0; r<ImageDimension; r++)
>> {
>> spacing[r] *= static_cast<double>(fixedImageSize[r] - 1) /
>> static_cast<double>(gridSizeOnImage[r] - 1);
>> }
>>
>> image_t::DirectionType gridDirection = itkFixImage->GetDirection();
>> SpacingType gridOriginOffset = gridDirection * spacing;
>>
>> OriginType gridOrigin = origin - gridOriginOffset;
>>
>> transform->SetGridSpacing( spacing );
>> transform->SetGridOrigin( gridOrigin );
>> transform->SetGridRegion( bsplineRegion );
>> transform->SetGridDirection( gridDirection );
>>
>>
>> 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;
>>
>> 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+1 );
>> optimizer->SetProjectedGradientTolerance( 1e-7 );
>> optimizer->SetMaximumNumberOfIterations( 20 );
>> optimizer->SetMaximumNumberOfEvaluations( 500 );
>> optimizer->SetMaximumNumberOfCorrections( 12 );
>>
>> //optimizer->MaximizeOn(); // MI needs to be maximized to align images
>> optimizer->MinimizeOn();
>>
>>
>> metric->SetNumberOfHistogramBins( 50 );
>> const unsigned int numberOfSamples =
>> static_cast<unsigned int>( fixedRegion.GetNumberOfPixels() / 10.0 );
>>
>> metric->SetNumberOfSpatialSamples( numberOfSamples );
>> //metric->SetNumberOfSpatialSamples( 50000 );
>> metric->ReinitializeSeed( 76926294 );
>> metric->SetUseCachingOfBSplineWeights( true );
>>
>>
>> //std::cout<< std::endl<< "Starting Registration"<< std::endl;
>>
>> try
>> {
>> registration->StartRegistration();
>> // std::cout<< "Optimizer stop condition ="
>> //<<
>> registration->GetOptimizer()->GetStopConditionDescription()
>> //<< std::endl;
>> }
>> catch( itk::ExceptionObject& err )
>> {
>> std::cerr<< "ExceptionObject caught !"<< std::endl;
>> std::cerr<< err<< std::endl;
>> raise(SIGABRT);
>> }
>>
>> OptimizerType::ParametersType finalParameters =
>> registration->GetLastTransformParameters();
>>
>> transform->SetParameters( finalParameters );
>>
>> typedef itk::ResampleImageFilter< image_t, image_t>
>> ResampleFilterType;
>>
>>
>> // record final cost function value
>> finalCostValue = optimizer->GetValue();
>>
>> //tst
>> std::cout<<"finalCostValue ="<<finalCostValue<<std::endl;
>> //tst//
>>
>> } //reg
>>
>>
>>
>> return 0;
>> }
>>
>>
>>
>> Thank you for any hint on what might be happening!
>>
>>
>> Best,
>> yi
>> _______________________________________________
>> slicer-devel mailing list
>> slicer-devel at bwh.harvard.edu
>> http://massmail.spl.harvard.edu/mailman/listinfo/slicer-devel
>> To unsubscribe: send email to
>> slicer-devel-request at massmail.spl.harvard.edu with unsubscribe as the
>> subject
>
More information about the Insight-developers
mailing list