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

  Program:   Visualization Toolkit
  Module:    TestIncrementalOctreePointLocator.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.

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


#include "vtkMath.h"
#include "vtkPoints.h"
#include "vtkIdList.h"
#include "vtkTestUtilities.h"
#include "vtkUnstructuredGrid.h"
#include "vtkUnstructuredGridReader.h"
#include "vtkIncrementalOctreePointLocator.h"


#define  _BRUTE_FORCE_VERIFICATION_


// The following epsilon value is needed to address numerical inaccuracy /
// precision issues on some platforms: Fast-Release-g++, Win32-cygwin344,
// Win32-minGW, Win32-FreeVC++, and Win32-VS71. In particular, two of the
// references of this value below, as marked, handle Win32-cygwin344 and
// Win32-minGW. The numerical inaccuracy problem has nothing to do with
// the incremental octree point locator or the associated incremental octree
// node itself. Instead it is just due to the multiple sub-tests themselves
// (combined in this single file) in which the brute-force mode employs many
// if-statements / comparisons that involve double values.
//
// For example, vtkMath::Distance2BetweenPoints( point1, point2 ) may not be
// directly used in if-statements (even for some platforms other than the above
// mentioned 5), even though the incremental octree point locator always uses
// double variables for computation and returning values. Another example is
// that the min SQUARED distance D between point A (a given point) and point B
// (the closest point to A) may not be directly used to test the number of points
// within the SQUARED radius D relative to point A herein, though it is supposed
// to be ok (and the number is expected to be 1). The fact is that an epsilon
// needs to be added to D for such a test. Otherwise the numerical inaccuracy
// issue would just cause 0 to be returned --- no any point exists within the
// SQUARED radius D relative to A. Please note that this problem is not caused by
// sqrt() at all because the incremental octree point locator offers an accurate
// function variant FindPointsWithinSquaredRadius() to avoid the obvious numerical
// inaccuracy related to sqrt().
//
// Given the numerical inaccuracy issues on some platforms, the rapid verification
// mode might not be used. Fortunately, this test is fast enough.
#define  INC_OCT_PNT_LOC_TESTS_ZERO  0.00000000000001


// ---------------------------------------------------------------------------
// Meta information of the test data
//
// number of grid points           = 2288
// number of unique points         = 2200 (in terms of zero tolerance)
// number of duplicate occurrences = 88
// bounding box: [ -2.839926, 2.862497 ]
//               [ -2.856848, 2.856848 ]
//               [  0.000000, 1.125546 ]
//
// min squared distance = 1.036624e-005 (for zero-tolerance unique points)
// max squared distance = 3.391319e+001


//----------------------------------------------------------------------------
// Swap an array of (or a single) int, double, or vtkIdType values when
// loading little-endian / Win-based disk files on big-endian platforms.
// Note the two data files for this test are in little-endian binary format.
void SwapForBigEndian( unsigned char * theArray, int tupleSiz, int numTuple )
{
  int      i, j;
  unsigned char * tmpChar0 = theArray;
  unsigned char * tmpChar1 = ( unsigned char * ) malloc( tupleSiz );

  for ( j = 0; j < numTuple; j ++, tmpChar0 += tupleSiz )
  {
    for ( i = 0; i < tupleSiz; i ++ )
    {
      tmpChar1[i] = tmpChar0[ tupleSiz - 1 - i ];
    }

    for ( i = 0; i < tupleSiz; i ++ )
    {
      tmpChar0[i] = tmpChar1[i];
    }
  }

  tmpChar0 = nullptr;
  free( tmpChar1 );  tmpChar1 = nullptr;
}


//----------------------------------------------------------------------------
int TestIncrementalOctreePointLocator( int argc, char * argv[] )
{
  // specifications
  int         nClzNpts    = 4;
  int         octreRes[3] = { 64, 128, 256 };
  double      tolerans[2] = { 0.0, 0.01 };


  // variables
  int         r, t, m, i, j, k;
  int         retValue = 0;
  int         numbPnts = 0;
  int         inserted = 0;
  int         numInsrt = 0;
  int         nLocPnts = 0;
  int         nUniques = 0;
  int         nDuplics = 0;
  int         arrayIdx = 0;
  int         numExPts = 0;
  char *      fileName = nullptr;
  FILE *      diskFile = nullptr;
  FILE *      pntsFile = nullptr;
  double      pntCoord[3] = { 0.0, 0.0, 0.0 };
  double      tmpDist2 = 0.0;
  double      fTempRad = 0.0;
  double *    pDataPts = nullptr;
  double *    xtentPts = nullptr;
  double *    pLocPnts = nullptr;
  double *    minDist2 = nullptr;
  double *    maxDist2 = nullptr;
  double *    clzNdst2 = nullptr;
  vtkPoints * dataPnts = nullptr;
  vtkPoints * insrtPts = nullptr;
  vtkIdType   pointIdx;
  vtkIdType * truthIds = nullptr;
  vtkIdType * resltIds = nullptr;
  vtkIdList * ptIdList = nullptr;
  vtkIdList * idxLists[3] = { nullptr, nullptr, nullptr };
  vtkUnstructuredGrid *              unstruct = nullptr;
  vtkUnstructuredGridReader *        ugReader = nullptr;
  vtkIncrementalOctreePointLocator * octLocat = nullptr;
  size_t      n;


  // load an unstructured grid dataset
  fileName = vtkTestUtilities::ExpandDataFileName
                               ( argc, argv, "Data/post.vtk" );
  ugReader = vtkUnstructuredGridReader::New();
  ugReader->SetFileName( fileName );
  delete []  fileName;
  ugReader->Update();
  unstruct = ugReader->GetOutput();
  dataPnts = unstruct->GetPoints();
  numbPnts = dataPnts->GetNumberOfPoints();


  // open a file for reading or writing the ground truth data
  fileName = vtkTestUtilities::ExpandDataFileName
             ( argc, argv, "Data/IncOctPntLocResult.dat" );
  #ifndef _BRUTE_FORCE_VERIFICATION_
  diskFile = fopen( fileName, "rb" );
  truthIds = ( vtkIdType * )
             realloc(  truthIds,  sizeof( vtkIdType ) * numbPnts  );
  #endif
  delete []  fileName;  fileName = nullptr;


  // obtain the 3D point coordinates
  pDataPts = ( double * )
             realloc(  pDataPts,  sizeof( double ) * 3 * numbPnts  );
  for ( i = 0; i < numbPnts; i ++ )
  {
    dataPnts->GetPoint(  i,  pDataPts + ( i << 1 ) + i  );
  }


  // allocate memory for exactly duplicate points and inherit some points
  nUniques = 4;
  nDuplics = 300;
  arrayIdx = numbPnts * 3;
  numExPts = numbPnts + nUniques * nDuplics;
  xtentPts = ( double * )
             realloc(  xtentPts, sizeof( double ) * 3 * numExPts  );
  memcpy(  xtentPts,  pDataPts,  sizeof( double ) * 3 * numbPnts  );

  // add an additional number of exactly duplicate points
  for ( j = 0; j < nUniques; j ++ )
  {
    i = (  numbPnts / ( nUniques + 2 )  )  *  ( j + 1 );
    i = ( i << 1 ) + i;
    pntCoord[0] = pDataPts[ i     ];
    pntCoord[1] = pDataPts[ i + 1 ];
    pntCoord[2] = pDataPts[ i + 2 ];
    for ( i = 0; i < nDuplics; i ++, arrayIdx += 3 )
    {
      xtentPts[ arrayIdx     ] = pntCoord[0];
      xtentPts[ arrayIdx + 1 ] = pntCoord[1];
      xtentPts[ arrayIdx + 2 ] = pntCoord[2];
    }
  }


  // memory allocation
  ptIdList = vtkIdList::New();
  ptIdList->Allocate( numbPnts, numbPnts >> 1 );
  insrtPts = vtkPoints::New();
  insrtPts->Allocate( numbPnts, numbPnts >> 1 );
  octLocat = vtkIncrementalOctreePointLocator::New();


  // =========================================================================
  // ============================ Point Insertion ============================
  // =========================================================================
  for (  r = 0;  ( r < 3 ) && ( retValue == 0 );  r ++  ) // three resolutions
  {

    // -----------------------------------------------------------------------
    // --------------------- check-based point insertion ---------------------
    // -----------------------------------------------------------------------
    for (  t = 0;  ( t < 2 ) && ( retValue == 0 );  t ++  ) // two  tolerances
    for (  m = 0;  ( m < 3 ) && ( retValue == 0 );  m ++  ) // three functions
    {

      ptIdList->Reset();  // indices of the inserted points: based on dataPnts
      insrtPts->Reset();
      octLocat->FreeSearchStructure();
      octLocat->SetMaxPointsPerLeaf( octreRes[r] );
      octLocat->SetTolerance( tolerans[t] );
      octLocat->InitPointInsertion
                ( insrtPts, dataPnts->GetBounds(), numbPnts );

      // ---------------------------------------------------------------------
      if ( m == 0 )   // vtkIncrementalOctreePointLocator::InsertUniquePoint()
      {
        for ( i = 0; i < numbPnts; i ++ )
        {
          inserted = octLocat->InsertUniquePoint
                               (  pDataPts + ( i << 1 ) + i,  pointIdx  );
          if ( inserted ) ptIdList->InsertNextId( i );
        }
      }
      else
      if ( m == 1 )   // vtkIncrementalOctreePointLocator::InsertNextPoint()
      {
        for ( i = 0; i < numbPnts; i ++ )
        {
          inserted = octLocat->IsInsertedPoint(  pDataPts + ( i << 1 ) + i  );
          if ( inserted == -1 )
          {
            octLocat->InsertNextPoint(  pDataPts + ( i << 1 ) + i  );
            ptIdList->InsertNextId( i );
          }
        }
      }
      else            // vtkIncrementalOctreePointLocator::InsertPoint()
      {
        numInsrt = 0;
        for ( i = 0; i < numbPnts; i ++ )
        {
          inserted = octLocat->IsInsertedPoint(  pDataPts + ( i << 1 ) + i  );
          if ( inserted == -1 )
          {
            octLocat->InsertPoint(  numInsrt ++,  pDataPts + ( i << 1 ) + i  );
            ptIdList->InsertNextId( i );
          }
        }
      }
      // -------------------------------------------------------------------//

      #ifdef _BRUTE_FORCE_VERIFICATION_
      // ---------------------------------------------------------------------
      // verify the point insertion result in brute force mode

      double     tempPnt0[3];
      double     tempPnt1[3];
      double     tolerns2[2] = { tolerans[0] * tolerans[0],
                                 tolerans[1] * tolerans[1]
                               };
      numInsrt = ptIdList->GetNumberOfIds();

      // check if the squared distance between any two inserted points is
      // less than the threshold
      for (  j = 0;  ( j < numInsrt - 1 ) && ( retValue == 0 );  j ++  )
      {
        insrtPts->GetPoint( j, tempPnt0 );
        for (  i = j + 1;  ( i < numInsrt ) && ( retValue == 0 );  i ++  )
        {
          insrtPts->GetPoint( i, tempPnt1 );
          tmpDist2 = vtkMath::Distance2BetweenPoints( tempPnt0, tempPnt1 );
          if ( tmpDist2 <= tolerns2[t]) retValue = 1;
        }
      }

      // check if there is any point whose distance to ALL inserted points
      // is greater than the threshold
      for (  j = 0;  ( j < numbPnts ) && ( retValue == 0 );  j ++  )
      {
        if (  ptIdList->IsId( j )  !=  -1  )  continue;    // already inserted
        int   bGreater = 1;
        for ( i = 0; i < numInsrt; i ++ )
        {
          insrtPts->GetPoint( i, tempPnt1 );
          tmpDist2 = vtkMath::Distance2BetweenPoints
                              (  pDataPts + ( j << 1 ) + j,  tempPnt1  );
          if ( tmpDist2 <= tolerns2[t] )  bGreater = 0; // No 'break' here !!!
        }
        retValue = bGreater;
      }

      // -------------------------------------------------------------------//
      #else
      // ---------------------------------------------------------------------
      // rapid point index-based verification
      n = fread(  &numInsrt,  sizeof( int ),  1,  diskFile  );
      if (n != 1)
      {
        cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
        return 1;
      }
      #ifdef VTK_WORDS_BIGENDIAN
      SwapForBigEndian
        (  ( unsigned char * ) ( &numInsrt ),  sizeof( int ),  1  );
      #endif

      if (  numInsrt  ==  ptIdList->GetNumberOfIds()  )
      {
        int  samePtId = 1;
        n = fread(  truthIds,  sizeof( vtkIdType ),  numInsrt,  diskFile  );
        if (n != static_cast<size_t>(numInsrt))
        {
          cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
          return 1;
        }
        #ifdef VTK_WORDS_BIGENDIAN
        SwapForBigEndian
          (  ( unsigned char * ) truthIds,  sizeof( vtkIdType ),  numInsrt  );
        #endif

        for (  i = 0;  ( i < numInsrt ) && ( samePtId == 1 );  i ++  )
        {
          samePtId = (  truthIds[i] == ptIdList->GetId( i )  )  ?  1  :  0;
        }
        retValue = 1 - samePtId;
      }
      else  retValue = 1;
      // -------------------------------------------------------------------//
      #endif


    } // end of two tolerances and three functions


    // -----------------------------------------------------------------------
    // ------------------ direct check-free point insertion ------------------
    // -----------------------------------------------------------------------
    insrtPts->Reset();
    octLocat->FreeSearchStructure();
    octLocat->SetMaxPointsPerLeaf( octreRes[r] );
    octLocat->InitPointInsertion
              ( insrtPts, dataPnts->GetBounds(), numbPnts );
    for ( i = 0; i < numbPnts; i ++ )
    {
      octLocat->InsertPointWithoutChecking
                ( pDataPts + ( i << 1 ) + i, pointIdx, 1 );
    }
    retValue = ( insrtPts->GetNumberOfPoints() == numbPnts ) ? 0 : 1;

  } // end of three resolutions


  // =========================================================================
  // direct check-free insertion of  a huge number of EXACTLY DUPLICATE points
  //           (number > the maximum number of points per leaf node)
  // =========================================================================
  if ( retValue == 0 )
  {
    // perform direct / check-free point insertion
    for (  r = 0;  ( r < 3 ) && ( retValue == 0 );  r ++ ) // three resolutions
    {
      insrtPts->Reset();
      octLocat->FreeSearchStructure();
      octLocat->SetMaxPointsPerLeaf( octreRes[r] );
      octLocat->InitPointInsertion
                ( insrtPts, dataPnts->GetBounds(), numExPts );
      for ( i = 0; i < numExPts; i ++ )
      {
        octLocat->InsertPointWithoutChecking
                  (  xtentPts + ( i << 1 ) + i,  pointIdx,  1  );
      }

      retValue = ( insrtPts->GetNumberOfPoints() == numExPts ) ? 0 : 1;
    }
  }
  if ( xtentPts )
  {
    free( xtentPts );
    xtentPts = nullptr;
  }
  // =======================================================================//
  // =======================================================================//


  // reclaim this vtkPoints as it will be never used again
  if ( insrtPts )
  {
    insrtPts->Delete();
    insrtPts = nullptr;
  }


  // =========================================================================
  // ============================ Point  Location ============================
  // =========================================================================


  // load points and radius data from a disk file for point location tasks
  fileName = vtkTestUtilities::ExpandDataFileName
                               ( argc, argv, "Data/IncOctPntLocData.dat" );
  pntsFile = fopen( fileName, "rb" );
  delete []  fileName;  fileName = nullptr;
  n = fread(  &nLocPnts,  sizeof( int ),  1,  pntsFile  );
  if (n != 1)
  {
    cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
    fclose(pntsFile);
    free(pDataPts);
    return 1;
  }
  #ifdef VTK_WORDS_BIGENDIAN
  SwapForBigEndian
    (  ( unsigned char * ) ( &nLocPnts ),  sizeof( int ),  1  );
  #endif
  pLocPnts = ( double * ) malloc( sizeof( double ) * nLocPnts * 3 );
  minDist2 = ( double * ) malloc( sizeof( double ) * nLocPnts     );
  maxDist2 = ( double * ) malloc( sizeof( double ) * nLocPnts     );
  n = fread( pLocPnts, sizeof( double ), nLocPnts * 3, pntsFile );
  if (n != static_cast<size_t>(nLocPnts * 3))
  {
    cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
    fclose(pntsFile);
    free(pLocPnts);
    free(minDist2);
    free(maxDist2);
    free(pDataPts);
    return 1;
  }
  //fread( minDist2, sizeof( double ), nLocPnts,     pntsFile );
  //fread( maxDist2, sizeof( double ), nLocPnts,     pntsFile );
  #ifdef VTK_WORDS_BIGENDIAN
  SwapForBigEndian
    (  ( unsigned char * ) pLocPnts,  sizeof( double ),  nLocPnts * 3  );
  //SwapForBigEndian
  //  (  ( unsigned char * ) minDist2,  sizeof( double ),  nLocPnts      );
  //SwapForBigEndian
  //  (  ( unsigned char * ) maxDist2,  sizeof( double ),  nLocPnts      );
  #endif
  fclose( pntsFile );                pntsFile = nullptr;


  // destroy the context of point insertion while attaching the dataset
  octLocat->FreeSearchStructure();
  octLocat->SetDataSet( unstruct );


  // memory allocation
  clzNdst2 = ( double * ) realloc(  clzNdst2,  sizeof( double ) * nClzNpts  );
  for ( i = 0; i < 3; i ++ )
  {
    idxLists[i] = vtkIdList::New();
    idxLists[i]->Allocate( nClzNpts, nClzNpts );
  }


  // the main component
  for (  r = 0;  ( r < 3 ) && ( retValue == 0 );  r ++  ) // three resolutions
  {

    // establish a new octree with the specified resolution
    octLocat->Modified();
    octLocat->SetMaxPointsPerLeaf( octreRes[r] );
    octLocat->BuildLocator();

    // memory allocation
    resltIds = ( vtkIdType * )
               realloc( resltIds, sizeof( vtkIdType ) * nLocPnts );
    #ifndef _BRUTE_FORCE_VERIFICATION_
    truthIds = ( vtkIdType * )
               realloc( truthIds, sizeof( vtkIdType ) * nLocPnts );
    #endif


    // -----------------------------------------------------------------------
    // -------------------- location of the closest point --------------------
    // -----------------------------------------------------------------------

    // find the closest point
    for ( i = 0; i < nLocPnts; i ++ )
    {
      resltIds[i] = octLocat->FindClosestPoint(  pLocPnts + ( i << 1 ) + i,
                                                 minDist2 + i  );
    }

    #ifdef _BRUTE_FORCE_VERIFICATION_
    // -----------------------------------------------------------------------
    // verify the result in brute force mode
    for (  j = 0;  ( j < nLocPnts ) && ( retValue == 0 );  j ++  )
    for (  i = 0;  ( i < numbPnts ) && ( retValue == 0 );  i ++  )
    {
      if ( i == resltIds[j] ) continue;     // just the selected closest point
      tmpDist2 = vtkMath::Distance2BetweenPoints
                 (  pLocPnts + ( j << 1 ) + j,  pDataPts + ( i << 1 ) + i  );
      if ( tmpDist2 + INC_OCT_PNT_LOC_TESTS_ZERO < minDist2[j] ) retValue = 1;
    }
    // ---------------------------------------------------------------------//
    #else
    // -----------------------------------------------------------------------
    // rapid point index-based verification
    n = fread( &nLocPnts,  sizeof( int       ),  1,         diskFile  );
    if (n != 1)
    {
      cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
      return 1;
    }
    #ifdef VTK_WORDS_BIGENDIAN
    SwapForBigEndian
      (  ( unsigned char * ) ( &nLocPnts ),  sizeof( int ),  1  );
    #endif

    n = fread(  truthIds,  sizeof( vtkIdType ),  nLocPnts,  diskFile  );
    if (n != static_cast<size_t>(nLocPnts))
    {
      cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
      return 1;
    }
    #ifdef VTK_WORDS_BIGENDIAN
    SwapForBigEndian
      (  ( unsigned char * ) truthIds,  sizeof( vtkIdType ),  nLocPnts  );
    #endif

    for (  i = 0;  ( i < nLocPnts ) && ( retValue == 0 );  i ++  )
    {
      retValue = ( resltIds[i] == truthIds[i] ) ? 0 : 1;
    }
    // ---------------------------------------------------------------------//
    #endif

    if ( retValue == 1 ) continue;


    // -----------------------------------------------------------------------
    // ------------------ location of the closest N points -------------------
    // -----------------------------------------------------------------------

    // memory allocation
    ptIdList->SetNumberOfIds( nClzNpts * 10 );  // to claim part of the memory
    resltIds = ( vtkIdType * )
               realloc( resltIds, sizeof( vtkIdType ) * nClzNpts * nLocPnts );
    #ifndef _BRUTE_FORCE_VERIFICATION_
    truthIds = ( vtkIdType * )
               realloc( truthIds, sizeof( vtkIdType ) * nClzNpts * nLocPnts );
    n = fread( &numInsrt,  sizeof( int       ),  1,         diskFile  );
    if (n != 1)
    {
      cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
      return 1;
    }
    #ifdef VTK_WORDS_BIGENDIAN
    SwapForBigEndian
      (  ( unsigned char * ) ( &numInsrt ),  sizeof( int ),  1  );
    #endif

    n = fread(  truthIds,  sizeof( vtkIdType ),  numInsrt,  diskFile  );
    if (n != static_cast<size_t>(numInsrt))
    {
      cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
      return 1;
    }
    #ifdef VTK_WORDS_BIGENDIAN
    SwapForBigEndian
      (  ( unsigned char * ) truthIds,  sizeof( vtkIdType ),  numInsrt  );
    #endif
    #endif

    // -----------------------------------------------------------------------
    // -----------------------------------------------------------------------
    // find the closest N points (with embedded brute-force verification)
    for ( i = 0; i < nLocPnts; i ++ )
    {
      ptIdList->Reset();
      octLocat->FindClosestNPoints
                (  nClzNpts,  pLocPnts + ( i << 1 ) + i,  ptIdList  );

      // ---------------------------------------------------------------------
      // verify the result in brute force mode
      #ifdef _BRUTE_FORCE_VERIFICATION_

      // check the order of the closest points
      for ( j = 0; j < nClzNpts; j ++ )
      {
        pointIdx    = ptIdList->GetId( j );
        clzNdst2[j] = vtkMath::Distance2BetweenPoints
                               (  pLocPnts + ( i        << 1 ) + i,
                                  pDataPts + ( pointIdx << 1 ) + pointIdx  );
      }
      for (  j = 0;  ( j < nClzNpts - 1 ) && ( retValue == 0 );  j ++  )
      {
        retValue = ( clzNdst2[j + 1] >= clzNdst2[j] ) ? 0 : 1;
      }

      // write some data for locating the closest point
      // within a radius and the points within a radius
      minDist2[i] = clzNdst2[0];
      maxDist2[i] = clzNdst2[ nClzNpts - 1 ];

      // check if there are any ignored but closer points
      for (  j = 0;  ( j < numbPnts ) && ( retValue == 0 );  j ++  )
      {
        tmpDist2 = vtkMath::Distance2BetweenPoints
                   (  pLocPnts + ( i << 1 ) + i,  pDataPts + ( j << 1 ) + j  );
         if (    tmpDist2 + INC_OCT_PNT_LOC_TESTS_ZERO // for cygwin and minGW
                 < clzNdst2[ nClzNpts - 1 ]   // Not "<=" here as there
              && ptIdList->IsId( j ) == -1    // may be other points that were
            )    retValue = 1;        // rejected simply due to the limit of N
      }

      if ( retValue == 1 ) break;

      #endif
      // -------------------------------------------------------------------//

      // transfer the point indices for subsequent file writing purposes
      memcpy(  resltIds + i * nClzNpts,  ptIdList->GetPointer( 0 ),
               sizeof( vtkIdType ) * nClzNpts  );
    }
    // ---------------------------------------------------------------------//
    // ---------------------------------------------------------------------//

    // -----------------------------------------------------------------------
    // rapid point index-based verification
    #ifndef _BRUTE_FORCE_VERIFICATION_
    numInsrt = nClzNpts * nLocPnts;     // data has been read at the beginning
    for (  i = 0;  ( i < numInsrt ) && ( retValue == 0 );  i ++  )
    {
      retValue = ( resltIds[i] == truthIds[i] ) ? 0 : 1;
    }
    #endif
    // ---------------------------------------------------------------------//


    if ( retValue == 1 ) continue;


    // -----------------------------------------------------------------------
    // ------------ location of the closest point within a radius ------------
    // -----------------------------------------------------------------------

    // memory allocation
    resltIds = ( vtkIdType * )
               realloc(  resltIds,  sizeof( vtkIdType ) * nLocPnts * 3  );

    // find the closest point within three radii
    for ( i = 0; i < nLocPnts; i ++ )
    {
      j = ( i << 1 ) + i;
      pointIdx = octLocat->FindClosestPoint( pLocPnts + j, minDist2 + i );

      // there are some points falling on the in-octree points
      fTempRad = ( minDist2[i] <= INC_OCT_PNT_LOC_TESTS_ZERO )
                 ? 0.000001 : minDist2[i];
      fTempRad = sqrt( fTempRad ); // causes inaccuracy if minDist2 is nonzero

      resltIds[ j     ] = octLocat->FindClosestPointWithinRadius
                                    ( fTempRad * 0.5, pLocPnts + j, tmpDist2 );
      if ( minDist2[i] <= INC_OCT_PNT_LOC_TESTS_ZERO )
      resltIds[ j + 1 ] = octLocat->FindClosestPointWithinRadius
                                    ( fTempRad * 1.0, pLocPnts + j, tmpDist2 );
      else  // for non-zero cases, use the original squared radius for accuracy
      resltIds[ j + 1 ] = octLocat->FindClosestPointWithinSquaredRadius
                                    ( minDist2[i],    pLocPnts + j, tmpDist2 );
      resltIds[ j + 2 ] = octLocat->FindClosestPointWithinRadius
                                    ( fTempRad * 1.5, pLocPnts + j, tmpDist2 );

      // -----------------------------------------------------------------------
      // verify the result in brute force mode
      #ifdef _BRUTE_FORCE_VERIFICATION_
      if (    (    ( minDist2[i] >  INC_OCT_PNT_LOC_TESTS_ZERO )
                && ( resltIds[j] != -1 )
              )
           || (    ( minDist2[i] <= INC_OCT_PNT_LOC_TESTS_ZERO )
                && ( resltIds[j] != pointIdx )
              )
           || ( resltIds[ j + 1 ] != pointIdx )
           || ( resltIds[ j + 2 ] != pointIdx )
         ) {  retValue = 1;  break;  }
      // ---------------------------------------------------------------------//
      #endif
    }

    // -----------------------------------------------------------------------
    // rapid point index-based verification
    #ifndef _BRUTE_FORCE_VERIFICATION_
    n = fread( &numInsrt,  sizeof( int ),  1,  diskFile  );
    if (n != 1)
    {
      cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
      return 1;
    }
    #ifdef VTK_WORDS_BIGENDIAN
    SwapForBigEndian
      (  ( unsigned char * ) ( &numInsrt ),  sizeof( int ),  1  );
    #endif

    truthIds = ( vtkIdType * )
               realloc(  truthIds,  sizeof( vtkIdType ) * numInsrt  );
    n = fread(  truthIds,  sizeof( vtkIdType ),  numInsrt,  diskFile  );
    if (n != static_cast<size_t>(numInsrt))
    {
      cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
      return 1;
    }
    #ifdef VTK_WORDS_BIGENDIAN
    SwapForBigEndian
      (  ( unsigned char * ) truthIds,  sizeof( vtkIdType ),  numInsrt  );
    #endif

    for (  i = 0;  ( i < numInsrt ) && ( retValue == 0 );  i ++  )
    {
      retValue = ( resltIds[i] == truthIds[i] ) ? 0 : 1;
    }
    #endif
    // ---------------------------------------------------------------------//


    if ( retValue == 1 ) continue;


    // -----------------------------------------------------------------------
    // --------------- location of the points within a radius ----------------
    // -----------------------------------------------------------------------

    ptIdList->Reset();
    for ( i = 0; i < nLocPnts; i ++ )
    {
      // set the coordinate index of the point under check
      j = ( i << 1 ) + i;
      minDist2[i] += INC_OCT_PNT_LOC_TESTS_ZERO; // for cygwin and minGW

      // obtain the points within three radii (use squared radii for test
      // as sqrt can incur inaccuracy that complicates the test)
      idxLists[0]->Reset();
      idxLists[1]->Reset();
      idxLists[2]->Reset();
      if ( minDist2[i] <= INC_OCT_PNT_LOC_TESTS_ZERO )
      {
        // each ( maxDist2[i] * 0.3 ) has been guaranteed to be >
        // INC_OCT_PNT_LOC_TESTS_ZERO
        octLocat->FindPointsWithinSquaredRadius
                  ( maxDist2[i] * 0.3, pLocPnts + j, idxLists[0] );
        octLocat->FindPointsWithinSquaredRadius
                  ( maxDist2[i] * 0.6, pLocPnts + j, idxLists[1] );
        octLocat->FindPointsWithinSquaredRadius
                  ( maxDist2[i],       pLocPnts + j, idxLists[2] );
      }
      else
      {
        octLocat->FindPointsWithinSquaredRadius
                  ( minDist2[i] * 0.5, pLocPnts + j, idxLists[0] );
        octLocat->FindPointsWithinSquaredRadius
                  ( minDist2[i],       pLocPnts + j, idxLists[1] );
        octLocat->FindPointsWithinSquaredRadius
                  ( maxDist2[i],       pLocPnts + j, idxLists[2] );
        if ( idxLists[0]->GetNumberOfIds() == 0 )
             idxLists[0]->InsertNextId( -1 );    // to handle an empty id list
      }

      // copy the point indices to a vtkIdList for comparison and file writing
      for ( m = 0; m < 3; m ++ )
      {
        numInsrt = idxLists[m]->GetNumberOfIds();
        for ( k = 0; k < numInsrt; k ++ )
        {
          ptIdList->InsertNextId(  idxLists[m]->GetId( k )  );
        }
      }

      // ---------------------------------------------------------------------
      // verify the result in brute force mode
      #ifdef _BRUTE_FORCE_VERIFICATION_

      // check if the monotonical property holds among the 3 point-index lists
      if ( minDist2[i] <= INC_OCT_PNT_LOC_TESTS_ZERO )
      {
        pointIdx = octLocat->FindClosestPoint( pLocPnts + j );
        if ( idxLists[0]->IsId( pointIdx ) == -1 ||
             idxLists[1]->IsId( pointIdx ) == -1 ||
             idxLists[2]->IsId( pointIdx ) == -1 ||
             idxLists[1]->GetNumberOfIds() < idxLists[0]->GetNumberOfIds() ||
             idxLists[2]->GetNumberOfIds() < idxLists[0]->GetNumberOfIds() ||
             idxLists[2]->GetNumberOfIds() < idxLists[1]->GetNumberOfIds()
           ) retValue = 1;
      }
      else
      {
        if ( idxLists[0]->GetNumberOfIds() !=  1 ||
             idxLists[0]->GetId( 0 )       != -1 ||
             idxLists[1]->GetNumberOfIds() <   1 ||
             idxLists[2]->GetNumberOfIds() <  idxLists[1]->GetNumberOfIds()
           ) retValue = 1;
      }

      // check the points within each of the three radii
      for (  m = 0;  ( m < 3 ) && ( retValue == 0 );  m ++  )
      {
        // get the squared radius actually used
        if ( minDist2[i] <= INC_OCT_PNT_LOC_TESTS_ZERO )
        {
          if ( m == 0 ) fTempRad = maxDist2[i] * 0.3;
          else
          if ( m == 1 ) fTempRad = maxDist2[i] * 0.6;
          else          fTempRad = maxDist2[i];
        }
        else
        {
          if ( m == 0 ) fTempRad = minDist2[i] * 0.5;
          else
          if ( m == 1 ) fTempRad = minDist2[i];
          else          fTempRad = maxDist2[i];
        }

        // check if there is any false insertion
        numInsrt = idxLists[m]->GetNumberOfIds();
        for (  k = 0;  ( k < numInsrt ) && ( retValue == 0 );  k ++  )
        {
          if (  m == 0 && idxLists[0]->GetId( 0 )  ==  -1  ) break;

          pointIdx = idxLists[m]->GetId( k );
          pointIdx = ( pointIdx << 1 ) + pointIdx;
          tmpDist2 = vtkMath::Distance2BetweenPoints
                              ( pLocPnts + j, pDataPts + pointIdx );

          if ( tmpDist2 > fTempRad + INC_OCT_PNT_LOC_TESTS_ZERO ) retValue = 1;
        }

        // check if there is any missed insertion
        numInsrt = 0;
        for ( k = 0; k < numbPnts; k ++ )
        {
          tmpDist2 = vtkMath::Distance2BetweenPoints
                              (  pLocPnts + j,  pDataPts + ( k << 1 ) + k  );
          if ( tmpDist2 + INC_OCT_PNT_LOC_TESTS_ZERO <= fTempRad ) numInsrt ++;
        }

        // get the actual size of the vtkIdList for comparison
        int  listSize = ( m == 0 && idxLists[0]->GetId( 0 ) == -1 )
                        ? 0       // for an actually nullptr vtkIdList
                        : idxLists[m]->GetNumberOfIds();
        if ( numInsrt > listSize ) retValue = 1;
      }
      #endif
      // -------------------------------------------------------------------//
    }

    // -----------------------------------------------------------------------
    // rapid point index-based verification
    #ifndef _BRUTE_FORCE_VERIFICATION_
    n = fread( &numInsrt,  sizeof( int ),  1,  diskFile  );
    if (n != 1)
    {
      cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
      return 1;
    }
    #ifdef VTK_WORDS_BIGENDIAN
    SwapForBigEndian
      (  ( unsigned char * ) ( &numInsrt ),  sizeof( int ),  1  );
    #endif

    truthIds = ( vtkIdType * )
               realloc(  truthIds,  sizeof( vtkIdType ) * numInsrt  );
    n = fread(  truthIds,  sizeof( vtkIdType ),  numInsrt,  diskFile  );
    if (n != static_cast<size_t>(numInsrt))
    {
      cerr << "IO error " << __FILE__ << ":" << __LINE__ << "\n";
      return 1;
    }
    #ifdef VTK_WORDS_BIGENDIAN
    SwapForBigEndian
      (  ( unsigned char * ) truthIds,  sizeof( vtkIdType ),  numInsrt  );
    #endif

    vtkIdType * tmpPtIds = ptIdList->GetPointer( 0 );
    for (  i = 0;  ( i < numInsrt ) && ( retValue == 0 );  i ++  )
    {
      retValue = ( tmpPtIds[i] == truthIds[i] ) ? 0 : 1;
    }
    tmpPtIds = nullptr;
    #endif
    // ---------------------------------------------------------------------//

  }


  // memory clearance
  dataPnts = nullptr;
  unstruct = nullptr;
  if ( ptIdList )
  {
    ptIdList->Delete();
    ptIdList = nullptr;
  }
  if ( octLocat )
  {
    octLocat->Delete();
    octLocat = nullptr;
  }
  if ( ugReader )
  {
    ugReader->Delete();
    ugReader = nullptr;
  }

  free( truthIds );
  truthIds = nullptr;

  free( resltIds );
  resltIds = nullptr;

  free( pDataPts );
  pDataPts = nullptr;

  free( pLocPnts );
  pLocPnts = nullptr;

  free( minDist2 );
  minDist2 = nullptr;

  free( maxDist2 );
  maxDist2 = nullptr;

  free( clzNdst2 );
  clzNdst2 = nullptr;

  if ( diskFile )
  {
    fclose( diskFile );
    diskFile = nullptr;
  }

  for ( i = 0; i < 3; i ++ )
  {
    if ( idxLists[i] )
    {
      idxLists[i]->Delete();
      idxLists[i] = nullptr;
    }
  }


  return retValue;
}