~ubuntu-branches/ubuntu/wily/pygts/wily

« back to all changes in this revision

Viewing changes to .pc/fix_examples.patch/examples/isosurface.py

  • Committer: Package Import Robot
  • Author(s): Anton Gladky, Václav Šmilauer
  • Date: 2013-01-26 21:22:37 UTC
  • Revision ID: package-import@ubuntu.com-20130126212237-5l783abxeg1k9iuh
Tags: 0.3.1-1
[ Václav Šmilauer ]
Initial debian package. (Closes: #691799)

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
#! /usr/bin/env python
 
2
 
 
3
"""isosurface.py -- plot isosurfaces
 
4
 
 
5
   Heavily based on iso.c in the GTS examples directory
 
6
 
 
7
   Copyright (C) 2009 Richard Everson
 
8
   All rights reserved.
 
9
 
 
10
   Richard Everson <R.M.Everson@exeter.ac.uk>
 
11
   School of Engineering, Computing and Mathematics,
 
12
   University of Exeter
 
13
   Exeter, EX4 4 QF.  UK. 
 
14
 
 
15
NOTICE
 
16
 
 
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.
 
21
 
 
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.
 
26
 
 
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.
 
31
"""
 
32
 
 
33
import sys, string
 
34
from optparse import OptionParser
 
35
import numpy
 
36
from enthought.mayavi import mlab
 
37
import gts
 
38
 
 
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.
 
45
 
 
46
 
 
47
 
 
48
def ellipsoid(N=50):
 
49
    Nj = N*(0+1j)
 
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
 
52
    return scalars
 
53
 
 
54
def clebsch(N=50):
 
55
    # The Clebsch diagonal surface: a smooth cubic surface admitting the
 
56
    # symmetry of the tetrahedron (courtesy Johannes Beigel via GTS).
 
57
    Nj = N*(0+1j)
 
58
    w2 = numpy.sqrt(2.0)
 
59
    x, y, z = numpy.ogrid[-10:10:Nj, -10:10:Nj, -10:10:Nj]
 
60
    p = 1 - z - w2*x
 
61
    q = 1 - z + w2*x
 
62
    r = 1 + z + w2*y
 
63
    s = 1 + z - w2*y
 
64
    c1 = p + q + r - s
 
65
    c2 = p*p*p + q*q*q + r*r*r - s*s*s
 
66
    return c2 - c1*c1*c1
 
67
 
 
68
 
 
69
 
 
70
functions = { "ellipsoid" : ellipsoid,
 
71
              "clebsch" : clebsch
 
72
            }
 
73
 
 
74
parser = OptionParser()
 
75
 
 
76
parser.add_option("-n", "--N", type="int",
 
77
                  dest="N", default=50,
 
78
                  help="Size of data cube")
 
79
 
 
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
 
84
                  """)
 
85
 
 
86
parser.add_option("--function", type="string", dest="function",
 
87
                  default="ellipsoid", 
 
88
                  help="Function defining surface: ['ellipsoid'|'clebsch']")
 
89
 
 
90
parser.add_option("--method", type="string",
 
91
                  dest="method", default="cube",
 
92
                  help="Iso-surface generation method: " \
 
93
                      "['cube'|'tetra'|'dual'|'bounded']")
 
94
 
 
95
(opt, args) = parser.parse_args()
 
96
 
 
97
extents = numpy.asarray([-10.0, 10.0, -10.0, 10.0, -10.0, 10.0])
 
98
func = functions[opt.function]
 
99
data = func(opt.N)
 
100
            
 
101
S = gts.Surface()
 
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)
 
106
mlab.clf()
 
107
mlab.triangular_mesh(x,y,z,t,color=(0.25,0.25,0.75))
 
108
 
 
109
mlab.show()