MeMeshUtils.cpp

Go to the documentation of this file.

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

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

// 1. Precompiled header

// 2. My own header
#include <xmsmesh/meshing/MeMeshUtils.h>

// 3. Standard library headers
#include <map>
#include <sstream>

// 4. External library headers

// 5. Shared code headers
#include <xmscore/math/math.h>
#include <xmscore/misc/XmError.h>
#include <xmscore/misc/DynBitset.h>
#include <xmsinterp/triangulate/TrTin.h>

// 6. Non-shared code headers

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

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

//----- Namespace declaration --------------------------------------------------
namespace xms
{
//----- Constants / Enumerations -----------------------------------------------

//----- Classes / Structs ------------------------------------------------------

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
class SmoothIo
{
public:
  SmoothIo()
  : m_tin()
  , m_sizes(nullptr)
  , m_anchorType(0)
  , m_ptsFlag()
  , m_smoothSize(nullptr)
  , m_checkMinSize(false)
  , m_sizeRatio(0)
  , m_minSize(0)
  , m_scaleFactor(0)
  , m_percentGrowth(0)
  , m_logPercentGrowth(0)
  , m_maxSize(0)
  {
  }

  bool CalcMaxSize(double a_length, float a_smoothVal, double& a_maxSize);
  double CalcMinSize(double a_length, float a_smoothVal, double a_calcMaxSize);

  BSHP<TrTin> m_tin;     
  const VecFlt* m_sizes; 
  int m_anchorType;      
  DynBitset m_ptsFlag;   
  VecFlt* m_smoothSize;  

  // used with size function smoothing
  bool m_checkMinSize;       
  double m_sizeRatio;        
  double m_minSize;          
  double m_scaleFactor;      
  double m_percentGrowth;    
  double m_logPercentGrowth; 

  double m_maxSize; 
};

//----- Internal functions -----------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool SmoothIo::CalcMaxSize(double a_length, float a_smoothVal, double& a_maxSize)
{
  if (m_checkMinSize)
  {
    double sizeOneAway = a_length / (m_scaleFactor * a_smoothVal);
    double negLogCheck = sizeOneAway * (m_percentGrowth - 1.0) + 1.0;
    XM_ENSURE_TRUE(negLogCheck > 0.0, false);
    double numElems = log(negLogCheck) / m_logPercentGrowth;
    a_maxSize = pow(m_percentGrowth, numElems) * a_smoothVal;
  }
  else
  {
    a_maxSize = (double)a_smoothVal + a_length * m_maxSize;
  }
  return true;
} // SmoothIo::CalcMaxSize
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
double SmoothIo::CalcMinSize(double a_length, float a_smoothVal, double a_calcMaxSize)
{
  double minSize(0);
  if (m_checkMinSize)
  {
    minSize = a_smoothVal - (a_calcMaxSize - a_smoothVal);
  }
  else
  {
    minSize = (double)a_smoothVal - a_length * m_maxSize;
  }
  return minSize;
} // SmoothIo::CalcMaxSize
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
static void meiDoSmooth(SmoothIo& a_)
{
  if (a_.m_ptsFlag.empty())
    a_.m_ptsFlag.resize(a_.m_tin->NumPoints(), true);
  XM_ENSURE_TRUE(a_.m_ptsFlag.size() == (size_t)a_.m_tin->NumPoints());

  DynBitset ptsFlag(a_.m_tin->NumPoints(), false);

  VecFlt& smoothSize(*a_.m_smoothSize);
  const VecFlt& sz(*a_.m_sizes);
  // set the smoothed size equal to the incoming size
  smoothSize = sz;
  // sort points based on size
  std::multimap<float, int> mapSizeIdx;
  std::vector<std::multimap<float, int>::iterator> vecIt;
  float val(1.0);
  // if we are anchoring to the max size then make the map key the negative
  // of the size
  if (1 == a_.m_anchorType)
    val = -1.0;
  for (int i = 0; i < (int)sz.size(); ++i)
  {
    vecIt.push_back(mapSizeIdx.insert(std::make_pair(val * sz[i], i)));
  }

  // get references to the Tin
  const VecPt3d& pts(a_.m_tin->Points());
  const VecInt2d& trisAdjToPts(a_.m_tin->TrisAdjToPts());
  const VecInt& tris(a_.m_tin->Triangles());

  // iterate through the points
  auto it = mapSizeIdx.begin();
  auto endIt = mapSizeIdx.end();
  for (; it != endIt; ++it)
  {
    int i = it->second; // index of the point
    if (!a_.m_ptsFlag[i])
      continue;
    if (a_.m_checkMinSize && (double)smoothSize[i] < a_.m_minSize)
    {
      smoothSize[i] = (float)a_.m_minSize;
    }
    const VecInt adjTris(trisAdjToPts[i]); // triangles adjacent to the point
    for (const auto& t : adjTris)
    {
      int tIdx = t * 3;
      int triPts[3];
      // points in the adjacent triangle
      triPts[0] = tris[tIdx];
      triPts[1] = tris[tIdx + 1];
      triPts[2] = tris[tIdx + 2];
      for (int t1 = 0; t1 < 3; ++t1)
      {
        int ix = triPts[t1]; // index of triangle point
        // make sure this is not the vertex in the "i" loop
        if (ix == i)
          continue;
        // make sure we have not already adjusted this point
        if (ptsFlag[ix])
          continue;

        const Pt3d &p0(pts[i]), &p1(pts[ix]);
        // calculate what the min or max size can be based on how close
        // the elements are
        double length = Mdist(p0.x, p0.y, p1.x, p1.y);
        double maxSize(0);
        XM_ENSURE_TRUE(a_.CalcMaxSize(length, smoothSize[i], maxSize));

        switch (a_.m_anchorType)
        {
        case 0: // anchor to the min size
        {
          if (a_.m_checkMinSize)
            XM_ENSURE_TRUE(maxSize > a_.m_minSize);
          XM_ENSURE_TRUE(maxSize >= (double)smoothSize[i]);
          if (maxSize < (double)smoothSize[ix])
          {
            smoothSize[ix] = (float)maxSize;
            ptsFlag.set(ix, true);
          }
          else if (a_.m_checkMinSize && (double)smoothSize[ix] < a_.m_minSize)
          {
            smoothSize[ix] = (float)a_.m_minSize;
            ptsFlag.set(ix, true);
          }
        }
        break;
        case 1: // anchor to the max size
        {
          double minSize = a_.CalcMinSize(length, smoothSize[i], maxSize);
          // double minSize = smoothSize[i] - (maxSize - smoothSize[i]);
          XM_ENSURE_TRUE(minSize < smoothSize[i]);
          if (minSize > (double)smoothSize[ix])
          {
            smoothSize[ix] = (float)minSize;
            ptsFlag.set(ix, true);
          }
          else if (a_.m_checkMinSize && (double)smoothSize[ix] < a_.m_minSize)
          {
            smoothSize[ix] = (float)a_.m_minSize;
            ptsFlag.set(ix, true);
          }
        }
        break;
        default:
          XM_ASSERT(0);
          break;
        }
        // resort the point if it changed size
        if (ptsFlag[ix])
        {
          ptsFlag.set(ix, false);
          mapSizeIdx.erase(vecIt[ix]);
          vecIt[ix] = mapSizeIdx.insert(std::make_pair(val * smoothSize[ix], ix));
        }
      }
    }
  }
} // meiDoSmooth

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void meSizeFunctionFromDepth(const VecDbl& a_depths,
                             VecDbl& a_size,
                             double a_minSize,
                             double a_maxSize)
{
  XM_ENSURE_TRUE(a_minSize > 0);
  XM_ENSURE_TRUE(a_maxSize > a_minSize);
  double minD = *std::min_element(a_depths.begin(), a_depths.end());
  double maxD = *std::max_element(a_depths.begin(), a_depths.end());
  double dDiff = maxD - minD;
  XM_ENSURE_TRUE(dDiff > 0);
  double sDiff = a_maxSize - a_minSize;
  XM_ENSURE_TRUE(sDiff > 0);

  a_size.assign(a_depths.size(), a_minSize);
  int i(0);
  for (auto& d : a_size)
  {
    d += ((a_depths[i++] - minD) / dDiff) * sDiff;
  }
} // meSizeFunctionFromDepth
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void meSmoothSizeFunction(BSHP<TrTin> a_tin,
                          const VecFlt& a_sizes,
                          double a_sizeRatio,
                          double a_minSize,
                          int a_anchorType,
                          const DynBitset& a_ptsFlag,
                          VecFlt& a_smoothSize)
{
  a_smoothSize.resize(0);

  // error checking
  XM_ENSURE_TRUE(!a_sizes.empty());
  XM_ENSURE_TRUE(a_sizeRatio > 0.0);
  XM_ENSURE_TRUE(a_minSize > 0.0);
  double percentGrowth = sqrt(1.0 / a_sizeRatio);
  XM_ENSURE_TRUE(percentGrowth > 0.0);
  double logPercentGrowth = log(percentGrowth);
  double scaleFactor = (percentGrowth - 1.0) / logPercentGrowth;
  XM_ENSURE_TRUE(scaleFactor != 0.0);

  SmoothIo io;
  io.m_tin = a_tin;
  io.m_sizes = &a_sizes;
  io.m_anchorType = a_anchorType;
  io.m_ptsFlag = a_ptsFlag;
  io.m_smoothSize = &a_smoothSize;

  io.m_checkMinSize = true;
  io.m_sizeRatio = a_sizeRatio;
  io.m_minSize = a_minSize;
  io.m_percentGrowth = percentGrowth;
  io.m_logPercentGrowth = logPercentGrowth;
  io.m_scaleFactor = scaleFactor;
  meiDoSmooth(io);
} // meSmoothSizeFunction
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void meSmoothElevBySlope(BSHP<TrTin> a_tin,
                         const VecFlt& a_elevs,
                         double a_maxSlope,
                         int a_anchorType,
                         const DynBitset& a_ptsFlag,
                         VecFlt& a_smoothElevs)
{
  a_smoothElevs.resize(0);

  // error checking
  XM_ENSURE_TRUE(!a_elevs.empty());
  XM_ENSURE_TRUE(a_maxSlope > 0.0);

  SmoothIo io;
  io.m_tin = a_tin;
  io.m_sizes = &a_elevs;
  io.m_anchorType = a_anchorType;
  io.m_ptsFlag = a_ptsFlag;
  io.m_smoothSize = &a_smoothElevs;

  io.m_maxSize = a_maxSlope;

  meiDoSmooth(io);
} // meSmoothSizeFunction
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void meModifyMessageWithPolygonId(int a_polyId, std::string& a_msg)
{
  if (a_polyId > -1)
  {
    std::stringstream ss;
    ss << "Error meshing polygon id: " << a_polyId << ". ";
    a_msg = ss.str() + a_msg;
  }
} // meModifyMessageWithPolygonId

} // namespace xms

#if CXX_TEST

#include <xmsmesh/meshing/MeMeshUtils.t.h>

#include <xmscore/testing/TestTools.h>
#include <xmsinterp/triangulate/TrTriangulatorPoints.h>

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeMeshUtilsUnitTests::testSizeFuncFromDepth()
{
  // array of depths
  xms::VecDbl depths = {0, 5, 10, 20, 25, 5, 0};
  // array for the computed sizes
  xms::VecDbl elemSize;
  // set the value of the min and max element size
  double minElem(2), maxElem(102);
  // generate the size array
  xms::meSizeFunctionFromDepth(depths, elemSize, minElem, maxElem);
  // verify that the sizes are as expected
  xms::VecDbl baseElemSize = {2, 22, 42, 82, 102, 22, 2};
  TS_ASSERT_DELTA_VEC(baseElemSize, elemSize, 1e-9);
} // MeMeshUtilsUnitTests::testSizeFuncFromDepth
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeMeshUtilsUnitTests::testSmoothSizeFunc()
{
  BSHP<xms::VecPt3d> pts(new xms::VecPt3d());
  *pts = {{0, 0},   {10, 0},  {20, 0}, {30, 0},  {0, 10},  {10, 10},
          {20, 10}, {30, 10}, {0, 20}, {10, 20}, {20, 20}, {30, 20}};
  xms::VecFlt sizes(pts->size(), 100);
  sizes[4] = 1;
  BSHP<xms::VecInt> tris(new xms::VecInt());
  BSHP<xms::VecInt2d> adjTris(new xms::VecInt2d());
  xms::TrTriangulatorPoints tr(*pts, *tris, &*adjTris);
  tr.Triangulate();
  BSHP<xms::TrTin> tin = xms::TrTin::New();
  tin->SetGeometry(pts, tris, adjTris);

  xms::VecFlt vSmooth;
  xms::DynBitset ptFlags;
  xms::meSmoothSizeFunction(tin, sizes, 0.5, 1.0, 0, ptFlags, vSmooth);
  xms::VecFlt baseSmooth = {4.46f, 5.90f,  9.36f, 12.83f, 1.0f,   4.46f,
                            7.93f, 11.39f, 4.46f, 7.93f,  11.39f, 14.86f};
  TS_ASSERT_DELTA_VEC(baseSmooth, vSmooth, .1);
} // MeMeshUtilsUnitTests::testSmoothSizeFunc
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeMeshUtilsUnitTests::testSmoothSizeFunc1()
{
  BSHP<xms::VecPt3d> pts(new xms::VecPt3d());
  *pts = {{0, 0},   {10, 0},  {20, 0}, {30, 0},  {0, 10},  {10, 10},
          {20, 10}, {30, 10}, {0, 20}, {10, 20}, {20, 20}, {30, 20}};
  xms::VecFlt sizes(pts->size(), 1);
  sizes[4] = 100;
  BSHP<xms::VecInt> tris(new xms::VecInt());
  BSHP<xms::VecInt2d> adjTris(new xms::VecInt2d());
  xms::TrTriangulatorPoints tr(*pts, *tris, &*adjTris);
  tr.Triangulate();
  BSHP<xms::TrTin> tin = xms::TrTin::New();
  tin->SetGeometry(pts, tris, adjTris);

  xms::VecFlt vSmooth;
  xms::DynBitset ptFlags;
  xms::meSmoothSizeFunction(tin, sizes, 0.5, 1.0, 1, ptFlags, vSmooth);
  xms::VecFlt baseSmooth = {96.53f, 95.10f, 91.63f, 88.17f, 100.0f, 96.53f,
                            93.07f, 89.60f, 96.53f, 93.07f, 89.60f, 86.14f};
  TS_ASSERT_DELTA_VEC(baseSmooth, vSmooth, .1);
} // MeMeshUtilsUnitTests::testSmoothSizeFunc1

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeMeshUtilsUnitTests::testSmoothSizeFunc2()
{
  BSHP<xms::VecPt3d> pts(new xms::VecPt3d());
  *pts = {{0, 0},   {10, 0},  {20, 0}, {30, 0},  {0, 10},  {10, 10},
          {20, 10}, {30, 10}, {0, 20}, {10, 20}, {20, 20}, {30, 20}};
  xms::VecFlt sizes(pts->size(), 100);
  sizes[4] = 1;
  BSHP<xms::VecInt> tris(new xms::VecInt());
  BSHP<xms::VecInt2d> adjTris(new xms::VecInt2d());
  xms::TrTriangulatorPoints tr(*pts, *tris, &*adjTris);
  tr.Triangulate();
  BSHP<xms::TrTin> tin = xms::TrTin::New();
  tin->SetGeometry(pts, tris, adjTris);

  xms::VecFlt vSmooth;
  xms::DynBitset ptsFlag;
  xms::meSmoothElevBySlope(tin, sizes, 0.5, 0, ptsFlag, vSmooth);
  xms::VecFlt baseSmooth = {6.00f,  8.07f,  13.07f, 18.07f, 1.0f,   6.00f,
                            11.00f, 16.00f, 6.00f,  11.00f, 16.00f, 21.00f};
  TS_ASSERT_DELTA_VEC(baseSmooth, vSmooth, .1);
} // MeMeshUtilsUnitTests::testSmoothSizeFunc2
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void MeMeshUtilsUnitTests::testSmoothSizeFunc3()
{
  BSHP<xms::VecPt3d> pts(new xms::VecPt3d());
  *pts = {{0, 0},   {10, 0},  {20, 0}, {30, 0},  {0, 10},  {10, 10},
          {20, 10}, {30, 10}, {0, 20}, {10, 20}, {20, 20}, {30, 20}};
  xms::VecFlt sizes(pts->size(), 1);
  sizes[4] = 100;
  BSHP<xms::VecInt> tris(new xms::VecInt());
  BSHP<xms::VecInt2d> adjTris(new xms::VecInt2d());
  xms::TrTriangulatorPoints tr(*pts, *tris, &*adjTris);
  tr.Triangulate();
  BSHP<xms::TrTin> tin = xms::TrTin::New();
  tin->SetGeometry(pts, tris, adjTris);

  xms::VecFlt vSmooth;
  xms::DynBitset ptsFlag;
  xms::meSmoothElevBySlope(tin, sizes, 0.5, 1, ptsFlag, vSmooth);
  xms::VecFlt baseSmooth = {95.00f, 92.92f, 87.92f, 82.92f, 100.0f, 95.00f,
                            90.00f, 85.00f, 95.00f, 90.00f, 85.00f, 80.00f};
  TS_ASSERT_DELTA_VEC(baseSmooth, vSmooth, .1);
} // MeMeshUtilsUnitTests::testSmoothSizeFunc3

#endif