Menu
▾
▴
SF.net SVN: matplotlib:[7293] trunk/matplotlib/lib/matplotlib/delaunay
|
From: <js...@us...> - 2009-07-24 02:09:29
|
Revision: 7293
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=7293&view=rev
Author: jswhit
Date: 2009-07-24 02:09:20 +0000 (Fri, 24 Jul 2009)
Log Message:
-----------
merge from scikits.delaunay trunk
Modified Paths:
--------------
trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.cpp
trunk/matplotlib/lib/matplotlib/delaunay/__init__.py
trunk/matplotlib/lib/matplotlib/delaunay/_delaunay.cpp
trunk/matplotlib/lib/matplotlib/delaunay/triangulate.py
Modified: trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.cpp
===================================================================
--- trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.cpp 2009-07-23 22:23:04 UTC (rev 7292)
+++ trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.cpp 2009-07-24 02:09:20 UTC (rev 7293)
@@ -12,9 +12,9 @@
* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
*/
-/*
- * This code was originally written by Stephan Fortune in C code. Shane O'Sullivan,
- * have since modified it, encapsulating it in a C++ class and, fixing memory leaks and
+/*
+ * This code was originally written by Stephan Fortune in C code. Shane O'Sullivan,
+ * have since modified it, encapsulating it in a C++ class and, fixing memory leaks and
* adding accessors to the Voronoi Edges.
* Permission to use, copy, modify, and distribute this software for any
* purpose without fee is hereby granted, provided that this entire notice
@@ -26,7 +26,7 @@
* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
*/
-
+
/*
* Subsequently, Robert Kern modified it to yield Python objects.
* Copyright 2005 Robert Kern <rob...@gm...>
@@ -78,9 +78,9 @@
nsites=numPoints;
plot = 0;
- triangulate = 0;
+ triangulate = 0;
debug = 1;
- sorted = 0;
+ sorted = 0;
freeinit(&sfl, sizeof (Site));
sites = (struct Site *) myalloc(nsites*sizeof( *sites));
@@ -112,9 +112,9 @@
//printf("\n%f %f\n",xValues[i],yValues[i]);
}
-
+
qsort(sites, nsites, sizeof (*sites), scomp);
-
+
siteidx = 0;
geominit();
double temp = 0;
@@ -134,9 +134,9 @@
borderMinY = minY;
borderMaxX = maxX;
borderMaxY = maxY;
+
+ siteidx = 0;
- siteidx = 0;
-
voronoi(triangulate);
return true;
@@ -191,25 +191,25 @@
struct Halfedge * VoronoiDiagramGenerator::ELgethash(int b)
{
struct Halfedge *he;
-
- if(b<0 || b>=ELhashsize)
+
+ if(b<0 || b>=ELhashsize)
return((struct Halfedge *) NULL);
- he = ELhash[b];
- if (he == (struct Halfedge *) NULL || he->ELedge != (struct Edge *) DELETED )
+ he = ELhash[b];
+ if (he == (struct Halfedge *) NULL || he->ELedge != (struct Edge *) DELETED )
return (he);
-
+
/* Hash table points to deleted half edge. Patch as necessary. */
ELhash[b] = (struct Halfedge *) NULL;
- if ((he -> ELrefcnt -= 1) == 0)
+ if ((he -> ELrefcnt -= 1) == 0)
makefree((Freenode*)he, &hfl);
return ((struct Halfedge *) NULL);
-}
+}
struct Halfedge * VoronoiDiagramGenerator::ELleftbnd(struct Point *p)
{
int i, bucket;
struct Halfedge *he;
-
+
/* Use hash table to get close to desired halfedge */
bucket = (int)((p->x - xmin)/deltax * ELhashsize); //use the hash function to find the place in the hash map that this HalfEdge should be
@@ -218,12 +218,12 @@
he = ELgethash(bucket);
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
- {
+ {
for(i=1; 1 ; i += 1)
- {
- if ((he=ELgethash(bucket-i)) != (struct Halfedge *) NULL)
+ {
+ if ((he=ELgethash(bucket-i)) != (struct Halfedge *) NULL)
break;
- if ((he=ELgethash(bucket+i)) != (struct Halfedge *) NULL)
+ if ((he=ELgethash(bucket+i)) != (struct Halfedge *) NULL)
break;
};
totalsearch += i;
@@ -232,22 +232,22 @@
/* Now search linear list of halfedges for the correct one */
if (he==ELleftend || (he != ELrightend && right_of(he,p)))
{
- do
+ do
{
he = he -> ELright;
} 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
he = he -> ELleft; //isn't to the right of
}
else //if the point is to the left of the HalfEdge, then search left for the HE just to the left of the point
- do
+ do
{
he = he -> ELleft;
} while (he!=ELleftend && !right_of(he,p));
-
+
/* Update hash table and reference counts */
if(bucket > 0 && bucket <ELhashsize-1)
- {
- if(ELhash[bucket] != (struct Halfedge *) NULL)
+ {
+ if(ELhash[bucket] != (struct Halfedge *) NULL)
{
ELhash[bucket] -> ELrefcnt -= 1;
}
@@ -281,9 +281,9 @@
struct Site * VoronoiDiagramGenerator::leftreg(struct Halfedge *he)
{
- if(he -> ELedge == (struct Edge *)NULL)
+ if(he -> ELedge == (struct Edge *)NULL)
return(bottomsite);
- return( he -> ELpm == le ?
+ return( he -> ELpm == le ?
he -> ELedge -> reg[le] : he -> ELedge -> reg[re]);
}
@@ -297,7 +297,7 @@
}
void VoronoiDiagramGenerator::geominit()
-{
+{
double sn;
freeinit(&efl, sizeof(Edge));
@@ -313,17 +313,17 @@
struct Edge * VoronoiDiagramGenerator::bisect(struct Site *s1, struct Site *s2)
{
double dx,dy,adx,ady;
- struct Edge *newedge;
+ struct Edge *newedge;
newedge = (struct Edge *) getfree(&efl);
-
+
newedge -> reg[0] = s1; //store the sites that this edge is bisecting
newedge -> reg[1] = s2;
- ref(s1);
+ ref(s1);
ref(s2);
newedge -> ep[0] = (struct Site *) NULL; //to begin with, there are no endpoints on the bisector - it goes to infinity
newedge -> ep[1] = (struct Site *) NULL;
-
+
dx = s2->coord.x - s1->coord.x; //get the difference in x dist between the sites
dy = s2->coord.y - s1->coord.y;
adx = dx>0 ? dx : -dx; //make sure that the difference in positive
@@ -331,18 +331,18 @@
newedge -> c = (double)(s1->coord.x * dx + s1->coord.y * dy + (dx*dx + dy*dy)*0.5);//get the slope of the line
if (adx>ady)
- {
+ {
newedge -> a = 1.0; newedge -> b = dy/dx; newedge -> c /= dx;//set formula of line, with x fixed to 1
}
else
- {
+ {
newedge -> b = 1.0; newedge -> a = dx/dy; newedge -> c /= dy;//set formula of line, with y fixed to 1
};
-
+
newedge -> edgenbr = nedges;
//printf("\nbisect(%d) ((%f,%f) and (%f,%f)",nedges,s1->coord.x,s1->coord.y,s2->coord.x,s2->coord.y);
-
+
nedges += 1;
return(newedge);
}
@@ -355,40 +355,40 @@
double d, xint, yint;
int right_of_site;
struct Site *v;
-
+
e1 = el1 -> ELedge;
e2 = el2 -> ELedge;
- if(e1 == (struct Edge*)NULL || e2 == (struct Edge*)NULL)
+ if(e1 == (struct Edge*)NULL || e2 == (struct Edge*)NULL)
return ((struct Site *) NULL);
//if the two edges bisect the same parent, return null
- if (e1->reg[1] == e2->reg[1])
+ if (e1->reg[1] == e2->reg[1])
return ((struct Site *) NULL);
-
+
d = e1->a * e2->b - e1->b * e2->a;
- if (-1.0e-10<d && d<1.0e-10)
+ if (-1.0e-10<d && d<1.0e-10)
return ((struct Site *) NULL);
-
+
xint = (e1->c*e2->b - e2->c*e1->b)/d;
yint = (e2->c*e1->a - e1->c*e2->a)/d;
-
+
if( (e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
(e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
e1->reg[1]->coord.x < e2->reg[1]->coord.x) )
- {
- el = el1;
+ {
+ el = el1;
e = e1;
}
else
- {
- el = el2;
+ {
+ el = el2;
e = e2;
};
-
+
right_of_site = xint >= e -> reg[1] -> coord.x;
- if ((right_of_site && el -> ELpm == le) || (!right_of_site && el -> ELpm == re))
+ if ((right_of_site && el -> ELpm == le) || (!right_of_site && el -> ELpm == re))
return ((struct Site *) NULL);
-
+
//create a new site at the point of intersection - this is a new vector event waiting to happen
v = (struct Site *) getfree(&sfl);
v -> refcnt = 0;
@@ -404,22 +404,22 @@
struct Site *topsite;
int right_of_site, above, fast;
double dxp, dyp, dxs, t1, t2, t3, yl;
-
+
e = el -> ELedge;
topsite = e -> reg[1];
right_of_site = p -> x > topsite -> coord.x;
if(right_of_site && el -> ELpm == le) return(1);
if(!right_of_site && el -> ELpm == re) return (0);
-
+
if (e->a == 1.0)
{ dyp = p->y - topsite->coord.y;
dxp = p->x - topsite->coord.x;
fast = 0;
if ((!right_of_site & (e->b<0.0)) | (right_of_site & (e->b>=0.0)) )
- { above = dyp>= e->b*dxp;
+ { above = dyp>= e->b*dxp;
fast = above;
}
- else
+ else
{ above = p->x + p->y*e->b > e-> c;
if(e->b<0.0) above = !above;
if (!above) fast = 1;
@@ -446,7 +446,7 @@
{
e -> ep[lr] = s;
ref(s);
- if(e -> ep[re-lr]== (struct Site *) NULL)
+ if(e -> ep[re-lr]== (struct Site *) NULL)
return;
clip_line(e);
@@ -477,7 +477,7 @@
void VoronoiDiagramGenerator::deref(struct Site *v)
{
v -> refcnt -= 1;
- if (v -> refcnt == 0 )
+ if (v -> refcnt == 0 )
makefree((Freenode*)v, &sfl);
}
@@ -490,7 +490,7 @@
void VoronoiDiagramGenerator::PQinsert(struct Halfedge *he,struct Site * v, double offset)
{
struct Halfedge *last, *next;
-
+
he -> vertex = v;
ref(v);
he -> ystar = (double)(v -> coord.y + offset);
@@ -498,23 +498,23 @@
while ((next = last -> PQnext) != (struct Halfedge *) NULL &&
(he -> ystar > next -> ystar ||
(he -> ystar == next -> ystar && v -> coord.x > next->vertex->coord.x)))
- {
+ {
last = next;
};
- he -> PQnext = last -> PQnext;
+ he -> PQnext = last -> PQnext;
last -> PQnext = he;
PQcount += 1;
}
-//remove the HalfEdge from the list of vertices
+//remove the HalfEdge from the list of vertices
void VoronoiDiagramGenerator::PQdelete(struct Halfedge *he)
{
struct Halfedge *last;
-
+
if(he -> vertex != (struct Site *) NULL)
- {
+ {
last = &PQhash[PQbucket(he)];
- while (last -> PQnext != he)
+ while (last -> PQnext != he)
last = last -> PQnext;
last -> PQnext = he -> PQnext;
@@ -527,7 +527,7 @@
int VoronoiDiagramGenerator::PQbucket(struct Halfedge *he)
{
int bucket;
-
+
bucket = (int)((he->ystar - ymin)/deltay * PQhashsize);
if (bucket<0) bucket = 0;
if (bucket>=PQhashsize) bucket = PQhashsize-1 ;
@@ -546,7 +546,7 @@
struct Point VoronoiDiagramGenerator::PQ_min()
{
struct Point answer;
-
+
while(PQhash[PQmin].PQnext == (struct Halfedge *)NULL) {PQmin += 1;};
answer.x = PQhash[PQmin].PQnext -> vertex -> coord.x;
answer.y = PQhash[PQmin].PQnext -> ystar;
@@ -556,7 +556,7 @@
struct Halfedge * VoronoiDiagramGenerator::PQextractmin()
{
struct Halfedge *curr;
-
+
curr = PQhash[PQmin].PQnext;
PQhash[PQmin].PQnext = curr -> PQnext;
PQcount -= 1;
@@ -566,8 +566,8 @@
bool VoronoiDiagramGenerator::PQinitialize()
{
- int i;
-
+ int i;
+
PQcount = 0;
PQmin = 0;
PQhashsize = 4 * sqrt_nsites;
@@ -590,23 +590,23 @@
char * VoronoiDiagramGenerator::getfree(struct Freelist *fl)
{
- int i;
+ int i;
struct Freenode *t;
if(fl->head == (struct Freenode *) NULL)
- {
+ {
t = (struct Freenode *) myalloc(sqrt_nsites * fl->nodesize);
if(t == 0)
return 0;
-
+
currentMemoryBlock->next = new FreeNodeArrayList;
currentMemoryBlock = currentMemoryBlock->next;
currentMemoryBlock->memory = t;
currentMemoryBlock->next = 0;
- for(i=0; i<sqrt_nsites; i+=1)
- makefree((struct Freenode *)((char *)t+i*fl->nodesize), fl);
+ for(i=0; i<sqrt_nsites; i+=1)
+ makefree((struct Freenode *)((char *)t+i*fl->nodesize), fl);
};
t = fl -> head;
fl -> head = (fl -> head) -> nextfree;
@@ -725,7 +725,7 @@
char * VoronoiDiagramGenerator::myalloc(unsigned n)
{
- char *t=0;
+ char *t=0;
t=(char*)malloc(n);
total_alloc += n;
return(t);
@@ -736,7 +736,7 @@
/* #include <plot.h> */
void VoronoiDiagramGenerator::openpl(){}
void VoronoiDiagramGenerator::line(double x1, double y1, double x2, double y2)
-{
+{
pushGraphEdge(x1,y1,x2,y2);
}
@@ -747,20 +747,20 @@
void VoronoiDiagramGenerator::out_bisector(struct Edge *e)
{
+
-
}
void VoronoiDiagramGenerator::out_ep(struct Edge *e)
{
-
-
+
+
}
void VoronoiDiagramGenerator::out_vertex(struct Site *v)
{
-
+
}
@@ -768,13 +768,13 @@
{
if(!triangulate & plot & !debug)
circle (s->coord.x, s->coord.y, cradius);
-
+
}
void VoronoiDiagramGenerator::out_triple(struct Site *s1, struct Site *s2,struct Site * s3)
{
-
+
}
@@ -782,7 +782,7 @@
void VoronoiDiagramGenerator::plotinit()
{
// double dx,dy,d;
-//
+//
// dy = ymax - ymin;
// dx = xmax - xmin;
// d = (double)(( dx > dy ? dx : dy) * 1.1);
@@ -808,7 +808,7 @@
// y1 = e->reg[0]->coord.y;
// y2 = e->reg[1]->coord.y;
//
-// //if the distance between the two points this line was created from is less than
+// //if the distance between the two points this line was created from is less than
// //the square root of 2, then ignore it
// if(sqrt(((x2 - x1) * (x2 - x1)) + ((y2 - y1) * (y2 - y1))) < minDistanceBetweenSites)
// {
@@ -820,16 +820,16 @@
// pymax = borderMaxY;
//
// if(e -> a == 1.0 && e ->b >= 0.0)
-// {
+// {
// s1 = e -> ep[1];
// s2 = e -> ep[0];
// }
-// else
+// else
// {
// s1 = e -> ep[0];
// s2 = e -> ep[1];
// };
-//
+//
// if(e -> a == 1.0)
// {
// y1 = pymin;
@@ -837,7 +837,7 @@
// {
// y1 = s1->coord.y;
// }
-// if(y1>pymax)
+// if(y1>pymax)
// {
// // printf("\nClipped (1) y1 = %f to %f",y1,pymax);
// y1 = pymax;
@@ -845,17 +845,17 @@
// }
// x1 = e -> c - e -> b * y1;
// y2 = pymax;
-// if (s2!=(struct Site *)NULL && s2->coord.y < pymax)
+// if (s2!=(struct Site *)NULL && s2->coord.y < pymax)
// y2 = s2->coord.y;
//
-// if(y2<pymin)
+// if(y2<pymin)
// {
// //printf("\nClipped (2) y2 = %f to %f",y2,pymin);
// y2 = pymin;
// //return;
// }
// x2 = (e->c) - (e->b) * y2;
-// if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin)))
+// if (((x1> pxmax) & (x2>pxmax)) | ((x1<pxmin)&(x2<pxmin)))
// {
// //printf("\nClipLine jumping out(3), x1 = %f, pxmin = %f, pxmax = %f",x1,pxmin,pxmax);
// return;
@@ -872,9 +872,9 @@
// else
// {
// x1 = pxmin;
-// if (s1!=(struct Site *)NULL && s1->coord.x > pxmin)
+// if (s1!=(struct Site *)NULL && s1->coord.x > pxmin)
// x1 = s1->coord.x;
-// if(x1>pxmax)
+// if(x1>pxmax)
// {
// //printf("\nClipped (3) x1 = %f to %f",x1,pxmin);
// //return;
@@ -882,16 +882,16 @@
// }
// y1 = e -> c - e -> a * x1;
// x2 = pxmax;
-// if (s2!=(struct Site *)NULL && s2->coord.x < pxmax)
+// if (s2!=(struct Site *)NULL && s2->coord.x < pxmax)
// x2 = s2->coord.x;
-// if(x2<pxmin)
+// if(x2<pxmin)
// {
// //printf("\nClipped (4) x2 = %f to %f",x2,pxmin);
// //return;
// x2 = pxmin;
// }
// y2 = e -> c - e -> a * x2;
-// if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin)))
+// if (((y1> pymax) & (y2>pymax)) | ((y1<pymin)&(y2<pymin)))
// {
// //printf("\nClipLine jumping out(6), y1 = %f, pymin = %f, pymax = %f",y2,pymin,pymax);
// return;
@@ -905,7 +905,7 @@
// if(y2<pymin)
// { y2 = pymin; x2 = (e -> c - y2)/e -> a;};
// };
-//
+//
// //printf("\nPushing line (%f,%f,%f,%f)",x1,y1,x2,y2);
// line(x1,y1,x2,y2);
}
@@ -924,7 +924,7 @@
int pm;
struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
struct Edge *e;
-
+
PQinitialize();
bottomsite = nextone();
out_site(bottomsite);
@@ -932,16 +932,16 @@
if(!retval)
return false;
-
+
newsite = nextone();
while(1)
{
- if(!PQempty())
+ if(!PQempty())
newintstar = PQ_min();
-
+
//if the lowest site has a smaller y value than the lowest vector intersection, process the site
- //otherwise process the vector intersection
+ //otherwise process the vector intersection
if (newsite != (struct Site *)NULL && (PQempty() || newsite -> coord.y < newintstar.y
|| (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x)))
@@ -950,31 +950,31 @@
lbnd = ELleftbnd(&(newsite->coord)); //get the first HalfEdge to the LEFT of the new site
rbnd = ELright(lbnd); //get the first HalfEdge to the RIGHT of the new site
bot = rightreg(lbnd); //if this halfedge has no edge, , bot = bottom site (whatever that is)
- e = bisect(bot, newsite); //create a new edge that bisects
- bisector = HEcreate(e, le); //create a new HalfEdge, setting its ELpm field to 0
- ELinsert(lbnd, bisector); //insert this new bisector edge between the left and right vectors in a linked list
+ e = bisect(bot, newsite); //create a new edge that bisects
+ bisector = HEcreate(e, le); //create a new HalfEdge, setting its ELpm field to 0
+ ELinsert(lbnd, bisector); //insert this new bisector edge between the left and right vectors in a linked list
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
- {
+ {
PQdelete(lbnd);
PQinsert(lbnd, p, dist(p,newsite));
};
- lbnd = bisector;
+ lbnd = bisector;
bisector = HEcreate(e, re); //create a new HalfEdge, setting its ELpm field to 1
ELinsert(lbnd, bisector); //insert the new HE to the right of the original bisector earlier in the IF stmt
if ((p = intersect(bisector, rbnd)) != (struct Site *) NULL) //if this new bisector intersects with the
- {
+ {
PQinsert(bisector, p, dist(p,newsite)); //push the HE into the ordered linked list of vertices
};
- newsite = nextone();
+ newsite = nextone();
}
- else if (!PQempty()) /* intersection is smallest - this is a vector event */
- {
- lbnd = PQextractmin(); //pop the HalfEdge with the lowest vector off the ordered list of vectors
+ else if (!PQempty()) /* intersection is smallest - this is a vector event */
+ {
+ lbnd = PQextractmin(); //pop the HalfEdge with the lowest vector off the ordered list of vectors
llbnd = ELleft(lbnd); //get the HalfEdge to the left of the above HE
rbnd = ELright(lbnd); //get the HalfEdge to the right of the above HE
- rrbnd = ELright(rbnd); //get the HalfEdge to the right of the HE to the right of the lowest HE
+ rrbnd = ELright(rbnd); //get the HalfEdge to the right of the HE to the right of the lowest HE
bot = leftreg(lbnd); //get the Site to the left of the left HE which it bisects
top = rightreg(rbnd); //get the Site to the right of the right HE which it bisects
@@ -984,16 +984,16 @@
makevertex(v); //set the vertex number - couldn't do this earlier since we didn't know when it would be processed
endpoint(lbnd->ELedge,lbnd->ELpm,v); //set the endpoint of the left HalfEdge to be this vector
endpoint(rbnd->ELedge,rbnd->ELpm,v); //set the endpoint of the right HalfEdge to be this vector
- ELdelete(lbnd); //mark the lowest HE for deletion - can't delete yet because there might be pointers to it in Hash Map
+ ELdelete(lbnd); //mark the lowest HE for deletion - can't delete yet because there might be pointers to it in Hash Map
PQdelete(rbnd); //remove all vertex events to do with the right HE
- ELdelete(rbnd); //mark the right HE for deletion - can't delete yet because there might be pointers to it in Hash Map
+ ELdelete(rbnd); //mark the right HE for deletion - can't delete yet because there might be pointers to it in Hash Map
pm = le; //set the pm variable to zero
-
+
if (bot->coord.y > top->coord.y) //if the site to the left of the event is higher than the Site
{ //to the right of it, then swap them and set the 'pm' variable to 1
- temp = bot;
- bot = top;
- top = temp;
+ temp = bot;
+ bot = top;
+ top = temp;
pm = re;
}
e = bisect(bot, top); //create an Edge (or line) that is between the two Sites. This creates
@@ -1005,27 +1005,27 @@
//Site, then this endpoint is put in position 0; otherwise in pos 1
deref(v); //delete the vector 'v'
- //if left HE and the new bisector don't intersect, then delete the left HE, and reinsert it
+ //if left HE and the new bisector don't intersect, then delete the left HE, and reinsert it
if((p = intersect(llbnd, bisector)) != (struct Site *) NULL)
- {
+ {
PQdelete(llbnd);
PQinsert(llbnd, p, dist(p,bot));
};
- //if right HE and the new bisector don't intersect, then reinsert it
+ //if right HE and the new bisector don't intersect, then reinsert it
if ((p = intersect(bisector, rrbnd)) != (struct Site *) NULL)
- {
+ {
PQinsert(bisector, p, dist(p,bot));
};
}
else break;
};
+
-
for(lbnd=ELright(ELleftend); lbnd != ELrightend; lbnd=ELright(lbnd))
- {
+ {
e = lbnd -> ELedge;
clip_line(e);
@@ -1035,7 +1035,7 @@
cleanup();
return true;
-
+
}
@@ -1054,16 +1054,16 @@
{
struct Site *s;
if(siteidx < nsites)
- {
+ {
s = &sites[siteidx];
siteidx += 1;
return(s);
}
- else
+ else
return( (struct Site *)NULL);
}
-bool VoronoiDiagramGenerator::getNextDelaunay(int& ep0, double& ep0x, double& ep0y,
+bool VoronoiDiagramGenerator::getNextDelaunay(int& ep0, double& ep0x, double& ep0y,
int& ep1, double& ep1x, double& ep1y,
int& reg0, int& reg1)
{
@@ -1080,7 +1080,7 @@
reg1 = iterEdgeList->reg1nbr;
iterEdgeList = iterEdgeList->next;
-
+
return true;
}
Modified: trunk/matplotlib/lib/matplotlib/delaunay/__init__.py
===================================================================
--- trunk/matplotlib/lib/matplotlib/delaunay/__init__.py 2009-07-23 22:23:04 UTC (rev 7292)
+++ trunk/matplotlib/lib/matplotlib/delaunay/__init__.py 2009-07-24 02:09:20 UTC (rev 7293)
@@ -8,4 +8,4 @@
from matplotlib._delaunay import delaunay
from triangulate import *
from interpolate import *
-__version__ = "0.1"
+__version__ = "0.2"
Modified: trunk/matplotlib/lib/matplotlib/delaunay/_delaunay.cpp
===================================================================
--- trunk/matplotlib/lib/matplotlib/delaunay/_delaunay.cpp 2009-07-23 22:23:04 UTC (rev 7292)
+++ trunk/matplotlib/lib/matplotlib/delaunay/_delaunay.cpp 2009-07-24 02:09:20 UTC (rev 7293)
@@ -12,8 +12,8 @@
extern "C" {
-static void reorder_edges(int npoints, int ntriangles,
- double *x, double *y,
+static void reorder_edges(int npoints, int ntriangles,
+ double *x, double *y,
int *edge_db, int *tri_edges, int *tri_nbrs)
{
int neighbors[3], nodes[3];
@@ -69,7 +69,7 @@
// Not trusting me? Okay, let's go through it:
// We have three edges to deal with and three nodes. Without loss
- // of generality, let's label the nodes A, B, and C with (A, B)
+ // of generality, let's label the nodes A, B, and C with (A, B)
// forming the first edge in the order they arrive on input.
// Then there are eight possibilities as to how the other edge-tuples
// may be labeled, but only two variations that are going to affect the
@@ -85,7 +85,7 @@
// The second test we need to perform is for counter-clockwiseness.
// Again, there are only two variations that will affect the outcome:
// either ABC is counter-clockwise, or it isn't. In the former case,
- // we're done setting the node order, we just need to associate the
+ // we're done setting the node order, we just need to associate the
// appropriate neighbor triangles with their opposite nodes, something
// which can be done by inspection. In the latter case, to order the
// nodes counter-clockwise, we only have to switch B and C to get
@@ -113,8 +113,8 @@
static PyObject* getMesh(int npoints, double *x, double *y)
{
- PyObject *vertices = NULL, *edge_db = NULL, *tri_edges = NULL, *tri_nbrs = NULL;
- PyObject *temp = NULL;
+ PyObject *vertices, *edge_db, *tri_edges, *tri_nbrs;
+ PyObject *temp;
int tri0, tri1, reg0, reg1;
double tri0x, tri0y, tri1x, tri1y;
int length, numtri, i, j;
@@ -136,21 +136,21 @@
dim[0] = length;
dim[1] = 2;
edge_db = PyArray_SimpleNew(2, dim, PyArray_INT);
- if (!edge_db) goto exit;
+ if (!edge_db) goto fail;
edge_db_ptr = (int*)PyArray_DATA(edge_db);
-
+
dim[0] = numtri;
vertices = PyArray_SimpleNew(2, dim, PyArray_DOUBLE);
- if (!vertices) goto exit;
+ if (!vertices) goto fail;
vertices_ptr = (double*)PyArray_DATA(vertices);
dim[1] = 3;
tri_edges = PyArray_SimpleNew(2, dim, PyArray_INT);
- if (!tri_edges) goto exit;
+ if (!tri_edges) goto fail;
tri_edges_ptr = (int*)PyArray_DATA(tri_edges);
tri_nbrs = PyArray_SimpleNew(2, dim, PyArray_INT);
- if (!tri_nbrs) goto exit;
+ if (!tri_nbrs) goto fail;
tri_nbrs_ptr = (int*)PyArray_DATA(tri_nbrs);
for (i=0; i<(3*numtri); i++) {
@@ -192,17 +192,25 @@
// tri_edges contains lists of edges; convert to lists of nodes in
// counterclockwise order and reorder tri_nbrs to match. Each node
// corresponds to the edge opposite it in the triangle.
- reorder_edges(npoints, numtri, x, y, edge_db_ptr, tri_edges_ptr,
+ reorder_edges(npoints, numtri, x, y, edge_db_ptr, tri_edges_ptr,
tri_nbrs_ptr);
temp = Py_BuildValue("(OOOO)", vertices, edge_db, tri_edges, tri_nbrs);
+ if (!temp) goto fail;
- exit:
+ Py_DECREF(vertices);
+ Py_DECREF(edge_db);
+ Py_DECREF(tri_edges);
+ Py_DECREF(tri_nbrs);
+
+ return temp;
+
+fail:
Py_XDECREF(vertices);
Py_XDECREF(edge_db);
Py_XDECREF(tri_edges);
Py_XDECREF(tri_nbrs);
- return temp;
+ return NULL;
}
static PyObject *linear_planes(int ntriangles, double *x, double *y, double *z,
@@ -213,7 +221,7 @@
int i;
double *planes_ptr;
double x02, y02, z02, x12, y12, z12, xy0212;
-
+
dims[0] = ntriangles;
dims[1] = 3;
planes = PyArray_SimpleNew(2, dims, PyArray_DOUBLE);
@@ -232,15 +240,15 @@
xy0212 = y02/y12;
INDEX3(planes_ptr,i,0) = (z02 - z12 * xy0212) / (x02 - x12 * xy0212);
INDEX3(planes_ptr,i,1) = (z12 - INDEX3(planes_ptr,i,0)*x12) / y12;
- INDEX3(planes_ptr,i,2) = (z[INDEX3(nodes,i,2)] -
- INDEX3(planes_ptr,i,0)*x[INDEX3(nodes,i,2)] -
+ INDEX3(planes_ptr,i,2) = (z[INDEX3(nodes,i,2)] -
+ INDEX3(planes_ptr,i,0)*x[INDEX3(nodes,i,2)] -
INDEX3(planes_ptr,i,1)*y[INDEX3(nodes,i,2)]);
} else {
xy0212 = x02/x12;
INDEX3(planes_ptr,i,1) = (z02 - z12 * xy0212) / (y02 - y12 * xy0212);
INDEX3(planes_ptr,i,0) = (z12 - INDEX3(planes_ptr,i,1)*y12) / x12;
- INDEX3(planes_ptr,i,2) = (z[INDEX3(nodes,i,2)] -
- INDEX3(planes_ptr,i,0)*x[INDEX3(nodes,i,2)] -
+ INDEX3(planes_ptr,i,2) = (z[INDEX3(nodes,i,2)] -
+ INDEX3(planes_ptr,i,0)*x[INDEX3(nodes,i,2)] -
INDEX3(planes_ptr,i,1)*y[INDEX3(nodes,i,2)]);
}
}
@@ -248,24 +256,24 @@
return (PyObject*)planes;
}
-static double linear_interpolate_single(double targetx, double targety,
+static double linear_interpolate_single(double targetx, double targety,
double *x, double *y, int *nodes, int *neighbors,
PyObject *planes, double defvalue, int start_triangle, int *end_triangle)
{
double *planes_ptr;
planes_ptr = (double*)PyArray_DATA(planes);
if (start_triangle == -1) start_triangle = 0;
- *end_triangle = walking_triangles(start_triangle, targetx, targety,
+ *end_triangle = walking_triangles(start_triangle, targetx, targety,
x, y, nodes, neighbors);
if (*end_triangle == -1) return defvalue;
- return (targetx*INDEX3(planes_ptr,*end_triangle,0) +
+ return (targetx*INDEX3(planes_ptr,*end_triangle,0) +
targety*INDEX3(planes_ptr,*end_triangle,1) +
INDEX3(planes_ptr,*end_triangle,2));
}
-static PyObject *linear_interpolate_grid(double x0, double x1, int xsteps,
+static PyObject *linear_interpolate_grid(double x0, double x1, int xsteps,
double y0, double y1, int ysteps,
- PyObject *planes, double defvalue,
+ PyObject *planes, double defvalue,
int npoints, double *x, double *y, int *nodes, int *neighbors)
{
int ix, iy;
@@ -304,10 +312,10 @@
static PyObject *compute_planes_method(PyObject *self, PyObject *args)
{
PyObject *pyx, *pyy, *pyz, *pynodes;
- PyObject *x = NULL, *y = NULL, *z = NULL, *nodes = NULL;
+ PyObject *x, *y, *z, *nodes;
int npoints, ntriangles;
- PyObject *planes = NULL;
+ PyObject *planes;
if (!PyArg_ParseTuple(args, "OOOO", &pyx, &pyy, &pyz, &pynodes)) {
return NULL;
@@ -315,52 +323,58 @@
x = PyArray_FROMANY(pyx, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
if (!x) {
PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
- goto exit;
+ goto fail;
}
y = PyArray_FROMANY(pyy, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
if (!y) {
PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
- goto exit;
+ goto fail;
}
z = PyArray_FROMANY(pyz, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
if (!z) {
PyErr_SetString(PyExc_ValueError, "z must be a 1-D array of floats");
- goto exit;
+ goto fail;
}
npoints = PyArray_DIM(x, 0);
if ((PyArray_DIM(y, 0) != npoints) || (PyArray_DIM(z, 0) != npoints)) {
PyErr_SetString(PyExc_ValueError, "x,y,z arrays must be of equal length");
- goto exit;
+ goto fail;
}
nodes = PyArray_FROMANY(pynodes, PyArray_INT, 2, 2, NPY_IN_ARRAY);
if (!nodes) {
PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
- goto exit;
+ goto fail;
}
ntriangles = PyArray_DIM(nodes, 0);
if (PyArray_DIM(nodes, 1) != 3) {
PyErr_SetString(PyExc_ValueError, "nodes must have shape (ntriangles, 3)");
- goto exit;
+ goto fail;
}
- planes = linear_planes(ntriangles, (double*)PyArray_DATA(x),
+ planes = linear_planes(ntriangles, (double*)PyArray_DATA(x),
(double*)PyArray_DATA(y), (double*)PyArray_DATA(z), (int*)PyArray_DATA(nodes));
-exit:
+ Py_DECREF(x);
+ Py_DECREF(y);
+ Py_DECREF(z);
+ Py_DECREF(nodes);
+
+ return planes;
+
+fail:
Py_XDECREF(x);
Py_XDECREF(y);
Py_XDECREF(z);
Py_XDECREF(nodes);
-
- return planes;
+ return NULL;
}
static PyObject *linear_interpolate_method(PyObject *self, PyObject *args)
{
double x0, x1, y0, y1, defvalue;
int xsteps, ysteps;
- PyObject *pyplanes, *pyx, *pyy, *pynodes, *pyneighbors, *grid = NULL;
- PyObject *planes = NULL, *x = NULL, *y = NULL, *nodes = NULL, *neighbors = NULL;
+ PyObject *pyplanes, *pyx, *pyy, *pynodes, *pyneighbors, *grid;
+ PyObject *planes, *x, *y, *nodes, *neighbors;
int npoints;
@@ -371,47 +385,54 @@
x = PyArray_FROMANY(pyx, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
if (!x) {
PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
- goto exit;
+ goto fail;
}
y = PyArray_FROMANY(pyy, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
if (!y) {
PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
- goto exit;
+ goto fail;
}
npoints = PyArray_DIM(x, 0);
if (PyArray_DIM(y, 0) != npoints) {
PyErr_SetString(PyExc_ValueError, "x,y arrays must be of equal length");
- goto exit;
+ goto fail;
}
planes = PyArray_FROMANY(pyplanes, PyArray_DOUBLE, 2, 2, NPY_IN_ARRAY);
if (!planes) {
PyErr_SetString(PyExc_ValueError, "planes must be a 2-D array of floats");
- goto exit;
+ goto fail;
}
nodes = PyArray_FROMANY(pynodes, PyArray_INT, 2, 2, NPY_IN_ARRAY);
if (!nodes) {
PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
- goto exit;
+ goto fail;
}
neighbors = PyArray_FROMANY(pyneighbors, PyArray_INT, 2, 2, NPY_IN_ARRAY);
if (!neighbors) {
PyErr_SetString(PyExc_ValueError, "neighbors must be a 2-D array of ints");
- goto exit;
+ goto fail;
}
grid = linear_interpolate_grid(x0, x1, xsteps, y0, y1, ysteps,
(PyObject*)planes, defvalue, npoints,
(double*)PyArray_DATA(x), (double*)PyArray_DATA(y),
(int*)PyArray_DATA(nodes), (int*)PyArray_DATA(neighbors));
+
+ Py_DECREF(x);
+ Py_DECREF(y);
+ Py_DECREF(planes);
+ Py_DECREF(nodes);
+ Py_DECREF(neighbors);
- exit:
+ return grid;
+
+fail:
Py_XDECREF(x);
Py_XDECREF(y);
Py_XDECREF(planes);
Py_XDECREF(nodes);
Py_XDECREF(neighbors);
-
- return grid;
+ return NULL;
}
// Thanks to C++'s memory rules, we can't use the usual "goto fail;" method of
@@ -434,7 +455,7 @@
{
PyObject *pyx, *pyy, *pyz, *pycenters, *pynodes, *pyneighbors, *pyintx, *pyinty;
PyObject *x = NULL, *y = NULL, *z = NULL, *centers = NULL, *nodes = NULL,
- *neighbors = NULL, *intx = NULL, *inty = NULL, *intz = NULL;
+ *neighbors = NULL, *intx = NULL, *inty = NULL, *intz;
double defvalue;
int size, npoints, ntriangles;
@@ -485,7 +506,7 @@
return NULL;
}
ntriangles = PyArray_DIM(neighbors, 0);
- if ((PyArray_DIM(nodes, 0) != ntriangles) ||
+ if ((PyArray_DIM(nodes, 0) != ntriangles) ||
(PyArray_DIM(centers, 0) != ntriangles)) {
PyErr_SetString(PyExc_ValueError, "centers,nodes,neighbors must be of equal length");
CLEANUP
@@ -521,13 +542,13 @@
return NULL;
}
- NaturalNeighbors nn(npoints, ntriangles,
+ NaturalNeighbors nn(npoints, ntriangles,
(double*)PyArray_DATA(x), (double*)PyArray_DATA(y),
- (double*)PyArray_DATA(centers), (int*)PyArray_DATA(nodes),
+ (double*)PyArray_DATA(centers), (int*)PyArray_DATA(nodes),
(int*)PyArray_DATA(neighbors));
size = PyArray_Size(intx);
- nn.interpolate_unstructured((double*)PyArray_DATA(z), size,
- (double*)PyArray_DATA(intx), (double*)PyArray_DATA(inty),
+ nn.interpolate_unstructured((double*)PyArray_DATA(z), size,
+ (double*)PyArray_DATA(intx), (double*)PyArray_DATA(inty),
(double*)PyArray_DATA(intz), defvalue);
Py_XDECREF(x);
@@ -554,13 +575,13 @@
static PyObject *nn_interpolate_method(PyObject *self, PyObject *args)
{
PyObject *pyx, *pyy, *pyz, *pycenters, *pynodes, *pyneighbors, *grid;
- PyObject *x = NULL, *y = NULL, *z = NULL, *centers = NULL, *nodes = NULL, *neighbors = NULL;
+ PyObject *x, *y, *z, *centers, *nodes, *neighbors;
double x0, x1, y0, y1, defvalue;
int xsteps, ysteps;
int npoints, ntriangles;
intp dims[2];
- if (!PyArg_ParseTuple(args, "ddiddidOOOOOO", &x0, &x1, &xsteps,
+ if (!PyArg_ParseTuple(args, "ddiddidOOOOOO", &x0, &x1, &xsteps,
&y0, &y1, &ysteps, &defvalue, &pyx, &pyy, &pyz, &pycenters, &pynodes,
&pyneighbors)) {
return NULL;
@@ -608,7 +629,7 @@
return NULL;
}
ntriangles = PyArray_DIM(neighbors, 0);
- if ((PyArray_DIM(nodes, 0) != ntriangles) ||
+ if ((PyArray_DIM(nodes, 0) != ntriangles) ||
(PyArray_DIM(centers, 0) != ntriangles)) {
PyErr_SetString(PyExc_ValueError, "centers,nodes,neighbors must be of equal length");
CLEANUP
@@ -623,11 +644,11 @@
return NULL;
}
- NaturalNeighbors nn(npoints, ntriangles,
+ NaturalNeighbors nn(npoints, ntriangles,
(double*)PyArray_DATA(x), (double*)PyArray_DATA(y),
- (double*)PyArray_DATA(centers), (int*)PyArray_DATA(nodes),
+ (double*)PyArray_DATA(centers), (int*)PyArray_DATA(nodes),
(int*)PyArray_DATA(neighbors));
- nn.interpolate_grid((double*)PyArray_DATA(z),
+ nn.interpolate_grid((double*)PyArray_DATA(z),
x0, x1, xsteps,
y0, y1, ysteps,
(double*)PyArray_DATA(grid),
@@ -642,8 +663,8 @@
static PyObject *delaunay_method(PyObject *self, PyObject *args)
{
- PyObject *pyx, *pyy, *mesh = NULL;
- PyObject *x = NULL, *y = NULL;
+ PyObject *pyx, *pyy, *mesh;
+ PyObject *x, *y;
int npoints;
if (!PyArg_ParseTuple(args, "OO", &pyx, &pyy)) {
@@ -652,31 +673,37 @@
x = PyArray_FROMANY(pyx, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
if (!x) {
PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
- goto exit;
+ goto fail;
}
y = PyArray_FROMANY(pyy, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
if (!y) {
PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
- goto exit;
+ goto fail;
}
npoints = PyArray_DIM(x, 0);
if (PyArray_DIM(y, 0) != npoints) {
PyErr_SetString(PyExc_ValueError, "x and y must have the same length");
- goto exit;
+ goto fail;
}
mesh = getMesh(npoints, (double*)PyArray_DATA(x), (double*)PyArray_DATA(y));
- exit:
+ if (!mesh) goto fail;
+
+ Py_DECREF(x);
+ Py_DECREF(y);
+
+ return mesh;
+
+fail:
Py_XDECREF(x);
Py_XDECREF(y);
-
- return mesh;
+ return NULL;
}
static PyMethodDef delaunay_methods[] = {
- {"delaunay", (PyCFunction)delaunay_method, METH_VARARGS,
+ {"delaunay", (PyCFunction)delaunay_method, METH_VARARGS,
"Compute the Delaunay triangulation of a cloud of 2-D points.\n\n"
"circumcenters, edges, tri_points, tri_neighbors = delaunay(x, y)\n\n"
"x, y -- shape-(npoints,) arrays of floats giving the X and Y coordinates of the points\n"
@@ -703,7 +730,7 @@
PyMODINIT_FUNC init_delaunay(void)
{
PyObject* m;
- m = Py_InitModule3("_delaunay", delaunay_methods,
+ m = Py_InitModule3("_delaunay", delaunay_methods,
"Tools for computing the Delaunay triangulation and some operations on it.\n"
);
if (m == NULL)
Modified: trunk/matplotlib/lib/matplotlib/delaunay/triangulate.py
===================================================================
--- trunk/matplotlib/lib/matplotlib/delaunay/triangulate.py 2009-07-23 22:23:04 UTC (rev 7292)
+++ trunk/matplotlib/lib/matplotlib/delaunay/triangulate.py 2009-07-24 02:09:20 UTC (rev 7293)
@@ -1,8 +1,10 @@
import warnings
+# 2.3 compatibility
try:
set
except NameError:
- from sets import Set as set
+ import sets
+ set = sets.Set
import numpy as np
@@ -96,7 +98,7 @@
# Find the indices of the unique entries
j_sorted = np.lexsort(keys=(self.x, self.y))
mask_unique = np.hstack([
- True,
+ True,
(np.diff(self.x[j_sorted]) != 0) | (np.diff(self.y[j_sorted]) != 0),
])
return j_sorted[mask_unique]
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
