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