[Insight-developers] [slicer-devel] bspline registration final cost value using lbfgs wrong
Stephen Aylward
stephen.aylward at kitware.com
Fri Nov 12 11:04:50 EST 2010
I think it is a great idea.
This is the "final" version of the ITKv3 series and should be
supported for a long time.
s
On Fri, Nov 12, 2010 at 9:16 AM, Steve Pieper <pieper at bwh.harvard.edu> wrote:
> Hi Yi -
>
> Interesting - sounds like a fix was added for ITK 3.20.
>
> Are there are objections to a move of the slicer3 trunk to use 3.20?
>
> -Steve
>
> On 11/11/2010 05:02 PM, Gao, Yi wrote:
>> 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
>>>
> _______________________________________________
> 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
>
--
==============================
Stephen R. Aylward, Ph.D.
Director of Medical Imaging Research
Kitware, Inc. - North Carolina Office
http://www.kitware.com
stephen.aylward (Skype)
(919) 969-6990 x300
More information about the Insight-developers
mailing list