1 /*
2 * Copyright (C) 2011 The Android Open Source Project
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 // Delaunay.cpp
18 // $Id: Delaunay.cpp,v 1.10 2011/06/17 13:35:48 mbansal Exp $
19
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <memory.h>
23 #include "Delaunay.h"
24
25 #define QQ 9 // Optimal value as determined by testing
26 #define DM 38 // 2^(1+DM/2) element sort capability. DM=38 for >10^6 elements
27 #define NYL (-1)
28 #define valid(l) ccw(orig(basel), dest(l), dest(basel))
29
30
CDelaunay()31 CDelaunay::CDelaunay()
32 {
33 }
34
~CDelaunay()35 CDelaunay::~CDelaunay()
36 {
37 }
38
39 // Allocate storage, construct triangulation, compute voronoi corners
triangulate(SEdgeVector ** edges,int n_sites,int width,int height)40 int CDelaunay::triangulate(SEdgeVector **edges, int n_sites, int width, int height)
41 {
42 EdgePointer cep;
43
44 deleteAllEdges();
45 buildTriangulation(n_sites);
46 cep = consolidateEdges();
47 *edges = ev;
48
49 // Note: construction_list will change ev
50 return constructList(cep, width, height);
51 }
52
53 // builds delaunay triangulation
buildTriangulation(int size)54 void CDelaunay::buildTriangulation(int size)
55 {
56 int i, rows;
57 EdgePointer lefte, righte;
58
59 rows = (int)( 0.5 + sqrt( (double) size / log( (double) size )));
60
61 // Sort the pointers by x-coordinate of site
62 for ( i=0 ; i < size ; i++ ) {
63 sp[i] = (SitePointer) i;
64 }
65
66 spsortx( sp, 0, size-1 );
67 build( 0, size-1, &lefte, &righte, rows );
68 oneBndryEdge = lefte;
69 }
70
71 // Recursive Delaunay Triangulation Procedure
72 // Contains modifications for axis-switching division.
build(int lo,int hi,EdgePointer * le,EdgePointer * re,int rows)73 void CDelaunay::build(int lo, int hi, EdgePointer *le, EdgePointer *re, int rows)
74 {
75 EdgePointer a, b, c, ldo, rdi, ldi, rdo, maxx, minx;
76 int split, lowrows;
77 int low, high;
78 SitePointer s1, s2, s3;
79 low = lo;
80 high = hi;
81
82 if ( low < (high-2) ) {
83 // more than three elements; do recursion
84 minx = sp[low];
85 maxx = sp[high];
86 if (rows == 1) { // time to switch axis of division
87 spsorty( sp, low, high);
88 rows = 65536;
89 }
90 lowrows = rows/2;
91 split = low - 1 + (int)
92 (0.5 + ((double)(high-low+1) * ((double)lowrows / (double)rows)));
93 build( low, split, &ldo, &ldi, lowrows );
94 build( split+1, high, &rdi, &rdo, (rows-lowrows) );
95 doMerge(&ldo, ldi, rdi, &rdo);
96 while (orig(ldo) != minx) {
97 ldo = rprev(ldo);
98 }
99 while (orig(rdo) != maxx) {
100 rdo = (SitePointer) lprev(rdo);
101 }
102 *le = ldo;
103 *re = rdo;
104 }
105 else if (low >= (high - 1)) { // two or one points
106 a = makeEdge(sp[low], sp[high]);
107 *le = a;
108 *re = (EdgePointer) sym(a);
109 } else { // three points
110 // 3 cases: triangles of 2 orientations, and 3 points on a line
111 a = makeEdge((s1 = sp[low]), (s2 = sp[low+1]));
112 b = makeEdge(s2, (s3 = sp[high]));
113 splice((EdgePointer) sym(a), b);
114 if (ccw(s1, s3, s2)) {
115 c = connectLeft(b, a);
116 *le = (EdgePointer) sym(c);
117 *re = c;
118 } else {
119 *le = a;
120 *re = (EdgePointer) sym(b);
121 if (ccw(s1, s2, s3)) {
122 // not colinear
123 c = connectLeft(b, a);
124 }
125 }
126 }
127 }
128
129 // Quad-edge manipulation primitives
makeEdge(SitePointer origin,SitePointer destination)130 EdgePointer CDelaunay::makeEdge(SitePointer origin, SitePointer destination)
131 {
132 EdgePointer temp, ans;
133 temp = allocEdge();
134 ans = temp;
135
136 onext(temp) = ans;
137 orig(temp) = origin;
138 onext(++temp) = (EdgePointer) (ans + 3);
139 onext(++temp) = (EdgePointer) (ans + 2);
140 orig(temp) = destination;
141 onext(++temp) = (EdgePointer) (ans + 1);
142
143 return(ans);
144 }
145
splice(EdgePointer a,EdgePointer b)146 void CDelaunay::splice(EdgePointer a, EdgePointer b)
147 {
148 EdgePointer alpha, beta, temp;
149 alpha = (EdgePointer) rot(onext(a));
150 beta = (EdgePointer) rot(onext(b));
151 temp = onext(alpha);
152 onext(alpha) = onext(beta);
153 onext(beta) = temp;
154 temp = onext(a);
155 onext(a) = onext(b);
156 onext(b) = temp;
157 }
158
connectLeft(EdgePointer a,EdgePointer b)159 EdgePointer CDelaunay::connectLeft(EdgePointer a, EdgePointer b)
160 {
161 EdgePointer ans;
162 ans = makeEdge(dest(a), orig(b));
163 splice(ans, (EdgePointer) lnext(a));
164 splice((EdgePointer) sym(ans), b);
165 return(ans);
166 }
167
connectRight(EdgePointer a,EdgePointer b)168 EdgePointer CDelaunay::connectRight(EdgePointer a, EdgePointer b)
169 {
170 EdgePointer ans;
171 ans = makeEdge(dest(a), orig(b));
172 splice(ans, (EdgePointer) sym(a));
173 splice((EdgePointer) sym(ans), (EdgePointer) oprev(b));
174 return(ans);
175 }
176
177 // disconnects e from the rest of the structure and destroys it
deleteEdge(EdgePointer e)178 void CDelaunay::deleteEdge(EdgePointer e)
179 {
180 splice(e, (EdgePointer) oprev(e));
181 splice((EdgePointer) sym(e), (EdgePointer) oprev(sym(e)));
182 freeEdge(e);
183 }
184
185 //
186 // Overall storage allocation
187 //
188
189 // Quad-edge storage allocation
allocMemory(int n)190 CSite *CDelaunay::allocMemory(int n)
191 {
192 unsigned int size;
193
194 size = ((sizeof(CSite) + sizeof(SitePointer)) * n +
195 (sizeof(SitePointer) + sizeof(EdgePointer)) * 12
196 ) * n;
197 if (!(sa = (CSite*) malloc(size))) {
198 return NULL;
199 }
200 sp = (SitePointer *) (sa + n);
201 ev = (SEdgeVector *) (org = sp + n);
202 next = (EdgePointer *) (org + 12 * n);
203 ei = (struct EDGE_INFO *) (next + 12 * n);
204 return sa;
205 }
206
freeMemory()207 void CDelaunay::freeMemory()
208 {
209 if (sa) {
210 free(sa);
211 sa = (CSite*)NULL;
212 }
213 }
214
215 //
216 // Edge storage management
217 //
218
deleteAllEdges()219 void CDelaunay::deleteAllEdges()
220 {
221 nextEdge = 0;
222 availEdge = NYL;
223 }
224
allocEdge()225 EdgePointer CDelaunay::allocEdge()
226 {
227 EdgePointer ans;
228
229 if (availEdge == NYL) {
230 ans = nextEdge, nextEdge += 4;
231 } else {
232 ans = availEdge, availEdge = onext(availEdge);
233 }
234 return(ans);
235 }
236
freeEdge(EdgePointer e)237 void CDelaunay::freeEdge(EdgePointer e)
238 {
239 e ^= e & 3;
240 onext(e) = availEdge;
241 availEdge = e;
242 }
243
consolidateEdges()244 EdgePointer CDelaunay::consolidateEdges()
245 {
246 EdgePointer e;
247 int i,j;
248
249 while (availEdge != NYL) {
250 nextEdge -= 4; e = availEdge; availEdge = onext(availEdge);
251
252 if (e==nextEdge) {
253 continue; // the one deleted was the last one anyway
254 }
255 if ((oneBndryEdge&~3) == nextEdge) {
256 oneBndryEdge = (EdgePointer) (e | (oneBndryEdge&3));
257 }
258 for (i=0,j=3; i<4; i++,j=rot(j)) {
259 onext(e+i) = onext(nextEdge+i);
260 onext(rot(onext(e+i))) = (EdgePointer) (e+j);
261 }
262 }
263 return nextEdge;
264 }
265
266 //
267 // Sorting Routines
268 //
269
xcmpsp(int i,int j)270 int CDelaunay::xcmpsp(int i, int j)
271 {
272 double d = sa[(i>=0)?sp[i]:sp1].X() - sa[(j>=0)?sp[j]:sp1].X();
273 if ( d > 0. ) {
274 return 1;
275 }
276 if ( d < 0. ) {
277 return -1;
278 }
279 d = sa[(i>=0)?sp[i]:sp1].Y() - sa[(j>=0)?sp[j]:sp1].Y();
280 if ( d > 0. ) {
281 return 1;
282 }
283 if ( d < 0. ) {
284 return -1;
285 }
286 return 0;
287 }
288
ycmpsp(int i,int j)289 int CDelaunay::ycmpsp(int i, int j)
290 {
291 double d = sa[(i>=0)?sp[i]:sp1].Y() - sa[(j>=0)?sp[j]:sp1].Y();
292 if ( d > 0. ) {
293 return 1;
294 }
295 if ( d < 0. ) {
296 return -1;
297 }
298 d = sa[(i>=0)?sp[i]:sp1].X() - sa[(j>=0)?sp[j]:sp1].X();
299 if ( d > 0. ) {
300 return 1;
301 }
302 if ( d < 0. ) {
303 return -1;
304 }
305 return 0;
306 }
307
cmpev(int i,int j)308 int CDelaunay::cmpev(int i, int j)
309 {
310 return (ev[i].first - ev[j].first);
311 }
312
swapsp(int i,int j)313 void CDelaunay::swapsp(int i, int j)
314 {
315 int t;
316 t = (i>=0) ? sp[i] : sp1;
317
318 if (i>=0) {
319 sp[i] = (j>=0)?sp[j]:sp1;
320 } else {
321 sp1 = (j>=0)?sp[j]:sp1;
322 }
323
324 if (j>=0) {
325 sp[j] = (SitePointer) t;
326 } else {
327 sp1 = (SitePointer) t;
328 }
329 }
330
swapev(int i,int j)331 void CDelaunay::swapev(int i, int j)
332 {
333 SEdgeVector temp;
334
335 temp = ev[i];
336 ev[i] = ev[j];
337 ev[j] = temp;
338 }
339
copysp(int i,int j)340 void CDelaunay::copysp(int i, int j)
341 {
342 if (j>=0) {
343 sp[j] = (i>=0)?sp[i]:sp1;
344 } else {
345 sp1 = (i>=0)?sp[i]:sp1;
346 }
347 }
348
copyev(int i,int j)349 void CDelaunay::copyev(int i, int j)
350 {
351 ev[j] = ev[i];
352 }
353
spsortx(SitePointer * sp_in,int low,int high)354 void CDelaunay::spsortx(SitePointer *sp_in, int low, int high)
355 {
356 sp = sp_in;
357 rcssort(low,high,-1,&CDelaunay::xcmpsp,&CDelaunay::swapsp,&CDelaunay::copysp);
358 }
359
spsorty(SitePointer * sp_in,int low,int high)360 void CDelaunay::spsorty(SitePointer *sp_in, int low, int high )
361 {
362 sp = sp_in;
363 rcssort(low,high,-1,&CDelaunay::ycmpsp,&CDelaunay::swapsp,&CDelaunay::copysp);
364 }
365
rcssort(int lowelt,int highelt,int temp,int (CDelaunay::* comparison)(int,int),void (CDelaunay::* swap)(int,int),void (CDelaunay::* copy)(int,int))366 void CDelaunay::rcssort(int lowelt, int highelt, int temp,
367 int (CDelaunay::*comparison)(int,int),
368 void (CDelaunay::*swap)(int,int),
369 void (CDelaunay::*copy)(int,int))
370 {
371 int m,sij,si,sj,sL,sk;
372 int stack[DM];
373
374 if (highelt-lowelt<=1) {
375 return;
376 }
377 if (highelt-lowelt>QQ) {
378 m = 0;
379 si = lowelt; sj = highelt;
380 for (;;) { // partition [si,sj] about median-of-3.
381 sij = (sj+si) >> 1;
382
383 // Now to sort elements si,sij,sj into order & set temp=their median
384 if ( (this->*comparison)( si,sij ) > 0 ) {
385 (this->*swap)( si,sij );
386 }
387 if ( (this->*comparison)( sij,sj ) > 0 ) {
388 (this->*swap)( sj,sij );
389 if ( (this->*comparison)( si,sij ) > 0 ) {
390 (this->*swap)( si,sij );
391 }
392 }
393 (this->*copy)( sij,temp );
394
395 // Now to partition into elements <=temp, >=temp, and ==temp.
396 sk = si; sL = sj;
397 do {
398 do {
399 sL--;
400 } while( (this->*comparison)( sL,temp ) > 0 );
401 do {
402 sk++;
403 } while( (this->*comparison)( temp,sk ) > 0 );
404 if ( sk < sL ) {
405 (this->*swap)( sL,sk );
406 }
407 } while(sk <= sL);
408
409 // Now to recurse on shorter partition, store longer partition on stack
410 if ( sL-si > sj-sk ) {
411 if ( sL-si < QQ ) {
412 if( m==0 ) {
413 break; // empty stack && both partitions < QQ so break
414 } else {
415 sj = stack[--m];
416 si = stack[--m];
417 }
418 }
419 else {
420 if ( sj-sk < QQ ) {
421 sj = sL;
422 } else {
423 stack[m++] = si;
424 stack[m++] = sL;
425 si = sk;
426 }
427 }
428 }
429 else {
430 if ( sj-sk < QQ ) {
431 if ( m==0 ) {
432 break; // empty stack && both partitions < QQ so break
433 } else {
434 sj = stack[--m];
435 si = stack[--m];
436 }
437 }
438 else {
439 if ( sL-si < QQ ) {
440 si = sk;
441 } else {
442 stack[m++] = sk;
443 stack[m++] = sj;
444 sj = sL;
445 }
446 }
447 }
448 }
449 }
450
451 // Now for 0 or Data bounded "straight insertion" sort of [0,nels-1]; if it is
452 // known that el[-1] = -INF, then can omit the "sk>=0" test and save time.
453 for (si=lowelt; si<highelt; si++) {
454 if ( (this->*comparison)( si,si+1 ) > 0 ) {
455 (this->*copy)( si+1,temp );
456 sj = sk = si;
457 sj++;
458 do {
459 (this->*copy)( sk,sj );
460 sj = sk;
461 sk--;
462 } while ( (this->*comparison)( sk,temp ) > 0 && sk>=lowelt );
463 (this->*copy)( temp,sj );
464 }
465 }
466 }
467
468 //
469 // Geometric primitives
470 //
471
472 // incircle, as in the Guibas-Stolfi paper.
incircle(SitePointer a,SitePointer b,SitePointer c,SitePointer d)473 int CDelaunay::incircle(SitePointer a, SitePointer b, SitePointer c, SitePointer d)
474 {
475 double adx, ady, bdx, bdy, cdx, cdy, dx, dy, nad, nbd, ncd;
476 dx = sa[d].X();
477 dy = sa[d].Y();
478 adx = sa[a].X() - dx;
479 ady = sa[a].Y() - dy;
480 bdx = sa[b].X() - dx;
481 bdy = sa[b].Y() - dy;
482 cdx = sa[c].X() - dx;
483 cdy = sa[c].Y() - dy;
484 nad = adx*adx+ady*ady;
485 nbd = bdx*bdx+bdy*bdy;
486 ncd = cdx*cdx+cdy*cdy;
487 return( (0.0 < (nad * (bdx * cdy - bdy * cdx)
488 + nbd * (cdx * ady - cdy * adx)
489 + ncd * (adx * bdy - ady * bdx))) ? TRUE : FALSE );
490 }
491
492 // TRUE iff A, B, C form a counterclockwise oriented triangle
ccw(SitePointer a,SitePointer b,SitePointer c)493 int CDelaunay::ccw(SitePointer a, SitePointer b, SitePointer c)
494 {
495 int result;
496
497 double ax = sa[a].X();
498 double bx = sa[b].X();
499 double cx = sa[c].X();
500 double ay = sa[a].Y();
501 double by = sa[b].Y();
502 double cy = sa[c].Y();
503
504 double val = (ax - cx)*(by - cy) - (bx - cx)*(ay - cy);
505 if ( val > 0.0) {
506 return true;
507 }
508
509 return false;
510 }
511
512 //
513 // The Merge Procedure.
514 //
515
doMerge(EdgePointer * ldo,EdgePointer ldi,EdgePointer rdi,EdgePointer * rdo)516 void CDelaunay::doMerge(EdgePointer *ldo, EdgePointer ldi, EdgePointer rdi, EdgePointer *rdo)
517 {
518 int rvalid, lvalid;
519 EdgePointer basel,lcand,rcand,t;
520
521 for (;;) {
522 while (ccw(orig(ldi), dest(ldi), orig(rdi))) {
523 ldi = (EdgePointer) lnext(ldi);
524 }
525 if (ccw(dest(rdi), orig(rdi), orig(ldi))) {
526 rdi = (EdgePointer)rprev(rdi);
527 } else {
528 break;
529 }
530 }
531
532 basel = connectLeft((EdgePointer) sym(rdi), ldi);
533 lcand = rprev(basel);
534 rcand = (EdgePointer) oprev(basel);
535 if (orig(basel) == orig(*rdo)) {
536 *rdo = basel;
537 }
538 if (dest(basel) == orig(*ldo)) {
539 *ldo = (EdgePointer) sym(basel);
540 }
541
542 for (;;) {
543 #if 1
544 if (valid(t=onext(lcand))) {
545 #else
546 t = (EdgePointer)onext(lcand);
547 if (valid(basel, t)) {
548 #endif
549 while (incircle(dest(lcand), dest(t), orig(lcand), orig(basel))) {
550 deleteEdge(lcand);
551 lcand = t;
552 t = onext(lcand);
553 }
554 }
555 #if 1
556 if (valid(t=(EdgePointer)oprev(rcand))) {
557 #else
558 t = (EdgePointer)oprev(rcand);
559 if (valid(basel, t)) {
560 #endif
561 while (incircle(dest(t), dest(rcand), orig(rcand), dest(basel))) {
562 deleteEdge(rcand);
563 rcand = t;
564 t = (EdgePointer)oprev(rcand);
565 }
566 }
567
568 #if 1
569 lvalid = valid(lcand);
570 rvalid = valid(rcand);
571 #else
572 lvalid = valid(basel, lcand);
573 rvalid = valid(basel, rcand);
574 #endif
575 if ((! lvalid) && (! rvalid)) {
576 return;
577 }
578
579 if (!lvalid ||
580 (rvalid && incircle(dest(lcand), orig(lcand), orig(rcand), dest(rcand)))) {
581 basel = connectLeft(rcand, (EdgePointer) sym(basel));
582 rcand = (EdgePointer) lnext(sym(basel));
583 } else {
584 basel = (EdgePointer) sym(connectRight(lcand, basel));
585 lcand = rprev(basel);
586 }
587 }
588 }
589
590 int CDelaunay::constructList(EdgePointer last, int width, int height)
591 {
592 int c, i;
593 EdgePointer curr, src, nex;
594 SEdgeVector *currv, *prevv;
595
596 c = (int) ((curr = (EdgePointer) ((last & ~3))) >> 1);
597
598 for (last -= 4; last >= 0; last -= 4) {
599 src = orig(last);
600 nex = dest(last);
601 orig(--curr) = src;
602 orig(--curr) = nex;
603 orig(--curr) = nex;
604 orig(--curr) = src;
605 }
606 rcssort(0, c - 1, -1, &CDelaunay::cmpev, &CDelaunay::swapev, &CDelaunay::copyev);
607
608 // Throw out any edges that are too far apart
609 currv = prevv = ev;
610 for (i = c; i--; currv++) {
611 if ((int) fabs(sa[currv->first].getVCenter().x - sa[currv->second].getVCenter().x) <= width &&
612 (int) fabs(sa[currv->first].getVCenter().y - sa[currv->second].getVCenter().y) <= height) {
613 *(prevv++) = *currv;
614 } else {
615 c--;
616 }
617 }
618 return c;
619 }
620
621 // Fill in site neighbor information
622 void CDelaunay::linkNeighbors(SEdgeVector *edge, int nedge, int nsite)
623 {
624 int i;
625
626 for (i = 0; i < nsite; i++) {
627 sa[i].setNeighbor(edge);
628 sa[i].setNumNeighbors(0);
629 for (; edge->first == i && nedge; edge++, nedge--) {
630 sa[i].incrNumNeighbors();
631 }
632 }
633 }
634