itkbsplineinterpolateimagefunctiontest.cxx

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

      passed = TestGeometricPoint<InterpolatorType2D, PointType2D>( interp, point, b_Inside[ii], truth[ii][splineOrder ]  );
  
      if( !passed ) flag += 1;
      }
    }  // end of splineOrder

  return (flag);
}

int test3DSpline()
{
  int flag = 0;

  /* Allocate a simple test image */
  ImageTypePtr3D image = ImageType3D::New();

  set3DInterpData<ImageType3D> (image);

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

  /* Create and initialize the interpolator */
  for (int splineOrder = 2; splineOrder<=5; splineOrder++)
    {
    InterpolatorType3D::Pointer interp = InterpolatorType3D::New();
    interp->SetSplineOrder(splineOrder);
    interp->SetInputImage(image);
    interp->Print( std::cout );

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

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

    double darray1[NPOINTS3][ImageDimension3D] = {{0.1, 20.1, 28.4}, {21.58, 34.5, 17.2 }, {10, 20, 12}, { 15, 20.2, 31}};
    double truth[NPOINTS3][4] = {{48.621593795, 48.651173138, 48.656914878, 48.662256571},
        { 73.280126903, 73.280816965, 73.282780615, 73.285315943},
        { 42.0, 42.0, 42.0, 42.0},
        {0,0,0,0}};
    bool b_Inside[NPOINTS3] = {true, true, true, false};
   // double darray1[2];

    // an integer position inside the image
    for (int ii=0; ii < NPOINTS3; ii++)
      {
     // darray1[0] = darray[ii][0];
     // darray1[1] = darray[ii][1];
      cindex = ContinuousIndexType3D(&darray1[ii][0]);
      passed = TestContinuousIndex<InterpolatorType3D, ContinuousIndexType3D >( interp, cindex, b_Inside[ii], truth[ii][splineOrder -2] );
  
      if( !passed ) flag += 1;
  
      image->TransformContinuousIndexToPhysicalPoint( cindex, point );
      passed = TestGeometricPoint<InterpolatorType3D, PointType3D>( interp, point, b_Inside[ii], truth[ii][splineOrder -2]  );
  
      if( !passed ) flag += 1;
      }
    }  // end of splineOrder

  return (flag);
}

int test3DSplineDerivative()
{
  int flag = 0;

  /* Allocate a simple test image */
  ImageTypePtr3D image = ImageType3D::New();

  set3DDerivativeData(image);

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

  /* Create and initialize the interpolator */
  for (int splineOrder = 1; splineOrder<=5; splineOrder++)
    {
    InterpolatorType3D::Pointer interp = InterpolatorType3D::New();
    interp->SetSplineOrder(splineOrder);
    interp->SetInputImage(image);
    interp->Print( std::cout );

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

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

    double darray1[NPOINTS4][ImageDimension3D] = {{25.3,26.8,24.5}, {21.0, 1.4, 0.6}, {18, 31, 10 }, { 4.3, 17.9, 42} };
    // Calculated Truth is: {19.4158,5,-24}, {0.9,5,71.6}, {-7.2, 5, 34}, {0,0,0}
    // TODO: Value near border is way off, is this an algorithm problem?  Also,
    //       Is error for 1st order splines in the expected range?
    double truth[5][NPOINTS4][ImageDimension3D] = { 
      { {23.6,   5,-24}, {0,        5,       72.0},    {-3.0,     5,       32},      {0,0,0} },
      { {19.345, 5,-24}, {0.875,    4.8873,  98.6607}, {-7.525,   5,       34},      {0,0,0} },
      { {19.399, 5,-24}, {0.9,      4.95411, 92.9006}, {-7.2,     5,       33.9999}, {0,0,0} },
      { {19.4164,5,-24}, {0.9,      4.9925,  94.5082}, {-7.2,     5.00044, 33.9976}, {0,0,0} },
      { {19.4223,5,-24}, {0.900157, 5.0544,  93.8607}, {-7.19929, 5.00189, 33.9879}, {0,0,0} }};
    bool b_Inside[NPOINTS4] = {true, true, true, false};
   // double darray1[2];

    // an integer position inside the image
    for (int ii=0; ii < NPOINTS4; ii++)
      {
     // darray1[0] = darray[ii][0];
     // darray1[1] = darray[ii][1];
      cindex = ContinuousIndexType3D(&darray1[ii][0]);
      passed = TestContinuousIndexDerivative<InterpolatorType3D, ContinuousIndexType3D >( interp, cindex, b_Inside[ii], &truth[splineOrder - 1][ii][0] );
  
      if( !passed ) flag += 1;
  
//      interp->ConvertContinuousIndexToPoint( cindex, point );
//      passed = TestGeometricPoint<InterpolatorType3D, PointType3D>( interp, point, b_Inside[ii], truth[ii][splineOrder -2]  );
  
      if( !passed ) flag += 1;
      }
    }  // end of splineOrder

  return (flag);
}

int testInteger3DSpline()
{
  int flag = 0;

  /* Allocate a simple test image */
  ImageIntegerType3D::Pointer image = ImageIntegerType3D::New();

  set3DInterpData<ImageIntegerType3D> (image);

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

  /* Create and initialize the interpolator */
  for (int splineOrder = 2; splineOrder<=5; splineOrder++)
    {
    InterpolatorIntegerType3D::Pointer interp = InterpolatorIntegerType3D::New();
    interp->SetSplineOrder(splineOrder);
    interp->SetInputImage(image);
    interp->Print( std::cout );

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

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

    // Note: the answers should be the same as for the test3DSpline
    double darray1[NPOINTS5][ImageDimension3D] = {{0.1, 20.1, 28.4}, {21.58, 34.5, 17.2 }, {10, 20, 12}, { 15, 20.2, 31}};
    double truth[NPOINTS5][4] = {{48.621593795, 48.651173138, 48.656914878, 48.662256571},
        { 73.280126903, 73.280816965, 73.282780615, 73.285315943},
        { 42.0, 42.0, 42.0, 42.0},
        {0,0,0,0}};
    bool b_Inside[NPOINTS3] = {true, true, true, false};
   // double darray1[2];

    // an integer position inside the image
    for (int ii=0; ii < NPOINTS3; ii++)
      {
     // darray1[0] = darray[ii][0];
     // darray1[1] = darray[ii][1];
      cindex = ContinuousIntegerIndexType3D(&darray1[ii][0]);
      passed = TestContinuousIndex<InterpolatorIntegerType3D, ContinuousIntegerIndexType3D >( interp, cindex, b_Inside[ii], truth[ii][splineOrder -2] );
  
      if( !passed ) flag += 1;
  
      image->TransformContinuousIndexToPhysicalPoint( cindex, point );
      passed = TestGeometricPoint<InterpolatorIntegerType3D, PointIntegerType3D>( interp, point, b_Inside[ii], truth[ii][splineOrder -2]  );
  
      if( !passed ) flag += 1;
      }
    }  // end of splineOrder

  return (flag);
}

int 
itkBSplineInterpolateImageFunctionTest(
    int itkNotUsed(argc),
    char * itkNotUsed(argv) [] )
{
  int flag = 0;           /* Did this test program work? */

  std::cout << "Testing B Spline interpolation methods:\n";

  flag += test1DCubicSpline();

  flag += test2DSpline();

  flag += test3DSpline();

  flag += test3DSplineDerivative();

  flag += testInteger3DSpline();




  /* Return results of test */
  if (flag != 0) {
    std::cout << "*** " << flag << " tests failed" << std::endl;
  
    return EXIT_FAILURE; }
  else {
    std::cout << "All tests successfully passed" << std::endl;
    return EXIT_SUCCESS; }

}

void set1DInterpData(ImageType1D::Pointer imgPtr)
{

  SizeType1D size = {{36}};
  double mydata[36] = {454.0000,  369.4000,  295.2000,  230.8000,  175.6000,  129.0000,   90.4000, 59.2000,   34.8000,   16.6000,    4.0000,   -3.6000,   -6.8000,   -6.2000,
    -2.4000,    4.0000,   12.4000,   22.2000,   32.8000,   43.6000,   54.0000, 63.4000,   71.2000,   76.8000,   79.6000,   79.0000,   74.4000,   65.2000,
    50.8000,   30.6000,    4.0000,  -29.6000,  -70.8000, -120.2000, -178.4000, -246.0000 };

  ImageType1D::RegionType region;
  region.SetSize( size );

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

  typedef itk::ImageRegionIterator<ImageType1D>    InputIterator;

  InputIterator inIter( imgPtr, region );

  int j = 0;
  while( !inIter.IsAtEnd() )
    {
    inIter.Set(mydata[j]);
    ++inIter;
    ++j;
    }

  
}

void set2DInterpData(ImageType2D::Pointer imgPtr)
{
  SizeType2D size = {{7,7}};
  double mydata[ 49 ] = {  154.5000,   82.4000,   30.9000,         0,  -10.3000,         0,   30.9000 ,
    117.0000,   62.4000,   23.4000,         0,   -7.8000,         0,   23.4000 ,
   18.0000,    9.6000,    3.6000,         0,   -1.2000,         0,    3.6000 ,
 -120.0000,  -64.0000,  -24.0000,         0,    8.0000,         0,  -24.0000 ,
 -274.5000, -146.4000,  -54.9000,         0,   18.3000,         0,  -54.9000 ,
 -423.0000, -225.6000,  -84.6000,         0,   28.2000,         0,  -84.6000 ,
 -543.0000, -289.6000, -108.6000,         0,   36.2000,         0, -108.6000  };

  ImageType2D::RegionType region;
  region.SetSize( size );

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

  typedef itk::ImageRegionIterator<ImageType2D>  InputIterator;

  InputIterator inIter( imgPtr, region );

  int j = 0;
  while( !inIter.IsAtEnd() )
    {
    inIter.Set(mydata[j]);
    ++inIter;
    ++j;
    }

  
}



void set3DDerivativeData(ImageType3D::Pointer imgPtr)
{
  SizeType3D size = {{41,41,41}};

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

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


  /* Set origin and spacing of physical coordinates */

  /* Initialize the image contents */
  // Note [x,y,z] = [x,y,z] - 20 so that image ranges from -20 + 20;
  // f(x) = 0.1x^4 - 0.5x^3 + 2x - 43
  // f(y) = 5y + 7
  // f(z) = -2z^2 - 6z + 10
  // df(x)/dx = 0.4x^3 - 1.5x^2 + 2
  // df(y)/dy = 5
  // df(z)/dz = -4z - 6
  double value;
  double slice1, row1, col1;
  IndexType3D index;
  for (unsigned int slice = 0; slice < size[2]; slice++) {
      index[2] = slice;
      slice1 = slice - 20.0;  // Center offset
      for (unsigned int row = 0; row < size[1]; row++) {
          index[1] = row;
          row1 = row - 20.0;  // Center
          for (unsigned int col = 0; col < size[0]; col++) {
              index[0] = col;
              col1 = col - 20.0; //Center
              value = 0.1*col1*col1*col1*col1 - 0.5* col1*col1*col1 + 2.0*col1 - 43.0;
              value += 5.0*row1 + 7.0;
              value += -2.0*slice1*slice1 - 6.0* slice1 + 10.0;
              imgPtr->SetPixel(index, value);
//              imgPtr->SetPixel(index, row);
          }
      }
  }

}