1# Pizza.py toolkit, www.cs.sandia.gov/~sjplimp/pizza.html 2# Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories 3# 4# Copyright (2005) Sandia Corporation. Under the terms of Contract 5# DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains 6# certain rights in this software. This software is distributed under 7# the GNU General Public License. 8 9# gl tool 10 11oneline = "3d interactive visualization via OpenGL" 12 13docstr = """ 14g = gl(d) create OpenGL display for data in d 15 16 d = atom snapshot object (dump, data) 17 18g.bg("black") set background color (def = "black") 19g.size(N) set image size to NxN 20g.size(N,M) set image size to NxM 21g.rotate(60,135) view from z theta and azimuthal phi (def = 60,30) 22g.shift(x,y) translate by x,y pixels in view window (def = 0,0) 23g.zoom(0.5) scale image by factor (def = 1) 24g.box(0/1/2) 0/1/2 = none/variable/fixed box 25g.box(0/1/2,"green") set box color 26g.box(0/1/2,"red",4) set box edge thickness 27g.file = "image" file prefix for created images (def = "image") 28 29g.show(N) show image of snapshot at timestep N 30 31g.all() make images of all selected snapshots 32g.all(P) images of all, start file label at P 33g.all(N,M,P) make M images of snapshot N, start label at P 34 35g.pan(60,135,1.0,40,135,1.5) pan during all() operation 36g.pan() no pan during all() (default) 37 38 args = z theta, azimuthal phi, zoom factor at beginning and end 39 values at each step are interpolated between beginning and end values 40 41g.select = "$x > %g*3.0" string to pass to d.aselect.test() during all() 42g.select = "" no extra aselect (default) 43 44 %g varies from 0.0 to 1.0 from beginning to end of all() 45 46g.acol(2,"green") set atom colors by atom type (1-N) 47g.acol([2,4],["red","blue"]) 1st arg = one type or list of types 48g.acol(0,"blue") 2nd arg = one color or list of colors 49g.acol(range(20),["red","blue"]) if list lengths unequal, interpolate 50g.acol(range(10),"loop") assign colors in loop, randomly ordered 51 52 if 1st arg is 0, set all types to 2nd arg 53 if list of types has a 0 (e.g. range(10)), +1 is added to each value 54 interpolate means colors blend smoothly from one value to the next 55 56g.arad([1,2],[0.5,0.3]) set atom radii, same rules as acol() 57 58g.bcol() set bond color, same args as acol() 59g.brad() set bond thickness, same args as arad() 60 61g.tcol() set triangle color, same args as acol() 62g.tfill() set triangle fill, 0 fill, 1 line, 2 both 63 64g.lcol() set line color, same args as acol() 65g.lrad() set line thickness, same args as arad() 66 67g.adef() set atom/bond/tri/line properties to default 68g.bdef() default = "loop" for colors, 0.45 for radii 69g.tdef() default = 0.25 for bond/line thickness 70g.ldef() default = 0 fill 71 72 by default 100 types are assigned 73 if atom/bond/tri/line has type > # defined properties, is an error 74 75from vizinfo import colors access color list 76print colors list defined color names and RGB values 77colors["nickname"] = [R,G,B] set new RGB values from 0 to 255 78 79 140 pre-defined colors: red, green, blue, purple, yellow, black, white, etc 80 81Settings specific to gl tool: 82 83g.q(10) set quality of image (def = 5) 84g.axis(0/1) turn xyz axes off/on 85g.ortho(0/1) perspective (0) vs orthographic (1) view 86g.clip('xlo',0.25) clip in xyz from lo/hi at box fraction (0-1) 87g.reload() force all data to be reloaded 88g.cache = 0/1 turn off/on GL cache lists (def = on) 89theta,phi,x,y,scale,up = g.gview() grab all current view parameters 90g.sview(theta,phi,x,y,scale,up) set all view parameters 91 92 data reload is necessary if dump selection is used to change the data 93 cache lists usually improve graphics performance 94 gview returns values to use in other commands: 95 theta,phi are args to rotate() 96 x,y are args to shift() 97 scale is arg to zoom() 98 up is a 3-vector arg to sview() 99""" 100 101# History 102# 9/05, Steve Plimpton (SNL): original version 103 104# ToDo list 105# when do aselect with select str while looping N times on same timestep 106# would not let you grow # of atoms selected 107 108# Variables 109# ztheta = vertical angle from z-azis of viewpoint 110# azphi = azimuthal angle of viewpoint 111# xshift,yshift = xy translation of scene (in pixels) 112# distance = size of simulation box (largest dim) 113# eye = viewpoint distance from center of scene 114# file = filename prefix to use for images produced 115# boxflag = 0/1/2 for drawing simulation box: none/variable/fixed 116# bxcol = color of box 117# bxthick = thickness of box lines 118# bgcol = color of background 119# vizinfo = scene attributes 120# center[3] = center point of simulation box 121# view[3] = direction towards eye in simulation box (unit vector) 122# up[3] = screen up direction in simulation box (unit vector) 123# right[3] = screen right direction in simulation box (unit vector) 124 125# Imports and external programs 126 127from math import sin,cos,sqrt,pi,acos 128from OpenGL.Tk import * 129from OpenGL.GLUT import * 130import Image 131from vizinfo import vizinfo 132 133# Class definition 134 135class gl: 136 137# -------------------------------------------------------------------- 138 139 def __init__(self,data): 140 self.data = data 141 self.root = None 142 self.xpixels = 512 143 self.ypixels = 512 144 self.ztheta = 60 145 self.azphi = 30 146 self.scale = 1.0 147 self.xshift = self.yshift = 0 148 149 self.file = "image" 150 self.boxflag = 0 151 self.bxcol = [1,1,0] 152 self.bxthick = 0.3 153 self.bgcol = [0,0,0] 154 self.labels = [] 155 self.panflag = 0 156 self.select = "" 157 158 self.axisflag = 0 159 self.orthoflag = 1 160 self.nslices = 5 161 self.nstacks = 5 162 self.nsides = 10 163 self.theta_amplify = 2 164 self.shiny = 2 165 166 self.clipflag = 0 167 self.clipxlo = self.clipylo = self.clipzlo = 0.0 168 self.clipxhi = self.clipyhi = self.clipzhi = 1.0 169 170 self.nclist = 0 171 self.calllist = [0] # indexed by 1-Ntype, so start with 0 index 172 self.cache = 1 173 self.cachelist = 0 174 175 self.boxdraw = [] 176 self.atomdraw = [] 177 self.bonddraw = [] 178 self.tridraw = [] 179 self.linedraw = [] 180 181 self.ready = 0 182 self.create_window() 183 184 self.vizinfo = vizinfo() 185 self.adef() 186 self.bdef() 187 self.tdef() 188 self.ldef() 189 190 self.center = 3*[0] 191 self.view = 3*[0] 192 self.up = 3*[0] 193 self.right = 3*[0] 194 self.viewupright() 195 196 # -------------------------------------------------------------------- 197 198 def bg(self,color): 199 from vizinfo import colors 200 self.bgcol = [colors[color][0]/255.0,colors[color][1]/255.0, 201 colors[color][2]/255.0] 202 self.w.tkRedraw() 203 204 # -------------------------------------------------------------------- 205 206 def size(self,xnew,ynew=None): 207 self.xpixels = xnew 208 if not ynew: self.ypixels = self.xpixels 209 else: self.ypixels = ynew 210 self.create_window() 211 212 # -------------------------------------------------------------------- 213 214 def axis(self,value): 215 self.axisflag = value 216 self.cachelist = -self.cachelist 217 self.w.tkRedraw() 218 219 # -------------------------------------------------------------------- 220 221 def create_window(self): 222 if self.root: self.root.destroy() 223 224 from __main__ import tkroot 225 self.root = Toplevel(tkroot) 226 self.root.title('Pizza.py gl tool') 227 228 self.w = MyOpengl(self.root,width=self.xpixels,height=self.ypixels, 229 double=1,depth=1) 230 self.w.pack(expand=YES) 231# self.w.pack(expand=YES,fill=BOTH) 232 233 glViewport(0,0,self.xpixels,self.ypixels) 234 glEnable(GL_LIGHTING); 235 glEnable(GL_LIGHT0); 236 glEnable(GL_DEPTH_TEST); 237 glLightModeli(GL_LIGHT_MODEL_TWO_SIDE,GL_TRUE); 238 glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) 239 240 self.rtrack = self.xpixels 241 if self.ypixels > self.xpixels: self.rtrack = self.ypixels 242 243 self.w.redraw = self.redraw 244 self.w.parent = self 245 self.w.tkRedraw() 246 tkroot.update_idletasks() # force window to appear 247 248 # -------------------------------------------------------------------- 249 250 def clip(self,which,value): 251 if which == "xlo": 252 self.clipxlo = value 253 if value > self.clipxhi: self.clipxlo = self.clipxhi 254 elif which == "xhi": 255 self.clipxhi = value 256 if value < self.clipxlo: self.clipxhi = self.clipxlo 257 elif which == "ylo": 258 self.clipylo = value 259 if value > self.clipyhi: self.clipylo = self.clipyhi 260 elif which == "yhi": 261 self.clipyhi = value 262 if value < self.clipylo: self.clipyhi = self.clipylo 263 elif which == "zlo": 264 self.clipzlo = value 265 if value > self.clipzhi: self.clipzlo = self.clipzhi 266 elif which == "zhi": 267 self.clipzhi = value 268 if value < self.clipzlo: self.clipzhi = self.clipzlo 269 270 oldflag = self.clipflag 271 if self.clipxlo > 0 or self.clipylo > 0 or self.clipzlo > 0 or \ 272 self.clipxhi < 1 or self.clipyhi < 1 or self.clipzhi < 1: 273 self.clipflag = 1 274 else: self.clipflag = 0 275 276 if oldflag == 0 and self.clipflag == 0: return 277 self.cachelist = -self.cachelist 278 self.w.tkRedraw() 279 280 # -------------------------------------------------------------------- 281 282 def q(self,value): 283 self.nslices = value 284 self.nstacks = value 285 self.make_atom_calllist() 286 self.cachelist = -self.cachelist 287 self.w.tkRedraw() 288 289 # -------------------------------------------------------------------- 290 291 def ortho(self,value): 292 self.orthoflag = value 293 self.w.tkRedraw() 294 295 # -------------------------------------------------------------------- 296 # set unit vectors for view,up,right from ztheta,azphi 297 # assume +z in scene should be up on screen (unless looking down z-axis) 298 # right = up x view 299 300 def viewupright(self): 301 self.view[0] = cos(pi*self.azphi/180) * sin(pi*self.ztheta/180) 302 self.view[1] = sin(pi*self.azphi/180) * sin(pi*self.ztheta/180) 303 self.view[2] = cos(pi*self.ztheta/180) 304 305 if self.ztheta == 0.0: 306 self.up[0] = cos(pi*self.azphi/180) 307 self.up[1] = -sin(pi*self.azphi/180) 308 self.up[2] = 0.0 309 elif self.ztheta == 180.0: 310 self.up[0] = cos(pi*self.azphi/180) 311 self.up[1] = sin(pi*self.azphi/180) 312 self.up[2] = 0.0 313 else: 314 dot = self.view[2] # dot = (0,0,1) . view 315 self.up[0] = -dot*self.view[0] # up projected onto v = dot * v 316 self.up[1] = -dot*self.view[1] # up perp to v = up - dot * v 317 self.up[2] = 1.0 - dot*self.view[2] 318 319 self.up = vecnorm(self.up) 320 self.right = veccross(self.up,self.view) 321 322 # -------------------------------------------------------------------- 323 # reset ztheta,azphi and thus view,up.right 324 # called as function from Pizza.py 325 326 def rotate(self,ztheta,azphi): 327 self.ztheta = ztheta 328 self.azphi = azphi 329 self.viewupright() 330 self.setview() 331 self.w.tkRedraw() 332 333 # -------------------------------------------------------------------- 334 # return all view params to reproduce current display via sview() 335 336 def gview(self): 337 return self.ztheta,self.azphi,self.xshift,self.yshift,self.scale,self.up 338 339 # -------------------------------------------------------------------- 340 # set current view, called by user with full set of view params 341 # up is not settable via any other call, all other params are 342 343 def sview(self,ztheta,azphi,xshift,yshift,scale,up): 344 self.ztheta = ztheta 345 self.azphi = azphi 346 self.xshift = xshift 347 self.yshift = yshift 348 self.scale = scale 349 self.up[0] = up[0] 350 self.up[1] = up[1] 351 self.up[2] = up[2] 352 self.up = vecnorm(self.up) 353 self.view[0] = cos(pi*self.azphi/180) * sin(pi*self.ztheta/180) 354 self.view[1] = sin(pi*self.azphi/180) * sin(pi*self.ztheta/180) 355 self.view[2] = cos(pi*self.ztheta/180) 356 self.right = veccross(self.up,self.view) 357 self.setview() 358 self.w.tkRedraw() 359 360 # -------------------------------------------------------------------- 361 # rotation triggered by mouse trackball 362 # project old,new onto unit trackball surf 363 # rotate view,up around axis of rotation = old x new 364 # right = up x view 365 # reset ztheta,azphi from view 366 367 def mouse_rotate(self,xnew,ynew,xold,yold): 368 369 # change y pixels to measure from bottom of window instead of top 370 371 yold = self.ypixels - yold 372 ynew = self.ypixels - ynew 373 374 # vold = unit vector to (xold,yold) projected onto trackball 375 # vnew = unit vector to (xnew,ynew) projected onto trackball 376 # return (no rotation) if either projection point is outside rtrack 377 378 vold = [0,0,0] 379 vold[0] = xold - (0.5*self.xpixels + self.xshift) 380 vold[1] = yold - (0.5*self.ypixels + self.yshift) 381 vold[2] = self.rtrack*self.rtrack - vold[0]*vold[0] - vold[1]*vold[1] 382 if vold[2] < 0: return 383 vold[2] = sqrt(vold[2]) 384 vold = vecnorm(vold) 385 386 vnew = [0,0,0] 387 vnew[0] = xnew - (0.5*self.xpixels + self.xshift) 388 vnew[1] = ynew - (0.5*self.ypixels + self.yshift) 389 vnew[2] = self.rtrack*self.rtrack - vnew[0]*vnew[0] - vnew[1]*vnew[1] 390 if vnew[2] < 0: return 391 vnew[2] = sqrt(vnew[2]) 392 vnew = vecnorm(vnew) 393 394 # rot = trackball rotation axis in screen ref frame = vold x vnew 395 # theta = angle of rotation = sin(theta) for small theta 396 # axis = rotation axis in body ref frame described by right,up,view 397 398 rot = veccross(vold,vnew) 399 theta = sqrt(rot[0]*rot[0] + rot[1]*rot[1] + rot[2]*rot[2]) 400 theta *= self.theta_amplify 401 402 axis = [0,0,0] 403 axis[0] = rot[0]*self.right[0] + rot[1]*self.up[0] + rot[2]*self.view[0] 404 axis[1] = rot[0]*self.right[1] + rot[1]*self.up[1] + rot[2]*self.view[1] 405 axis[2] = rot[0]*self.right[2] + rot[1]*self.up[2] + rot[2]*self.view[2] 406 axis = vecnorm(axis) 407 408 # view is changed by (axis x view) scaled by theta 409 # up is changed by (axis x up) scaled by theta 410 # force up to be perp to view via up_perp = up - (up . view) view 411 # right = up x view 412 413 delta = veccross(axis,self.view) 414 self.view[0] -= theta*delta[0] 415 self.view[1] -= theta*delta[1] 416 self.view[2] -= theta*delta[2] 417 self.view = vecnorm(self.view) 418 419 delta = veccross(axis,self.up) 420 self.up[0] -= theta*delta[0] 421 self.up[1] -= theta*delta[1] 422 self.up[2] -= theta*delta[2] 423 424 dot = vecdot(self.up,self.view) 425 self.up[0] -= dot*self.view[0] 426 self.up[1] -= dot*self.view[1] 427 self.up[2] -= dot*self.view[2] 428 self.up = vecnorm(self.up) 429 430 self.right = veccross(self.up,self.view) 431 432 # convert new view to ztheta,azphi 433 434 self.ztheta = acos(self.view[2])/pi * 180.0 435 if (self.ztheta == 0.0): self.azphi = 0.0 436 else: self.azphi = acos(self.view[0]/sin(pi*self.ztheta/180.0))/pi * 180.0 437 if self.view[1] < 0: self.azphi = 360.0 - self.azphi 438 self.setview() 439 self.w.tkRedraw() 440 441 # -------------------------------------------------------------------- 442 443 def shift(self,x,y): 444 self.xshift = x; 445 self.yshift = y; 446 self.setview() 447 self.w.tkRedraw() 448 449 # -------------------------------------------------------------------- 450 451 def zoom(self,scale): 452 self.scale = scale 453 self.setview() 454 self.w.tkRedraw() 455 456 # -------------------------------------------------------------------- 457 # set view params needed by redraw 458 # input: center = center of box 459 # distance = size of scene (longest box length) 460 # scale = zoom factor (1.0 = no zoom) 461 # xshift,yshift = translation factor in pixels 462 # view = unit vector from center to viewpoint 463 # up = unit vector in up direction in scene 464 # right = unit vector in right direction in scene 465 # output: eye = distance to view scene from 466 # xto,yto,zto = point to look to 467 # xfrom,yfrom,zfrom = point to look from 468 469 def setview(self): 470 if not self.ready: return # no distance since no scene yet 471 472 self.eye = 3 * self.distance / self.scale 473 xfactor = 0.5*self.eye*self.xshift/self.xpixels 474 yfactor = 0.5*self.eye*self.yshift/self.ypixels 475 476 self.xto = self.center[0] - xfactor*self.right[0] - yfactor*self.up[0] 477 self.yto = self.center[1] - xfactor*self.right[1] - yfactor*self.up[1] 478 self.zto = self.center[2] - xfactor*self.right[2] - yfactor*self.up[2] 479 480 self.xfrom = self.xto + self.eye*self.view[0] 481 self.yfrom = self.yto + self.eye*self.view[1] 482 self.zfrom = self.zto + self.eye*self.view[2] 483 484 # -------------------------------------------------------------------- 485 # box attributes, also used for triangle lines 486 487 def box(self,*args): 488 self.boxflag = args[0] 489 if len(args) > 1: 490 from vizinfo import colors 491 self.bxcol = [colors[args[1]][0]/255.0,colors[args[1]][1]/255.0, 492 colors[args[1]][2]/255.0] 493 if len(args) > 2: self.bxthick = args[2] 494 self.cachelist = -self.cachelist 495 self.w.tkRedraw() 496 497 # -------------------------------------------------------------------- 498 # grab all selected snapshots from data object 499 # add GL-specific info to each bond 500 501 def reload(self): 502 print "Loading data into gl tool ..." 503 data = self.data 504 505 self.timeframes = [] 506 self.boxframes = [] 507 self.atomframes = [] 508 self.bondframes = [] 509 self.triframes = [] 510 self.lineframes = [] 511 512 box = [] 513 if self.boxflag == 2: box = data.maxbox() 514 515 flag = 0 516 while 1: 517 which,time,flag = data.iterator(flag) 518 if flag == -1: break 519 time,boxone,atoms,bonds,tris,lines = data.viz(which) 520 if self.boxflag < 2: box = boxone 521 if bonds: self.bonds_augment(bonds) 522 523 self.timeframes.append(time) 524 self.boxframes.append(box) 525 self.atomframes.append(atoms) 526 self.bondframes.append(bonds) 527 self.triframes.append(tris) 528 self.lineframes.append(lines) 529 530 print time, 531 sys.stdout.flush() 532 print 533 534 self.nframes = len(self.timeframes) 535 self.distance = compute_distance(self.boxframes[0]) 536 self.center = compute_center(self.boxframes[0]) 537 self.ready = 1 538 self.setview() 539 540 # -------------------------------------------------------------------- 541 542 def nolabel(self): 543 self.cachelist = -self.cachelist 544 self.labels = [] 545 546 # -------------------------------------------------------------------- 547 # show a single snapshot 548 # distance from snapshot box or max box for all selected steps 549 550 def show(self,ntime): 551 data = self.data 552 which = data.findtime(ntime) 553 time,box,atoms,bonds,tris,lines = data.viz(which) 554 if self.boxflag == 2: box = data.maxbox() 555 self.distance = compute_distance(box) 556 self.center = compute_center(box) 557 558 if bonds: self.bonds_augment(bonds) 559 560 self.boxdraw = box 561 self.atomdraw = atoms 562 self.bonddraw = bonds 563 self.tridraw = tris 564 self.linedraw = lines 565 566 self.ready = 1 567 self.setview() 568 self.cachelist = -self.cachelist 569 self.w.tkRedraw() 570 self.save() 571 572 # -------------------------------------------------------------------- 573 574 def pan(self,*list): 575 if len(list) == 0: self.panflag = 0 576 else: 577 self.panflag = 1 578 self.ztheta_start = list[0] 579 self.azphi_start = list[1] 580 self.scale_start = list[2] 581 self.ztheta_stop = list[3] 582 self.azphi_stop = list[4] 583 self.scale_stop = list[5] 584 585 # -------------------------------------------------------------------- 586 587 def all(self,*list): 588 data = self.data 589 if len(list) == 0: 590 nstart = 0 591 ncount = data.nselect 592 elif len(list) == 1: 593 nstart = list[0] 594 ncount = data.nselect 595 else: 596 ntime = list[0] 597 nstart = list[2] 598 ncount = list[1] 599 600 if self.boxflag == 2: box = data.maxbox() 601 602 # loop over all selected steps 603 # distance from 1st snapshot box or max box for all selected steps 604 # recompute box center on 1st step or if panning 605 606 if len(list) <= 1: 607 608 n = nstart 609 i = flag = 0 610 while 1: 611 which,time,flag = data.iterator(flag) 612 if flag == -1: break 613 614 fraction = float(i) / (ncount-1) 615 616 if self.select != "": 617 newstr = self.select % fraction 618 data.aselect.test(newstr,time) 619 time,boxone,atoms,bonds,tris,lines = data.viz(which) 620 621 if self.boxflag < 2: box = boxone 622 if n == nstart: self.distance = compute_distance(box) 623 624 if n < 10: file = self.file + "000" + str(n) 625 elif n < 100: file = self.file + "00" + str(n) 626 elif n < 1000: file = self.file + "0" + str(n) 627 else: file = self.file + str(n) 628 629 if self.panflag: 630 self.ztheta = self.ztheta_start + \ 631 fraction*(self.ztheta_stop - self.ztheta_start) 632 self.azphi = self.azphi_start + \ 633 fraction*(self.azphi_stop - self.azphi_start) 634 self.scale = self.scale_start + \ 635 fraction*(self.scale_stop - self.scale_start) 636 self.viewupright() 637 638 if n == nstart or self.panflag: self.center = compute_center(box) 639 640 if bonds: self.bonds_augment(bonds) 641 642 self.boxdraw = box 643 self.atomdraw = atoms 644 self.bonddraw = bonds 645 self.tridraw = tris 646 self.linedraw = lines 647 648 self.ready = 1 649 self.setview() 650 self.cachelist = -self.cachelist 651 self.w.tkRedraw() 652 self.save(file) 653 654 print time, 655 sys.stdout.flush() 656 i += 1 657 n += 1 658 659 # loop ncount times on same step 660 # distance from 1st snapshot box or max box for all selected steps 661 # recompute box center on 1st step or if panning 662 663 else: 664 665 which = data.findtime(ntime) 666 667 n = nstart 668 for i in range(ncount): 669 fraction = float(i) / (ncount-1) 670 671 if self.select != "": 672 newstr = self.select % fraction 673 data.aselect.test(newstr,ntime) 674 time,boxone,atoms,bonds,tris,lines = data.viz(which) 675 676 if self.boxflag < 2: box = boxone 677 if n == nstart: self.distance = compute_distance(box) 678 679 if n < 10: file = self.file + "000" + str(n) 680 elif n < 100: file = self.file + "00" + str(n) 681 elif n < 1000: file = self.file + "0" + str(n) 682 else: file = self.file + str(n) 683 684 if self.panflag: 685 self.ztheta = self.ztheta_start + \ 686 fraction*(self.ztheta_stop - self.ztheta_start) 687 self.azphi = self.azphi_start + \ 688 fraction*(self.azphi_stop - self.azphi_start) 689 self.scale = self.scale_start + \ 690 fraction*(self.scale_stop - self.scale_start) 691 self.viewupright() 692 693 if n == nstart or self.panflag: self.center = compute_center(box) 694 695 if bonds: self.bonds_augment(bonds) 696 697 self.boxdraw = box 698 self.atomdraw = atoms 699 self.bonddraw = bonds 700 self.tridraw = tris 701 self.linedraw = lines 702 703 self.ready = 1 704 self.setview() 705 self.cachelist = -self.cachelist 706 self.w.tkRedraw() 707 self.save(file) 708 709 print n, 710 sys.stdout.flush() 711 n += 1 712 713 print "\n%d images" % ncount 714 715 # -------------------------------------------------------------------- 716 717 def display(self,index): 718 self.boxdraw = self.boxframes[index] 719 self.atomdraw = self.atomframes[index] 720 self.bonddraw = self.bondframes[index] 721 self.tridraw = self.triframes[index] 722 self.linedraw = self.lineframes[index] 723 724 self.ready = 1 725 self.cachelist = -self.cachelist 726 self.w.tkRedraw() 727 return (self.timeframes[index],len(self.atomdraw)) 728 729 # -------------------------------------------------------------------- 730 # draw the GL scene 731 732 def redraw(self,o): 733 # clear window to background color 734 735 glClearColor(self.bgcol[0],self.bgcol[1],self.bgcol[2],0) 736 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) 737 738 # not ready if no scene yet 739 740 if not self.ready: return 741 742 # set view from eye, distance, 3 lookat vectors (from,to,up) 743 744 glMatrixMode(GL_PROJECTION) 745 glLoadIdentity() 746 if self.orthoflag: 747 glOrtho(-0.25*self.eye,0.25*self.eye,-0.25*self.eye,0.25*self.eye, 748 self.eye-2*self.distance,self.eye+2*self.distance) 749 else: 750 gluPerspective(30.0,1.0,0.01,10000.0) 751 752 glMatrixMode(GL_MODELVIEW) 753 glLoadIdentity() 754 gluLookAt(self.xfrom,self.yfrom,self.zfrom,self.xto,self.yto,self.zto, 755 self.up[0],self.up[1],self.up[2]) 756 757 # draw scene from display list if caching allowed and list hasn't changed 758 # else redraw and store as new display list if caching allowed 759 760 if self.cache and self.cachelist > 0: glCallList(self.cachelist); 761 else: 762 if self.cache: 763 if self.cachelist < 0: glDeleteLists(-self.cachelist,1) 764 self.cachelist = glGenLists(1) 765 glNewList(self.cachelist,GL_COMPILE_AND_EXECUTE) 766 767 # draw box, clip-box, xyz axes, lines 768 769 glDisable(GL_LIGHTING) 770 771 if self.boxflag: 772 self.draw_box(0) 773 if self.clipflag: self.draw_box(1) 774 if self.axisflag: self.draw_axes() 775 776 ncolor = self.vizinfo.nlcolor 777 for line in self.linedraw: 778 itype = int(line[1]) 779 if itype > ncolor: raise StandardError,"line type too big" 780 red,green,blue = self.vizinfo.lcolor[itype] 781 glColor3f(red,green,blue) 782 thick = self.vizinfo.lrad[itype] 783 glLineWidth(thick) 784 glBegin(GL_LINES) 785 glVertex3f(line[2],line[3],line[4]) 786 glVertex3f(line[5],line[6],line[7]) 787 glEnd() 788 789 glEnable(GL_LIGHTING) 790 791 # draw non-clipped scene = atoms, bonds, triangles 792 793# draw atoms as collection of points 794# cannot put PointSize inside glBegin 795# so probably need to group atoms by type for best performance 796# or just allow one radius 797# need to scale radius appropriately with box size 798# or could leave it at absolute value 799# use POINT_SMOOTH to enable anti-aliasing and round points 800# multiple timesteps via vcr::play() is still not fast 801# caching makes it fast for single frame, but multiple frames is slow 802# need to enable clipping 803 804# if not self.clipflag: 805# glDisable(GL_LIGHTING) 806# glEnable(GL_POINT_SMOOTH) 807# glPointSize(self.vizinfo.arad[int(self.atomdraw[0][1])]) 808# glBegin(GL_POINTS) 809# for atom in self.atomdraw: 810# red,green,blue = self.vizinfo.acolor[int(atom[1])] 811# glColor(red,green,blue) 812# glVertex3d(atom[2],atom[3],atom[4]) 813# glEnd() 814# glEnable(GL_LIGHTING) 815 816 if not self.clipflag: 817 for atom in self.atomdraw: 818 glTranslatef(atom[2],atom[3],atom[4]); 819 glCallList(self.calllist[int(atom[1])]); 820 glTranslatef(-atom[2],-atom[3],-atom[4]); 821 822 if self.bonddraw: 823 bound = 0.25 * self.distance 824 ncolor = self.vizinfo.nbcolor 825 for bond in self.bonddraw: 826 if bond[10] > bound: continue 827 itype = int(bond[1]) 828 if itype > ncolor: raise StandardError,"bond type too big" 829 red,green,blue = self.vizinfo.bcolor[itype] 830 rad = self.vizinfo.brad[itype] 831 glPushMatrix() 832 glTranslatef(bond[2],bond[3],bond[4]) 833 glRotatef(bond[11],bond[12],bond[13],0.0) 834 glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,[red,green,blue,1.0]); 835 glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); 836 obj = gluNewQuadric() 837 gluCylinder(obj,rad,rad,bond[10],self.nsides,self.nsides) 838 glPopMatrix() 839 840 if self.tridraw: 841 fillflag = self.vizinfo.tfill[int(self.tridraw[0][1])] 842 843 if fillflag != 1: 844 if fillflag: 845 glEnable(GL_POLYGON_OFFSET_FILL) 846 glPolygonOffset(1.0,1.0) 847 glBegin(GL_TRIANGLES) 848 ncolor = self.vizinfo.ntcolor 849 for tri in self.tridraw: 850 itype = int(tri[1]) 851 if itype > ncolor: raise StandardError,"tri type too big" 852 red,green,blue = self.vizinfo.tcolor[itype] 853 glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,[red,green,blue,1.0]); 854 glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); 855 glNormal3f(tri[11],tri[12],tri[13]) 856 glVertex3f(tri[2],tri[3],tri[4]) 857 glVertex3f(tri[5],tri[6],tri[7]) 858 glVertex3f(tri[8],tri[9],tri[10]) 859 glEnd() 860 if fillflag: glDisable(GL_POLYGON_OFFSET_FILL) 861 862 if fillflag: 863 glDisable(GL_LIGHTING) 864 glPolygonMode(GL_FRONT_AND_BACK,GL_LINE) 865 glLineWidth(self.bxthick) 866 glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) 867 glBegin(GL_TRIANGLES) 868 for tri in self.tridraw: 869 glVertex3f(tri[2],tri[3],tri[4]) 870 glVertex3f(tri[5],tri[6],tri[7]) 871 glVertex3f(tri[8],tri[9],tri[10]) 872 glEnd() 873 glEnable(GL_LIGHTING) 874 glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) 875 876 # draw clipped scene = atoms, bonds, triangles 877 878 else: 879 box = self.boxdraw 880 xlo = box[0] + self.clipxlo*(box[3] - box[0]) 881 xhi = box[0] + self.clipxhi*(box[3] - box[0]) 882 ylo = box[1] + self.clipylo*(box[4] - box[1]) 883 yhi = box[1] + self.clipyhi*(box[4] - box[1]) 884 zlo = box[2] + self.clipzlo*(box[5] - box[2]) 885 zhi = box[2] + self.clipzhi*(box[5] - box[2]) 886 887 for atom in self.atomdraw: 888 x,y,z = atom[2],atom[3],atom[4] 889 if x >= xlo and x <= xhi and y >= ylo and y <= yhi and \ 890 z >= zlo and z <= zhi: 891 glTranslatef(x,y,z); 892 glCallList(self.calllist[int(atom[1])]); 893 glTranslatef(-x,-y,-z); 894 895 if self.bonddraw: 896 bound = 0.25 * self.distance 897 ncolor = self.vizinfo.nbcolor 898 for bond in self.bonddraw: 899 xmin = min2(bond[2],bond[5]) 900 xmax = max2(bond[2],bond[5]) 901 ymin = min2(bond[3],bond[6]) 902 ymax = max2(bond[3],bond[6]) 903 zmin = min2(bond[4],bond[7]) 904 zmax = max2(bond[4],bond[7]) 905 if xmin >= xlo and xmax <= xhi and \ 906 ymin >= ylo and ymax <= yhi and zmin >= zlo and zmax <= zhi: 907 if bond[10] > bound: continue 908 itype = int(bond[1]) 909 if itype > ncolor: raise StandardError,"bond type too big" 910 red,green,blue = self.vizinfo.bcolor[itype] 911 rad = self.vizinfo.brad[itype] 912 glPushMatrix() 913 glTranslatef(bond[2],bond[3],bond[4]) 914 glRotatef(bond[11],bond[12],bond[13],0.0) 915 glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,[red,green,blue,1.0]); 916 glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); 917 obj = gluNewQuadric() 918 gluCylinder(obj,rad,rad,bond[10],self.nsides,self.nsides) 919 glPopMatrix() 920 921 if self.tridraw: 922 fillflag = self.vizinfo.tfill[int(self.tridraw[0][1])] 923 924 if fillflag != 1: 925 if fillflag: 926 glEnable(GL_POLYGON_OFFSET_FILL) 927 glPolygonOffset(1.0,1.0) 928 glBegin(GL_TRIANGLES) 929 ncolor = self.vizinfo.ntcolor 930 for tri in self.tridraw: 931 xmin = min3(tri[2],tri[5],tri[8]) 932 xmax = max3(tri[2],tri[5],tri[8]) 933 ymin = min3(tri[3],tri[6],tri[9]) 934 ymax = max3(tri[3],tri[6],tri[9]) 935 zmin = min3(tri[4],tri[7],tri[10]) 936 zmax = max3(tri[4],tri[7],tri[10]) 937 if xmin >= xlo and xmax <= xhi and \ 938 ymin >= ylo and ymax <= yhi and \ 939 zmin >= zlo and zmax <= zhi: 940 itype = int(tri[1]) 941 if itype > ncolor: raise StandardError,"tri type too big" 942 red,green,blue = self.vizinfo.tcolor[itype] 943 glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION, 944 [red,green,blue,1.0]); 945 glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,self.shiny); 946 glNormal3f(tri[11],tri[12],tri[13]) 947 glVertex3f(tri[2],tri[3],tri[4]) 948 glVertex3f(tri[5],tri[6],tri[7]) 949 glVertex3f(tri[8],tri[9],tri[10]) 950 glEnd() 951 if fillflag: glDisable(GL_POLYGON_OFFSET_FILL) 952 953 if fillflag: 954 glDisable(GL_LIGHTING) 955 glPolygonMode(GL_FRONT_AND_BACK,GL_LINE) 956 glLineWidth(self.bxthick) 957 glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) 958 glBegin(GL_TRIANGLES) 959 for tri in self.tridraw: 960 xmin = min3(tri[2],tri[5],tri[8]) 961 xmax = max3(tri[2],tri[5],tri[8]) 962 ymin = min3(tri[3],tri[6],tri[9]) 963 ymax = max3(tri[3],tri[6],tri[9]) 964 zmin = min3(tri[4],tri[7],tri[10]) 965 zmax = max3(tri[4],tri[7],tri[10]) 966 if xmin >= xlo and xmax <= xhi and \ 967 ymin >= ylo and ymax <= yhi and \ 968 zmin >= zlo and zmax <= zhi: 969 glVertex3f(tri[2],tri[3],tri[4]) 970 glVertex3f(tri[5],tri[6],tri[7]) 971 glVertex3f(tri[8],tri[9],tri[10]) 972 glEnd() 973 glEnable(GL_LIGHTING) 974 glPolygonMode(GL_FRONT_AND_BACK,GL_FILL) 975 976 if self.cache: glEndList() 977 978 glFlush() 979 980 # -------------------------------------------------------------------- 981 # make new call list for each atom type 982 # called when atom color/rad/quality is changed 983 984 def make_atom_calllist(self): 985 # extend calllist array if necessary 986 987 if self.vizinfo.nacolor > self.nclist: 988 for i in range(self.vizinfo.nacolor-self.nclist): self.calllist.append(0) 989 self.nclist = self.vizinfo.nacolor 990 991 # create new calllist for each atom type 992 993 for itype in xrange(1,self.vizinfo.nacolor+1): 994 if self.calllist[itype]: glDeleteLists(self.calllist[itype],1) 995 ilist = glGenLists(1) 996 self.calllist[itype] = ilist 997 glNewList(ilist,GL_COMPILE) 998 red,green,blue = self.vizinfo.acolor[itype] 999 rad = self.vizinfo.arad[itype] 1000 glColor3f(red,green,blue); 1001 1002# glPointSize(10.0*rad) 1003# glBegin(GL_POINTS) 1004# glVertex3f(0.0,0.0,0.0) 1005# glEnd() 1006 1007 glMaterialfv(GL_FRONT,GL_EMISSION,[red,green,blue,1.0]); 1008 glMaterialf(GL_FRONT,GL_SHININESS,self.shiny); 1009 glutSolidSphere(rad,self.nslices,self.nstacks) 1010 glEndList() 1011 1012 # -------------------------------------------------------------------- 1013 # augment bond info returned by viz() with info needed for GL draw 1014 # info = length, theta, -dy, dx for bond orientation 1015 1016 def bonds_augment(self,bonds): 1017 for bond in bonds: 1018 dx = bond[5] - bond[2] 1019 dy = bond[6] - bond[3] 1020 dz = bond[7] - bond[4] 1021 length = sqrt(dx*dx + dy*dy + dz*dz) 1022 dx /= length 1023 dy /= length 1024 dz /= length 1025 theta = acos(dz)*180.0/pi 1026 bond += [length,theta,-dy,dx] 1027 1028 # -------------------------------------------------------------------- 1029 1030 def draw_box(self,flag): 1031 xlo,ylo,zlo,xhi,yhi,zhi = self.boxdraw 1032 1033 if flag: 1034 tmp = xlo + self.clipxlo*(xhi - xlo) 1035 xhi = xlo + self.clipxhi*(xhi - xlo) 1036 xlo = tmp 1037 tmp = ylo + self.clipylo*(yhi - ylo) 1038 yhi = ylo + self.clipyhi*(yhi - ylo) 1039 ylo = tmp 1040 tmp = zlo + self.clipzlo*(zhi - zlo) 1041 zhi = zlo + self.clipzhi*(zhi - zlo) 1042 zlo = tmp 1043 1044 glLineWidth(self.bxthick) 1045 glColor3f(self.bxcol[0],self.bxcol[1],self.bxcol[2]) 1046 1047 glBegin(GL_LINE_LOOP) 1048 glVertex3f(xlo,ylo,zlo) 1049 glVertex3f(xhi,ylo,zlo) 1050 glVertex3f(xhi,yhi,zlo) 1051 glVertex3f(xlo,yhi,zlo) 1052 glEnd() 1053 1054 glBegin(GL_LINE_LOOP) 1055 glVertex3f(xlo,ylo,zhi) 1056 glVertex3f(xhi,ylo,zhi) 1057 glVertex3f(xhi,yhi,zhi) 1058 glVertex3f(xlo,yhi,zhi) 1059 glEnd() 1060 1061 glBegin(GL_LINES) 1062 glVertex3f(xlo,ylo,zlo) 1063 glVertex3f(xlo,ylo,zhi) 1064 glVertex3f(xhi,ylo,zlo) 1065 glVertex3f(xhi,ylo,zhi) 1066 glVertex3f(xhi,yhi,zlo) 1067 glVertex3f(xhi,yhi,zhi) 1068 glVertex3f(xlo,yhi,zlo) 1069 glVertex3f(xlo,yhi,zhi) 1070 glEnd() 1071 1072 # -------------------------------------------------------------------- 1073 1074 def draw_axes(self): 1075 xlo,ylo,zlo,xhi,yhi,zhi = self.boxdraw 1076 1077 delta = xhi-xlo 1078 if yhi-ylo > delta: delta = yhi-ylo 1079 if zhi-zlo > delta: delta = zhi-zlo 1080 delta *= 0.1 1081 1082 glLineWidth(self.bxthick) 1083 1084 glBegin(GL_LINES) 1085 glColor3f(1,0,0) 1086 glVertex3f(xlo-delta,ylo-delta,zlo-delta) 1087 glVertex3f(xhi-delta,ylo-delta,zlo-delta) 1088 glColor3f(0,1,0) 1089 glVertex3f(xlo-delta,ylo-delta,zlo-delta) 1090 glVertex3f(xlo-delta,yhi-delta,zlo-delta) 1091 glColor3f(0,0,1) 1092 glVertex3f(xlo-delta,ylo-delta,zlo-delta) 1093 glVertex3f(xlo-delta,ylo-delta,zhi-delta) 1094 glEnd() 1095 1096 # -------------------------------------------------------------------- 1097 1098 def save(self,file=None): 1099 self.w.update() # force image on screen to be current before saving it 1100 1101 pstring = glReadPixels(0,0,self.xpixels,self.ypixels, 1102 GL_RGBA,GL_UNSIGNED_BYTE) 1103 snapshot = Image.fromstring("RGBA",(self.xpixels,self.ypixels),pstring) 1104 snapshot = snapshot.transpose(Image.FLIP_TOP_BOTTOM) 1105 1106 if not file: file = self.file 1107 snapshot.save(file + ".png") 1108 1109 # -------------------------------------------------------------------- 1110 1111 def adef(self): 1112 self.vizinfo.setcolors("atom",range(100),"loop") 1113 self.vizinfo.setradii("atom",range(100),0.45) 1114 self.make_atom_calllist() 1115 self.cachelist = -self.cachelist 1116 self.w.tkRedraw() 1117 1118 # -------------------------------------------------------------------- 1119 1120 def bdef(self): 1121 self.vizinfo.setcolors("bond",range(100),"loop") 1122 self.vizinfo.setradii("bond",range(100),0.25) 1123 self.cachelist = -self.cachelist 1124 self.w.tkRedraw() 1125 1126 # -------------------------------------------------------------------- 1127 1128 def tdef(self): 1129 self.vizinfo.setcolors("tri",range(100),"loop") 1130 self.vizinfo.setfills("tri",range(100),0) 1131 self.cachelist = -self.cachelist 1132 self.w.tkRedraw() 1133 1134 # -------------------------------------------------------------------- 1135 1136 def ldef(self): 1137 self.vizinfo.setcolors("line",range(100),"loop") 1138 self.vizinfo.setradii("line",range(100),0.25) 1139 self.cachelist = -self.cachelist 1140 self.w.tkRedraw() 1141 1142 # -------------------------------------------------------------------- 1143 1144 def acol(self,atypes,colors): 1145 self.vizinfo.setcolors("atom",atypes,colors) 1146 self.make_atom_calllist() 1147 self.cachelist = -self.cachelist 1148 self.w.tkRedraw() 1149 1150 # -------------------------------------------------------------------- 1151 1152 def arad(self,atypes,radii): 1153 self.vizinfo.setradii("atom",atypes,radii) 1154 self.make_atom_calllist() 1155 self.cachelist = -self.cachelist 1156 self.w.tkRedraw() 1157 1158 # -------------------------------------------------------------------- 1159 1160 def bcol(self,btypes,colors): 1161 self.vizinfo.setcolors("bond",btypes,colors) 1162 self.cachelist = -self.cachelist 1163 self.w.tkRedraw() 1164 1165 # -------------------------------------------------------------------- 1166 1167 def brad(self,btypes,radii): 1168 self.vizinfo.setradii("bond",btypes,radii) 1169 self.cachelist = -self.cachelist 1170 self.w.tkRedraw() 1171 1172 # -------------------------------------------------------------------- 1173 1174 def tcol(self,ttypes,colors): 1175 self.vizinfo.setcolors("tri",ttypes,colors) 1176 self.cachelist = -self.cachelist 1177 self.w.tkRedraw() 1178 1179 # -------------------------------------------------------------------- 1180 1181 def tfill(self,ttypes,flags): 1182 self.vizinfo.setfills("tri",ttypes,flags) 1183 self.cachelist = -self.cachelist 1184 self.w.tkRedraw() 1185 1186 # -------------------------------------------------------------------- 1187 1188 def lcol(self,ltypes,colors): 1189 self.vizinfo.setcolors("line",ltypes,colors) 1190 self.cachelist = -self.cachelist 1191 self.w.tkRedraw() 1192 1193 # -------------------------------------------------------------------- 1194 1195 def lrad(self,ltypes,radii): 1196 self.vizinfo.setradii("line",ltypes,radii) 1197 self.cachelist = -self.cachelist 1198 self.w.tkRedraw() 1199 1200# -------------------------------------------------------------------- 1201# derived class from Togl's Opengl 1202# overwrite redraw, translate, rotate, scale methods 1203# latter 3 are mouse-motion methods 1204 1205class MyOpengl(Opengl): 1206 def __init__(self, master, cnf={}, **kw): 1207 args = (self,master,cnf) 1208 Opengl.__init__(*args,**kw) 1209 Opengl.autospin_allowed = 0 1210 1211 # redraw Opengl scene 1212 # call parent redraw() method 1213 1214 def tkRedraw(self,*dummy): 1215 if not self.initialised: return 1216 self.tk.call(self._w,'makecurrent') 1217 self.redraw(self) 1218 self.tk.call(self._w,'swapbuffers') 1219 1220 # left button translate 1221 # access parent xshift/yshift and call parent trans() method 1222 1223 def tkTranslate(self,event): 1224 dx = event.x - self.xmouse 1225 dy = event.y - self.ymouse 1226 x = self.parent.xshift + dx 1227 y = self.parent.yshift - dy 1228 self.parent.shift(x,y) 1229 self.tkRedraw() 1230 self.tkRecordMouse(event) 1231 1232 # middle button trackball 1233 # call parent mouse_rotate() method 1234 1235 def tkRotate(self,event): 1236 self.parent.mouse_rotate(event.x,event.y,self.xmouse,self.ymouse) 1237 self.tkRedraw() 1238 self.tkRecordMouse(event) 1239 1240 # right button zoom 1241 # access parent scale and call parent zoom() method 1242 1243 def tkScale(self,event): 1244 scale = 1 - 0.01 * (event.y - self.ymouse) 1245 if scale < 0.001: scale = 0.001 1246 elif scale > 1000: scale = 1000 1247 scale *= self.parent.scale 1248 self.parent.zoom(scale) 1249 self.tkRedraw() 1250 self.tkRecordMouse(event) 1251 1252# -------------------------------------------------------------------- 1253# draw a line segment 1254 1255def segment(p1,p2): 1256 glVertex3f(p1[0],p1[1],p1[2]) 1257 glVertex3f(p2[0],p2[1],p2[2]) 1258 1259# -------------------------------------------------------------------- 1260# normalize a 3-vector to unit length 1261 1262def vecnorm(v): 1263 length = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]) 1264 return [v[0]/length,v[1]/length,v[2]/length] 1265 1266# -------------------------------------------------------------------- 1267# dot product of two 3-vectors 1268 1269def vecdot(v1,v2): 1270 return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] 1271 1272# -------------------------------------------------------------------- 1273# cross product of two 3-vectors 1274 1275def veccross(v1,v2): 1276 v = [0,0,0] 1277 v[0] = v1[1]*v2[2] - v1[2]*v2[1] 1278 v[1] = v1[2]*v2[0] - v1[0]*v2[2] 1279 v[2] = v1[0]*v2[1] - v1[1]*v2[0] 1280 return v 1281 1282# -------------------------------------------------------------------- 1283# return characteristic distance of simulation domain = max dimension 1284 1285def compute_distance(box): 1286 distance = box[3]-box[0] 1287 if box[4]-box[1] > distance: distance = box[4]-box[1] 1288 if box[5]-box[2] > distance: distance = box[5]-box[2] 1289 return distance 1290 1291# -------------------------------------------------------------------- 1292# return center of box as 3 vector 1293 1294def compute_center(box): 1295 c = [0,0,0] 1296 c[0] = 0.5 * (box[0] + box[3]) 1297 c[1] = 0.5 * (box[1] + box[4]) 1298 c[2] = 0.5 * (box[2] + box[5]) 1299 return c 1300 1301# -------------------------------------------------------------------- 1302# return min of 2 values 1303 1304def min2(a,b): 1305 if b < a: a = b 1306 return a 1307 1308# -------------------------------------------------------------------- 1309# return max of 2 values 1310 1311def max2(a,b): 1312 if b > a: a = b 1313 return a 1314 1315# -------------------------------------------------------------------- 1316# return min of 3 values 1317 1318def min3(a,b,c): 1319 if b < a: a = b 1320 if c < a: a = c 1321 return a 1322 1323# -------------------------------------------------------------------- 1324# return max of 3 values 1325 1326def max3(a,b,c): 1327 if b > a: a = b 1328 if c > a: a = c 1329 return a 1330