star-cpr

commit e81c4569e9b09f3499df965c7819d60bae8589ac
Author: Vincent Forest <vincent.forest@meso-star.com>
Date:   Fri, 19 Aug 2016 18:07:02 +0200

First draft of the Weiler Atherton clipping algorithm

Diffstat:
A
main.c
 | 
361
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

1 file changed, 361 insertions(+), 0 deletions(-)
diff --git a/main.c b/main.c
@@ -0,0 +1,361 @@
+#include <rsys/double2.h>
+#include <rsys/stretchy_array.h>
+
+#include <stdio.h>
+#include <math.h>
+
+#define NDISKS 10
+#define N 50
+#define LEN 5.0
+#define DELTA (LEN / NDISKS)
+
+struct mesh {
+  double* coords;
+  size_t npoints;
+  int* tris;
+  size_t ntris;
+};
+
+struct hit {
+  double pos[2]; /* Position */
+  double t0; /* Parametric intersection on the polygon */
+  double t1; /* Parametric intersection on the clipping polygon */
+};
+
+struct node {
+  int prev;
+  int next;
+  int link;
+  double pos[2];
+  double t;
+};
+
+#define MESH_NULL { NULL, 0, NULL, 0 }
+
+static void
+dump(FILE* stream, const struct mesh* poly)
+{
+  size_t i;
+  for(i = 0; i < poly->npoints; ++i) {
+    fprintf(stream, "v %g %g %g\n",
+      poly->coords[i*3+0],
+      poly->coords[i*3+1],
+      poly->coords[i*3+2]);
+  }
+
+  for(i = 0; i < poly->ntris; ++i) {
+    fprintf(stream, "f %d %d %d\n",
+      poly->tris[i*3+0] + 1,
+      poly->tris[i*3+1] + 1,
+      poly->tris[i*3+2] + 1);
+  }
+}
+
+static int
+intersect
+  (const double a0[2], const double a1[2], /* 1st segment */
+   const double b0[2], const double b1[2], /* 2nd segment */
+   struct hit* hit)
+{
+  double p0[2], p1[2];
+  double d0[2], d1[2];
+  double tmp[2];
+  double a, b, t;
+
+  d2_set(p0, a0);
+  d2_set(p1, b0);
+  d2_sub(d0, a1, a0);
+  d2_sub(d1, b1, b0);
+
+  d2_sub(tmp, p1, p0);
+  a = d2_cross(tmp, d1);
+  b = d2_cross(d0, d1);
+  if(b == 0) return 0;
+
+  hit->t0 = a / b;
+  d2_add(hit->pos, p0, d2_muld(hit->pos, d0, hit->t0));
+  if(d1[0] != 0) {
+    hit->t1 = (hit->pos[0] - p1[0]) / d1[0];
+  } else {
+    hit->t1 = (hit->pos[1] - p1[1]) / d1[1];
+  }
+
+  if(hit->t0 < 0 || hit->t0 > 1
+  || hit->t1 < 0 || hit->t1 > 1)
+    return 0; /* No intersection */
+
+#ifndef NDEBUG
+  {
+    double pt1[2];
+    d2_add(pt1, p1, d2_muld(pt1, d1, hit->t1));
+    ASSERT(d2_eq_eps(hit->pos, pt1, 1.e-6));
+  }
+#endif
+  return 1;
+}
+
+static void
+clip_weiler_atherton
+  (double v[3][2], /* CW order */
+   const double* clip, /* CCW order */
+   const size_t nclips)
+{
+  struct node* poly_nodes = NULL;
+  struct node* clip_nodes = NULL;
+  double* output = NULL;
+  int* is_inside = NULL;
+  struct node* list[2];
+  int prev, next;
+  int i, j, n;
+
+  /* Register the polygon nodes */
+  for(i=0, prev=2, next=1, n=3; i < n; ++i) {
+    struct node* node = sa_add(poly_nodes, 1);
+    node->prev = prev;
+    node->next = next;
+    node->link = -1;
+    node->pos[0] = v[i][0];
+    node->pos[1] = v[i][1];
+
+    prev = (prev + 1) % n;
+    next = (next + 1) % n;
+  }
+
+  /* Register the clipping nodes */
+  for(i=0, prev=nclips-1, next=1, n=nclips; i < n; ++i) {
+    struct node* node = sa_add(clip_nodes, 1);
+    node->prev = prev;
+    node->next = next;
+    node->link = -1;
+    node->pos[0] = clip[i*2 + 0];
+    node->pos[1] = clip[i*2 + 1];
+
+    prev = (prev + 1) % n;
+    next = (next + 1) % n;
+  }
+
+  /* Add the intersection nodes into the polygon and clipping nodes */
+  for(i=0; i < 3; ++i) {
+    double poly_edge[2][2];
+
+    /* Fetch the edge vertices of the triangle */
+    d2_set(poly_edge[0], v[i]);
+    d2_set(poly_edge[1], v[(i+1) % 3]);
+
+    for(j=0; j < nclips; ++j) {
+      struct node* node;
+      struct node* clip_node;
+      struct node* poly_node;
+      struct hit hit;
+      size_t ipoly_node;
+      size_t iclip_node;
+      size_t inode;
+      double clip_edge[2][2];
+      int has_intersection;
+
+      /* Fetch the vertices of the clip polygon edge */
+      d2_set(clip_edge[0], clip + j*2);
+      d2_set(clip_edge[1], clip + ((j+1)%nclips)*2);
+
+      /* Find the intersection between the triangle and the clip polygon edges */
+      has_intersection = intersect
+        (poly_edge[0], poly_edge[1], clip_edge[0], clip_edge[1], &hit);
+      if(!has_intersection) continue; /* No edge intersection */
+
+      /* Register the intersection into the list of poly/clip nodes */
+      poly_node = sa_add(poly_nodes, 1);
+      clip_node = sa_add(clip_nodes, 1);
+      ipoly_node = poly_node - poly_nodes;
+      iclip_node = clip_node - clip_nodes;
+
+      /* Insert the new poly node into the list of triangle nodes */
+      node = poly_nodes + poly_nodes[i].next;
+      while(node->link >= 0) {
+        if(node->t > hit.t0) break;
+        node = poly_nodes + node->next;
+      }
+      inode = node - poly_nodes;
+      poly_node->next = inode;
+      poly_node->prev = node->prev;
+      poly_node->link = iclip_node; /* Id of the corresponding clip node */
+      poly_node->pos[0] = hit.pos[0];
+      poly_node->pos[1] = hit.pos[1];
+      poly_node->t = hit.t0;
+      poly_nodes[node->prev].next = ipoly_node;
+      node->prev = ipoly_node;
+
+      /* Insert the new clip node into the list of clip nodes */
+      node = clip_nodes + clip_nodes[j].next;
+      while(node->link >= 0) {
+        if(node->t > hit.t1) break;
+        node = clip_nodes + node->next;
+      }
+      inode = node - clip_nodes;
+      clip_node->next = inode;
+      clip_node->prev = node->prev;
+      clip_node->link = ipoly_node; /* Id of the corresponding poly node */
+      clip_node->pos[0] = hit.pos[0];
+      clip_node->pos[1] = hit.pos[1];
+      clip_node->t = hit.t1;
+      clip_nodes[node->prev].next = iclip_node;
+      node->prev = iclip_node;
+    }
+  }
+
+  /* Define which triangle vertices are inside the clip polygon */
+  for(i=0; i < 3; ++i) {
+    int inside = 1;
+    for(j=0; inside && j < nclips; ++j) {
+      const double* c0 = clip + j*2/*#coords per vertex*/;
+      const double* c1 = clip + ((j+1)%nclips)*2/*#coords per vertex*/;
+      double e0[2], e1[2];
+
+      d2_sub(e0, c0, v[i]);
+      d2_sub(e1, c1, v[i]);
+      inside = d2_cross(e1, e0) < 0;
+    }
+    sa_push(is_inside, inside);
+  }
+
+  list[0] = poly_nodes;
+  list[1] = clip_nodes;
+
+  /* Note that hit nodes are registered *after* the nodes of the polygon */
+  if(sa_size(poly_nodes) == 3 && !is_inside[0]) {
+    for(i=0; i < 3; ++i) {
+      double x = v[i][0];
+      double y = v[i][1];
+      double z = x*x + y*y;
+      printf("%g %g %g\n", x, y, z);
+    }
+  } else {
+    for(i=3, n=sa_size(poly_nodes); i < n; ++i) {
+      struct node* node_start = poly_nodes + poly_nodes[i].next;
+      struct node* node = node_start;
+      int ilist = 0;
+
+      if(poly_nodes[poly_nodes[i].next].link >= 0 || is_inside[poly_nodes[i].next])
+        continue;
+
+      sa_clear(output);
+      node_start = poly_nodes + poly_nodes[i].next;
+      node = node_start;
+      ilist = 0;
+
+      do {
+        d2_set(sa_add(output, 2), node->pos);
+        node = list[ilist] + node->next;
+
+        if(node->link >= 0) {
+          ilist = !ilist;
+          node = list[ilist] + node->link;
+        }
+      } while(node != node_start);
+
+      /* Print the output polygon */
+      for(j = 0; j < sa_size(output)/2; ++j) {
+        double x = output[j*2+0];
+        double y = output[j*2+1];
+        double z = x*x + y*y;
+        printf("%g %g %g\n", x, y, z);
+      }
+      printf("\n");
+    }
+  }
+
+  sa_release(is_inside);
+  sa_release(output);
+  sa_release(poly_nodes);
+  sa_release(clip_nodes);
+}
+
+static void
+test_paraboloid(void)
+{
+  struct mesh paraboloid = MESH_NULL;
+  const double xlen;
+  double* clip = NULL;
+  int i, j;
+
+  for(i = 0; i < NDISKS+1; ++i) {
+    const double r = (i+1/*FIXME*/)*DELTA;
+    for(j = 0; j < N; ++j) {
+      const double theta = (double)j / (double)N * 2 * M_PI;
+      const double x = r * cos(theta);
+      const double y = r * sin(theta);
+      const double z = x*x + y*y;
+      sa_push(paraboloid.coords, x);
+      sa_push(paraboloid.coords, y);
+      sa_push(paraboloid.coords, z);
+    }
+  }
+  /* Polar point */
+  sa_push(paraboloid.coords, 0);
+  sa_push(paraboloid.coords, 0);
+  sa_push(paraboloid.coords, 0);
+  paraboloid.npoints = sa_size(paraboloid.coords) / 3 /* #coords per point */;
+
+  for(i = 0; i < NDISKS; ++i) {
+    const int offset = (i+1)*N;
+    for(j = 0; j < N; ++j) {
+      const int id0 = j + offset;
+      const int id1 = ((j + 1) % N) + offset;
+      const int id2 = id0 - N;
+      const int id3 = id1 - N;
+
+      sa_push(paraboloid.tris, id0);
+      sa_push(paraboloid.tris, id1);
+      sa_push(paraboloid.tris, id2);
+
+      sa_push(paraboloid.tris, id2);
+      sa_push(paraboloid.tris, id1);
+      sa_push(paraboloid.tris, id3);
+    }
+  }
+
+  /* Polar triangles */
+  for(j = 0; j < N; ++j) {
+    const int id0 = j;
+    const int id1 = (j + 1) % N;
+    const int id2 = sa_size(paraboloid.coords)/3 - 1;
+
+    sa_push(paraboloid.tris, id0);
+    sa_push(paraboloid.tris, id1);
+    sa_push(paraboloid.tris, id2);
+  }
+  paraboloid.ntris = sa_size(paraboloid.tris) / 3 /* #ids per triangle */;
+  /*dump(stdout, &paraboloid);*/
+
+  /* Build the clip polygon */
+  d2(sa_add(clip, 2), 0, -2.5);
+  d2(sa_add(clip, 2), 4, 2);
+  d2(sa_add(clip, 2), -3.5, 0);
+
+  for(i = 0; i < paraboloid.ntris; ++i) {
+    double tri[3][2];
+    d2_set(tri[0], paraboloid.coords + paraboloid.tris[i*3+0] * 3);
+    d2_set(tri[1], paraboloid.coords + paraboloid.tris[i*3+1] * 3);
+    d2_set(tri[2], paraboloid.coords + paraboloid.tris[i*3+2] * 3);
+    clip_weiler_atherton(tri, clip, sa_size(clip)/2);
+  }
+
+  sa_release(paraboloid.coords);
+  sa_release(paraboloid.tris);
+}
+
+static void
+test_simple(void)
+{
+  double triangle[3][2] = { {0.0, 0.0}, {0.0, 1.0}, {1.0, 0.0} };
+  double clip[6] = { -1.0, 0.25, 1.0, 0.25, 1, 0.75 };
+  clip_weiler_atherton(triangle, clip, 3);
+}
+
+int
+main(int argc, char** argv)
+{
+  test_paraboloid();
+  /*test_simple();*/
+  return 0;
+}
+