Une implémentation Java de l'algorithme de Triangulation de Delaunay incrémentale inventé par Leonidas Guibas et Jorge Stolfi.



La classe QuadEdge qui décrit la structure de données duale (face/segment):
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
 
import java.awt.Point;
 
/**
 * Quad-Edge data structure
 *
 * @author Code by X.Philippeau - Structure by Guibas and Stolfi
 *
 * @see Primitives for the Manipulation of General Subdivisions
 *      and the Computation of Voronoi Diagrams (Leonidas Guibas,Jorge Stolfi)
 */
public class QuadEdge {
 
	// pointer to the next (direct order) QuadEdge
	private QuadEdge onext;
 
	// pointer to the dual QuadEdge (faces graph <-> edges graph)
	private QuadEdge rot;
 
	// origin point of the edge/face
	private Point orig;
 
	// marker for triangle generation
	public boolean mark=false;
 
	/**
         * (private) constructor. Use makeEdge() to create a new QuadEdge
         *
         * @param onext pointer to the next QuadEdge on the ring
         * @param rot   pointer to the next (direct order) crossing edge
         * @param orig  Origin point
         */
	private QuadEdge(QuadEdge Onext, QuadEdge rot, Point orig) {
		this.onext = Onext;
		this.rot = rot;
		this.orig = orig;
	}
 
	// ----------------------------------------------------------------
	//                             Getter/Setter
	// ----------------------------------------------------------------
 
	public QuadEdge onext() {
		return onext;
	}
 
	public QuadEdge rot() {
		return rot;
	}
 
	public Point orig() {
		return orig;
	}
 
	public void setOnext(QuadEdge next) {
		onext = next;
	}
 
	public void setRot(QuadEdge rot) {
		this.rot = rot;
	}
 
	public void setOrig(Point p) {
		this.orig = p;
	}
 
	// ----------------------------------------------------------------
	//                      QuadEdge Navigation
	// ----------------------------------------------------------------
 
	/**
         * @return the symetric (reverse) QuadEdge
         */
	public QuadEdge sym() {
		return rot.rot();
	}
 
	/**
         * @return the other extremity point
         */
	public Point dest() {
		return sym().orig();
	}
 
	/**
         * @return the symetric dual QuadEdge
         */
	public QuadEdge rotSym() {
		return rot.sym();
	}
 
	/**
         * @return the previous QuadEdge (pointing to this.orig)
         */
	public QuadEdge oprev() {
		return rot.onext().rot();
	}
 
	/**
         * @return the previous QuadEdge starting from dest()
         */
	public QuadEdge dprev() {
		return rotSym().onext().rotSym();
	}
 
	/**
         * @return the next QuadEdge on left Face
         */
	public QuadEdge lnext() {
		return rotSym().onext().rot();
	}
 
	/**
         * @return the previous QuadEdge on left Face
         */
	public QuadEdge lprev() {
		return onext().sym();
	}
 
 
	// ************************** STATIC ******************************
 
 
	/**
         * Create a new edge (i.e. a segment)
         *
         * @param  orig origin of the segment
         * @param  dest end of the segment
         * @return the QuadEdge of the origin point
         */
	public static QuadEdge makeEdge(Point orig, Point dest) {
		QuadEdge q0 = new QuadEdge(null, null, orig);
		QuadEdge q1 = new QuadEdge(null, null, null);
		QuadEdge q2 = new QuadEdge(null, null, dest);
		QuadEdge q3 = new QuadEdge(null, null, null);
 
		// create the segment
		q0.onext = q0; q2.onext = q2; // lonely segment: no "next" quadedge
		q1.onext = q3; q3.onext = q1; // in the dual: 2 communicating facets
 
		// dual switch
		q0.rot = q1; q1.rot = q2;
		q2.rot = q3; q3.rot = q0;
 
		return q0;
	}
 
	/**
         * attach/detach the two edges = combine/split the two rings in the dual space
         *
         * @param q1,q2 the 2 QuadEdge to attach/detach
         */
	public static void splice(QuadEdge a, QuadEdge b) {
		QuadEdge alpha = a.onext().rot();
		QuadEdge beta  = b.onext().rot();
 
		QuadEdge t1 = b.onext();
		QuadEdge t2 = a.onext();
		QuadEdge t3 = beta.onext();
		QuadEdge t4 = alpha.onext();
 
		a.setOnext(t1);
		b.setOnext(t2);
		alpha.setOnext(t3);
		beta.setOnext(t4);
	}
 
	/**
         * Create a new QuadEdge by connecting 2 QuadEdges
     *
         * @param e1,e2 the 2 QuadEdges to connect
         * @return the new QuadEdge
         */
	public static QuadEdge connect(QuadEdge e1, QuadEdge e2) {
		QuadEdge q = makeEdge(e1.dest(), e2.orig());
		splice(q, e1.lnext());
		splice(q.sym(), e2);
		return q;
	}
 
	public static void swapEdge(QuadEdge e) {
		QuadEdge a = e.oprev();
		QuadEdge b = e.sym().oprev();
		splice(e, a);
		splice(e.sym(), b);
		splice(e, a.lnext());
		splice(e.sym(), b.lnext());
		e.orig = a.dest();
		e.sym().orig = b.dest();
	}
 
	/**
         * Delete a QuadEdge
         *
         * @param q the QuadEdge to delete
         */
	public static void deleteEdge(QuadEdge q) {
		splice(q, q.oprev());
		splice(q.sym(), q.sym().oprev());
	}
 
	// ----------------------------------------------------------------
	//                      Geometric computation
	// ----------------------------------------------------------------
 
	/**
         * Test if the Point p is on the line porting the edge
         *
         * @param e QuadEdge
         * @param p Point to test
         * @return true/false
         */
	public static boolean isOnLine(QuadEdge e, Point p) {
		// test if the vector product is zero
		if ((p.x-e.orig().x)*(p.y-e.dest().y)==(p.y-e.orig().y)*(p.x-e.dest().x))
			return true;
		return false;
	}
 
	/**
         * Test if the Point p is at the right of the QuadEdge q.
         *
         * @param QuadEdge reference
         * @param p Point to test
         * @return true/false
         */
	public static boolean isAtRightOf(QuadEdge q, Point p) {
		return isCounterClockwise(p, q.dest(), q.orig());
	}
 
	/** return true if a, b and c turn in Counter Clockwise direction
         *
         * @param a,b,c the 3 points to test
         * @return true if a, b and c turn in Counter Clockwise direction
         */
	public static boolean isCounterClockwise(Point a, Point b, Point c) {
		// test the sign of the determinant of ab x cb
		if ( (a.x - b.x)*(b.y - c.y) > (a.y - b.y)*(b.x - c.x) ) return true;
		return false;
	}
 
	/**
         * The Delaunay criteria:
         *
         *   test if the point d is inside the circumscribed circle of triangle a,b,c
         *
         * @param a,b,c triangle
         * @param d point to test
         * @return  true/false
         */
	public static boolean inCircle(Point a, Point b, Point c, Point d) {
		/*
		 if "d" is strictly INSIDE the circle, then
 
		     |d² dx dy 1|
             |a² ax ay 1|
		 det |b² bx by 1| < 0
		     |c² cx cy 1|
 
		*/
		long a2 = a.x*a.x + a.y*a.y;
		long b2 = b.x*b.x + b.y*b.y;
		long c2 = c.x*c.x + c.y*c.y;
		long d2 = d.x*d.x + d.y*d.y;
 
		long det44 = 0;
		det44 += d2  * det33( a.x,a.y, 1,  b.x,b.y, 1,  c.x,c.y, 1 );
		det44 -= d.x * det33( a2 ,a.y, 1,  b2 ,b.y, 1,  c2 ,c.y, 1 );
		det44 += d.y * det33( a2 ,a.x, 1,  b2 ,b.x, 1,  c2 ,c.x, 1 );
		det44 -= 1   * det33( a2,a.x,a.y,  b2,b.x,b.y,  c2,c.x,c.y );
 
		if (det44<0) return true;
		return false;
	}
 
	private static long det33( long... m ) {
		long det33=0;
		det33 += m[0] * (m[4]*m[8] - m[5]*m[7]);
		det33 -= m[1] * (m[3]*m[8] - m[5]*m[6]);
		det33 += m[2] * (m[3]*m[7] - m[4]*m[6]);
		return det33;
	}
}

La classe Delaunay qui effectue la triangulation:
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
 
import java.awt.Point;
import java.util.ArrayList;
import java.util.List;
 
/**
 * Incremental Delaunay Triangulation
 *
 * @author Java-code by X.Philippeau - Pseudo-code by Guibas and Stolfi
 *
 * @see Primitives for the Manipulation of General Subdivisions
 *      and the Computation of Voronoi Diagrams (Leonidas Guibas,Jorge Stolfi)
 */
public class Delaunay {
 
	// starting edge for walk (see locate() method)
	private QuadEdge startingEdge = null;
 
	// list of quadEdge belonging to Delaunay triangulation
	private List<QuadEdge> quadEdge = new ArrayList<QuadEdge>();
 
	// Bounding box of the triangulation
	class BoundingBox {
		int minx,miny,maxx,maxy;
		Point a = new Point(); // lower left
		Point b = new Point(); // lower right
		Point c = new Point(); // upper right
		Point d = new Point(); // upper left
	}
	private BoundingBox bbox = new BoundingBox();
 
	/**
         * Constuctor:
         */
	public Delaunay() {
 
		bbox.minx = Integer.MAX_VALUE;
		bbox.maxx = Integer.MIN_VALUE;
		bbox.miny = Integer.MAX_VALUE;
		bbox.maxy = Integer.MIN_VALUE;
 
		// create the QuadEdge graph of the bounding box
		QuadEdge ab = QuadEdge.makeEdge(bbox.a,bbox.b);
		QuadEdge bc = QuadEdge.makeEdge(bbox.b,bbox.c);
		QuadEdge cd = QuadEdge.makeEdge(bbox.c,bbox.d);
		QuadEdge da = QuadEdge.makeEdge(bbox.d,bbox.a);
		QuadEdge.splice(ab.sym(), bc);
		QuadEdge.splice(bc.sym(), cd);
		QuadEdge.splice(cd.sym(), da);
		QuadEdge.splice(da.sym(), ab);
 
		this.startingEdge = ab;
	}
 
	/**
         * update the dimension of the bounding box
         *
         * @param minx,miny,maxx,maxy summits of the rectangle
         */
	public void setBoundigBox(int minx,int miny,int maxx,int maxy) {
		// update saved values
		bbox.minx=minx; bbox.maxx=maxx;
		bbox.miny=miny; bbox.maxy=maxy;
 
		// extend the bounding-box to surround min/max
		int centerx = (minx+maxx)/2;
		int centery = (miny+maxy)/2;
		int x_min = (int)((minx-centerx-1)*10+centerx);
		int x_max = (int)((maxx-centerx+1)*10+centerx);
		int y_min = (int)((miny-centery-1)*10+centery);
		int y_max = (int)((maxy-centery+1)*10+centery);
 
		// set new positions
		bbox.a.x = x_min;	bbox.a.y = y_min;
		bbox.b.x = x_max;	bbox.b.y = y_min;
		bbox.c.x = x_max;	bbox.c.y = y_max;
		bbox.d.x = x_min;	bbox.d.y = y_max;
	}
 
	// update the size of the bounding box (cf locate() method)
	private void updateBoundigBox(Point p) {
		int minx = Math.min(bbox.minx, p.x);
		int maxx = Math.max(bbox.maxx, p.x);
		int miny = Math.min(bbox.miny, p.y);
		int maxy = Math.max(bbox.maxy, p.y);
		setBoundigBox(minx, miny, maxx, maxy);
		//System.out.println("resizing bounding-box: "+minx+" "+miny+" "+maxx+" "+maxy);
	}
 
	/**
         * Returns an edge e of the triangle containing the point p
         * (Guibas and Stolfi)
         *
         * @param p the point to localte
         * @return the edge of the triangle
         */
	private QuadEdge locate(Point p) {
 
		/* outside the bounding box ? */
		if ( p.x<bbox.minx || p.x>bbox.maxx || p.y<bbox.miny || p.y>bbox.maxy ) {
			updateBoundigBox(p);
		}
 
		QuadEdge e = startingEdge;
		while (true) {
			/* duplicate point ? */
			if(p.x==e.orig().x && p.y==e.orig().y) return e;
			if(p.x==e.dest().x && p.y==e.dest().y) return e;
 
			/* walk */
			if (QuadEdge.isAtRightOf(e,p))
				e = e.sym();
			else if (!QuadEdge.isAtRightOf(e.onext(),p))
				e = e.onext();
			else if (!QuadEdge.isAtRightOf(e.dprev(),p))
				e = e.dprev();
			else
				return e;
		}
	}
 
	/**
         *  Inserts a new point into a Delaunay triangulation
         *  (Guibas and Stolfi)
         *
         * @param p the point to insert
         */
	public void insertPoint(Point p) {
		QuadEdge e = locate(p);
 
		// point is a duplicate -> nothing to do
		if(p.x==e.orig().x && p.y==e.orig().y) return;
		if(p.x==e.dest().x && p.y==e.dest().y) return;
 
		// point is on an existing edge -> remove the edge
		if (QuadEdge.isOnLine(e,p)) {
			e = e.oprev();
			this.quadEdge.remove(e.onext().sym());
			this.quadEdge.remove(e.onext());
			QuadEdge.deleteEdge(e.onext());
		}
 
		// Connect the new point to the vertices of the containing triangle
		// (or quadrilateral in case of the point is on an existing edge)
		QuadEdge base = QuadEdge.makeEdge(e.orig(),p);
		this.quadEdge.add(base);
 
		QuadEdge.splice(base, e);
		this.startingEdge = base;
		do {
			base = QuadEdge.connect(e, base.sym());
			this.quadEdge.add(base);
			e = base.oprev();
		} while (e.lnext() != startingEdge);
 
		// Examine suspect edges to ensure that the Delaunay condition is satisfied.
		do {
			QuadEdge t = e.oprev();
 
			if (QuadEdge.isAtRightOf(e,t.dest()) &&
					QuadEdge.inCircle(e.orig(), t.dest(), e.dest(), p)) {
				// flip triangles
				QuadEdge.swapEdge(e);
				e = e.oprev();
			}
			else if (e.onext() == startingEdge)
				return; // no more suspect edges
			else
				e = e.onext().lprev();  // next suspect edge
		} while (true);
	}
 
	/**
         *  compute and return the list of edges
         */
	public List<Point[]> computeEdges() {
		List<Point[]> edges = new ArrayList<Point[]>();
		// do not return edges pointing to/from surrouding triangle
		for(QuadEdge q:this.quadEdge) {
			if ( q.orig()==bbox.a || q.orig()==bbox.b || q.orig()==bbox.c || q.orig()==bbox.d ) continue;
			if ( q.dest()==bbox.a || q.dest()==bbox.b || q.dest()==bbox.c || q.dest()==bbox.d ) continue;
			edges.add(new Point[]{q.orig(),q.dest()});
		}
		return edges;
	}
 
	/**
         *  compute and return the list of triangles
         */
	public List<Point[]> computeTriangles() {
		List<Point[]> triangles = new ArrayList<Point[]>();
 
		// do not process edges pointing to/from surrouding triangle
		// --> mark them as already computed
		for(QuadEdge q:this.quadEdge) {
			q.mark = false; q.sym().mark=false;
			if ( q.orig()==bbox.a || q.orig()==bbox.b || q.orig()==bbox.c || q.orig()==bbox.d ) q.mark = true;
			if ( q.dest()==bbox.a || q.dest()==bbox.b || q.dest()==bbox.c || q.dest()==bbox.d ) q.sym().mark=true;
		}
 
		// compute the 2 triangles associated to each quadEdge
		for(QuadEdge qe:quadEdge) {
			// first triangle
			QuadEdge q1 = qe;
			QuadEdge q2 = q1.lnext();
			QuadEdge q3 = q2.lnext();
			if (!q1.mark && !q2.mark && !q3.mark) {
				triangles.add(new Point[] { q1.orig(),q2.orig(),q3.orig() });
			}
 
			// second triangle
			QuadEdge qsym1 = qe.sym();
			QuadEdge qsym2 = qsym1.lnext();
			QuadEdge qsym3 = qsym2.lnext();
			if (!qsym1.mark && !qsym2.mark && !qsym3.mark) {
				triangles.add(new Point[] { qsym1.orig(),qsym2.orig(),qsym3.orig() });
			}
 
			// mark as used
			qe.mark=true;
			qe.sym().mark=true;
		}
 
		return triangles;
	}
 
	public List<Point[]> computeVoronoi() {
		List<Point[]> voronoi = new ArrayList<Point[]>();
 
		// do not process edges pointing to/from surrouding triangle
		// --> mark them as already computed
		for(QuadEdge q:this.quadEdge) {
			q.mark = false; q.sym().mark=false;
			if ( q.orig()==bbox.a || q.orig()==bbox.b || q.orig()==bbox.c || q.orig()==bbox.d ) q.mark = true;
			if ( q.dest()==bbox.a || q.dest()==bbox.b || q.dest()==bbox.c || q.dest()==bbox.d ) q.sym().mark=true;
		}
 
		for(QuadEdge qe:quadEdge) {
 
			// walk throught left and right region
			for(int b=0;b<=1;b++) {
				QuadEdge qstart = (b==0)?qe:qe.sym();
				if (qstart.mark) continue;
 
				// new region start
				List<Point> poly = new ArrayList<Point>();
 
				// walk around region
				QuadEdge qregion = qstart;
				while(true) {
					qregion.mark=true;
 
					// compute CircumCenter if needed
					if (qregion.rot().orig()==null) {
						QuadEdge q1 = qregion;		Point p0=q1.orig();
						QuadEdge q2 = q1.lnext();	Point p1=q2.orig();
						QuadEdge q3 = q2.lnext();	Point p2=q3.orig();
 
						double ex=p1.x-p0.x, ey=p1.y-p0.y;
						double nx=p2.y-p1.y, ny=p1.x-p2.x;
						double dx=(p0.x-p2.x)*0.5, dy=(p0.y-p2.y)*0.5;
						double s=(ex*dx+ey*dy)/(ex*nx+ey*ny);
						double cx=(p1.x+p2.x)*0.5+s*nx;
						double cy=(p1.y+p2.y)*0.5+s*ny;
 
						Point p = new Point( (int)cx,(int)cy );
						qregion.rot().setOrig(p);
					}
 
					poly.add(qregion.rot().orig());
 
					qregion = qregion.onext();
					if (qregion==qstart) break;
				}
 
				// add region to output list
				voronoi.add(poly.toArray(new Point[0]));
			}
		}
		return voronoi;
	}
}

et enfin un exemple d'utilisation:
Code java : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 
/* creation de l'instance */
Delaunay delaunay = new Delaunay();
 
/* ajout des points (calcul incrémental) */
delaunay.insertPoint( new Point(12,98) );
delaunay.insertPoint( new Point(34,76) );
delaunay.insertPoint( new Point(56,54) );
delaunay.insertPoint( new Point(78,32) );
/* ... */
 
/* retourne la liste des segments */
List<Point[]> edges = delaunay.computeEdges();
 
/* retourne la liste des triangles */
List<Point[]> triangles = delaunay.computeTriangles();
 
/* retourne la liste des régions */
List<Point[]> voronoi = delaunay.computeVoronoi();


Edit : Modif de la méthode QuadEdge.IsOnLine(), cf. post #38 de Commander Keen.