[SciPy-User] Scipy.spatial.Delaunay
Dan Richards
Hi,
I am working on a project that requires Delaunay triangulation of a set of points in 3D space.
I am trying to use ‘Scipy.spatial.Delaunay’ to convert a set of [x,y,z] points into vertices that define the triangulated simplexes/faces, however I am struggling to work out which attributes give me this data?
#So far I simply have a set of points:
Points = [[0,0,0],[10,91,23],[53,4,66],[49,392,49],[39,20,0]]
#I am then using:
x = scipy.spatial.Delaunay(Points)
I am sure for here it is very simple, however I have been struggling to understand which attributes I need access to get the vertices of connecting lines?
· x.vertices?
· x.neighbors?
· x.vertex_to_simplex?
· x.convex_hull...?
If anyone can help or point me in the right direction that would be very much appreciated.
Thanks,
Dan
Pauli Virtanen
[clip]
> I am sure for here it is very simple, however I have been struggling to
> understand which attributes I need access to get the vertices of
> connecting lines?
>
> * x.neighbors?
> * x.vertex_to_simplex?
> * x.convex_hull...?
>
> If anyone can help or point me in the right direction that would be very
> much appreciated.
The edges are recorded in the `vertices` array, which contains indices of
the points making up each triangle. The overall structure is recorded
`neighbors`.
This is maybe easiest to explain in code. The set of edges is:
edges = []
for i in xrange(x.nsimplex):
edges.append((x.vertices[i,0], x.vertices[i,1]))
edges.append((x.vertices[i,1], x.vertices[i,2]))
edges.append((x.vertices[i,2], x.vertices[i,0]))
This however counts each edge multiple times. To get around, that:
edges = []
for i in xrange(x.nsimplex):
if i > x.neighbors[i,2]:
edges.append((x.vertices[i,0], x.vertices[i,1]))
if i > x.neighbors[i,0]:
edges.append((x.vertices[i,1], x.vertices[i,2]))
if i > x.neighbors[i,1]:
edges.append((x.vertices[i,2], x.vertices[i,0]))
This counts each edge only once. Note how the `neighbors` array relates
to `vertices`: its j-th entry gives the neighboring triangle on the
other side of the edge formed that remains after the j-th vertex is
removed from the triangle.
--
Pauli Virtanen
_______________________________________________
SciPy-User mailing list
SciPy...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user
Daniel Lepage
provides a 3D triangulation where each element is a tetrahedron. Each
tetrahedron has four vertices and up to four neighboring tetrahedra.
Pauli's example has fewer connections than you expect because it
assumes a 2D triangulation. In the 3D case, you should have six lines
of code connect each of the four vertices to each other:
edges = []
for i in xrange(x.nsimplex):
edges.append((x.vertices[i,0], x.vertices[i,1]))
edges.append((x.vertices[i,0], x.vertices[i,2]))
edges.append((x.vertices[i,0], x.vertices[i,3]))
edges.append((x.vertices[i,1], x.vertices[i,2]))
edges.append((x.vertices[i,1], x.vertices[i,3]))
edges.append((x.vertices[i,2], x.vertices[i,3]))
--
Dan Lepage
On Wed, Mar 16, 2011 at 3:22 PM, Dan Richards <D.Ric...@mmu.ac.uk> wrote:
> Hi Pauli,
>
> Thanks for your quick reply, I really appreciate the help.
>
> I am still a little confused as to how the points, vertices and neighbors
> relate to one another. Perhaps I can explain how I understand them and you
> can correct me?
>
> When I type x.vertices I get an array that has values for each index:
>>>>x.vertices
> array([[6, 4, 5, 9],
> [8, 6, 4, 5],
> [8, 1, 4, 7],
> [8, 1, 6, 4],
> [3, 6, 4, 9]...])
>
> Do these numbers [w,x,y,z] represent a triangulation whereby the connections
> are as follows?:
>
> w-x
> x-y
> y-z
> w-y
> w-z
>
> Your code did seem to work well, although I added an extra line which I
> assume should have been there?
>
> edges = []
> for i in xrange(x.nsimplex):
> edges.append((x.vertices[i,0], x.vertices[i,1]))
> edges.append((x.vertices[i,1], x.vertices[i,2]))
> edges.append((x.vertices[i,2], x.vertices[i,3])) # New line here
> edges.append((x.vertices[i,3], x.vertices[i,0]))
>
> The confusion on my part is that I expected the vertices to hold three
> indices relating to the points of a triangle so I am confused as to how to
> interpret the four values?
>
> Equally with the neighbors, could you tell me what the four indices define?
>>>>x.neighbors
> array([[-1, -1, 1, 6],
> [ 0, 2, 54, 9],
> [-1, 1, 19, 4],
> [12, 27, 31, 10],
> [ 2, 5, 13, 44]...])
>
> Many Thanks,
> Dan
> "Before acting on this email or opening any attachments you should read the Manchester Metropolitan University email disclaimer available on its website http://www.mmu.ac.uk/emaildisclaimer "
Dan Richards
Thanks for your quick reply, I really appreciate the help.
I am still a little confused as to how the points, vertices and neighbors
relate to one another. Perhaps I can explain how I understand them and you
can correct me?
When I type x.vertices I get an array that has values for each index:
>>>x.vertices
array([[6, 4, 5, 9],
[8, 6, 4, 5],
[8, 1, 4, 7],
[8, 1, 6, 4],
[3, 6, 4, 9]...])
Do these numbers [w,x,y,z] represent a triangulation whereby the connections
are as follows?:
w-x
x-y
y-z
w-y
w-z
Your code did seem to work well, although I added an extra line which I
assume should have been there?
edges = []
for i in xrange(x.nsimplex):
edges.append((x.vertices[i,0], x.vertices[i,1]))
edges.append((x.vertices[i,1], x.vertices[i,2]))
edges.append((x.vertices[i,2], x.vertices[i,3])) # New line here
edges.append((x.vertices[i,3], x.vertices[i,0]))
The confusion on my part is that I expected the vertices to hold three
indices relating to the points of a triangle so I am confused as to how to
interpret the four values?
Equally with the neighbors, could you tell me what the four indices define?
>>>x.neighbors
array([[-1, -1, 1, 6],
[ 0, 2, 54, 9],
[-1, 1, 19, 4],
[12, 27, 31, 10],
[ 2, 5, 13, 44]...])
Many Thanks,
Dan
The edges are recorded in the `vertices` array, which contains indices of
--
Pauli Virtanen
"Before acting on this email or opening any attachments you should read the Manchester Metropolitan University email disclaimer available on its website http://www.mmu.ac.uk/emaildisclaimer "
Pauli Virtanen
[clip]
> Pauli's example has fewer connections than you expect because it assumes
> a 2D triangulation. In the 3D case, you should have six lines of code
> connect each of the four vertices to each other:
Yes, sorry, I didn't read the original mail carefully enough. For 3D,
each tetrahedron has 6 edges, and they can be found as Daniel shows:
> edges = []
> for i in xrange(x.nsimplex):
> edges.append((x.vertices[i,0], x.vertices[i,1]))
> edges.append((x.vertices[i,0], x.vertices[i,2]))
> edges.append((x.vertices[i,0], x.vertices[i,3]))
> edges.append((x.vertices[i,1], x.vertices[i,2]))
> edges.append((x.vertices[i,1], x.vertices[i,3]))
> edges.append((x.vertices[i,2], x.vertices[i,3]))
If you want each edge to appear only once, the trick of using
x.neighbours will not work in 3D, since ndim - 1 != 1. So a brute
force algorithm needs to be used, for example detecting duplicates
while processing, or by removing them afterwards. (This information
is not available from Qhull.)
If you want to write something yourself using Cython, this may give
some ideas on how to proceed:
https://github.com/pv/scipy-work/blob/enh%2Finterpnd-smooth/scipy/spatial/qhull.pyx#L969
Note that it'll need some adaptation.
Robert Kern
> Hi Pauli,
>
> Thanks for your quick reply, I really appreciate the help.
>
> I am still a little confused as to how the points, vertices and neighbors
> relate to one another. Perhaps I can explain how I understand them and you
> can correct me?
>
> When I type x.vertices I get an array that has values for each index:
>>>>x.vertices
> array([[6, 4, 5, 9],
> [8, 6, 4, 5],
> [8, 1, 4, 7],
> [8, 1, 6, 4],
> [3, 6, 4, 9]...])
>
> Do these numbers [w,x,y,z] represent a triangulation whereby the connections
> are as follows?:
>
> w-x
> x-y
> y-z
> w-y
> w-z
And x-z. It's a tetrahedralization, technically. But you probably
don't want to deal with edges. Rather, you usually want to deal with
faces.
w-x-y
x-z-y
w-z-x
w-y-z
> Your code did seem to work well, although I added an extra line which I
> assume should have been there?
>
> edges = []
> for i in xrange(x.nsimplex):
> edges.append((x.vertices[i,0], x.vertices[i,1]))
> edges.append((x.vertices[i,1], x.vertices[i,2]))
> edges.append((x.vertices[i,2], x.vertices[i,3])) # New line here
> edges.append((x.vertices[i,3], x.vertices[i,0]))
>
> The confusion on my part is that I expected the vertices to hold three
> indices relating to the points of a triangle so I am confused as to how to
> interpret the four values?
He was showing you the 2D case for triangles. For the 3D case of
tetrahedra, you probably want faces rather than edges.
faces = []
v = x.vertices
for i in xrange(x.nsimplex):
faces.extend([
(v[i,0], v[i,1], v[i,2]),
(v[i,1], v[i,3], v[i,2]),
(v[i,0], v[i,3], v[i,1]),
(v[i,0], v[i,2], v[i,3]),
])
> Equally with the neighbors, could you tell me what the four indices define?
>>>>x.neighbors
> array([[-1, -1, 1, 6],
> [ 0, 2, 54, 9],
> [-1, 1, 19, 4],
> [12, 27, 31, 10],
> [ 2, 5, 13, 44]...])
Each is the index into x.vertices of the simplex that is adjacent to
the face opposite the given point. This index is -1 when that face is
on the outer boundary. The first simplex has these neighbors:
[-1, -1, 1, 6]
This means that the face opposite point w (x-y-z) is on the outer
boundary. The face opposite point y (w-x-z) adjoins the simplex given
by x.vertices[1], and so on.
--
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
Dan Richards
I will look into writing my own brute force algorithm to remove duplicate
edges/faces as your example.
>https://github.com/pv/scipy-work/blob/enh%2Finterpnd-smooth/scipy/spatial/q
hull.pyx#L969
Thanks again, much appreciated.
Dan
"Before acting on this email or opening any attachments you should read the Manchester Metropolitan University email disclaimer available on its website http://www.mmu.ac.uk/emaildisclaimer "
Daniel Lepage
make sure to sort the indices so that each edge always has the same
representation):
edges = set()
def add_edge(v1, v2):
edges.add((min(v1, v2), max(v1, v2))
for i in xrange(x.nsimplex):
add_edge(x.vertices[i,0], x.vertices[i,1])
add_edge(x.vertices[i,0], x.vertices[i,2])
add_edge(x.vertices[i,0], x.vertices[i,3])
add_edge(x.vertices[i,1], x.vertices[i,2])
add_edge(x.vertices[i,1], x.vertices[i,3])
add_edge(x.vertices[i,2], x.vertices[i,3])
You can then iterate over edges directly, or call `list(edges)` if you
need them in an ordered list.
--
Dan Lepage
On Thu, Mar 17, 2011 at 5:20 AM, Dan Richards <D.Ric...@mmu.ac.uk> wrote:
> Thanks Pauli, Dan and Rob, this is excellent!
>
> I will look into writing my own brute force algorithm to remove duplicate
> edges/faces as your example.
>
>>https://github.com/pv/scipy-work/blob/enh%2Finterpnd-smooth/scipy/spatial/q
> hull.pyx#L969
>
> Thanks again, much appreciated.
> Dan
>
> "Before acting on this email or opening any attachments you should read the Manchester Metropolitan University email disclaimer available on its website http://www.mmu.ac.uk/emaildisclaimer "