LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
dodecahedron.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 import math
3 import itertools
4 import numpy
5 
6 class Dodecahedron(object):
7  """A dodecahedron
8 
9  Contains positions of faces and associated vertices
10  """
11  def __init__(self, withFacesOnPoles=False):
12  """Construct a Dodecahedron
13 
14  @param[in] withFacesOnPoles: if True center a face on each pole, else put a vertex on each pole
15  """
16  self._withFacesOnPoles = bool(withFacesOnPoles)
17 
18  # Basis cartesian vectors describing the faces of a dodecahedron; the full set of vectors is obtained
19  # by choosing both signs of each nonzero component of each vector.
20  # The orientation of the resulting dodecahedron, while very convenient
21  # for specifying face vectors, is not an orientation we want so it must be rotated.
22  g = (1.0 + math.sqrt(5.0)) / 2.0
23  faceBases = (
24  (0, 1, g),
25  (1, g, 0),
26  (g, 0, 1),
27  )
28  unrotFaceVecList = _computeFullVecList(faceBases)
29  unrotVertexVecList = _computeDodecahedronVertices(unrotFaceVecList)
30 
31  if self._withFacesOnPoles:
32  # one face is centered on each pole
33  vec0, vec1 = _findClosePair(unrotFaceVecList, 0)
34  rotMat = _computeCoordTransform(vec0, vec1)
35  else:
36  # one vertex is on each pole
37  vec0, vec1 = _findClosePair(unrotVertexVecList, 0)
38  rotMat = _computeCoordTransform(vec0, vec1, vec1NegativeX=True)
39  self.vertexVecList = [numpy.dot(rotMat, unrotVertexVec) for unrotVertexVec in unrotVertexVecList]
40  unsortedFaceList = [numpy.dot(rotMat, unrotFaceVec) for unrotFaceVec in unrotFaceVecList]
41  self.faceVecList = _sortedVectorList(unsortedFaceList)
42 
43  def getFaceCtrList(self):
44  """Return a list of face centers
45 
46  @return a list of face centers (in index order); each a unit vector (numpy array)
47  """
48  return self.faceVecList[:]
49 
50  def getFaceCtr(self, ind):
51  """Return the center of the specified face
52 
53  @param[in] ind: face index
54  @return face center as a unit vector (numpy array)
55  """
56  return self.faceVecList[ind][:]
57 
58  def getVertices(self, ind):
59  """Return the vertices for a given face
60 
61  @param[in] ind: face index
62  @return a list of vertices, each a unit vector (numpy array)
63  """
64  faceVec = self.getFaceCtr(ind)
65  vertexList, indList = _findCloseList(self.vertexVecList, faceVec)
66 
67  # sort vertex list about face vector (direction is random)
68  sortedVertexList = [vertexList[0]]
69  vertexList = list(vertexList[1:])
70  while len(vertexList) != 0:
71  nearVertexList, nearInd = _findCloseList(vertexList, sortedVertexList[-1])
72  sortedVertexList.append(nearVertexList[0])
73  vertexList.pop(nearInd[0])
74  return sortedVertexList
75 
76  def getFaceInd(self, vec):
77  """Return the index of the face containing the cartesian vector
78 
79  @param[in] vec: cartesian vector (length is ignored)
80  @return index of face containing vec
81  """
82  return numpy.argmax(numpy.dot(self.faceVecList, vec))
83 
85  """Get withFacesOnPoles parameter
86  """
87  return self._withFacesOnPoles
88 
89 def computeRotationMatrix(angle, axis):
90  """Return a 3D rotation matrix for rotation by a specified amount around a specified axis
91 
92  Inputs:
93  - angle: amount of rotation (rad)
94  - axis: axis of rotation; one of 0, 1 or 2 for x, y or z
95  """
96  cosAng = math.cos(angle)
97  sinAng = math.sin(angle)
98  rotMat = numpy.zeros((3,3), dtype=float)
99  rotMat[axis, axis] = 1
100  rotMat[(axis + 1) % 3, (axis + 1) % 3] = cosAng
101  rotMat[(axis + 2) % 3, (axis + 1) % 3] = sinAng
102  rotMat[(axis + 1) % 3, (axis + 2) % 3] = -sinAng
103  rotMat[(axis + 2) % 3, (axis + 2) % 3] = cosAng
104  return rotMat
105 
106 def _computeCoordTransform(vec0, vec1, vec1NegativeX=False):
107  """Compute a rotation matrix that puts vec0 along z and vec1 along +x in the xz plane
108 
109  Inputs:
110  - vec0: vector 0
111  - vec1: vector 1
112  - vec1NegativeX: if True then vec1 is rotated to face negative x
113  """
114  # rotate around x by angle of vec0 from z to y
115  xAng = math.atan2(vec0[1], vec0[2])
116  xRotMat = computeRotationMatrix(xAng, 0)
117 
118  # rotate around y by -angle of rotated vec0 from z to x
119  vec0RotX = numpy.dot(xRotMat, vec0)
120  yAng = -math.atan2(vec0RotX[0], vec0RotX[2])
121  yRotMat = computeRotationMatrix(yAng, 1)
122  xyRotMat = numpy.dot(yRotMat, xRotMat)
123 
124  # rotate around z by -angle of rotated vec1 from +/-x to y
125  vec1RotXY = numpy.dot(xyRotMat, vec1)
126  xVal = vec1RotXY[0]
127  if vec1NegativeX:
128  xVal = -xVal
129  zAng = -math.atan2(vec1RotXY[1], xVal)
130  zRotMat = computeRotationMatrix(zAng, 2)
131  xyzRotMat = numpy.dot(zRotMat, xyRotMat)
132  return xyzRotMat
133 
135  """Given a vector of face positions of a Dodecahedron compute the vertices
136  """
137  closeIndSetList = []
138  vertexDict = {}
139  for i in range(len(faceVecList)):
140  closeIndSet = _findCloseIndexSet(faceVecList, i)
141  if len(closeIndSet) != 5:
142  raise RuntimeError("Found %s vertices instead of 5 near %s: %s" % \
143  (len(closeIndSet), faceVecList[i], closeIndSet))
144  closeIndSetList.append(closeIndSet)
145  for i, iCloseIndSet in enumerate(closeIndSetList):
146  for j in iCloseIndSet:
147  jCloseIndSet = closeIndSetList[j]
148  sharedCloseIndSet = iCloseIndSet.intersection(jCloseIndSet)
149  if len(sharedCloseIndSet) != 2:
150  raise RuntimeError("Found %s vertices instead of 2 near %s and %s: %s" % \
151  (len(sharedCloseIndSet), faceVecList[i], faceVecList[j], sharedCloseIndSet))
152  for k in sharedCloseIndSet:
153  key = frozenset((i, j, k))
154  if key in vertexDict:
155  continue
156  vertexVec = faceVecList[i] + faceVecList[j] + faceVecList[k]
157  vertexVec /= numpy.sqrt(numpy.sum(vertexVec**2))
158  vertexDict[key] = vertexVec
159  return vertexDict.values()
160 
161 def _computeFullVecList(basisSet):
162  """Given a collection of basis vectors, compute all permutations with both signs of all nonzero values
163 
164  For example: [(0, 1, 2)] -> [(0, 1, 2), (0, -1, 2), (0, 1, -2), (0, -1, -2)]
165  """
166  fullSet = []
167  for basisVec in basisSet:
168  vecLen = math.sqrt(numpy.sum(numpy.array(basisVec)**2))
169  valueList = []
170  for basisValue in basisVec:
171  if basisValue == 0:
172  valueList.append((0,))
173  else:
174  valueList.append((basisValue, -basisValue))
175  fullSet += list(numpy.array((x, y, z))/vecLen
176  for z in valueList[2]
177  for y in valueList[1]
178  for x in valueList[0]
179  )
180  return fullSet
181 
182 def _findCloseIndexSet(vecList, ind):
183  """Given a list of cartesian vectors, return a set of indices of those closest to one of them
184 
185  This is intended for regular grids where distances are quantized
186 
187  Inputs:
188  - vecList: list of cartesian vectors
189  - ind: index of vector to be nearest
190  """
191  dotProductList = numpy.round(numpy.dot(vecList, vecList[ind]), 2)
192  dotProductList[ind] = -9e99
193  minDist = numpy.max(dotProductList)
194  indList = numpy.arange(len(dotProductList))[dotProductList==minDist]
195  return set(indList)
196 
197 def _findCloseList(vecList, vec):
198  """Given a list of cartesian vectors, return all those closest to a specified position
199 
200  This is intended for regular grids where distances are quantized
201 
202  Inputs:
203  - vecList: list of cartesian vectors
204  - vec: vector to be near
205 
206  @return two items:
207  - list of closest vectors
208  - list if indices of those vectors
209  """
210  dotProductList = numpy.round(numpy.dot(vecList, vec), 2)
211  minDist = numpy.max(dotProductList)
212  indList = numpy.arange(len(dotProductList))[dotProductList==minDist]
213  retList = numpy.take(vecList, indList, 0)
214  return retList, indList
215 
216 def _findClosePair(vecList, ind=0):
217  """Given a list of cartesian vectors and an index, return the vector and one of its closest neighbors
218 
219  Inputs:
220  - vecList: list of cartesian vectors
221  - ind: index of first vector
222  """
223  vec = vecList[ind]
224  otherVecList = vecList[0:ind] + vecList[ind+1:]
225  ind1 = numpy.argmax(numpy.dot(otherVecList, vec))
226  return vec, otherVecList[ind1]
227 
228 def _sortedVectorList(vecList):
229  """Return a list of cartesian vectors sorted by decreasing latitude and increasing longitude
230  """
231  def vecToSort(vec):
232  ang = round(math.atan2(vec[1], vec[0]), 2)
233  if ang < 0:
234  ang += 2.0 * math.pi
235  return (-round(vec[2], 1), ang, vec)
236 
237  decoratedList = [vecToSort(v) for v in vecList]
238  decoratedList.sort()
239  return [d[2] for d in decoratedList]
240 
241 if __name__ == "__main__":
242  numpy.set_printoptions(precision=2, suppress=True, linewidth=120)
243  import lsst.afw.coord as afwCoord
244  import lsst.afw.geom as afwGeom
245 
246  print "Dodecahedron with vertices on poles"
247  vertexDodec = Dodecahedron(withFacesOnPoles=False)
248  for i in range(12):
249  faceVec = vertexDodec.getFaceCtr(i)
250  print "Face %2d: %s" % (i, faceVec)