[Insight-users] bspline deformable registration problem

[Samuel Gerber]

Hi

I am trying to adapt the DeformableRegistration6 example so that i can 
use a Mask for the MovingImage and reduce computation by setting the 
FixedImageRegion. But the registration doesn't work proberly anymore as 
soon as I set a FixedImageRegion. The error I get is:

Starting Registration with low resolution transform
 F = 0, GNORM = 0
*************************************************
   I   NFN    FUNC        GNORM       STEPLENGTH
THE SEARCH DIRECTION IS NOT A DESCENT DIRECTIONIFLAG= -1. LINE SEARCH 
FAILED. SEE DOCUMENTATION OF ROUTINE MCSRCH ERROR RETURN  OF LINE SEAR
CH: INFO=0 POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE  INCORRECT OR 
INCORRECT TOLERANCESvnl_lbfgs: ** EEEK **
ExceptionObject caught !

I think i don't set the GridOrigin, GridSpacing and GridRegion of the 
transformation right. I am not sure if I fully understand what they are 
for. But to my understanding the code below should setup the Grid right. 
Or is there something wrong? Any other hints what might be wong?



      //Image Region to use for registration
      //Find bounding "box" of mask image
      FixedImageType::IndexType minIndex;
        minIndex.Fill( vnl_numeric_traits< 
FixedImageType::IndexType::IndexValueType >::maxval );
        FixedImageType::IndexType maxIndex;
        maxIndex.Fill(0);

        typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> 
MaskImageIterator;
        MaskImageIterator maskIt(maskImage, 
maskImage->GetLargestPossibleRegion());
        for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt){
            if(maskIt.Get() != 0){
                FixedImageType::IndexType tmp = maskIt.GetIndex();
                for(int i = 0; i < ImageDimension; i++){
                    if(maxIndex[i] < tmp[i]){
                        maxIndex[i] = tmp[i];
                    }
                    if(minIndex[i] > tmp[i]){
                        minIndex[i] = tmp[i];
                    }
                }
            }
        }


        FixedImageType::RegionType::SizeType size;
        for(int i = 0; i < ImageDimension; i++){
            size[i] = 5*(maxIndex[i]    - minIndex[i]);
            minIndex[i] = minIndex[i] - 2*( maxIndex[i] - minIndex[i] );
        }

        FixedImageType::RegionType maxRegion = 
fixedImage->GetLargestPossibleRegion();
        FixedImageType::RegionType::SizeType maxSize = maxRegion.GetSize();
        for(int i = 0; i < ImageDimension; i++){
            if(minIndex[i] < 0){
                minIndex[i] = 0;
            }
            if(minIndex[i] + size[i] > maxSize[i]){
                size[i] = maxSize[i] - minIndex[i];
            }
        }

        FixedImageType::RegionType fixedRegion(minIndex, size);
        std::cout << fixedRegion << std::endl;
        registration->SetFixedImageRegion( fixedRegion );


        //Moving Image mask for metric
        ImageMaskSpatialObjectType::Pointer mask = 
ImageMaskSpatialObjectType::New();
        mask->SetImage(maskImage);
        metric->SetMovingImageMask(mask);


        //Low resolution transform
        typedef TransformType::RegionType RegionType;
        RegionType::SizeType   gridBorderSize;
        RegionType::SizeType   totalGridSize;

        gridBorderSize.Fill( 3 );    // Border for spline order = 3 ( 1 
lower, 2 upper )

        RegionType::SizeType   gridLowSizeOnImage;
        gridLowSizeOnImage.Fill( 5 );
        totalGridSize = gridLowSizeOnImage + gridBorderSize;

        RegionType bsplineRegion;
        bsplineRegion.SetSize( totalGridSize );

        typedef TransformType::SpacingType SpacingType;
        SpacingType spacingLow = fixedImage->GetSpacing();
       
        typedef TransformType::OriginType OriginType;
        OriginType originLow;
        fixedImage->TransformIndexToPhysicalPoint(minIndex, originLow);

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

        for(unsigned int r=0; r<ImageDimension; r++){
            spacingLow[r] = floor( 
static_cast<double>(fixedImageSize[r])  / 
static_cast<double>(gridLowSizeOnImage[r]) );
            originLow[r]  -=  spacingLow[r];
        }
       

        transformLow->SetGridSpacing( spacingLow );
        transformLow->SetGridOrigin( originLow );
        transformLow->SetGridRegion( bsplineRegion );


Kind regards
Samuel