InterpLinear.cpp

Go to the documentation of this file.

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

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

// 1. Precompiled header

// 2. My header
#include <xmsinterp/interpolate/InterpLinear.h>

// 3. Standard Library Headers
#include <cfloat>
#include <sstream>

// 4. External Library Headers

// 5. Shared Headers
#include <xmscore/math/math.h>
#include <xmsgrid/geometry/GmPtSearch.h>
#include <xmsgrid/geometry/GmTriSearch.h>
#include <xmsinterp/interpolate/detail/InterpCt.h>
#include <xmsinterp/interpolate/detail/InterpNatNeigh.h>
#include <xmsinterp/interpolate/detail/NodalFunc.h>
#include <xmsinterp/interpolate/InterpUtil.h>
#include <xmscore/stl/utility.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 InterpLinearImpl : public InterpLinear
{
public:
  InterpLinearImpl();
  virtual ~InterpLinearImpl();

  virtual void SetPtsTris(BSHP<VecPt3d> a_, BSHP<VecInt> a_tris) override;

  virtual void SetScalars(const float* a_scalar, size_t a_n) override;
  virtual void SetScalars(BSHP<VecFlt> a_scalar) override;
  virtual void SetPtActivity(DynBitset& a_activity) override;
  // bitset is number of triangles in length not numtri*3 like the tris array
  virtual void SetTriActivity(DynBitset& a_activity) 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 VecInt GetExtrapolationPointIndexes() const override;

  virtual float InterpToPt(const Pt3d& a_pt) override;
          float InterpToPtImpl(const Pt3d& a_pt, int a_ptIdx);
  virtual void InterpToPts(const VecPt3d& a_pts, VecFlt& a_scalars) override;

  // find triangle
  virtual int TriContainingPt(const Pt3d& a_pt) override;
  virtual void TriEnvelopsContainingPt(const Pt3d& a_pt, VecInt& a_tris) override;
  // interpolation
  virtual bool InterpWeights(const Pt3d& a_pt, VecInt& a_idxs, VecDbl& a_wts) override;

  float CloughTocherInterp(VecInt& a_idxs, const Pt3d& a_pt);
  float NatNeighInterp(const Pt3d& a_pt);
  // interpolation options
  virtual void SetExtrapVal(double a_val) override;
  virtual void SetTrunc(double a_sMax, double a_sMin) override;
  virtual void SetUseCloughTocher(bool a_, BSHP<Observer> a_prog) override;
  virtual void SetUseNatNeigh(bool a_,
                              int a_nodalFunc,
                              int a_ndFuncOpt,
                              int a_ndFuncNumNearestPts,
                              bool a_blendWeights,
                              BSHP<Observer> a_prog) override;

  virtual double GetExtrapVal() const override { return m_extrap; }
  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 bool GetUseCloughTocher() const override { return m_cloughTocher; }
  virtual bool GetUseNatNeigh() const override { return m_natNeigh; }
  virtual int GetNatNeighNodalFunc() const override
  {
    if (!m_nodal)
      return 0;
    return m_nodal->GetType();
  }
  virtual int GetNatNeighNodalFuncNearestPtsOption() const override
  {
    if (!m_nodal)
      return 0;
    return m_nodal->GetNearestPointsOption();
  }
  virtual int GetNatNeighNodalFuncNumNearestPts() const override
  {
    if (!m_nodal)
      return 0;
    return m_nodal->GetNumNearestPoints();
  }
  virtual bool GetNatNeighBlendWeights() const override
  {
    if (!m_nodal)
      return false;
    return m_nodal->GetUseModifiedShepardWeights();
  }


  BSHP<GmPtSearch> CreatePtSearch();
  void RecalcNodalFunc();
  virtual std::string ToString() const override;

protected:
  BSHP<GmTriSearch> m_triSearch; 
  BSHP<VecPt3d> m_pts;           
  BSHP<VecInt> m_tris;           
  BSHP<VecFlt> m_scalar;         
  double m_extrap;               
  bool m_trunc;                  
  double m_truncMax;             
  double m_truncMin;             
  bool m_cloughTocher;           
  bool m_natNeigh;           
  BSHP<NodalFunc> m_nodal;   
  BSHP<InterpCt> m_ct;       
  BSHP<InterpNatNeigh> m_nn; 
  VecInt m_extrapolationPointIndexes; 
};

//----- Internal functions -----------------------------------------------------

//----- Class / Function definitions -------------------------------------------

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

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
InterpLinearImpl::InterpLinearImpl()
: m_triSearch(GmTriSearch::New())
, m_pts(new VecPt3d())
, m_tris(new VecInt())
, m_scalar(new VecFlt())
, m_extrap(XM_NODATA)
, m_trunc(false)
, m_truncMax(0.0)
, m_truncMin(0.0)
, m_cloughTocher(false)
, m_natNeigh(false)
{
} // InterpLinearImpl::InterpLinearImpl
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
InterpLinearImpl::~InterpLinearImpl()
{
} // InterpLinearImpl::~InterpLinearImpl
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::SetPtsTris(BSHP<VecPt3d> a_pts, BSHP<VecInt> a_tris)
{
  m_pts = a_pts;
  m_tris = a_tris;
  m_scalar->assign(m_pts->size(), 0);
  for (size_t i = 0; i < m_pts->size(); ++i)
    (*m_scalar)[i] = (float)(*m_pts)[i].z;
  m_triSearch->TrisToSearch(a_pts, a_tris);
} // InterpLinearImpl::SetPtsTris
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::SetScalars(const float* a_scalar, size_t a_n)
{
  XM_ENSURE_TRUE(a_n == m_pts->size());
  m_scalar->assign(&a_scalar[0], &a_scalar[a_n]);
  RecalcNodalFunc();
} // InterpLinearImpl::SetScalars
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::SetScalars(BSHP<VecFlt> a_scalar)
{
  XM_ENSURE_TRUE(a_scalar->size() == m_pts->size());
  m_scalar = a_scalar;
  RecalcNodalFunc();
} // InterpLinear::SetScalars
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::SetPtActivity(DynBitset& a_activity)
{
  m_triSearch->SetPtActivity(a_activity);
  RecalcNodalFunc();
} // InterpLinearImpl::SetPtActivity
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::SetTriActivity(DynBitset& a_activity)
{
  m_triSearch->SetTriActivity(a_activity);
  RecalcNodalFunc();
} // InterpLinearImpl::SetTriActivity
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const BSHP<VecPt3d> InterpLinearImpl::GetPts() const
{
  return m_pts;
} // InterpLinearImpl::GetPts
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const BSHP<VecInt> InterpLinearImpl::GetTris() const
{
  return m_tris;
} // InterpLinearImpl::GetTris
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const BSHP<VecFlt> InterpLinearImpl::GetScalars() const
{
  return m_scalar;
} // InterpLinearImpl::GetScalars
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
DynBitset InterpLinearImpl::GetPtActivity() const
{
  return m_triSearch->GetPtActivity();
} // InterpLinearImpl::GetPtActivity
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
DynBitset InterpLinearImpl::GetTriActivity() const
{
  return m_triSearch->GetTriActivity();
} // InterpLinearImpl::GetTriActivity
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
VecInt InterpLinearImpl::GetExtrapolationPointIndexes() const
{
  return m_extrapolationPointIndexes;
} // InterpLinearImpl::GetExtrapolationPointIndexes
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
float InterpLinearImpl::InterpToPt(const Pt3d& a_pt)
{
  return InterpToPtImpl(a_pt, 0);
} // InterpLinearImpl::InterpToPt
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
float InterpLinearImpl::InterpToPtImpl(const Pt3d& a_pt, int a_ptIdx)
{
  // return m_triSearch->InterpToPt(a_pt);
  float rval((float)m_extrap);
  VecInt idxs(3, 0);
  VecDbl wts(3, 0);
  if (m_triSearch->InterpWeights(a_pt, idxs, wts))
  {
    float* s = &(*m_scalar)[0];
    rval = (float)((s[idxs[0]] * wts[0]) + (s[idxs[1]] * wts[1]) + (s[idxs[2]] * wts[2]));
    if (m_cloughTocher || m_natNeigh)
    {
      if (EQ_TOL(wts[0], 1, FLT_EPSILON) || EQ_TOL(wts[1], 1, FLT_EPSILON) ||
          EQ_TOL(wts[2], 1, FLT_EPSILON))
      {
      }
      else
      {
        if (m_cloughTocher)
        {
          rval = CloughTocherInterp(idxs, a_pt);
        }
        else
          rval = NatNeighInterp(a_pt);
      }
    }
  }
  else
  {
    m_extrapolationPointIndexes.push_back(a_ptIdx);
  }
  if (m_trunc)
  {
    if (rval < m_truncMin)
      rval = (float)m_truncMin;
    if (rval > m_truncMax)
      rval = (float)m_truncMax;
  }
  return rval;
} // InterpLinearImpl::InterpToPt
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
void InterpLinearImpl::InterpToPts(const VecPt3d& a_pts, VecFlt& a_scalars)
{
  // m_triSearch->InterpToPts(a_pts, a_scalars);
  m_extrapolationPointIndexes.clear();
  a_scalars.assign(a_pts.size(), 0);
  for (size_t i = 0; i < a_pts.size(); i++)
  {
    a_scalars[i] = InterpToPtImpl(a_pts[i], static_cast<int>(i));
  }
} // InterpLinearImpl::InterpToPts
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int InterpLinearImpl::TriContainingPt(const Pt3d& a_pt)
{
  return m_triSearch->TriContainingPt(a_pt);
} // InterpLinearImpl::TriContainingPt
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::TriEnvelopsContainingPt(const Pt3d& a_pt, VecInt& a_tris)
{
  m_triSearch->TriEnvelopsContainingPt(a_pt, a_tris);
} // InterpLinearImpl::TriContainingPt
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool InterpLinearImpl::InterpWeights(const Pt3d& a_pt, VecInt& a_idxs, VecDbl& a_wts)
{
  return m_triSearch->InterpWeights(a_pt, a_idxs, a_wts);
} // InterpLinearImpl::InterpWeights
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
float InterpLinearImpl::CloughTocherInterp(VecInt& a_idxs, const Pt3d& a_pt)
{
  m_ct->ComputeCtCoeff(&a_idxs[0]);
  return (float)(m_ct->InterpToPt(a_pt));
} // InterpLinearImpl::CloughTocherInterp
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
float InterpLinearImpl::NatNeighInterp(const Pt3d& a_pt)
{
  return m_nn->InterpToPt(a_pt);
} // InterpLinearImpl::NatNeighInterp
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::SetExtrapVal(double a_val)
{
  m_extrap = a_val;
} // InterpLinearImpl::SetExtrapVal
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::SetTrunc(double a_sMax, double a_sMin)
{
  m_trunc = true;
  m_truncMax = a_sMax;
  m_truncMin = a_sMin;
} // InterpLinearImpl::SetTrunc
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::SetUseCloughTocher(bool a_, BSHP<Observer> a_prog)
{
  // m_triSearch->SetUseCloughTocher(a_, a_prog);
  m_cloughTocher = a_;
  if (!m_cloughTocher)
    return;

  BSHP<GmPtSearch> pts = CreatePtSearch();
  int nodalfunc(1), nearest(16);
  bool twod(true), quad_oct(false), modifiedShepardWeights(true);
  double dPower(2);
  m_nodal.reset();
  if (!inAllScalarsEqual(*m_scalar, m_triSearch->GetPtActivity()) && m_pts->size() > 6)
  {
    m_nodal = NodalFunc::New(nodalfunc, twod, pts, *m_pts, *m_scalar, nearest, quad_oct, dPower,
                             modifiedShepardWeights, a_prog, nullptr);
  }
  m_ct.reset(new InterpCt(*m_pts, m_nodal));
} // InterpLinearImpl::SetUseCloughTocher
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::SetUseNatNeigh(bool a_,
                                      int a_nodalFunc,
                                      int a_ndFuncOpt,
                                      int a_ndFuncNumNearestPts,
                                      bool a_blendWeights,
                                      BSHP<Observer> a_prog)
{
  m_natNeigh = a_;
  if (!m_natNeigh)
    return;

  m_nn = InterpNatNeigh::New(*m_pts, *m_tris, *m_scalar, m_triSearch.get());
  m_nn->SetBlendWeights(a_blendWeights);

  if (0 != a_nodalFunc && !inAllScalarsEqual(*m_scalar, m_triSearch->GetPtActivity()))
  {
    int nearest(a_ndFuncNumNearestPts);
    bool twod(true), quad_oct(false), modifiedShepardWeights(true);
    double dPower(2);
    InterpNatNeigh* nnPtr(nullptr);
    BSHP<GmPtSearch> pts;
    if (0 != a_ndFuncOpt)
      pts = CreatePtSearch();
    else
      nnPtr = m_nn.get();

    m_nodal = NodalFunc::New(a_nodalFunc, twod, pts, *m_pts, *m_scalar, nearest, quad_oct, dPower,
                             modifiedShepardWeights, a_prog, nnPtr);
    m_nn->SetNodalFunc(m_nodal);
  }
} // InterpLinearImpl::SetUseNatNeigh
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
BSHP<GmPtSearch> InterpLinearImpl::CreatePtSearch()
{
  BSHP<GmPtSearch> pts = GmPtSearch::New(true);
  pts->PtsToSearch(m_pts);
  DynBitset b = m_triSearch->GetPtActivity();
  pts->SetActivity(b);
  return pts;
} // InterpLinearImpl::CreatePtSearch
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearImpl::RecalcNodalFunc()
{
  if (m_nn)
    m_nn->RecalcNodalFunc();
  else if (m_ct)
    m_ct->RecalcNodalFunc();
} // InterpLinearImpl::RecalcNodalFunc
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
std::string InterpLinearImpl::ToString() const
{
  std::stringstream ss;
  ss << m_cloughTocher << "=cloughTocher " << m_natNeigh << "=natNeigh " << m_trunc << "=trunc "
     << m_truncMax << "=truncMax " << m_truncMin << "=truncMin " << m_extrap << "=extrap "
     << "\n";
  if (m_triSearch)
  {
    ss << m_triSearch->ToString();
    ss << "=triSearch "
       << "\n";
  }
  if (m_pts)
    VecToStream(ss, *m_pts, "pts");
  if (m_tris)
    VecToStream(ss, *m_tris, "tris");
  if (m_scalar)
    VecToStream(ss, *m_scalar, "scalar");
  if (m_nodal)
  {
    ss << m_nodal->ToString();
    ss << "=nodal "
       << "\n";
  }
  if (m_ct)
  {
    ss << m_ct->ToString();
    ss << "=ct "
       << "\n";
  }
  if (m_nn)
  {
    ss << m_nn->ToString();
    ss << "=nn "
       << "\n";
  }
  ss << "\n";
  return ss.str();
} // InterpLinearImpl::ToString

} // namespace xms

#ifdef CXX_TEST

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

// namespace xms {
using namespace xms;

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearUnitTests::testToString()
{
#if 0
  BSHP<InterpLinear> interpolator = InterpLinear::New();
  std::string s = interpolator->ToString();
  std::string expected = "0=cloughTocher 0=natNeigh 0=trunc 0=truncMax 0=truncMin 3.40282e+038,3.40282e+038,3.40282e+038=min -3.40282e+038,-3.40282e+038,-3.40282e+038=max -1e+007=extrap ";
  TS_ASSERT_EQUALS(expected, s);
#endif
} // InterpLinearUnitTests::testToString
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearUnitTests::testOneTriangle()
{
  InterpLinearImpl linear;
  BSHP<VecPt3d> p(new VecPt3d());
  *p = {{0, 0, 1}, {1, 1, 1}, {1, 0, 1}};
  BSHP<VecInt> t(new VecInt());
  *t = {0, 2, 1};
  linear.SetPtsTris(p, t);
  const float EXTRAPVAL = (float)-3.12345678;
  linear.SetExtrapVal((double)EXTRAPVAL);

  Pt3d pt(-.1, 0, 0);
  // double val(0);
  // outside the triangle
  TS_ASSERT_DELTA(EXTRAPVAL, linear.InterpToPt(pt), FLT_EPSILON);
  pt = Pt3d(-1e-6, 0, 0);
  TS_ASSERT_DELTA(EXTRAPVAL, linear.InterpToPt(pt), FLT_EPSILON);
  // inside the triangle
  pt = Pt3d(.5, .25, 0);
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);
  pt = Pt3d(-1e-13, 0, 0);
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);
  // on each edge of the triangle
  pt = Pt3d(.5, .5, 0);
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);
  pt = Pt3d(.5, 0, 0);
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);
  pt = Pt3d(1, .5, 0);
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);
  // at the points of the triangle
  pt = Pt3d();
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);
  pt = Pt3d(1, 1, 0);
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);
  pt = Pt3d(1, 0, 0);
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);

  VecInt idxes = linear.GetExtrapolationPointIndexes();
  VecInt baseIdxes = { 0, 0 };
  TS_ASSERT_EQUALS(baseIdxes, idxes);

  VecPt3d vPts = {{ .5, .5, 0 }, { -1, 0, 0 }, { -.5, 0, 0 }};
  VecFlt scalar;
  linear.InterpToPts(vPts, scalar);
  idxes = linear.GetExtrapolationPointIndexes();
  baseIdxes = { 1, 2 };
  TS_ASSERT_EQUALS(baseIdxes, idxes);
} // InterpLinearUnitTests::testOneTriangle
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearUnitTests::testInterp()
{
  InterpLinearImpl linear;
  BSHP<VecPt3d> p(new VecPt3d());
  *p = {{0, 0, 0}, {1, 1, 1}, {1, 0, 2}};
  BSHP<VecInt> t(new VecInt());
  *t = {0, 2, 1};
  linear.SetPtsTris(p, t);

  Pt3d pt(.5, .5, 0);
  // double val(0);
  TS_ASSERT_DELTA(.5, linear.InterpToPt(pt), FLT_EPSILON);
  pt = Pt3d(.5, 0, 0);
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);
  pt = Pt3d(1, .5, 0);
  TS_ASSERT_DELTA(1.5, linear.InterpToPt(pt), FLT_EPSILON);
  pt = Pt3d(.66666666666666666, .3333333333333333333, 0);
  TS_ASSERT_DELTA(1, linear.InterpToPt(pt), FLT_EPSILON);
} // InterpLinearUnitTests::testInterp
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void InterpLinearUnitTests::testCt()
{
  InterpLinearImpl linear;
  BSHP<VecPt3d> p(new VecPt3d());
  *p = {{0, 0, 0}, {1, 1, 1}, {1, 0, 2}};
  BSHP<VecInt> t(new VecInt());
  *t = {0, 2, 1};
  linear.SetPtsTris(p, t);

  Pt3d pt(.5, .5, 0);
  linear.SetUseCloughTocher(true, BSHP<Observer>());
  TS_ASSERT_DELTA(.5f, linear.InterpToPt(pt), 1e-7);
  linear.SetTrunc(-1, 0);
  TS_ASSERT_DELTA(-1.0f, linear.InterpToPt(pt), FLT_EPSILON);
  pt.x = .5;
  pt.y = .1;
  TS_ASSERT_DELTA(-1.0f, linear.InterpToPt(pt), FLT_EPSILON);
  pt.x = .2;
  pt.y = .15;
  TS_ASSERT_DELTA(-1.0f, linear.InterpToPt(pt), FLT_EPSILON);
  pt.x = .9;
  pt.y = .5;
  TS_ASSERT_DELTA(-1.0f, linear.InterpToPt(pt), FLT_EPSILON);
} // InterpLinearUnitTests::testCt

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