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

  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 "vtkEndian.h"
#include "vtkIdList.h"
#include "vtkIncrementalOctreePointLocator.h"
#include "vtkMath.h"
#include "vtkPoints.h"
#include "vtkTestUtilities.h"
#include "vtkUnstructuredGrid.h"
#include "vtkUnstructuredGridReader.h"
#include <vtksys/SystemTools.hxx>

#define VTK_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 VTK_BRUTE_FORCE_VERIFICATION
  diskFile = vtksys::SystemTools::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 VTK_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 = vtksys::SystemTools::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 VTK_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 VTK_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 VTK_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 VTK_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 VTK_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 VTK_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 VTK_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 VTK_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 VTK_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;
}