[Insight-users] How to deform segmentation to compare registration result??

Fri Nov 13 08:13:24 EST 2009

Hi, I just started to use ITK and am comparing the performance of the
registration algorithms.
When I deform a brain label segmentation, I have many artifacts in the
segmentation.

The procedure is as follows.
Let’s say I have 5, 3D brain images with ground truth segmentations.
So I register all the brains to one specific brain then using the
deformation field found, I transform all the segmentation to the same
specific brain. Then I measure how well the each segmented region is
overlapping each other.
But problem I have now is that after I transform the segmentation, I observe
lots of salt and pepper like noise here there inside each region.

One more question is...how long does 3D BSplineTransform usually take??? I
have 7*7*7 landmarks and it takes around 2 minutes in intel 2.4 mobile cpu
2GB memory.

Thank you!!!

Here is the code for Bspline deformation that I am using .......

    const unsigned int ImageDimension = 3;

    typedef itk::ImageFileWriter< ITKIMAGE3D > MovingWriterType;
    MovingWriterType::Pointer movingWriter = MovingWriterType::New();

    typedef itk::ResampleImageFilter< ITKIMAGE3D,
    ITKIMAGE3D > FilterType;
        ITKIMAGE3D::Pointer fixedImage=itkimage_src;

    FilterType::Pointer resampler = FilterType::New();
    typedef itk::LinearInterpolateImageFunction<
    ITKIMAGE3D, double > InterpolatorType;
    InterpolatorType::Pointer interpolator = InterpolatorType::New();
    resampler->SetInterpolator( interpolator );
        ITKIMAGE3D::SpacingType fixedSpacing = fixedImage->GetSpacing();
        ITKIMAGE3D::PointType fixedOrigin = fixedImage->GetOrigin();
    ITKIMAGE3D::DirectionType fixedDirection  = fixedImage->GetDirection();

    resampler->SetOutputSpacing( fixedSpacing );
    resampler->SetOutputOrigin( fixedOrigin );
    resampler->SetOutputDirection( fixedDirection );

        ITKIMAGE3D::RegionType fixedRegion =
fixedImage->GetLargestPossibleRegion();
    ITKIMAGE3D::SizeType fixedSize = fixedRegion.GetSize();
    resampler->SetSize( fixedSize );
    resampler->SetOutputStartIndex( fixedRegion.GetIndex() );
    resampler->SetInput(itkimage_dst );
    movingWriter->SetInput( resampler->GetOutput() );
    movingWriter->SetFileName(outputPath);
    /*
     We instantiate now the type of the BSplineDeformableTransform using as
template param- 
     eters the type for coordinates representation, the dimension of the
space, and the order of the
     B-spline.*/

    TransformType3D::Pointer bsplineTransform = TransformType3D::New();
    /*Since we are using a B-spline of order 3, the coverage of the BSpling
grid should exceed by one
     the spatial extent of the image on the lower region of image indices,
and by two grid points on
     the upper region of image indices. We choose here to use a 8
     × 8 B-spline grid, from which only
     a 5 by 5 sub-grid will be covering the input image. If we use an input
image of size 500 by 500
     pixels, and pixel spacing 2.0 × 2.0 then we need the 5 × 5 B-spline
grid to cover a physical
     extent of 1000 × 1000 mm. This can be achieved by setting the pixel
spacing of the B-spline
     grid to 250.0 × 250.0 mm. The origin of the B-spline grid must be set
at one grid position away
     from the origin of the desired output image. All this is done with the
following lines of code. */

    typedef TransformType3D::RegionType RegionType;
    RegionType bsplineRegion; //itk::BSplineDeformableTransform::RegionType
    RegionType::SizeType size;
//itk::BSplineDeformableTransform::RegionType::SizeType
    const unsigned int numberOfGridNodesOutsideTheImageSupport = 3;
    const unsigned int numberOfGridNodesInsideTheImageSupportH = nH;
    const unsigned int numberOfGridNodesInsideTheImageSupportW = nW;
    const unsigned int numberOfGridNodesInsideTheImageSupportD = nD;
    const unsigned int numberOfGridNodesH =
    numberOfGridNodesInsideTheImageSupportH +
    numberOfGridNodesOutsideTheImageSupport;
    const unsigned int numberOfGridNodesW =
    numberOfGridNodesInsideTheImageSupportW +
    numberOfGridNodesOutsideTheImageSupport;
    const unsigned int numberOfGridNodesD =
    numberOfGridNodesInsideTheImageSupportD +
    numberOfGridNodesOutsideTheImageSupport;

    //We need to modify this size[0]=m_nW size[1]=m_nH
    //itk::BSplineDeformableTransform::RegionType::SizeType.Fill
    size[0]=numberOfGridNodesW;
    size[1]=numberOfGridNodesH;
    size[2]=numberOfGridNodesD;
    bsplineRegion.SetSize( size );
//itk::BSplineDeformableTransform::RegionType::SetSize
    typedef TransformType3D::SpacingType SpacingType;
//itk::BSplineDeformableTransform::SpacingType
    SpacingType spacing;
    spacing[0] = nHorSpace*fixedSpacing[0];//floor( fixedSpacing[0] *
(fixedSize[0] - 1) / numberOfGridCellsW ); //convert pixel spacing to
metric....
    spacing[1] = nVerSpace*fixedSpacing[1];//floor( fixedSpacing[1] *
(fixedSize[1] - 1) / numberOfGridCellsH );
    spacing[2] = nSliSpace*fixedSpacing[2];//floor( fixedSpacing[1] *
(fixedSize[1] - 1) / numberOfGridCellsH );

    typedef TransformType3D::OriginType OriginType;
//itk::BSplineDeformableTransform::OriginType
    OriginType origin;
    origin = fixedOrigin - fixedDirection*spacing; //due to boundary...

    bsplineTransform->SetGridSpacing( spacing );
//itk::BSplineDeformableTransform->SetGridSpacing
    bsplineTransform->SetGridOrigin( origin );
//itk::BSplineDeformableTransform->SetGridOrigin
    bsplineTransform->SetGridRegion( bsplineRegion
);//itk::BSplineDeformableTransform->SetGridRegion
    bsplineTransform->SetGridDirection( fixedDirection );

    typedef TransformType3D::ParametersType ParametersType;
//itk::BSplineDeformableTransform::ParametersType
    const unsigned int numberOfParameters =
    bsplineTransform->GetNumberOfParameters();
//itk::BSplineDeformableTransform->GetNumberOfParameters();
    const unsigned int numberOfNodes = numberOfParameters / ImageDimension;
    ParametersType parameters( numberOfParameters );

    /*
     The B-spline grid should now be fed with coeﬁcients at each node. Since
this is a two dimen-
     sional grid, each node should receive two coefﬁcients. Each coefﬁcient
pair is representing a
     displacement vector at this node. The coefﬁcients can be passed to the
B-spline in the form of
     an array where the ﬁrst set of elements are the ﬁrst component of the
displacements for all the
     nodes, and the second set of elemets is formed by the second component
of the displacements
     for all the nodes.
     In this example we read such displacements from a ﬁle, but for
convinience we have written this
     ﬁle using the pairs of (x, y) displacement for every node. The elements
read from the ﬁle should
     therefore be reorganized when assigned to the elements of the array. We
do this by storing all
     the odd elements from the ﬁle in the ﬁrst block of the array, and all
the even elements from the
     ﬁle in the second block of the array. Finally the array is passed to
the B-spline transform using
     the SetParameters(). */

    int nPos=0;
        for(int z=-1;z<nD+2;z++)
    {
        for(int y=-1;y<nH+2;y++)
        {
            for(int x=-1;x<nW+2;x++)
            {
                if(z>=0 &&z< nD && y>=0 && y< nH && x>=0 && x<nW)
                {
                    parameters[nPos]=disVector[z*nW*nH+y*nW+x].x;

parameters[nPos+numberOfNodes]=disVector[z*nW*nH+y*nW+x].y;

parameters[nPos+numberOfNodes*2]=disVector[z*nW*nH+y*nW+x].z;
                }
                else
                {
                    parameters[nPos]=0;//disVector[y*nW+x].x;
                    parameters[nPos+numberOfNodes]=0;//disVector[y*nW+x].y;
                    parameters[nPos+numberOfNodes*2]=0;
                }
                nPos++;
            }
        }
    }
    bsplineTransform->SetParameters( parameters );

    resampler->SetDefaultPixelValue(0);
    resampler->SetTransform( bsplineTransform );
    try 
    { 
        movingWriter->Update();
    } 
    catch( itk::ExceptionObject & excp )
    { 
        std::cerr << "Exception thrown " << std::endl;
        std::cerr << excp << std::endl;
        return EXIT_FAILURE;
    }

    //We just transform segmentation image to evaluate.....
    if(m_strSegPath!="")
    {
        if(nLevel==0)
        {//input is fixed segmentation image...
            ITK_IMG_READER_TYPE3D::Pointer
seg_reader=ITK_IMG_READER_TYPE3D::New();
            seg_reader->SetFileName(m_strSegPath);

            //we first compute maximum and minimum of the label...and
remember it....
            typedef itk::MinimumMaximumImageCalculator<ITKIMAGE3D> calType;
            calType::Pointer cal=calType::New();
            cal->SetImage(seg_reader->GetOutput());
            seg_reader->Update();
                        cal->ComputeMaximum();
                        cal->ComputeMinimum();
            ITKIMAGE3D::PixelType maxval=cal->GetMaximum();
            ITKIMAGE3D::PixelType minval=cal->GetMinimum();

            typedef itk::RescaleIntensityImageFilter<
            ITKIMAGE3D, ITKIMAGE3D> RescaleFilterType;
            RescaleFilterType::Pointer enlarge_rescaler =
RescaleFilterType::New();
            //This rescaler was used to avoid interpolation ambiguity
between adjacent labels.
            // in hope for avoiding salt and pepper like noise....
            enlarge_rescaler->SetOutputMinimum(   0 );
            enlarge_rescaler->SetOutputMaximum( 30000 );
            enlarge_rescaler->SetInput(seg_reader->GetOutput());

            resampler->SetInput(enlarge_rescaler->GetOutput());
            ITK_IMG_WRITER_TYPE3D::Pointer
seg_writer=ITK_IMG_WRITER_TYPE3D::New();
            QString strFileName;

strFileName.sprintf("%s%s_seg_deformed%.2d.hdr",qPrintable(m_strOutPath),qPr
intable(m_strFileName),nLevel);
            seg_writer->SetFileName(strFileName);

            RescaleFilterType::Pointer shorten_rescaler =
RescaleFilterType::New();

            shorten_rescaler->SetOutputMinimum(minval);
            shorten_rescaler->SetOutputMaximum(maxval);
            shorten_rescaler->SetInput(resampler->GetOutput());
            seg_writer->SetInput(shorten_rescaler->GetOutput());
            seg_writer->Update();
        }
        else
        {//input is transformed segmentation image....
            QString strFileName;

strFileName.sprintf("%s%s_seg_deformed%.2d.hdr",qPrintable(m_strOutPath),qPr
intable(m_strFileName),nLevel-1);

            ITK_IMG_READER_TYPE3D::Pointer
seg_reader=ITK_IMG_READER_TYPE3D::New();

            //we first compute maximum and minimum of the label...and
remember it....
            typedef itk::MinimumMaximumImageCalculator<ITKIMAGE3D> calType;
            calType::Pointer cal=calType::New();
            cal->SetImage(seg_reader->GetOutput());
            seg_reader->Update();
            cal->ComputeMaximum();
            cal->ComputeMinimum();
            ITKIMAGE3D::PixelType maxval=cal->GetMaximum();
            ITKIMAGE3D::PixelType minval=cal->GetMinimum();

            typedef itk::RescaleIntensityImageFilter<
            ITKIMAGE3D, ITKIMAGE3D> RescaleFilterType;
            RescaleFilterType::Pointer enlarge_rescaler =
RescaleFilterType::New();

            enlarge_rescaler->SetOutputMinimum(   0 );
                        enlarge_rescaler->SetOutputMaximum( 30000 );
            enlarge_rescaler->SetInput(seg_reader->GetOutput());

            resampler->SetInput(enlarge_rescaler->GetOutput());

            seg_reader->SetFileName(strFileName);

            ITK_IMG_WRITER_TYPE3D::Pointer
seg_writer=ITK_IMG_WRITER_TYPE3D::New();

strFileName.sprintf("%s%s_seg_deformed%.2d.hdr",qPrintable(m_strOutPath),qPr
intable(m_strFileName),nLevel);
            seg_writer->SetFileName(strFileName);

            RescaleFilterType::Pointer shorten_rescaler =
RescaleFilterType::New();

            shorten_rescaler->SetOutputMinimum(minval);
            shorten_rescaler->SetOutputMaximum(maxval);
            shorten_rescaler->SetInput(resampler->GetOutput());
            seg_writer->SetInput(shorten_rescaler->GetOutput());
            seg_writer->Update();
        }
    } 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091113/6e16f024/attachment-0001.htm>