src/common/triangulate.c

#ifndef lint
static const char RCSid[] = "$Id: triangulate.c,v 2.6 2021/04/19 19:40:03 greg Exp $";
#endif
/*
 *  triangulate.c
 *  
 *  Adapted by Greg Ward on 1/23/14.
 *  Fixes for polygons with seams/holes and co-linear vertices added
 *  by Nathaniel Jones on 12/21/16.
 *  Copyright 2016 Anyhere Software. All rights reserved.
 *
 */

/* COTD Entry submitted by John W. Ratcliff [[email protected]]

// ** THIS IS A CODE SNIPPET WHICH WILL EFFICIEINTLY TRIANGULATE ANY
// ** POLYGON/CONTOUR (without holes) AS A STATIC CLASS.  THIS SNIPPET
// ** IS COMPRISED OF 3 FILES, TRIANGULATE.H, THE HEADER FILE FOR THE
// ** TRIANGULATE BASE CLASS, TRIANGULATE.CPP, THE IMPLEMENTATION OF
// ** THE TRIANGULATE BASE CLASS, AND TEST.CPP, A SMALL TEST PROGRAM
// ** DEMONSTRATING THE USAGE OF THE TRIANGULATOR.  THE TRIANGULATE
// ** BASE CLASS ALSO PROVIDES TWO USEFUL HELPER METHODS, ONE WHICH
// ** COMPUTES THE AREA OF A POLYGON, AND ANOTHER WHICH DOES AN EFFICENT
// ** POINT IN A TRIANGLE TEST.
// ** SUBMITTED BY JOHN W. RATCLIFF ([email protected]) July 22, 2000
*/

#include <stdio.h>
#include <stdlib.h>
#include "triangulate.h"

#ifndef true
#define true    1
#define false   0
#endif

#define EPSILON         0.0000000001

static int
polySnip(const Vert2_list *contour, int u, int v, int w, int n, int *V)
{
  int p;
  double Ax, Ay, Bx, By, Cx, Cy, Px, Py, cross;

  Ax = contour->v[V[u]].mX;
  Ay = contour->v[V[u]].mY;

  Bx = contour->v[V[v]].mX;
  By = contour->v[V[v]].mY;

  Cx = contour->v[V[w]].mX;
  Cy = contour->v[V[w]].mY;

  cross = ((Bx - Ax)*(Cy - Ay)) - ((By - Ay)*(Cx - Ax));
  if (cross < EPSILON)
        return cross > -EPSILON ? -1 : false; /* Negative if colinear points */

  for (p=0;p<n;p++)
  {
    if( (p == u) | (p == v) | (p == w) ) continue;
    Px = contour->v[V[p]].mX;
    Py = contour->v[V[p]].mY;
    if ((Px == Ax) & (Py == Ay) || (Px == Bx) & (Py == By) ||
                (Px == Cx) & (Py == Cy)) continue; /* Handle donuts */
    if (insideTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,Py)) return false;
  }

  return true;
}

Vert2_list *
polyAlloc(int nv)
{
        Vert2_list      *pnew;

        if (nv < 3) return NULL;
        
        pnew = (Vert2_list *)malloc(sizeof(Vert2_list) + sizeof(Vert2)*(nv-3));
        if (pnew == NULL) return NULL;
        pnew->nv = nv;
        pnew->p = NULL;
        
        return pnew;
}

/*
  Area is positive if vertices listed counter-clockwise, negative if clockwise
 */
double
polyArea(const Vert2_list *contour)
{
  double A=0.0;
  int   p, q;

  for(p = contour->nv-1, q = 0; q < contour->nv; p=q++)
  {
    A += contour->v[p].mX*contour->v[q].mY - contour->v[q].mX*contour->v[p].mY;
  }
  return A*0.5;
}

/*
  InsideTriangle decides if a point P is Inside of the triangle
  defined by A, B, C.
*/
int
insideTriangle(double Ax, double Ay,
                      double Bx, double By,
                      double Cx, double Cy,
                      double Px, double Py)
{
  double ax, ay, bx, by, cx, cy, apx, apy, bpx, bpy, cpx, cpy;
  double cCROSSap, bCROSScp, aCROSSbp;

  ax = Cx - Bx;  ay = Cy - By;
  bx = Ax - Cx;  by = Ay - Cy;
  cx = Bx - Ax;  cy = By - Ay;
  apx= Px - Ax;  apy= Py - Ay;
  bpx= Px - Bx;  bpy= Py - By;
  cpx= Px - Cx;  cpy= Py - Cy;

  aCROSSbp = ax*bpy - ay*bpx;
  cCROSSap = cx*apy - cy*apx;
  bCROSScp = bx*cpy - by*cpx;

  return ((aCROSSbp >= 0.0) & (bCROSScp >= 0.0) & (cCROSSap >= 0.0));
};

int
polyTriangulate(const Vert2_list *contour, tri_out_t *cb)
{
  /* allocate and initialize list of Vertices in polygon */

  int nv, u, v, w, count, result;
  int *V;

  if ( contour->nv < 3 ) return false;

  V = (int *)malloc(sizeof(int)*contour->nv);
  if (V == NULL) return false;

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

  if ( polyArea(contour) > 0.0 )
    for (v=0; v<contour->nv; v++) V[v] = v;
  else
    for (v=0; v<contour->nv; v++) V[v] = (contour->nv-1)-v;

  nv = contour->nv;

  /*  remove nv-2 Vertices, creating 1 triangle every time */
  count = 2*nv;   /* error detection */

  v = nv-1;
  while (nv > 2)
  {
    /* if we loop, it is probably a non-simple polygon */
    if (count-- <= 0)
    {
      /* Triangulate: ERROR - probable bad polygon */
      free(V);
      return false;
    }

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

    result = polySnip(contour, u, v, w, nv, V);
    if (result > 0) /* successfully found a triangle */
    {
      /* output Triangle */
      if (!(*cb)(contour, V[u], V[v], V[w])) {
        free(V);
        return false;
      }
    }
    if (result) /* successfully found a triangle or three consecutive colinear points */
    {
      int s,t;

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

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

  free(V);

  return true;
}
Revision:	2.7
Committed:	Mon Apr 19 20:37:04 2021 UTC (3 years ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, HEAD
Changes since 2.6:	+2 -2 lines
Log Message:	fix: Repaired issue in last change with colinear vertices
#	Content
1	#ifndef lint
2	static const char RCSid[] = "$Id: triangulate.c,v 2.6 2021/04/19 19:40:03 greg Exp $";
3	#endif
4	/*
5	* triangulate.c
6	*
7	* Adapted by Greg Ward on 1/23/14.
8	* Fixes for polygons with seams/holes and co-linear vertices added
9	* by Nathaniel Jones on 12/21/16.
10	* Copyright 2016 Anyhere Software. All rights reserved.
11	*
12	*/
13
14	/* COTD Entry submitted by John W. Ratcliff [[email protected]]
15
16	// ** THIS IS A CODE SNIPPET WHICH WILL EFFICIEINTLY TRIANGULATE ANY
17	// ** POLYGON/CONTOUR (without holes) AS A STATIC CLASS. THIS SNIPPET
18	// ** IS COMPRISED OF 3 FILES, TRIANGULATE.H, THE HEADER FILE FOR THE
19	// ** TRIANGULATE BASE CLASS, TRIANGULATE.CPP, THE IMPLEMENTATION OF
20	// ** THE TRIANGULATE BASE CLASS, AND TEST.CPP, A SMALL TEST PROGRAM
21	// ** DEMONSTRATING THE USAGE OF THE TRIANGULATOR. THE TRIANGULATE
22	// ** BASE CLASS ALSO PROVIDES TWO USEFUL HELPER METHODS, ONE WHICH
23	// ** COMPUTES THE AREA OF A POLYGON, AND ANOTHER WHICH DOES AN EFFICENT
24	// ** POINT IN A TRIANGLE TEST.
25	// ** SUBMITTED BY JOHN W. RATCLIFF ([email protected]) July 22, 2000
26	*/
27
28	#include <stdio.h>
29	#include <stdlib.h>
30	#include "triangulate.h"
31
32	#ifndef true
33	#define true 1
34	#define false 0
35	#endif
36
37	#define EPSILON 0.0000000001
38
39	static int
40	polySnip(const Vert2_list contour, int u, int v, int w, int n, int V)
41	{
42	int p;
43	double Ax, Ay, Bx, By, Cx, Cy, Px, Py, cross;
44
45	Ax = contour->v[V[u]].mX;
46	Ay = contour->v[V[u]].mY;
47
48	Bx = contour->v[V[v]].mX;
49	By = contour->v[V[v]].mY;
50
51	Cx = contour->v[V[w]].mX;
52	Cy = contour->v[V[w]].mY;
53
54	cross = ((Bx - Ax)(Cy - Ay)) - ((By - Ay)(Cx - Ax));
55	if (cross < EPSILON)
56	return cross > -EPSILON ? -1 : false; /* Negative if colinear points */
57
58	for (p=0;p<n;p++)
59	{
60	if( (p == u) \| (p == v) \| (p == w) ) continue;
61	Px = contour->v[V[p]].mX;
62	Py = contour->v[V[p]].mY;
63	if ((Px == Ax) & (Py == Ay) \|\| (Px == Bx) & (Py == By) \|\|
64	(Px == Cx) & (Py == Cy)) continue; /* Handle donuts */
65	if (insideTriangle(Ax,Ay,Bx,By,Cx,Cy,Px,Py)) return false;
66	}
67
68	return true;
69	}
70
71	Vert2_list *
72	polyAlloc(int nv)
73	{
74	Vert2_list *pnew;
75
76	if (nv < 3) return NULL;
77
78	pnew = (Vert2_list )malloc(sizeof(Vert2_list) + sizeof(Vert2)(nv-3));
79	if (pnew == NULL) return NULL;
80	pnew->nv = nv;
81	pnew->p = NULL;
82
83	return pnew;
84	}
85
86	/*
87	Area is positive if vertices listed counter-clockwise, negative if clockwise
88	*/
89	double
90	polyArea(const Vert2_list *contour)
91	{
92	double A=0.0;
93	int p, q;
94
95	for(p = contour->nv-1, q = 0; q < contour->nv; p=q++)
96	{
97	A += contour->v[p].mXcontour->v[q].mY - contour->v[q].mXcontour->v[p].mY;
98	}
99	return A*0.5;
100	}
101
102	/*
103	InsideTriangle decides if a point P is Inside of the triangle
104	defined by A, B, C.
105	*/
106	int
107	insideTriangle(double Ax, double Ay,
108	double Bx, double By,
109	double Cx, double Cy,
110	double Px, double Py)
111	{
112	double ax, ay, bx, by, cx, cy, apx, apy, bpx, bpy, cpx, cpy;
113	double cCROSSap, bCROSScp, aCROSSbp;
114
115	ax = Cx - Bx; ay = Cy - By;
116	bx = Ax - Cx; by = Ay - Cy;
117	cx = Bx - Ax; cy = By - Ay;
118	apx= Px - Ax; apy= Py - Ay;
119	bpx= Px - Bx; bpy= Py - By;
120	cpx= Px - Cx; cpy= Py - Cy;
121
122	aCROSSbp = axbpy - aybpx;
123	cCROSSap = cxapy - cyapx;
124	bCROSScp = bxcpy - bycpx;
125
126	return ((aCROSSbp >= 0.0) & (bCROSScp >= 0.0) & (cCROSSap >= 0.0));
127	};
128
129	int
130	polyTriangulate(const Vert2_list contour, tri_out_t cb)
131	{
132	/* allocate and initialize list of Vertices in polygon */
133
134	int nv, u, v, w, count, result;
135	int *V;
136
137	if ( contour->nv < 3 ) return false;
138
139	V = (int )malloc(sizeof(int)contour->nv);
140	if (V == NULL) return false;
141
142	/* we want a counter-clockwise polygon in V */
143
144	if ( polyArea(contour) > 0.0 )
145	for (v=0; v<contour->nv; v++) V[v] = v;
146	else
147	for (v=0; v<contour->nv; v++) V[v] = (contour->nv-1)-v;
148
149	nv = contour->nv;
150
151	/* remove nv-2 Vertices, creating 1 triangle every time */
152	count = 2nv; / error detection */
153
154	v = nv-1;
155	while (nv > 2)
156	{
157	/* if we loop, it is probably a non-simple polygon */
158	if (count-- <= 0)
159	{
160	/* Triangulate: ERROR - probable bad polygon */
161	free(V);
162	return false;
163	}
164
165	/* three consecutive vertices in current polygon, <u,v,w> */
166	u = v ; u = (nv > u); / previous */
167	v = u+1; v = (nv > v); / new v */
168	w = v+1; w = (nv > w); / next */
169
170	result = polySnip(contour, u, v, w, nv, V);
171	if (result > 0) /* successfully found a triangle */
172	{
173	/* output Triangle */
174	if (!(*cb)(contour, V[u], V[v], V[w])) {
175	free(V);
176	return false;
177	}
178	}
179	if (result) /* successfully found a triangle or three consecutive colinear points */
180	{
181	int s,t;
182
183	/* remove v from remaining polygon */
184	for(s=v,t=v+1;t<nv;s++,t++) V[s] = V[t]; nv--;
185
186	/* reset error detection counter */
187	count = 2*nv;
188	}
189	}
190
191	free(V);
192
193	return true;
194	}