src/common/triangulate.c

#ifndef lint
static const char RCSid[] = "$Id: triangulate.c,v 2.2 2014/01/24 01:26:44 greg Exp $";
#endif
/*
 *  triangulate.c
 *  
 *  Adapted by Greg Ward on 1/23/14.
 *  Copyright 2014 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

static const double 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;

  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;

  if ( EPSILON > (((Bx-Ax)*(Cy-Ay)) - ((By-Ay)*(Cx-Ax))) ) return false;

  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 (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;
}

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, m, u, v, w, count;
  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 ( 0.0 < polyArea(contour) )
    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 */

  for(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 */
      return false;
    }

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

    if ( polySnip(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 */
      if (!(*cb)(contour, a, b, c)) return false;
 
      m++;

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