[Yt-svn] yt-commit r432 - in trunk/yt/raven: . delaunay
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Fri May 2 16:44:31 PDT 2008
Author: mturk
Date: Fri May 2 16:44:22 2008
New Revision: 432
URL: http://yt.spacepope.org/changeset/432
Log:
We are now going to include the Delaunay Triangulation scikit from the SciPy
project. This module was written by Robert Kern and is released under a
BSD-like license, which is included in LICENSE.txt. Now contours (and,
hopefully, other things) will be available without the installation of the
SciKits or SciPy repositories.
This officially closes #66. Woohoo!
Added:
trunk/yt/raven/delaunay/
trunk/yt/raven/delaunay/LICENSE.txt
trunk/yt/raven/delaunay/README.txt
trunk/yt/raven/delaunay/VoronoiDiagramGenerator.cpp
trunk/yt/raven/delaunay/VoronoiDiagramGenerator.h
trunk/yt/raven/delaunay/__init__.py
trunk/yt/raven/delaunay/_delaunay.cpp
trunk/yt/raven/delaunay/delaunay_utils.cpp
trunk/yt/raven/delaunay/delaunay_utils.h
trunk/yt/raven/delaunay/interpolate.py
trunk/yt/raven/delaunay/natneighbors.cpp
trunk/yt/raven/delaunay/natneighbors.h
trunk/yt/raven/delaunay/setup.py
trunk/yt/raven/delaunay/testfuncs.py
trunk/yt/raven/delaunay/triangulate.py
Modified:
trunk/yt/raven/PlotTypes.py
trunk/yt/raven/setup.py
Modified: trunk/yt/raven/PlotTypes.py
==============================================================================
--- trunk/yt/raven/PlotTypes.py (original)
+++ trunk/yt/raven/PlotTypes.py Fri May 2 16:44:22 2008
@@ -755,7 +755,7 @@
def contourCallback(field, ncont=5, factor=4, take_log=False, clim=None):
try:
- import scipy.sandbox.delaunay as de
+ import delaunay as de
except ImportError:
mylog.warning("Callback failed; no delaunay module")
return lambda a: None
Added: trunk/yt/raven/delaunay/LICENSE.txt
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/LICENSE.txt Fri May 2 16:44:22 2008
@@ -0,0 +1,29 @@
+Copyright (c) 2001, 2002 Enthought, Inc.
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+ a. Redistributions of source code must retain the above copyright notice,
+ this list of conditions and the following disclaimer.
+ b. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+ c. Neither the name of the Enthought nor the names of its contributors
+ may be used to endorse or promote products derived from this software
+ without specific prior written permission.
+
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
+DAMAGE.
+
Added: trunk/yt/raven/delaunay/README.txt
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/README.txt Fri May 2 16:44:22 2008
@@ -0,0 +1,6 @@
+The code in this directory is subject to the license as specified in
+LICENSE.txt. It was primarily written by Robert Kern, and is formally part of
+the SciKits project, a sub-project of SciPy.
+
+For more information, please see __init__.py, or the SciPy website at
+www.scipy.org .
Added: trunk/yt/raven/delaunay/VoronoiDiagramGenerator.cpp
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/VoronoiDiagramGenerator.cpp Fri May 2 16:44:22 2008
@@ -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 <robert.kern at gmail.com>
+ * 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/yt/raven/delaunay/VoronoiDiagramGenerator.h
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/VoronoiDiagramGenerator.h Fri May 2 16:44:22 2008
@@ -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/yt/raven/delaunay/__init__.py
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/__init__.py Fri May 2 16:44:22 2008
@@ -0,0 +1,10 @@
+"""Delaunay triangulation and interpolation tools.
+
+:Author: Robert Kern <robert.kern at gmail.com>
+:Copyright: Copyright 2005 Robert Kern.
+:License: BSD-style license. See LICENSE.txt
+"""
+
+from _delaunay import delaunay
+from triangulate import *
+from interpolate import *
Added: trunk/yt/raven/delaunay/_delaunay.cpp
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/_delaunay.cpp Fri May 2 16:44:22 2008
@@ -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, PyArray_INT, 2, 2, NPY_IN_ARRAY);
+ if (!nodes) {
+ PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
+ goto fail;
+ }
+ ntriangles = PyArray_DIM(nodes, 0);
+ if (PyArray_DIM(nodes, 1) != 3) {
+ PyErr_SetString(PyExc_ValueError, "nodes must have shape (ntriangles, 3)");
+ goto fail;
+ }
+
+ planes = linear_planes(ntriangles, (double*)PyArray_DATA(x),
+ (double*)PyArray_DATA(y), (double*)PyArray_DATA(z), (int*)PyArray_DATA(nodes));
+
+ 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 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;
+ PyObject *planes, *x, *y, *nodes, *neighbors;
+ int npoints;
+
+
+ if (!PyArg_ParseTuple(args, "ddiddidOOOOO", &x0, &x1, &xsteps, &y0, &y1, &ysteps,
+ &defvalue, &pyplanes, &pyx, &pyy, &pynodes, &pyneighbors)) {
+ 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;
+ }
+ npoints = PyArray_DIM(x, 0);
+ if (PyArray_DIM(y, 0) != npoints) {
+ PyErr_SetString(PyExc_ValueError, "x,y arrays must be of equal length");
+ 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 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 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 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);
+
+ return grid;
+
+fail:
+ Py_XDECREF(x);
+ Py_XDECREF(y);
+ Py_XDECREF(planes);
+ Py_XDECREF(nodes);
+ Py_XDECREF(neighbors);
+ return NULL;
+}
+
+// Thanks to C++'s memory rules, we can't use the usual "goto fail;" method of
+// error handling.
+
+#define CLEANUP \
+ Py_XDECREF(x);\
+ Py_XDECREF(y);\
+ Py_XDECREF(z);\
+ Py_XDECREF(intx);\
+ Py_XDECREF(inty);\
+ Py_XDECREF(centers);\
+ Py_XDECREF(nodes);\
+ Py_XDECREF(neighbors);\
+ Py_XDECREF(intz);
+
+#define PyArray_ND(arr) (((PyArrayObject*)arr)->nd)
+
+static PyObject *nn_interpolate_unstructured_method(PyObject *self, PyObject *args)
+{
+ 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;
+ double defvalue;
+ int size, npoints, ntriangles;
+
+ if (!PyArg_ParseTuple(args, "OOdOOOOOO", &pyintx, &pyinty, &defvalue,
+ &pyx, &pyy, &pyz, &pycenters, &pynodes, &pyneighbors)) {
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ centers = PyArray_FROMANY(pycenters, PyArray_DOUBLE, 2, 2, NPY_IN_ARRAY);
+ if (!centers) {
+ PyErr_SetString(PyExc_ValueError, "centers must be a 2-D array of ints");
+ CLEANUP
+ return NULL;
+ }
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ ntriangles = PyArray_DIM(neighbors, 0);
+ if ((PyArray_DIM(nodes, 0) != ntriangles) ||
+ (PyArray_DIM(centers, 0) != ntriangles)) {
+ PyErr_SetString(PyExc_ValueError, "centers,nodes,neighbors must be of equal length");
+ CLEANUP
+ return NULL;
+ }
+ intx = PyArray_FROM_OTF(pyintx, PyArray_DOUBLE, NPY_IN_ARRAY);
+ if (!intx) {
+ PyErr_SetString(PyExc_ValueError, "intx must be an array of floats");
+ CLEANUP
+ return NULL;
+ }
+ inty = PyArray_FROM_OTF(pyinty, PyArray_DOUBLE, NPY_IN_ARRAY);
+ if (!inty) {
+ PyErr_SetString(PyExc_ValueError, "inty must be an array of floats");
+ CLEANUP
+ return NULL;
+ }
+ if (PyArray_ND(intx) != PyArray_ND(inty)) {
+ PyErr_SetString(PyExc_ValueError, "intx,inty must have same shapes");
+ CLEANUP
+ return NULL;
+ }
+ for (int i=0; i<PyArray_ND(intx); i++) {
+ if (PyArray_DIM(intx, i) != PyArray_DIM(inty, i)) {
+ PyErr_SetString(PyExc_ValueError, "intx,inty must have same shapes");
+ CLEANUP
+ return NULL;
+ }
+ }
+ intz = PyArray_SimpleNew(PyArray_ND(intx), PyArray_DIMS(intx), PyArray_DOUBLE);
+ if (!intz) {
+ CLEANUP
+ return NULL;
+ }
+
+ NaturalNeighbors nn(npoints, ntriangles,
+ (double*)PyArray_DATA(x), (double*)PyArray_DATA(y),
+ (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),
+ (double*)PyArray_DATA(intz), defvalue);
+
+ Py_XDECREF(x);
+ Py_XDECREF(y);
+ Py_XDECREF(z);
+ Py_XDECREF(intx);
+ Py_XDECREF(inty);
+ Py_XDECREF(centers);
+ Py_XDECREF(nodes);
+ Py_XDECREF(neighbors);
+ return intz;
+}
+
+#undef CLEANUP
+
+#define CLEANUP \
+ Py_XDECREF(x);\
+ Py_XDECREF(y);\
+ Py_XDECREF(z);\
+ Py_XDECREF(centers);\
+ Py_XDECREF(nodes);\
+ Py_XDECREF(neighbors);
+
+static PyObject *nn_interpolate_method(PyObject *self, PyObject *args)
+{
+ PyObject *pyx, *pyy, *pyz, *pycenters, *pynodes, *pyneighbors, *grid;
+ 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,
+ &y0, &y1, &ysteps, &defvalue, &pyx, &pyy, &pyz, &pycenters, &pynodes,
+ &pyneighbors)) {
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ npoints = PyArray_DIM(x, 0);
+ if (PyArray_DIM(y, 0) != npoints) {
+ PyErr_SetString(PyExc_ValueError, "x,y arrays must be of equal length");
+ CLEANUP
+ return NULL;
+ }
+ centers = PyArray_FROMANY(pycenters, PyArray_DOUBLE, 2, 2, NPY_IN_ARRAY);
+ if (!centers) {
+ PyErr_SetString(PyExc_ValueError, "centers must be a 2-D array of ints");
+ CLEANUP
+ return NULL;
+ }
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ 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");
+ CLEANUP
+ return NULL;
+ }
+ ntriangles = PyArray_DIM(neighbors, 0);
+ if ((PyArray_DIM(nodes, 0) != ntriangles) ||
+ (PyArray_DIM(centers, 0) != ntriangles)) {
+ PyErr_SetString(PyExc_ValueError, "centers,nodes,neighbors must be of equal length");
+ CLEANUP
+ return NULL;
+ }
+
+ dims[0] = ysteps;
+ dims[1] = xsteps;
+ grid = PyArray_SimpleNew(2, dims, PyArray_DOUBLE);
+ if (!grid) {
+ CLEANUP
+ return NULL;
+ }
+
+ NaturalNeighbors nn(npoints, ntriangles,
+ (double*)PyArray_DATA(x), (double*)PyArray_DATA(y),
+ (double*)PyArray_DATA(centers), (int*)PyArray_DATA(nodes),
+ (int*)PyArray_DATA(neighbors));
+ nn.interpolate_grid((double*)PyArray_DATA(z),
+ x0, x1, xsteps,
+ y0, y1, ysteps,
+ (double*)PyArray_DATA(grid),
+ defvalue, 0);
+
+ CLEANUP
+
+ return grid;
+
+}
+#undef CLEANUP
+
+static PyObject *delaunay_method(PyObject *self, PyObject *args)
+{
+ PyObject *pyx, *pyy, *mesh;
+ PyObject *x, *y;
+ int npoints;
+
+ if (!PyArg_ParseTuple(args, "OO", &pyx, &pyy)) {
+ 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;
+ }
+
+ npoints = PyArray_DIM(x, 0);
+ if (PyArray_DIM(y, 0) != npoints) {
+ PyErr_SetString(PyExc_ValueError, "x and y must have the same length");
+ goto fail;
+ }
+
+ mesh = getMesh(npoints, (double*)PyArray_DATA(x), (double*)PyArray_DATA(y));
+
+ if (!mesh) goto fail;
+
+ Py_DECREF(x);
+ Py_DECREF(y);
+
+ return mesh;
+
+fail:
+ Py_XDECREF(x);
+ Py_XDECREF(y);
+ return NULL;
+}
+
+static PyMethodDef delaunay_methods[] = {
+ {"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"
+ "circumcenters -- shape-(numtri,2) array of floats giving the coordinates of the\n"
+ " circumcenters of each triangle (numtri being the number of triangles)\n"
+ "edges -- shape-(nedges,2) array of integers giving the indices into x and y\n"
+ " of each edge in the triangulation\n"
+ "tri_points -- shape-(numtri,3) array of integers giving the indices into x and y\n"
+ " of each node in each triangle\n"
+ "tri_neighbors -- shape-(numtri,3) array of integers giving the indices into circumcenters\n"
+ " tri_points, and tri_neighbors of the neighbors of each triangle\n"},
+ {"compute_planes", (PyCFunction)compute_planes_method, METH_VARARGS,
+ ""},
+ {"linear_interpolate_grid", (PyCFunction)linear_interpolate_method, METH_VARARGS,
+ ""},
+ {"nn_interpolate_grid", (PyCFunction)nn_interpolate_method, METH_VARARGS,
+ ""},
+ {"nn_interpolate_unstructured", (PyCFunction)nn_interpolate_unstructured_method, METH_VARARGS,
+ ""},
+ {NULL, NULL, 0, NULL}
+};
+
+
+PyMODINIT_FUNC init_delaunay(void)
+{
+ PyObject* m;
+ m = Py_InitModule3("_delaunay", delaunay_methods,
+ "Tools for computing the Delaunay triangulation and some operations on it.\n"
+ );
+ if (m == NULL)
+ return;
+ import_array();
+}
+
+} // extern "C"
Added: trunk/yt/raven/delaunay/delaunay_utils.cpp
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/delaunay_utils.cpp Fri May 2 16:44:22 2008
@@ -0,0 +1,134 @@
+#include "delaunay_utils.h"
+#include <algorithm>
+#include <queue>
+#include <vector>
+#include <iostream>
+
+using namespace std;
+
+int walking_triangles(int start, double targetx, double targety,
+ double *x, double *y, int *nodes, int *neighbors)
+{
+ int i, j, k, t;
+
+ if (start == -1) start = 0;
+ t = start;
+ while (1) {
+ for (i=0; i<3; i++) {
+ j = EDGE0(i);
+ k = EDGE1(i);
+ if (ONRIGHT(x[INDEX3(nodes,t,j)], y[INDEX3(nodes,t,j)],
+ x[INDEX3(nodes,t,k)], y[INDEX3(nodes,t,k)],
+ targetx, targety)) {
+ t = INDEX3(neighbors, t, i);
+ if (t < 0) return t; // outside the convex hull
+ break;
+ }
+ }
+ if (i == 3) break;
+ }
+
+ return t;
+}
+
+void getminmax(double *arr, int n, double& minimum, double& maximum)
+{
+ int i;
+ minimum = arr[0];
+ maximum = arr[0];
+ for (i=1; i<n; i++) {
+ if (arr[i] < minimum) {
+ minimum = arr[i];
+ } else if (arr[i] > maximum) {
+ maximum = arr[i];
+ }
+ }
+}
+
+// Find the circumcenter of three 2-D points by Cramer's Rule to find
+// the intersection of two perpendicular bisectors of the triangle's
+// edges.
+// http://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html
+//
+// Return true if successful; return false if points are collinear
+bool circumcenter(double x0, double y0,
+ double x1, double y1,
+ double x2, double y2,
+ double& centerx, double& centery)
+{
+ double D;
+ double x0m2, y1m2, x1m2, y0m2;
+ double x0p2, y1p2, x1p2, y0p2;
+ x0m2 = x0 - x2;
+ y1m2 = y1 - y2;
+ x1m2 = x1 - x2;
+ y0m2 = y0 - y2;
+ x0p2 = x0 + x2;
+ y1p2 = y1 + y2;
+ x1p2 = x1 + x2;
+ y0p2 = y0 + y2;
+
+ D = x0m2*y1m2 - x1m2*y0m2;
+ if ((D < TOLERANCE_EPS) && (D > -TOLERANCE_EPS)) return false;
+
+ centerx = (((x0m2*x0p2 + y0m2*y0p2)/2*y1m2)
+ - (x1m2*x1p2 + y1m2*y1p2)/2*y0m2) / D;
+ centery = (((x1m2*x1p2 + y1m2*y1p2)/2*x0m2)
+ - (x0m2*x0p2 + y0m2*y0p2)/2*x1m2) / D;
+
+ return true;
+}
+
+double signed_area(double x0, double y0,
+ double x1, double y1,
+ double x2, double y2)
+{
+ return 0.5*(x0 * (y1 - y2)
+ + x1 * (y2 - y0)
+ + x2 * (y0 - y1));
+}
+
+ConvexPolygon::ConvexPolygon(){
+ seeded = false;
+}
+
+ConvexPolygon::~ConvexPolygon() {
+}
+
+void ConvexPolygon::seed(double x0c, double y0c) {
+ this->x0 = x0c;
+ this->y0 = y0c;
+}
+
+void ConvexPolygon::push(double x, double y) {
+ if (!seeded) {
+ seed(x, y);
+ seeded = true;
+ return;
+ }
+ SeededPoint xy(this->x0, this->y0, x, y);
+
+ this->points.push_back(xy);
+}
+
+double ConvexPolygon::area() {
+ double A=0.0;
+ int i;
+
+ sort(this->points.begin(), this->points.end());
+
+ this->points.push_back(SeededPoint(x0, y0, x0, y0));
+ int n = (int)this->points.size();
+
+ for (i=0; i<n; i++) {
+ int j = i-1;
+ if (j<0) j = n-1;
+ int k = i+1;
+ if (k>=n) k = 0;
+
+ A += points[i].x*(points[k].y - points[j].y);
+ }
+
+ A *= 0.5;
+ return A;
+}
Added: trunk/yt/raven/delaunay/delaunay_utils.h
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/delaunay_utils.h Fri May 2 16:44:22 2008
@@ -0,0 +1,77 @@
+
+#ifndef _DELAUNAY_UTILS_H
+#define _DELAUNAY_UTILS_H
+
+#include <vector>
+#include <queue>
+#include <functional>
+
+using namespace std;
+
+#define ONRIGHT(x0, y0, x1, y1, x, y) ((y0-y)*(x1-x) > (x0-x)*(y1-y))
+#define EDGE0(node) ((node + 1) % 3)
+#define EDGE1(node) ((node + 2) % 3)
+#define INDEX2(arr,ix,jx) (arr[2*ix+jx])
+#define INDEX3(arr,ix,jx) (arr[3*ix+jx])
+#define INDEXN(arr,N,ix,jx) (arr[N*ix+jx])
+#define SQ(a) ((a)*(a))
+
+#define TOLERANCE_EPS (4e-13)
+#define PERTURB_EPS (1e-3)
+#define GINORMOUS (1e100)
+
+extern int walking_triangles(int start, double targetx, double targety,
+ double *x, double *y, int *nodes, int *neighbors);
+extern void getminmax(double *arr, int n, double& minimum, double& maximum);
+extern bool circumcenter(double x0, double y0,
+ double x1, double y1,
+ double x2, double y2,
+ double& centerx, double& centery);
+extern double signed_area(double x0, double y0,
+ double x1, double y1,
+ double x2, double y2);
+
+class SeededPoint {
+public:
+ SeededPoint() {};
+ SeededPoint(double x0c, double y0c, double xc, double yc) {
+ this->x0 = x0c;
+ this->y0 = y0c;
+ this->x = xc;
+ this->y = yc;
+ };
+ ~SeededPoint() {};
+
+ double x0, y0;
+ double x, y;
+
+ bool operator<(const SeededPoint& p2) const {
+ double test = (this->y0-p2.y)*(this->x-p2.x) - (this->x0-p2.x)*(this->y-p2.y);
+ if (test == 0) {
+ double length1 = SQ(this->x-this->x0) + SQ(this->y-this->y0);
+ double length2 = SQ(p2.x-this->x0) + SQ(p2.y-this->y0);
+
+ return (length2 > length1);
+ } else return (test < 0);
+ }
+
+};
+
+class ConvexPolygon {
+public:
+ ConvexPolygon();
+ ~ConvexPolygon();
+
+ void seed(double x0c, double y0c);
+ void push(double x, double y);
+
+ double area();
+
+// private: // I don't care much for data-hiding
+ double x0, y0;
+ vector<SeededPoint> points;
+ bool seeded;
+};
+
+
+#endif // _DELAUNAY_UTILS_H
Added: trunk/yt/raven/delaunay/interpolate.py
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/interpolate.py Fri May 2 16:44:22 2008
@@ -0,0 +1,160 @@
+import numpy as np
+
+from _delaunay import compute_planes, linear_interpolate_grid, nn_interpolate_grid
+from _delaunay import nn_interpolate_unstructured
+
+__all__ = ['LinearInterpolator', 'NNInterpolator']
+
+def slice2gridspec(key):
+ """Convert a 2-tuple of slices to start,stop,steps for x and y.
+
+ key -- (slice(ystart,ystop,ystep), slice(xtart, xstop, xstep))
+
+ For now, the only accepted step values are imaginary integers (interpreted
+ in the same way numpy.mgrid, etc. do).
+ """
+ if ((len(key) != 2) or
+ (not isinstance(key[0], slice)) or
+ (not isinstance(key[1], slice))):
+ raise ValueError("only 2-D slices, please")
+
+ x0 = key[1].start
+ x1 = key[1].stop
+ xstep = key[1].step
+ if not isinstance(xstep, complex) or int(xstep.real) != xstep.real:
+ raise ValueError("only the [start:stop:numsteps*1j] form supported")
+ xstep = int(xstep.imag)
+ y0 = key[0].start
+ y1 = key[0].stop
+ ystep = key[0].step
+ if not isinstance(ystep, complex) or int(ystep.real) != ystep.real:
+ raise ValueError("only the [start:stop:numsteps*1j] form supported")
+ ystep = int(ystep.imag)
+
+ return x0, x1, xstep, y0, y1, ystep
+
+class LinearInterpolator(object):
+ """Interpolate a function defined on the nodes of a triangulation by
+ using the planes defined by the three function values at each corner of
+ the triangles.
+
+ LinearInterpolator(triangulation, z, default_value=numpy.nan)
+
+ triangulation -- Triangulation instance
+ z -- the function values at each node of the triangulation
+ default_value -- a float giving the default value should the interpolating
+ point happen to fall outside of the convex hull of the triangulation
+
+ At the moment, the only regular rectangular grids are supported for
+ interpolation.
+
+ vals = interp[ystart:ystop:ysteps*1j, xstart:xstop:xsteps*1j]
+
+ vals would then be a (ysteps, xsteps) array containing the interpolated
+ values. These arguments are interpreted the same way as numpy.mgrid.
+
+ Attributes:
+ planes -- (ntriangles, 3) array of floats specifying the plane for each
+ triangle.
+
+ Linear Interpolation
+ --------------------
+ Given the Delauany triangulation (or indeed *any* complete triangulation) we
+ can interpolate values inside the convex hull by locating the enclosing
+ triangle of the interpolation point and returning the value at that point of
+ the plane defined by the three node values.
+
+ f = planes[tri,0]*x + planes[tri,1]*y + planes[tri,2]
+
+ The interpolated function is C0 continuous across the convex hull of the
+ input points. It is C1 continuous across the convex hull except for the
+ nodes and the edges of the triangulation.
+ """
+ def __init__(self, triangulation, z, default_value=np.nan):
+ self.triangulation = triangulation
+ self.z = np.asarray(z, dtype=np.float64)
+ self.default_value = default_value
+
+ self.planes = compute_planes(triangulation.x, triangulation.y, self.z,
+ triangulation.triangle_nodes)
+
+ def __getitem__(self, key):
+ x0, x1, xstep, y0, y1, ystep = slice2gridspec(key)
+ grid = linear_interpolate_grid(x0, x1, xstep, y0, y1, ystep, self.default_value,
+ self.planes, self.triangulation.x, self.triangulation.y,
+ self.triangulation.triangle_nodes, self.triangulation.triangle_neighbors)
+ return grid
+
+class NNInterpolator(object):
+ """Interpolate a function defined on the nodes of a triangulation by
+ the natural neighbors method.
+
+ NNInterpolator(triangulation, z, default_value=numpy.nan)
+
+ triangulation -- Triangulation instance
+ z -- the function values at each node of the triangulation
+ default_value -- a float giving the default value should the interpolating
+ point happen to fall outside of the convex hull of the triangulation
+
+ At the moment, the only regular rectangular grids are supported for
+ interpolation.
+
+ vals = interp[ystart:ystop:ysteps*1j, xstart:xstop:xsteps*1j]
+
+ vals would then be a (ysteps, xsteps) array containing the interpolated
+ values. These arguments are interpreted the same way as numpy.mgrid.
+
+ Natural Neighbors Interpolation
+ -------------------------------
+ One feature of the Delaunay triangulation is that for each triangle, its
+ circumcircle contains no other point (although in degenerate cases, like
+ squares, other points may be *on* the circumcircle). One can also construct
+ what is called the Voronoi diagram from a Delaunay triangulation by
+ connecting the circumcenters of the triangles to those of their neighbors to
+ form a tesselation of irregular polygons covering the plane and containing
+ only one node from the triangulation. Each point in one node's Voronoi
+ polygon is closer to that node than any other node.
+
+ To compute the Natural Neighbors interpolant, we consider adding the
+ interpolation point to the triangulation. We define the natural neighbors of
+ this point as the set of nodes participating in Delaunay triangles whose
+ circumcircles contain the point. To restore the Delaunay-ness of the
+ triangulation, one would only have to alter those triangles and Voronoi
+ polygons. The new Voronooi diagram would have a polygon around the inserted
+ point. This polygon would "steal" area from the original Voronoi polygons.
+ For each node i in the natural neighbors set, we compute the area stolen
+ from its original Voronoi polygon, stolen[i]. We define the natural
+ neighbors coordinates
+
+ phi[i] = stolen[i] / sum(stolen,axis=0)
+
+ We then use these phi[i] to weight the corresponding function values from
+ the input data z to compute the interpolated value.
+
+ The interpolated surface is C1-continuous except at the nodes themselves
+ across the convex hull of the input points. One can find the set of points
+ that a given node will affect by computing the union of the areas covered by
+ the circumcircles of each Delaunay triangle that node participates in.
+ """
+
+ def __init__(self, triangulation, z, default_value=np.nan):
+ self.triangulation = triangulation
+ self.z = np.asarray(z, dtype=np.float64)
+ self.default_value = default_value
+
+ def __getitem__(self, key):
+ x0, x1, xstep, y0, y1, ystep = slice2gridspec(key)
+ grid = nn_interpolate_grid(x0, x1, xstep, y0, y1, ystep, self.default_value,
+ self.triangulation.x, self.triangulation.y, self.z,
+ self.triangulation.circumcenters,
+ self.triangulation.triangle_nodes,
+ self.triangulation.triangle_neighbors)
+ return grid
+
+ def __call__(self, intx, inty):
+ intz = nn_interpolate_unstructured(intx, inty, self.default_value,
+ self.triangulation.x, self.triangulation.y, self.z,
+ self.triangulation.circumcenters,
+ self.triangulation.triangle_nodes,
+ self.triangulation.triangle_neighbors)
+ return intz
Added: trunk/yt/raven/delaunay/natneighbors.cpp
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/natneighbors.cpp Fri May 2 16:44:22 2008
@@ -0,0 +1,299 @@
+
+#include "delaunay_utils.h"
+#include "natneighbors.h"
+
+#include <stack>
+#include <list>
+#include <vector>
+#include <set>
+#include <math.h>
+#include <iostream>
+#include <iterator>
+
+using namespace std;
+
+NaturalNeighbors::NaturalNeighbors(int npoints, int ntriangles, double *x, double *y,
+ double *centers, int *nodes, int *neighbors)
+{
+ this->npoints = npoints;
+ this->ntriangles = ntriangles;
+ this->x = x;
+ this->y = y;
+ this->centers = centers;
+ this->nodes = nodes;
+ this->neighbors = neighbors;
+
+ this->radii2 = new double[ntriangles];
+ for (int i=0; i<ntriangles; i++) {
+ double x2 = x[INDEX3(nodes,i,0)] - INDEX2(centers,i,0);
+ x2 = x2*x2;
+ double y2 = y[INDEX3(nodes,i,0)] - INDEX2(centers,i,1);
+ y2 = y2*y2;
+ this->radii2[i] = x2 + y2;
+ }
+
+}
+
+NaturalNeighbors::~NaturalNeighbors()
+{
+ delete[] this->radii2;
+}
+
+int NaturalNeighbors::find_containing_triangle(double targetx, double targety, int start_triangle)
+{
+ int final_triangle;
+ final_triangle = walking_triangles(start_triangle, targetx, targety,
+ x, y, nodes, neighbors);
+ return final_triangle;
+}
+
+double NaturalNeighbors::interpolate_one(double *z, double targetx, double targety,
+ double defvalue, int &start_triangle)
+{
+ int t = find_containing_triangle(targetx, targety, start_triangle);
+ if (t == -1) return defvalue;
+
+ start_triangle = t;
+ vector<int> circumtri;
+
+ circumtri.push_back(t);
+ stack<int> stackA;
+ stack<int> stackB;
+ int tnew, i;
+
+ for (i=0; i<3; i++) {
+ tnew = INDEX3(this->neighbors, t, i);
+ if (tnew != -1) {
+ stackA.push(tnew);
+ stackB.push(t);
+ }
+ }
+ while (!stackA.empty()) {
+ tnew = stackA.top();
+ stackA.pop();
+ t = stackB.top();
+ stackB.pop();
+ double d2 = (SQ(targetx - INDEX2(this->centers,tnew,0))
+ + SQ(targety - INDEX2(this->centers,tnew,1)));
+ if ((this->radii2[tnew]-d2) > TOLERANCE_EPS) {
+ // tnew is a circumtriangle of the target
+ circumtri.push_back(tnew);
+ for (i=0; i<3; i++) {
+ int ti = INDEX3(this->neighbors, tnew, i);
+ if ((ti != -1) && (ti != t)) {
+ stackA.push(ti);
+ stackB.push(tnew);
+ }
+ }
+ }
+ }
+
+ vector<int>::iterator it;
+ double f = 0.0;
+ double A = 0.0;
+ double tA=0.0, yA=0.0, cA=0.0; // Kahan summation temps for A
+ double tf=0.0, yf=0.0, cf=0.0; // Kahan summation temps for f
+
+ vector<int> edge;
+ bool onedge = false;
+ bool onhull = false;
+
+ for (it = circumtri.begin(); it != circumtri.end(); it++) {
+ int t = *it;
+ double vx = INDEX2(this->centers, t, 0);
+ double vy = INDEX2(this->centers, t, 1);
+ vector<double> c(6);
+ for (int i=0; i<3; i++) {
+ int j = EDGE0(i);
+ int k = EDGE1(i);
+
+ if (!circumcenter(
+ this->x[INDEX3(this->nodes, t, j)],
+ this->y[INDEX3(this->nodes, t, j)],
+ this->x[INDEX3(this->nodes, t, k)],
+ this->y[INDEX3(this->nodes, t, k)],
+ targetx, targety,
+ INDEX2(c, i, 0), INDEX2(c, i, 1))) {
+
+ // bail out with the appropriate values if we're actually on a
+ // node
+ if ((fabs(targetx - this->x[INDEX3(this->nodes, t, j)]) < TOLERANCE_EPS)
+ && (fabs(targety - this->y[INDEX3(this->nodes, t, j)]) < TOLERANCE_EPS)) {
+ return z[INDEX3(this->nodes, t, j)];
+ } else if ((fabs(targetx - this->x[INDEX3(this->nodes, t, k)]) < TOLERANCE_EPS)
+ && (fabs(targety - this->y[INDEX3(this->nodes, t, k)]) < TOLERANCE_EPS)) {
+ return z[INDEX3(this->nodes, t, k)];
+ } else if (!onedge) {
+ onedge = true;
+ edge.push_back(INDEX3(this->nodes, t, j));
+ edge.push_back(INDEX3(this->nodes, t, k));
+ onhull = (INDEX3(neighbors, t, i) == -1);
+ }
+ }
+ }
+ for (int i=0; i<3; i++) {
+ int j = EDGE0(i);
+ int k = EDGE1(i);
+ int q = INDEX3(this->nodes, t, i);
+ double ati = 0.0;
+
+ if (!onedge || ((edge[0] != q) && edge[1] != q)) {
+ ati = signed_area(vx, vy,
+ INDEX2(c, j, 0), INDEX2(c, j, 1),
+ INDEX2(c, k, 0), INDEX2(c, k, 1));
+
+
+ yA = ati - cA;
+ tA = A + yA;
+ cA = (tA - A) - yA;
+ A = tA;
+
+ yf = ati*z[q] - cf;
+ tf = f + yf;
+ cf = (tf - f) - yf;
+ f = tf;
+ }
+ }
+ }
+
+ // If we're on an edge, then the scheme of adding up triangles as above
+ // doesn't work so well. We'll take care of these two nodes here.
+ if (onedge) {
+
+ // If we're on the convex hull, then the other nodes don't actually
+ // contribute anything, just the nodes for the edge we're on. The
+ // Voronoi "polygons" are infinite in extent.
+ if (onhull) {
+ double a = (hypot(targetx-x[edge[0]], targety-y[edge[0]]) /
+ hypot(x[edge[1]]-x[edge[0]], y[edge[1]]-y[edge[0]]));
+ return (1-a) * z[edge[0]] + a*z[edge[1]];
+ }
+
+ set<int> T(circumtri.begin(), circumtri.end());
+ vector<int> newedges0; // the two nodes that edge[0] still connect to
+ vector<int> newedges1; // the two nodes that edge[1] still connect to
+ set<int> alltri0; // all of the circumtriangle edge[0] participates in
+ set<int> alltri1; // all of the circumtriangle edge[1] participates in
+ for (it = circumtri.begin(); it != circumtri.end(); it++) {
+ for (int i=0; i<3; i++) {
+ int ti = INDEX3(this->neighbors, *it, i);
+ int j = EDGE0(i);
+ int k = EDGE1(i);
+ int q0 = INDEX3(this->nodes, *it, j);
+ int q1 = INDEX3(this->nodes, *it, k);
+
+ if ((q0 == edge[0]) || (q1 == edge[0])) alltri0.insert(*it);
+ if ((q0 == edge[1]) || (q1 == edge[1])) alltri1.insert(*it);
+
+ if (!T.count(ti)) {
+ // neighbor is not in the set of circumtriangles
+ if (q0 == edge[0]) newedges0.push_back(q1);
+ if (q1 == edge[0]) newedges0.push_back(q0);
+ if (q0 == edge[1]) newedges1.push_back(q1);
+ if (q1 == edge[1]) newedges1.push_back(q0);
+ }
+ }
+ }
+
+ set<int>::iterator sit;
+
+ double cx, cy;
+ ConvexPolygon poly0;
+ ConvexPolygon poly1;
+
+ if (edge[1] != newedges0[0]) {
+ circumcenter(this->x[edge[0]], this->y[edge[0]],
+ this->x[newedges0[0]], this->y[newedges0[0]],
+ targetx, targety,
+ cx, cy);
+ poly0.push(cx, cy);
+ }
+ if (edge[1] != newedges0[1]) {
+ circumcenter(this->x[edge[0]], this->y[edge[0]],
+ this->x[newedges0[1]], this->y[newedges0[1]],
+ targetx, targety,
+ cx, cy);
+ poly0.push(cx, cy);
+ }
+
+ if (edge[0] != newedges1[0]) {
+ circumcenter(this->x[edge[1]], this->y[edge[1]],
+ this->x[newedges1[0]], this->y[newedges1[0]],
+ targetx, targety,
+ cx, cy);
+ poly1.push(cx, cy);
+ }
+ if (edge[0] != newedges1[1]) {
+ circumcenter(this->x[edge[1]], this->y[edge[1]],
+ this->x[newedges1[1]], this->y[newedges1[1]],
+ targetx, targety,
+ cx, cy);
+ poly1.push(cx, cy);
+ }
+
+ for (sit = alltri0.begin(); sit != alltri0.end(); sit++) {
+ poly0.push(INDEX2(this->centers, *sit, 0),
+ INDEX2(this->centers, *sit, 1));
+ }
+ for (sit = alltri1.begin(); sit != alltri1.end(); sit++) {
+ poly1.push(INDEX2(this->centers, *sit, 0),
+ INDEX2(this->centers, *sit, 1));
+ }
+
+ double a0 = poly0.area();
+ double a1 = poly1.area();
+
+
+ f += a0*z[edge[0]];
+ A += a0;
+ f += a1*z[edge[1]];
+ A += a1;
+
+ // Anticlimactic, isn't it?
+ }
+
+ f /= A;
+ return f;
+}
+
+void NaturalNeighbors::interpolate_grid(double *z,
+ double x0, double x1, int xsteps,
+ double y0, double y1, int ysteps,
+ double *output,
+ double defvalue, int start_triangle)
+{
+ int i, ix, iy, rowtri, coltri, tri;
+ double dx, dy, targetx, targety;
+
+ dx = (x1 - x0) / (xsteps-1);
+ dy = (y1 - y0) / (ysteps-1);
+
+ rowtri = 0;
+ i = 0;
+ for (iy=0; iy<ysteps; iy++) {
+ targety = y0 + dy*iy;
+ rowtri = find_containing_triangle(x0, targety, rowtri);
+ tri = rowtri;
+ for (ix=0; ix<xsteps; ix++) {
+ targetx = x0 + dx*ix;
+ coltri = tri;
+ INDEXN(output, xsteps, iy, ix) = interpolate_one(z, targetx, targety,
+ defvalue, coltri);
+ if (coltri != -1) tri = coltri;
+ }
+ }
+}
+
+void NaturalNeighbors::interpolate_unstructured(double *z, int size,
+ double *intx, double *inty, double *output, double defvalue)
+{
+ int i, tri1, tri2;
+
+ tri1 = 0;
+ tri2 = 0;
+ for (i=0; i<size; i++) {
+ tri2 = tri1;
+ output[i] = interpolate_one(z, intx[i], inty[i], defvalue, tri2);
+ if (tri2 != -1) tri1 = tri2;
+ }
+}
Added: trunk/yt/raven/delaunay/natneighbors.h
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/natneighbors.h Fri May 2 16:44:22 2008
@@ -0,0 +1,33 @@
+#ifndef _NATNEIGHBORS_H
+#define _NATNEIGHBORS_H
+
+#include <list>
+using namespace std;
+
+class NaturalNeighbors
+{
+public:
+ NaturalNeighbors(int npoints, int ntriangles, double *x, double *y,
+ double *centers, int *nodes, int *neighbors);
+ ~NaturalNeighbors();
+
+ double interpolate_one(double *z, double targetx, double targety,
+ double defvalue, int &start_triangle);
+
+ void interpolate_grid(double *z,
+ double x0, double x1, int xsteps,
+ double y0, double y1, int ysteps,
+ double *output, double defvalue, int start_triangle);
+
+ void interpolate_unstructured(double *z, int size,
+ double *intx, double *inty, double *output, double defvalue);
+
+private:
+ int npoints, ntriangles;
+ double *x, *y, *centers, *radii2;
+ int *nodes, *neighbors;
+
+ int find_containing_triangle(double targetx, double targety, int start_triangle);
+};
+
+#endif // _NATNEIGHBORS_H
Added: trunk/yt/raven/delaunay/setup.py
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/setup.py Fri May 2 16:44:22 2008
@@ -0,0 +1,17 @@
+from numpy.distutils.core import setup
+from numpy.distutils.misc_util import Configuration
+
+def configuration(parent_package='', top_path=None):
+
+ config = Configuration('delaunay', parent_package, top_path)
+
+ config.add_extension("_delaunay",
+ sources=["_delaunay.cpp", "VoronoiDiagramGenerator.cpp",
+ "delaunay_utils.cpp", "natneighbors.cpp"],
+ include_dirs=['.'],
+ )
+
+ return config
+
+if __name__ == '__main__':
+ setup(**configuration(top_path='').todict())
Added: trunk/yt/raven/delaunay/testfuncs.py
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/testfuncs.py Fri May 2 16:44:22 2008
@@ -0,0 +1,461 @@
+"""Some test functions for bivariate interpolation.
+
+Most of these have been yoinked from ACM TOMS 792.
+http://netlib.org/toms/792
+"""
+
+import numpy as np
+from triangulate import Triangulation
+
+class TestData(dict):
+ def __init__(self, *args, **kwds):
+ dict.__init__(self, *args, **kwds)
+ self.__dict__ = self
+
+class TestDataSet(object):
+ def __init__(self, **kwds):
+ self.__dict__.update(kwds)
+
+data = TestData(
+franke100=TestDataSet(
+ x=np.array([ 0.0227035, 0.0539888, 0.0217008, 0.0175129, 0.0019029,
+ -0.0509685, 0.0395408, -0.0487061, 0.0315828, -0.0418785,
+ 0.1324189, 0.1090271, 0.1254439, 0.093454 , 0.0767578,
+ 0.1451874, 0.0626494, 0.1452734, 0.0958668, 0.0695559,
+ 0.2645602, 0.2391645, 0.208899 , 0.2767329, 0.1714726,
+ 0.2266781, 0.1909212, 0.1867647, 0.2304634, 0.2426219,
+ 0.3663168, 0.3857662, 0.3832392, 0.3179087, 0.3466321,
+ 0.3776591, 0.3873159, 0.3812917, 0.3795364, 0.2803515,
+ 0.4149771, 0.4277679, 0.420001 , 0.4663631, 0.4855658,
+ 0.4092026, 0.4792578, 0.4812279, 0.3977761, 0.4027321,
+ 0.5848691, 0.5730076, 0.6063893, 0.5013894, 0.5741311,
+ 0.6106955, 0.5990105, 0.5380621, 0.6096967, 0.5026188,
+ 0.6616928, 0.6427836, 0.6396475, 0.6703963, 0.7001181,
+ 0.633359 , 0.6908947, 0.6895638, 0.6718889, 0.6837675,
+ 0.7736939, 0.7635332, 0.7410424, 0.8258981, 0.7306034,
+ 0.8086609, 0.8214531, 0.729064 , 0.8076643, 0.8170951,
+ 0.8424572, 0.8684053, 0.8366923, 0.9418461, 0.8478122,
+ 0.8599583, 0.91757 , 0.8596328, 0.9279871, 0.8512805,
+ 1.044982 , 0.9670631, 0.9857884, 0.9676313, 1.0129299,
+ 0.965704 , 1.0019855, 1.0359297, 1.0414677, 0.9471506]),
+ y=np.array([-0.0310206, 0.1586742, 0.2576924, 0.3414014, 0.4943596,
+ 0.5782854, 0.6993418, 0.7470194, 0.9107649, 0.996289 ,
+ 0.050133 , 0.0918555, 0.2592973, 0.3381592, 0.4171125,
+ 0.5615563, 0.6552235, 0.7524066, 0.9146523, 0.9632421,
+ 0.0292939, 0.0602303, 0.2668783, 0.3696044, 0.4801738,
+ 0.5940595, 0.6878797, 0.8185576, 0.9046507, 0.9805412,
+ 0.0396955, 0.0684484, 0.2389548, 0.3124129, 0.4902989,
+ 0.5199303, 0.6445227, 0.8203789, 0.8938079, 0.9711719,
+ -0.0284618, 0.1560965, 0.2262471, 0.3175094, 0.3891417,
+ 0.5084949, 0.6324247, 0.7511007, 0.8489712, 0.9978728,
+ -0.0271948, 0.127243 , 0.2709269, 0.3477728, 0.4259422,
+ 0.6084711, 0.6733781, 0.7235242, 0.9242411, 1.0308762,
+ 0.0255959, 0.0707835, 0.2008336, 0.3259843, 0.4890704,
+ 0.5096324, 0.669788 , 0.7759569, 0.9366096, 1.0064516,
+ 0.0285374, 0.1021403, 0.1936581, 0.3235775, 0.4714228,
+ 0.6091595, 0.6685053, 0.8022808, 0.847679 , 1.0512371,
+ 0.0380499, 0.0902048, 0.2083092, 0.3318491, 0.4335632,
+ 0.5910139, 0.6307383, 0.8144841, 0.904231 , 0.969603 ,
+ -0.01209 , 0.1334114, 0.2695844, 0.3795281, 0.4396054,
+ 0.5044425, 0.6941519, 0.7459923, 0.8682081, 0.9801409])),
+franke33=TestDataSet(
+ x=np.array([ 5.00000000e-02, 0.00000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 1.00000000e-01, 1.00000000e-01,
+ 1.50000000e-01, 2.00000000e-01, 2.50000000e-01,
+ 3.00000000e-01, 3.50000000e-01, 5.00000000e-01,
+ 5.00000000e-01, 5.50000000e-01, 6.00000000e-01,
+ 6.00000000e-01, 6.00000000e-01, 6.50000000e-01,
+ 7.00000000e-01, 7.00000000e-01, 7.00000000e-01,
+ 7.50000000e-01, 7.50000000e-01, 7.50000000e-01,
+ 8.00000000e-01, 8.00000000e-01, 8.50000000e-01,
+ 9.00000000e-01, 9.00000000e-01, 9.50000000e-01,
+ 1.00000000e+00, 1.00000000e+00, 1.00000000e+00]),
+ y=np.array([ 4.50000000e-01, 5.00000000e-01, 1.00000000e+00,
+ 0.00000000e+00, 1.50000000e-01, 7.50000000e-01,
+ 3.00000000e-01, 1.00000000e-01, 2.00000000e-01,
+ 3.50000000e-01, 8.50000000e-01, 0.00000000e+00,
+ 1.00000000e+00, 9.50000000e-01, 2.50000000e-01,
+ 6.50000000e-01, 8.50000000e-01, 7.00000000e-01,
+ 2.00000000e-01, 6.50000000e-01, 9.00000000e-01,
+ 1.00000000e-01, 3.50000000e-01, 8.50000000e-01,
+ 4.00000000e-01, 6.50000000e-01, 2.50000000e-01,
+ 3.50000000e-01, 8.00000000e-01, 9.00000000e-01,
+ 0.00000000e+00, 5.00000000e-01, 1.00000000e+00])),
+lawson25=TestDataSet(
+ x=np.array([ 0.1375, 0.9125, 0.7125, 0.225 , -0.05 , 0.475 , 0.05 ,
+ 0.45 , 1.0875, 0.5375, -0.0375, 0.1875, 0.7125, 0.85 ,
+ 0.7 , 0.275 , 0.45 , 0.8125, 0.45 , 1. , 0.5 ,
+ 0.1875, 0.5875, 1.05 , 0.1 ]),
+ y=np.array([ 0.975 , 0.9875 , 0.7625 , 0.8375 , 0.4125 , 0.6375 ,
+ -0.05 , 1.0375 , 0.55 , 0.8 , 0.75 , 0.575 ,
+ 0.55 , 0.4375 , 0.3125 , 0.425 , 0.2875 , 0.1875 ,
+ -0.0375 , 0.2625 , 0.4625 , 0.2625 , 0.125 , -0.06125, 0.1125 ])),
+random100=TestDataSet(
+ x=np.array([ 0.0096326, 0.0216348, 0.029836 , 0.0417447, 0.0470462,
+ 0.0562965, 0.0646857, 0.0740377, 0.0873907, 0.0934832,
+ 0.1032216, 0.1110176, 0.1181193, 0.1251704, 0.132733 ,
+ 0.1439536, 0.1564861, 0.1651043, 0.1786039, 0.1886405,
+ 0.2016706, 0.2099886, 0.2147003, 0.2204141, 0.2343715,
+ 0.240966 , 0.252774 , 0.2570839, 0.2733365, 0.2853833,
+ 0.2901755, 0.2964854, 0.3019725, 0.3125695, 0.3307163,
+ 0.3378504, 0.3439061, 0.3529922, 0.3635507, 0.3766172,
+ 0.3822429, 0.3869838, 0.3973137, 0.4170708, 0.4255588,
+ 0.4299218, 0.4372839, 0.4705033, 0.4736655, 0.4879299,
+ 0.494026 , 0.5055324, 0.5162593, 0.5219219, 0.5348529,
+ 0.5483213, 0.5569571, 0.5638611, 0.5784908, 0.586395 ,
+ 0.5929148, 0.5987839, 0.6117561, 0.6252296, 0.6331381,
+ 0.6399048, 0.6488972, 0.6558537, 0.6677405, 0.6814074,
+ 0.6887812, 0.6940896, 0.7061687, 0.7160957, 0.7317445,
+ 0.7370798, 0.746203 , 0.7566957, 0.7699998, 0.7879347,
+ 0.7944014, 0.8164468, 0.8192794, 0.8368405, 0.8500993,
+ 0.8588255, 0.8646496, 0.8792329, 0.8837536, 0.8900077,
+ 0.8969894, 0.9044917, 0.9083947, 0.9203972, 0.9347906,
+ 0.9434519, 0.9490328, 0.9569571, 0.9772067, 0.9983493]),
+ y=np.array([ 0.3083158, 0.2450434, 0.8613847, 0.0977864, 0.3648355,
+ 0.7156339, 0.5311312, 0.9755672, 0.1781117, 0.5452797,
+ 0.1603881, 0.7837139, 0.9982015, 0.6910589, 0.104958 ,
+ 0.8184662, 0.7086405, 0.4456593, 0.1178342, 0.3189021,
+ 0.9668446, 0.7571834, 0.2016598, 0.3232444, 0.4368583,
+ 0.8907869, 0.064726 , 0.5692618, 0.2947027, 0.4332426,
+ 0.3347464, 0.7436284, 0.1066265, 0.8845357, 0.515873 ,
+ 0.9425637, 0.4799701, 0.1783069, 0.114676 , 0.8225797,
+ 0.2270688, 0.4073598, 0.887508 , 0.7631616, 0.9972804,
+ 0.4959884, 0.3410421, 0.249812 , 0.6409007, 0.105869 ,
+ 0.5411969, 0.0089792, 0.8784268, 0.5515874, 0.4038952,
+ 0.1654023, 0.2965158, 0.3660356, 0.0366554, 0.950242 ,
+ 0.2638101, 0.9277386, 0.5377694, 0.7374676, 0.4674627,
+ 0.9186109, 0.0416884, 0.1291029, 0.6763676, 0.8444238,
+ 0.3273328, 0.1893879, 0.0645923, 0.0180147, 0.8904992,
+ 0.4160648, 0.4688995, 0.2174508, 0.5734231, 0.8853319,
+ 0.8018436, 0.6388941, 0.8931002, 0.1000558, 0.2789506,
+ 0.9082948, 0.3259159, 0.8318747, 0.0508513, 0.970845 ,
+ 0.5120548, 0.2859716, 0.9581641, 0.6183429, 0.3779934,
+ 0.4010423, 0.9478657, 0.7425486, 0.8883287, 0.549675 ])),
+uniform9=TestDataSet(
+ x=np.array([ 1.25000000e-01, 0.00000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
+ 0.00000000e+00, 1.25000000e-01, 1.25000000e-01,
+ 1.25000000e-01, 1.25000000e-01, 1.25000000e-01,
+ 1.25000000e-01, 1.25000000e-01, 1.25000000e-01,
+ 2.50000000e-01, 2.50000000e-01, 2.50000000e-01,
+ 2.50000000e-01, 2.50000000e-01, 2.50000000e-01,
+ 2.50000000e-01, 2.50000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 3.75000000e-01, 3.75000000e-01,
+ 3.75000000e-01, 3.75000000e-01, 3.75000000e-01,
+ 3.75000000e-01, 3.75000000e-01, 3.75000000e-01,
+ 5.00000000e-01, 5.00000000e-01, 5.00000000e-01,
+ 5.00000000e-01, 5.00000000e-01, 5.00000000e-01,
+ 5.00000000e-01, 5.00000000e-01, 5.00000000e-01,
+ 6.25000000e-01, 6.25000000e-01, 6.25000000e-01,
+ 6.25000000e-01, 6.25000000e-01, 6.25000000e-01,
+ 6.25000000e-01, 6.25000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 7.50000000e-01, 7.50000000e-01,
+ 7.50000000e-01, 7.50000000e-01, 7.50000000e-01,
+ 7.50000000e-01, 7.50000000e-01, 7.50000000e-01,
+ 8.75000000e-01, 8.75000000e-01, 8.75000000e-01,
+ 8.75000000e-01, 8.75000000e-01, 8.75000000e-01,
+ 8.75000000e-01, 8.75000000e-01, 8.75000000e-01,
+ 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,
+ 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,
+ 1.00000000e+00, 1.00000000e+00, 1.00000000e+00]),
+ y=np.array([ 0.00000000e+00, 1.25000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 5.00000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 8.75000000e-01, 1.00000000e+00,
+ 0.00000000e+00, 1.25000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 5.00000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 8.75000000e-01, 1.00000000e+00,
+ 0.00000000e+00, 1.25000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 5.00000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 8.75000000e-01, 1.00000000e+00,
+ 0.00000000e+00, 1.25000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 5.00000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 8.75000000e-01, 1.00000000e+00,
+ 0.00000000e+00, 1.25000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 5.00000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 8.75000000e-01, 1.00000000e+00,
+ 0.00000000e+00, 1.25000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 5.00000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 8.75000000e-01, 1.00000000e+00,
+ 0.00000000e+00, 1.25000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 5.00000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 8.75000000e-01, 1.00000000e+00,
+ 0.00000000e+00, 1.25000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 5.00000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 8.75000000e-01, 1.00000000e+00,
+ 0.00000000e+00, 1.25000000e-01, 2.50000000e-01,
+ 3.75000000e-01, 5.00000000e-01, 6.25000000e-01,
+ 7.50000000e-01, 8.75000000e-01, 1.00000000e+00])),
+)
+
+
+
+
+
+def constant(x, y):
+ return np.ones(x.shape, x.dtype)
+constant.title = 'Constant'
+
+def xramp(x, y):
+ return x
+xramp.title = 'X Ramp'
+
+def yramp(x, y):
+ return y
+yramp.title = 'Y Ramp'
+
+def exponential(x, y):
+ x = x*9
+ y = y*9
+ x1 = x+1.0
+ x2 = x-2.0
+ x4 = x-4.0
+ x7 = x-7.0
+ y1 = x+1.0
+ y2 = y-2.0
+ y3 = y-3.0
+ y7 = y-7.0
+ f = (0.75 * np.exp(-(x2*x2+y2*y2)/4.0) +
+ 0.75 * np.exp(-x1*x1/49.0 - y1/10.0) +
+ 0.5 * np.exp(-(x7*x7 + y3*y3)/4.0) -
+ 0.2 * np.exp(-x4*x4 -y7*y7))
+ return f
+exponential.title = 'Exponential and Some Gaussians'
+
+def cliff(x, y):
+ f = np.tanh(9.0*(y-x) + 1.0)/9.0
+ return f
+cliff.title = 'Cliff'
+
+def saddle(x, y):
+ f = (1.25 + np.cos(5.4*y))/(6.0 + 6.0*(3*x-1.0)**2)
+ return f
+saddle.title = 'Saddle'
+
+def gentle(x, y):
+ f = np.exp(-5.0625*((x-0.5)**2+(y-0.5)**2))/3.0
+ return f
+gentle.title = 'Gentle Peak'
+
+def steep(x, y):
+ f = np.exp(-20.25*((x-0.5)**2+(y-0.5)**2))/3.0
+ return f
+steep.title = 'Steep Peak'
+
+def sphere(x, y):
+ circle = 64-81*((x-0.5)**2 + (y-0.5)**2)
+ f = np.where(circle >= 0, np.sqrt(np.clip(circle,0,100)) - 0.5, 0.0)
+ return f
+sphere.title = 'Sphere'
+
+def trig(x, y):
+ f = 2.0*np.cos(10.0*x)*np.sin(10.0*y) + np.sin(10.0*x*y)
+ return f
+trig.title = 'Cosines and Sines'
+
+def gauss(x, y):
+ x = 5.0-10.0*x
+ y = 5.0-10.0*y
+ g1 = np.exp(-x*x/2)
+ g2 = np.exp(-y*y/2)
+ f = g1 + 0.75*g2*(1 + g1)
+ return f
+gauss.title = 'Gaussian Peak and Gaussian Ridges'
+
+def cloverleaf(x, y):
+ ex = np.exp((10.0-20.0*x)/3.0)
+ ey = np.exp((10.0-20.0*y)/3.0)
+ logitx = 1.0/(1.0+ex)
+ logity = 1.0/(1.0+ey)
+ f = (((20.0/3.0)**3 * ex*ey)**2 * (logitx*logity)**5 *
+ (ex-2.0*logitx)*(ey-2.0*logity))
+ return f
+cloverleaf.title = 'Cloverleaf'
+
+def cosine_peak(x, y):
+ circle = np.hypot(80*x-40.0, 90*y-45.)
+ f = np.exp(-0.04*circle) * np.cos(0.15*circle)
+ return f
+cosine_peak.title = 'Cosine Peak'
+
+allfuncs = [exponential, cliff, saddle, gentle, steep, sphere, trig, gauss, cloverleaf, cosine_peak]
+
+
+class LinearTester(object):
+ name = 'Linear'
+ def __init__(self, xrange=(0.0, 1.0), yrange=(0.0, 1.0), nrange=101, npoints=250):
+ self.xrange = xrange
+ self.yrange = yrange
+ self.nrange = nrange
+ self.npoints = npoints
+
+ rng = np.random.RandomState(1234567890)
+ self.x = rng.uniform(xrange[0], xrange[1], size=npoints)
+ self.y = rng.uniform(yrange[0], yrange[1], size=npoints)
+ self.tri = Triangulation(self.x, self.y)
+
+ def replace_data(self, dataset):
+ self.x = dataset.x
+ self.y = dataset.y
+ self.tri = Triangulation(self.x, self.y)
+
+ def interpolator(self, func):
+ z = func(self.x, self.y)
+ return self.tri.linear_extrapolator(z, bbox=self.xrange+self.yrange)
+
+ def plot(self, func, interp=True, plotter='imshow'):
+ import matplotlib as mpl
+ from matplotlib import pylab as pl
+ if interp:
+ lpi = self.interpolator(func)
+ z = lpi[self.yrange[0]:self.yrange[1]:complex(0,self.nrange),
+ self.xrange[0]:self.xrange[1]:complex(0,self.nrange)]
+ else:
+ y, x = np.mgrid[self.yrange[0]:self.yrange[1]:complex(0,self.nrange),
+ self.xrange[0]:self.xrange[1]:complex(0,self.nrange)]
+ z = func(x, y)
+
+ z = np.where(np.isinf(z), 0.0, z)
+
+ extent = (self.xrange[0], self.xrange[1],
+ self.yrange[0], self.yrange[1])
+ pl.ioff()
+ pl.clf()
+ pl.hot() # Some like it hot
+ if plotter == 'imshow':
+ pl.imshow(np.nan_to_num(z), interpolation='nearest', extent=extent, origin='lower')
+ elif plotter == 'contour':
+ Y, X = np.ogrid[self.yrange[0]:self.yrange[1]:complex(0,self.nrange),
+ self.xrange[0]:self.xrange[1]:complex(0,self.nrange)]
+ pl.contour(np.ravel(X), np.ravel(Y), z, 20)
+ x = self.x
+ y = self.y
+ lc = mpl.collections.LineCollection(np.array([((x[i], y[i]), (x[j], y[j]))
+ for i, j in self.tri.edge_db]), colors=[(0,0,0,0.2)])
+ ax = pl.gca()
+ ax.add_collection(lc)
+
+ if interp:
+ title = '%s Interpolant' % self.name
+ else:
+ title = 'Reference'
+ if hasattr(func, 'title'):
+ pl.title('%s: %s' % (func.title, title))
+ else:
+ pl.title(title)
+
+ pl.show()
+ pl.ion()
+
+class NNTester(LinearTester):
+ name = 'Natural Neighbors'
+ def interpolator(self, func):
+ z = func(self.x, self.y)
+ return self.tri.nn_extrapolator(z, bbox=self.xrange+self.yrange)
+
+def plotallfuncs(allfuncs=allfuncs):
+ from matplotlib import pylab as pl
+ pl.ioff()
+ nnt = NNTester(npoints=1000)
+ lpt = LinearTester(npoints=1000)
+ for func in allfuncs:
+ print func.title
+ nnt.plot(func, interp=False, plotter='imshow')
+ pl.savefig('%s-ref-img.png' % func.func_name)
+ nnt.plot(func, interp=True, plotter='imshow')
+ pl.savefig('%s-nn-img.png' % func.func_name)
+ lpt.plot(func, interp=True, plotter='imshow')
+ pl.savefig('%s-lin-img.png' % func.func_name)
+ nnt.plot(func, interp=False, plotter='contour')
+ pl.savefig('%s-ref-con.png' % func.func_name)
+ nnt.plot(func, interp=True, plotter='contour')
+ pl.savefig('%s-nn-con.png' % func.func_name)
+ lpt.plot(func, interp=True, plotter='contour')
+ pl.savefig('%s-lin-con.png' % func.func_name)
+ pl.ion()
+
+def plot_dt(tri, colors=None):
+ import matplotlib as mpl
+ from matplotlib import pylab as pl
+ if colors is None:
+ colors = [(0,0,0,0.2)]
+ lc = mpl.collections.LineCollection(np.array([((tri.x[i], tri.y[i]), (tri.x[j], tri.y[j]))
+ for i, j in tri.edge_db]), colors=colors)
+ ax = pl.gca()
+ ax.add_collection(lc)
+ pl.draw_if_interactive()
+
+def plot_vo(tri, colors=None):
+ import matplotlib as mpl
+ from matplotlib import pylab as pl
+ if colors is None:
+ colors = [(0,1,0,0.2)]
+ lc = mpl.collections.LineCollection(np.array(
+ [(tri.circumcenters[i], tri.circumcenters[j])
+ for i in xrange(len(tri.circumcenters))
+ for j in tri.triangle_neighbors[i] if j != -1]),
+ colors=colors)
+ ax = pl.gca()
+ ax.add_collection(lc)
+ pl.draw_if_interactive()
+
+def plot_cc(tri, edgecolor=None):
+ import matplotlib as mpl
+ from matplotlib import pylab as pl
+ if edgecolor is None:
+ edgecolor = (0,0,1,0.2)
+ dxy = (np.array([(tri.x[i], tri.y[i]) for i,j,k in tri.triangle_nodes])
+ - tri.circumcenters)
+ r = np.hypot(dxy[:,0], dxy[:,1])
+ ax = pl.gca()
+ for i in xrange(len(r)):
+ p = mpl.patches.Circle(tri.circumcenters[i], r[i], resolution=100, edgecolor=edgecolor,
+ facecolor=(1,1,1,0), linewidth=0.2)
+ ax.add_patch(p)
+ pl.draw_if_interactive()
+
+def quality(func, mesh, interpolator='nn', n=33):
+ """Compute a quality factor (the quantity r**2 from TOMS792).
+
+ interpolator must be in ('linear', 'nn').
+ """
+ fz = func(mesh.x, mesh.y)
+ tri = Triangulation(mesh.x, mesh.y)
+ intp = getattr(tri, interpolator+'_extrapolator')(fz, bbox=(0.,1.,0.,1.))
+ Y, X = np.mgrid[0:1:complex(0,n),0:1:complex(0,n)]
+ Z = func(X, Y)
+ iz = intp[0:1:complex(0,n),0:1:complex(0,n)]
+ #nans = np.isnan(iz)
+ #numgood = n*n - np.sum(np.array(nans.flat, np.int32))
+ numgood = n*n
+
+ SE = (Z - iz)**2
+ SSE = np.sum(SE.flat)
+ meanZ = np.sum(Z.flat) / numgood
+ SM = (Z - meanZ)**2
+ SSM = np.sum(SM.flat)
+
+
+ r2 = 1.0 - SSE/SSM
+ print func.func_name, r2, SSE, SSM, numgood
+ return r2
+
+def allquality(interpolator='nn', allfuncs=allfuncs, data=data, n=33):
+ results = {}
+ kv = data.items()
+ kv.sort()
+ for name, mesh in kv:
+ reslist = results.setdefault(name, [])
+ for func in allfuncs:
+ reslist.append(quality(func, mesh, interpolator, n))
+ return results
+
+
+def funky():
+ x0 = np.array([0.25, 0.3, 0.5, 0.6, 0.6])
+ y0 = np.array([0.2, 0.35, 0.0, 0.25, 0.65])
+ tx = 0.46
+ ty = 0.23
+ t0 = Triangulation(x0, y0)
+ t1 = Triangulation(np.hstack((x0, [tx])), np.hstack((y0, [ty])))
+ return t0, t1
Added: trunk/yt/raven/delaunay/triangulate.py
==============================================================================
--- (empty file)
+++ trunk/yt/raven/delaunay/triangulate.py Fri May 2 16:44:22 2008
@@ -0,0 +1,211 @@
+import warnings
+# 2.3 compatibility
+try:
+ set
+except NameError:
+ import sets
+ set = sets.Set
+
+import numpy as np
+
+from _delaunay import delaunay
+from interpolate import LinearInterpolator, NNInterpolator
+
+__all__ = ['Triangulation', 'DuplicatePointWarning']
+
+
+class DuplicatePointWarning(RuntimeWarning):
+ """Duplicate points were passed in to the triangulation routine.
+ """
+
+
+class Triangulation(object):
+ """A Delaunay triangulation of points in a plane.
+
+ Triangulation(x, y)
+ x, y -- the coordinates of the points as 1-D arrays of floats
+
+ Let us make the following definitions:
+ npoints = number of points input
+ nedges = number of edges in the triangulation
+ ntriangles = number of triangles in the triangulation
+
+ point_id = an integer identifying a particular point (specifically, an
+ index into x and y), range(0, npoints)
+ edge_id = an integer identifying a particular edge, range(0, nedges)
+ triangle_id = an integer identifying a particular triangle
+ range(0, ntriangles)
+
+ Attributes: (all should be treated as read-only to maintain consistency)
+ x, y -- the coordinates of the points as 1-D arrays of floats.
+
+ circumcenters -- (ntriangles, 2) array of floats giving the (x,y)
+ coordinates of the circumcenters of each triangle (indexed by a
+ triangle_id).
+
+ edge_db -- (nedges, 2) array of point_id's giving the points forming
+ each edge in no particular order; indexed by an edge_id.
+
+ triangle_nodes -- (ntriangles, 3) array of point_id's giving the points
+ forming each triangle in counter-clockwise order; indexed by a
+ triangle_id.
+
+ triangle_neighbors -- (ntriangles, 3) array of triangle_id's giving the
+ neighboring triangle; indexed by a triangle_id.
+
+ The value can also be -1 meaning that that edge is on the convex hull of
+ the points and there is no neighbor on that edge. The values are ordered
+ such that triangle_neighbors[tri, i] corresponds with the edge
+ *opposite* triangle_nodes[tri, i]. As such, these neighbors are also in
+ counter-clockwise order.
+
+ hull -- list of point_id's giving the nodes which form the convex hull
+ of the point set. This list is sorted in counter-clockwise order.
+ """
+ def __init__(self, x, y):
+ self.x = np.asarray(x, dtype=np.float64)
+ self.y = np.asarray(y, dtype=np.float64)
+
+ if self.x.shape != self.y.shape or len(self.x.shape) != 1:
+ raise ValueError("x,y must be equal-length 1-D arrays")
+
+ self.old_shape = self.x.shape
+ j_unique = self._collapse_duplicate_points()
+
+ if j_unique.shape != self.x.shape:
+ warnings.warn(
+ "Input data contains duplicate x,y points; some values are ignored.",
+ DuplicatePointWarning,
+ )
+ self.j_unique = j_unique
+ self.x = self.x[self.j_unique]
+ self.y = self.y[self.j_unique]
+ else:
+ self.j_unique = None
+
+
+ self.circumcenters, self.edge_db, self.triangle_nodes, \
+ self.triangle_neighbors = delaunay(self.x, self.y)
+
+ self.hull = self._compute_convex_hull()
+
+ def _collapse_duplicate_points(self):
+ """Generate index array that picks out unique x,y points.
+
+ This appears to be required by the underlying delaunay triangulation
+ code.
+ """
+ # Find the indices of the unique entries
+ j_sorted = np.lexsort(keys=(self.x, self.y))
+ mask_unique = np.hstack([
+ True,
+ (np.diff(self.x[j_sorted]) != 0) | (np.diff(self.y[j_sorted]) != 0),
+ ])
+ return j_sorted[mask_unique]
+
+ def _compute_convex_hull(self):
+ """Extract the convex hull from the triangulation information.
+
+ The output will be a list of point_id's in counter-clockwise order
+ forming the convex hull of the data set.
+ """
+ border = (self.triangle_neighbors == -1)
+
+ edges = {}
+ edges.update(dict(zip(self.triangle_nodes[border[:,0]][:,1],
+ self.triangle_nodes[border[:,0]][:,2])))
+ edges.update(dict(zip(self.triangle_nodes[border[:,1]][:,2],
+ self.triangle_nodes[border[:,1]][:,0])))
+ edges.update(dict(zip(self.triangle_nodes[border[:,2]][:,0],
+ self.triangle_nodes[border[:,2]][:,1])))
+
+ # Take an arbitrary starting point and its subsequent node
+ hull = list(edges.popitem())
+ while edges:
+ hull.append(edges.pop(hull[-1]))
+
+ # hull[-1] == hull[0], so remove hull[-1]
+ hull.pop()
+
+ return hull
+
+ def linear_interpolator(self, z, default_value=np.nan):
+ """Get an object which can interpolate within the convex hull by
+ assigning a plane to each triangle.
+
+ z -- an array of floats giving the known function values at each point
+ in the triangulation.
+ """
+ z = np.asarray(z, dtype=np.float64)
+ if z.shape != self.old_shape:
+ raise ValueError("z must be the same shape as x and y")
+ if self.j_unique is not None:
+ z = z[self.j_unique]
+
+ return LinearInterpolator(self, z, default_value)
+
+ def nn_interpolator(self, z, default_value=np.nan):
+ """Get an object which can interpolate within the convex hull by
+ the natural neighbors method.
+
+ z -- an array of floats giving the known function values at each point
+ in the triangulation.
+ """
+ z = np.asarray(z, dtype=np.float64)
+ if z.shape != self.old_shape:
+ raise ValueError("z must be the same shape as x and y")
+ if self.j_unique is not None:
+ z = z[self.j_unique]
+
+ return NNInterpolator(self, z, default_value)
+
+ def prep_extrapolator(self, z, bbox=None):
+ if bbox is None:
+ bbox = (self.x[0], self.x[0], self.y[0], self.y[0])
+ minx, maxx, miny, maxy = np.asarray(bbox, np.float64)
+ minx = min(minx, np.minimum.reduce(self.x))
+ miny = min(miny, np.minimum.reduce(self.y))
+ maxx = max(maxx, np.maximum.reduce(self.x))
+ maxy = max(maxy, np.maximum.reduce(self.y))
+ M = max((maxx-minx)/2, (maxy-miny)/2)
+ midx = (minx + maxx)/2.0
+ midy = (miny + maxy)/2.0
+
+ xp, yp= np.array([[midx+3*M, midx, midx-3*M],
+ [midy, midy+3*M, midy-3*M]])
+ x1 = np.hstack((self.x, xp))
+ y1 = np.hstack((self.y, yp))
+ newtri = self.__class__(x1, y1)
+
+ # do a least-squares fit to a plane to make pseudo-data
+ xy1 = np.ones((len(self.x), 3), np.float64)
+ xy1[:,0] = self.x
+ xy1[:,1] = self.y
+ from numpy.dual import lstsq
+ c, res, rank, s = lstsq(xy1, z)
+ zp = np.hstack((z, xp*c[0] + yp*c[1] + c[2]))
+
+ return newtri, zp
+
+ def nn_extrapolator(self, z, bbox=None, default_value=np.nan):
+ newtri, zp = self.prep_extrapolator(z, bbox)
+ return newtri.nn_interpolator(zp, default_value)
+
+ def linear_extrapolator(self, z, bbox=None, default_value=np.nan):
+ newtri, zp = self.prep_extrapolator(z, bbox)
+ return newtri.linear_interpolator(zp, default_value)
+
+ def node_graph(self):
+ """Return a graph of node_id's pointing to node_id's.
+
+ The arcs of the graph correspond to the edges in the triangulation.
+
+ {node_id: set([node_id, ...]), ...}
+ """
+ g = {}
+ for i, j in self.edge_db:
+ s = g.setdefault(i, set())
+ s.add(j)
+ s = g.setdefault(j, set())
+ s.add(i)
+ return g
Modified: trunk/yt/raven/setup.py
==============================================================================
--- trunk/yt/raven/setup.py (original)
+++ trunk/yt/raven/setup.py Fri May 2 16:44:22 2008
@@ -4,6 +4,7 @@
from numpy.distutils.misc_util import Configuration
config = Configuration('raven',parent_package,top_path)
config.add_subpackage("deliveration")
+ config.add_subpackage("delaunay") # From SciPy, primarily written by Robert Kern
config.make_config_py() # installs __config__.py
config.add_extension("_MPL", "_MPL.c", libraries=["m"])
return config
More information about the yt-svn
mailing list