text: /*=========================================================================
  
    Program:   Visualization Toolkit
    Module:    TestImageDataInterpolation.cxx
  
    Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
    All rights reserved.
    See Copyright.txt or http://www.kitware.com/Copyright.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 notice for more information.
  
  =========================================================================*/
  // .NAME TestImageDataInterpolation.cxx -- Test interpolation from image data
  //
  // .SECTION Description
  //  This test applies a function, F(x,y,z), to the image data and tests the
  //  interpolation.
  
  #include "vtkDebugLeaks.h"
  #include "vtkImageData.h"
  #include "vtkSmartPointer.h"
  #include "vtkPoints.h"
  #include "vtkDoubleArray.h"
  #include "vtkIdList.h"
  #include "vtkPointData.h"
  #include "vtkCell.h"
  
  #include <cassert>
  #include <cmath>
  #include <limits>
  
  // Description:
  // Performs safe division a/b which also checks for underflow & overflow
  double SafeDiv( const double a, const double b )
  {
    // Catch overflow
    if( (b < 1) && (a > ( b*std::numeric_limits<double>::max() ) ) )
      return std::numeric_limits<double>::max();
  
    // Catch underflow
    if( (a == static_cast<double>(0.0)) ||
        ( (b > 1) && (a < b*std::numeric_limits<double>::max() ) ) )
        return( static_cast<double>( 0.0 ) );
  
    return( a/b );
  }
  
  // Description:
  // Checks if two given floating numbers are equivalent.
  // This algorithm is based on Knuth, The of Computer Programming (vol II).
  bool eq( double a, double b, double TOL=1e-9 )
  {
    double adiff = std::abs( a-b );
    double d1    = SafeDiv( adiff,std::abs(a) );
    double d2    = SafeDiv( adiff,std::abs(b) );
    if( (d1 <= TOL) || (d2 <= TOL) )
      return true;
    return false;
  }
  
  
  // Description:
  // Applies the function to the point with the given coordinates.
  // The function is defined as F(x,y,z)= x + y + z
  inline double F( const double x, const double y, const double z )
  {
    return( (x+y+z) );
  }
  
  // Description:
  // Returns the grid with the specified extent, origin and spacing.
  inline vtkImageData* GetGrid( int dims[3], double origin[3], double h[3] )
  {
    vtkImageData *image = vtkImageData::New();
    assert( "pre: image data is NULL!" && (image != NULL) );
  
    image->SetDimensions( dims );
    image->SetOrigin( origin );
    image->SetSpacing( h );
  
    vtkPointData *PD = image->GetPointData();
    assert( "pre: point-data is NULL" && (PD != NULL) );
  
    vtkDoubleArray *dataArray = vtkDoubleArray::New();
    dataArray->SetName( "Fx" );
    dataArray->SetNumberOfTuples( image->GetNumberOfPoints() );
    dataArray->SetNumberOfComponents( 1 );
  
    vtkIdType idx = 0;
    for( ; idx < image->GetNumberOfPoints(); ++idx )
    {
      double pnt[3];
      image->GetPoint(idx, pnt);
      dataArray->SetComponent( idx, 0, F(pnt[0],pnt[1],pnt[2]) );
    } // END for all points
  
    PD->AddArray( dataArray );
    dataArray->Delete();
  
    return( image );
  }
  
  // Description:
  // Given the image data this method returns a list of test points for
  // interpolation at each cell center.
  inline vtkPoints* GetReceivePoints(
      vtkImageData *img, vtkIdList* donorCellList)
  {
    vtkPoints *rcvPoints = vtkPoints::New();
    rcvPoints->SetNumberOfPoints( img->GetNumberOfCells() );
    donorCellList->SetNumberOfIds( img->GetNumberOfCells() );
  
    vtkIdType cellIdx = 0;
    for( ; cellIdx < img->GetNumberOfCells(); ++cellIdx  )
    {
      // Get cell
      vtkCell *myCell = img->GetCell( cellIdx  );
      assert( "pre: myCell != NULL" && (myCell != NULL) );
  
      // Compute center
      double c[3];
      double pCenter[3];
      double *weights = new double[ myCell->GetNumberOfPoints() ];
      int subId       = myCell->GetParametricCenter( pCenter );
      myCell->EvaluateLocation( subId,pCenter,c,weights );
      delete [] weights;
  
      // Insert center to point list
      donorCellList->SetId( cellIdx, cellIdx );
      rcvPoints->SetPoint( cellIdx, c );
    } // END for all cells
  
    return( rcvPoints );
  }
  
  // Description:
  // Given the mesh data, the donor cell and the interpolation weights, this
  // method returns the interpolated value at the corresponding point location.
  inline double InterpolateValue(vtkImageData *img, vtkCell *c, double* w)
  {
    assert( "pre: image is NULL!" && (img != NULL) );
    assert( "pre: donor cell is NULL" && (c != NULL) );
    assert( "pre: weights is NULL" && (w != NULL) );
  
    vtkPointData *PD = img->GetPointData();
    assert( "pre: point data is NULL" && (PD != NULL) );
    vtkDataArray *dataArray = PD->GetArray( "Fx" );
    assert( "pre: data array is NULL" && (dataArray != NULL) );
  
    std::cout << "W:[";
    std::cout.flush();
    double val         = 0.0;
    vtkIdType numNodes = c->GetNumberOfPoints();
    for( vtkIdType nodeIdx=0; nodeIdx < numNodes; ++nodeIdx )
    {
      std::cout << w[ nodeIdx ] << " ";
      std::cout.flush();
      vtkIdType meshNodeIdx = c->GetPointId( nodeIdx );
      val += w[nodeIdx]*dataArray->GetComponent( meshNodeIdx, 0 );
    } // END for all cells nodes
    std::cout << "]\n";
    std::cout.flush();
    return( val );
  }
  
  // Description:
  // Main test routine for testing the interpolation
  inline int TestInterpolation( int dims[3], double origin[3], double h[3] )
  {
    vtkImageData *grid = GetGrid( dims, origin, h );
    assert( "pre: grid is NULL" && (grid != NULL) );
  
    std::cout << "NUMBER OF CELLS:  " << grid->GetNumberOfCells() << std::endl;
    std::cout.flush();
  
    vtkIdList *donors = vtkIdList::New();
    vtkPoints *pnts = GetReceivePoints( grid, donors );
  
    int subId = 0;
    double pcoords[3];
    double weights[8];
    double x[3];
  
    int interpErrors = 0;
  
    vtkIdType idx = 0;
    for( ; idx < pnts->GetNumberOfPoints(); ++idx )
    {
      pnts->GetPoint( idx, x );
      vtkIdType cellIdx =
       grid->FindCell(x,NULL,0,0.0,subId,pcoords,weights);
  
      if( cellIdx < 0 )
      {
        std::cerr << "point (" << x[0] << ", " << x[1] << ", " << x[2];
        std::cerr << ") is out-of-bounds!\n";
        return 1;
      }
  
      vtkCell *donorCell = grid->GetCell( cellIdx );
      assert( "pre: donor cell is NULL" && (donorCell != NULL));
  
      std::cout << "N:  [" << pcoords[0] << " ";
      std::cout << pcoords[1] << " " << pcoords[2] << "]\n";
      std::cout.flush();
  
      double f         = InterpolateValue( grid, donorCell, weights );
      double fExpected = F( x[0],x[1],x[2] );
  
      if( !eq( f, fExpected ) )
      {
        std::cout << "INTERPOLATION ERROR: ";
        std::cout << "f_expeted=" << fExpected << " ";
        std::cout << "f_interp=" << f << std::endl;
        std::cout.flush();
        ++ interpErrors;
      }
  
    } // END for all points
  
    grid->Delete();
    donors->Delete();
    pnts->Delete();
  
    return interpErrors;
  }
  
  // Description:
  // Test main function
  int TestImageDataInterpolation(int, char *[])
  {
    static int dims[4][3]={
        {3,3,1}, // XY Plane
        {3,1,3}, // XZ Plane
        {1,3,3}, // YZ Plane
        {3,3,3}  // XYZ Volume
    };
  
    static const char* TestName[4]={
        "XY-Plane",
        "XZ-Plane",
        "YZ-Plane",
        "XYZ-Grid"
    };
  
    double origin[3];
    origin[0] = origin[1] = origin[2] = 0.0;
  
    double h[3];
    h[0] = h[1] = h[2] = 0.2;
  
    int testStatus = 0;
    for( int i=0; i < 4; ++i )
    {
      std::cout << "Testing " << TestName[ i ] << "...\n";
      std::cout.flush();
      int rc = 0;
      rc     = TestInterpolation( dims[i], origin, h );
      testStatus += rc;
      std::cout << "[DONE]\n";
      std::cout.flush();
  
      if( rc != 0 )
      {
        std::cout << rc << " failures detected!\n";
        std::cout << "TEST FAILED!\n\n";
        std::cout.flush();
      }
      else
      {
        std::cout << "TEST PASSED!\n";
        std::cout << std::endl;
      }
    } // END for all tests
  
    return( testStatus );
  
  }
  