1
## Automatically adapted for scipy Oct 31, 2005 by
3
# $Id: pl3d.py 2183 2006-08-29 10:30:44Z oliphant $
4
# Copyright (c) 1996, 1997, The Regents of the University of California.
5
# All rights reserved. See Legal.htm for full text and disclaimer.
8
from shapetest import *
11
from gistfuncs import *
15
# Viewing transforms and other aids for 3D plotting.
17
# $Id: pl3d.py 2183 2006-08-29 10:30:44Z oliphant $
19
# Copyright (c) 1997. The Regents of the University of California.
20
# All rights reserved.
23
General overview of module pl3d:
25
(1) Viewing transform machinery. Arguably the simplest model
26
is the CAD/CAM notion that the object you see is oriented
27
as you see it in the current picture. You can then move
28
it left, right, up, down, or toward or away from you,
29
or you can rotate it about any of the three axes (horizontal,
30
vertical, or out of the screen). The xyz coordinates of the
31
object remains unchanged throughout all of this, but this
32
object coordinate system changes relative to the fixed
33
xyz of the viewer, in which x is always to the right, y is
34
up, and z is directed out of the screen. Initially, the
35
two coordinate systems coincide.
37
rot3 (xangle,yangle,zangle)
38
Rotate the object about viewer's x-axis by xangle, then
39
about viewer's y-axis by yangle, then about viewer's
41
mov3 (xchange,ychange,zchange)
42
Move the object by the specified amounts.
45
The "camera" is located at (0,0,zcamera) in the viewer's
46
coordinate system, looking in the minus-z direction.
47
Initially, zcamera is very large, and the magnification
48
factor is correspondingly large, giving an isometric view.
49
Decreasing zcamera makes the perspective more extreme.
50
If parts of the object are behind the camera, strange things
55
Undo the last N (default 1) viewpoint commands (rot3, mov3,
56
or setz3). Up to 100 viewpoint changes are remembered.
60
The current viewpoint transformation can be saved and later
64
Toggle the gnomon (a simple display showing the orientation
65
of the xyz axes of the object).
68
# ------------------------------------------------------------------------
71
def set_draw3_ ( n ) :
74
set_draw3_ ( 0 | 1 ) is used to set the global draw3_,
75
which controls whether the function draw3 actually shows a drawing.
83
# ZCM 2/21/97 change reflects the fact that I hadn't realized
84
# that car and cdr, as functions, return the item replaced.
87
oldx = _draw3_list [0]
89
undo3_set_ (setrot3_, oldx)
91
def rot3 (xa = 0., ya = 0., za = 0.) :
95
rotate the current 3D plot by XA about viewer's x-axis,
96
YA about viewer's y-axis, and ZA about viewer's z-axis.
97
SEE ALSO: orient3, mov3, aim3, setz3, undo3, save3, restore3, light3
100
x = array ([1.,0.,0.], Float)
101
y = array ([0.,1.,0.], Float)
102
z = array ([0.,0.,1.], Float)
103
[x, y] = rot3_ (za, x, y)
104
[z, x] = rot3_ (ya, z, x)
105
[y, z] = rot3_ (xa, y, z)
106
# n. b. matrixMultiply has the unfortunate effect of destroying
107
# the matrix that calls it.
108
gr3 = array (getrot3_ (), copy = 1)
109
setrot3_ (transpose (dot (transpose (gr3), array ( [x, y, z]))))
111
def rot3_ (a, x, y) :
114
return [multiply (ca, x) + multiply (sa, y), multiply (-sa, x) + multiply (ca, y)]
116
def mov3 ( xa = 0., ya = 0., za = 0. ) :
119
mov3 ( [xa [, ya [, za]]])
120
move the current 3D plot by XA along the viewer's x axis,
121
YA along the viewer's y axis, and ZA along the viewer's z axis.
122
SEE ALSO: rot3, orient3, setz3, undo3, save3, restore3, light3
125
gr = dot (transpose (gr), transpose (xa))
126
setorg3_ ( getorg3_ () - gr)
128
def aim3 ( xa = 0., ya = 0., za = 0. ) :
131
aim3 ( [xa [, ya [, za]]])
132
move the current 3D plot to put the point (XA, YA, ZA) in object
133
coordinates at the point (0, 0, 0) -- the aim point -- in the
134
viewer's coordinates. If any of the XA, YA, or ZA is nil, it defaults
135
SEE ALSO: mov3, rot3, orient3, setz3, undo3, save3, restore3, light3
142
def setz3 ( zc = None ) :
146
Set the camera position to z = ZC (x = y = 0) in the viewer's coordinate
147
system. If zc is None, set the camera to infinity (default).
148
SEE ALSO: rot3, orient3, undo3, save3, restore3, light3
151
if not is_scalar (zc) :
152
raise _ZcError, "camera position must be scalar."
156
def orient3 ( ** kw ) :
159
orient3 ( [phi = val1, theta = val2] )
160
Set the orientation of the object to (PHI, THETA). Orientations
161
are a subset of the possible rotation matrices in which the z axis
162
of the object appears vertical on the screen (that is, the object
163
z axis projects onto the viewer y axis). The THETA angle is the
164
angle from the viewer y axis to the object z axis, positive if
165
the object z axis is tilted towards you (toward viewer +z). PHI is
166
zero when the object x axis coincides with the viewer x axis. If
167
neither PHI nor THETA is specified, PHI defaults to - pi / 4 and
168
THETA defaults to pi / 6. If only PHI is specified, THETA remains
169
unchanged, unless the current THETA is near pi / 2, in which case
170
THETA returns to pi / 6, or unless the current orientation does
171
not have a vertical z axis, in which case THETA returns to its
173
Unlike rot3, orient3 is not a cumulative operation.
175
# Notes with regard to global variables: (ZCM 2/21/97)
176
# _orient3_phi, _orient3_theta, the default orientation angles,
177
# are known and referred to only in this routine. I have started
178
# them with an underscore, too, to make them inaccessible
179
# from outside this module.
180
# phi and theta need not be global here since they are recalculated
181
# each time this routine is called.
183
global _orient3_phi, _orient3_theta
185
dummy = _orient3_theta
187
_orient3_theta = pi / 6.
192
_orient3_phi = - pi / 4.
194
if kw.has_key ("phi") and kw ["phi"] == None :
195
kw ["phi"] = _orient3_phi
196
if kw.has_key ("theta") and kw ["theta"] == None :
197
kw ["theta"] = _orient3_theta
198
if not kw.has_key ("phi") and not kw.has_key ("theta") :
200
theta = _orient3_theta
201
elif not kw.has_key ("phi") or not kw.has_key ("theta") :
202
gr3 = array (getrot3_ (), copy = 1)
203
z = dot (transpose (gr3), array ( [0., 0., 1.]))
204
if abs (z [0]) > 1.e-6 :
205
# object z-axis not aligned with viewer y-axis
206
if not kw.has_key ("theta") :
207
theta = _orient3_theta
212
elif not kw.has_key ("theta") :
214
if (abs (z [1]) < 1.e-6) :
215
theta = _orient3_theta
217
theta = arctan2 (z [2], z [1])
220
y = array ( [0., z [2], -z [1]])
221
x = dot (transpose (gr3), array ( [1., 0., 0.]))
222
phi = arctan2 (sum (y * x,axis=0), x [0])
227
x = array ( [1., 0., 0.], Float)
228
y = array ( [0., 1., 0.], Float)
229
z = array ( [0., 0., 1.], Float)
230
[y, z] = rot3_ (theta, y, z)
231
[z, x] = rot3_ (phi, z, x)
232
setrot3_ (array ( [x, -z, y], Float))
240
Save the current 3D viewing transformation and lighting.
241
Actually, this doesn't save anything; it returns a copy
242
of the current 3D viewing transformation and lighting, so
243
that the user can put it aside somewhere.
244
SEE ALSO: restore3, rot3, mov3, aim3, light3
247
return _draw3_list [0:_draw3_n]
249
def restore3 ( view = None ) :
253
Restore a previously saved 3D viewing transformation and lighting.
254
If view is missing, rotate object to viewer's coordinate system.
255
SEE ALSO: restore3, rot3, mov3, aim3, light3
258
global _draw3_list, _draw3_view, _light3_list, _draw3_n
261
view = view [0:len (view)] # Copies view
263
view = _draw3_view + _light3_list
264
old = _draw3_list [0:_draw3_n]
265
_draw3_list = view [0:_draw3_n] + _draw3_list [_draw3_n:]
266
undo3_set_ (restore3, old)
268
_AmbientError = "AmbientError"
269
_DiffuseError = "DiffuseError"
270
_LightingError = "LightingError"
272
def light3 ( * kw, ** kwds ) :
275
light3 (ambient=a_level,
280
Sets lighting properties for 3D shading effects.
281
A surface will be shaded according to its to its orientation
282
relative to the viewing direction.
284
The ambient level A_LEVEL is a light level (arbitrary units)
285
that is added to every surface independent of its orientation.
287
The diffuse level D_LEVEL is a light level which is proportional
288
to cos(theta), where theta is the angle between the surface
289
normal and the viewing direction, so that surfaces directly
290
facing the viewer are bright, while surfaces viewed edge on are
291
unlit (and surfaces facing away, if drawn, are shaded as if they
294
The specular level S_LEVEL is a light level proportional to a high
295
power spower=N of 1+cos(alpha), where alpha is the angle between
296
the specular reflection angle and the viewing direction. The light
297
source for the calculation of alpha lies in the direction XYZ (a
298
3 element vector) in the viewer's coordinate system at infinite
299
distance. You can have ns light sources by making S_LEVEL, N, and
300
XYZ (or any combination) be vectors of length ns (3-by-ns in the
301
case of XYZ). (See source code for specular_hook function
302
definition if powers of 1+cos(alpha) aren't good enough for you.)
304
With no arguments, return to the default lighting.
307
light3 ( diffuse=.1, specular=1., sdir=[0,0,-1])
308
(dramatic "tail lighting" effect)
309
light3 ( diffuse=.5, specular=1., sdir=[1,.5,1])
310
(classic "over your right shoulder" lighting)
311
light3 ( ambient=.1,diffuse=.1,specular=1.,
312
sdir=[[0,0,-1],[1,.5,1]],spower=[4,2])
313
(two light sources combining previous effects)
314
SEE ALSO: rot3, save3, restore3
317
global _draw3_list, _draw3_nv
318
if len (kw) > 0 : kwds = kw [0]
319
old = _draw3_list [_draw3_nv:] [0:5]
321
if kwds.has_key ("ambient") and kwds ["ambient"] != None :
322
ambient = kwds ["ambient"]
323
if not is_scalar (ambient) :
324
raise _AmbientError, "ambient light level must be scalar."
326
_draw3_list [_draw3_nv] = ambient
327
if kwds.has_key ("diffuse") and kwds ["diffuse"] != None :
328
diffuse = kwds ["diffuse"]
329
if not is_scalar (diffuse) :
330
raise _DiffuseError, "diffuse light level must be scalar."
332
_draw3_list [_draw3_nv + 1 ] = diffuse
334
if kwds.has_key ("specular") and kwds ["specular"] != None :
335
specular = kwds ["specular"]
338
specular = _draw3_list [_draw3_nv + 2]
339
if kwds.has_key ("spower") and kwds ["spower"] != None :
340
spower = kwds ["spower"]
343
spower = _draw3_list [_draw3_nv + 3]
344
if kwds.has_key ("sdir") and kwds ["sdir"] != None :
347
if dims == 0 or len (dims) == 2 and dims [1] != 3 :
348
raise _LightingError, \
349
"lighting direction must be 3 vector or ns-by-3 array."
352
sdir = _draw3_list [_draw3_nv + 4]
354
if flags & 4 : _draw3_list [_draw3_nv + 2] = specular
355
if flags & 8 : _draw3_list [_draw3_nv + 3] = spower
356
if flags & 16 : _draw3_list [_draw3_nv + 4] = sdir
358
_draw3_list [_draw3_nv: _draw3_nv + 5] = _light3_list [0:5]
359
undo3_set_ (light3_, old)
362
global _draw3_list, _draw3_nv
364
_draw3_list [_draw3_nv:_draw3_nv + 5] = arg [0:5]
366
def get3_light (xyz, * nxyz) :
369
get3_light(xyz, nxyz)
372
return 3D lighting for polygons with vertices XYZ. If NXYZ is
373
specified, XYZ should be sum(nxyz,axis=0)-by-3, with NXYZ being the
374
list of numbers of vertices for each polygon (as for the plfp
375
function). If NXYZ is not specified, XYZ should be a quadrilateral
376
mesh, ni-by-nj-by-3 (as for the plf function). In the first case,
377
the return value is len (NXYZ) long; in the second case, the
378
return value is (ni-1)-by-(nj-1).
380
The parameters of the lighting calculation are set by the
383
SEE ALSO: light3, set3_object, get3_normal, get3_centroid
386
global _draw3_list, _draw3_nv
387
list = _draw3_list [_draw3_nv:]
395
normal = get3_normal (xyz)
397
normal = get3_normal (xyz, nxyz [0])
401
view = array ( [0., 0., 1.], Float)
402
elif len (nxyz) == 0 :
403
view = array ( [0., 0., zc], Float) - get3_centroid (xyz)
405
view = array ( [0., 0., zc], Float) - get3_centroid (xyz, nxyz [0])
407
sqrt ( sum (view * view,axis=0))
408
if m1 == 0. : m1 = 1.
411
nv = normal [0, ...] * view [0] + normal [1, ...] * view [1] + \
412
normal [2, ...] * view [2]
413
light = ambient + diffuse * abs (nv)
415
sv = transpose (transpose (sdir) / sqrt (sum (transpose (sdir*sdir),axis=0)))
417
if len (shape (sdir)) == 1 :
418
sn = sum(array([sdir[0]*normal[0],sdir[1]*normal[1],
419
sdir[2]*normal[2]]),axis=0)
420
####### I left out the specular_hook stuff.
421
m1 = maximum (sn * nv -0.5 * sv + 0.5, 1.e-30)
423
light = light + (specular * m1)
424
elif len (shape (sdir)) >= 2 :
425
# multiple light sources
426
nsrc = len (shape (sdir))
427
for i in range (nsrc) :
428
sn = sum(array([sdir[i,0]*normal[0],sdir[i,1]*normal[1],
429
sdir[i,2]*normal[2]]),axis=0)
430
m1 = maximum (sn * nv -0.5 * sv [i] + 0.5, 1.e-30) ** spower [i]
431
light = light + specular * m1
434
def get3_normal (xyz, *nxyz) :
437
get3_normal(xyz, nxyz)
440
return 3D normals for polygons with vertices XYZ. If NXYZ is
441
specified, XYZ should be sum(nxyz,axis=0)-by-3, with NXYZ being the
442
list of numbers of vertices for each polygon (as for the plfp
443
function). If NXYZ is not specified, XYZ should be a quadrilateral
444
mesh, ni-by-nj-by-3 (as for the plf function). In the first case,
445
the return value is len(NXYZ)-by-3; in the second case, the
446
return value is (ni-1)-by-(nj-1)-by-3.
448
The normals are constructed from the cross product of the lines
449
joining the midpoints of two edges which as nearly quarter the
450
polygon as possible (the medians for a quadrilateral). No check
451
is made that these not be parallel; the returned "normal" is
452
[0,0,0] in that case. Also, if the polygon vertices are not
453
coplanar, the "normal" has no precisely definable meaning.
455
SEE ALSO: get3_centroid, get3_light
459
# if no polygon list is given, assume xyz is 2D mesh
460
# form normal as cross product of medians
461
m1 = dif_ (zcen_ (xyz, 1), 2)
462
m2 = zcen_ (dif_ (xyz, 1), 2)
464
# with polygon list, more elaborate calculation required
465
# (1) frst subscripts the first vertex of each polygon
466
frst = cumsum (nxyz [0],axis=0) - nxyz [0]
468
# form normal by getting two approximate diameters
469
# (reduces to above medians for quads)
470
# (2) compute midpoints of first three sides
471
n2 = (nxyz [0] + 1) / 2
472
c0 = (take(xyz, frst, 0) + take(xyz, frst + 1, 0)) / 2.
474
c1 = (take(xyz, i, 0) + take(xyz, i + 1, 0)) / 2.
476
c2 = (take(xyz, frst + i, 0) + take(xyz, frst + (i + 1) % nxyz [0], 0)) / 2.
477
i = minimum (i + n2, nxyz [0]) - 1
478
c3 = (take(xyz, frst + i, 0) + take(xyz, frst + (i + 1) % nxyz [0], 0)) / 2.
482
# poly normal is cross product of two medians (or diameters)
483
# normal = m1; I had to reverse the sign.
484
if len (shape (xyz)) == 3 :
485
n1 = m1 [2, :] * m2 [1, :] - \
486
m1 [1, :] * m2 [2, :]
487
n2 = m1 [0, :] * m2 [2, :] - \
488
m1 [2, :] * m2 [0, :]
489
n3 = m1 [1, :] * m2 [0, :] - \
490
m1 [0, :] * m2 [1, :]
492
n1 = m1 [:, 2] * m2 [:, 1] - \
493
m1 [:, 1] * m2 [:, 2]
494
n2 = m1 [:, 0] * m2 [:, 2] - \
495
m1 [:, 2] * m2 [:, 0]
496
n3 = m1 [:, 1] * m2 [:, 0] - \
497
m1 [:, 0] * m2 [:, 1]
498
m1 = sqrt (n1 ** 2 + n2 **2 + n3 **2)
499
m1 = m1 + equal (m1, 0.0)
500
normal = array([n1 / m1, n2 / m1, n3 / m1])
504
def get3_centroid (xyz, * nxyz) :
507
get3_centroid(xyz, *nxyz)
508
or get3_centroid(xyz)
510
return 3D centroids for polygons with vertices XYZ. If NXYZ is
511
specified, XYZ should be sum(nxyz,axis=0)-by-3, with NXYZ being the
512
list of numbers of vertices for each polygon (as for the plfp
513
function). If NXYZ is not specified, XYZ should be a quadrilateral
514
mesh, ni-by-nj-by-3 (as for the plf function). In the first case,
515
the return value is len(NXYZ) in length; in the second case, the
516
return value is (ni-1)-by-(nj-1)-by-3.
518
The centroids are constructed as the mean value of all vertices
521
SEE ALSO: get3_normal, get3_light
525
# if no polygon list is given, assume xyz is 2D mesh
526
centroid = zcen_ (zcen_ (xyz, 1), 0)
528
# with polygon list, more elaborate calculation required
529
last = cumsum (nxyz [0],axis=0)
530
list = histogram (1 + last) [0:-1]
531
list = cumsum (list,axis=0)
534
centroid = zeros ( (k, 3))
535
centroid [0:k, 0] = histogram (list, xyz [0:l,0])
536
centroid [0:k, 1] = histogram (list, xyz [0:l,1])
537
centroid [0:k, 2] = histogram (list, xyz [0:l,2])
538
fnxyz = array (nxyz [0], Float )
539
centroid = centroid / fnxyz
542
_Get3Error = "Get3Error"
544
def get3_xy (xyz, *flg) :
550
Given anything-by-3 coordinates XYZ, return X and Y in viewer's
551
coordinate system (set by rot3, mov3, orient3, etc.). If the
552
second argument is present and non-zero, also return Z (for use
553
in sort3d or get3_light, for example). If the camera position
554
has been set to a finite distance with setz3, the returned
555
coordinates will be tangents of angles for a perspective
556
drawing (and Z will be scaled by 1/zc).
557
Unlike the Yorick version, this function returns a 3-by-anything
558
array of coordinates.
559
Actually, what it returns is a 3-by-anything python array, whose
560
0th element is the x array, whose 1th element is the y array, and
561
whose 2th element is the z array if asked for.
562
I believe that x, y, and z can be either 1d or 2d, so this
563
routine is written in two cases.
567
# rotate and translate to viewer's coordinate system
570
# 2d mesh case is much more complex than in Yorick
577
tmpxyz = array ( [xx, yy, zz])
578
gr3 = array (getrot3_ (), copy = 1)
579
tmpxyz = dot (transpose (gr3),
580
tmpxyz - array ( [ [go3_ [0]], [go3_ [1]], [go3_ [2]]]))
581
## xx = transpose (reshape (ravel (tmpxyz [0]), (k,l)))
582
## yy = transpose (reshape (ravel (tmpxyz [1]), (k,l)))
583
## zz = transpose (reshape (ravel (tmpxyz [2]), (k,l)))
584
xx = (reshape (ravel (tmpxyz [0]), (k,l)))
585
yy = (reshape (ravel (tmpxyz [1]), (k,l)))
586
zz = (reshape (ravel (tmpxyz [2]), (k,l)))
587
tmpxyz = array ( [xx, yy, zz])
590
lm = array (getrot3_ (), copy = 1)
591
rm = (xyz - array ( [ go3_ [0], go3_ [1], go3_ [2]]))
592
tmpxyz = dot (rm, lm)
594
raise _Get3Error, "xyz has a bad shape: " + `shp`
596
# do optional perspective projection
601
zc = maximum (zc - z, 1.e-35) # protect behind camera, avoid zero divide
602
tmpxyz [:, 0] = tmpxyz [:, 0] / zc
603
tmpxyz [:, 1] = tmpxyz [:, 1] / zc
604
if len (flg) != 0 and flg [0] != 0 :
605
tmpxyz [:, 2] = tmpxyz [:, 2] / zc
606
elif len (shp) == 3 :
608
zc = maximum (zc - z, 1.e-35) # protect behind camera, avoid zero divide
609
tmpxyz [:,:, 0] = tmpxyz [:,:, 0] / zc
610
tmpxyz [:,:, 1] = tmpxyz [:,:, 1] / zc
611
if len (flg) != 0 and flg [0] != 0 :
612
tmpxyz [:,:, 2] = tmpxyz [:,:, 2] / zc
615
_UndoError = "UndoError"
625
Undo the effects of the last N (default 1) rot3, orient3, mov3, aim3,
626
setz3, or light3 commands.
629
global _in_undo3, _undo3_list
631
if n < 0 or n > len (_undo3_list) :
632
raise _UndoError, "not that many items in undo list"
633
_in_undo3 = 1 # flag to skip undo3_set_
634
# perhaps should save discarded items in a redo list?
635
use_list = undo3_list [-n:]
636
undo3_list = undo3_list [:-n]
647
def set3_object (fnc, arg) :
650
set3_object (drawing_function, [arg1,arg2,...])
652
set up to trigger a call to draw3, adding a call to the
653
3D display list of the form:
655
DRAWING_FUNCTION ( [ARG1, ARG2, ...]))
657
When draw3 calls DRAWING_FUNCTION, the external variable draw3_
658
will be non-zero, so DRAWING_FUNCTION can be written like this:
660
def drawing_function(arg) :
666
...<calls to get3_xy, sort3d, get3_light, etc.>...
667
...<calls to graphics functions plfp, plf, etc.>...
671
...<do orientation and lighting independent calcs>...
672
set3_object (drawing_function, [arg1,arg2,...])
674
SEE ALSO: get3_xy, get3_light, sort3d
678
_draw3_list = _draw3_list + [fnc, arg]
682
# ZCM 2/21/97 change reflects the fact that I hadn't realized
683
# that car and cdr, as functions, return the item replaced.
685
oldx = _draw3_list [1]
687
undo3_set_ ( setorg3_, oldx)
690
# ZCM 2/21/97 change reflects the fact that I hadn't realized
691
# that car and cdr, as functions, return the item replaced.
693
oldx = _draw3_list [2]
695
undo3_set_ ( setzc3_, oldx)
698
return _draw3_list [0]
701
return _draw3_list [1]
704
return _draw3_list [2]
706
def undo3_set_ (fnc, arg) :
707
global _undo3_list, _in_undo3, _undo3_limit
708
# arg = copy.deepcopy (arg)
710
if len (_undo3_list) >= 2 * _undo3_limit :
711
_undo3_list = _undo3_list [0:2 * _undo3_limit - 2]
712
_undo3_list = [fnc, arg] + _undo3_list
715
_in_undo3 = 0 # ??????????????
722
def clear_idler ( ) :
723
_idler = do_nothing ( )
725
def set_idler ( fnc ) :
733
def _draw3_idler ( ) :
734
# I have added orientation and limits to this because they may not
735
# have been set by a previous command. If the user doesn't like this,
736
# he/she can write his/her own idler. (ZCM 7/1/97)
737
global _default_gnomon
739
if current_window () == -1 :
742
window3 (current_window ())
743
gnomon (_default_gnomon)
748
limits (lims [0], lims [1], lims [2], lims [3])
750
def set_default_idler ( ) :
751
set_idler (_draw3_idler)
753
set_default_idler ( )
755
_draw3_changes = None
757
def set_multiple_components ( n = 0 ) :
758
global _multiple_components
759
_multiple_components = n
761
set_multiple_components (0)
763
def has_multiple_components () :
764
global _multiple_components
765
return _multiple_components
767
def draw3_trigger ( ) :
768
"arrange to call draw3 when everything else is finished"
769
global _draw3_changes
771
set_idler ( _draw3_idler )
775
"clear3 ( ) : Clear the current 3D display list."
776
global _draw3_list, _draw3_n
777
_draw3_list [_draw3_n:] = []
778
set_multiple_components (0)
780
def window3 ( * n , **kw ) :
783
window3 ( ) or window3 (n)
784
initialize style="nobox.gs" window for 3D graphics
787
if kw.has_key ("dump") :
791
if kw.has_key ("hcp") :
793
window (wait=1, style="nobox.gs", legends=0, hcp=kw ["hcp"],
797
window (n [0], wait=1, style="nobox.gs", legends=0, hcp=kw ["hcp"],
802
window (wait=1, style="nobox.gs", legends=0)
804
window (n [0], wait=1, style="nobox.gs", legends=0)
806
def sort3d (z, npolys) :
810
given Z and NPOLYS, with len(Z)==sum(npolys,axis=0), return
811
a 2-element list [LIST, VLIST] such that Z[VLIST] and NPOLYS[LIST] are
812
sorted from smallest average Z to largest average Z, where
813
the averages are taken over the clusters of length NPOLYS.
814
Within each cluster (polygon), the cyclic order of Z[VLIST]
815
remains unchanged, but the absolute order may change.
817
This sorting order produces correct or nearly correct order
818
for a plfp command to make a plot involving hidden or partially
819
hidden surfaces in three dimensions. It works best when the
820
polys form a set of disjoint closed, convex surfaces, and when
821
the surface normal changes only very little between neighboring
822
polys. (If the latter condition holds, then even if sort3d
823
mis-orders two neighboring polys, their colors will be very
824
nearly the same, and the mistake won't be noticeable.) A truly
825
correct 3D sorting routine is impossible, since there may be no
826
rendering order which produces correct surface hiding (some polys
827
may need to be split into pieces in order to do that). There
828
are more nearly correct algorithms than this, but they are much
833
# first compute z, the z-centroid of every poly
834
# get a list the same length as x, y, or z which is 1 for each
835
# vertex of poly 1, 2 for each vertex of poly2, etc.
836
# the goal is to make nlist with histogram(nlist)==npolys
837
nlist = histogram(cumsum (npolys,axis=0)) [0:-1]
838
nlist = cumsum (nlist,axis=0)
839
# now sum the vertex values and divide by the number of vertices
840
z = histogram (nlist, z) / npolys
842
# sort the polygons from smallest z to largest z
843
list = index_sort (z)
844
# next, find the list which sorts the polygon vertices
845
# first, find a list vlist such that sort(vlist) is above list
846
vlist = zeros (len (list), Int)
847
array_set (vlist, list, arange (len (list), dtype = Int))
848
# then reset the nlist values to that pre-sorted order, so that
849
# sort(nlist) will be the required vertex sorting list
850
nlist = take(vlist, nlist, 0)
851
# the final hitch is to ensure that the vertices within each polygon
852
# remain in their initial order (sort scrambles equal values)
853
# since the vertices of a polygon can be cyclically permuted,
854
# it suffices to add a sawtooth function to a scaled nlist to
855
# produce a list in which each cluster of equal values will retain
856
# the same cyclic order after the sort
857
# (note that the more complicated msort routine would leave the
858
# clusters without even a cyclic permutation, if that were
860
n1max = max (npolys) # this must never be so large that
861
# numberof(npolys)*nmax > 2e9
862
nmax = n1max * ones (len (nlist), Int)
863
vlist = index_sort (nmax * nlist +
864
arange (len (nlist), dtype = Int) % n1max)
865
# primary sort key ^ secondary key ^
868
_square = 1 # Global variable which tells whether to force equal axes
870
_yfactor = 1. # These globals enable one to scale one or both axes up or down
872
def get_factors_ ( ) :
873
return [_xfactor, _yfactor]
875
def get_square_ ( ) :
879
def limits_ (square = 0, yfactor = 1., xfactor = 1.) :
880
global _square, _xfactor, _yfactor
885
def draw3 (called_as_idler = 0, lims = None) :
888
draw3 (called_as_idler = 0, lims = None):
889
Draw the current 3d display list.
890
Ordinarily triggered automatically when the drawing changes.
893
global _draw3, _draw3_changes, _draw3_list, _draw3_n, _gnomon
898
# the first _draw3_n elements of _draw3_list are the viewing
899
# transforms, lighting, etc.
900
# thereafter, elements are (function, argument-list) pairs
901
# the _draw3 flag alerts the functions that these are the draw
902
# calls rather than the interactive setup calls
904
list = _draw3_list [_draw3_n:]
905
no_lims = lims == None
907
# ZCM Feb. 1997: Because Gist command 'limits' seems to
908
# misbehave and be timing dependent, I have added the kludge
909
# below, which seems to make things work.
914
lims = fnc (list [1])
918
if fv != None and lims != None :
919
lims = [min (fv [0], lims [0]),
920
max (fv [1], lims [1]),
921
min (fv [2], lims [2]),
922
max (fv [3], lims [3])]
930
_draw3_changes = None
938
_draw3_view = [array ([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), [0., 0., 0.], None]
939
_draw3_nv = len (_draw3_view)
951
dummy = _light3_ambient
953
_light3_ambient = 0.2
956
dummy = _light3_diffuse
958
_light3_diffuse = 1.0
961
dummy = _light3_specular
963
_light3_specular = 0.0
966
dummy = _light3_spower
973
_light3_sdir = array ( [1.0, 0.5, 1.0]) / sqrt(2.25)
975
_light3_list = [_light3_ambient, _light3_diffuse, _light3_specular,
976
_light3_spower, _light3_sdir]
979
_draw3_list = _draw3_view + _light3_list
980
_draw3_n = len (_draw3_list)
982
def get_draw3_list_ ( ) :
986
def get_draw3_n_ ( ) :
995
def set_default_gnomon ( * n ) :
996
# The default gnomon value is used when _draw3 is nonzero, i. e.,
997
# when a plot is actually done after every plot call.
998
global _default_gnomon
1004
set_default_gnomon (0)
1006
def gnomon (* on, ** kw) :
1011
Toggle the gnomon display. If on is present and non-zero,
1012
turn on the gnomon. If zero, turn it off.
1014
The gnomon shows the X, Y, and Z axis directions in the
1015
object coordinate system. The directions are labeled.
1016
The gnomon is always infinitely far behind the object
1017
(away from the camera).
1019
There is a mirror-through-the-screen-plane ambiguity in the
1020
display which is resolved in two ways: (1) the (X, Y, Z)
1021
coordinate system is right-handed, and (2) If the tip of an
1022
axis projects into the screen, its label is drawn in opposite
1023
polarity to the other text in the screen.
1026
# (ZCM 4/4/97) Add keyword argument chr to allow specification
1027
# of the axis labels.
1032
_gnomon = 1 - _gnomon
1039
if kw.has_key ("chr") :
1042
chr = ["X", "Y", "Z"]
1044
def _gnomon_draw ( ) :
1046
o = array ( [0., 0., 0.], Float)
1047
x1 = array ( [1., 0., 0.], Float)
1048
y1 = array ( [0., 1., 0.], Float)
1049
z1 = array ( [0., 0., 1.], Float)
1050
xyz1 = array (getrot3_ ( ), copy = 1)
1051
xyz2 = array([[o,x1],[o,y1],[o,z1]])
1054
xyz = zeros ( (s2 [1], s2 [0], s1 [1] ), Float)
1055
xyz [0, :, :] = dot (transpose (xyz1), xyz2 [:, 0, :])
1056
xyz [1, :, :] = dot (transpose (xyz1), xyz2 [:, 1, :])
1057
xyz = .0013 * _gnomon_scale * xyz
1058
x1 = xyz [0:2, 0, 0:3]
1059
y1 = xyz [0:2, 1, 0:3]
1060
z1 = xyz [1, 2, 0:3]
1065
wid = min (_gnomon_scale / 18., 6.)
1066
if ( wid < 0.5 ) : wid = 0.
1068
pldj (x0 + _gnomon_x, y0 + _gnomon_y, x1 + _gnomon_x, y1 + _gnomon_y,
1069
width = wid, type = 1, legend = "")
1072
# Compute point size of labels (1/3 of axis length)
1073
pts = [8, 10, 12, 14, 18, 24] [digitize (_gnomon_scale / 3.0,
1074
array ([9, 11, 13, 16, 21], Int))]
1076
if _gnomon_scale < 21.0 :
1077
x1 = x1 * 21. / _gnomon_scale
1078
y1 = y1 * 21. / _gnomon_scale
1079
# label positions: first find shortest axis
1080
xy = sqrt (x1 * x1 + y1 * y1)
1081
xysum = add.reduce (xy)
1082
i = argmin (xy,axis=-1) # mnx (xy)
1083
jk = [ [1, 2], [2, 0], [0, 1]] [i]
1086
if xy [i] < 1.e-7 * xysum : # guarantee not exactly zero
1087
x1 [i] = -1.e-6 * (x1 [j] + x1 [k] )
1088
y1 [i] = -1.e-6 * (y1 [j] + y1 [k] )
1089
xy [i] = sqrt (x1 [i] * x1 [i] + y1 [i] * y1 [i])
1091
# next find axis nearest to shortest
1092
if abs (x1 [j] * y1 [i] - y1 [j] * x1 [i]) * xy [k] > \
1093
abs (x1 [k] * y1 [i] - y1 [k] * x1 [i]) * xy [j] :
1097
# furthest axis first--move perpendicular to nearest axis
1100
xy = sqrt (xk * xk + yk * yk)
1103
if (xk * x1 [k] + yk * y1 [k] < 0.0 ) :
1106
# nearer axis next--move perpendicular to furthest axis
1109
xy = sqrt (xj * xj + yj * yj)
1112
if (xj * x1[j] + yj * y1 [j] < 0.0 ) :
1115
# shortest axis last -- move perpendicular to nearer
1118
xy = sqrt (xi * xi + yi * yi)
1121
if (xi *x1 [i] + yi * y1 [i] < 0.0) :
1125
# shortest axis label may need adjustment
1128
# just center it in correct quadrant
1129
jk = sign_ (xi * xj + yi * yj)
1130
yi = sign_ (xi * xk + yi * yk)
1131
xi = jk * xj + yi * xk
1132
yi = jk * yj + yi * yk
1133
jk = sqrt (xi * xi + yi * yi)
1136
x = zeros (3, Float)
1137
y = zeros (3, Float)
1146
x = x + x1 + _gnomon_x
1147
y = y + y1 + _gnomon_y
1151
chr = ["X", "Y", "Z"]
1152
gnomon_text_ (chr [i], x [i], y [i], pts, z1 [i] < 1.e-6)
1153
gnomon_text_ (chr [j], x [j], y [j], pts, z1 [j] < 1.e-6)
1154
gnomon_text_ (chr [k], x [k], y [k], pts, z1 [k] < 1.e-6)
1157
dummy = _gnomon_scale
1159
_gnomon_scale = 30. # axes lengths in points
1163
_gnomon_x = 0.18 # gnomon origin in system 0 coordinates
1169
def gnomon_text_ (chr, x, y, pts, invert) :
1170
# pts = 8, 10, 12, 14, 18, or 24
1174
plg (array ( [y, y]), array ( [x, x]), type = 1, width = 2.2 * pts,
1175
color = col, marks = 0, legend = "")
1178
plt (chr, x, y, justify = "CH", color = col, height = pts,
1179
font = "helvetica", opaque = 0)
1185
def spin3 (nframes = 30, axis = array ([-1, 1, 0], Float), tlimit = 60.,
1186
dtmin = 0.0, bracket_time = array ([2., 2.], Float), lims = None,
1187
timing = 0, angle = 2. * pi) :
1190
spin3 ( ) or spin3 (nframes) os spin3 (nframes, axis)
1191
Spin the current 3D display list about AXIS over NFRAMES. Keywords
1192
tlimit= the total time allowed for the movie in seconds (default 60),
1193
dtmin= the minimum allowed interframe time in seconds (default 0.0),
1194
bracket_time= (as for movie function in movie.i), timing = 1 if
1195
you want timing measured and printed out, 0 if not.
1197
The default AXIS is [-1,1,0] and the default NFRAMES is 30.
1201
# Note on global variables (ZCM 2/21/97):
1202
# I see no better way of sharing these between spin3 and _spin3
1203
# than making them global. Otherwise one would have to pass
1204
# them to movie, which would then send them as arguments to
1205
# _spin3. But because movie may call other routines, every one
1206
# of them would have to have these values, necessary or not.
1207
# So I have started their names with underscores; at least
1208
# this makes them inaccessible outside this module.
1209
global _phi, _theta, _dtheta
1211
_g_nframes = nframes
1212
_dtheta = angle / (nframes - 1)
1213
_theta = arccos (axis [2] / sqrt (axis [0] * axis [0] + axis [1] * axis [1] +
1214
axis [2] * axis [2]))
1215
inc = axis [0] == axis [1] == 0
1216
_phi = arctan2 (axis [1], axis [0] + inc)
1218
movie (_spin3, tlimit, dtmin, bracket_time, lims, timing = 0)
1223
global _phi, _theta, _dtheta
1227
rot3 (ya = -_theta, za = _dtheta)
1228
rot3 (ya = _theta, za = _phi)
1230
limits (lims [0], lims [1], lims [2], lims [3])
added
removed
scipy/sandbox/xplt/pl3d.py
