# Copyright 2004-2007 Nanorex, Inc. See LICENSE file for details. """ geometry.py -- miscellaneous purely geometric routines. @author: Josh @version: $Id$ @copyright: 2004-2007 Nanorex, Inc. See LICENSE file for details. History: - Made by bruce 060119 from code (mostly by Josh) split out of other files. """ # Note: this module should not import from model modules # (such as those defining class Atom or Chunk). Callers # should handle model or ui dependent features themselves, # e.g. extract geometric info from model objects before # passing it here, or generate history warnings, or make up # default axes based on the screen direction. The code in # this file should remain purely geometric. import math from Numeric import transpose, minimum, maximum, remainder, size, add from Numeric import Float, zeros, multiply, sign, dot from LinearAlgebra import solve_linear_equations, eigenvectors from utilities import debug_flags # for atom_debug from geometry.VQT import V, A, cat, norm, cross, X_AXIS, Y_AXIS def selection_polyhedron(basepos, borderwidth = 1.8): """ Given basepos (a Numeric array of 3d positions), compute and return (as a list of vertices, each pair being an edge to draw) a simple bounding polyhedron, convenient for designating the approximate extent of this set of points. (This is used on the array of atom and open bond positions in a chunk to designate a selected chunk.) Make every face farther out than it needs to be to enclose all points in basepos, by borderwidth (default 1.8). Negative values of borderwidth might sometimes work, but are likely to fail drastically if the polyhedron is too small. The orientation of the faces matches the coordinate axes used in basepos (in typical calls, these are chunk-relative). [#e Someday we might permit caller-specified orientation axes. The axes of inertia would be interesting to try.] """ #bruce 060226/060605 made borderwidth an argument (non-default values untested); documented retval format #bruce 060119 split this out of shakedown_poly_evals_evecs_axis() in chunk.py. # Note, it has always had some bad bugs for certain cases, like long diagonal rods. if not len(basepos): return [] # a guess # find extrema in many directions xtab = dot(basepos, polyXmat) mins = minimum.reduce(xtab) - borderwidth maxs = maximum.reduce(xtab) + borderwidth polyhedron = makePolyList(cat(maxs,mins)) # apparently polyhedron is just a long list of vertices [v0,v1,v2,v3,...] to be passed to GL_LINES # (via drawlinelist), which will divide them into pairs and draw lines v0-v1, v2-v3, etc. # Maybe we should generalize the format, so some callers could also draw faces. # [bruce 060226/060605 comment] return polyhedron # == helper definitions for selection_polyhedron [moved here from drawer.py by bruce 060119] # mat mult by rowvector list and max/min reduce to find extrema D = math.sqrt(2.0)/2.0 T = math.sqrt(3.0)/3.0 # 0 1 2 3 4 5 6 7 8 9 10 11 12 # 13 14 15 16 17 18 19 20 21 22 23 24 25 polyXmat = A([[1, 0, 0, D, D, D, D, 0, 0, T, T, T, T], [0, 1, 0, D, -D, 0, 0, D, D, T, T, -T, -T], [0, 0, 1, 0, 0, D, -D, D, -D, T, -T, T, -T]]) del D, T polyMat = cat(transpose(polyXmat),transpose(polyXmat)) polyTab = [( 9, (0,7,3), [3,0,5,2,7,1,3]), (10, (0,1,4), [3,1,8,15,6,0,3]),#(10, (0,4,1), [3,0,6,15,8,1,3]), (11, (8,11,7), [4,14,21,2,5,0,4]), (12, (8,4,9), [4,0,6,15,20,14,4]), (22, (5,10,9), [18,13,16,14,20,15,18]), (23, (10,6,11), [16,13,19,2,21,14,16]), (24, (1,2,5), [8,1,17,13,18,15,8]), (25, (3,6,2), [7,2,19,13,17,1,7])] #( 9, (0,7,3), [3,0,5,2,7,1,3]), #(10, (0,1,4), [3,1,8,15,6,0,3]), #(11, (8,11,7), [4,14,21,2,5,0,4]), #(12, (8,4,9), [4,0,6,15,20,14,4]), #(22, (5,10,9), [18,13,16,14,20,15,18]), #(23, (10,6,11), [16,13,19,2,21,14,16]), #(24, (1,2,5), [8,1,17,13,18,15,8]), #(25, (3,6,2), [7,2,19,13,17,1,7])] ##(22, (10,5,9), [16,13,18,15,20,14,16]), #(23, (10,6,11), [16,13,19,2,21,14,16]), #(24, (2, 5, 1), [17,13,18,15,8,1,17]), #(25, (2,3,6), [17,1,7,2,19,13,17])] def planepoint(v,i,j,k): mat = V(polyMat[i],polyMat[j],polyMat[k]) vec = V(v[i],v[j],v[k]) return solve_linear_equations(mat, vec) def makePolyList(v): xlines = [[],[],[],[],[],[],[],[],[],[],[],[]] segs = [] for corner, edges, planes in polyTab: linx = [] for i in range(3): l,s,r = planes[2*i:2*i+3] e = remainder(i+1,3) p1 = planepoint(v,corner,l,r) if abs(dot(p1,polyMat[s])) <= abs(v[s]): p2 = planepoint(v,l,s,r) linx += [p1] xlines[edges[i]] += [p2] xlines[edges[e]] += [p2] segs += [p1,p2] else: p1 = planepoint(v,corner,l,s) p2 = planepoint(v,corner,r,s) linx += [p1,p2] xlines[edges[i]] += [p1] xlines[edges[e]] += [p2] e=edges[0] xlines[e] = xlines[e][:-2] + [xlines[e][-1],xlines[e][-2]] for p1,p2 in zip(linx, linx[1:]+[linx[0]]): segs += [p1,p2] ctl = 12 for lis in xlines[:ctl]: segs += [lis[0],lis[3],lis[1],lis[2]] assert type(segs) == type([]) #bruce 041119 return segs def trialMakePolyList(v): # [i think this is experimental code by Huaicai, never called, perhaps never tested. -- bruce 051117] pMat = [] for i in range(size(v)): pMat += [polyMat[i] * v[i]] segs = [] for corner, edges, planes in polyTab: pts = size(planes) for i in range(pts): segs += [pMat[corner], pMat[planes[i]]] segs += [pMat[planes[i]], pMat[planes[(i+1)%pts]]] return segs # == def inertia_eigenvectors(basepos, already_centered = False): """ Given basepos (an array of positions), compute and return (as a 2-tuple) the lists of eigenvalues and eigenvectors of the inertia tensor (computed as if all points had the same mass). These lists are always length 3, even for len(basepos) of 0,1, or 2, overlapping or colinear points, etc, though some evals will be 0 in these cases. Optional small speedup: if caller knows basepos is centered at the origin, it can say so. """ #bruce 060119 split this out of shakedown_poly_evals_evecs_axis() in chunk.py basepos = A(basepos) # make sure it's a Numeric array if not already_centered and len(basepos): center = add.reduce(basepos)/len(basepos) basepos = basepos - center # compute inertia tensor tensor = zeros((3,3),Float) for p in basepos: rsq = dot(p, p) m= - multiply.outer(p, p) m[0,0] += rsq m[1,1] += rsq m[2,2] += rsq tensor += m evals, evecs = eigenvectors(tensor) assert len(evals) == len(evecs) == 3 return evals, evecs if 0: # self-test; works fine as of 060119 def testie(points): print "ie( %r ) is %r" % ( points, inertia_eigenvectors(points) ) map( testie, [ [], [V(0,0,0)], [V(0,0,1)], # not centered [V(0,0,-1), V(0,0,1)], [V(0,0,-1), V(0,0,0), V(0,0,1)], [V(0,1,1), V(0,-1,-1), V(0,1,-1), V(0,-1,1)], [V(1,4,9), V(1,4,-9), V(1,-4,9), V(1,-4,-9), V(-1,4,9), V(-1,4,-9), V(-1,-4,9), V(-1,-4,-9)], ]) # == def unzip(pairs): """ inverse of zip, for a list of pairs. [#e should generalize.] [#e probably there's a simple general implem - transpose?] """ return [pair[0] for pair in pairs], [pair[1] for pair in pairs] def compute_heuristic_axis( basepos, type, evals_evecs = None, already_centered = False, aspect_threshhold = 0.95, numeric_threshhold = 0.0001, near1 = None, near2 = None, nears = (), dflt = None ): #bruce 060120 adding nears; will doc this and let it fully replace near1 and near2, but not yet, # since Mark is right now adding code that calls this with near1 and near2 ###@@@ """ Given basepos (an array of positions), compute and return an axis in one of various ways according to 'type' (choices are listed below), optionally adjusting the algorithm using aspect_threshhold (when are two dimensions close enough to count as circular), and numeric_threshhold (roughly, precision of coordinate values) (numeric_threshhold might be partly nim ###). Optional speed optimizations: if you know basepos is already_centered, you can say so; if (even better) you have the value of inertia_eigenvectors(basepos, already_centered = already_centered), you can pass that value so it doesn't have to be recomputed. (Computing it is the only way basepos is used in this routine.) The answer is always arbitrary in sign, so try to disambiguate it using (in order) whichever of near1, near2, and dflt is provided (by making it have positive dot product with the first one possible); but if none of those are provided and helpful (not perpendicular), return the answer with an arbitrary sign. If the positions require an arbitrary choice within a plane, and near1 is provided, try to choose the direction in that plane closest to near1; if they're all equally close (using numeric_threshhold), try near2 in the same way (ideally near2 should be perpendicular to near1). If the positions require a completely arbitrary choice, return dflt. Typical calls from the UI pass near1 = out of screen, near2 = up (parallel to screen), dflt = out of screen or None (depending on whether caller wants to detect the need to use dflt in order to print a warning). [#k guesses] Typical calls with no UI pass near1 = V(0,0,1), near2 = V(0,1,0), dflt = V(0,0,1). Choices of type (required argument) are: 'normal' - try to return the shortest axis; if ambiguous, use near directions to choose within a plane, or return dflt for a blob. Works for 2 or more points; treats single point (or no points) as a blob. This should be used by rotary and linear motors, since they want to connect to atoms on a surface (intention as of 060119), and by viewNormalTo. 'parallel' - try to return the longest axis; otherwise like 'normal'. This should be used by viewParallelTo. 'chunk' - for slabs close enough to circles or squares (using aspect_threshhold), like 'normal', otherwise like 'parallel'. This is used for computing axes of chunks for purposes of how they interactively slide in certain UI modes. (Note that viewParallelTo and viewNormalTo, applied to a chunk or perhaps to jigs containing sets of atoms, should compute their own values as if applied to the same atoms, since the chunk or jig will compute its axis in a different way, both in choice of type argument, and in values of near1, near2, and dflt.) """ #bruce 060119 split this out of shakedown_poly_evals_evecs_axis() in chunk.py, added some options. # Could clean it up by requiring caller to pass evals_evecs (no longer needing args basepos and already_centered). # According to tested behavior of inertia_eigenvectors, we don't need any special cases for # len(basepos) in [0,1,2], or points overlapping or colinear, etc! if evals_evecs is None: evals_evecs = inertia_eigenvectors( basepos, already_centered = already_centered) del basepos evals, evecs = evals_evecs # evals[i] tells you moment of inertia when rotating around axis evecs[i], which is larger if points are farther from axis. # So for a slab of dims 1x4x9, for axis '1' we see dims '4' and '9' and eval is the highest in that case. valvecs = zip(evals, evecs) valvecs.sort() evals, evecs = unzip(valvecs) # (custom helper function, inverse of zip in this case) assert evals[0] <= evals[1] <= evals[2] # evals[0] is now the lowest-valued evals[i] (i.e. longest axis, heuristically), evals[2] the highest. # (The heuristic relating this to axis-length assumes a uniform density of points # over some volume or its smooth-enough surface.) if type == 'chunk': if aspect_too_close( evals[0], evals[1], aspect_threshhold ): # sufficiently circular-or-square pancake/disk/slab/plane type = 'normal' # return the axle of a wheel else: # anything else type = 'parallel' # return the long axis pass if type == 'normal': # we want the shortest axis == largest eval; try axis (evec) with largest eval first: evals.reverse() evecs.reverse() elif type == 'parallel': # we want longest axis == shortest eval; order is already correct pass else: assert 0, "unrecognized type %r" % (type,) axes = [ evecs[0] ] for i in [1,2]: # order of these matters, I think if aspect_too_close( evals[0], evals[i], aspect_threshhold ): axes.append( evecs[i] ) del evals, evecs #bruce 060120 adding nears; will replace near1, near2, but for now, accept them all, use in that order: goodvecs = [near1, near2] + list(nears) # len(axes) now tells us how ambiguous the answer is, and axes are the vectors to make it from. if len(axes) == 1: # we know it up to a sign. answer = best_sign_on_vector( axes[0], goodvecs + [dflt], numeric_threshhold ) elif len(axes) == 2: # the answer is arbitrary within a plane. answer = best_vector_in_plane( axes, goodvecs, numeric_threshhold ) # (this could be None if caller didn't provide orthogonal near1 and near2) else: assert len(axes) == 3 # the answer is completely arbitrary answer = dflt return answer # end of compute_heuristic_axis # == helper functions for compute_heuristic_axis (likely to also be generally useful) def aspect_too_close( dim1, dim2, aspect_threshhold ): """ Are dim1 and dim2 (positive or zero real numbers) as close to 1:1 ratio as aspect_threshhold is to 1.0? """ # make sure it doesn't matter whether aspect_threshhold or 1/aspect_threshhold is passed aspect_threshhold = float(aspect_threshhold) if aspect_threshhold > 1: aspect_threshhold = 1.0 / aspect_threshhold if dim2 < dim1: dim2, dim1 = dim1, dim2 return dim1 >= (dim2 * aspect_threshhold) def best_sign_on_vector(vec, goodvecs, numeric_threshhold): """ vec is an arbitrary vector, and goodvecs is a list of unit vectors or (ignored) zero vectors or Nones; return vec or - vec, whichever is more in the same direction as the first goodvec which helps determine this (i.e. which is not zero or None, and is not perpendicular to vec, using numeric_threshhold to determine what's too close to call). If none of the goodvecs help, just return vec (this always happens if vec itself is too small, so it's safest if vec is a unit vector, though it's not required). """ for good in goodvecs: if good is not None: d = dot(vec, good) s = sign_with_threshhold(d, numeric_threshhold) if s: # it helps! return s * vec return vec def sign_with_threshhold( num, thresh ): """ Return -1, 0, or 1 as num is << 0, close to 0, or >> 0, where "close to 0" means abs(num) <= thresh. """ if abs(num) <= thresh: return 0 return sign(num) def best_vector_in_plane( axes, goodvecs, numeric_threshhold ): """ axes is a list of two orthonormal vectors defining a plane, and goodvecs is a list of unit vectors or (ignored) zero vectors or Nones; return whichever unit vector in the plane defined by axes is closest in direction to the first goodvec which helps determine this (i.e. which is not zero or None, and is not perpendicular to the plane, using numeric_threshhold to determine what's too close to call). If none of the goodvecs help, return None. """ x,y = axes for good in goodvecs: if good is not None: dx = dot(x,good) # will be 0 for good = V(0,0,0); not sensible for vlen(good) other than 0 or 1 dy = dot(y,good) if abs(dx) < numeric_threshhold and abs(dy) < numeric_threshhold: continue # good is perpendicular to the plane (or is zero) return norm(dx * x + dy * y) return None def arbitrary_perpendicular( vec, nicevecs = [] ): #bruce 060608, probably duplicates other code somewhere (check in VQT?) """ Return an arbitrary unit vector perpendicular to vec (a Numeric array, 3 floats), making it from vec and as-early-as-possible nicevecs (if any are passed). """ nicevecs = map(norm, nicevecs) + [X_AXIS, Y_AXIS] # we'll look at vec and each nicevec until they span a subspace which includes a perpendicular. # if we're not done, it means our subspace so far is spanned by just vec itself. vec = norm(vec) if not vec: return nicevecs[0] # the best we can do for nice in nicevecs: res = norm( nice - vec * dot(nice, vec) ) if res: return res return nicevecs[0] # should never happen def matrix_putting_axis_at_z(axis): #bruce 060608 """ Return an orthonormal matrix which can be used (via dot(points, matrix) for a Numeric array of 3d points) to transform points so that axis transforms to the z axis. (Not sure if it has determinant 1 or -1. If you're sure, change it to always be determinant 1.) """ # this code was modified from ops_select.py's findAtomUnderMouse, as it prepares to call Chunk.findAtomUnderMouse z = norm(axis) x = arbitrary_perpendicular(axis) y = cross(z,x) matrix = transpose(V(x,y,z)) return matrix """ notes, 060119, probably can be removed soon: motor axis - use 'view normal to' algorithm chunk axis - weirdest case - ### - for 1x4x9 take 9, but for 1x9x9 take 1. view normal to - for 1x4x9, take the short '1' dim; if answer ambiguous w/in a plane (atoms are in a sausage; aspect_threshhold?), prefer axis w/in that plane near the current one view parallel to - for 1x4x9, take the long '9' dim; if answer is ambiguous w/in a plane (atom plane is circular) (use aspect_threshhold to decide it's ambiguous) prefer an axis w/in that plane near the current one those work fine for 2 atoms (in different posns), as if 1x1x9 slab; for single atom or blob or 1x1x1 slab, in those cases, just return None (let caller do noop), or could return exact current axis... so the options are: - normal, parallel, or chunk - aspect_threshhold 0.95 - numeric_threshhold 0.0001 """ #end