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