← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/delaunay/VoronoiDiagramGenerator.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sandbox/delaunay/VoronoiDiagramGenerator.cpp
1
 
/*
2
 
 * The author of this software is Steven Fortune.  Copyright (c) 1994 by AT&T
3
 
 * Bell Laboratories.
4
 
 * Permission to use, copy, modify, and distribute this software for any
5
 
 * purpose without fee is hereby granted, provided that this entire notice
6
 
 * is included in all copies of any software which is or includes a copy
7
 
 * or modification of this software and in all copies of the supporting
8
 
 * documentation for such software.
9
 
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
10
 
 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
11
 
 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
12
 
 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
13
 
 */
14
 
 
15
 
/* 
16
 
 * This code was originally written by Stephan Fortune in C code.  Shane O'Sullivan, 
17
 
 * have since modified it, encapsulating it in a C++ class and, fixing memory leaks and 
18
 
 * adding accessors to the Voronoi Edges.
19
 
 * Permission to use, copy, modify, and distribute this software for any
20
 
 * purpose without fee is hereby granted, provided that this entire notice
21
 
 * is included in all copies of any software which is or includes a copy
22
 
 * or modification of this software and in all copies of the supporting
23
 
 * documentation for such software.
24
 
 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
25
 
 * WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
26
 
 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
27
 
 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
28
 
 */
29
 
 
30
 
/*
31
 
 * Subsequently, Robert Kern modified it to yield Python objects.
32
 
 * Copyright 2005 Robert Kern <robert.kern@gmail.com>
33
 
 * See LICENSE.txt in the scipy source directory.
34
 
 */
35
 
 
36
 
#include "VoronoiDiagramGenerator.h"
37
 
 
38
 
VoronoiDiagramGenerator::VoronoiDiagramGenerator()
39
 
{
40
 
    siteidx = 0;
41
 
    sites = 0;
42
 
 
43
 
    allMemoryList = new FreeNodeArrayList;
44
 
    allMemoryList->memory = 0;
45
 
    allMemoryList->next = 0;
46
 
    currentMemoryBlock = allMemoryList;
47
 
    allEdges = 0;
48
 
    allEdgeList = 0;
49
 
    iteratorEdges = 0;
50
 
    iterEdgeList = 0;
51
 
    minDistanceBetweenSites = 0;
52
 
}
53
 
 
54
 
VoronoiDiagramGenerator::~VoronoiDiagramGenerator()
55
 
{
56
 
    cleanupEdgeList();
57
 
    cleanup();
58
 
    cleanupEdges();
59
 
 
60
 
    if(allMemoryList != 0)
61
 
        delete allMemoryList;
62
 
}
63
 
 
64
 
 
65
 
 
66
 
bool VoronoiDiagramGenerator::generateVoronoi(double *xValues, double *yValues, int numPoints, double minX, double maxX, double minY, double maxY, double minDist)
67
 
{
68
 
    cleanupEdgeList();
69
 
    cleanup();
70
 
    cleanupEdges();
71
 
    int i;
72
 
 
73
 
    minDistanceBetweenSites = minDist;
74
 
 
75
 
    nsites=numPoints;
76
 
    plot = 0;
77
 
    triangulate = 0;    
78
 
    debug = 1;
79
 
    sorted = 0; 
80
 
    freeinit(&sfl, sizeof (Site));
81
 
 
82
 
    sites = (struct Site *) myalloc(nsites*sizeof( *sites));
83
 
 
84
 
    if(sites == 0)
85
 
        return false;
86
 
 
87
 
    xmin = xValues[0];
88
 
    ymin = yValues[0];
89
 
    xmax = xValues[0];
90
 
    ymax = yValues[0];
91
 
 
92
 
    for(i = 0; i< nsites; i++)
93
 
    {
94
 
        sites[i].coord.x = xValues[i];
95
 
        sites[i].coord.y = yValues[i];
96
 
        sites[i].sitenbr = i;
97
 
        sites[i].refcnt = 0;
98
 
 
99
 
        if(xValues[i] < xmin)
100
 
            xmin = xValues[i];
101
 
        else if(xValues[i] > xmax)
102
 
            xmax = xValues[i];
103
 
 
104
 
        if(yValues[i] < ymin)
105
 
            ymin = yValues[i];
106
 
        else if(yValues[i] > ymax)
107
 
            ymax = yValues[i];
108
 
 
109
 
        //printf("\n%f %f\n",xValues[i],yValues[i]);
110
 
    }
111
 
    
112
 
    qsort(sites, nsites, sizeof (*sites), scomp);
113
 
    
114
 
    siteidx = 0;
115
 
    geominit();
116
 
    double temp = 0;
117
 
    if(minX > maxX)
118
 
    {
119
 
        temp = minX;
120
 
        minX = maxX;
121
 
        maxX = temp;
122
 
    }
123
 
    if(minY > maxY)
124
 
    {
125
 
        temp = minY;
126
 
        minY = maxY;
127
 
        maxY = temp;
128
 
    }
129
 
    borderMinX = minX;
130
 
    borderMinY = minY;
131
 
    borderMaxX = maxX;
132
 
    borderMaxY = maxY;
133
 
    
134
 
    siteidx = 0;    
135
 
 
136
 
    voronoi(triangulate);
137
 
 
138
 
    return true;
139
 
}
140
 
 
141
 
bool VoronoiDiagramGenerator::ELinitialize()
142
 
{
143
 
    int i;
144
 
    freeinit(&hfl, sizeof **ELhash);
145
 
    ELhashsize = 2 * sqrt_nsites;
146
 
    ELhash = (struct Halfedge **) myalloc ( sizeof *ELhash * ELhashsize);
147
 
 
148
 
    if(ELhash == 0)
149
 
        return false;
150
 
 
151
 
    for(i=0; i<ELhashsize; i +=1) ELhash[i] = (struct Halfedge *)NULL;
152
 
    ELleftend = HEcreate( (struct Edge *)NULL, 0);
153
 
    ELrightend = HEcreate( (struct Edge *)NULL, 0);
154
 
    ELleftend -> ELleft = (struct Halfedge *)NULL;
155
 
    ELleftend -> ELright = ELrightend;
156
 
    ELrightend -> ELleft = ELleftend;
157
 
    ELrightend -> ELright = (struct Halfedge *)NULL;
158
 
    ELhash[0] = ELleftend;
159
 
    ELhash[ELhashsize-1] = ELrightend;
160
 
 
161
 
    return true;
162
 
}
163
 
 
164
 
 
165
 
struct Halfedge* VoronoiDiagramGenerator::HEcreate(struct Edge *e,int pm)
166
 
{
167
 
    struct Halfedge *answer;
168
 
    answer = (struct Halfedge *) getfree(&hfl);
169
 
    answer -> ELedge = e;
170
 
    answer -> ELpm = pm;
171
 
    answer -> PQnext = (struct Halfedge *) NULL;
172
 
    answer -> vertex = (struct Site *) NULL;
173
 
    answer -> ELrefcnt = 0;
174
 
    return(answer);
175
 
}
176
 
 
177
 
 
178
 
void VoronoiDiagramGenerator::ELinsert(struct    Halfedge *lb, struct Halfedge *newHe)
179
 
{
180
 
    newHe -> ELleft = lb;
181
 
    newHe -> ELright = lb -> ELright;
182
 
    (lb -> ELright) -> ELleft = newHe;
183
 
    lb -> ELright = newHe;
184
 
}
185
 
 
186
 
/* Get entry from hash table, pruning any deleted nodes */
187
 
struct Halfedge * VoronoiDiagramGenerator::ELgethash(int b)
188
 
{
189
 
    struct Halfedge *he;
190
 
    
191
 
    if(b<0 || b>=ELhashsize) 
192
 
        return((struct Halfedge *) NULL);
193
 
    he = ELhash[b]; 
194
 
    if (he == (struct Halfedge *) NULL || he->ELedge != (struct Edge *) DELETED ) 
195
 
        return (he);
196
 
    
197
 
    /* Hash table points to deleted half edge.  Patch as necessary. */
198
 
    ELhash[b] = (struct Halfedge *) NULL;
199
 
    if ((he -> ELrefcnt -= 1) == 0) 
200
 
        makefree((Freenode*)he, &hfl);
201
 
    return ((struct Halfedge *) NULL);
202
 
}    
203
 
 
204
 
struct Halfedge * VoronoiDiagramGenerator::ELleftbnd(struct Point *p)
205
 
{
206
 
    int i, bucket;
207
 
    struct Halfedge *he;
208
 
    
209
 
    /* Use hash table to get close to desired halfedge */
210
 
    bucket = (int)((p->x - xmin)/deltax * ELhashsize);    //use the hash function to find the place in the hash map that this HalfEdge should be
211
 
 
212
 
    if(bucket<0) bucket =0;                    //make sure that the bucket position in within the range of the hash array
213
 
    if(bucket>=ELhashsize) bucket = ELhashsize - 1;
214
 
 
215
 
    he = ELgethash(bucket);
216
 
    if(he == (struct Halfedge *) NULL)            //if the HE isn't found, search backwards and forwards in the hash map for the first non-null entry
217
 
    {   
218
 
        for(i=1; 1 ; i += 1)
219
 
        {    
220
 
            if ((he=ELgethash(bucket-i)) != (struct Halfedge *) NULL) 
221
 
                break;
222
 
            if ((he=ELgethash(bucket+i)) != (struct Halfedge *) NULL) 
223
 
                break;
224
 
        };
225
 
        totalsearch += i;
226
 
    };
227
 
    ntry += 1;
228
 
    /* Now search linear list of halfedges for the correct one */
229
 
    if (he==ELleftend  || (he != ELrightend && right_of(he,p)))
230
 
    {
231
 
        do 
232
 
        {
233
 
            he = he -> ELright;
234
 
        } while (he!=ELrightend && right_of(he,p));    //keep going right on the list until either the end is reached, or you find the 1st edge which the point
235
 
        he = he -> ELleft;                //isn't to the right of
236
 
    }
237
 
    else                             //if the point is to the left of the HalfEdge, then search left for the HE just to the left of the point
238
 
        do 
239
 
        {
240
 
            he = he -> ELleft;
241
 
        } while (he!=ELleftend && !right_of(he,p));
242
 
        
243
 
    /* Update hash table and reference counts */
244
 
    if(bucket > 0 && bucket <ELhashsize-1)
245
 
    {    
246
 
        if(ELhash[bucket] != (struct Halfedge *) NULL) 
247
 
        {
248
 
            ELhash[bucket] -> ELrefcnt -= 1;
249
 
        }
250
 
        ELhash[bucket] = he;
251
 
        ELhash[bucket] -> ELrefcnt += 1;
252
 
    };
253
 
    return (he);
254
 
}
255
 
 
256
 
 
257
 
/* This delete routine can't reclaim node, since pointers from hash
258
 
table may be present.   */
259
 
void VoronoiDiagramGenerator::ELdelete(struct Halfedge *he)
260
 
{
261
 
    (he -> ELleft) -> ELright = he -> ELright;
262
 
    (he -> ELright) -> ELleft = he -> ELleft;
263
 
    he -> ELedge = (struct Edge *)DELETED;
264
 
}
265
 
 
266
 
 
267
 
struct Halfedge * VoronoiDiagramGenerator::ELright(struct Halfedge *he)
268
 
{
269
 
    return (he -> ELright);
270
 
}
271
 
 
272
 
struct Halfedge * VoronoiDiagramGenerator::ELleft(struct Halfedge *he)
273
 
{
274
 
    return (he -> ELleft);
275
 
}
276
 
 
277
 
 
278
 
struct Site * VoronoiDiagramGenerator::leftreg(struct Halfedge *he)
279
 
{
280
 
    if(he -> ELedge == (struct Edge *)NULL) 
281
 
        return(bottomsite);
282
 
    return( he -> ELpm == le ? 
283
 
        he -> ELedge -> reg[le] : he -> ELedge -> reg[re]);
284
 
}
285
 
 
286
 
struct Site * VoronoiDiagramGenerator::rightreg(struct Halfedge *he)
287
 
{
288
 
    if(he -> ELedge == (struct Edge *)NULL) //if this halfedge has no edge, return the bottom site (whatever that is)
289
 
        return(bottomsite);
290
 
 
291
 
    //if the ELpm field is zero, return the site 0 that this edge bisects, otherwise return site number 1
292
 
    return( he -> ELpm == le ? he -> ELedge -> reg[re] : he -> ELedge -> reg[le]);
293
 
}
294
 
 
295
 
void VoronoiDiagramGenerator::geominit()
296
 
{    
297
 
    double sn;
298
 
 
299
 
    freeinit(&efl, sizeof(Edge));
300
 
    nvertices = 0;
301
 
    nedges = 0;
302
 
    sn = (double)nsites+4;
303
 
    sqrt_nsites = (int)sqrt(sn);
304
 
    deltay = ymax - ymin;
305
 
    deltax = xmax - xmin;
306
 
}
307
 
 
308
 
 
309
 
struct Edge * VoronoiDiagramGenerator::bisect(struct Site *s1, struct Site *s2)
310
 
{
311
 
    double dx,dy,adx,ady;
312
 
    struct Edge *newedge;    
313
 
 
314
 
    newedge = (struct Edge *) getfree(&efl);
315
 
    
316
 
    newedge -> reg[0] = s1; //store the sites that this edge is bisecting
317
 
    newedge -> reg[1] = s2;
318
 
    ref(s1); 
319
 
    ref(s2);
320
 
    newedge -> ep[0] = (struct Site *) NULL; //to begin with, there are no endpoints on the bisector - it goes to infinity
321
 
    newedge -> ep[1] = (struct Site *) NULL;
322
 
    
323
 
    dx = s2->coord.x - s1->coord.x;            //get the difference in x dist between the sites
324
 
    dy = s2->coord.y - s1->coord.y;
325
 
    adx = dx>0 ? dx : -dx;                    //make sure that the difference in positive
326
 
    ady = dy>0 ? dy : -dy;
327
 
    newedge -> c = (double)(s1->coord.x * dx + s1->coord.y * dy + (dx*dx + dy*dy)*0.5);//get the slope of the line
328
 
 
329
 
    if (adx>ady)
330
 
    {    
331
 
        newedge -> a = 1.0; newedge -> b = dy/dx; newedge -> c /= dx;//set formula of line, with x fixed to 1
332
 
    }
333
 
    else
334
 
    {    
335
 
        newedge -> b = 1.0; newedge -> a = dx/dy; newedge -> c /= dy;//set formula of line, with y fixed to 1
336
 
    };
337
 
    
338
 
    newedge -> edgenbr = nedges;
339
 
 
340
 
    //printf("\nbisect(%d) ((%f,%f) and (%f,%f)",nedges,s1->coord.x,s1->coord.y,s2->coord.x,s2->coord.y);
341
 
    
342
 
    nedges += 1;
343
 
    return(newedge);
344
 
}
345
 
 
346
 
//create a new site where the HalfEdges el1 and el2 intersect - note that the Point in the argument list is not used, don't know why it's there
347
 
struct Site * VoronoiDiagramGenerator::intersect(struct Halfedge *el1, struct Halfedge *el2, struct Point *p)
348
 
{
349
 
    struct    Edge *e1,*e2, *e;
350
 
    struct  Halfedge *el;
351
 
    double d, xint, yint;
352
 
    int right_of_site;
353
 
    struct Site *v;
354
 
    
355
 
    e1 = el1 -> ELedge;
356
 
    e2 = el2 -> ELedge;
357
 
    if(e1 == (struct Edge*)NULL || e2 == (struct Edge*)NULL) 
358
 
        return ((struct Site *) NULL);
359
 
 
360
 
    //if the two edges bisect the same parent, return null
361
 
    if (e1->reg[1] == e2->reg[1]) 
362
 
        return ((struct Site *) NULL);
363
 
    
364
 
    d = e1->a * e2->b - e1->b * e2->a;
365
 
    if (-1.0e-10<d && d<1.0e-10) 
366
 
        return ((struct Site *) NULL);
367
 
    
368
 
    xint = (e1->c*e2->b - e2->c*e1->b)/d;
369
 
    yint = (e2->c*e1->a - e1->c*e2->a)/d;
370
 
    
371
 
    if( (e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
372
 
        (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
373
 
        e1->reg[1]->coord.x < e2->reg[1]->coord.x) )
374
 
    {    
375
 
        el = el1; 
376
 
        e = e1;
377
 
    }
378
 
    else
379
 
    {    
380
 
        el = el2; 
381
 
        e = e2;
382
 
    };
383
 
    
384
 
    right_of_site = xint >= e -> reg[1] -> coord.x;
385
 
    if ((right_of_site && el -> ELpm == le) || (!right_of_site && el -> ELpm == re)) 
386
 
        return ((struct Site *) NULL);
387
 
    
388
 
    //create a new site at the point of intersection - this is a new vector event waiting to happen
389
 
    v = (struct Site *) getfree(&sfl);
390
 
    v -> refcnt = 0;
391
 
    v -> coord.x = xint;
392
 
    v -> coord.y = yint;
393
 
    return(v);
394
 
}
395
 
 
396
 
/* returns 1 if p is to right of halfedge e */
397
 
int VoronoiDiagramGenerator::right_of(struct Halfedge *el,struct Point *p)
398
 
{
399
 
    struct Edge *e;
400
 
    struct Site *topsite;
401
 
    int right_of_site, above, fast;
402
 
    double dxp, dyp, dxs, t1, t2, t3, yl;
403
 
    
404
 
    e = el -> ELedge;
405
 
    topsite = e -> reg[1];
406
 
    right_of_site = p -> x > topsite -> coord.x;
407
 
    if(right_of_site && el -> ELpm == le) return(1);
408
 
    if(!right_of_site && el -> ELpm == re) return (0);
409
 
    
410
 
    if (e->a == 1.0)
411
 
    {    dyp = p->y - topsite->coord.y;
412
 
    dxp = p->x - topsite->coord.x;
413
 
    fast = 0;
414
 
    if ((!right_of_site & (e->b<0.0)) | (right_of_site & (e->b>=0.0)) )
415
 
    {    above = dyp>= e->b*dxp;    
416
 
    fast = above;
417
 
    }
418
 
    else 
419
 
    {    above = p->x + p->y*e->b > e-> c;
420
 
    if(e->b<0.0) above = !above;
421
 
    if (!above) fast = 1;
422
 
    };
423
 
    if (!fast)
424
 
    {    dxs = topsite->coord.x - (e->reg[0])->coord.x;
425
 
    above = e->b * (dxp*dxp - dyp*dyp) <
426
 
        dxs*dyp*(1.0+2.0*dxp/dxs + e->b*e->b);
427
 
    if(e->b<0.0) above = !above;
428
 
    };
429
 
    }
430
 
    else  /*e->b==1.0 */
431
 
    {    yl = e->c - e->a*p->x;
432
 
    t1 = p->y - yl;
433
 
    t2 = p->x - topsite->coord.x;
434
 
    t3 = yl - topsite->coord.y;
435
 
    above = t1*t1 > t2*t2 + t3*t3;
436
 
    };
437
 
    return (el->ELpm==le ? above : !above);
438
 
}
439
 
 
440
 
 
441
 
void VoronoiDiagramGenerator::endpoint(struct Edge *e,int lr,struct Site * s)
442
 
{
443
 
    e -> ep[lr] = s;
444
 
    ref(s);
445
 
    if(e -> ep[re-lr]== (struct Site *) NULL) 
446
 
        return;
447
 
 
448
 
    clip_line(e);
449
 
 
450
 
    deref(e->reg[le]);
451
 
    deref(e->reg[re]);
452
 
    makefree((Freenode*)e, &efl);
453
 
}
454
 
 
455
 
 
456
 
double VoronoiDiagramGenerator::dist(struct Site *s,struct Site *t)
457
 
{
458
 
    double dx,dy;
459
 
    dx = s->coord.x - t->coord.x;
460
 
    dy = s->coord.y - t->coord.y;
461
 
    return (double)(sqrt(dx*dx + dy*dy));
462
 
}
463
 
 
464
 
 
465
 
void VoronoiDiagramGenerator::makevertex(struct Site *v)
466
 
{
467
 
    v -> sitenbr = nvertices;
468
 
    nvertices += 1;
469
 
    out_vertex(v);
470
 
}
471
 
 
472
 
 
473
 
void VoronoiDiagramGenerator::deref(struct Site *v)
474
 
{
475
 
    v -> refcnt -= 1;
476
 
    if (v -> refcnt == 0 ) 
477
 
        makefree((Freenode*)v, &sfl);
478
 
}
479
 
 
480
 
void VoronoiDiagramGenerator::ref(struct Site *v)
481
 
{
482
 
    v -> refcnt += 1;
483
 
}
484
 
 
485
 
//push the HalfEdge into the ordered linked list of vertices
486
 
void VoronoiDiagramGenerator::PQinsert(struct Halfedge *he,struct Site * v, double offset)
487
 
{
488
 
    struct Halfedge *last, *next;
489
 
    
490
 
    he -> vertex = v;
491
 
    ref(v);
492
 
    he -> ystar = (double)(v -> coord.y + offset);
493
 
    last = &PQhash[PQbucket(he)];
494
 
    while ((next = last -> PQnext) != (struct Halfedge *) NULL &&
495
 
        (he -> ystar  > next -> ystar  ||
496
 
        (he -> ystar == next -> ystar && v -> coord.x > next->vertex->coord.x)))
497
 
    {    
498
 
        last = next;
499
 
    };
500
 
    he -> PQnext = last -> PQnext; 
501
 
    last -> PQnext = he;
502
 
    PQcount += 1;
503
 
}
504
 
 
505
 
//remove the HalfEdge from the list of vertices 
506
 
void VoronoiDiagramGenerator::PQdelete(struct Halfedge *he)
507
 
{
508
 
    struct Halfedge *last;
509
 
    
510
 
    if(he -> vertex != (struct Site *) NULL)
511
 
    {    
512
 
        last = &PQhash[PQbucket(he)];
513
 
        while (last -> PQnext != he) 
514
 
            last = last -> PQnext;
515
 
 
516
 
        last -> PQnext = he -> PQnext;
517
 
        PQcount -= 1;
518
 
        deref(he -> vertex);
519
 
        he -> vertex = (struct Site *) NULL;
520
 
    };
521
 
}
522
 
 
523
 
int VoronoiDiagramGenerator::PQbucket(struct Halfedge *he)
524
 
{
525
 
    int bucket;
526
 
    
527
 
    bucket = (int)((he->ystar - ymin)/deltay * PQhashsize);
528
 
    if (bucket<0) bucket = 0;
529
 
    if (bucket>=PQhashsize) bucket = PQhashsize-1 ;
530
 
    if (bucket < PQmin) PQmin = bucket;
531
 
    return(bucket);
532
 
}
533
 
 
534
 
 
535
 
 
536
 
int VoronoiDiagramGenerator::PQempty()
537
 
{
538
 
    return(PQcount==0);
539
 
}
540
 
 
541
 
 
542
 
struct Point VoronoiDiagramGenerator::PQ_min()
543
 
{
544
 
    struct Point answer;
545
 
    
546
 
    while(PQhash[PQmin].PQnext == (struct Halfedge *)NULL) {PQmin += 1;};
547
 
    answer.x = PQhash[PQmin].PQnext -> vertex -> coord.x;
548
 
    answer.y = PQhash[PQmin].PQnext -> ystar;
549
 
    return (answer);
550
 
}
551
 
 
552
 
struct Halfedge * VoronoiDiagramGenerator::PQextractmin()
553
 
{
554
 
    struct Halfedge *curr;
555
 
    
556
 
    curr = PQhash[PQmin].PQnext;
557
 
    PQhash[PQmin].PQnext = curr -> PQnext;
558
 
    PQcount -= 1;
559
 
    return(curr);
560
 
}
561
 
 
562
 
 
563
 
bool VoronoiDiagramGenerator::PQinitialize()
564
 
{
565
 
    int i; 
566
 
    
567
 
    PQcount = 0;
568
 
    PQmin = 0;
569
 
    PQhashsize = 4 * sqrt_nsites;
570
 
    PQhash = (struct Halfedge *) myalloc(PQhashsize * sizeof *PQhash);
571
 
 
572
 
    if(PQhash == 0)
573
 
        return false;
574
 
 
575
 
    for(i=0; i<PQhashsize; i+=1) PQhash[i].PQnext = (struct Halfedge *)NULL;
576
 
 
577
 
    return true;
578
 
}
579
 
 
580
 
 
581
 
void VoronoiDiagramGenerator::freeinit(struct Freelist *fl,int size)
582
 
{
583
 
    fl -> head = (struct Freenode *) NULL;
584
 
    fl -> nodesize = size;
585
 
}
586
 
 
587
 
char * VoronoiDiagramGenerator::getfree(struct Freelist *fl)
588
 
{
589
 
    int i; 
590
 
    struct Freenode *t;
591
 
 
592
 
    if(fl->head == (struct Freenode *) NULL)
593
 
    {    
594
 
        t =  (struct Freenode *) myalloc(sqrt_nsites * fl->nodesize);
595
 
 
596
 
        if(t == 0)
597
 
            return 0;
598
 
        
599
 
        currentMemoryBlock->next = new FreeNodeArrayList;
600
 
        currentMemoryBlock = currentMemoryBlock->next;
601
 
        currentMemoryBlock->memory = t;
602
 
        currentMemoryBlock->next = 0;
603
 
 
604
 
        for(i=0; i<sqrt_nsites; i+=1)     
605
 
            makefree((struct Freenode *)((char *)t+i*fl->nodesize), fl);        
606
 
    };
607
 
    t = fl -> head;
608
 
    fl -> head = (fl -> head) -> nextfree;
609
 
    return((char *)t);
610
 
}
611
 
 
612
 
 
613
 
 
614
 
void VoronoiDiagramGenerator::makefree(struct Freenode *curr,struct Freelist *fl)
615
 
{
616
 
    curr -> nextfree = fl -> head;
617
 
    fl -> head = curr;
618
 
}
619
 
 
620
 
void VoronoiDiagramGenerator::cleanup()
621
 
{
622
 
    if(sites != 0)
623
 
    {
624
 
        free(sites);
625
 
        sites = 0;
626
 
    }
627
 
 
628
 
    FreeNodeArrayList* current=0, *prev = 0;
629
 
 
630
 
    current = prev = allMemoryList;
631
 
 
632
 
    while(current->next != 0)
633
 
    {
634
 
        prev = current;
635
 
        current = current->next;
636
 
        free(prev->memory);
637
 
        delete prev;
638
 
        prev = 0;
639
 
    }
640
 
 
641
 
    if(current != 0 && current->memory != 0)
642
 
    {
643
 
        free(current->memory);
644
 
        delete current;
645
 
    }
646
 
 
647
 
    allMemoryList = new FreeNodeArrayList;
648
 
    allMemoryList->next = 0;
649
 
    allMemoryList->memory = 0;
650
 
    currentMemoryBlock = allMemoryList;
651
 
}
652
 
 
653
 
void VoronoiDiagramGenerator::cleanupEdges()
654
 
{
655
 
    GraphEdge* geCurrent = 0, *gePrev = 0;
656
 
    geCurrent = gePrev = allEdges;
657
 
 
658
 
    while(geCurrent != 0 && geCurrent->next != 0)
659
 
    {
660
 
        gePrev = geCurrent;
661
 
        geCurrent = geCurrent->next;
662
 
        delete gePrev;
663
 
    }
664
 
 
665
 
    allEdges = 0;
666
 
 
667
 
}
668
 
 
669
 
void VoronoiDiagramGenerator::cleanupEdgeList()
670
 
{
671
 
    EdgeList* elCurrent = 0, *elPrev = 0;
672
 
    elCurrent = elPrev = allEdgeList;
673
 
 
674
 
    while (elCurrent != 0 && elCurrent->next != 0)
675
 
    {
676
 
        elPrev = elCurrent;
677
 
        elCurrent = elCurrent->next;
678
 
        delete elPrev;
679
 
    }
680
 
 
681
 
    allEdgeList = 0;
682
 
}
683
 
 
684
 
void VoronoiDiagramGenerator::pushGraphEdge(double x1, double y1, double x2, double y2)
685
 
{
686
 
    GraphEdge* newEdge = new GraphEdge;
687
 
    newEdge->next = allEdges;
688
 
    allEdges = newEdge;
689
 
    newEdge->x1 = x1;
690
 
    newEdge->y1 = y1;
691
 
    newEdge->x2 = x2;
692
 
    newEdge->y2 = y2;
693
 
}
694
 
 
695
 
void VoronoiDiagramGenerator::pushEdgeList(Edge *e)
696
 
{
697
 
    EdgeList* newEdge = new EdgeList;
698
 
    newEdge->next = allEdgeList;
699
 
    allEdgeList = newEdge;
700
 
    newEdge->a = e->a;
701
 
    newEdge->b = e->b;
702
 
    newEdge->c = e->c;
703
 
    if (e->ep[0]) {
704
 
        newEdge->ep0nbr = e->ep[0]->sitenbr;
705
 
        newEdge->ep0x = e->ep[0]->coord.x;
706
 
        newEdge->ep0y = e->ep[0]->coord.y;
707
 
    } else {
708
 
        newEdge->ep0nbr = -1;
709
 
    }
710
 
    if (e->ep[1]) {
711
 
        newEdge->ep1nbr = e->ep[1]->sitenbr;
712
 
        newEdge->ep1x = e->ep[1]->coord.x;
713
 
        newEdge->ep1y = e->ep[1]->coord.y;
714
 
    } else {
715
 
        newEdge->ep1nbr = -1;
716
 
    }
717
 
    newEdge->reg0nbr = e->reg[0]->sitenbr;
718
 
    newEdge->reg1nbr = e->reg[1]->sitenbr;
719
 
    newEdge->edgenbr = e->edgenbr;
720
 
}
721
 
 
722
 
char * VoronoiDiagramGenerator::myalloc(unsigned n)
723
 
{
724
 
    char *t=0;    
725
 
    t=(char*)malloc(n);
726
 
    total_alloc += n;
727
 
    return(t);
728
 
}
729
 
 
730
 
 
731
 
/* for those who don't have Cherry's plot */
732
 
/* #include <plot.h> */
733
 
void VoronoiDiagramGenerator::openpl(){}
734
 
void VoronoiDiagramGenerator::line(double x1, double y1, double x2, double y2)
735
 
{    
736
 
    pushGraphEdge(x1,y1,x2,y2);
737
 
 
738
 
}
739
 
void VoronoiDiagramGenerator::circle(double x, double y, double radius){}
740
 
void VoronoiDiagramGenerator::range(double minX, double minY, double maxX, double maxY){}
741
 
 
742
 
 
743
 
 
744
 
void VoronoiDiagramGenerator::out_bisector(struct Edge *e)
745
 
{
746
 
    
747
 
 
748
 
}
749
 
 
750
 
 
751
 
void VoronoiDiagramGenerator::out_ep(struct Edge *e)
752
 
{
753
 
    
754
 
    
755
 
}
756
 
 
757
 
void VoronoiDiagramGenerator::out_vertex(struct Site *v)
758
 
{
759
 
    
760
 
}
761
 
 
762
 
 
763
 
void VoronoiDiagramGenerator::out_site(struct Site *s)
764
 
{
765
 
    if(!triangulate & plot & !debug)
766
 
        circle (s->coord.x, s->coord.y, cradius);
767
 
    
768
 
}
769
 
 
770
 
 
771
 
void VoronoiDiagramGenerator::out_triple(struct Site *s1, struct Site *s2,struct Site * s3)
772
 
{
773
 
    
774
 
}
775
 
 
776
 
 
777
 
 
778
 
void VoronoiDiagramGenerator::plotinit()
779
 
{
780
 
//    double dx,dy,d;
781
 
//    
782
 
//    dy = ymax - ymin;
783
 
//    dx = xmax - xmin;
784
 
//    d = (double)(( dx > dy ? dx : dy) * 1.1);
785
 
//    pxmin = (double)(xmin - (d-dx)/2.0);
786
 
//    pxmax = (double)(xmax + (d-dx)/2.0);
787
 
//    pymin = (double)(ymin - (d-dy)/2.0);
788
 
//    pymax = (double)(ymax + (d-dy)/2.0);
789
 
//    cradius = (double)((pxmax - pxmin)/350.0);
790
 
//    openpl();
791
 
//    range(pxmin, pymin, pxmax, pymax);
792
 
}
793
 
 
794
 
 
795
 
void VoronoiDiagramGenerator::clip_line(struct Edge *e)
796
 
{
797
 
//    struct Site *s1, *s2;
798
 
//    double x1=0,x2=0,y1=0,y2=0;
799
 
 
800
 
    pushEdgeList(e);
801
 
 
802
 
//    x1 = e->reg[0]->coord.x;
803
 
//    x2 = e->reg[1]->coord.x;
804
 
//    y1 = e->reg[0]->coord.y;
805
 
//    y2 = e->reg[1]->coord.y;
806
 
//
807
 
//    //if the distance between the two points this line was created from is less than 
808
 
//    //the square root of 2, then ignore it
809
 
//    if(sqrt(((x2 - x1) * (x2 - x1)) + ((y2 - y1) * (y2 - y1))) < minDistanceBetweenSites)
810
 
//    {
811
 
//        return;
812
 
//    }
813
 
//    pxmin = borderMinX;
814
 
//    pxmax = borderMaxX;
815
 
//    pymin = borderMinY;
816
 
//    pymax = borderMaxY;
817
 
//
818
 
//    if(e -> a == 1.0 && e ->b >= 0.0)
819
 
//    {    
820
 
//        s1 = e -> ep[1];
821
 
//        s2 = e -> ep[0];
822
 
//    }
823
 
//    else 
824
 
//    {
825
 
//        s1 = e -> ep[0];
826
 
//        s2 = e -> ep[1];
827
 
//    };
828
 
//    
829
 
//    if(e -> a == 1.0)
830
 
//    {
831
 
//        y1 = pymin;
832
 
//        if (s1!=(struct Site *)NULL && s1->coord.y > pymin)
833
 
//        {
834
 
//            y1 = s1->coord.y;
835
 
//        }
836
 
//        if(y1>pymax) 
837
 
//        {
838
 
//        //    printf("\nClipped (1) y1 = %f to %f",y1,pymax);
839
 
//            y1 = pymax;
840
 
//            //return;
841
 
//        }
842
 
//        x1 = e -> c - e -> b * y1;
843
 
//        y2 = pymax;
844
 
//        if (s2!=(struct Site *)NULL && s2->coord.y < pymax) 
845
 
//            y2 = s2->coord.y;
846
 
//
847
 
//        if(y2<pymin) 
848
 
//        {
849
 
//            //printf("\nClipped (2) y2 = %f to %f",y2,pymin);
850
 
//            y2 = pymin;
851
 
//            //return;
852
 
//        }
853
 
//        x2 = (e->c) - (e->b) * y2;
854
 
//        if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin))) 
855
 
//        {
856
 
//            //printf("\nClipLine jumping out(3), x1 = %f, pxmin = %f, pxmax = %f",x1,pxmin,pxmax);
857
 
//            return;
858
 
//        }
859
 
//        if(x1> pxmax)
860
 
//        {    x1 = pxmax; y1 = (e -> c - x1)/e -> b;};
861
 
//        if(x1<pxmin)
862
 
//        {    x1 = pxmin; y1 = (e -> c - x1)/e -> b;};
863
 
//        if(x2>pxmax)
864
 
//        {    x2 = pxmax; y2 = (e -> c - x2)/e -> b;};
865
 
//        if(x2<pxmin)
866
 
//        {    x2 = pxmin; y2 = (e -> c - x2)/e -> b;};
867
 
//    }
868
 
//    else
869
 
//    {
870
 
//        x1 = pxmin;
871
 
//        if (s1!=(struct Site *)NULL && s1->coord.x > pxmin) 
872
 
//            x1 = s1->coord.x;
873
 
//        if(x1>pxmax) 
874
 
//        {
875
 
//            //printf("\nClipped (3) x1 = %f to %f",x1,pxmin);
876
 
//            //return;
877
 
//            x1 = pxmax;
878
 
//        }
879
 
//        y1 = e -> c - e -> a * x1;
880
 
//        x2 = pxmax;
881
 
//        if (s2!=(struct Site *)NULL && s2->coord.x < pxmax) 
882
 
//            x2 = s2->coord.x;
883
 
//        if(x2<pxmin) 
884
 
//        {
885
 
//            //printf("\nClipped (4) x2 = %f to %f",x2,pxmin);
886
 
//            //return;
887
 
//            x2 = pxmin;
888
 
//        }
889
 
//        y2 = e -> c - e -> a * x2;
890
 
//        if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin))) 
891
 
//        {
892
 
//            //printf("\nClipLine jumping out(6), y1 = %f, pymin = %f, pymax = %f",y2,pymin,pymax);
893
 
//            return;
894
 
//        }
895
 
//        if(y1> pymax)
896
 
//        {    y1 = pymax; x1 = (e -> c - y1)/e -> a;};
897
 
//        if(y1<pymin)
898
 
//        {    y1 = pymin; x1 = (e -> c - y1)/e -> a;};
899
 
//        if(y2>pymax)
900
 
//        {    y2 = pymax; x2 = (e -> c - y2)/e -> a;};
901
 
//        if(y2<pymin)
902
 
//        {    y2 = pymin; x2 = (e -> c - y2)/e -> a;};
903
 
//    };
904
 
//    
905
 
//    //printf("\nPushing line (%f,%f,%f,%f)",x1,y1,x2,y2);
906
 
//    line(x1,y1,x2,y2);
907
 
}
908
 
 
909
 
 
910
 
/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
911
 
deltax, deltay (can all be estimates).
912
 
Performance suffers if they are wrong; better to make nsites,
913
 
deltax, and deltay too big than too small.  (?) */
914
 
 
915
 
bool VoronoiDiagramGenerator::voronoi(int triangulate)
916
 
{
917
 
    struct Site *newsite, *bot, *top, *temp, *p;
918
 
    struct Site *v;
919
 
    struct Point newintstar;
920
 
    int pm;
921
 
    struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
922
 
    struct Edge *e;
923
 
    
924
 
    PQinitialize();
925
 
    bottomsite = nextone();
926
 
    out_site(bottomsite);
927
 
    bool retval = ELinitialize();
928
 
 
929
 
    if(!retval)
930
 
        return false;
931
 
    
932
 
    newsite = nextone();
933
 
    while(1)
934
 
    {
935
 
 
936
 
        if(!PQempty()) 
937
 
            newintstar = PQ_min();
938
 
        
939
 
        //if the lowest site has a smaller y value than the lowest vector intersection, process the site
940
 
        //otherwise process the vector intersection        
941
 
 
942
 
        if (newsite != (struct Site *)NULL     && (PQempty() || newsite -> coord.y < newintstar.y
943
 
            || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x)))
944
 
        {/* new site is smallest - this is a site event*/
945
 
            out_site(newsite);                        //output the site
946
 
            lbnd = ELleftbnd(&(newsite->coord));                //get the first HalfEdge to the LEFT of the new site
947
 
            rbnd = ELright(lbnd);                        //get the first HalfEdge to the RIGHT of the new site
948
 
            bot = rightreg(lbnd);                        //if this halfedge has no edge, , bot = bottom site (whatever that is)
949
 
            e = bisect(bot, newsite);                    //create a new edge that bisects 
950
 
            bisector = HEcreate(e, le);                    //create a new HalfEdge, setting its ELpm field to 0            
951
 
            ELinsert(lbnd, bisector);                    //insert this new bisector edge between the left and right vectors in a linked list    
952
 
 
953
 
            if ((p = intersect(lbnd, bisector)) != (struct Site *) NULL)     //if the new bisector intersects with the left edge, remove the left edge's vertex, and put in the new one
954
 
            {    
955
 
                PQdelete(lbnd);
956
 
                PQinsert(lbnd, p, dist(p,newsite));
957
 
            };
958
 
            lbnd = bisector;                        
959
 
            bisector = HEcreate(e, re);                    //create a new HalfEdge, setting its ELpm field to 1
960
 
            ELinsert(lbnd, bisector);                    //insert the new HE to the right of the original bisector earlier in the IF stmt
961
 
 
962
 
            if ((p = intersect(bisector, rbnd)) != (struct Site *) NULL)    //if this new bisector intersects with the
963
 
            {    
964
 
                PQinsert(bisector, p, dist(p,newsite));            //push the HE into the ordered linked list of vertices
965
 
            };
966
 
            newsite = nextone();    
967
 
        }
968
 
        else if (!PQempty()) /* intersection is smallest - this is a vector event */            
969
 
        {    
970
 
            lbnd = PQextractmin();                        //pop the HalfEdge with the lowest vector off the ordered list of vectors                
971
 
            llbnd = ELleft(lbnd);                        //get the HalfEdge to the left of the above HE
972
 
            rbnd = ELright(lbnd);                        //get the HalfEdge to the right of the above HE
973
 
            rrbnd = ELright(rbnd);                        //get the HalfEdge to the right of the HE to the right of the lowest HE 
974
 
            bot = leftreg(lbnd);                        //get the Site to the left of the left HE which it bisects
975
 
            top = rightreg(rbnd);                        //get the Site to the right of the right HE which it bisects
976
 
 
977
 
            out_triple(bot, top, rightreg(lbnd));        //output the triple of sites, stating that a circle goes through them
978
 
 
979
 
            v = lbnd->vertex;                        //get the vertex that caused this event
980
 
            makevertex(v);                            //set the vertex number - couldn't do this earlier since we didn't know when it would be processed
981
 
            endpoint(lbnd->ELedge,lbnd->ELpm,v);    //set the endpoint of the left HalfEdge to be this vector
982
 
            endpoint(rbnd->ELedge,rbnd->ELpm,v);    //set the endpoint of the right HalfEdge to be this vector
983
 
            ELdelete(lbnd);                            //mark the lowest HE for deletion - can't delete yet because there might be pointers to it in Hash Map    
984
 
            PQdelete(rbnd);                            //remove all vertex events to do with the  right HE
985
 
            ELdelete(rbnd);                            //mark the right HE for deletion - can't delete yet because there might be pointers to it in Hash Map    
986
 
            pm = le;                                //set the pm variable to zero
987
 
            
988
 
            if (bot->coord.y > top->coord.y)        //if the site to the left of the event is higher than the Site
989
 
            {                                        //to the right of it, then swap them and set the 'pm' variable to 1
990
 
                temp = bot; 
991
 
                bot = top; 
992
 
                top = temp; 
993
 
                pm = re;
994
 
            }
995
 
            e = bisect(bot, top);                    //create an Edge (or line) that is between the two Sites. This creates
996
 
                                                    //the formula of the line, and assigns a line number to it
997
 
            bisector = HEcreate(e, pm);                //create a HE from the Edge 'e', and make it point to that edge with its ELedge field
998
 
            ELinsert(llbnd, bisector);                //insert the new bisector to the right of the left HE
999
 
            endpoint(e, re-pm, v);                    //set one endpoint to the new edge to be the vector point 'v'.
1000
 
                                                    //If the site to the left of this bisector is higher than the right
1001
 
                                                    //Site, then this endpoint is put in position 0; otherwise in pos 1
1002
 
            deref(v);                                //delete the vector 'v'
1003
 
 
1004
 
            //if left HE and the new bisector don't intersect, then delete the left HE, and reinsert it 
1005
 
            if((p = intersect(llbnd, bisector)) != (struct Site *) NULL)
1006
 
            {    
1007
 
                PQdelete(llbnd);
1008
 
                PQinsert(llbnd, p, dist(p,bot));
1009
 
            };
1010
 
 
1011
 
            //if right HE and the new bisector don't intersect, then reinsert it 
1012
 
            if ((p = intersect(bisector, rrbnd)) != (struct Site *) NULL)
1013
 
            {    
1014
 
                PQinsert(bisector, p, dist(p,bot));
1015
 
            };
1016
 
        }
1017
 
        else break;
1018
 
    };
1019
 
 
1020
 
    
1021
 
 
1022
 
 
1023
 
    for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd))
1024
 
    {    
1025
 
        e = lbnd -> ELedge;
1026
 
 
1027
 
        clip_line(e);
1028
 
 
1029
 
    };
1030
 
 
1031
 
    cleanup();
1032
 
 
1033
 
    return true;
1034
 
    
1035
 
}
1036
 
 
1037
 
 
1038
 
int scomp(const void *p1,const void *p2)
1039
 
{
1040
 
    struct Point *s1 = (Point*)p1, *s2=(Point*)p2;
1041
 
    if(s1 -> y < s2 -> y) return(-1);
1042
 
    if(s1 -> y > s2 -> y) return(1);
1043
 
    if(s1 -> x < s2 -> x) return(-1);
1044
 
    if(s1 -> x > s2 -> x) return(1);
1045
 
    return(0);
1046
 
}
1047
 
 
1048
 
/* return a single in-storage site */
1049
 
struct Site * VoronoiDiagramGenerator::nextone()
1050
 
{
1051
 
    struct Site *s;
1052
 
    if(siteidx < nsites)
1053
 
    {    
1054
 
        s = &sites[siteidx];
1055
 
        siteidx += 1;
1056
 
        return(s);
1057
 
    }
1058
 
    else    
1059
 
        return( (struct Site *)NULL);
1060
 
}
1061
 
 
1062
 
bool VoronoiDiagramGenerator::getNextDelaunay(int& ep0, double& ep0x, double& ep0y, 
1063
 
                     int& ep1, double& ep1x, double& ep1y,
1064
 
                     int& reg0, int& reg1)
1065
 
{
1066
 
    if (iterEdgeList == 0)
1067
 
        return false;
1068
 
 
1069
 
    ep0 = iterEdgeList->ep0nbr;
1070
 
    ep0x = iterEdgeList->ep0x;
1071
 
    ep0y = iterEdgeList->ep0y;
1072
 
    ep1 = iterEdgeList->ep1nbr;
1073
 
    ep1x = iterEdgeList->ep1x;
1074
 
    ep1y = iterEdgeList->ep1y;
1075
 
    reg0 = iterEdgeList->reg0nbr;
1076
 
    reg1 = iterEdgeList->reg1nbr;
1077
 
 
1078
 
    iterEdgeList = iterEdgeList->next;
1079
 
    
1080
 
    return true;
1081
 
}
1082
 
 
1083
 
//PyObject* VoronoiDiagramGenerator::_getMesh()
1084
 
//{
1085
 
//    PyObject *vlist, *dlist, *tlist;
1086
 
//    PyObject *temp, *faces, *face;
1087
 
//    int tri0, tri1, reg0, reg1;
1088
 
//    double tri0x, tri0y, tri1x, tri1y;
1089
 
//    int length, numtri, i;
1090
 
//
1091
 
//    length = nedges;
1092
 
//    numtri = nvertices;
1093
 
//
1094
 
//    dlist = PyList_New(length);
1095
 
//    if (!dlist) goto fail;
1096
 
//    vlist = PyList_New(numtri);
1097
 
//    if (!vlist) goto fail;
1098
 
//    tlist = PyList_New(numtri);
1099
 
//    if (!tlist) goto fail;
1100
 
//
1101
 
//    for (i=0; i<numtri; i++) {
1102
 
//        faces = PyList_New(0);
1103
 
//        if (!faces) goto fail;
1104
 
//        PyList_SET_ITEM(tlist, i, faces);
1105
 
//    }
1106
 
//
1107
 
//    resetEdgeListIter();
1108
 
//    i = -1;
1109
 
//    while (getNextDelaunay(tri0, tri0x, tri0y, tri1, tri1x, tri1y, reg0, reg1)) {
1110
 
//        i++;
1111
 
//        face = Py_BuildValue("(ii)", reg0, reg1);
1112
 
//        if (!face) goto fail;
1113
 
//        PyList_SET_ITEM(dlist, i, face);
1114
 
//        if (tri0 > -1) {
1115
 
//            temp = PyList_GET_ITEM(vlist, tri0);
1116
 
//            if (!temp) {
1117
 
//                temp = Py_BuildValue("(dd)", tri0x, tri0y);
1118
 
//                PyList_SET_ITEM(vlist, tri0, temp);
1119
 
//            }
1120
 
//            faces = PyList_GET_ITEM(tlist, tri0);
1121
 
//            if (PyList_Append(faces, face) < 0) goto fail;
1122
 
//        }
1123
 
//        if (tri1 > -1) {
1124
 
//            temp = PyList_GET_ITEM(vlist, tri1);
1125
 
//            if (!temp) {
1126
 
//                temp = Py_BuildValue("(dd)", tri1x, tri1y);
1127
 
//                PyList_SET_ITEM(vlist, tri1, temp);
1128
 
//            }
1129
 
//            faces = PyList_GET_ITEM(tlist, tri1);
1130
 
//            if (PyList_Append(faces, face) < 0) goto fail;
1131
 
//        }
1132
 
//    }
1133
 
//
1134
 
//    temp = PyTuple_Pack(3, vlist, dlist, tlist);
1135
 
//    if (!temp) goto fail;
1136
 
//
1137
 
//    Py_DECREF(vlist);
1138
 
//    Py_DECREF(dlist);
1139
 
//    Py_DECREF(tlist);
1140
 
//
1141
 
//    return temp;
1142
 
//
1143
 
//fail:
1144
 
//    Py_XDECREF(vlist);
1145
 
//    Py_XDECREF(dlist);
1146
 
//    Py_XDECREF(temp);
1147
 
//    Py_XDECREF(faces);
1148
 
//    Py_XDECREF(face);
1149
 
//    return NULL;
1150
 
//}
1151
 
 
1152