## Automatically adapted for scipy Oct 31, 2005 by # $Id: pl3d.py 2183 2006-08-29 10:30:44Z oliphant $ # Copyright (c) 1996, 1997, The Regents of the University of California. # All rights reserved. See Legal.htm for full text and disclaimer. from scipy import * from shapetest import * from yorick import * from numpy import * from gistfuncs import * # PL3D.PY # Viewing transforms and other aids for 3D plotting. # # $Id: pl3d.py 2183 2006-08-29 10:30:44Z oliphant $ # Copyright (c) 1997. The Regents of the University of California. # All rights reserved. """ General overview of module pl3d: (1) Viewing transform machinery. Arguably the simplest model is the CAD/CAM notion that the object you see is oriented as you see it in the current picture. You can then move it left, right, up, down, or toward or away from you, or you can rotate it about any of the three axes (horizontal, vertical, or out of the screen). The xyz coordinates of the object remains unchanged throughout all of this, but this object coordinate system changes relative to the fixed xyz of the viewer, in which x is always to the right, y is up, and z is directed out of the screen. Initially, the two coordinate systems coincide. rot3 (xangle,yangle,zangle) Rotate the object about viewer's x-axis by xangle, then about viewer's y-axis by yangle, then about viewer's z-axis by zangle mov3 (xchange,ychange,zchange) Move the object by the specified amounts. setz3 (zcamera) The "camera" is located at (0,0,zcamera) in the viewer's coordinate system, looking in the minus-z direction. Initially, zcamera is very large, and the magnification factor is correspondingly large, giving an isometric view. Decreasing zcamera makes the perspective more extreme. If parts of the object are behind the camera, strange things may happen. undo3 () undo3 (n) Undo the last N (default 1) viewpoint commands (rot3, mov3, or setz3). Up to 100 viewpoint changes are remembered. viewpoint= save3() ... restore3 (viewpoint) The current viewpoint transformation can be saved and later restored. gnomon (on_off) Toggle the gnomon (a simple display showing the orientation of the xyz axes of the object). """ # ------------------------------------------------------------------------ def set_draw3_ ( n ) : """ set_draw3_ ( 0 | 1 ) is used to set the global draw3_, which controls whether the function draw3 actually shows a drawing. """ global _draw3 _draw3 = n def setrot3_ (x) : # ZCM 2/21/97 change reflects the fact that I hadn't realized # that car and cdr, as functions, return the item replaced. global _draw3_list oldx = _draw3_list [0] _draw3_list [0] = x undo3_set_ (setrot3_, oldx) def rot3 (xa = 0., ya = 0., za = 0.) : """ rot3 (xa, ya, za) rotate the current 3D plot by XA about viewer's x-axis, YA about viewer's y-axis, and ZA about viewer's z-axis. SEE ALSO: orient3, mov3, aim3, setz3, undo3, save3, restore3, light3 """ x = array ([1.,0.,0.], Float) y = array ([0.,1.,0.], Float) z = array ([0.,0.,1.], Float) [x, y] = rot3_ (za, x, y) [z, x] = rot3_ (ya, z, x) [y, z] = rot3_ (xa, y, z) # n. b. matrixMultiply has the unfortunate effect of destroying # the matrix that calls it. gr3 = array (getrot3_ (), copy = 1) setrot3_ (transpose (dot (transpose (gr3), array ( [x, y, z])))) def rot3_ (a, x, y) : ca = cos (a) sa = sin (a) return [multiply (ca, x) + multiply (sa, y), multiply (-sa, x) + multiply (ca, y)] def mov3 ( xa = 0., ya = 0., za = 0. ) : """ mov3 ( [xa [, ya [, za]]]) move the current 3D plot by XA along the viewer's x axis, YA along the viewer's y axis, and ZA along the viewer's z axis. SEE ALSO: rot3, orient3, setz3, undo3, save3, restore3, light3 """ gr = dot (transpose (gr), transpose (xa)) setorg3_ ( getorg3_ () - gr) def aim3 ( xa = 0., ya = 0., za = 0. ) : """ aim3 ( [xa [, ya [, za]]]) move the current 3D plot to put the point (XA, YA, ZA) in object coordinates at the point (0, 0, 0) -- the aim point -- in the viewer's coordinates. If any of the XA, YA, or ZA is nil, it defaults SEE ALSO: mov3, rot3, orient3, setz3, undo3, save3, restore3, light3 """ setorg3_ (x) _ZcError = "ZcError" def setz3 ( zc = None ) : """ setz3 ( [zc] ) Set the camera position to z = ZC (x = y = 0) in the viewer's coordinate system. If zc is None, set the camera to infinity (default). SEE ALSO: rot3, orient3, undo3, save3, restore3, light3 """ if not is_scalar (zc) : raise _ZcError, "camera position must be scalar." setzc3_ (zc) def orient3 ( ** kw ) : """ orient3 ( [phi = val1, theta = val2] ) Set the orientation of the object to (PHI, THETA). Orientations are a subset of the possible rotation matrices in which the z axis of the object appears vertical on the screen (that is, the object z axis projects onto the viewer y axis). The THETA angle is the angle from the viewer y axis to the object z axis, positive if the object z axis is tilted towards you (toward viewer +z). PHI is zero when the object x axis coincides with the viewer x axis. If neither PHI nor THETA is specified, PHI defaults to - pi / 4 and THETA defaults to pi / 6. If only PHI is specified, THETA remains unchanged, unless the current THETA is near pi / 2, in which case THETA returns to pi / 6, or unless the current orientation does not have a vertical z axis, in which case THETA returns to its default. Unlike rot3, orient3 is not a cumulative operation. """ # Notes with regard to global variables: (ZCM 2/21/97) # _orient3_phi, _orient3_theta, the default orientation angles, # are known and referred to only in this routine. I have started # them with an underscore, too, to make them inaccessible # from outside this module. # phi and theta need not be global here since they are recalculated # each time this routine is called. global _orient3_phi, _orient3_theta try : dummy = _orient3_theta except : _orient3_theta = pi / 6. try : dummy = _orient3_phi except : _orient3_phi = - pi / 4. if kw.has_key ("phi") and kw ["phi"] == None : kw ["phi"] = _orient3_phi if kw.has_key ("theta") and kw ["theta"] == None : kw ["theta"] = _orient3_theta if not kw.has_key ("phi") and not kw.has_key ("theta") : phi = _orient3_phi theta = _orient3_theta elif not kw.has_key ("phi") or not kw.has_key ("theta") : gr3 = array (getrot3_ (), copy = 1) z = dot (transpose (gr3), array ( [0., 0., 1.])) if abs (z [0]) > 1.e-6 : # object z-axis not aligned with viewer y-axis if not kw.has_key ("theta") : theta = _orient3_theta phi = kw ["phi"] else : phi = _orient3_phi theta = kw ["theta"] elif not kw.has_key ("theta") : phi = kw ["phi"] if (abs (z [1]) < 1.e-6) : theta = _orient3_theta else : theta = arctan2 (z [2], z [1]) else : theta = kw ["theta"] y = array ( [0., z [2], -z [1]]) x = dot (transpose (gr3), array ( [1., 0., 0.])) phi = arctan2 (sum (y * x,axis=0), x [0]) else : phi = kw ["phi"] theta = kw ["theta"] x = array ( [1., 0., 0.], Float) y = array ( [0., 1., 0.], Float) z = array ( [0., 0., 1.], Float) [y, z] = rot3_ (theta, y, z) [z, x] = rot3_ (phi, z, x) setrot3_ (array ( [x, -z, y], Float)) import copy def save3 ( ) : """ view = save3 ( ) Save the current 3D viewing transformation and lighting. Actually, this doesn't save anything; it returns a copy of the current 3D viewing transformation and lighting, so that the user can put it aside somewhere. SEE ALSO: restore3, rot3, mov3, aim3, light3 """ return _draw3_list [0:_draw3_n] def restore3 ( view = None ) : """ restore3 ( view ) Restore a previously saved 3D viewing transformation and lighting. If view is missing, rotate object to viewer's coordinate system. SEE ALSO: restore3, rot3, mov3, aim3, light3 """ global _draw3_list, _draw3_view, _light3_list, _draw3_n if view != None : view = view [0:len (view)] # Copies view else : view = _draw3_view + _light3_list old = _draw3_list [0:_draw3_n] _draw3_list = view [0:_draw3_n] + _draw3_list [_draw3_n:] undo3_set_ (restore3, old) _AmbientError = "AmbientError" _DiffuseError = "DiffuseError" _LightingError = "LightingError" def light3 ( * kw, ** kwds ) : """ light3 (ambient=a_level, diffuse=d_level, specular=s_level, spower=n, sdir=xyz) Sets lighting properties for 3D shading effects. A surface will be shaded according to its to its orientation relative to the viewing direction. The ambient level A_LEVEL is a light level (arbitrary units) that is added to every surface independent of its orientation. The diffuse level D_LEVEL is a light level which is proportional to cos(theta), where theta is the angle between the surface normal and the viewing direction, so that surfaces directly facing the viewer are bright, while surfaces viewed edge on are unlit (and surfaces facing away, if drawn, are shaded as if they faced the viewer). The specular level S_LEVEL is a light level proportional to a high power spower=N of 1+cos(alpha), where alpha is the angle between the specular reflection angle and the viewing direction. The light source for the calculation of alpha lies in the direction XYZ (a 3 element vector) in the viewer's coordinate system at infinite distance. You can have ns light sources by making S_LEVEL, N, and XYZ (or any combination) be vectors of length ns (3-by-ns in the case of XYZ). (See source code for specular_hook function definition if powers of 1+cos(alpha) aren't good enough for you.) With no arguments, return to the default lighting. EXAMPLES: light3 ( diffuse=.1, specular=1., sdir=[0,0,-1]) (dramatic "tail lighting" effect) light3 ( diffuse=.5, specular=1., sdir=[1,.5,1]) (classic "over your right shoulder" lighting) light3 ( ambient=.1,diffuse=.1,specular=1., sdir=[[0,0,-1],[1,.5,1]],spower=[4,2]) (two light sources combining previous effects) SEE ALSO: rot3, save3, restore3 """ global _draw3_list, _draw3_nv if len (kw) > 0 : kwds = kw [0] old = _draw3_list [_draw3_nv:] [0:5] flags = 0 if kwds.has_key ("ambient") and kwds ["ambient"] != None : ambient = kwds ["ambient"] if not is_scalar (ambient) : raise _AmbientError, "ambient light level must be scalar." flags = flags | 1 _draw3_list [_draw3_nv] = ambient if kwds.has_key ("diffuse") and kwds ["diffuse"] != None : diffuse = kwds ["diffuse"] if not is_scalar (diffuse) : raise _DiffuseError, "diffuse light level must be scalar." flags = flags | 2 _draw3_list [_draw3_nv + 1 ] = diffuse if kwds.has_key ("specular") and kwds ["specular"] != None : specular = kwds ["specular"] flags = flags | 4 else : specular = _draw3_list [_draw3_nv + 2] if kwds.has_key ("spower") and kwds ["spower"] != None : spower = kwds ["spower"] flags = flags | 8 else : spower = _draw3_list [_draw3_nv + 3] if kwds.has_key ("sdir") and kwds ["sdir"] != None : sdir = kwds ["sdir"] dims = shape (sdir) if dims == 0 or len (dims) == 2 and dims [1] != 3 : raise _LightingError, \ "lighting direction must be 3 vector or ns-by-3 array." flags = flags | 16 else : sdir = _draw3_list [_draw3_nv + 4] if flags & 28 : if flags & 4 : _draw3_list [_draw3_nv + 2] = specular if flags & 8 : _draw3_list [_draw3_nv + 3] = spower if flags & 16 : _draw3_list [_draw3_nv + 4] = sdir if not flags : _draw3_list [_draw3_nv: _draw3_nv + 5] = _light3_list [0:5] undo3_set_ (light3_, old) def light3_ (arg) : global _draw3_list, _draw3_nv _draw3_list [_draw3_nv:_draw3_nv + 5] = arg [0:5] def get3_light (xyz, * nxyz) : """ get3_light(xyz, nxyz) or get3_light(xyz) return 3D lighting for polygons with vertices XYZ. If NXYZ is specified, XYZ should be sum(nxyz,axis=0)-by-3, with NXYZ being the list of numbers of vertices for each polygon (as for the plfp function). If NXYZ is not specified, XYZ should be a quadrilateral mesh, ni-by-nj-by-3 (as for the plf function). In the first case, the return value is len (NXYZ) long; in the second case, the return value is (ni-1)-by-(nj-1). The parameters of the lighting calculation are set by the light3 function. SEE ALSO: light3, set3_object, get3_normal, get3_centroid """ global _draw3_list, _draw3_nv list = _draw3_list [_draw3_nv:] ambient = list [0] diffuse = list [1] specular = list [2] spower = list [3] sdir = list [4] if len (nxyz) == 0 : normal = get3_normal (xyz) else : normal = get3_normal (xyz, nxyz [0]) zc = getzc3_ ( ) if ( not zc ) : view = array ( [0., 0., 1.], Float) elif len (nxyz) == 0 : view = array ( [0., 0., zc], Float) - get3_centroid (xyz) else : view = array ( [0., 0., zc], Float) - get3_centroid (xyz, nxyz [0]) m1 = \ sqrt ( sum (view * view,axis=0)) if m1 == 0. : m1 = 1. view = view / m1 nv = normal [0, ...] * view [0] + normal [1, ...] * view [1] + \ normal [2, ...] * view [2] light = ambient + diffuse * abs (nv) if specular != 0. : sv = transpose (transpose (sdir) / sqrt (sum (transpose (sdir*sdir),axis=0))) sv = dot (sv, view) if len (shape (sdir)) == 1 : sn = sum(array([sdir[0]*normal[0],sdir[1]*normal[1], sdir[2]*normal[2]]),axis=0) ####### I left out the specular_hook stuff. m1 = maximum (sn * nv -0.5 * sv + 0.5, 1.e-30) m1 = m1 ** spower light = light + (specular * m1) elif len (shape (sdir)) >= 2 : # multiple light sources nsrc = len (shape (sdir)) for i in range (nsrc) : sn = sum(array([sdir[i,0]*normal[0],sdir[i,1]*normal[1], sdir[i,2]*normal[2]]),axis=0) m1 = maximum (sn * nv -0.5 * sv [i] + 0.5, 1.e-30) ** spower [i] light = light + specular * m1 return light def get3_normal (xyz, *nxyz) : """ get3_normal(xyz, nxyz) or get3_normal(xyz) return 3D normals for polygons with vertices XYZ. If NXYZ is specified, XYZ should be sum(nxyz,axis=0)-by-3, with NXYZ being the list of numbers of vertices for each polygon (as for the plfp function). If NXYZ is not specified, XYZ should be a quadrilateral mesh, ni-by-nj-by-3 (as for the plf function). In the first case, the return value is len(NXYZ)-by-3; in the second case, the return value is (ni-1)-by-(nj-1)-by-3. The normals are constructed from the cross product of the lines joining the midpoints of two edges which as nearly quarter the polygon as possible (the medians for a quadrilateral). No check is made that these not be parallel; the returned "normal" is [0,0,0] in that case. Also, if the polygon vertices are not coplanar, the "normal" has no precisely definable meaning. SEE ALSO: get3_centroid, get3_light """ if len (nxyz) == 0 : # if no polygon list is given, assume xyz is 2D mesh # form normal as cross product of medians m1 = dif_ (zcen_ (xyz, 1), 2) m2 = zcen_ (dif_ (xyz, 1), 2) else : # with polygon list, more elaborate calculation required # (1) frst subscripts the first vertex of each polygon frst = cumsum (nxyz [0],axis=0) - nxyz [0] # form normal by getting two approximate diameters # (reduces to above medians for quads) # (2) compute midpoints of first three sides n2 = (nxyz [0] + 1) / 2 c0 = (take(xyz, frst, 0) + take(xyz, frst + 1, 0)) / 2. i = frst + n2 - 1 c1 = (take(xyz, i, 0) + take(xyz, i + 1, 0)) / 2. i = n2 / 2 c2 = (take(xyz, frst + i, 0) + take(xyz, frst + (i + 1) % nxyz [0], 0)) / 2. i = minimum (i + n2, nxyz [0]) - 1 c3 = (take(xyz, frst + i, 0) + take(xyz, frst + (i + 1) % nxyz [0], 0)) / 2. m1 = c1 - c0 m2 = c3 - c2 # poly normal is cross product of two medians (or diameters) # normal = m1; I had to reverse the sign. if len (shape (xyz)) == 3 : n1 = m1 [2, :] * m2 [1, :] - \ m1 [1, :] * m2 [2, :] n2 = m1 [0, :] * m2 [2, :] - \ m1 [2, :] * m2 [0, :] n3 = m1 [1, :] * m2 [0, :] - \ m1 [0, :] * m2 [1, :] else : n1 = m1 [:, 2] * m2 [:, 1] - \ m1 [:, 1] * m2 [:, 2] n2 = m1 [:, 0] * m2 [:, 2] - \ m1 [:, 2] * m2 [:, 0] n3 = m1 [:, 1] * m2 [:, 0] - \ m1 [:, 0] * m2 [:, 1] m1 = sqrt (n1 ** 2 + n2 **2 + n3 **2) m1 = m1 + equal (m1, 0.0) normal = array([n1 / m1, n2 / m1, n3 / m1]) return normal def get3_centroid (xyz, * nxyz) : """ get3_centroid(xyz, *nxyz) or get3_centroid(xyz) return 3D centroids for polygons with vertices XYZ. If NXYZ is specified, XYZ should be sum(nxyz,axis=0)-by-3, with NXYZ being the list of numbers of vertices for each polygon (as for the plfp function). If NXYZ is not specified, XYZ should be a quadrilateral mesh, ni-by-nj-by-3 (as for the plf function). In the first case, the return value is len(NXYZ) in length; in the second case, the return value is (ni-1)-by-(nj-1)-by-3. The centroids are constructed as the mean value of all vertices of each polygon. SEE ALSO: get3_normal, get3_light """ if len (nxyz) == 0 : # if no polygon list is given, assume xyz is 2D mesh centroid = zcen_ (zcen_ (xyz, 1), 0) else : # with polygon list, more elaborate calculation required last = cumsum (nxyz [0],axis=0) list = histogram (1 + last) [0:-1] list = cumsum (list,axis=0) k = len (nxyz [0]) l = shape (xyz) [0] centroid = zeros ( (k, 3)) centroid [0:k, 0] = histogram (list, xyz [0:l,0]) centroid [0:k, 1] = histogram (list, xyz [0:l,1]) centroid [0:k, 2] = histogram (list, xyz [0:l,2]) fnxyz = array (nxyz [0], Float ) centroid = centroid / fnxyz return centroid _Get3Error = "Get3Error" def get3_xy (xyz, *flg) : """ get3_xy (xyz) or get3_xy(xyz, 1) Given anything-by-3 coordinates XYZ, return X and Y in viewer's coordinate system (set by rot3, mov3, orient3, etc.). If the second argument is present and non-zero, also return Z (for use in sort3d or get3_light, for example). If the camera position has been set to a finite distance with setz3, the returned coordinates will be tangents of angles for a perspective drawing (and Z will be scaled by 1/zc). Unlike the Yorick version, this function returns a 3-by-anything array of coordinates. Actually, what it returns is a 3-by-anything python array, whose 0th element is the x array, whose 1th element is the y array, and whose 2th element is the z array if asked for. I believe that x, y, and z can be either 1d or 2d, so this routine is written in two cases. """ # rotate and translate to viewer's coordinate system shp = shape (xyz) if len (shp) == 3: # 2d mesh case is much more complex than in Yorick (k, l) = shp [1:3] go3_ = getorg3_ () # Unwind xyz xx = ravel (xyz [0]) yy = ravel (xyz [1]) zz = ravel (xyz [2]) tmpxyz = array ( [xx, yy, zz]) gr3 = array (getrot3_ (), copy = 1) tmpxyz = dot (transpose (gr3), tmpxyz - array ( [ [go3_ [0]], [go3_ [1]], [go3_ [2]]])) ## xx = transpose (reshape (ravel (tmpxyz [0]), (k,l))) ## yy = transpose (reshape (ravel (tmpxyz [1]), (k,l))) ## zz = transpose (reshape (ravel (tmpxyz [2]), (k,l))) xx = (reshape (ravel (tmpxyz [0]), (k,l))) yy = (reshape (ravel (tmpxyz [1]), (k,l))) zz = (reshape (ravel (tmpxyz [2]), (k,l))) tmpxyz = array ( [xx, yy, zz]) elif len (shp) == 2: go3_ = getorg3_ () lm = array (getrot3_ (), copy = 1) rm = (xyz - array ( [ go3_ [0], go3_ [1], go3_ [2]])) tmpxyz = dot (rm, lm) else: raise _Get3Error, "xyz has a bad shape: " + `shp` # do optional perspective projection zc = getzc3_ () if zc != None : if len (shp) == 2 : z = tmpxyz [:, 2] zc = maximum (zc - z, 1.e-35) # protect behind camera, avoid zero divide tmpxyz [:, 0] = tmpxyz [:, 0] / zc tmpxyz [:, 1] = tmpxyz [:, 1] / zc if len (flg) != 0 and flg [0] != 0 : tmpxyz [:, 2] = tmpxyz [:, 2] / zc elif len (shp) == 3 : z = tmpxyz [:,:, 2] zc = maximum (zc - z, 1.e-35) # protect behind camera, avoid zero divide tmpxyz [:,:, 0] = tmpxyz [:,:, 0] / zc tmpxyz [:,:, 1] = tmpxyz [:,:, 1] / zc if len (flg) != 0 and flg [0] != 0 : tmpxyz [:,:, 2] = tmpxyz [:,:, 2] / zc return tmpxyz _UndoError = "UndoError" _in_undo3 = 0 _undo3_list = [] def undo3 (n = 1) : """ undo3 () or undo3 (n) Undo the effects of the last N (default 1) rot3, orient3, mov3, aim3, setz3, or light3 commands. """ global _in_undo3, _undo3_list n = 2 * n if n < 0 or n > len (_undo3_list) : raise _UndoError, "not that many items in undo list" _in_undo3 = 1 # flag to skip undo3_set_ # perhaps should save discarded items in a redo list? use_list = undo3_list [-n:] undo3_list = undo3_list [:-n] while n > 0 : fnc = use_list_ [0] del use_list_ [0] arg = use_list_ [0] del use_list_ [0] fnc (arg) n = n - 2 _in_undo3 = 0 draw3_trigger ( ) def set3_object (fnc, arg) : """ set3_object (drawing_function, [arg1,arg2,...]) set up to trigger a call to draw3, adding a call to the 3D display list of the form: DRAWING_FUNCTION ( [ARG1, ARG2, ...])) When draw3 calls DRAWING_FUNCTION, the external variable draw3_ will be non-zero, so DRAWING_FUNCTION can be written like this: def drawing_function(arg) : if (draw3_) : arg1= arg [0] arg1= arg [1] ... ...... ...... return ...... ...... set3_object (drawing_function, [arg1,arg2,...]) SEE ALSO: get3_xy, get3_light, sort3d """ global _draw3_list _draw3_list = _draw3_list + [fnc, arg] draw3_trigger () def setorg3_ ( x ) : # ZCM 2/21/97 change reflects the fact that I hadn't realized # that car and cdr, as functions, return the item replaced. global _draw3_list oldx = _draw3_list [1] _draw3_list [1] = x undo3_set_ ( setorg3_, oldx) def setzc3_ (x) : # ZCM 2/21/97 change reflects the fact that I hadn't realized # that car and cdr, as functions, return the item replaced. global _draw3_list oldx = _draw3_list [2] _draw3_list [2] = x undo3_set_ ( setzc3_, oldx) def getrot3_ () : return _draw3_list [0] def getorg3_ () : return _draw3_list [1] def getzc3_ () : return _draw3_list [2] def undo3_set_ (fnc, arg) : global _undo3_list, _in_undo3, _undo3_limit # arg = copy.deepcopy (arg) if not _in_undo3 : if len (_undo3_list) >= 2 * _undo3_limit : _undo3_list = _undo3_list [0:2 * _undo3_limit - 2] _undo3_list = [fnc, arg] + _undo3_list draw3_trigger ( ) _in_undo3 = 0 # ?????????????? _in_undo3 = 100 def do_nothing ( ) : pass return def clear_idler ( ) : _idler = do_nothing ( ) def set_idler ( fnc ) : global _idler _idler = fnc def call_idler ( ) : global _idler _idler ( ) def _draw3_idler ( ) : # I have added orientation and limits to this because they may not # have been set by a previous command. If the user doesn't like this, # he/she can write his/her own idler. (ZCM 7/1/97) global _default_gnomon orient3 () if current_window () == -1 : window3 (0) else : window3 (current_window ()) gnomon (_default_gnomon) lims = draw3 (1) if lims == None : return else : limits (lims [0], lims [1], lims [2], lims [3]) def set_default_idler ( ) : set_idler (_draw3_idler) set_default_idler ( ) _draw3_changes = None def set_multiple_components ( n = 0 ) : global _multiple_components _multiple_components = n set_multiple_components (0) def has_multiple_components () : global _multiple_components return _multiple_components def draw3_trigger ( ) : "arrange to call draw3 when everything else is finished" global _draw3_changes global _draw3_idler set_idler ( _draw3_idler ) _draw3_changes = 1 def clear3 ( ) : "clear3 ( ) : Clear the current 3D display list." global _draw3_list, _draw3_n _draw3_list [_draw3_n:] = [] set_multiple_components (0) def window3 ( * n , **kw ) : """ window3 ( ) or window3 (n) initialize style="nobox.gs" window for 3D graphics """ if kw.has_key ("dump") : dump = kw ["dump"] else : dump = 0 if kw.has_key ("hcp") : if len (n) == 0 : window (wait=1, style="nobox.gs", legends=0, hcp=kw ["hcp"], dump = dump) hcpon () else : window (n [0], wait=1, style="nobox.gs", legends=0, hcp=kw ["hcp"], dump = dump) hcpon () else : if len (n) == 0 : window (wait=1, style="nobox.gs", legends=0) else : window (n [0], wait=1, style="nobox.gs", legends=0) def sort3d (z, npolys) : """ sort3d(z, npolys) given Z and NPOLYS, with len(Z)==sum(npolys,axis=0), return a 2-element list [LIST, VLIST] such that Z[VLIST] and NPOLYS[LIST] are sorted from smallest average Z to largest average Z, where the averages are taken over the clusters of length NPOLYS. Within each cluster (polygon), the cyclic order of Z[VLIST] remains unchanged, but the absolute order may change. This sorting order produces correct or nearly correct order for a plfp command to make a plot involving hidden or partially hidden surfaces in three dimensions. It works best when the polys form a set of disjoint closed, convex surfaces, and when the surface normal changes only very little between neighboring polys. (If the latter condition holds, then even if sort3d mis-orders two neighboring polys, their colors will be very nearly the same, and the mistake won't be noticeable.) A truly correct 3D sorting routine is impossible, since there may be no rendering order which produces correct surface hiding (some polys may need to be split into pieces in order to do that). There are more nearly correct algorithms than this, but they are much slower. SEE ALSO: get3_xy """ # first compute z, the z-centroid of every poly # get a list the same length as x, y, or z which is 1 for each # vertex of poly 1, 2 for each vertex of poly2, etc. # the goal is to make nlist with histogram(nlist)==npolys nlist = histogram(cumsum (npolys,axis=0)) [0:-1] nlist = cumsum (nlist,axis=0) # now sum the vertex values and divide by the number of vertices z = histogram (nlist, z) / npolys # sort the polygons from smallest z to largest z list = index_sort (z) # next, find the list which sorts the polygon vertices # first, find a list vlist such that sort(vlist) is above list vlist = zeros (len (list), Int) array_set (vlist, list, arange (len (list), dtype = Int)) # then reset the nlist values to that pre-sorted order, so that # sort(nlist) will be the required vertex sorting list nlist = take(vlist, nlist, 0) # the final hitch is to ensure that the vertices within each polygon # remain in their initial order (sort scrambles equal values) # since the vertices of a polygon can be cyclically permuted, # it suffices to add a sawtooth function to a scaled nlist to # produce a list in which each cluster of equal values will retain # the same cyclic order after the sort # (note that the more complicated msort routine would leave the # clusters without even a cyclic permutation, if that were # necessary) n1max = max (npolys) # this must never be so large that # numberof(npolys)*nmax > 2e9 nmax = n1max * ones (len (nlist), Int) vlist = index_sort (nmax * nlist + arange (len (nlist), dtype = Int) % n1max) # primary sort key ^ secondary key ^ return [list, vlist] _square = 1 # Global variable which tells whether to force equal axes _xfactor = 1. _yfactor = 1. # These globals enable one to scale one or both axes up or down def get_factors_ ( ) : return [_xfactor, _yfactor] def get_square_ ( ) : global _square return _square def limits_ (square = 0, yfactor = 1., xfactor = 1.) : global _square, _xfactor, _yfactor _square = square _xfactor = xfactor _yfactor = yfactor def draw3 (called_as_idler = 0, lims = None) : """ draw3 (called_as_idler = 0, lims = None): Draw the current 3d display list. Ordinarily triggered automatically when the drawing changes. """ global _draw3, _draw3_changes, _draw3_list, _draw3_n, _gnomon if _draw3_changes : if called_as_idler : fma ( ) # the first _draw3_n elements of _draw3_list are the viewing # transforms, lighting, etc. # thereafter, elements are (function, argument-list) pairs # the _draw3 flag alerts the functions that these are the draw # calls rather than the interactive setup calls set_draw3_ (1) list = _draw3_list [_draw3_n:] no_lims = lims == None first = 1 # ZCM Feb. 1997: Because Gist command 'limits' seems to # misbehave and be timing dependent, I have added the kludge # below, which seems to make things work. while list != [] : fnc = list [0] if no_lims : if (first) : lims = fnc (list [1]) first = 0 else : fv = fnc (list [1]) if fv != None and lims != None : lims = [min (fv [0], lims [0]), max (fv [1], lims [1]), min (fv [2], lims [2]), max (fv [3], lims [3])] elif fv != None : lims = fv else : fnc (list [1]) list = list [2:] if _gnomon : _gnomon_draw ( ) _draw3_changes = None set_draw3_ (0) return lims # _draw3 = 0 try : dummy = _draw3_view except : _draw3_view = [array ([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), [0., 0., 0.], None] _draw3_nv = len (_draw3_view) try : dummy = _draw3 except : set_draw3_ (0) def get_draw3_ ( ) : global _draw3 return _draw3 try : dummy = _light3_ambient except : _light3_ambient = 0.2 try : dummy = _light3_diffuse except : _light3_diffuse = 1.0 try : dummy = _light3_specular except : _light3_specular = 0.0 try : dummy = _light3_spower except : _light3_spower = 2 try : dummy = _light3_sdir except : _light3_sdir = array ( [1.0, 0.5, 1.0]) / sqrt(2.25) _light3_list = [_light3_ambient, _light3_diffuse, _light3_specular, _light3_spower, _light3_sdir] _draw3_list = _draw3_view + _light3_list _draw3_n = len (_draw3_list) def get_draw3_list_ ( ) : global _draw3_list return _draw3_list def get_draw3_n_ ( ) : global _draw3_n return _draw3_n try : dummy = _gnomon except : _gnomon = 0 def set_default_gnomon ( * n ) : # The default gnomon value is used when _draw3 is nonzero, i. e., # when a plot is actually done after every plot call. global _default_gnomon if len (n) > 0 : _default_gnomon = n else : _default_gnomon = 0 set_default_gnomon (0) def gnomon (* on, ** kw) : """ gnomon () or gnomon (onoff) Toggle the gnomon display. If on is present and non-zero, turn on the gnomon. If zero, turn it off. The gnomon shows the X, Y, and Z axis directions in the object coordinate system. The directions are labeled. The gnomon is always infinitely far behind the object (away from the camera). There is a mirror-through-the-screen-plane ambiguity in the display which is resolved in two ways: (1) the (X, Y, Z) coordinate system is right-handed, and (2) If the tip of an axis projects into the screen, its label is drawn in opposite polarity to the other text in the screen. """ # (ZCM 4/4/97) Add keyword argument chr to allow specification # of the axis labels. global _gnomon, chr old = _gnomon if len (on) == 0 : _gnomon = 1 - _gnomon elif (on [0]) : _gnomon = 1 else : _gnomon = 0 if old != _gnomon : draw3_trigger () if kw.has_key ("chr") : chr = kw ["chr"] else : chr = ["X", "Y", "Z"] def _gnomon_draw ( ) : global chr o = array ( [0., 0., 0.], Float) x1 = array ( [1., 0., 0.], Float) y1 = array ( [0., 1., 0.], Float) z1 = array ( [0., 0., 1.], Float) xyz1 = array (getrot3_ ( ), copy = 1) xyz2 = array([[o,x1],[o,y1],[o,z1]]) s1 = shape ( xyz1 ) s2 = shape ( xyz2 ) xyz = zeros ( (s2 [1], s2 [0], s1 [1] ), Float) xyz [0, :, :] = dot (transpose (xyz1), xyz2 [:, 0, :]) xyz [1, :, :] = dot (transpose (xyz1), xyz2 [:, 1, :]) xyz = .0013 * _gnomon_scale * xyz x1 = xyz [0:2, 0, 0:3] y1 = xyz [0:2, 1, 0:3] z1 = xyz [1, 2, 0:3] x0 = x1 [0] x1 = x1 [1] y0 = y1 [0] y1 = y1 [1] wid = min (_gnomon_scale / 18., 6.) if ( wid < 0.5 ) : wid = 0. plsys (0) pldj (x0 + _gnomon_x, y0 + _gnomon_y, x1 + _gnomon_x, y1 + _gnomon_y, width = wid, type = 1, legend = "") plsys (1) # Compute point size of labels (1/3 of axis length) pts = [8, 10, 12, 14, 18, 24] [digitize (_gnomon_scale / 3.0, array ([9, 11, 13, 16, 21], Int))] if _gnomon_scale < 21.0 : x1 = x1 * 21. / _gnomon_scale y1 = y1 * 21. / _gnomon_scale # label positions: first find shortest axis xy = sqrt (x1 * x1 + y1 * y1) xysum = add.reduce (xy) i = argmin (xy,axis=-1) # mnx (xy) jk = [ [1, 2], [2, 0], [0, 1]] [i] j = jk [0] k = jk [1] if xy [i] < 1.e-7 * xysum : # guarantee not exactly zero x1 [i] = -1.e-6 * (x1 [j] + x1 [k] ) y1 [i] = -1.e-6 * (y1 [j] + y1 [k] ) xy [i] = sqrt (x1 [i] * x1 [i] + y1 [i] * y1 [i]) xyi = xy [i] # next find axis nearest to shortest if abs (x1 [j] * y1 [i] - y1 [j] * x1 [i]) * xy [k] > \ abs (x1 [k] * y1 [i] - y1 [k] * x1 [i]) * xy [j] : jk = j j = k k = jk # furthest axis first--move perpendicular to nearest axis xk = - y1 [j] yk = x1 [j] xy = sqrt (xk * xk + yk * yk) xk = xk / xy yk = yk / xy if (xk * x1 [k] + yk * y1 [k] < 0.0 ) : xk = - xk yk = - yk # nearer axis next--move perpendicular to furthest axis xj = - y1 [k] yj = x1 [k] xy = sqrt (xj * xj + yj * yj) xj = xj / xy yj = yj / xy if (xj * x1[j] + yj * y1 [j] < 0.0 ) : xj = - xj yj = - yj # shortest axis last -- move perpendicular to nearer xi = - y1 [j] yi = x1 [j] xy = sqrt (xi * xi + yi * yi) xi = xi / xy yi = yi / xy if (xi *x1 [i] + yi * y1 [i] < 0.0) : xi = - xi yi = - yi # shortest axis label may need adjustment d = 0.0013 * pts if xyi < d : # just center it in correct quadrant jk = sign_ (xi * xj + yi * yj) yi = sign_ (xi * xk + yi * yk) xi = jk * xj + yi * xk yi = jk * yj + yi * yk jk = sqrt (xi * xi + yi * yi) xi = xi / jk yi = yi / jk x = zeros (3, Float) y = zeros (3, Float) x [i] = xi x [j] = xj x [k] = xk y [i] = yi y [j] = yj y [k] = yk x = x * d y = y * d x = x + x1 + _gnomon_x y = y + y1 + _gnomon_y try : dum = chr except : chr = ["X", "Y", "Z"] gnomon_text_ (chr [i], x [i], y [i], pts, z1 [i] < 1.e-6) gnomon_text_ (chr [j], x [j], y [j], pts, z1 [j] < 1.e-6) gnomon_text_ (chr [k], x [k], y [k], pts, z1 [k] < 1.e-6) try : dummy = _gnomon_scale except : _gnomon_scale = 30. # axes lengths in points try : dummy = _gnomon_x except : _gnomon_x = 0.18 # gnomon origin in system 0 coordinates try : dummy = _gnomon_y except : _gnomon_y = 0.42 def gnomon_text_ (chr, x, y, pts, invert) : # pts = 8, 10, 12, 14, 18, or 24 col = "fg" if invert : plsys (0) plg (array ( [y, y]), array ( [x, x]), type = 1, width = 2.2 * pts, color = col, marks = 0, legend = "") plsys (1) col = "bg" plt (chr, x, y, justify = "CH", color = col, height = pts, font = "helvetica", opaque = 0) from movie import * g_nframes = 30 def spin3 (nframes = 30, axis = array ([-1, 1, 0], Float), tlimit = 60., dtmin = 0.0, bracket_time = array ([2., 2.], Float), lims = None, timing = 0, angle = 2. * pi) : """ spin3 ( ) or spin3 (nframes) os spin3 (nframes, axis) Spin the current 3D display list about AXIS over NFRAMES. Keywords tlimit= the total time allowed for the movie in seconds (default 60), dtmin= the minimum allowed interframe time in seconds (default 0.0), bracket_time= (as for movie function in movie.i), timing = 1 if you want timing measured and printed out, 0 if not. The default AXIS is [-1,1,0] and the default NFRAMES is 30. SEE ALSO: rot3 """ # Note on global variables (ZCM 2/21/97): # I see no better way of sharing these between spin3 and _spin3 # than making them global. Otherwise one would have to pass # them to movie, which would then send them as arguments to # _spin3. But because movie may call other routines, every one # of them would have to have these values, necessary or not. # So I have started their names with underscores; at least # this makes them inaccessible outside this module. global _phi, _theta, _dtheta global _g_nframes _g_nframes = nframes _dtheta = angle / (nframes - 1) _theta = arccos (axis [2] / sqrt (axis [0] * axis [0] + axis [1] * axis [1] + axis [2] * axis [2])) inc = axis [0] == axis [1] == 0 _phi = arctan2 (axis [1], axis [0] + inc) orig = save3 ( ) movie (_spin3, tlimit, dtmin, bracket_time, lims, timing = 0) restore3 (orig) def _spin3 (i) : global _g_nframes global _phi, _theta, _dtheta if i >= _g_nframes: return 0 rot3 (za = -_phi) rot3 (ya = -_theta, za = _dtheta) rot3 (ya = _theta, za = _phi) lims = draw3 ( ) limits (lims [0], lims [1], lims [2], lims [3]) return 1