vtkPolygon.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#ifndef vtkPolygon_h
#define vtkPolygon_h

#include "vtkCell.h"
#include "vtkCellStatus.h"            // For return type
#include "vtkCommonDataModelModule.h" // For export macro

#include <cmath>   // For std::sqrt in the quality ear clip
#include <utility> // For std::swap in the quality ear clip
#include <vector>  // For ear-clip scratch buffers

VTK_ABI_NAMESPACE_BEGIN
class vtkDoubleArray;
class vtkIdTypeArray;
class vtkLine;
class vtkPoints;
class vtkQuad;
class vtkTriangle;
class vtkPriorityQueue;
class vtkIncrementalPointLocator;

class VTKCOMMONDATAMODEL_EXPORT vtkPolygon : public vtkCell
{
public:
  static vtkPolygon* New();
  vtkTypeMacro(vtkPolygon, vtkCell);
  void PrintSelf(ostream& os, vtkIndent indent) override;


  int GetCellType() override { return VTK_POLYGON; }
  int GetCellDimension() override { return 2; }
  int GetNumberOfEdges() override { return this->GetNumberOfPoints(); }
  int GetNumberOfFaces() override { return 0; }
  vtkCell* GetEdge(int edgeId) override;
  vtkCell* GetFace(int) override { return nullptr; }
  int CellBoundary(int subId, const double pcoords[3], vtkIdList* pts) override;
  void Contour(double value, vtkDataArray* cellScalars, vtkIncrementalPointLocator* locator,
    vtkCellArray* verts, vtkCellArray* lines, vtkCellArray* polys, vtkPointData* inPd,
    vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId, vtkCellData* outCd) override;
  void Clip(double value, vtkDataArray* cellScalars, vtkIncrementalPointLocator* locator,
    vtkCellArray* tris, vtkPointData* inPd, vtkPointData* outPd, vtkCellData* inCd,
    vtkIdType cellId, vtkCellData* outCd, int insideOut) override;
  int EvaluatePosition(const double x[3], double closestPoint[3], int& subId, double pcoords[3],
    double& dist2, double weights[]) override;
  void EvaluateLocation(int& subId, const double pcoords[3], double x[3], double* weights) override;
  int IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t, double x[3],
    double pcoords[3], int& subId) override;
  int TriangulateLocalIds(int index, vtkIdList* ptIds) override;
  void Derivatives(
    int subId, const double pcoords[3], const double* values, int dim, double* derivs) override;
  int IsPrimaryCell() VTK_FUTURE_CONST override { return 0; }

  double ComputeArea();

  void InterpolateFunctions(const double x[3], double* sf) override;


  static vtkCellStatus ComputeNormal(vtkPoints* p, int numPts, const vtkIdType* pts, double n[3]);
  static vtkCellStatus ComputeNormal(vtkPoints* p, double n[3]);
  static vtkCellStatus ComputeNormal(vtkIdTypeArray* ids, vtkPoints* pts, double n[3]);

  static vtkCellStatus ComputeNormal(int numPts, double* pts, double n[3]);

  bool IsConvex();


  static bool IsConvex(vtkPoints* p, int numPts, const vtkIdType* pts);
  static bool IsConvex(vtkIdTypeArray* ids, vtkPoints* p);
  static bool IsConvex(vtkPoints* p);


  static vtkCellStatus ComputeCentroid(
    vtkPoints* p, int numPts, const vtkIdType* pts, double centroid[3], double tolerance);
  static bool ComputeCentroid(vtkPoints* p, int numPts, const vtkIdType* pts, double centroid[3]);
  static bool ComputeCentroid(vtkIdTypeArray* ids, vtkPoints* pts, double centroid[3]);

  static bool ComputeInteriorCircle(
    vtkPoints* p, int numPts, const vtkIdType* ids, double center[3], double& radius2);

  static double ComputeArea(vtkPoints* p, vtkIdType numPts, const vtkIdType* pts, double normal[3]);

  int ParameterizePolygon(
    double p0[3], double p10[3], double& l10, double p20[3], double& l20, double n[3]);

  static int PointInPolygon(double x[3], int numPts, double* pts, double bounds[6], double n[3]);

  // Needed to remove warning "member function does not override any
  // base class virtual member function"
  int Triangulate(int index, vtkIdList* ptIds, vtkPoints* pts) override
  {
    return vtkCell::Triangulate(index, ptIds, pts);
  }

  int NonDegenerateTriangulate(vtkIdList* outTris);

  int BoundedTriangulate(vtkIdList* outTris, double tol);

  static double DistanceToPolygon(
    double x[3], int numPts, double* pts, double bounds[6], double closest[3]);

  static int IntersectPolygonWithPolygon(int npts, double* pts, double bounds[6], int npts2,
    double* pts2, double bounds2[6], double tol, double x[3]);

  static int IntersectConvex2DCells(
    vtkCell* cell1, vtkCell* cell2, double tol, double p0[3], double p1[3]);


  vtkGetMacro(UseMVCInterpolation, bool);
  vtkSetMacro(UseMVCInterpolation, bool);


  vtkSetClampMacro(Tolerance, double, 0.0, 1.0);
  vtkGetMacro(Tolerance, double);

protected:
  vtkPolygon();
  ~vtkPolygon() override = default;

  // Compute the interpolation functions using Mean Value Coordinate.
  void InterpolateFunctionsUsingMVC(const double x[3], double* weights);

  // variables used by instances of this class
  double Tolerance;        // Intersection tolerance set by public API
  double Tol;              // Internal tolerance set by ComputeBounds()
  void ComputeTolerance(); // Compute the internal tolerance Tol

  int SuccessfulTriangulation;     // Stops recursive triangulation if necessary
  vtkSmartPointer<vtkIdList> Tris; // Output triangulation placed here

  // These are used for internal computation.
  vtkSmartPointer<vtkTriangle> Triangle;
  vtkSmartPointer<vtkQuad> Quad;
  vtkSmartPointer<vtkDoubleArray> TriScalars;
  vtkSmartPointer<vtkLine> Line;
  vtkSmartPointer<vtkPriorityQueue> EarClipQueue; // reused ear-clip removal queue

  // Parameter indicating whether to use Mean Value Coordinate algorithm
  // for interpolation. The parameter is false by default.
  bool UseMVCInterpolation;

  // Helper methods for triangulation------------------------------
  // Made public for external access
public:
  // Ear cut triangulation options. The order in which vertices are
  // removed are controlled by different measures. Changing this can
  // make subtle differences in some cases. Historically the
  // PERIMETER2_TO_AREA_RATIO has been used.
  enum EarCutMeasureTypes
  {
    PERIMETER2_TO_AREA_RATIO = 0,
    DOT_PRODUCT = 1,
    BEST_QUALITY = 2
  };


  int EarCutTriangulation(int measure = PERIMETER2_TO_AREA_RATIO);
  int EarCutTriangulation(vtkIdList* outTris, int measure = PERIMETER2_TO_AREA_RATIO);


  int UnbiasedEarCutTriangulation(int seed, int measure = PERIMETER2_TO_AREA_RATIO);
  int UnbiasedEarCutTriangulation(
    int seed, vtkIdList* outTris, int measure = PERIMETER2_TO_AREA_RATIO);


  template <typename PointsRange, typename CellIter>
  static int CompactPolygonRing(
    const PointsRange& points, CellIter cell, int cellSize, std::vector<int>& ring);
  template <typename PointsRange, typename CellIter>
  static vtkIdType EarClipTriangleCount(
    const PointsRange& points, CellIter cell, int cellSize, std::vector<int>& ring);
  template <typename PointsRange, typename CellIter, typename EmitFn>
  static void EarClipPolygon3D(const PointsRange& points, CellIter cell, int cellSize,
    std::vector<int>& prevBuf, std::vector<int>& nextBuf, std::vector<int>& ring, EmitFn&& emit);

private:
  vtkPolygon(const vtkPolygon&) = delete;
  void operator=(const vtkPolygon&) = delete;
};

//------------------------------------------------------------------------------
template <typename PointsRange, typename CellIter>
int vtkPolygon::CompactPolygonRing(
  const PointsRange& points, CellIter cell, int cellSize, std::vector<int>& ring)
{
  // Build the ring of distinct vertices, dropping consecutive coincident
  // vertices (degenerate edges). Banded contouring and clipping routinely emit
  // polygons with coincident consecutive vertices when a scalar lands exactly
  // on a clip value; the historical fan triangulation tolerated these by
  // skipping zero-area triangles, and the ear clip must do the same or it
  // computes a bogus normal / ear test at the duplicated vertex.
  ring.clear();
  for (int i = 0; i < cellSize; ++i)
  {
    if (!ring.empty())
    {
      auto pPrev = points[cell[ring.back()]];
      auto pCur = points[cell[i]];
      if (pPrev[0] == pCur[0] && pPrev[1] == pCur[1] && pPrev[2] == pCur[2])
      {
        continue; // coincident with previous kept vertex
      }
    }
    ring.push_back(i);
  }
  // Drop the last vertex if it coincides with the first (wrap-around duplicate).
  if (ring.size() >= 2)
  {
    auto pFirst = points[cell[ring.front()]];
    auto pLast = points[cell[ring.back()]];
    if (pFirst[0] == pLast[0] && pFirst[1] == pLast[1] && pFirst[2] == pLast[2])
    {
      ring.pop_back();
    }
  }
  return static_cast<int>(ring.size());
}

//------------------------------------------------------------------------------
template <typename PointsRange, typename CellIter>
vtkIdType vtkPolygon::EarClipTriangleCount(
  const PointsRange& points, CellIter cell, int cellSize, std::vector<int>& ring)
{
  const int m = vtkPolygon::CompactPolygonRing(points, cell, cellSize, ring);
  return m >= 3 ? static_cast<vtkIdType>(m - 2) : 0;
}

//------------------------------------------------------------------------------
template <typename PointsRange, typename CellIter, typename EmitFn>
void vtkPolygon::EarClipPolygon3D(const PointsRange& points, CellIter cell, int cellSize,
  std::vector<int>& prevBuf, std::vector<int>& nextBuf, std::vector<int>& ring, EmitFn&& emit)
{
  const int m = vtkPolygon::CompactPolygonRing(points, cell, cellSize, ring);
  if (m < 3)
  {
    return; // fully degenerate (collapses to a point or segment): no triangles
  }

  // ring[k] is the polygon-local index of the k-th distinct vertex. Work in
  // compacted space 0..m-1; map back through ring[] when emitting and looking
  // up coordinates.
  auto P = [&](int k) { return points[cell[ring[k]]]; };

  // Compute polygon normal via Newell's method over the compacted ring.
  double normal[3] = { 0.0, 0.0, 0.0 };
  {
    auto pLast = P(m - 1);
    double xp = pLast[0], yp = pLast[1], zp = pLast[2];
    for (int i = 0; i < m; ++i)
    {
      auto pi = P(i);
      double x = pi[0], y = pi[1], z = pi[2];
      normal[0] += (yp - y) * (zp + z);
      normal[1] += (zp - z) * (xp + x);
      normal[2] += (xp - x) * (yp + y);
      xp = x;
      yp = y;
      zp = z;
    }
  }

  // Emit a triangle given compacted-ring indices, computing the polygon-
  // boundary edge mask from compacted-ring adjacency, then mapping to local.
  auto emitRing = [&](int a, int b, int c)
  {
    auto isBoundary = [m](int x, int y) -> int
    {
      const int d = (y - x + m) % m;
      return (d == 1 || d == m - 1) ? 1 : 0;
    };
    int mask = isBoundary(a, b);
    mask |= isBoundary(b, c) << 1;
    mask |= isBoundary(c, a) << 2;
    emit(ring[a], ring[b], ring[c], mask);
  };

  const double normLen2 = normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2];
  if (normLen2 == 0.0)
  {
    // Zero-area (collinear) polygon. No meaningful triangulation; fall back to
    // a fan over the compacted ring so the emit count matches m - 2.
    for (int i = 1; i < m - 1; ++i)
    {
      emitRing(0, i, i + 1);
    }
    return;
  }

  // Doubly-linked circular list over the compacted ring 0..m-1.
  prevBuf.resize(m);
  nextBuf.resize(m);
  for (int i = 0; i < m; ++i)
  {
    prevBuf[i] = (i + m - 1) % m;
    nextBuf[i] = (i + 1) % m;
  }

  auto isEar = [&](int b) -> bool
  {
    const int a = prevBuf[b];
    const int c = nextBuf[b];
    auto pa = P(a);
    auto pb = P(b);
    auto pc = P(c);
    const double ax = pa[0], ay = pa[1], az = pa[2];
    const double bx = pb[0], by = pb[1], bz = pb[2];
    const double cx = pc[0], cy = pc[1], cz = pc[2];

    // Convexity at b: ((b - a) x (c - b)) . normal > 0
    const double e1x = bx - ax, e1y = by - ay, e1z = bz - az;
    const double e2x = cx - bx, e2y = cy - by, e2z = cz - bz;
    const double crx = e1y * e2z - e1z * e2y;
    const double cry = e1z * e2x - e1x * e2z;
    const double crz = e1x * e2y - e1y * e2x;
    if (crx * normal[0] + cry * normal[1] + crz * normal[2] <= 0.0)
    {
      return false; // reflex or zero-area
    }

    // No other vertex in the remaining polygon may lie strictly inside
    // triangle (a, b, c). Use the polygon normal as the projection axis for
    // a same-side test against each triangle edge.
    for (int q = nextBuf[c]; q != a; q = nextBuf[q])
    {
      auto pq = P(q);
      const double qx = pq[0], qy = pq[1], qz = pq[2];
      auto sideOf = [&](double sx, double sy, double sz, double ex, double ey, double ez) -> double
      {
        const double dx = ex - sx, dy = ey - sy, dz = ez - sz;
        const double rx = qx - sx, ry = qy - sy, rz = qz - sz;
        const double tcrx = dy * rz - dz * ry;
        const double tcry = dz * rx - dx * rz;
        const double tcrz = dx * ry - dy * rx;
        return tcrx * normal[0] + tcry * normal[1] + tcrz * normal[2];
      };
      const double s1 = sideOf(ax, ay, az, bx, by, bz);
      const double s2 = sideOf(bx, by, bz, cx, cy, cz);
      const double s3 = sideOf(cx, cy, cz, ax, ay, az);
      if (s1 > 0.0 && s2 > 0.0 && s3 > 0.0)
      {
        return false;
      }
    }
    return true;
  };

  // Clip ears walking forward from vertex 1, keeping the clipped vertex's
  // successor as the next candidate. For a convex polygon this clips vertices
  // 1, 2, 3, ... in order while vertex 0 remains the common apex, reproducing
  // exactly the fan-from-vertex-0 triangulation VTK has always emitted (and
  // that existing image baselines were generated with). Triangulation is not
  // invariant under interior diagonals: per-triangle texture-coordinate and
  // Gouraud interpolation depend on the diagonals chosen, so a convex polygon
  // must keep producing the fan. For non-convex polygons reflex vertices are
  // skipped and the walk still finds valid ears.
  int remaining = m;
  int current = 1 % m;
  int safetyBudget = 2 * m;
  while (remaining > 3 && safetyBudget > 0)
  {
    if (isEar(current))
    {
      const int a = prevBuf[current];
      const int c = nextBuf[current];
      emitRing(a, current, c);
      nextBuf[a] = c;
      prevBuf[c] = a;
      current = c;
      --remaining;
      safetyBudget = 2 * m;
    }
    else
    {
      current = nextBuf[current];
      --safetyBudget;
    }
  }

  if (remaining > 3)
  {
    // Ear-clip stalled. Emit remaining vertices as a fan from `current` to
    // satisfy the emit-pass triangle count. Result will be visually wrong for
    // a non-convex remainder, but degenerate input has no correct triangulation.
    const int head = current;
    int v = nextBuf[head];
    int vn = nextBuf[v];
    while (vn != head)
    {
      emitRing(head, v, vn);
      v = vn;
      vn = nextBuf[vn];
    }
    return;
  }

  // Final triangle: the three vertices still in the list. Emit it anchored at
  // prev(current) so that for a convex polygon (where current has advanced to
  // the second-to-last vertex with vertex 0 still the apex) the triangle is
  // (0, m-2, m-1) - the same vertex order and boundary-edge mask the fan emits.
  const int b = current;
  const int a = prevBuf[b];
  const int c = nextBuf[b];
  emitRing(a, b, c);
}

VTK_ABI_NAMESPACE_END
#endif