3
"""isosurface.py -- plot isosurfaces
5
Heavily based on iso.c in the GTS examples directory
7
Copyright (C) 2009 Richard Everson
10
Richard Everson <R.M.Everson@exeter.ac.uk>
11
School of Engineering, Computing and Mathematics,
17
This library is free software; you can redistribute it and/or
18
modify it under the terms of the GNU Library General Public
19
License as published by the Free Software Foundation; either
20
version 2 of the License, or (at your option) any later version.
22
This library is distributed in the hope that it will be useful,
23
but WITHOUT ANY WARRANTY; without even the implied warranty of
24
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
25
Library General Public License for more details.
27
You should have received a copy of the GNU Library General Public
28
License along with this library; if not, write to the
29
Free Software Foundation, Inc., 59 Temple Place - Suite 330,
30
Boston, MA 02111-1307, USA.
34
from optparse import OptionParser
36
from enthought.mayavi import mlab
39
# Visualise a bug in GTS by running this with args --function=ellipsoid
40
# -c80.0 --method=dual
41
# In addition to the expected ellipsoid, there is an open surface in the
42
# z=-10 plane which appears to be the "shadow" of the ellipsoid intersected
43
# with the data grid. Running the GTS/examples/iso -d 50 50 50 80 yields
44
# similar, confirming that this is a GTS bug.
50
x, y, z = numpy.ogrid[-10:10:Nj, -10:10:Nj, -10:10:Nj]
51
scalars = x*x + 2.0*y*y + z*z/2.0
55
# The Clebsch diagonal surface: a smooth cubic surface admitting the
56
# symmetry of the tetrahedron (courtesy Johannes Beigel via GTS).
59
x, y, z = numpy.ogrid[-10:10:Nj, -10:10:Nj, -10:10:Nj]
65
c2 = p*p*p + q*q*q + r*r*r - s*s*s
70
functions = { "ellipsoid" : ellipsoid,
74
parser = OptionParser()
76
parser.add_option("-n", "--N", type="int",
78
help="Size of data cube")
80
parser.add_option("-c", "--isovalue", type="string",
81
dest="c", default="80.0",
82
help="""Comma separated list of values defining isosurfaces.
83
Try: --isovalue=25.0,120.0 --function=ellipsoid
86
parser.add_option("--function", type="string", dest="function",
88
help="Function defining surface: ['ellipsoid'|'clebsch']")
90
parser.add_option("--method", type="string",
91
dest="method", default="cube",
92
help="Iso-surface generation method: " \
93
"['cube'|'tetra'|'dual'|'bounded']")
95
(opt, args) = parser.parse_args()
97
extents = numpy.asarray([-10.0, 10.0, -10.0, 10.0, -10.0, 10.0])
98
func = functions[opt.function]
102
for c in [string.atof(x) for x in opt.c.split(',')]:
103
S.add(gts.isosurface(data, c, extents=extents, method=opt.method))
104
print S.stats()['n_faces'], 'facets'
105
x,y,z,t = gts.get_coords_and_face_indices(S,True)
107
mlab.triangular_mesh(x,y,z,t,color=(0.25,0.25,0.75))