0

I came over the code in this link: 2D Triangulation by John Ratcliff. The code presents a triangulation algorithm in C++ for 2D polygons. I've tried to port it to 3D polygons but unsuccessfully. How can I modify this code to make it use 3D planar polygons? I know I can use a projection, but is there a way to make a 3D version of the code linked above? Thanks. NB: Below is my porting attempt.

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;
}
N. Wells
  • 123
  • 10
  • As polygons are planar you can transfer your coordinates to 2D before calculation use "2D Triangulation by John Ratcliff" as it is and then transfer them back. So you do not even need to modify that code, just add pre and post processing. – Slava Oct 03 '18 at 14:31
  • I think 12 parameters for a function is a little bit too much. Why don't you use 4 `std::array` instead? – mch Oct 03 '18 at 14:40
  • @Slava, pre-projection comes with an additional computational cost, I cannot afford, as this piece is designed to process hundreds of thousands of polygons. – N. Wells Oct 03 '18 at 14:49
  • @mch, the code above is for testing. I intend to beautify it once it is ready and working. – N. Wells Oct 03 '18 at 14:50
  • If you do want to use algo from John Ratcliff I would recommend to rewrite it more generic way - do not pass coordinates as separate variables etc but only using Vector2D and operations on them. Then test it. After you make it work correct make that code templated or just replace `Vector2D` to `Vector3D`. Unfortunately code shown is not written in generic way. – Slava Oct 03 '18 at 15:22
  • For example you should not have parameters `double Ax` `double Ay` and `double Az` but `Vector2D A` and then replaced with `Vector3D A`. And you should not access separate coordinates of `A` you should use methods that work with 2 or 3 points. – Slava Oct 03 '18 at 15:24
  • @Slava, that's somehow, what I tried to do, especially identifying which vector function in 2D was used. For example, Ratcliff, doesn't explicitly use the cross product in finding areas. I had to google it out. – N. Wells Oct 03 '18 at 16:00
  • "*comes with an additional computational cost, I cannot afford*" – the cost of performing a simple projection 2D is likely to be significantly, if not *much*, lower than that of the triangulation algorithm itself, the former being only two or three vector operations per vertex. – meowgoesthedog Oct 03 '18 at 16:13
  • @meowgoesthedog, how do you find the plane of projection, for every polygon in a very complex mesh? How do you make sure the projection section (2D triangle area) is significant enough for the triangulation? – N. Wells Oct 04 '18 at 12:16
  • If your polygon is approximately planar then a good choice of normal would simply be the mean normal over all polygons, or even a small random sample of polygons. Again, this is also only a few vector operations per polygon. – meowgoesthedog Oct 04 '18 at 12:23
  • @meowgoesthedog, your proposition only works when the mesh is convex, right? Otherwise, I'll have to take as projection plane, the polygon plane for each face of the mesh. – N. Wells Oct 04 '18 at 12:27
  • 1
    Wasn't the very **premise** of your question that the polygons are planar? – meowgoesthedog Oct 04 '18 at 12:28
  • Nope, I didn't say so. I said that each polygon seats on its own plane. – N. Wells Oct 04 '18 at 12:30
  • "*How can I modify this code to make it use 3D planar polygons?*" Also, if each polygon seats on its own plane, why would it be difficult to obtain said plane? – meowgoesthedog Oct 04 '18 at 12:34
  • @meowgoesthedog, sorry for the misunderstanding in my post. But some bad 3D meshes do have non-planar polygons. Why does it matter? Ok, for each polygon, I have to build the projection matrix (5 to six vector operations, including square roots and inverse matrix calculations). This looks to me like a lot of overhead, doesn't it? – N. Wells Oct 04 '18 at 12:46
  • @N.Wells 1) the polygons don't necessarily have to be *exactly* planar, as long as they don't "fold back". 2) Still significantly cheaper than the triangulation routine, not to mention that you only need to normalise the normal once at the end (i.e. area-weighted mean). 3) You don't need to do any explicit matrix inversions; just note that the inverse of a rotation matrix is equal to its transpose (orthogonal). – meowgoesthedog Oct 04 '18 at 12:56
  • @meowgoesthedog, "the normal at end"? Any sample code? Mine is similar to this post: https://stackoverflow.com/questions/38588842/coordinate-system-transformation-3d-projection-to-2d-plane – N. Wells Oct 04 '18 at 13:01
  • Oops ignore that bit; but a few square roots for a matrix, which you only need to compute once per polygon, are still much cheaper than the work needed to triangulate it. The inverse in that post is a generalisation; for your case the matrix `M` is just a pure rotation matrix. Moral of the story: don't prematurely optimise – benchmark instead. – meowgoesthedog Oct 04 '18 at 13:04
  • @N. Wells, I am sorry but I still don't understand what your input is. From your question it seems that the input is a set of polygons and each of them is planar (although the plane might be different for each of those). Then, in a comment you state that "some bad 3D meshes do have non-planar polygons". Would you mind clarifying this, please? – Dominik Mokriš Oct 12 '18 at 18:57
  • Sorry, @DominikMokriš, I meant to say that my input meshes are "Good" well-formed meshes. Meshes with non-planar polygons are usually from a bad modeling. But you're right, a good triangulation program must take those "extreme" cases into account. – N. Wells Oct 15 '18 at 06:09

0 Answers0