InterpIdw.cpp

Go to the documentation of this file.

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------

//----- Included files ---------------------------------------------------------

// 1. Precompiled header
#include <boost/thread.hpp>
// 2. My header
#include <xmsinterp/interpolate/InterpIdw.h>

// 3. Standard Library Headers

// 4. External Library Headers

// 5. Shared Headers
#include <xmscore/math/math.h>
#include <xmscore/stl/utility.h>
#include <xmscore/stl/vector.h>
#include <xmsgrid/geometry/GmPtSearch.h>
#include <xmsinterp/interpolate/InterpSignals.h>
#include <xmsinterp/interpolate/InterpUtil.h>
#include <xmsinterp/interpolate/detail/NodalFunc.h>
#include <xmsgrid/geometry/geoms.h>
#include <xmsinterp/thread/ThreadLoop.h>
#include <xmsinterp/thread/ThreadMgr.h>
#include <xmscore/misc/XmError.h>
#include <xmscore/misc/xmstype.h>

// 6. Non-shared Headers

//----- Forward declarations ---------------------------------------------------

//----- External globals -------------------------------------------------------

//----- Namespace declaration --------------------------------------------------

namespace xms
{
//----- Constants / Enumerations -----------------------------------------------

//----- Classes / Structs ------------------------------------------------------
class InterpIdwImpl : public InterpIdw
{
public:
  class InterpThread : public ThreadLoop
  {
  public:
    InterpThread(InterpIdwImpl& a_, VecFlt& a_s, const VecPt3d& a_p)
    : m_interp(a_)
    , m_scalarTo(a_s)
    , m_pts(a_p)
    {
    } // InterpThread

    InterpIdwImpl& m_interp;        
    VecFlt& m_scalarTo; 
    const VecPt3d& m_pts; 

    BSHP<ThreadLoop> CreateForNewThread() override
    {
      BSHP<InterpThread> r(new InterpThread(m_interp, m_scalarTo, m_pts));
      BSHP<ThreadLoop> q = BDPC<ThreadLoop>(r);
      return q;
    }

  private:
    void Worker() override
    {
      int i = (int)CurrIdx();
      m_scalarTo[i] = m_interp.InterpToPt(m_pts[i], i);
    }
  };

  InterpIdwImpl();
  virtual ~InterpIdwImpl();

  virtual void SetPtsTris(BSHP<VecPt3d> a_pts, BSHP<VecInt> a_tris) override;
  void SetScalars(const float* a_scalar, size_t a_n) override;
  void SetScalars(BSHP<VecFlt> a_scalar) override;
  virtual void SetPts(BSHP<VecPt3d> a_pts, bool a_2d) override;
  // this method will run parallel
  virtual float InterpToPt(const Pt3d& a_pt) override;
  virtual void InterpToPts(const VecPt3d& a_pts, VecFlt& a_scalars) override;
  virtual void SetPtActivity(DynBitset& a_activity) override;
  virtual void SetTriActivity(DynBitset&) override {}
  virtual const BSHP<VecPt3d> GetPts() const override;
  virtual const BSHP<VecInt> GetTris() const override;
  virtual const BSHP<VecFlt> GetScalars() const override;
  virtual DynBitset GetPtActivity() const override;
  virtual DynBitset GetTriActivity() const override;
  virtual void SetTrunc(double a_sMax, double a_sMin) override;

  virtual bool GetTruncateInterpolatedValues() const override { return m_trunc; }
  virtual double GetTruncMin() const override { return m_truncMin; }
  virtual double GetTruncMax() const override { return m_truncMax; }
  virtual double GetPower() const override { return m_power; }
  virtual int GetSearchOptsNumNearestPts() const override { return m_nNearestPts; }
  virtual bool GetSearchOptsUseQuadrantSearch() const override { return m_quadOctSearch; }
  virtual WeightEnum GetWeightCalcMethod() const override
  {
    if (!m_modifiedShepardWeights)
      return InterpIdw::CLASSIC;
    return InterpIdw::MODIFIED;
  }
  virtual NodalFuncEnum GetNodalFunctionType() const override
  {
    if (!m_nodalFunc)
      return InterpIdw::CONSTANT;
    return static_cast<InterpIdw::NodalFuncEnum>(m_nodalFunc->GetType());
  }
  virtual int GetNodalFunctionNumNearestPts() const override
  {
    if (!m_nodalFunc)
      return 0;
    return m_nodalFunc->GetNumNearestPoints();
  }
  virtual bool GetNodalFunctionUseQuadrantSearch() const override
  {
    if (!m_nodalFunc)
      return false;
    return m_nodalFunc->GetUseQuadrantSearch();
  }


  virtual void SetSaveWeights(bool a_) override { m_saveWeights = a_; }
  virtual void InterpWeights(const Pt3d& a_pt,
                             VecInt& a_idx,
                             VecDbl& a_wt) const override;

  virtual void SetMultiThreading(bool a_) override { m_multiThread = a_; }

  virtual std::string ToString() const override;

  void SetPower(double a_power) override;
  void RecalcNodalFunc();
  void SetSearchOpts(int a_nNearestPoints, bool a_quad_oct_Search) override;
  float InterpToPt(const Pt3d& a_pt, int a_idx) const;
  void SetObserver(BSHP<Observer> a_prog) override { m_prog = a_prog; }

  void ValFromWeights(VecDbl& a_w,
                      VecInt& a_nPts,
                      float& a_val,
                      int a_ptIdx,
                      const Pt3d& a_pt) const;
  void SetWeightCalcMethod(InterpIdw::WeightEnum a_) override;
  void SetNodalFunction(NodalFuncEnum a_,
                        int a_nNearest,
                        bool a_quad_oct,
                        BSHP<Observer> a_p) override;
  void CreateNodalFunctionPtr(NodalFuncEnum a_,
                              int a_nNearest,
                              bool a_quad_oct,
                              BSHP<Observer> a_p);
  double ScalarFromNodalFunc(int a_ptIdx, const Pt3d& a_loc) const;

  bool m_2d; 
  bool m_quadOctSearch;
  bool m_modifiedShepardWeights;
  int m_nNearestPts;             
  double m_power;
  BSHP<GmPtSearch> m_ptSearch; 
  BSHP<VecPt3d> m_pts;         
  BSHP<VecInt> m_tris;         
  BSHP<VecFlt> m_scalarFrom;   
  BSHP<NodalFunc> m_nodalFunc; 
  BSHP<Observer> m_prog;       
  bool m_saveWeights;          
  VecInt2d m_ptIdx;            
  VecDbl2d m_weights;          
  bool m_trunc;                
  double m_truncMax;      
  double m_truncMin;      
  bool m_multiThread;     
  std::string m_idString; 
};

//----- Internal functions -----------------------------------------------------
namespace
{
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
static VecInt2d& iPtIdx()
{
  static VecInt2d m_;
  return m_;
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
static VecDbl2d& iWeights()
{
  static VecDbl2d m_weights;
  return m_weights;
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
static void iSaveWeight(int a_ptIdx, int a_pt, double a_val)
{
  XM_ENSURE_TRUE_VOID_NO_ASSERT(a_ptIdx > XM_NONE);
  iPtIdx()[a_ptIdx].push_back(a_pt);
  iWeights()[a_ptIdx].push_back(a_val);
} // iSaveWeight
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
static void iCallSignals(VecInt2d& a_ptIdx,
                         VecDbl2d& a_weights)
{
  for (size_t i = 0; i < a_ptIdx.size(); ++i)
  {
    inSetDataIndex((int)i);
    for (size_t j = 0; j < a_ptIdx[i].size(); ++j)
    {
      inFillWtArray(a_ptIdx[i][j] + 1, a_weights[i][j]);
    }
  }
  a_ptIdx.clear();
  a_weights.clear();
} // iCallSignals
} // unnamed namespace
//----- Class / Function definitions -------------------------------------------

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
BSHP<InterpIdw> InterpIdw::New()
{
  BSHP<InterpIdw> ptr(new InterpIdwImpl());
  return ptr;
} // InterpIdw::New
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
InterpIdw::InterpIdw()
{
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
InterpIdw::~InterpIdw()
{
}

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
InterpIdwImpl::InterpIdwImpl()
: m_2d(true)
, m_quadOctSearch(false)
, m_modifiedShepardWeights(true)
, m_nNearestPts(16)
, m_power(2)
, m_ptSearch()
, m_pts(new VecPt3d())
, m_tris(new VecInt())
, m_scalarFrom(new VecFlt())
, m_nodalFunc()
, m_prog()
, m_saveWeights(false)
, m_trunc(false)
, m_truncMax(0.0)
, m_truncMin(0.0)
, m_multiThread(true)
{
} // InterpIdwImpl::InterpIdwImpl
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
InterpIdwImpl::~InterpIdwImpl()
{
} // InterpIdwImpl::~InterpIdwImpl
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetPtsTris(BSHP<VecPt3d> a_pts, BSHP<VecInt> a_tris)
{
  SetPts(a_pts, true);
  m_tris = a_tris;
} // InterpIdwImpl::SetPtsTris
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetPts(BSHP<VecPt3d> a_pts, bool a_2d)
{
  m_2d = a_2d;
  m_pts = a_pts;
  m_ptSearch = GmPtSearch::New(m_2d);
  m_ptSearch->PtsToSearch(m_pts);
  if (m_2d && m_scalarFrom->empty())
  {
    m_scalarFrom->resize(0);
    m_scalarFrom->reserve(m_pts->size());
    VecFlt& s(*m_scalarFrom);
    VecPt3d& pts(*m_pts);
    for (Pt3d& p : pts)
      s.push_back((float)p.z);
  }
} // InterpIdwImpl::SetPts
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetScalars(const float* a_scalar, size_t a_n)
{
  m_scalarFrom->resize(0);
  m_scalarFrom->insert(m_scalarFrom->begin(), a_scalar, a_scalar + a_n);
  RecalcNodalFunc();
} // InterpIdwImpl::SetScalars
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetScalars(BSHP<VecFlt> a_scalar)
{
  m_scalarFrom = a_scalar;
  RecalcNodalFunc();
} // InterpIdwImpl::SetScalars
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetPower(double a_power)
{
  m_power = a_power;
  XM_ENSURE_TRUE(!m_nodalFunc); // Not handled: SetPower after creating nodal func
} // InterpIdwImpl::SetPower
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::RecalcNodalFunc()
{
  if (m_nodalFunc)
  {
    m_nodalFunc->ComputeNodalFuncs();
  }
} // InterpIdwImpl::RecalcNodalFunc
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetSearchOpts(int a_nNearestPoints, bool a_quad_oct_Search)
{
  m_nNearestPts = a_nNearestPoints;
  m_quadOctSearch = a_quad_oct_Search;
} // InterpIdwImpl::SetPower
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
float InterpIdwImpl::InterpToPt(const Pt3d& a_pt)
{
  return InterpToPt(a_pt, XM_NONE);
} // InterpIdwImpl::InterpToPt
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
float InterpIdwImpl::InterpToPt(const Pt3d& a_pt, int a_idx) const
{
  float rval(0);
  XM_ENSURE_FALSE(m_scalarFrom->empty(), rval);
  VecInt nPts;
  VecDbl w;
  InterpWeights(a_pt, nPts, w);
  if (1 == nPts.size())
  {
    rval = (*m_scalarFrom)[nPts[0]];
    if (m_saveWeights)
      iSaveWeight(a_idx, nPts[0], 1.0);
    return rval;
  }
  ValFromWeights(w, nPts, rval, a_idx, a_pt);
  if (m_trunc)
  {
    if (rval < m_truncMin)
      rval = (float)m_truncMin;
    if (rval > m_truncMax)
      rval = (float)m_truncMax;
  }
  return rval;
} // InterpIdwImpl::InterpToPt
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::InterpToPts(const VecPt3d& a_pts, VecFlt& a_scalars)
{
  a_scalars.assign(a_pts.size(), 0);
  if (m_saveWeights)
  {
    iPtIdx().assign(a_scalars.size(), VecInt());
    iWeights().assign(a_scalars.size(), VecDbl());
  }

  if (!m_multiThread)
  {
    for (size_t i = 0; i < a_pts.size(); ++i)
    {
      a_scalars[i] = InterpToPt(a_pts[i], (int)i);
    }
  }
  else
  {
    BSHP<ThreadMgr> threadMgr = ThreadMgr::New();
    threadMgr->SetThreadLoopClass(InterpThread(*this, a_scalars, a_pts).CreateForNewThread());
    threadMgr->SetObserver(m_prog);
    if (a_pts.size() < 10000)
      threadMgr->ExplicitlySetNumThreads(2);
    threadMgr->RunThreads((int)a_scalars.size());
  }

  if (m_saveWeights)
    iCallSignals(iPtIdx(), iWeights());
} // InterpIdwImpl::InterpToPts
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetPtActivity(DynBitset& a_activity)
{
  m_ptSearch->SetActivity(a_activity);
} // InterpIdwImpl::SetPtActivity
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const BSHP<VecPt3d> InterpIdwImpl::GetPts() const
{
  return m_pts;
} // InterpIdwImpl::GetPts
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const BSHP<VecInt> InterpIdwImpl::GetTris() const
{
  return m_tris;
} // InterpIdwImpl::GetTris
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const BSHP<VecFlt> InterpIdwImpl::GetScalars() const
{
  return m_scalarFrom;
} // InterpIdwImpl::GetScalars
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
DynBitset InterpIdwImpl::GetPtActivity() const
{
  if (!m_ptSearch)
    return DynBitset();
  return m_ptSearch->GetActivity();
} // InterpIdwImpl::GetPtActivity
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
DynBitset InterpIdwImpl::GetTriActivity() const
{
  return DynBitset();
} // InterpIdwImpl::GetTriActivity
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetTrunc(double a_sMax, double a_sMin)
{
  m_trunc = true;
  m_truncMax = a_sMax;
  m_truncMin = a_sMin;
} // InterpIdwImpl::SetTrunc
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::InterpWeights(const Pt3d& a_pt,
                                  VecInt& a_nPts,
                                  VecDbl& a_w) const
{
  VecDbl d2;
  // first see if we are at this location
  m_ptSearch->NearestPtsToPt(a_pt, 1, false, a_nPts);
  // get the distance between our location and the nearest point
  inDistanceSquared(a_pt, a_nPts, *m_pts, m_2d, d2);
  if (m_nNearestPts == 1 || d2[0] < gmXyTol())
  {
    a_w.resize(1);
    a_w[0] = 1;
    return;
  }
  // see if we are doing all points
  if (m_nNearestPts < 0)
  {
    a_nPts.assign(m_pts->size(), 0);
    for (size_t i = 0; i < m_pts->size(); ++i)
      a_nPts[i] = (int)i;
  }
  else
  {
    m_ptSearch->NearestPtsToPt(a_pt, m_nNearestPts, m_quadOctSearch, a_nPts);
  }
  // compute weights for each point
  inDistanceSquared(a_pt, a_nPts, *m_pts, m_2d, d2);
  inIdwWeights(d2, m_power, m_modifiedShepardWeights, a_w);
} // InterpIdwImpl::InterpWeights
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::string InterpIdwImpl::ToString() const
{
  std::stringstream ss;
  ss << m_2d << "=2d " << m_quadOctSearch << "=quadOctSearch " << m_modifiedShepardWeights
     << "=modifiedShepardWeights " << m_nNearestPts << "=nNearestPts " << m_power << "=power "
     << m_saveWeights << "=saveWeights " << m_multiThread << "=multiThread "
     << "\n";

  if (m_ptSearch)
    ss << m_ptSearch->ToString() << "=ptSearch "
       << "\n";

  if (m_pts)
    VecToStream(ss, *m_pts, "pts");
  if (m_tris)
    VecToStream(ss, *m_tris, "tris");
  if (m_scalarFrom)
    VecToStream(ss, *m_scalarFrom, "scalarFrom");

  if (m_nodalFunc)
    ss << m_nodalFunc->ToString() << "=nodalFunc "
       << "\n";

  for (size_t i = 0; i < m_ptIdx.size(); ++i)
    VecToStream(ss, m_ptIdx[i], "ptIdx");
  ss << "=ptIdx "
     << "\n";

  for (size_t i = 0; i < m_weights.size(); ++i)
    VecToStream(ss, m_weights[i], "weights");
  ss << "=weights "
     << "\n";

  return ss.str();
} // InterpIdwImpl::ToString
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::ValFromWeights(VecDbl& a_w,
                                   VecInt& a_nPts,
                                   float& a_val,
                                   int a_ptIdx,
                                   const Pt3d& a_pt) const
{
  a_val = 0;
  double val(0);
  for (size_t i = 0; i < a_w.size(); ++i)
  {
    int j = a_nPts[i];
    // val += a_w[i] * m_scalarFrom[j];
    val += a_w[i] * ScalarFromNodalFunc(j, a_pt);
    if (m_saveWeights)
      iSaveWeight(a_ptIdx, j, a_w[i]);
  }
  a_val = (float)val;
} // InterpIdwImpl::ValFromWeights
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetWeightCalcMethod(InterpIdw::WeightEnum a_)
{
  if (a_ == InterpIdw::CLASSIC)
    m_modifiedShepardWeights = false;
  else if (a_ == InterpIdw::MODIFIED)
    m_modifiedShepardWeights = true;
  else
    m_modifiedShepardWeights = true;
} // InterpIdwImpl::SetWeightCalcMethod
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::SetNodalFunction(NodalFuncEnum a_,
                                     int a_nNearest,
                                     bool a_quad_oct,
                                     BSHP<Observer> a_p)
{
  XM_ENSURE_TRUE_VOID_NO_ASSERT(a_ != CONSTANT);
  XM_ENSURE_TRUE_VOID_NO_ASSERT(!inAllScalarsEqual(*m_scalarFrom, m_ptSearch->GetActivity()));
  XM_ENSURE_TRUE(!m_pts->empty());

  CreateNodalFunctionPtr(a_, a_nNearest, a_quad_oct, a_p);
} // InterpIdwImpl::SetNodalFunction
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwImpl::CreateNodalFunctionPtr(NodalFuncEnum a_,
                                           int a_nNearest,
                                           bool a_quad_oct,
                                           BSHP<Observer> a_p)
{
  m_nodalFunc = NodalFunc::New((int)a_, m_2d, m_ptSearch, *m_pts, *m_scalarFrom, a_nNearest,
                               a_quad_oct, m_power, m_modifiedShepardWeights, a_p, nullptr);
} // InterpIdwImpl::CreateNodalFunctionPtr
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
double InterpIdwImpl::ScalarFromNodalFunc(int a_ptIdx, const Pt3d& a_loc) const
{
  double rval((double)(*m_scalarFrom)[a_ptIdx]);
  if (m_nodalFunc)
  {
    rval = m_nodalFunc->ScalarFromPtIdx(a_ptIdx, a_loc);
  }
  return rval;
} // InterpIdwImpl::ScalarFromNodalFunc

} // namespace xms

#ifdef CXX_TEST

#include <xmsinterp/interpolate/InterpIdw.t.h>

#include <xmscore/testing/TestTools.h>
#include <xmscore/misc/XmLog.h>

// namespace xms {
using namespace xms;
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
static void iGetPoints(VecPt3d& a_pts,
                       VecFlt& a_scalar,
                       VecPt3d& a_iPts)
{
  a_pts = {{-75.0, -16.0, 0.0}, {-60.0, 32.0, 0.0},   {-34.0, 50.0, 0.0},   {-34.0, 27.0, 0.0},
           {-8.0, 40.0, 0.0},   {16.0, 38.0, 0.0},    {-25.0, 14.0, 43.64}, {10.0, 18.0, 44.16},
           {27.0, 26.0, 0.0},   {63.0, 35.0, 0.0},    {-32.0, 0.0, 59.04},  {-7.0, 7.0, 90.2},
           {26.0, 6.0, 67.2},   {75.0, 7.0, 0.0},     {-37.0, -15.0, 9.24}, {-7.0, -13.0, 71.0},
           {2.0, -3.0, 98.4},   {31.0, -15.0, 25.56}, {60.0, -13.0, 0.0},   {-50.0, -30.0, 0.0},
           {-30.0, -28.0, 0.0}, {43.0, -22.0, 0.0},   {-32.0, -50.0, 0.0},  {27.0, -37.0, 0.0},
           {60.0, -33.0, 0.0}};
  a_scalar.resize(a_pts.size());
  for (size_t i = 0; i < a_pts.size(); ++i)
    a_scalar[i] = (float)a_pts[i].z;

  a_iPts = {{-90, 60, 0}, {90, 60, 0}, {-90, -60, 0}, {90, -60, 0}, {60, -33, 0}};
} // iGetPoints

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testCreateClass()
{
  BSHP<InterpIdw> pts = InterpIdw::New();
  TS_ASSERT(pts);
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testInterpToPts()
{
  BSHP<InterpIdw> idw = InterpIdw::New();
  BSHP<VecPt3d> pts(new VecPt3d());
  pts->push_back(Pt3d());
  VecFlt scalar(1, 2);
  idw->SetPts(pts, true);
  idw->SetScalars(&scalar[0], scalar.size());

  pts->resize(0);
  pts->push_back(Pt3d(1, 1, 1));
  scalar.resize(0);
  idw->SetMultiThreading(false);
  idw->InterpToPts(*pts, scalar);
  TS_ASSERT_EQUALS(2.0f, scalar[0]);
} // InterpIdwUnitTests::testInterpToPts
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testInterp2d_a()
{
  BSHP<VecPt3d> pts(new VecPt3d());
  VecPt3d ipts;
  BSHP<VecFlt> scalar(new VecFlt());
  VecFlt s, vBase;
  iGetPoints(*pts, *scalar, ipts);

  BSHP<InterpIdw> idw = InterpIdw::New();
  idw->SetPts(pts, true);
  idw->SetScalars(scalar);
  idw->SetMultiThreading(false);
  idw->InterpToPts(ipts, s);
  vBase = {4.5606198f, 2.9842966f, 4.8022637f, 2.8352659f, 0.0f};
  TS_ASSERT_DELTA_VEC(vBase, s, 1e-7);
} // InterpIdwUnitTests::testInterp2d_a
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testInterp2d_b()
{
  BSHP<VecPt3d> pts(new VecPt3d());
  VecPt3d ipts;
  VecFlt scalar, s(5), vBase;
  iGetPoints(*pts, scalar, ipts);
  BSHP<InterpIdw> idw = InterpIdw::New();
  idw->SetPts(pts, true);
  idw->SetScalars(&scalar[0], scalar.size());
  idw->SetSearchOpts(-1, false);
  for (int i = 0; i < 5; ++i)
    s[i] = idw->InterpToPt(ipts[i]);
  vBase = {10.172515f, 7.9736009f, 9.5391436f, 7.2670202f, 0.0f};
  TS_ASSERT_DELTA_VEC(vBase, s, 1e-7);
} // InterpIdwUnitTests::testInterp2d_b
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testInterp2d_c()
{
  BSHP<VecPt3d> pts(new VecPt3d());
  VecPt3d ipts;
  VecFlt scalar;
  iGetPoints(*pts, scalar, ipts);
  BSHP<InterpIdw> idw = InterpIdw::New();
  idw->SetPts(pts, true);
  idw->SetScalars(&scalar[0], scalar.size());
  idw->SetSearchOpts(5, false);
  VecInt idx;
  VecDbl wt, wtBase;
  Pt3d myPt;
  idw->InterpWeights(myPt, idx, wt);
  wtBase = {0.0019633969662346436, 0.064476394365326525, 0.0, 0.014616102731534851,
            0.91894410593690401};
  TS_ASSERT_DELTA_VEC(wtBase, wt, 1e-7);
} // InterpIdwUnitTests::testInterp2d_c
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testInterp2d_d()
{
  BSHP<VecPt3d> pts(new VecPt3d());
  VecPt3d ipts;
  VecFlt scalar, s(5), vBase;
  iGetPoints(*pts, scalar, ipts);
  BSHP<InterpIdw> idw = InterpIdw::New();
  idw->SetPts(pts, true);
  idw->SetScalars(&scalar[0], scalar.size());
  idw->SetSearchOpts(6, true);
  idw->SetPower(3);
  idw->SetNodalFunction(InterpIdw::GRAD_PLANE, 6, true, NULL);
  for (int i = 0; i < 5; ++i)
    s[i] = idw->InterpToPt(ipts[i]);
  vBase = {6.1976342f, 0.60451090f, -0.80145311f, 13.252287f, 0.0f};
  TS_ASSERT_DELTA_VEC(vBase, s, 1e-7);
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testInterp2d_e()
{
  BSHP<VecPt3d> pts(new VecPt3d());
  VecPt3d ipts;
  VecFlt scalar, s(5), vBase;
  iGetPoints(*pts, scalar, ipts);
  BSHP<InterpIdw> idw = InterpIdw::New();
  idw->SetPts(pts, true);
  idw->SetScalars(&scalar[0], scalar.size());
  idw->SetWeightCalcMethod(InterpIdw::CLASSIC);
  for (int i = 0; i < 5; ++i)
    s[i] = idw->InterpToPt(ipts[i]);
  vBase = {16.884138f, 15.516511f, 16.895596f, 14.160349f, 0.0f};
  TS_ASSERT_DELTA_VEC(vBase, s, 1e-7);
  idw->SetPower(1);
  for (int i = 0; i < 5; ++i)
    s[i] = idw->InterpToPt(ipts[i]);
  vBase = {21.907907f, 21.933628f, 22.670679f, 19.839111f, 0.0f};
  TS_ASSERT_DELTA_VEC(vBase, s, 1e-7);
} // InterpIdwUnitTests::testInterp2d_e
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testInterp2d_f()
{
  BSHP<VecPt3d> pts(new VecPt3d());
  VecPt3d ipts;
  VecFlt scalar, s(5), vBase;
  iGetPoints(*pts, scalar, ipts);
  BSHP<InterpIdw> idw = InterpIdw::New();
  idw->SetPts(pts, true);
  idw->SetScalars(&scalar[0], scalar.size());
  idw->SetNodalFunction(InterpIdw::QUADRATIC, XM_NONE, true, NULL);
  for (int i = 0; i < 5; ++i)
    s[i] = idw->InterpToPt(ipts[i]);
  vBase = {-46.129936f, -65.856499f, -71.997749f, -47.643417f, 0.0f};
  TS_ASSERT_DELTA_VEC(vBase, s, 1e-7);
} // InterpIdwUnitTests::testInterp2d_f
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testInterp3d()
{
  BSHP<VecPt3d> pts(new VecPt3d());
  VecPt3d ipts;
  VecFlt scalar, s(5), vBase;
  iGetPoints(*pts, scalar, ipts);
  BSHP<InterpIdw> idw = InterpIdw::New();
  idw->SetPts(pts, false);
  idw->SetScalars(&scalar[0], scalar.size());
  idw->SetNodalFunction(InterpIdw::GRAD_PLANE, XM_NONE, true, NULL);
  for (int i = 0; i < 5; ++i)
    s[i] = idw->InterpToPt(ipts[i]);
  vBase = {1.5025957e-008f, -7.7562149e-008f, 1.4442284e-007f, -6.1540938e-008f, 0.0f};
  TS_ASSERT_DELTA_VEC(vBase, s, 1e-7);
  idw->SetNodalFunction(InterpIdw::QUADRATIC, XM_NONE, true, NULL);
  for (int i = 0; i < 5; ++i)
    s[i] = idw->InterpToPt(ipts[i]);
  vBase = {5.2535050e-008f, -1.2138121e-007f, 4.1972626e-007f, -1.1271275e-007f, 0.0f};
  TS_ASSERT_DELTA_VEC(vBase, s, 1e-7);
} // InterpIdwUnitTests::testInterp2d_f
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testErrors()
{
  BSHP<VecPt3d> pts(new VecPt3d());
  VecPt3d ipts;
  VecFlt scalar, s(5), vBase;
  *pts = {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {3, 0, 0}};
  scalar.assign(pts->size(), 0);
  ipts.assign(5, Pt3d(-1, 0, 0));

  XmLog::Instance().GetAndClearStackStr();
  InterpIdwImpl idw;
  idw.SetPts(pts, false);
  idw.SetScalars(&scalar[0], scalar.size());
  idw.SetNodalFunction(InterpIdw::QUADRATIC, 1, true, NULL);
  // no nodal function created on previous line because all scalars are equal

  idw.CreateNodalFunctionPtr(InterpIdw::QUADRATIC, 1, true, NULL);
  TS_ASSERT_EQUALS(true, (XmLog::Instance().ErrCount() != 0));
  XmLog::Instance().GetAndClearStackStr();
} // InterpIdwUnitTests::testErrors
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpIdwUnitTests::testErrors2()
{
  BSHP<VecPt3d> pts(new VecPt3d());
  VecPt3d ipts;
  VecFlt scalar, s(5), vBase;
  *pts = {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {3, 0, 0}};
  scalar.assign(pts->size(), 0);
  ipts.assign(5, Pt3d(-1, 0, 0));

  XmLog::Instance().GetAndClearStackStr();
  InterpIdwImpl idw;
  idw.SetPts(pts, true);
  idw.SetScalars(&scalar[0], scalar.size());
  idw.SetNodalFunction(InterpIdw::GRAD_PLANE, 2, true, NULL);
  // no nodal function created on previous line because all scalars are equal

  idw.CreateNodalFunctionPtr(InterpIdw::GRAD_PLANE, 2, true, NULL);
  for (int i = 0; i < 5; ++i)
    s[i] = idw.InterpToPt(ipts[i]);
  TS_ASSERT_EQUALS(true, (XmLog::Instance().ErrCount() != 0));
  XmLog::Instance().GetAndClearStackStr();
} // InterpIdwUnitTests::testErrors

//} // namespace xms
#endif // CXX_TEST