Nils Wagner wrote:

John Hunter wrote:

...

...

...

...

...

> "Nils" == Nils Wagner <nwagner@iam.uni-stuttgart.de> writes: > Nils> Hi all, I have installed delaunay from the sandbox. How can Nils> I triangulate L-shaped domains ?

Nils> My first try is somewhat unsatisfactory (delaunay.png) ? Nils> How can I remove the unwanted triangles down to the right ?

delaunay assumes a convex shape, which your domain is not. I think you'll need a more sophisticated mesh algorithm.

A short note in the docstring would be very helpful. **delaunay assumes a convex shape**

Actually, it doesn't "assume" any shape at all. It computes a Delaunay triangulation of a set of points irrespective of any boundary edges that you might have wanted. That's what an (unqualified) Delaunay triangulation is. Docstrings are not the place for tutorials on basic concepts. Google works quite well. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco

...

...

...

...

"Rob" == Rob Hetland <hetland@tamu.edu> writes:

Rob> I have a feeling that there *might* be a way to tease out Rob> which triangles are inside the set of points you give, and Rob> which are outside. However, I have to say, that I thought Rob> about this for a while, and then gave up. Rob> More generally: Is there a good utility for finding points Rob> inside an arbitrary polygon? I, for one, could really use Rob> such a utility. Here is something I wrote for matplotlib -- someone more clever than I might be able to do it w/o the loop over the vertices by making better use of numpy .... import matplotlib.numerix as nx def inside_poly(points, verts): """ points is a sequence of x,y points verts is a sequence of x,y vertices of a poygon return value is a sequence on indices into points for the points that are inside the polygon """ xys = nx.asarray(points) Nxy = xys.shape[0] Nv = len(verts) def angle(x1, y1, x2, y2): twopi = 2*nx.pi theta1 = nx.arctan2(y1, x1) theta2 = nx.arctan2(y2, x2) dtheta = theta2-theta1 d = dtheta%twopi d = nx.where(nx.less(d, 0), twopi + d, d) return nx.where(nx.greater(d,nx.pi), d-twopi, d) angles = nx.zeros((Nxy,), nx.Float) x1 = nx.zeros((Nxy,), nx.Float) y1 = nx.zeros((Nxy,), nx.Float) x2 = nx.zeros((Nxy,), nx.Float) y2 = nx.zeros((Nxy,), nx.Float) x = xys[:,0] y = xys[:,1] for i in range(Nv): thisx, thisy = verts[i] x1 = thisx - x y1 = thisy - y thisx, thisy = verts[(i+1)%Nv] x2 = thisx - x y2 = thisy - y a = angle(x1, y1, x2, y2) angles += a return nx.nonzero(nx.greater_equal(nx.absolute(angles), nx.pi))

On 7/7/06, Rob Hetland <hetland@tamu.edu> wrote:

I have a feeling that there *might* be a way to tease out which triangles are inside the set of points you give, and which are outside. However, I have to say, that I thought about this for a while, and then gave up.

More generally: Is there a good utility for finding points inside an arbitrary polygon? I, for one, could really use such a utility.

The easiest thing is to do a ray-cast to infinity and count how many edges of the polygon you intersect. If you cross an odd number then the point was inside, if even then it was outside. So to implement it basically you only need a 2D ray-edge intersection routine, and the simplest thing is to just check the ray against every edge of the polygon one by one. If you have a lot of edges you can speed up the search in a variety of ways, like by using a spatial data strucure for your edges, e.g. a quadtree, k-d tree, bsp tree, etc. --bb

A very simple implementation (not heavily tested) of this method is attached. Again, unfortunately, not vectorised.

Cheers Stéfan <poly.py>

[I sent this message a while ago, but it appears to have been lost... so again:] I modified Stéfan's code (he may not recognize it anymore -- now it's ugly), so that the area and centroid methods are vectorized (using numpy arrays). I also simplified the inside method (was call contains, but I think inside is more common). Finally, I renamed the entire class from Poly to Polygon to prevent confusion with polynomial classes. Everything seems to work, but I have not had a chance to test it extensively. The inside class is slow when there are many points. The code available here: http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html seems to work great. Unfortunately, I am not very good at wrapping C into python, and I also not good at putting compiled code into distutils. If anyone else has the skills and is so motivated, could we make a small class (with a compiled inside method) and throw it in the scipy sandbox? -Rob  ---- Rob Hetland, Assistant Professor Dept. of Oceanography, Texas A&M University http://pong.tamue.edu/~rob phone: 979-458-0096, fax: 979-845-6331

If your really set on making this work, you should probably have a look at the wonderful project of wrapping cgal to python: http://cgal-python.gforge.inria.fr/ A very, very exciting project in terms of computational geometry! Cgal has a constrained version of the Delaunay triangulation, not sure if its as of yet wrapped in cgal-python... The project is moving rapidly though... -jelle -----Original Message----- From: scipy-user-bounces@scipy.org [mailto:scipy-user-bounces@scipy.org] On Behalf Of Robert Kern Sent: Friday, July 07, 2006 8:21 PM To: SciPy Users List Subject: Re: [SciPy-user] Triangulation of L-shaped domains Rob Hetland wrote:

I have a feeling that there *might* be a way to tease out which triangles are inside the set of points you give, and which are outside. However, I have to say, that I thought about this for a while, and then gave up.

It wouldn't necessarily work. The Delaunay triangulation performed by scipy.sandbox.delaunay is not constrained to require the edges of your non-convex boundary. But supposing you did have a constrained Delaunay triangulation, you could pick a triangle that is known to be outside and then push its neighbors onto a stack if the edge that joins them is not a boundary edge. Pop a triangle from the stack and repeat. (More or less; there are some details you need to fill in). You will have to put multiple "seed" triangles to initialize the stack if there are multiple "outside" regions (including holes!).

More generally: Is there a good utility for finding points inside an arbitrary polygon? I, for one, could really use such a utility.

Raycast with simulation of simplicity to handle degeneracy is probably your best bet. http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html John's "add up the angles" approach is not really a good one. I frequently find it referred to in the literature as "the worst thing you could possibly do." :-) -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user

...

...

...

...

"Robert" == Robert Kern <robert.kern@gmail.com> writes:

Robert> Raycast with simulation of simplicity to handle degeneracy Robert> is probably your best bet. Robert> Robert> http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html Robert> John's "add up the angles" approach is not really a good Robert> one. I frequently find it referred to in the literature as Robert> "the worst thing you could possibly do." :-) OK, I could only stand Robert's taunt for so long. I added the pnpoly routine to matplotlib svn extension code: there is a function to test whether a a single x, y point is inside a polygon (matplotlib.nxutils.pnpoly) and a function which takes a list of points and returns a mask for the points inside (matplotlib.nxutils.points_inside_poly). When profiling against my original "worst thing you could possibly do" implementation, this new version is 15-35 times faster, in addition to being a better algorithm in terms of handling the degenerate cases. It is also considerably faster than the pure numpy pnpoly implementation, since it avoids all the temporaries and extra passes through the data of doing things at the numeric level. For 50 vertices and 1000 candidate inclusion points: nxutils.pnpoly vs pure numpy pnpoly: 50x speedup nxutils.points_inside_poly vs pure numpy pnpoly: 250x speedup #!/usr/bin/env python import time import matplotlib.nxutils as nxutils import matplotlib.numerix as nx import numpy as N def pnpoly(x, y, verts): """Check whether point is in the polygon defined by verts. verts - 2xN array point - (2,) array See http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html """ verts = verts.astype(float) xpi = verts[:,0] ypi = verts[:,1] # shift xpj = xpi[N.arange(xpi.size)-1] ypj = ypi[N.arange(ypi.size)-1] possible_crossings = ((ypi <= y) & (y < ypj)) | ((ypj <= y) & (y < ypi)) xpi = xpi[possible_crossings] ypi = ypi[possible_crossings] xpj = xpj[possible_crossings] ypj = ypj[possible_crossings] crossings = x < (xpj-xpi)*(y - ypi) / (ypj - ypi) + xpi return sum(crossings) % 2 numtrials, numverts, numpoints = 10, 50, 1000 verts = nx.mlab.rand(numverts, 2) points = nx.mlab.rand(numpoints, 2) mask1 = nx.zeros((numpoints, )) mask2 = nx.zeros((numpoints, )) t0 = time.time() for i in range(numtrials): for j in range(numpoints): x, y = points[j] b = pnpoly(x, y, verts) mask1[j] = b told = time.time() - t0 t0 = time.time() for i in range(numtrials): for j in range(numpoints): x, y = points[j] b = nxutils.pnpoly(x, y, verts) mask2[j] = b tnew = time.time() - t0 t0 = time.time() for i in range(numtrials): mask3 = nxutils.points_inside_poly(points, verts) tmany = time.time() - t0 print told, tnew, told/tnew print told, tmany, told/tmany for v0, v1, v2 in zip(mask1, mask2, mask3): assert(v0==v1) assert(v1==v2) JDH

Hi John On Mon, Sep 04, 2006 at 09:31:51PM -0500, John Hunter wrote:

Robert> Raycast with simulation of simplicity to handle degeneracy Robert> is probably your best bet.

Robert> Robert> http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

Robert> John's "add up the angles" approach is not really a good Robert> one. I frequently find it referred to in the literature as Robert> "the worst thing you could possibly do." :-)

OK, I could only stand Robert's taunt for so long. I added the pnpoly routine to matplotlib svn extension code: there is a function to test whether a a single x, y point is inside a polygon (matplotlib.nxutils.pnpoly) and a function which takes a list of points and returns a mask for the points inside (matplotlib.nxutils.points_inside_poly). When profiling against my original "worst thing you could possibly do" implementation, this new version is 15-35 times faster, in addition to being a better algorithm in terms of handling the degenerate cases. It is also considerably faster than the pure numpy pnpoly implementation, since it avoids all the temporaries and extra passes through the data of doing things at the numeric level.

I am glad to see this included. The python implementation is my first choice, but if you are looking for speed specifically, I also wrote weave and ctype versions. Regards Stéfan

...

...

...

...

"John" == John Hunter <jdhunter@ace.bsd.uchicago.edu> writes:

John> need maximum portability, so generally I can't count on the John> presence of weave (no scipy dependency) or ctypes (no python John> 2.4 dependency). But having your python implementation was Hmm, it now occurs to me that ctypes is only in python 2.5. Note to self: never post before the morning coffee has kicked in. JDH

On Wed, Sep 06, 2006 at 02:08:00AM +0200, Stefan van der Walt wrote:

On Tue, Sep 05, 2006 at 09:55:53AM -0500, John Hunter wrote:

...

...

...

...

...

> "John" == John Hunter <jdhunter@ace.bsd.uchicago.edu> writes:

John> need maximum portability, so generally I can't count on the John> presence of weave (no scipy dependency) or ctypes (no python John> 2.4 dependency). But having your python implementation was

Hmm, it now occurs to me that ctypes is only in python 2.5. Note to self: never post before the morning coffee has kicked in.

Either way, the code can be found at

http://mentat.za.net

I guess a Numpy C-API wrapper would be easy to write -- I just never got round to it, since I had the ctypes and weave versions working.

You guys don't waste much time! I see this has already been done in matplotlib SVN. Stéfan