Home | History | Annotate | Download | only in mosaic
      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 
     31 CDelaunay::CDelaunay()
     32 {
     33 }
     34 
     35 CDelaunay::~CDelaunay()
     36 {
     37 }
     38 
     39 // Allocate storage, construct triangulation, compute voronoi corners
     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
     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.
     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
    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 
    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 
    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 
    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
    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
    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 
    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 
    219 void CDelaunay::deleteAllEdges()
    220 {
    221     nextEdge = 0;
    222     availEdge = NYL;
    223 }
    224 
    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 
    237 void CDelaunay::freeEdge(EdgePointer e)
    238 {
    239   e ^= e & 3;
    240   onext(e) = availEdge;
    241   availEdge = e;
    242 }
    243 
    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 
    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 
    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 
    308 int CDelaunay::cmpev(int i, int j)
    309 {
    310   return (ev[i].first - ev[j].first);
    311 }
    312 
    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 
    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 
    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 
    349 void CDelaunay::copyev(int i, int j)
    350 {
    351   ev[j] = ev[i];
    352 }
    353 
    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 
    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 
    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.
    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
    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 
    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