Как перенести алгоритм триангуляции 2D-полигонов в 3D?

Question

Я пришел по коду в этой ссылке: 2D триангуляция Джона Рэтклиффа .Код представляет алгоритм триангуляции в C ++ для 2D полигонов.Я пытался перенести его на 3D полигоны, но безуспешно.Как я могу изменить этот код, чтобы он использовал трехмерные плоские многоугольники?Я знаю, что могу использовать проекцию, но есть ли способ сделать 3D-версию кода, указанного выше?Благодарю.NB: Ниже моя попытка портирования.

double Triangulate::Area(const std::vector<int> &poly)
{

  int n = poly.size();

  if (n < 3) // not a plane - no area
        return 0;

  std::vector<double> total = {0, 0, 0};
  int vi1, vi2;
  std::vector<double> prod;
  std::vector<double> v1, v2, v3;
  double* point = NULL;
  for(int i = 0; i < n; i++){
        v1.clear();
        v2.clear();
        vi1 = poly[i];
        if (i == n-1)
            vi2 = poly[0];
        else
            vi2 = poly[i+1];
        point = input.GetPoint(vi1);
        v1.insert(v1.begin(), point, point+3);

        point = input.GetPoint(vi2);
        v2.insert(v2.begin(), point, point+3);

        prod = Cross(v1, v2);
        total[0] += prod[0];
        total[1] += prod[1];
        total[2] += prod[2];
    }
    v1.clear();
    v2.clear();
    v3.clear();
    point = input.GetPoint(poly[0]);
    v1.insert(v1.begin(), point, point+3);

    point = input.GetPoint(poly[1]);
    v2.insert(v2.begin(), point, point+3);

    point = input.GetPoint(poly[2]);
    v3.insert(v3.begin(), point, point+3);

    double result = Dot(total, UnitNormal(v1, v2, v3));
    return result/2;
}

bool Triangulate::SameSide(double P1x,  double P1y, double P1z, 
                    double P2x,  double P2y, double P2z,
                    double Ax,  double Ay, double Az,
                    double Bx,  double By, double Bz){
    double ABx, ABy, ABz, AP1x, AP1y, AP1z,
        AP2x, AP2y, AP2z;

    double cp1x, cp1y, cp1z, cp2x, cp2y, cp2z;

    double dot;
    std::vector<double> AB = {Bx - Ax, By - Ay, Bz - Az};

    std::vector<double> AP1 = {P1x - Ax,P1y - Ay,  P1z - Az};

    std::vector<double> AP2 = {P2x - Ax, P2y - Ay, P2z - Az};

    std::vector<double> cp1 = Cross(AB, AP1);
    std::vector<double> cp2 = Cross(AB, AP2);

    dot = Dot(cp1, cp2);

   return (dot >= 0);
}

   /*
     InsideTriangle decides if a point P is Inside of the triangle
     defined by A, B, C.
   */
bool Triangulate::InsideTriangle(double Ax, double Ay, double Az,
                      double Bx, double By, double Bz,
                      double Cx, double Cy, double Cz,
                      double Px, double Py, double Pz)

{
    return SameSide(Px, Py, Pz, Ax, Ay, Az, Bx, By, Bz, Cx, Cy, Cz) 
      || SameSide(Px, Py, Pz, Bx, By, Bz, Ax, Ay, Az, Cx, Cy, Cz) 
      || SameSide(Px, Py, Pz, Cx, Cy, Cz, Ax, Ay, Az, Bx, By, Bz) ;
}

bool Triangulate::Snip(const std::vector<int> &contour,
   int u,int v,int w,int n,int *V)
{
  int p;
  double Ax, Ay,Az, Bx, By,Bz, Cx, Cy, Cz, Px, Py, Pz;
  double BAx, BAy, BAz, CAx, CAy, CAz;

  Ax = input.GetPoint(V[u])[X];
  Ay = input.GetPoint(V[u])[Y];
  Az = input.GetPoint(V[u])[Z];

  Bx = input.GetPoint(V[v])[X];
  By = input.GetPoint(V[v])[Y];
  Bz = input.GetPoint(V[v])[Z];

  Cx = input.GetPoint(V[w])[X];
  Cy = input.GetPoint(V[w])[Y];
  Cz = input.GetPoint(V[w])[Z];

  BAx = Bx-Ax;
  BAy = By - Ay;
  BAz = Bz - Az;

  CAx = Cx-Ax;
  CAy = Cy - Ay;
  CAz = Cz - Az;

  if ( EPSILON > ( (BAy*CAz-BAz*CAy)*(BAy*CAz-BAz*CAy)+ 
   (BAx*CAz-BAz*CAx)*(BAx*CAz-BAz*CAx) +
   (BAx*CAy-BAy*CAx)*(BAx*CAy-BAy*CAx))  ) return false;

  for (p=0;p<n;p++)
  {
    if( (p == u) || (p == v) || (p == w) ) continue;
    Px = input.GetPoint(V[p])[X];
    Py = input.GetPoint(V[p])[Y];
    Pz = input.GetPoint(V[p])[Z];

    if (InsideTriangle(Ax,Ay,Az,Bx,By,Bz,Cx,Cy,Cz,Px,Py,Pz)) return false;
  }

  return true;
}

int Triangulate::Process(const std::vector<int> &contour,
std::vector<std::vector<int>> &result)
{
  /* allocate and initialize list of Vertices in polygon */

  int n = contour.size();
  if ( n < 3 ) {
    LOG_VAL("Cannot happen: contour is a line", "")
   return 0;
  }

  int *V = new int[n];

  /* we want a counter-clockwise polygon in V */

  if ( 0.0f < Area(contour) )
    for (int v=0; v<n; v++) V[v] = v;
  else
    for(int v=0; v<n; v++) V[v] = (n-1)-v;

  int nv = n;

  /*  remove nv-2 Vertices, creating 1 triangle every time */
  int count = 2*nv;   /* error detection */
  std::vector<int> tri;
  int trinum = 0;

  for(int m=0, v=nv-1; nv>2; )
  {
    /* if we loop, it is probably a non-simple polygon */
    if (0 >= (count--))
    {
      //** Triangulate: ERROR - probable bad polygon!
      LOG_VAL("Bad polygon:", trinum)
      return trinum;
    }

    /* three consecutive vertices in current polygon, <u,v,w> */
    int u = v  ; if (nv <= u) u = 0;     /* previous */
    v = u+1; if (nv <= v) v = 0;     /* new v    */
    int w = v+1; if (nv <= w) w = 0;     /* next     */

    if ( Snip(contour,u,v,w,nv,V) )
    {
      int a,b,c,s,t;

      /* true names of the vertices */
      a = V[u]; b = V[v]; c = V[w];

      /* output Triangle */
      tri.clear();
      tri.push_back( contour[a] );
      tri.push_back( contour[b] );
      tri.push_back( contour[c] );
      result.push_back(tri);
      trinum++;

      m++;

      /* remove v from remaining polygon */
      for(s=v,t=v+1;t<nv;s++,t++) V[s] = V[t]; nv--;

      /* resest error detection counter */
      count = 2*nv;
    }
  }



  delete V;

  if(result.size() < 1){
      LOG_VAL("Result Is Empty:", "")
      exit(1);
  }

  return trinum;
}