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
31CDelaunay::CDelaunay()
32{
33}
34
35CDelaunay::~CDelaunay()
36{
37}
38
39// Allocate storage, construct triangulation, compute voronoi corners
40int 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
54void 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.
73void 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
130EdgePointer 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
146void 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
159EdgePointer 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
168EdgePointer 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
178void 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
190CSite *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
207void CDelaunay::freeMemory()
208{
209  if (sa) {
210    free(sa);
211    sa = (CSite*)NULL;
212  }
213}
214
215//
216// Edge storage management
217//
218
219void CDelaunay::deleteAllEdges()
220{
221    nextEdge = 0;
222    availEdge = NYL;
223}
224
225EdgePointer 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
237void CDelaunay::freeEdge(EdgePointer e)
238{
239  e ^= e & 3;
240  onext(e) = availEdge;
241  availEdge = e;
242}
243
244EdgePointer 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
270int 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
289int 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
308int CDelaunay::cmpev(int i, int j)
309{
310  return (ev[i].first - ev[j].first);
311}
312
313void 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
331void 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
340void 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
349void CDelaunay::copyev(int i, int j)
350{
351  ev[j] = ev[i];
352}
353
354void 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
360void 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
366void 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.
473int 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
493int 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
516void 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
590int 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
622void 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