XmUGridTriangles2d.cpp

Go to the documentation of this file.

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

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

// 1. Precompiled header

// 2. My own header
#include <xmsextractor/ugrid/XmUGridTriangles2d.h>

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

// 4. External library headers
#include <boost/container/flat_map.hpp>

// 5. Shared code headers
#include <xmscore/misc/XmError.h>
#include <xmscore/misc/XmLog.h>
#include <xmsgrid/geometry/geoms.h>
#include <xmsgrid/geometry/GmTriSearch.h>
#include <xmsgrid/ugrid/XmUGrid.h>

// 6. Non-shared code headers

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

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

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

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

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

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

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

namespace
{
typedef std::tuple<int, int, int> Triangle; 
typedef boost::container::flat_map<Triangle, double>
  RatioCache; 
typedef boost::container::flat_map<Triangle, bool>
  ValidCache; 

class XmUGridTriangles2dImpl : public XmUGridTriangles2d
{
public:
  XmUGridTriangles2dImpl();

  virtual void BuildTriangles(const XmUGrid& a_ugrid, bool a_addTriangleCenters) override;
  virtual void BuildEarcutTriangles(const XmUGrid& a_ugrid) override;
  virtual void SetCellActivity(const DynBitset& a_cellActivity) override;

  virtual const VecPt3d& GetPoints() const override;
  virtual const VecInt& GetTriangles() const override;
  virtual BSHP<VecPt3d> GetPointsPtr() override;
  virtual BSHP<VecInt> GetTrianglesPtr() override;

  int AddCellCentroid(int a_cellIdx, const Pt3d& a_point);
  void AddCellTriangle(int a_cellIdx, int a_idx1, int a_idx2, int a_idx3);

  virtual int GetCellCentroid(int a_cellIdx) const override;

  virtual int GetIntersectedCell(const Pt3d& a_point, VecInt& a_idxs, VecDbl& a_weights) override;

private:
  void Initialize(const XmUGrid& a_ugrid);
  BSHP<GmTriSearch> GetTriSearch();

  BSHP<VecPt3d> m_points;                
  BSHP<VecInt> m_triangles;              
  VecInt m_centroidIdxs;                 
  VecInt m_triangleToCellIdx;            
  mutable BSHP<GmTriSearch> m_triSearch; 
};

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
double iMagnitude(const Pt3d& a_vec)
{
  double magnitude = sqrt(a_vec.x * a_vec.x + a_vec.y * a_vec.y + a_vec.z * a_vec.z);
  return magnitude;
} // iMagnitude
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
double iGetEarcutTriangleRatio(const VecPt3d& a_points,
                               int a_idx1,
                               int a_idx2,
                               int a_idx3,
                               RatioCache& a_ratioCache)
{
  Triangle triangle(a_idx1, a_idx2, a_idx3);
  auto ratioIter = a_ratioCache.find(triangle);
  if (ratioIter != a_ratioCache.end())
    return ratioIter->second;

  Pt3d pt1 = a_points[a_idx1];
  Pt3d pt2 = a_points[a_idx2];
  Pt3d pt3 = a_points[a_idx3];
  Pt3d v1 = pt1 - pt2;
  Pt3d v2 = pt3 - pt2;
  Pt3d v3 = pt3 - pt1;

  // get 2*area
  Pt3d cross;
  gmCross3D(v2, v1, &cross);
  double ratio;
  if (cross.z <= 0.0)
  {
    // inverted
    ratio = -1.0;
  }
  else
  {
    double area2 = iMagnitude(cross);
    if (area2 == 0.0)
    {
      // degenerate triangle
      ratio = -2.0;
    }
    else
    {
      double perimeter = iMagnitude(v1) + iMagnitude(v2) + iMagnitude(v3);
      ratio = perimeter * perimeter / area2;
    }
  }

  a_ratioCache[triangle] = ratio;
  return ratio;
} // iGetEarcutTriangleRatio
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool iValidTriangle(const VecPt3d& a_points,
                    const VecInt a_polygon,
                    int a_idx1,
                    int a_idx2,
                    int a_idx3,
                    ValidCache& a_validCache)
{
  Triangle triangle(a_idx1, a_idx2, a_idx3);
  auto validIter = a_validCache.find(triangle);
  if (validIter != a_validCache.end())
    return validIter->second;

  const Pt3d& pt1 = a_points[a_idx1];
  const Pt3d& pt2 = a_points[a_idx2];
  const Pt3d& pt3 = a_points[a_idx3];
  for (size_t pointIdx = 0; pointIdx < a_polygon.size(); ++pointIdx)
  {
    int idx = a_polygon[pointIdx];
    if (idx != a_idx1 && idx != a_idx2 && idx != a_idx3)
    {
      const Pt3d& pt = a_points[idx];
      if (gmTurn(pt1, pt2, pt, 0.0) == TURN_LEFT && gmTurn(pt2, pt3, pt) == TURN_LEFT &&
          gmTurn(pt3, pt1, pt, 0.0) == TURN_LEFT)
      {
        a_validCache[triangle] = false;
        return false;
      }
    }
  }

  a_validCache[triangle] = true;
  return true;
} // iValidTriangle
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
bool iGenerateCentroidTriangles(XmUGridTriangles2dImpl& a_ugridTris,
                                int a_cellIdx,
                                const VecInt& a_polygonIdxs)
{
  const VecPt3d& points = a_ugridTris.GetPoints();
  size_t numPoints = a_polygonIdxs.size();

  VecPt3d polygon;
  for (size_t pointIdx = 0; pointIdx < numPoints; ++pointIdx)
    polygon.push_back(points[a_polygonIdxs[pointIdx]]);

  Pt3d centroid = gmComputeCentroid(polygon);

  // make sure none of the triangles connected to the centroid are inverted
  for (size_t pointIdx = 0; pointIdx < polygon.size(); ++pointIdx)
  {
    const Pt3d& pt1 = polygon[pointIdx];
    const Pt3d& pt2 = polygon[(pointIdx + 1) % numPoints];

    // centroid should be to the left of each edge
    if (gmTurn(pt1, pt2, centroid, 0.0) != TURN_LEFT)
      return false;
  }

  // add centroid to list of points
  int centroidIdx = a_ugridTris.AddCellCentroid(a_cellIdx, centroid);

  // add triangles
  for (size_t pointIdx = 0; pointIdx < polygon.size(); ++pointIdx)
  {
    int idx1 = a_polygonIdxs[pointIdx];
    int idx2 = a_polygonIdxs[(pointIdx + 1) % numPoints];
    a_ugridTris.AddCellTriangle(a_cellIdx, idx1, idx2, centroidIdx);
  }

  return true;
} // iGenerateCentroidTriangles
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void iBuildEarcutTriangles(XmUGridTriangles2dImpl& a_ugridTris,
                           int a_cellIdx,
                           const VecInt& a_polygonIdxs)
{
  VecInt polygonIdxs = a_polygonIdxs;
  const VecPt3d& points = a_ugridTris.GetPoints();

  // continually find best triangle on adjacent edges and cut it off polygon
  RatioCache ratioCache;
  ValidCache validCache;
  while (polygonIdxs.size() >= 4)
  {
    int bestIdx = -1;
    int secondBestIdx = -1;

    int numPoints = (int)polygonIdxs.size();
    double bestRatio = std::numeric_limits<double>::max();
    double secondBestRatio = std::numeric_limits<double>::max();
    for (int pointIdx = 0; pointIdx < numPoints; ++pointIdx)
    {
      int idx1 = polygonIdxs[(pointIdx + numPoints - 1) % numPoints];
      int idx2 = polygonIdxs[pointIdx];
      int idx3 = polygonIdxs[(pointIdx + 1) % numPoints];

      // make sure triangle is valid (not inverted and doesn't have other points in it)
      double ratio = iGetEarcutTriangleRatio(points, idx1, idx2, idx3, ratioCache);
      if (ratio > 0.0)
      {
        if (iValidTriangle(points, polygonIdxs, idx1, idx2, idx3, validCache))
        {
          if (ratio < bestRatio)
          {
            secondBestRatio = bestRatio;
            bestRatio = ratio;
            secondBestIdx = bestIdx;
            bestIdx = pointIdx;
          }
          else if (ratio < secondBestRatio)
          {
            secondBestRatio = ratio;
            secondBestIdx = pointIdx;
          }
        }
      }
    }

    if (bestIdx >= 0)
    {
      if (numPoints == 4 && secondBestIdx >= 0)
        bestIdx = secondBestIdx;

      // cut off the ear triangle
      int polygonIdx1 = (bestIdx + numPoints - 1) % numPoints;
      int polygonIdx2 = bestIdx;
      int polygonIdx3 = (bestIdx + 1) % numPoints;
      int idx1 = polygonIdxs[polygonIdx1];
      int idx2 = polygonIdxs[polygonIdx2];
      int idx3 = polygonIdxs[polygonIdx3];
      a_ugridTris.AddCellTriangle(a_cellIdx, idx1, idx2, idx3);
      polygonIdxs.erase(polygonIdxs.begin() + polygonIdx2);
    }
    else
    {
      std::ostringstream ss;
      ss << "Unable to split cell number " << a_cellIdx + 1 << " into triangles.";
      XM_LOG(xmlog::error, ss.str());
      return;
    }
  }

  // push on remaining triangle
  a_ugridTris.AddCellTriangle(a_cellIdx, polygonIdxs[0], polygonIdxs[1], polygonIdxs[2]);
} // iBuildEarcutTriangles

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridTriangles2dImpl::XmUGridTriangles2dImpl()
: m_points(new VecPt3d)
, m_triangles(new VecInt)
, m_centroidIdxs()
, m_triangleToCellIdx()
, m_triSearch()
{
} // XmUGridTriangles2dImpl::XmUGridTriangles2dImpl
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dImpl::BuildTriangles(const XmUGrid& a_ugrid, bool a_addTriangleCenters)
{
  Initialize(a_ugrid);

  int numCells = a_ugrid.GetCellCount();
  VecInt cellPoints;
  for (int cellIdx = 0; cellIdx < numCells; ++cellIdx)
  {
    a_ugrid.GetCellPoints(cellIdx, cellPoints);
    bool builtTriangles = false;
    if (a_addTriangleCenters)
      builtTriangles = iGenerateCentroidTriangles(*this, cellIdx, cellPoints);
    if (!builtTriangles)
      iBuildEarcutTriangles(*this, cellIdx, cellPoints);
  }
} // XmUGridTriangles2dImpl::BuildTriangles
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dImpl::BuildEarcutTriangles(const XmUGrid& a_ugrid)
{
  Initialize(a_ugrid);

  int numCells = a_ugrid.GetCellCount();
  VecInt cellPoints;
  for (int cellIdx = 0; cellIdx < numCells; ++cellIdx)
  {
    a_ugrid.GetCellPoints(cellIdx, cellPoints);
    iBuildEarcutTriangles(*this, cellIdx, cellPoints);
  }
} // XmUGridTriangles2dImpl::BuildEarcutTriangles
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dImpl::SetCellActivity(const DynBitset& a_cellActivity)
{
  if (a_cellActivity.empty())
  {
    DynBitset emptyActivity;
    GetTriSearch()->SetTriActivity(emptyActivity);
    return;
  }

  DynBitset triangleActivity;
  int numTriangles = (int)m_triangleToCellIdx.size();
  triangleActivity.resize(numTriangles);
  for (size_t triangleIdx = 0; triangleIdx < numTriangles; ++triangleIdx)
  {
    int cellIdx = m_triangleToCellIdx[triangleIdx];
    triangleActivity[triangleIdx] = cellIdx >= a_cellActivity.size() || a_cellActivity[cellIdx];
  }
  GetTriSearch()->SetTriActivity(triangleActivity);
} // XmUGridTriangles2dImpl::SetCellActivity
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const VecPt3d& XmUGridTriangles2dImpl::GetPoints() const
{
  return *m_points;
} // XmUGridTriangles2dImpl::GetPoints
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
const VecInt& XmUGridTriangles2dImpl::GetTriangles() const
{
  return *m_triangles;
} // XmUGridTriangles2dImpl::GetTriangles
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
BSHP<VecPt3d> XmUGridTriangles2dImpl::GetPointsPtr()
{
  return m_points;
} // XmUGridTriangles2dImpl::GetPointsPtr
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
BSHP<VecInt> XmUGridTriangles2dImpl::GetTrianglesPtr()
{
  return m_triangles;
} // XmUGridTriangles2dImpl::GetTrianglesPtr
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGridTriangles2dImpl::AddCellCentroid(int a_cellIdx, const Pt3d& a_point)
{
  int centroidIdx = (int)m_points->size();
  m_points->push_back(a_point);
  m_centroidIdxs[a_cellIdx] = centroidIdx;
  return centroidIdx;
} // XmUGridTriangles2dImpl::AddCellCentroid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dImpl::AddCellTriangle(int a_cellIdx, int a_idx1, int a_idx2, int a_idx3)
{
  m_triangles->push_back(a_idx1);
  m_triangles->push_back(a_idx2);
  m_triangles->push_back(a_idx3);
  m_triangleToCellIdx.push_back(a_cellIdx);
} // XmUGridTriangles2dImpl::AddCellTriangle
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGridTriangles2dImpl::GetCellCentroid(int a_cellIdx) const
{
  int pointIdx = -1;
  if (a_cellIdx >= 0 && a_cellIdx < (int)m_centroidIdxs.size())
    pointIdx = m_centroidIdxs[a_cellIdx];
  return pointIdx;
} // XmUGridTriangles2dImpl::GetCellCentroid
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
int XmUGridTriangles2dImpl::GetIntersectedCell(const Pt3d& a_point,
                                               VecInt& a_idxs,
                                               VecDbl& a_weights)
{
  if (!m_triSearch)
    GetTriSearch();
  int cellIdx = -1;
  int triangleLocation;
  if (m_triSearch->InterpWeightsTriangleIdx(a_point, triangleLocation, a_idxs, a_weights))
  {
    int triangleIdx = triangleLocation / 3;
    cellIdx = m_triangleToCellIdx[triangleIdx];
  }
  return cellIdx;
} // XmUGridTriangles2dImpl::GetIntersectedCell
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dImpl::Initialize(const XmUGrid& a_ugrid)
{
  *m_points = a_ugrid.GetLocations();
  m_triangles->clear();
  m_centroidIdxs.assign(a_ugrid.GetCellCount(), -1);
  m_triangleToCellIdx.clear();
  m_triSearch.reset();
} // XmUGridTriangles2dImpl::Initialize
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
BSHP<GmTriSearch> XmUGridTriangles2dImpl::GetTriSearch()
{
  if (!m_triSearch)
  {
    m_triSearch = GmTriSearch::New();
    m_triSearch->TrisToSearch(m_points, m_triangles);
  }
  return m_triSearch;
} // XmUGridTriangles2dImpl::GetTriSearch

} // namespace

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
BSHP<XmUGridTriangles2d> XmUGridTriangles2d::New()
{
  BSHP<XmUGridTriangles2d> triangles(new XmUGridTriangles2dImpl);
  return triangles;
} // XmUGridTriangles2d::New
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridTriangles2d::XmUGridTriangles2d()
{
} // XmUGridTriangles2d::XmUGridTriangles2d
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
XmUGridTriangles2d::~XmUGridTriangles2d()
{
} // XmUGridTriangles2d::XmUGridTriangles2d

} // namespace xms

#ifdef CXX_TEST
//------------------------------------------------------------------------------
// Unit Tests
//------------------------------------------------------------------------------
using namespace xms;
#include <xmsextractor/ugrid/XmUGridTriangles2d.t.h>

#include <xmscore/testing/TestTools.h>

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dUnitTests::testBuildCentroidTrianglesOnTriangle()
{
  VecPt3d points = {{0, 0, 0}, {1, 0, 0}, {0.5, 1, 0}};
  VecInt cells = {XMU_TRIANGLE, 3, 0, 1, 2};
  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);
  TS_REQUIRE_NOT_NULL(ugrid);
  XmUGridTriangles2dImpl triangles;

  // test building cetroid in triangle cell
  triangles.BuildTriangles(*ugrid, true);

  VecPt3d triPointsOut = triangles.GetPoints();
  VecPt3d triPointsExpected = {{0, 0, 0}, {1, 0, 0}, {0.5, 1, 0}, {0.5, 1 / 3.0, 0}};
  TS_ASSERT_EQUALS(triPointsExpected, triPointsOut);

  VecInt trianglesOut = triangles.GetTriangles();
  VecInt trianglesExpected = {0, 1, 3, 1, 2, 3, 2, 0, 3};
  TS_ASSERT_EQUALS(trianglesExpected, trianglesOut);

  TS_ASSERT_EQUALS(3, triangles.GetCellCentroid(0));

  VecInt idxs;
  VecDbl weights;
  // point intersecting last triangle
  int cellIdx = triangles.GetIntersectedCell(Pt3d(0.25, 0.25, 0), idxs, weights);
  TS_ASSERT_EQUALS(0, cellIdx);
  VecInt idxsExpected = {2, 0, 3};
  TS_ASSERT_EQUALS(idxsExpected, idxs);
  VecDbl weightsExpected = {0.125, 0.5, 0.375};
  TS_ASSERT_EQUALS(weightsExpected, weights);
  cellIdx = triangles.GetIntersectedCell(Pt3d(0.75, 0.25, 0), idxs, weights);
  TS_ASSERT_EQUALS(0, cellIdx);
  idxsExpected = {1, 2, 3};
  TS_ASSERT_EQUALS(idxsExpected, idxs);
  weightsExpected = {0.5, 0.125, 0.375};
  TS_ASSERT_EQUALS(weightsExpected, weights);
  cellIdx = triangles.GetIntersectedCell(Pt3d(0.5, 0.25, 0), idxs, weights);
  TS_ASSERT_EQUALS(0, cellIdx);
  idxsExpected = {0, 1, 3};
  TS_ASSERT_EQUALS(idxsExpected, idxs);
  weightsExpected = {0.125, 0.125, 0.75};
  TS_ASSERT_EQUALS(weightsExpected, weights);

  // test without building centroid in triangle cell
  triangles.BuildTriangles(*ugrid, false);

  triPointsOut = triangles.GetPoints();
  TS_ASSERT_EQUALS(points, triPointsOut);

  trianglesOut = triangles.GetTriangles();
  trianglesExpected = {0, 1, 2};
  TS_ASSERT_EQUALS(trianglesExpected, trianglesOut);

  TS_ASSERT_EQUALS(-1, triangles.GetCellCentroid(0));

  cellIdx = triangles.GetIntersectedCell(Pt3d(0.5, 0.25, 0), idxs, weights);
  TS_ASSERT_EQUALS(0, cellIdx);
  idxsExpected = {0, 1, 2};
  TS_ASSERT_EQUALS(idxsExpected, idxs);
  weightsExpected = {0.375, 0.375, 0.25};
  TS_ASSERT_EQUALS(weightsExpected, weights);

} // XmUGridTriangles2dUnitTests::testBuildCentroidTrianglesOnTriangle
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dUnitTests::testBuildCentroidTrianglesOnQuad()
{
  VecPt3d points = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}};
  VecInt cells = {XMU_QUAD, 4, 0, 1, 2, 3};
  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);
  TS_REQUIRE_NOT_NULL(ugrid);
  XmUGridTriangles2dImpl triangles;
  triangles.BuildTriangles(*ugrid, true);
  VecPt3d triPointsOut = triangles.GetPoints();
  VecPt3d triPointsExpected = {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {0.5, 0.5, 0}};
  TS_ASSERT_EQUALS(triPointsExpected, triPointsOut);

  VecInt trianglesOut = triangles.GetTriangles();
  VecInt trianglesExpected = {0, 1, 4, 1, 2, 4, 2, 3, 4, 3, 0, 4};
  TS_ASSERT_EQUALS(trianglesExpected, trianglesOut);

  TS_ASSERT_EQUALS(4, triangles.GetCellCentroid(0));
} // XmUGridTriangles2dUnitTests::testBuildCentroidTrianglesOnQuad
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dUnitTests::testBuildCentroidTriangles2dCellTypes()
{
  // clang-format off
  VecPt3d points = {
    {0, 0, 0},   // 0
    {10, 0, 0},  // 1
    {20, 0, 0},  // 2
    {30, 0, 0},  // 3
    {40, 0, 0},  // 4
    {0, 10, 0},  // 5
    {10, 10, 0}, // 6
    {20, 10, 0}, // 7
    {40, 10, 0}, // 8
    {0, 20, 0},  // 9
    {10, 20, 0}, // 10
    {20, 20, 0}, // 11
    {30, 20, 0}, // 12
    {40, 20, 0}  // 13
  };

  // Cell type (5), number of points (3), point numbers, counterclockwise
  std::vector<int> cells = {(int)XMU_QUAD, 4, 0, 1, 6, 5,             // 0
                            (int)XMU_PIXEL, 4, 1, 2, 6, 7,            // 1
                            (int)XMU_TRIANGLE, 3, 2, 3, 7,            // 2
                            (int)XMU_POLYGON, 6, 3, 4, 8, 13, 12, 7}; // 3
  // clang-format on

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);
  TS_REQUIRE_NOT_NULL(ugrid);
  XmUGridTriangles2dImpl triangles;
  triangles.BuildTriangles(*ugrid, true);

  VecPt3d triPointsOut = triangles.GetPoints();
  // clang-format off
  VecPt3d triPointsExpected = {
    {0.0, 0.0, 0.0},
    {10.0, 0.0, 0.0},
    {20.0, 0.0, 0.0},
    {30.0, 0.0, 0.0},
    {40.0, 0.0, 0.0},
    {0.0, 10.0, 0.0},
    {10.0, 10.0, 0.0},
    {20.0, 10.0, 0.0},
    {40.0, 10.0, 0.0},
    {0.0, 20.0, 0.0},
    {10.0, 20.0, 0.0},
    {20.0, 20.0, 0.0},
    {30.0, 20.0, 0.0},
    {40.0, 20.0, 0.0},

    // added centroids
    {5.0, 5.0, 0.0},                               // 14 - quad centroid
    {15.0, 5.0, 0.0},                              // 15 - pixel centroid
    {23.333333333333332, 3.3333333333333335, 0.0}, // 16 - triangle centroid
    {33.333333333333336, 10.0, 0.0}                // 17 - polygon centroid
  };
  // clang-format on
  TS_ASSERT_EQUALS(triPointsExpected, triPointsOut);

  VecInt trianglesOut = triangles.GetTriangles();
  // clang-format off
  VecInt trianglesExpected = {
    0, 1, 14, 1, 6, 14, 6, 5, 14, 5, 0, 14, // quad
    1, 2, 15, 2, 7, 15, 7, 6, 15, 6, 1, 15, // pixel
    2, 3, 16, 3, 7, 16, 7, 2, 16,           // triangle
    3, 4, 17, 4, 8, 17, 8, 13, 17, 13, 12, 17, 12, 7, 17, 7, 3, 17}; // polygon
  // clang-format on
  TS_ASSERT_EQUALS(trianglesExpected, trianglesOut);
} // XmUGridTriangles2dUnitTests::testBuildCentroidTriangles2dCellTypes
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dUnitTests::testBuildEarcutTriangles()
{
  {
    // V shaped cell
    // clang-format off
    VecPt3d points = {
      {0, 0, 0},   // 0
      {10, 0, 0},  // 1
      {1, 1, 0},   // 2
      {0, 10, 0},  // 3
      {10, 10, 0}  // 4
    };

    std::vector<int> cells = {
      XMU_QUAD, 4, 0, 1, 2, 3,
    };
    // clang-format on

    std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);
    XmUGridTriangles2dImpl ugridTris;

    ugridTris.BuildTriangles(*ugrid, true);

    VecPt3d triPointsOut = ugridTris.GetPoints();
    VecPt3d triPointsExpected = points;
    TS_ASSERT_EQUALS(triPointsExpected, triPointsOut);

    VecInt trianglesOut = ugridTris.GetTriangles();
    VecInt trianglesExpected = {2, 3, 0, 0, 1, 2};
    TS_ASSERT_EQUALS(trianglesExpected, trianglesOut);
  }

  {
    // house shaped cell
    // clang-format off
    VecPt3d points = {
      {0, 0, 0},   // 0
      {10, 0, 0},  // 1
      {10, 10, 0}, // 2
      {5, 11, 0},  // 3
      {0, 10, 0},  // 4
    };

    std::vector<int> cells = {(int)XMU_POLYGON, 5, 0, 1, 2, 3, 4};
    // clang-format on

    std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);
    XmUGridTriangles2dImpl ugridTris;
    ugridTris.BuildEarcutTriangles(*ugrid);

    VecPt3d triPointsOut = ugridTris.GetPoints();
    VecPt3d triPointsExpected = points;
    TS_ASSERT_EQUALS(triPointsExpected, triPointsOut);

    VecInt trianglesOut = ugridTris.GetTriangles();
    VecInt trianglesExpected = {4, 0, 1, 1, 2, 3, 1, 3, 4};
    TS_ASSERT_EQUALS(trianglesExpected, trianglesOut);
  }
} // XmUGridTriangles2dUnitTests::testBuildEarcutTriangles
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void XmUGridTriangles2dUnitTests::testBuildCentroidAndEarcutTriangles()
{
  // 7----6---------3----2
  // |    |         |    |
  // |    |         |    |
  // |    |         |    |
  // |    |   (8)   |    |
  // |    |         |    |
  // |    |         |    |
  // |    |         |    |
  // |    5---------4    |
  // |                   |
  // |                   |
  // |                   |
  // 0-------------------1
  // clang-format off
  VecPt3d points = {
    { 0,  0, 0}, // 0
    {15,  0, 0}, // 1
    {15, 15, 0}, // 2
    {10, 15, 0}, // 3
    {10,  5, 0}, // 4
    { 5,  5, 0}, // 5
    { 5, 15, 0}, // 6
    { 0, 15, 0}  // 7
  };

  // Cell type (5), number of points (3), point numbers, counterclockwise
  std::vector<int> cells = {
    XMU_POLYGON, 8, 0, 1, 2, 3, 4, 5, 6, 7,
    XMU_QUAD, 4, 3, 6, 5, 4
  };
  // clang-format on

  std::shared_ptr<XmUGrid> ugrid = XmUGrid::New(points, cells);
  XmUGridTriangles2dImpl ugridTris;

  ugridTris.BuildTriangles(*ugrid, true);

  VecPt3d triPointsOut = ugridTris.GetPoints();
  VecPt3d triPointsExpected = points;
  triPointsExpected.push_back({7.5, 10, 0});
  TS_ASSERT_EQUALS(triPointsExpected, triPointsOut);

  VecInt trianglesOut = ugridTris.GetTriangles();
  VecInt trianglesExpected = {
    // clang-format off
    2, 3, 4,   5, 6, 7,   1, 2, 4,   0, 1, 4,   0, 4, 5,   0, 5, 7, // earcut
    3, 6, 8,   6, 5, 8,   5, 4, 8,   4, 3, 8 // centroid
    // clang-format off
  };
  TS_ASSERT_EQUALS(trianglesExpected, trianglesOut);
} // XmUGridTriangles2dUnitTests::testBuildCentroidAndEarcutTriangles

#endif