itkbsplineinterpolateimagefunctiontest.cxx

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkBSplineInterpolateImageFunctionTest.cxx,v $
  Language:  C++
  Date:      $Date: 2007-08-10 14:34:01 $
  Version:   $Revision: 1.10 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

  Portions of this code are covered under the VTK copyright.
  See VTKCopyright.txt or http://www.kitware.com/VTKCopyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include <iostream>

#include "itkImage.h"
#include "itkSize.h"
#include "itkBSplineInterpolateImageFunction.h"
#include "itkImageRegionIterator.h"


#include "vnl/vnl_math.h"

  typedef double InputPixelType;
  typedef double CoordRepType;

// Set up for 1D Images
  enum { ImageDimension1D = 1 };

  typedef itk::Image< InputPixelType, ImageDimension1D > ImageType1D;
  typedef ImageType1D::Pointer ImageTypePtr1D;
  typedef ImageType1D::SizeType SizeType1D;
  typedef itk::BSplineInterpolateImageFunction<ImageType1D,CoordRepType> InterpolatorType1D;
//  typedef InterpolatorType1D::IndexType                 IndexType1D;
  typedef InterpolatorType1D::PointType                 PointType1D;
  typedef InterpolatorType1D::ContinuousIndexType       ContinuousIndexType1D;

  void set1DInterpData(ImageType1D::Pointer);

  // Set up for 2D Images
  enum { ImageDimension2D = 2 };

  typedef itk::Image< InputPixelType, ImageDimension2D > ImageType2D;
  typedef ImageType2D::Pointer ImageTypePtr2D;
  typedef ImageType2D::SizeType SizeType2D;
  typedef itk::BSplineInterpolateImageFunction<ImageType2D,CoordRepType> InterpolatorType2D;
//  typedef InterpolatorType2D::IndexType                 IndexType2D;
  typedef InterpolatorType2D::PointType                 PointType2D;
  typedef InterpolatorType2D::ContinuousIndexType       ContinuousIndexType2D;

  void set2DInterpData(ImageType2D::Pointer);

  // Set up for 3D Images
  enum { ImageDimension3D = 3 };

  typedef itk::Image< InputPixelType, ImageDimension3D > ImageType3D;
  typedef ImageType3D::Pointer ImageTypePtr3D;
  typedef ImageType3D::SizeType SizeType3D;
  typedef itk::BSplineInterpolateImageFunction<ImageType3D,CoordRepType> InterpolatorType3D;
  typedef InterpolatorType3D::IndexType                 IndexType3D;
  typedef InterpolatorType3D::PointType                 PointType3D;
  typedef InterpolatorType3D::ContinuousIndexType       ContinuousIndexType3D;

  typedef itk::Image< unsigned int, ImageDimension3D > ImageIntegerType3D;
  typedef itk::BSplineInterpolateImageFunction<ImageIntegerType3D,CoordRepType> InterpolatorIntegerType3D;
    typedef InterpolatorIntegerType3D::IndexType                 IndexIntegerType3D;
  typedef InterpolatorIntegerType3D::PointType                 PointIntegerType3D;
  typedef InterpolatorIntegerType3D::ContinuousIndexType       ContinuousIntegerIndexType3D;

//  template <class TImage>
//  void set3DInterpData(TImage::Pointer);
  void set3DDerivativeData(ImageType3D::Pointer);

template<class TImage>
void set3DInterpData(typename TImage::Pointer imgPtr)
{
  SizeType3D size = {{80,40,30}};

  /* Allocate a simple test image */
  typename TImage::RegionType region;
  region.SetSize( size );

  imgPtr->SetLargestPossibleRegion( region );
  imgPtr->SetBufferedRegion( region );
  imgPtr->Allocate();


  /* Set origin and spacing of physical coordinates */

  /* Initialize the image contents */
  IndexType3D index;
  for (unsigned int slice = 0; slice < size[2]; slice++) {
      index[2] = slice;
      for (unsigned int row = 0; row < size[1]; row++) {
          index[1] = row;
          for (unsigned int col = 0; col < size[0]; col++) {
              index[0] = col;
              imgPtr->SetPixel(index, slice+row+col);
          }
      }
  }
}

/**
 * Test a geometric point. Returns true if test has passed,
 * returns false otherwise
 */
template <class TInterpolator, class PointType>
bool TestGeometricPoint(
const TInterpolator * interp,
const PointType& point,
bool isInside,
double trueValue )
{

  std::cout << " Point: " << point;

  bool bvalue = interp->IsInsideBuffer( point );
  std::cout << " Inside: " << bvalue << " ";

  if( bvalue != isInside )
    {
    std::cout << "*** Error: inside should be " << isInside << std::endl;
    return false;
    }

  if( isInside )
    {
    double value = interp->Evaluate( point );
    std::cout << " Value: " << value;

    if( vnl_math_abs( value - trueValue ) > 1e-9 )
      {
      std::cout << "*** Error: value should be " << trueValue << std::endl;
      return false;
      }
    }

  std::cout << std::endl;
  return true;

}

/**
 * Test a continuous index. Returns true if test has passed,
 * returns false otherwise
 */
template<class TInterpolator, class ContinuousIndexType>
bool TestContinuousIndex(
const TInterpolator * interp,
const ContinuousIndexType& index,
bool isInside,
double trueValue )
{

  std::cout << " Index: " << index;

  bool bvalue = interp->IsInsideBuffer( index );
  std::cout << " Inside: " << bvalue;

  if( bvalue != isInside )
    {
    std::cout << "*** Error: inside should be " << isInside << std::endl;
    return false;
    }

  if( isInside )
    {
    double value = interp->EvaluateAtContinuousIndex( index );
    std::cout << " Value: " << value;

    if( vnl_math_abs( value - trueValue ) > 1e-4 )
      {
      std::cout << "*** Error: value should be " << trueValue << std::endl;
      return false;
      }
    }

  std::cout << std::endl;
  return true;

}
/**
 * Test a continuous index Derivative. Returns true if test has passed,
 * returns false otherwise
 */
template<class TInterpolator, class ContinuousIndexType>
bool TestContinuousIndexDerivative(
const TInterpolator * interp,
const ContinuousIndexType& index,
bool isInside,
double * trueValue )
{

  std::cout << " Index: " << index;

  bool bvalue = interp->IsInsideBuffer( index );
  std::cout << " Inside: " << bvalue << "\n";

  if( bvalue != isInside )
    {
    std::cout << "*** Error: inside should be " << isInside << std::endl;
    return false;
    }

  if( isInside )
    {
    typename TInterpolator::CovariantVectorType value;
    double value2 = interp->EvaluateAtContinuousIndex( index );
    std::cout << "Interpolated Value: " << value2 << "\n";
    value = interp->EvaluateDerivativeAtContinuousIndex( index );
    std::cout << " Value: "; 
    for (int i=0; i < ImageDimension3D; i++)
        {
        if (i != 0)
          {
          std::cout << ", ";
          }
        std::cout << value[i] ;
        if ( vnl_math_abs( value[i] - trueValue[i] ) > 1e-4 )
          {
          std::cout << "*** Error: value should be " << trueValue[i] << std::endl;
          return false;
          }
        }

    }

  std::cout << std::endl;
  return true;

}


// Run a series of tests to validate the 1D 
// cubic spline implementation.
int test1DCubicSpline()
{
  int flag = 0;

  // Allocate a simple test image 
  ImageTypePtr1D image = ImageType1D::New();

  set1DInterpData(image);

  // Set origin and spacing of physical coordinates 
  double origin [] = { 0.5  };
  double spacing[] = { 0.1  };
  image->SetOrigin(origin);
  image->SetSpacing(spacing);

  // Create and initialize the interpolator 
  InterpolatorType1D::Pointer interp = InterpolatorType1D::New();
//  interp->SetSplineOrder(1);
  interp->SetInputImage(image);
  interp->Print( std::cout );

  // Test evaluation at continuous indices and corresponding
  //gemetric points 
  std::cout << "Testing 1D Cubic B-Spline:\n";
  std::cout << "Evaluate at: " << std::endl;
  ContinuousIndexType1D cindex;
  PointType1D point;
  bool passed;

  // These values test 1) near border, 
  //    2) inside
  //    3) integer value
  //    4) outside image
#define NPOINTS 4  // number of points 
  double darray1[NPOINTS] = {1.4, 8.9, 10.0, 40.0};
  double truth[NPOINTS] = {334.41265437584, 18.158173426944, 4.0000, 0};
  bool b_Inside[NPOINTS] = {true, true, true, false};

  // an integer position inside the image
  for (int ii=0; ii < NPOINTS; ii++)
    {
  
    cindex = ContinuousIndexType1D(&darray1[ii]);
    passed = TestContinuousIndex<InterpolatorType1D, ContinuousIndexType1D >( interp, cindex, b_Inside[ii], truth[ii] );
  
    if( !passed ) flag += 1;
  
    image->TransformContinuousIndexToPhysicalPoint( cindex, point );
    passed = TestGeometricPoint<InterpolatorType1D, PointType1D>( interp, point, b_Inside[ii], truth[ii]  );
  
    if( !passed ) flag += 1;
    }

  return (flag);
}


int test2DSpline()
{
  int flag = 0;

  /* Allocate a simple test image */
  ImageTypePtr2D image = ImageType2D::New();

  set2DInterpData(image);

  /* Set origin and spacing of physical coordinates */
  double origin [] = { 0.5, 1.0 };
  double spacing[] = { 0.1, 0.5  };
  image->SetOrigin(origin);
  image->SetSpacing(spacing);

  /* Create and initialize the interpolator */
  for (unsigned int splineOrder = 0; splineOrder<=5; splineOrder++)
    {
    InterpolatorType2D::Pointer interp = InterpolatorType2D::New();
    interp->SetSplineOrder(splineOrder);

    std::cout << "SplineOrder: " << interp->GetSplineOrder() << std::endl;

    interp->SetInputImage(image);
    interp->Print( std::cout );

    /* Test evaluation at continuous indices and corresponding
    gemetric points */
    std::cout << "Testing 2D B-Spline of Order "<< splineOrder << ":\n";
    std::cout << "Evaluate at: " << std::endl;
    ContinuousIndexType2D cindex;
    PointType2D point;
    bool passed;

    // These values test 1) near border, 
    //    2) inside
    //    3) integer value
    //    4) outside image
#define NPOINTS2 4  // number of points 

    double darray1[NPOINTS2][2] = {{0.1, 0.2}, {3.4, 5.8}, {4.0, 6.0}, { 2.1, 8.0}};
    double truth[NPOINTS2][6] = {{154.5, 140.14, 151.86429192392, 151.650316034, 151.865916515, 151.882483111},
        { 0, 13.84, 22.688125812495, 22.411473093, 22.606968306, 22.908345604},
        { 36.2, 36.2, 36.2, 36.2, 36.2, 36.2 },
        {0, 0, 0,0,0,0}};
    bool b_Inside[NPOINTS2] = {true, true, true, false};
   // double darray1[2];

    // an integer position inside the image
    for (int ii=0; ii < NPOINTS2; ii++)
      {
     // darray1[0] = darray[ii][0];
     // darray1[1] = darray[ii][1];
      cindex = ContinuousIndexType2D(&darray1[ii][0]);
      passed = TestContinuousIndex<InterpolatorType2D, ContinuousIndexType2D >( interp, cindex, b_Inside[ii], truth[ii][splineOrder] );
  
      if( !passed ) flag += 1;
  
      image->TransformContinuousIndexToPhysicalPoint( cindex, point );