Menu
▾
▴
SF.net SVN: matplotlib:[5805] trunk/matplotlib
|
From: <js...@us...> - 2008-07-22 01:52:15
|
Revision: 5805
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=5805&view=rev
Author: jswhit
Date: 2008-07-22 01:52:12 +0000 (Tue, 22 Jul 2008)
Log Message:
-----------
added scikits.delaunay as matplotlib.delaunay, added griddata function to mlab.
Modified Paths:
--------------
trunk/matplotlib/CHANGELOG
trunk/matplotlib/examples/tests/backend_driver.py
trunk/matplotlib/setup.py
trunk/matplotlib/setupext.py
Added Paths:
-----------
trunk/matplotlib/examples/pylab_examples/griddata_demo.py
trunk/matplotlib/lib/matplotlib/delaunay/
trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.cpp
trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.h
trunk/matplotlib/lib/matplotlib/delaunay/__init__.py
trunk/matplotlib/lib/matplotlib/delaunay/_delaunay.cpp
trunk/matplotlib/lib/matplotlib/delaunay/delaunay_utils.cpp
trunk/matplotlib/lib/matplotlib/delaunay/delaunay_utils.h
trunk/matplotlib/lib/matplotlib/delaunay/interpolate.py
trunk/matplotlib/lib/matplotlib/delaunay/natneighbors.cpp
trunk/matplotlib/lib/matplotlib/delaunay/natneighbors.h
trunk/matplotlib/lib/matplotlib/delaunay/testfuncs.py
trunk/matplotlib/lib/matplotlib/delaunay/triangulate.py
Modified: trunk/matplotlib/CHANGELOG
===================================================================
--- trunk/matplotlib/CHANGELOG 2008-07-21 22:42:52 UTC (rev 5804)
+++ trunk/matplotlib/CHANGELOG 2008-07-22 01:52:12 UTC (rev 5805)
@@ -1,3 +1,8 @@
+2008-07-21 Added scikits.delaunay as matplotlib.delaunay. Added griddata
+ function in matplotlib.mlab, with example (griddata_demo.py) in
+ pylab_examples. griddata function will use mpl_toolkits._natgrid
+ if installed (haven't yet created the toolkit). - JSW
+
2008-07-21 Re-introduced offset_copy that works in the context of the
new transforms. - MGD
Added: trunk/matplotlib/examples/pylab_examples/griddata_demo.py
===================================================================
--- trunk/matplotlib/examples/pylab_examples/griddata_demo.py (rev 0)
+++ trunk/matplotlib/examples/pylab_examples/griddata_demo.py 2008-07-22 01:52:12 UTC (rev 5805)
@@ -0,0 +1,40 @@
+from numpy.random import uniform, seed
+from matplotlib.mlab import griddata
+import matplotlib.pyplot as plt
+import numpy as np
+# make up data.
+#npts = int(raw_input('enter # of random points to plot:'))
+seed(-1)
+npts = 200
+x = uniform(-2,2,npts)
+y = uniform(-2,2,npts)
+z = x*np.exp(-x**2-y**2)
+# define grid.
+xi = np.linspace(-2.1,2.1,100)
+yi = np.linspace(-2.1,2.1,100)
+# grid the data.
+zi = griddata(x,y,z,xi,yi)
+# contour the gridded data, plotting dots at the nonuniform data points.
+CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
+CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)
+plt.colorbar() # draw colorbar
+# plot data points.
+plt.scatter(x,y,marker='o',c='b',s=5)
+plt.xlim(-2,2)
+plt.ylim(-2,2)
+plt.title('griddata test (%d points)' % npts)
+plt.show()
+
+# test case that scikits.delaunay fails on, but natgrid passes..
+#data = np.array([[-1, -1], [-1, 0], [-1, 1],
+# [ 0, -1], [ 0, 0], [ 0, 1],
+# [ 1, -1 - np.finfo(np.float_).eps], [ 1, 0], [ 1, 1],
+# ])
+#x = data[:,0]
+#y = data[:,1]
+#z = x*np.exp(-x**2-y**2)
+## define grid.
+#xi = np.linspace(-1.1,1.1,100)
+#yi = np.linspace(-1.1,1.1,100)
+## grid the data.
+#zi = griddata(x,y,z,xi,yi)
Modified: trunk/matplotlib/examples/tests/backend_driver.py
===================================================================
--- trunk/matplotlib/examples/tests/backend_driver.py 2008-07-21 22:42:52 UTC (rev 5804)
+++ trunk/matplotlib/examples/tests/backend_driver.py 2008-07-22 01:52:12 UTC (rev 5805)
@@ -41,6 +41,7 @@
'contour_demo.py',
'contour_label_demo.py',
'contourf_demo.py',
+ 'griddata_demo.py',
'csd_demo.py',
'custom_ticker1.py',
'customize_rc.py',
Added: trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.cpp
===================================================================
--- trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.cpp (rev 0)
+++ trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.cpp 2008-07-22 01:52:12 UTC (rev 5805)
@@ -0,0 +1,1152 @@
+/*
+ * The author of this software is Steven Fortune. Copyright (c) 1994 by AT&T
+ * Bell Laboratories.
+ * Permission to use, copy, modify, and distribute this software for any
+ * purpose without fee is hereby granted, provided that this entire notice
+ * is included in all copies of any software which is or includes a copy
+ * or modification of this software and in all copies of the supporting
+ * documentation for such software.
+ * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
+ * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
+ * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
+ * 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
+ * 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
+ * is included in all copies of any software which is or includes a copy
+ * or modification of this software and in all copies of the supporting
+ * documentation for such software.
+ * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
+ * WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
+ * 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...>
+ * See LICENSE.txt in the scipy source directory.
+ */
+
+#include "VoronoiDiagramGenerator.h"
+
+VoronoiDiagramGenerator::VoronoiDiagramGenerator()
+{
+ siteidx = 0;
+ sites = 0;
+
+ allMemoryList = new FreeNodeArrayList;
+ allMemoryList->memory = 0;
+ allMemoryList->next = 0;
+ currentMemoryBlock = allMemoryList;
+ allEdges = 0;
+ allEdgeList = 0;
+ iteratorEdges = 0;
+ iterEdgeList = 0;
+ minDistanceBetweenSites = 0;
+}
+
+VoronoiDiagramGenerator::~VoronoiDiagramGenerator()
+{
+ cleanupEdgeList();
+ cleanup();
+ cleanupEdges();
+
+ if(allMemoryList != 0)
+ delete allMemoryList;
+}
+
+
+
+bool VoronoiDiagramGenerator::generateVoronoi(double *xValues, double *yValues, int numPoints, double minX, double maxX, double minY, double maxY, double minDist)
+{
+ cleanupEdgeList();
+ cleanup();
+ cleanupEdges();
+ int i;
+
+ minDistanceBetweenSites = minDist;
+
+ nsites=numPoints;
+ plot = 0;
+ triangulate = 0;
+ debug = 1;
+ sorted = 0;
+ freeinit(&sfl, sizeof (Site));
+
+ sites = (struct Site *) myalloc(nsites*sizeof( *sites));
+
+ if(sites == 0)
+ return false;
+
+ xmin = xValues[0];
+ ymin = yValues[0];
+ xmax = xValues[0];
+ ymax = yValues[0];
+
+ for(i = 0; i< nsites; i++)
+ {
+ sites[i].coord.x = xValues[i];
+ sites[i].coord.y = yValues[i];
+ sites[i].sitenbr = i;
+ sites[i].refcnt = 0;
+
+ if(xValues[i] < xmin)
+ xmin = xValues[i];
+ else if(xValues[i] > xmax)
+ xmax = xValues[i];
+
+ if(yValues[i] < ymin)
+ ymin = yValues[i];
+ else if(yValues[i] > ymax)
+ ymax = yValues[i];
+
+ //printf("\n%f %f\n",xValues[i],yValues[i]);
+ }
+
+ qsort(sites, nsites, sizeof (*sites), scomp);
+
+ siteidx = 0;
+ geominit();
+ double temp = 0;
+ if(minX > maxX)
+ {
+ temp = minX;
+ minX = maxX;
+ maxX = temp;
+ }
+ if(minY > maxY)
+ {
+ temp = minY;
+ minY = maxY;
+ maxY = temp;
+ }
+ borderMinX = minX;
+ borderMinY = minY;
+ borderMaxX = maxX;
+ borderMaxY = maxY;
+
+ siteidx = 0;
+
+ voronoi(triangulate);
+
+ return true;
+}
+
+bool VoronoiDiagramGenerator::ELinitialize()
+{
+ int i;
+ freeinit(&hfl, sizeof **ELhash);
+ ELhashsize = 2 * sqrt_nsites;
+ ELhash = (struct Halfedge **) myalloc ( sizeof *ELhash * ELhashsize);
+
+ if(ELhash == 0)
+ return false;
+
+ for(i=0; i<ELhashsize; i +=1) ELhash[i] = (struct Halfedge *)NULL;
+ ELleftend = HEcreate( (struct Edge *)NULL, 0);
+ ELrightend = HEcreate( (struct Edge *)NULL, 0);
+ ELleftend -> ELleft = (struct Halfedge *)NULL;
+ ELleftend -> ELright = ELrightend;
+ ELrightend -> ELleft = ELleftend;
+ ELrightend -> ELright = (struct Halfedge *)NULL;
+ ELhash[0] = ELleftend;
+ ELhash[ELhashsize-1] = ELrightend;
+
+ return true;
+}
+
+
+struct Halfedge* VoronoiDiagramGenerator::HEcreate(struct Edge *e,int pm)
+{
+ struct Halfedge *answer;
+ answer = (struct Halfedge *) getfree(&hfl);
+ answer -> ELedge = e;
+ answer -> ELpm = pm;
+ answer -> PQnext = (struct Halfedge *) NULL;
+ answer -> vertex = (struct Site *) NULL;
+ answer -> ELrefcnt = 0;
+ return(answer);
+}
+
+
+void VoronoiDiagramGenerator::ELinsert(struct Halfedge *lb, struct Halfedge *newHe)
+{
+ newHe -> ELleft = lb;
+ newHe -> ELright = lb -> ELright;
+ (lb -> ELright) -> ELleft = newHe;
+ lb -> ELright = newHe;
+}
+
+/* Get entry from hash table, pruning any deleted nodes */
+struct Halfedge * VoronoiDiagramGenerator::ELgethash(int b)
+{
+ struct Halfedge *he;
+
+ if(b<0 || b>=ELhashsize)
+ return((struct Halfedge *) NULL);
+ 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)
+ 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
+
+ if(bucket<0) bucket =0; //make sure that the bucket position in within the range of the hash array
+ if(bucket>=ELhashsize) bucket = ELhashsize - 1;
+
+ 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)
+ break;
+ if ((he=ELgethash(bucket+i)) != (struct Halfedge *) NULL)
+ break;
+ };
+ totalsearch += i;
+ };
+ ntry += 1;
+ /* Now search linear list of halfedges for the correct one */
+ if (he==ELleftend || (he != ELrightend && right_of(he,p)))
+ {
+ 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
+ {
+ 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)
+ {
+ ELhash[bucket] -> ELrefcnt -= 1;
+ }
+ ELhash[bucket] = he;
+ ELhash[bucket] -> ELrefcnt += 1;
+ };
+ return (he);
+}
+
+
+/* This delete routine can't reclaim node, since pointers from hash
+table may be present. */
+void VoronoiDiagramGenerator::ELdelete(struct Halfedge *he)
+{
+ (he -> ELleft) -> ELright = he -> ELright;
+ (he -> ELright) -> ELleft = he -> ELleft;
+ he -> ELedge = (struct Edge *)DELETED;
+}
+
+
+struct Halfedge * VoronoiDiagramGenerator::ELright(struct Halfedge *he)
+{
+ return (he -> ELright);
+}
+
+struct Halfedge * VoronoiDiagramGenerator::ELleft(struct Halfedge *he)
+{
+ return (he -> ELleft);
+}
+
+
+struct Site * VoronoiDiagramGenerator::leftreg(struct Halfedge *he)
+{
+ if(he -> ELedge == (struct Edge *)NULL)
+ return(bottomsite);
+ return( he -> ELpm == le ?
+ he -> ELedge -> reg[le] : he -> ELedge -> reg[re]);
+}
+
+struct Site * VoronoiDiagramGenerator::rightreg(struct Halfedge *he)
+{
+ if(he -> ELedge == (struct Edge *)NULL) //if this halfedge has no edge, return the bottom site (whatever that is)
+ return(bottomsite);
+
+ //if the ELpm field is zero, return the site 0 that this edge bisects, otherwise return site number 1
+ return( he -> ELpm == le ? he -> ELedge -> reg[re] : he -> ELedge -> reg[le]);
+}
+
+void VoronoiDiagramGenerator::geominit()
+{
+ double sn;
+
+ freeinit(&efl, sizeof(Edge));
+ nvertices = 0;
+ nedges = 0;
+ sn = (double)nsites+4;
+ sqrt_nsites = (int)sqrt(sn);
+ deltay = ymax - ymin;
+ deltax = xmax - xmin;
+}
+
+
+struct Edge * VoronoiDiagramGenerator::bisect(struct Site *s1, struct Site *s2)
+{
+ double dx,dy,adx,ady;
+ 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(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
+ ady = dy>0 ? dy : -dy;
+ 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);
+}
+
+//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
+struct Site * VoronoiDiagramGenerator::intersect(struct Halfedge *el1, struct Halfedge *el2, struct Point *p)
+{
+ struct Edge *e1,*e2, *e;
+ struct Halfedge *el;
+ 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)
+ return ((struct Site *) NULL);
+
+ //if the two edges bisect the same parent, return null
+ 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)
+ 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;
+ e = e1;
+ }
+ else
+ {
+ 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))
+ 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;
+ v -> coord.x = xint;
+ v -> coord.y = yint;
+ return(v);
+}
+
+/* returns 1 if p is to right of halfedge e */
+int VoronoiDiagramGenerator::right_of(struct Halfedge *el,struct Point *p)
+{
+ struct Edge *e;
+ 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;
+ fast = above;
+ }
+ else
+ { above = p->x + p->y*e->b > e-> c;
+ if(e->b<0.0) above = !above;
+ if (!above) fast = 1;
+ };
+ if (!fast)
+ { dxs = topsite->coord.x - (e->reg[0])->coord.x;
+ above = e->b * (dxp*dxp - dyp*dyp) <
+ dxs*dyp*(1.0+2.0*dxp/dxs + e->b*e->b);
+ if(e->b<0.0) above = !above;
+ };
+ }
+ else /*e->b==1.0 */
+ { yl = e->c - e->a*p->x;
+ t1 = p->y - yl;
+ t2 = p->x - topsite->coord.x;
+ t3 = yl - topsite->coord.y;
+ above = t1*t1 > t2*t2 + t3*t3;
+ };
+ return (el->ELpm==le ? above : !above);
+}
+
+
+void VoronoiDiagramGenerator::endpoint(struct Edge *e,int lr,struct Site * s)
+{
+ e -> ep[lr] = s;
+ ref(s);
+ if(e -> ep[re-lr]== (struct Site *) NULL)
+ return;
+
+ clip_line(e);
+
+ deref(e->reg[le]);
+ deref(e->reg[re]);
+ makefree((Freenode*)e, &efl);
+}
+
+
+double VoronoiDiagramGenerator::dist(struct Site *s,struct Site *t)
+{
+ double dx,dy;
+ dx = s->coord.x - t->coord.x;
+ dy = s->coord.y - t->coord.y;
+ return (double)(sqrt(dx*dx + dy*dy));
+}
+
+
+void VoronoiDiagramGenerator::makevertex(struct Site *v)
+{
+ v -> sitenbr = nvertices;
+ nvertices += 1;
+ out_vertex(v);
+}
+
+
+void VoronoiDiagramGenerator::deref(struct Site *v)
+{
+ v -> refcnt -= 1;
+ if (v -> refcnt == 0 )
+ makefree((Freenode*)v, &sfl);
+}
+
+void VoronoiDiagramGenerator::ref(struct Site *v)
+{
+ v -> refcnt += 1;
+}
+
+//push the HalfEdge into the ordered linked list of vertices
+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);
+ last = &PQhash[PQbucket(he)];
+ 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;
+ last -> PQnext = he;
+ PQcount += 1;
+}
+
+//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)
+ last = last -> PQnext;
+
+ last -> PQnext = he -> PQnext;
+ PQcount -= 1;
+ deref(he -> vertex);
+ he -> vertex = (struct Site *) NULL;
+ };
+}
+
+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 ;
+ if (bucket < PQmin) PQmin = bucket;
+ return(bucket);
+}
+
+
+
+int VoronoiDiagramGenerator::PQempty()
+{
+ return(PQcount==0);
+}
+
+
+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;
+ return (answer);
+}
+
+struct Halfedge * VoronoiDiagramGenerator::PQextractmin()
+{
+ struct Halfedge *curr;
+
+ curr = PQhash[PQmin].PQnext;
+ PQhash[PQmin].PQnext = curr -> PQnext;
+ PQcount -= 1;
+ return(curr);
+}
+
+
+bool VoronoiDiagramGenerator::PQinitialize()
+{
+ int i;
+
+ PQcount = 0;
+ PQmin = 0;
+ PQhashsize = 4 * sqrt_nsites;
+ PQhash = (struct Halfedge *) myalloc(PQhashsize * sizeof *PQhash);
+
+ if(PQhash == 0)
+ return false;
+
+ for(i=0; i<PQhashsize; i+=1) PQhash[i].PQnext = (struct Halfedge *)NULL;
+
+ return true;
+}
+
+
+void VoronoiDiagramGenerator::freeinit(struct Freelist *fl,int size)
+{
+ fl -> head = (struct Freenode *) NULL;
+ fl -> nodesize = size;
+}
+
+char * VoronoiDiagramGenerator::getfree(struct Freelist *fl)
+{
+ 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);
+ };
+ t = fl -> head;
+ fl -> head = (fl -> head) -> nextfree;
+ return((char *)t);
+}
+
+
+
+void VoronoiDiagramGenerator::makefree(struct Freenode *curr,struct Freelist *fl)
+{
+ curr -> nextfree = fl -> head;
+ fl -> head = curr;
+}
+
+void VoronoiDiagramGenerator::cleanup()
+{
+ if(sites != 0)
+ {
+ free(sites);
+ sites = 0;
+ }
+
+ FreeNodeArrayList* current=0, *prev = 0;
+
+ current = prev = allMemoryList;
+
+ while(current->next != 0)
+ {
+ prev = current;
+ current = current->next;
+ free(prev->memory);
+ delete prev;
+ prev = 0;
+ }
+
+ if(current != 0 && current->memory != 0)
+ {
+ free(current->memory);
+ delete current;
+ }
+
+ allMemoryList = new FreeNodeArrayList;
+ allMemoryList->next = 0;
+ allMemoryList->memory = 0;
+ currentMemoryBlock = allMemoryList;
+}
+
+void VoronoiDiagramGenerator::cleanupEdges()
+{
+ GraphEdge* geCurrent = 0, *gePrev = 0;
+ geCurrent = gePrev = allEdges;
+
+ while(geCurrent != 0 && geCurrent->next != 0)
+ {
+ gePrev = geCurrent;
+ geCurrent = geCurrent->next;
+ delete gePrev;
+ }
+
+ allEdges = 0;
+
+}
+
+void VoronoiDiagramGenerator::cleanupEdgeList()
+{
+ EdgeList* elCurrent = 0, *elPrev = 0;
+ elCurrent = elPrev = allEdgeList;
+
+ while (elCurrent != 0 && elCurrent->next != 0)
+ {
+ elPrev = elCurrent;
+ elCurrent = elCurrent->next;
+ delete elPrev;
+ }
+
+ allEdgeList = 0;
+}
+
+void VoronoiDiagramGenerator::pushGraphEdge(double x1, double y1, double x2, double y2)
+{
+ GraphEdge* newEdge = new GraphEdge;
+ newEdge->next = allEdges;
+ allEdges = newEdge;
+ newEdge->x1 = x1;
+ newEdge->y1 = y1;
+ newEdge->x2 = x2;
+ newEdge->y2 = y2;
+}
+
+void VoronoiDiagramGenerator::pushEdgeList(Edge *e)
+{
+ EdgeList* newEdge = new EdgeList;
+ newEdge->next = allEdgeList;
+ allEdgeList = newEdge;
+ newEdge->a = e->a;
+ newEdge->b = e->b;
+ newEdge->c = e->c;
+ if (e->ep[0]) {
+ newEdge->ep0nbr = e->ep[0]->sitenbr;
+ newEdge->ep0x = e->ep[0]->coord.x;
+ newEdge->ep0y = e->ep[0]->coord.y;
+ } else {
+ newEdge->ep0nbr = -1;
+ }
+ if (e->ep[1]) {
+ newEdge->ep1nbr = e->ep[1]->sitenbr;
+ newEdge->ep1x = e->ep[1]->coord.x;
+ newEdge->ep1y = e->ep[1]->coord.y;
+ } else {
+ newEdge->ep1nbr = -1;
+ }
+ newEdge->reg0nbr = e->reg[0]->sitenbr;
+ newEdge->reg1nbr = e->reg[1]->sitenbr;
+ newEdge->edgenbr = e->edgenbr;
+}
+
+char * VoronoiDiagramGenerator::myalloc(unsigned n)
+{
+ char *t=0;
+ t=(char*)malloc(n);
+ total_alloc += n;
+ return(t);
+}
+
+
+/* for those who don't have Cherry's plot */
+/* #include <plot.h> */
+void VoronoiDiagramGenerator::openpl(){}
+void VoronoiDiagramGenerator::line(double x1, double y1, double x2, double y2)
+{
+ pushGraphEdge(x1,y1,x2,y2);
+
+}
+void VoronoiDiagramGenerator::circle(double x, double y, double radius){}
+void VoronoiDiagramGenerator::range(double minX, double minY, double maxX, double maxY){}
+
+
+
+void VoronoiDiagramGenerator::out_bisector(struct Edge *e)
+{
+
+
+}
+
+
+void VoronoiDiagramGenerator::out_ep(struct Edge *e)
+{
+
+
+}
+
+void VoronoiDiagramGenerator::out_vertex(struct Site *v)
+{
+
+}
+
+
+void VoronoiDiagramGenerator::out_site(struct Site *s)
+{
+ 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)
+{
+
+}
+
+
+
+void VoronoiDiagramGenerator::plotinit()
+{
+// double dx,dy,d;
+//
+// dy = ymax - ymin;
+// dx = xmax - xmin;
+// d = (double)(( dx > dy ? dx : dy) * 1.1);
+// pxmin = (double)(xmin - (d-dx)/2.0);
+// pxmax = (double)(xmax + (d-dx)/2.0);
+// pymin = (double)(ymin - (d-dy)/2.0);
+// pymax = (double)(ymax + (d-dy)/2.0);
+// cradius = (double)((pxmax - pxmin)/350.0);
+// openpl();
+// range(pxmin, pymin, pxmax, pymax);
+}
+
+
+void VoronoiDiagramGenerator::clip_line(struct Edge *e)
+{
+// struct Site *s1, *s2;
+// double x1=0,x2=0,y1=0,y2=0;
+
+ pushEdgeList(e);
+
+// x1 = e->reg[0]->coord.x;
+// x2 = e->reg[1]->coord.x;
+// 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
+// //the square root of 2, then ignore it
+// if(sqrt(((x2 - x1) * (x2 - x1)) + ((y2 - y1) * (y2 - y1))) < minDistanceBetweenSites)
+// {
+// return;
+// }
+// pxmin = borderMinX;
+// pxmax = borderMaxX;
+// pymin = borderMinY;
+// pymax = borderMaxY;
+//
+// if(e -> a == 1.0 && e ->b >= 0.0)
+// {
+// s1 = e -> ep[1];
+// s2 = e -> ep[0];
+// }
+// else
+// {
+// s1 = e -> ep[0];
+// s2 = e -> ep[1];
+// };
+//
+// if(e -> a == 1.0)
+// {
+// y1 = pymin;
+// if (s1!=(struct Site *)NULL && s1->coord.y > pymin)
+// {
+// y1 = s1->coord.y;
+// }
+// if(y1>pymax)
+// {
+// // printf("\nClipped (1) y1 = %f to %f",y1,pymax);
+// y1 = pymax;
+// //return;
+// }
+// x1 = e -> c - e -> b * y1;
+// y2 = pymax;
+// if (s2!=(struct Site *)NULL && s2->coord.y < pymax)
+// y2 = s2->coord.y;
+//
+// 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)))
+// {
+// //printf("\nClipLine jumping out(3), x1 = %f, pxmin = %f, pxmax = %f",x1,pxmin,pxmax);
+// return;
+// }
+// if(x1> pxmax)
+// { x1 = pxmax; y1 = (e -> c - x1)/e -> b;};
+// if(x1<pxmin)
+// { x1 = pxmin; y1 = (e -> c - x1)/e -> b;};
+// if(x2>pxmax)
+// { x2 = pxmax; y2 = (e -> c - x2)/e -> b;};
+// if(x2<pxmin)
+// { x2 = pxmin; y2 = (e -> c - x2)/e -> b;};
+// }
+// else
+// {
+// x1 = pxmin;
+// if (s1!=(struct Site *)NULL && s1->coord.x > pxmin)
+// x1 = s1->coord.x;
+// if(x1>pxmax)
+// {
+// //printf("\nClipped (3) x1 = %f to %f",x1,pxmin);
+// //return;
+// x1 = pxmax;
+// }
+// y1 = e -> c - e -> a * x1;
+// x2 = pxmax;
+// if (s2!=(struct Site *)NULL && s2->coord.x < pxmax)
+// x2 = s2->coord.x;
+// 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)))
+// {
+// //printf("\nClipLine jumping out(6), y1 = %f, pymin = %f, pymax = %f",y2,pymin,pymax);
+// return;
+// }
+// if(y1> pymax)
+// { y1 = pymax; x1 = (e -> c - y1)/e -> a;};
+// if(y1<pymin)
+// { y1 = pymin; x1 = (e -> c - y1)/e -> a;};
+// if(y2>pymax)
+// { y2 = pymax; x2 = (e -> c - y2)/e -> a;};
+// 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);
+}
+
+
+/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
+deltax, deltay (can all be estimates).
+Performance suffers if they are wrong; better to make nsites,
+deltax, and deltay too big than too small. (?) */
+
+bool VoronoiDiagramGenerator::voronoi(int triangulate)
+{
+ struct Site *newsite, *bot, *top, *temp, *p;
+ struct Site *v;
+ struct Point newintstar;
+ int pm;
+ struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
+ struct Edge *e;
+
+ PQinitialize();
+ bottomsite = nextone();
+ out_site(bottomsite);
+ bool retval = ELinitialize();
+
+ if(!retval)
+ return false;
+
+ newsite = nextone();
+ while(1)
+ {
+
+ 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
+
+ if (newsite != (struct Site *)NULL && (PQempty() || newsite -> coord.y < newintstar.y
+ || (newsite->coord.y == newintstar.y && newsite->coord.x < newintstar.x)))
+ {/* new site is smallest - this is a site event*/
+ out_site(newsite); //output the site
+ 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
+
+ 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;
+ 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();
+ }
+ 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
+ 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
+
+ out_triple(bot, top, rightreg(lbnd)); //output the triple of sites, stating that a circle goes through them
+
+ v = lbnd->vertex; //get the vertex that caused this event
+ 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
+ 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
+ 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;
+ pm = re;
+ }
+ e = bisect(bot, top); //create an Edge (or line) that is between the two Sites. This creates
+ //the formula of the line, and assigns a line number to it
+ bisector = HEcreate(e, pm); //create a HE from the Edge 'e', and make it point to that edge with its ELedge field
+ ELinsert(llbnd, bisector); //insert the new bisector to the right of the left HE
+ endpoint(e, re-pm, v); //set one endpoint to the new edge to be the vector point 'v'.
+ //If the site to the left of this bisector is higher than the right
+ //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((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 ((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);
+
+ };
+
+ cleanup();
+
+ return true;
+
+}
+
+
+int scomp(const void *p1,const void *p2)
+{
+ struct Point *s1 = (Point*)p1, *s2=(Point*)p2;
+ if(s1 -> y < s2 -> y) return(-1);
+ if(s1 -> y > s2 -> y) return(1);
+ if(s1 -> x < s2 -> x) return(-1);
+ if(s1 -> x > s2 -> x) return(1);
+ return(0);
+}
+
+/* return a single in-storage site */
+struct Site * VoronoiDiagramGenerator::nextone()
+{
+ struct Site *s;
+ if(siteidx < nsites)
+ {
+ s = &sites[siteidx];
+ siteidx += 1;
+ return(s);
+ }
+ else
+ return( (struct Site *)NULL);
+}
+
+bool VoronoiDiagramGenerator::getNextDelaunay(int& ep0, double& ep0x, double& ep0y,
+ int& ep1, double& ep1x, double& ep1y,
+ int& reg0, int& reg1)
+{
+ if (iterEdgeList == 0)
+ return false;
+
+ ep0 = iterEdgeList->ep0nbr;
+ ep0x = iterEdgeList->ep0x;
+ ep0y = iterEdgeList->ep0y;
+ ep1 = iterEdgeList->ep1nbr;
+ ep1x = iterEdgeList->ep1x;
+ ep1y = iterEdgeList->ep1y;
+ reg0 = iterEdgeList->reg0nbr;
+ reg1 = iterEdgeList->reg1nbr;
+
+ iterEdgeList = iterEdgeList->next;
+
+ return true;
+}
+
+//PyObject* VoronoiDiagramGenerator::_getMesh()
+//{
+// PyObject *vlist, *dlist, *tlist;
+// PyObject *temp, *faces, *face;
+// int tri0, tri1, reg0, reg1;
+// double tri0x, tri0y, tri1x, tri1y;
+// int length, numtri, i;
+//
+// length = nedges;
+// numtri = nvertices;
+//
+// dlist = PyList_New(length);
+// if (!dlist) goto fail;
+// vlist = PyList_New(numtri);
+// if (!vlist) goto fail;
+// tlist = PyList_New(numtri);
+// if (!tlist) goto fail;
+//
+// for (i=0; i<numtri; i++) {
+// faces = PyList_New(0);
+// if (!faces) goto fail;
+// PyList_SET_ITEM(tlist, i, faces);
+// }
+//
+// resetEdgeListIter();
+// i = -1;
+// while (getNextDelaunay(tri0, tri0x, tri0y, tri1, tri1x, tri1y, reg0, reg1)) {
+// i++;
+// face = Py_BuildValue("(ii)", reg0, reg1);
+// if (!face) goto fail;
+// PyList_SET_ITEM(dlist, i, face);
+// if (tri0 > -1) {
+// temp = PyList_GET_ITEM(vlist, tri0);
+// if (!temp) {
+// temp = Py_BuildValue("(dd)", tri0x, tri0y);
+// PyList_SET_ITEM(vlist, tri0, temp);
+// }
+// faces = PyList_GET_ITEM(tlist, tri0);
+// if (PyList_Append(faces, face) < 0) goto fail;
+// }
+// if (tri1 > -1) {
+// temp = PyList_GET_ITEM(vlist, tri1);
+// if (!temp) {
+// temp = Py_BuildValue("(dd)", tri1x, tri1y);
+// PyList_SET_ITEM(vlist, tri1, temp);
+// }
+// faces = PyList_GET_ITEM(tlist, tri1);
+// if (PyList_Append(faces, face) < 0) goto fail;
+// }
+// }
+//
+// temp = PyTuple_Pack(3, vlist, dlist, tlist);
+// if (!temp) goto fail;
+//
+// Py_DECREF(vlist);
+// Py_DECREF(dlist);
+// Py_DECREF(tlist);
+//
+// return temp;
+//
+//fail:
+// Py_XDECREF(vlist);
+// Py_XDECREF(dlist);
+// Py_XDECREF(temp);
+// Py_XDECREF(faces);
+// Py_XDECREF(face);
+// return NULL;
+//}
+
+
Added: trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.h
===================================================================
--- trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.h (rev 0)
+++ trunk/matplotlib/lib/matplotlib/delaunay/VoronoiDiagramGenerator.h 2008-07-22 01:52:12 UTC (rev 5805)
@@ -0,0 +1,283 @@
+/*
+* The author of this software is Steven Fortune. Copyright (c) 1994 by AT&T
+* Bell Laboratories.
+* Permission to use, copy, modify, and distribute this software for any
+* purpose without fee is hereby granted, provided that this entire notice
+* is included in all copies of any software which is or includes a copy
+* or modification of this software and in all copies of the supporting
+* documentation for such software.
+* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
+* WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
+* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
+* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
+*/
+
+/*
+* This code was originally written by Stephan Fortune in C code. I, 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
+* is included in all copies of any software which is or includes a copy
+* or modification of this software and in all copies of the supporting
+* documentation for such software.
+* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
+* WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
+* REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
+* OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
+*/
+
+#ifndef VORONOI_DIAGRAM_GENERATOR
+#define VORONOI_DIAGRAM_GENERATOR
+
+#include "Python.h"
+#include "numpy/arrayobject.h"
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <iostream>
+
+
+#ifndef NULL
+#define NULL 0
+#endif
+#define DELETED -2
+
+#define le 0
+#define re 1
+
+#ifndef MAX
+#define MAX(x, y) (x > y ? x: y)
+#endif
+
+struct Freenode
+{
+ struct Freenode *nextfree;
+};
+
+struct FreeNodeArrayList
+{
+ struct Freenode* memory;
+ struct FreeNodeArrayList* next;
+
+};
+
+struct Freelist
+{
+ struct Freenode *head;
+ int nodesize;
+};
+
+struct Point
+{
+ double x,y;
+};
+
+// structure used both for sites and for vertices
+struct Site
+{
+ struct Point coord;
+ int sitenbr;
+ int refcnt;
+};
+
+
+
+struct Edge
+{
+ double a,b,c;
+ struct Site *ep[2];
+ struct Site *reg[2];
+ int edgenbr;
+};
+
+struct EdgeList
+{
+ double a,b,c;
+ int ep0nbr;
+ double ep0x, ep0y;
+ int ep1nbr;
+ double ep1x, ep1y;
+ int reg0nbr;
+ int reg1nbr;
+ int edgenbr;
+ struct EdgeList *next;
+};
+
+struct GraphEdge
+{
+ double x1,y1,x2,y2;
+ struct GraphEdge* next;
+};
+
+
+
+
+struct Halfedge
+{
+ struct Halfedge *ELleft, *ELright;
+ struct Edge *ELedge;
+ int ELrefcnt;
+ char ELpm;
+ struct Site *vertex;
+ double ystar;
+ struct Halfedge *PQnext;
+};
+
+
+
+
+class VoronoiDiagramGenerator
+{
+public:
+ VoronoiDiagramGenerator();
+ ~VoronoiDiagramGenerator();
+
+ bool generateVoronoi(double *xValues, double *yValues, int numPoints, double minX, double maxX, double minY, double maxY, double minDist=0);
+
+ void resetIterator()
+ {
+ iteratorEdges = allEdges;
+ }
+
+ bool getNext(double& x1, double& y1, double& x2, double& y2)
+ {
+ if(iteratorEdges == 0)
+ return false;
+
+ x1 = iteratorEdges->x1;
+ x2 = iteratorEdges->x2;
+ y1 = iteratorEdges->y1;
+ y2 = iteratorEdges->y2;
+
+ iteratorEdges = iteratorEdges->next;
+
+ return true;
+ }
+
+ void resetEdgeListIter()
+ {
+ iterEdgeList = allEdgeList;
+ }
+
+ bool getNextDelaunay(int& ep0, double& ep0x, double& ep0y,
+ int& ep1, double& ep1x, double& ep1y,
+ int& reg0, int& reg1);
+
+ void getNumbers(int& edges, int& vertices) {
+ edges = nedges;
+ vertices = nvertices;
+ }
+
+private:
+ void cleanup();
+ void cleanupEdgeList();
+ void cleanupEdges();
+ char *getfree(struct Freelist *fl);
+ struct Halfedge *PQfind();
+ int PQempty();
+
+
+ struct Halfedge **ELhash;
+ struct Halfedge *HEcreate(), *ELleft(), *ELright(), *ELleftbnd();
+ struct Halfedge *HEcreate(struct Edge *e,int pm);
+
+
+ struct Point PQ_min();
+ struct Halfedge *PQextractmin();
+ void freeinit(struct Freelist *fl,int size);
+ void makefree(struct Freenode *curr,struct Freelist *fl);
+ void geominit();
+ void plotinit();
+ bool voronoi(int triangulate);
+ void ref(struct Site *v);
+ void deref(struct Site *v);
+ void endpoint(struct Edge *e,int lr,struct Site * s);
+
+ void ELdelete(struct Halfedge *he);
+ struct Halfedge *ELleftbnd(struct Point *p);
+ struct Halfedge *ELright(struct Halfedge *he);
+ void makevertex(struct Site *v);
+ void out_triple(struct Site *s1, struct Site *s2,struct Site * s3);
+
+ void PQinsert(struct Halfedge *he,struct Site * v, double offset);
+ void PQdelete(struct Halfedge *he);
+ bool ELinitialize();
+ void ELinsert(struct Halfedge *lb, struct Halfedge *newHe);
+ struct Halfedge *ELgethash(int b);
+ struct Halfedge *ELleft(struct Halfedge *he);
+ struct Site *leftreg(struct Halfedge *he);
+ void out_site(struct Site *s);
+ bool PQinitialize();
+ int PQbucket(struct Halfedge *he);
+ void clip_line(struct Edge *e);
+ char *myalloc(unsigned n);
+ int right_of(struct Halfedge *el,struct Point *p);
+
+ struct Site *rightreg(struct Halfedge *he);
+ struct Edge *bisect(struct Site *s1, struct Site *s2);
+ double dist(struct Site *s,struct Site *t);
+ struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2, struct Point *p=0);
+
+ void out_bisector(struct Edge *e);
+ void out_ep(struct Edge *e);
+ void out_vertex(struct Site *v);
+ struct Site *nextone();
+
+ void pushGraphEdge(double x1, double y1, double x2, double y2);
+ void pushEdgeList(Edge *e);
+
+ void openpl();
+ void line(double x1, double y1, double x2, double y2);
+ void circle(double x, double y, double radius);
+ void range(double minX, double minY, double maxX, double maxY);
+
+
+ struct Freelist hfl;
+ struct Halfedge *ELleftend, *ELrightend;
+ int ELhashsize;
+
+ int triangulate, sorted, plot, debug;
+ double xmin, xmax, ymin, ymax, deltax, deltay;
+
+ struct Site *sites;
+ int nsites;
+ int siteidx;
+ int sqrt_nsites;
+ int nvertices;
+ struct Freelist sfl;
+ struct Site *bottomsite;
+
+ int nedges;
+ struct Freelist efl;
+ int PQhashsize;
+ struct Halfedge *PQhash;
+ int PQcount;
+ int PQmin;
+
+ int ntry, totalsearch;
+ double pxmin, pxmax, pymin, pymax, cradius;
+ int total_alloc;
+
+ double borderMinX, borderMaxX, borderMinY, borderMaxY;
+
+ FreeNodeArrayList* allMemoryList;
+ FreeNodeArrayList* currentMemoryBlock;
+
+ GraphEdge* allEdges;
+ GraphEdge* iteratorEdges;
+
+ EdgeList* allEdgeList;
+ EdgeList* iterEdgeList;
+
+ double minDistanceBetweenSites;
+
+};
+
+int scomp(const void *p1, const void *p2);
+
+
+#endif
+
+
Added: trunk/matplotlib/lib/matplotlib/delaunay/__init__.py
===================================================================
--- trunk/matplotlib/lib/matplotlib/delaunay/__init__.py (rev 0)
+++ trunk/matplotlib/lib/matplotlib/delaunay/__init__.py 2008-07-22 01:52:12 UTC (rev 5805)
@@ -0,0 +1,10 @@
+"""Delaunay triangulation and interpolation tools.
+
+:Author: Robert Kern <rob...@gm...>
+:Copyright: Copyright 2005 Robert Kern.
+:License: BSD-style license. See LICENSE.txt in the scipy source directory.
+"""
+
+from matplotlib._delaunay import delaunay
+from triangulate import *
+from interpolate import *
Added: trunk/matplotlib/lib/matplotlib/delaunay/_delaunay.cpp
===================================================================
--- trunk/matplotlib/lib/matplotlib/delaunay/_delaunay.cpp (rev 0)
+++ trunk/matplotlib/lib/matplotlib/delaunay/_delaunay.cpp 2008-07-22 01:52:12 UTC (rev 5805)
@@ -0,0 +1,741 @@
+#include "Python.h"
+#include <stdlib.h>
+#include <map>
+#include <iostream>
+
+#include "VoronoiDiagramGenerator.h"
+#include "delaunay_utils.h"
+#include "natneighbors.h"
+#include "numpy/noprefix.h"
+
+using namespace std;
+
+extern "C" {
+
+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];
+ int i, tmp;
+ int case1, case2;
+
+ for (i=0; i<ntriangles; i++) {
+ nodes[0] = INDEX2(edge_db, INDEX3(tri_edges,i,0), 0);
+ nodes[1] = INDEX2(edge_db, INDEX3(tri_edges,i,0), 1);
+ tmp = INDEX2(edge_db, INDEX3(tri_edges,i,1), 0);
+ if (tmp == nodes[0]) {
+ case1 = 1;
+ nodes[2] = INDEX2(edge_db, INDEX3(tri_edges,i,1), 1);
+ } else if (tmp == nodes[1]) {
+ case1 = 0;
+ nodes[2] = INDEX2(edge_db, INDEX3(tri_edges,i,1), 1);
+ } else if (INDEX2(edge_db, INDEX3(tri_edges,i,1), 1) == nodes[0]) {
+ case1 = 1;
+ nodes[2] = tmp;
+ } else {
+ case1 = 0;
+ nodes[2] = tmp;
+ }
+
+ if (ONRIGHT(x[nodes[0]], y[nodes[0]],
+ x[nodes[1]], y[nodes[1]],
+ x[nodes[2]], y[nodes[2]])) {
+ // flip to make counter-clockwise
+ tmp = nodes[2];
+ nodes[2] = nodes[1];
+ nodes[1] = tmp;
+ case2 = 1;
+ } else case2 = 0;
+
+ // I worked it out on paper. You're just gonna have to trust me on this.
+ if (!case1 && !case2) {
+ neighbors[0] = INDEX3(tri_nbrs, i, 1);
+ neighbors[1] = INDEX3(tri_nbrs, i, 2);
+ neighbors[2] = INDEX3(tri_nbrs, i, 0);
+ } else if (case1 && !case2) {
+ neighbors[0] = INDEX3(tri_nbrs, i, 2);
+ neighbors[1] = INDEX3(tri_nbrs, i, 1);
+ neighbors[2] = INDEX3(tri_nbrs, i, 0);
+ } else if (!case1 && case2) {
+ neighbors[0] = INDEX3(tri_nbrs, i, 1);
+ neighbors[1] = INDEX3(tri_nbrs, i, 0);
+ neighbors[2] = INDEX3(tri_nbrs, i, 2);
+ } else {
+ neighbors[0] = INDEX3(tri_nbrs, i, 2);
+ neighbors[1] = INDEX3(tri_nbrs, i, 0);
+ neighbors[2] = INDEX3(tri_nbrs, i, 1);
+ }
+
+ // 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)
+ // 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
+ // output:
+ //
+ // AB AB
+ // BC (CB) AC (CA)
+ // CA (AC) BC (CB)
+ //
+ // The distinction is whether A is in the second edge or B is.
+ // This is the test "case1" above.
+ //
+ // 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
+ // 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
+ // nodes ACB. Then we simply set the neighbor list by inspection again.
+ //
+ // CCW CW
+ // AB
+ // BC 120 102 -+
+ // CA |
+ // +- neighbor order
+ // AB |
+ // AC 210 201 -+
+ // BC
+ // ABC ACB -+- node order
+
+
+ INDEX3(tri_edges,i,0) = nodes[0];
+ INDEX3(tri_edges,i,1) = nodes[1];
+ INDEX3(tri_edges,i,2) = nodes[2];
+ INDEX3(tri_nbrs,i,0) = neighbors[0];
+ INDEX3(tri_nbrs,i,1) = neighbors[1];
+ INDEX3(tri_nbrs,i,2) = neighbors[2];
+ }
+}
+
+static PyObject* getMesh(int npoints, double *x, double *y)
+{
+ 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;
+ intp dim[MAX_DIMS];
+ int *edge_db_ptr, *tri_edges_ptr, *tri_nbrs_ptr;
+ double *vertices_ptr;
+ VoronoiDiagramGenerator vdg;
+
+ vdg.generateVoronoi(x, y, npoints, -100, 100, -100, 100, 0);
+ vdg.getNumbers(length, numtri);
+
+ // Count the actual number of edges
+ i = 0;
+ vdg.resetEdgeListIter();
+ while (vdg.getNextDelaunay(tri0, tri0x, tri0y, tri1, tri1x, tri1y, reg0, reg1))
+ i++;
+ length = i;
+
+ dim[0] = length;
+ dim[1] = 2;
+ edge_db = PyArray_SimpleNew(2, dim, PyArray_INT);
+ 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 fail;
+ vertices_ptr = (double*)PyArray_DATA(vertices);
+
+ dim[1] = 3;
+ tri_edges = PyArray_SimpleNew(2, dim, PyArray_INT);
+ 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 fail;
+ tri_nbrs_ptr = (int*)PyArray_DATA(tri_nbrs);
+
+ for (i=0; i<(3*numtri); i++) {
+ tri_edges_ptr[i] = tri_nbrs_ptr[i] = -1;
+ }
+
+ vdg.resetEdgeListIter();
+ i = -1;
+ while (vdg.getNextDelaunay(tri0, tri0x, tri0y, tri1, tri1x, tri1y, reg0, reg1)) {
+ i++;
+ INDEX2(edge_db_ptr,i,0) = reg0;
+ INDEX2(edge_db_ptr,i,1) = reg1;
+ if (tri0 > -1) {
+ INDEX2(vertices_ptr,tri0,0) = tri0x;
+ INDEX2(vertices_ptr,tri0,1) = tri0y;
+ for (j=0; j<3; j++) {
+ if (INDEX3(tri_edges_ptr,tri0,j) == i) break;
+ if (INDEX3(tri_edges_ptr,tri0,j) == -1) {
+ INDEX3(tri_edges_ptr,tri0,j) = i;
+ INDEX3(tri_nbrs_ptr,tri0,j) = tri1;
+ break;
+ }
+ }
+ }
+ if (tri1 > -1) {
+ INDEX2(vertices_ptr,tri1,0) = tri1x;
+ INDEX2(vertices_ptr,tri1,1) = tri1y;
+ for (j=0; j<3; j++) {
+ if (INDEX3(tri_edges_ptr,tri1,j) == i) break;
+ if (INDEX3(tri_edges_ptr,tri1,j) == -1) {
+ INDEX3(tri_edges_ptr,tri1,j) = i;
+ INDEX3(tri_nbrs_ptr,tri1,j) = tri0;
+ break;
+ }
+ }
+ }
+ }
+
+ // 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,
+ tri_nbrs_ptr);
+
+ temp = Py_BuildValue("(OOOO)", vertices, edge_db, tri_edges, tri_nbrs);
+ if (!temp) goto fail;
+
+ 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 NULL;
+}
+
+static PyObject *linear_planes(int ntriangles, double *x, double *y, double *z,
+ int *nodes)
+{
+ intp dims[2];
+ PyObject *planes;
+ 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);
+ if (!planes) return NULL;
+ planes_ptr = (double *)PyArray_DATA(planes);
+
+ for (i=0; i<ntriangles; i++) {
+ x02 = x[INDEX3(nodes,i,0)] - x[INDEX3(nodes,i,2)];
+ y02 = y[INDEX3(nodes,i,0)] - y[INDEX3(nodes,i,2)];
+ z02 = z[INDEX3(nodes,i,0)] - z[INDEX3(nodes,i,2)];
+ x12 = x[INDEX3(nodes,i,1)] - x[INDEX3(nodes,i,2)];
+ y12 = y[INDEX3(nodes,i,1)] - y[INDEX3(nodes,i,2)];
+ z12 = z[INDEX3(nodes,i,1)] - z[INDEX3(nodes,i,2)];
+
+ if (y12 != 0.0) {
+ 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,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,1)*y[INDEX3(nodes,i,2)]);
+ }
+ }
+
+ return (PyObject*)planes;
+}
+
+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,
+ x, y, nodes, neighbors);
+ if (*end_triangle == -1) return defvalue;
+ 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,
+ double y0, double y1, int ysteps,
+ PyObject *planes, double defvalue,
+ int npoints, double *x, double *y, int *nodes, int *neighbors)
+{
+ int ix, iy;
+ double dx, dy, targetx, targety;
+ int rowtri, coltri, tri;
+ PyObject *z;
+ double *z_ptr;
+ intp dims[2];
+
+ dims[0] = ysteps;
+ dims[1] = xsteps;
+ z = PyArray_SimpleNew(2, dims, PyArray_DOUBLE);
+ if (!z) return NULL;
+ z_ptr = (double*)PyArray_DATA(z);
+
+ dx = (x1 - x0) / (xsteps-1);
+ dy = (y1 - y0) / (ysteps-1);
+
+ rowtri = 0;
+ for (iy=0; iy<ysteps; iy++) {
+ targety = y0 + dy*iy;
+ rowtri = walking_triangles(rowtri, x0, targety, x, y, nodes, neighbors);
+ tri = rowtri;
+ for (ix=0; ix<xsteps; ix++) {
+ targetx = x0 + dx*ix;
+ INDEXN(z_ptr, xsteps, iy, ix) = linear_interpolate_single(
+ targetx, targety,
+ x, y, nodes, neighbors, planes, defvalue, tri, &coltri);
+ if (coltri != -1) tri = coltri;
+ }
+ }
+
+ return z;
+}
+
+static PyObject *compute_planes_method(PyObject *self, PyObject *args)
+{
+ PyObject *pyx, *pyy, *pyz, *pynodes;
+ PyObject *x, *y, *z, *nodes;
+ int npoints, ntriangles;
+
+ PyObject *planes;
+
+ if (!PyArg_ParseTuple(args, "OOOO", &pyx, &pyy, &pyz, &pynodes)) {
+ return NULL;
+ }
+ 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 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 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 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 fail;
+ }
+ nodes = PyArray_FROMANY(pynodes, PyA...
[truncated message content] |