LSSTApplications  10.0+286,10.0+36,10.0+46,10.0-2-g4f67435,10.1+152,10.1+37,11.0,11.0+1,11.0-1-g47edd16,11.0-1-g60db491,11.0-1-g7418c06,11.0-2-g04d2804,11.0-2-g68503cd,11.0-2-g818369d,11.0-2-gb8b8ce7
LSSTDataManagementBasePackage
quadsphere.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 
23 import math
24 
25 import lsst.pex.policy as pexPolicy
26 import lsst.geom as geom
27 
28 
29 # Methods for rotating vectors around the x, y, z and -x, -y, -z axes
30 def _rotX(v, sin, cos):
31  return (v[0], v[1] * cos - v[2] * sin, v[1] * sin + v[2] * cos)
32 def _rotY(v, sin, cos):
33  return (v[0] * cos + v[2] * sin, v[1], -v[0] * sin + v[2] * cos)
34 def _rotZ(v, sin, cos):
35  return (v[0] * cos - v[1] * sin, v[0] * sin + v[1] * cos, v[2])
36 def _rotNX(v, sin, cos):
37  return _rotX(v, -sin, cos)
38 def _rotNY(v, sin, cos):
39  return _rotY(v, -sin, cos)
40 def _rotNZ(v, sin, cos):
41  return _rotZ(v, -sin, cos)
42 
43 
44 class QuadSpherePixelization(object):
45  """A quad-sphere based pixelization of the unit sphere. Each of the 6
46  root pixels are divided into an R by R grid of sky-pixels, where R
47  is the resolution (a positive integer). Sky-pixel identifiers are
48  contiguous integers in range [0, 6 * R**2), and are ordered such that
49  enumerating the pixels 0, 1, ..., 6 * R**2 results in a roughly spiral
50  traversal (from south to north pole) of the unit sphere. In particular,
51  iterating over pixels in this order guarantees that when arriving at
52  pixel P:
53 
54  - all pixels with id P' < P have been visited
55  - all neighbors of pixels with id P' < P - 8R have been visited
56 
57  Each pixel may optionally be expanded by a padding angle A, such that
58  points on the unit sphere outside of a pixel but within angular separation
59  A of the fiducial pixel boundaries are also considered to be part of that
60  pixel. Pixel (and padded pixel) edges are great circles on the unit
61  sphere.
62 
63  This class provides methods for computing fiducial or padded geometric
64  representations of a sky-pixel - spherical convex polygons in both
65  cases. It also supports finding the list of sky-pixels intersecting
66  an input spherical convex polygon, as well as determining the neighbors
67  and center of a sky-pixel.
68 
69  TODO: Right now, pixels are defined by a simple tangent plane projection
70  of each cube face onto the sphere, but pixel boundaries could be adjusted
71  to produce equal area pixels.
72  """
73  def __init__(self, resolution, paddingRad):
74  """Creates a new quad-sphere sky pixelisation.
75 
76  resolution: the number of pixels along the x and y axes of each root
77  pixel.
78 
79  paddingRad: the angular separation (rad) by which fiducial sky-pixels
80  are to be padded.
81  """
82  if not isinstance(resolution, (int, long)):
83  raise TypeError(
84  'Quad-sphere resolution parameter must be an int or long')
85  if resolution < 3 or resolution > 16384:
86  raise RuntimeError(
87  'Quad-sphere resolution must be in range [3, 16384]')
88  if not isinstance(paddingRad, float):
89  raise TypeError(
90  'Quad-sphere pixel padding angle must be a float')
91  if paddingRad < 0.0 or paddingRad >= math.pi * 0.25:
92  raise RuntimeError(
93  'Quad-sphere pixel padding radius must be in range [0, 45) deg')
94  R = resolution
95  self.resolution = R
96  self.padding = paddingRad
97  x, y, z = (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0)
98  nx, ny, nz = (-1.0, 0.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, -1.0)
99  # center of each root pixel
100  self.center = [z, x, y, nx, ny, nz]
101  # x basis vector of each root pixel
102  self.x = [ny, y, nx, ny, x, ny]
103  # y basis vector of each root pixel
104  self.y = [x, z, z, z, z, nx]
105  # functions for rotating vectors in direction of increasing x (per root pixel)
106  self.xrot = [_rotX, _rotZ, _rotZ, _rotZ, _rotZ, _rotNX]
107  # functions for rotating vectors in direction of increasing y (per root pixel)
108  self.yrot = [_rotY, _rotNY, _rotX, _rotY, _rotNX, _rotY]
109  # Compute fiducial and splitting x/y planes for each root pixel
110  self.xplane = []
111  self.yplane = []
112  sp = math.sin(self.padding)
113  for root in xrange(6):
114  xplanes = []
115  yplanes = []
116  c, x, y = self.center[root], self.x[root], self.y[root]
117  for i in xrange(R + 1):
118  xfp = self._fiducialXPlane(root, i)
119  yfp = self._fiducialYPlane(root, i)
120  f = 2.0 * float(i) / float(R) - 1.0
121  thetaX = math.radians(0.5 * geom.cartesianAngularSep(
122  (c[0] + f * x[0] + y[0],
123  c[1] + f * x[1] + y[1],
124  c[2] + f * x[2] + y[2]),
125  (c[0] + f * x[0] - y[0],
126  c[1] + f * x[1] - y[1],
127  c[2] + f * x[2] - y[2])))
128  thetaY = math.radians(0.5 * geom.cartesianAngularSep(
129  (c[0] + x[0] + f * y[0],
130  c[1] + x[1] + f * y[1],
131  c[2] + x[2] + f * y[2]),
132  (c[0] - x[0] + f * y[0],
133  c[1] - x[1] + f * y[1],
134  c[2] - x[2] + f * y[2])))
135  sinrx = sp / math.cos(thetaX)
136  sinry = sp / math.cos(thetaY)
137  cosrx = math.sqrt(1.0 - sinrx * sinrx)
138  cosry = math.sqrt(1.0 - sinry * sinry)
139  if i == 0:
140  xlp = None
141  ybp = None
142  else:
143  xlp = self.xrot[root](xfp, sinrx, cosrx)
144  ybp = self.yrot[root](yfp, sinry, cosry)
145  xlp = (-xlp[0], -xlp[1], -xlp[2])
146  ybp = (-ybp[0], -ybp[1], -ybp[2])
147  if i == R:
148  xrp = None
149  ytp = None
150  else:
151  xrp = self.xrot[root](xfp, -sinrx, cosrx)
152  ytp = self.yrot[root](yfp, -sinry, cosry)
153  xplanes.append((xfp, xlp, xrp))
154  yplanes.append((yfp, ybp, ytp))
155  self.xplane.append(xplanes)
156  self.yplane.append(yplanes)
157  # Corner pixel neighbors.
158  self.cornerNeighbors = (
159  # root pixel 0
160  (
161  # x, y = 0, 0
162  (self.id(0, 1, 0),
163  self.id(0, 1, 1),
164  self.id(0, 0, 1),
165  self.id(2, R - 1, R - 1),
166  self.id(2, R - 2, R - 1),
167  self.id(3, 0, R - 1),
168  self.id(3, 1, R - 1)),
169  # x, y = 0, R-1
170  (self.id(0, 1, R - 1),
171  self.id(0, 1, R - 2),
172  self.id(0, 0, R - 2),
173  self.id(1, R - 1, R - 1),
174  self.id(1, R - 2, R - 1),
175  self.id(2, 0, R - 1),
176  self.id(2, 1, R - 1)),
177  # x, y = R-1, 0
178  (self.id(0, R - 2, 0),
179  self.id(0, R - 2, 1),
180  self.id(0, R - 1, 1),
181  self.id(3, R - 2, R - 1),
182  self.id(3, R - 1, R - 1),
183  self.id(4, 0, R - 1),
184  self.id(4, 1, R - 1)),
185  # x, y = R-1, R-1
186  (self.id(0, R - 2, R - 1),
187  self.id(0, R - 2, R - 2),
188  self.id(0, R - 1, R - 2),
189  self.id(1, 0, R - 1),
190  self.id(1, 1, R - 1),
191  self.id(4, R - 2, R - 1),
192  self.id(4, R - 1, R - 1))
193  ),
194  # root pixel 1
195  (
196  # x, y = 0, 0
197  (self.id(1, 1, 0),
198  self.id(1, 1, 1),
199  self.id(1, 0, 1),
200  self.id(4, R - 1, 0),
201  self.id(4, R - 1, 1),
202  self.id(5, R - 2, 0),
203  self.id(5, R - 1, 0)),
204  # x, y = 0, R-1
205  (self.id(1, 1, R - 1),
206  self.id(1, 1, R - 2),
207  self.id(1, 0, R - 2),
208  self.id(4, R - 1, R - 2),
209  self.id(4, R - 1, R - 1),
210  self.id(0, R - 2, R - 1),
211  self.id(0, R - 1, R - 1)),
212  # x, y = R-1, 0
213  (self.id(1, R - 2, 0),
214  self.id(1, R - 2, 1),
215  self.id(1, R - 1, 1),
216  self.id(2, 0, 0),
217  self.id(2, 0, 1),
218  self.id(5, 0, 0),
219  self.id(5, 1, 0)),
220  # x, y = R-1, R-1
221  (self.id(1, R - 2, R - 1),
222  self.id(1, R - 2, R - 2),
223  self.id(1, R - 1, R - 2),
224  self.id(2, 0, R - 2),
225  self.id(2, 0, R - 1),
226  self.id(0, 0, R - 1),
227  self.id(0, 1, R - 1))
228  ),
229  # root pixel 2
230  (
231  # x, y = 0, 0
232  (self.id(2, 1, 0),
233  self.id(2, 1, 1),
234  self.id(2, 0, 1),
235  self.id(1, R - 1, 0),
236  self.id(1, R - 1, 1),
237  self.id(5, 0, 0),
238  self.id(5, 0, 1)),
239  # x, y = 0, R-1
240  (self.id(2, 1, R - 1),
241  self.id(2, 1, R - 2),
242  self.id(2, 0, R - 2),
243  self.id(1, R - 1, R - 2),
244  self.id(1, R - 1, R - 1),
245  self.id(0, 0, R - 2),
246  self.id(0, 0, R - 1)),
247  # x, y = R-1, 0
248  (self.id(2, R - 2, 0),
249  self.id(2, R - 2, 1),
250  self.id(2, R - 1, 1),
251  self.id(3, 0, 0),
252  self.id(3, 0, 1),
253  self.id(5, 0, R - 2),
254  self.id(5, 0, R - 1)),
255  # x, y = R-1, R-1
256  (self.id(2, R - 2, R - 1),
257  self.id(2, R - 2, R - 2),
258  self.id(2, R - 1, R - 2),
259  self.id(3, 0, R - 2),
260  self.id(3, 0, R - 1),
261  self.id(0, 0, 0),
262  self.id(0, 0, 1))
263  ),
264  # root pixel 3
265  (
266  # x, y = 0, 0
267  (self.id(3, 1, 0),
268  self.id(3, 1, 1),
269  self.id(3, 0, 1),
270  self.id(2, R - 1, 0),
271  self.id(2, R - 1, 1),
272  self.id(5, 0, R - 1),
273  self.id(5, 1, R - 1)),
274  # x, y = 0, R-1
275  (self.id(3, 1, R - 1),
276  self.id(3, 1, R - 2),
277  self.id(3, 0, R - 2),
278  self.id(2, R - 1, R - 2),
279  self.id(2, R - 1, R - 1),
280  self.id(0, 0, 0),
281  self.id(0, 1, 0)),
282  # x, y = R-1, 0
283  (self.id(3, R - 2, 0),
284  self.id(3, R - 2, 1),
285  self.id(3, R - 1, 1),
286  self.id(4, 0, 0),
287  self.id(4, 0, 1),
288  self.id(5, R - 2, R - 1),
289  self.id(5, R - 1, R - 1)),
290  # x, y = R-1, R-1
291  (self.id(3, R - 2, R - 1),
292  self.id(3, R - 2, R - 2),
293  self.id(3, R - 1, R - 2),
294  self.id(4, 0, R - 2),
295  self.id(4, 0, R - 1),
296  self.id(0, R - 2, 0),
297  self.id(0, R - 1, 0))
298  ),
299  # root pixel 4
300  (
301  # x, y = 0, 0
302  (self.id(4, 1, 0),
303  self.id(4, 1, 1),
304  self.id(4, 0, 1),
305  self.id(3, R - 1, 0),
306  self.id(3, R - 1, 1),
307  self.id(5, R - 1, R - 2),
308  self.id(5, R - 1, R - 1)),
309  # x, y = 0, R-1
310  (self.id(4, 1, R - 1),
311  self.id(4, 1, R - 2),
312  self.id(4, 0, R - 2),
313  self.id(3, R - 1, R - 2),
314  self.id(3, R - 1, R - 1),
315  self.id(0, R - 1, 0),
316  self.id(0, R - 1, 1)),
317  # x, y = R-1, 0
318  (self.id(4, R - 2, 0),
319  self.id(4, R - 2, 1),
320  self.id(4, R - 1, 1),
321  self.id(1, 0, 0),
322  self.id(1, 0, 1),
323  self.id(5, R - 1, 0),
324  self.id(5, R - 1, 1)),
325  # x, y = R-1, R-1
326  (self.id(4, R - 2, R - 1),
327  self.id(4, R - 2, R - 2),
328  self.id(4, R - 1, R - 2),
329  self.id(1, 0, R - 2),
330  self.id(1, 0, R - 1),
331  self.id(0, R - 1, R - 2),
332  self.id(0, R - 1, R - 1))
333  ),
334  # root pixel 5
335  (
336  # x, y = 0, 0
337  (self.id(5, 1, 0),
338  self.id(5, 1, 1),
339  self.id(5, 0, 1),
340  self.id(1, R - 2, 0),
341  self.id(1, R - 1, 0),
342  self.id(2, 0, 0),
343  self.id(2, 1, 0)),
344  # x, y = 0, R-1
345  (self.id(5, 1, R - 1),
346  self.id(5, 1, R - 2),
347  self.id(5, 0, R - 2),
348  self.id(2, R - 2, 0),
349  self.id(2, R - 1, 0),
350  self.id(3, 0, 0),
351  self.id(3, 1, 0)),
352  # x, y = R-1, 0
353  (self.id(5, R - 2, 0),
354  self.id(5, R - 2, 1),
355  self.id(5, R - 1, 1),
356  self.id(1, 0, 0),
357  self.id(1, 1, 0),
358  self.id(4, R - 2, 0),
359  self.id(4, R - 1, 0)),
360  # x, y = R-1, R-1
361  (self.id(5, R - 2, R - 1),
362  self.id(5, R - 2, R - 2),
363  self.id(5, R - 1, R - 2),
364  self.id(3, R - 2, 0),
365  self.id(3, R - 1, 0),
366  self.id(4, 0, 0),
367  self.id(4, 1, 0))
368  ),
369  )
370 
371  def getGeometry(self, pixelId, fiducial=False):
372  """Returns a spherical convex polygon corresponding to the fiducial
373  or padded boundaries of the sky-pixel with the specified id.
374  """
375  root, ix, iy = self.coords(pixelId)
376  xl = 2.0 * float(ix) / float(self.resolution) - 1.0
377  xr = 2.0 * float(ix + 1) / float(self.resolution) - 1.0
378  yb = 2.0 * float(iy) / float(self.resolution) - 1.0
379  yt = 2.0 * float(iy + 1) / float(self.resolution) - 1.0
380  c, x, y = self.center[root], self.x[root], self.y[root]
381  v = map(geom.normalize, [(c[0] + xl * x[0] + yb * y[0],
382  c[1] + xl * x[1] + yb * y[1],
383  c[2] + xl * x[2] + yb * y[2]),
384  (c[0] + xr * x[0] + yb * y[0],
385  c[1] + xr * x[1] + yb * y[1],
386  c[2] + xr * x[2] + yb * y[2]),
387  (c[0] + xr * x[0] + yt * y[0],
388  c[1] + xr * x[1] + yt * y[1],
389  c[2] + xr * x[2] + yt * y[2]),
390  (c[0] + xl * x[0] + yt * y[0],
391  c[1] + xl * x[1] + yt * y[1],
392  c[2] + xl * x[2] + yt * y[2])])
393  if not fiducial and self.padding > 0.0:
394  # Determine angles by which edge planes must be rotated outwards
395  sp = math.sin(self.padding)
396  theta = map(lambda x: 0.5 * geom.cartesianAngularSep(x[0], x[1]),
397  ((v[0],v[3]), (v[1],v[0]), (v[2],v[1]), (v[3],v[2])))
398  sina = map(lambda x: sp / math.cos(math.radians(x)), theta)
399  cosa = map(lambda x: math.sqrt(1.0 - x * x), sina)
400  # find plane equations of fiducial pixel boundaries
401  xlp = self.xplane[root][ix][0]
402  ybp = self.yplane[root][iy][0]
403  xrp = self.xplane[root][ix + 1][0]
404  ytp = self.yplane[root][iy + 1][0]
405  # rotate edge planes outwards
406  xlp = self.xrot[root](xlp, -sina[0], cosa[0])
407  ybp = self.yrot[root](ybp, -sina[1], cosa[1])
408  xrp = self.xrot[root](xrp, sina[2], cosa[2])
409  ytp = self.yrot[root](ytp, sina[3], cosa[3])
410  # intersect rotated planes to find vertices of padded sky-pixel
411  v = map(geom.normalize, [geom.cross(xlp, ybp),
412  geom.cross(xrp, ybp),
413  geom.cross(xrp, ytp),
414  geom.cross(xlp, ytp)])
415  return geom.SphericalConvexPolygon(v)
416 
417  def getCenter(self, pixelId):
418  """Returns the center of a sky-pixel as a unit cartesian 3-vector.
419  """
420  root, ix, iy = self.coords(pixelId)
421  xc = 2.0 * (float(ix) + 0.5) / float(self.resolution) - 1.0
422  yc = 2.0 * (float(iy) + 0.5) / float(self.resolution) - 1.0
423  c, x, y = self.center[root], self.x[root], self.y[root]
424  return geom.normalize((c[0] + xc * x[0] + yc * y[0],
425  c[1] + xc * x[1] + yc * y[1],
426  c[2] + xc * x[2] + yc * y[2]))
427 
428  def getNeighbors(self, pixelId):
429  """Returns a list of ids for sky-pixels adjacent to specified pixel.
430  """
431  root, ix, iy = self.coords(pixelId)
432  R = self.resolution
433  if ix == 0:
434  # left edge
435  if iy == 0:
436  return self.cornerNeighbors[root][0]
437  elif iy == R - 1:
438  return self.cornerNeighbors[root][1]
439  elif root == 0:
440  neighbors = (self.id(2, R - iy - 2, R - 1),
441  self.id(2, R - iy - 1, R - 1),
442  self.id(2, R - iy, R - 1))
443  elif root == 5:
444  neighbors = (self.id(2, iy - 1, 0),
445  self.id(2, iy, 0),
446  self.id(2, iy + 1, 0))
447  else:
448  r = root - 1
449  if r == 0: r = 4
450  neighbors = (self.id(r, R - 1, iy - 1),
451  self.id(r, R - 1, iy ),
452  self.id(r, R - 1, iy + 1))
453  return neighbors + (self.id(root, 1, iy - 1),
454  self.id(root, 1, iy ),
455  self.id(root, 1, iy + 1),
456  self.id(root, 0, iy - 1),
457  self.id(root, 0, iy + 1))
458  elif ix == R - 1:
459  # right edge
460  if iy == 0:
461  return self.cornerNeighbors[root][2]
462  elif iy == R - 1:
463  return self.cornerNeighbors[root][3]
464  elif root == 0:
465  neighbors = (self.id(4, iy - 1, R - 1),
466  self.id(4, iy, R - 1),
467  self.id(4, iy + 1, R - 1))
468  elif root == 5:
469  neighbors = (self.id(4, R - iy - 2, 0),
470  self.id(4, R - iy - 1, 0),
471  self.id(4, R - iy, 0))
472  else:
473  r = root + 1
474  if r == 5: r = 1
475  neighbors = (self.id(r, 0, iy - 1),
476  self.id(r, 0, iy ),
477  self.id(r, 0, iy + 1))
478  return neighbors + (self.id(root, R - 2, iy - 1),
479  self.id(root, R - 2, iy),
480  self.id(root, R - 2, iy + 1),
481  self.id(root, R - 1, iy - 1),
482  self.id(root, R - 1, iy + 1))
483  elif iy == 0:
484  # bottom edge
485  if root == 0:
486  neighbors = (self.id(3, ix - 1, R - 1),
487  self.id(3, ix , R - 1),
488  self.id(3, ix + 1, R - 1))
489  elif root == 1:
490  neighbors = (self.id(5, R - ix - 2, 0),
491  self.id(5, R - ix - 1, 0),
492  self.id(5, R - ix, 0))
493  elif root == 2:
494  neighbors = (self.id(5, 0, ix - 1),
495  self.id(5, 0, ix ),
496  self.id(5, 0, ix + 1))
497  elif root == 3:
498  neighbors = (self.id(5, ix - 1, R - 1),
499  self.id(5, ix , R - 1),
500  self.id(5, ix + 1, R - 1))
501  elif root == 4:
502  neighbors = (self.id(5, R - 1, R - ix - 2),
503  self.id(5, R - 1, R - ix - 1),
504  self.id(5, R - 1, R - ix ))
505  else:
506  neighbors = (self.id(1, R - ix - 2, 0),
507  self.id(1, R - ix - 1, 0),
508  self.id(1, R - ix, 0))
509  return neighbors + (self.id(root, ix - 1, 1),
510  self.id(root, ix, 1),
511  self.id(root, ix + 1, 1),
512  self.id(root, ix - 1, 0),
513  self.id(root, ix + 1, 0))
514  elif iy == R - 1:
515  # top edge
516  if root == 0:
517  neighbors = (self.id(1, R - ix - 2, R - 1),
518  self.id(1, R - ix - 1, R - 1),
519  self.id(1, R - ix, R - 1))
520  elif root == 1:
521  neighbors = (self.id(0, R - ix - 2, R - 1),
522  self.id(0, R - ix - 1, R - 1),
523  self.id(0, R - ix, R - 1))
524  elif root == 2:
525  neighbors = (self.id(0, 0, R - ix - 2),
526  self.id(0, 0, R - ix - 1),
527  self.id(0, 0, R - ix ))
528  elif root == 3:
529  neighbors = (self.id(0, ix - 1, 0),
530  self.id(0, ix, 0),
531  self.id(0, ix + 1, 0))
532  elif root == 4:
533  neighbors = (self.id(0, R - 1, ix - 1),
534  self.id(0, R - 1, ix ),
535  self.id(0, R - 1, ix + 1))
536  else:
537  neighbors = (self.id(3, ix - 1, 0),
538  self.id(3, ix, 0),
539  self.id(3, ix + 1, 0))
540  return neighbors + (self.id(root, ix - 1, R - 2),
541  self.id(root, ix, R - 2),
542  self.id(root, ix + 1, R - 2),
543  self.id(root, ix - 1, R - 1),
544  self.id(root, ix + 1, R - 1))
545  # interior pixel
546  return (self.id(root, ix - 1, iy - 1),
547  self.id(root, ix , iy - 1),
548  self.id(root, ix + 1, iy - 1),
549  self.id(root, ix - 1, iy),
550  self.id(root, ix + 1, iy),
551  self.id(root, ix - 1, iy + 1),
552  self.id(root, ix , iy + 1),
553  self.id(root, ix + 1, iy + 1))
554 
555  def intersect(self, polygon):
556  """Returns a list of ids for sky-pixels intersecting the given
557  polygon. Intersection is computed against padded sky-pixels.
558  """
559  if not isinstance(polygon, geom.SphericalConvexPolygon):
560  raise TypeError('Expecting a SphericalConvexPolygon as input ' +
561  'to the sky-pixel intersection computation')
562  pixels = []
563  R = self.resolution
564  for root in xrange(6):
565  # clip against padded root pixel edge planes
566  poly = polygon
567  for p in (self.xplane[root][0][2], self.xplane[root][R][1],
568  self.yplane[root][0][2], self.yplane[root][R][1]):
569  poly = poly.clip(p)
570  if poly == None:
571  break
572  if poly == None:
573  # polygon doesn't intersect this root triangle
574  continue
575  self._intersect(pixels, poly, root, (0, R, 0, R))
576  return pixels
577 
578  def __len__(self):
579  """Returns the total number of sky-pixels in this pixelization.
580  """
581  return 6 * self.resolution * self.resolution
582 
583  def __iter__(self):
584  """Returns an iterator over the sky-pixel ids of all the pixels
585  in this pixelization.
586  """
587  return iter(xrange(len(self)))
588 
589  def __repr__(self):
590  return ''.join([self.__class__.__name__, '(',
591  repr(self.resolution), ', ', repr(self.padding), ')'])
592 
593  def pixel(self, theta, phi):
594  """Maps from spherical coordinates (longitude and latitude
595  angles theta and phi, both in radians) to a sky-pixel id.
596  """
597  R = self.resolution
598  root = int(math.fmod(0.5 + 2.0 * theta / math.pi, 4.0))
599  theta1 = theta - 0.5 * math.pi * root
600  root += 1
601  tanPhi = math.tan(phi)
602  y = tanPhi / math.cos(theta1)
603  if y > 1:
604  root = 0
605  x = -math.sin(theta) / tanPhi
606  y = math.cos(theta) / tanPhi
607  elif y < -1:
608  root = 5
609  x = math.sin(theta) / tanPhi
610  y = math.cos(theta) / tanPhi
611  else:
612  x = math.tan(theta1)
613  x = int(R * 0.5 * (x + 1.0))
614  y = int(R * 0.5 * (y + 1.0))
615  if x >= R:
616  x = R - 1
617  if y >= R:
618  y = R - 1
619  return self.id(root, x, y)
620 
621  def id(self, root, x, y):
622  """Maps from a root pixel number and x, y pixel coordinates
623  to a sky-pixel id.
624  """
625  if not all(isinstance(p, (int, long)) for p in (root, x, y)):
626  raise TypeError('Pixel coordinates must all be of type int or long')
627  if root < 0 or root > 5:
628  raise RuntimeError('root pixel number must be in range [0,6)')
629  R = self.resolution
630  if x < 0 or x >= R or y < 0 or y >= R:
631  raise RuntimeError('x and y pixel coordinates must ' +
632  'be in range [0,%d)' % R)
633  if root > 0 and root < 5:
634  # easy equatorial pixel case
635  return R ** 2 + 4 * R * y + (root - 1) * R + x
636  # assign ids to pixels in expanding concentric squares
637  dx = x - 0.5 * (R - 1)
638  dy = y - 0.5 * (R - 1)
639  ring = max(abs(dx), abs(dy))
640  r = 2.0 * ring - 1.0
641  if r < 0:
642  if root == 5:
643  return 0
644  else:
645  return 6 * R ** 2 - 1
646  if dy == -ring:
647  i = int(r * r + ring - dx)
648  elif dx == -ring:
649  i = int(r * r + 3.0 * ring + dy)
650  elif dy == ring:
651  i = int(r * r + 5.0 * ring + dx)
652  else: # dx == ring
653  i = int(r * r + 7.0 * ring - dy)
654  if root == 5:
655  # south pole - done
656  return i
657  # north pole: use shrinking concentric squares and wind
658  # in the opposite direction
659  return 6 * R ** 2 - i - 1
660 
661  def coords(self, pixelId):
662  """Maps from a sky-pixel id to a root pixel number and x, y
663  pixel coordinates.
664  """
665  if not isinstance(pixelId, (int, long)):
666  raise TypeError(
667  'Sky-pixel id must be an int or long')
668  R = self.resolution
669  R2 = R ** 2
670  if pixelId < 0 or pixelId >= 6 * R2:
671  raise RuntimeError(
672  'Invalid sky-pixel id %d - value must be in range [0, %d)' %
673  (pixelId, 6 * R2))
674  r = pixelId / R2
675  if r > 0 and r < 5:
676  # equatorial pixel case
677  pixelId -= R2
678  y = pixelId / (4 * R)
679  pixelId -= 4 * R * y
680  root = 1 + pixelId / R
681  x = pixelId - (root - 1) * R
682  return (root, x, y)
683  # polar root pixels
684  if r == 5:
685  root = 0
686  i = 6 * R2 - pixelId - 1
687  else:
688  root = 5
689  i = pixelId
690  s = int(math.floor(math.sqrt(i)))
691  if R & 1 == 0:
692  ring = 0.5 * (s + ((s & 1) ^ 1))
693  else:
694  ring = 0.5 * (s + (s & 1))
695  r = 2.0 * ring - 1.0
696  if r < 0.0:
697  return (root, R / 2, R / 2)
698  i -= r * r
699  if i <= 2.0 * ring:
700  dx = ring - i
701  dy = -ring
702  elif i <= 4.0 * ring:
703  dx = -ring
704  dy = i - 3.0 * ring
705  elif i <= 6.0 * ring:
706  dx = i - 5.0 * ring
707  dy = ring
708  else:
709  dx = ring
710  dy = 7.0 * ring - i
711  return (root, int(dx + 0.5 * (R - 1)), int(dy + 0.5 * (R - 1)))
712 
713  def _fiducialXPlane(self, root, ix):
714  assert isinstance(ix, (int,long)) and ix >= 0 and ix <= self.resolution
715  x = 2.0 * float(ix) / float(self.resolution) - 1.0
716  c, b = self.center[root], self.x[root]
717  v = (c[0] + x * b[0], c[1] + x * b[1], c[2] + x * b[2])
718  return geom.normalize(self.xrot[root](v, 1.0, 0.0))
719 
720  def _fiducialYPlane(self, root, iy):
721  assert isinstance(iy, (int,long)) and iy >= 0 and iy <= self.resolution
722  y = 2.0 * float(iy) / float(self.resolution) - 1.0
723  c, b = self.center[root], self.y[root]
724  v = (c[0] + y * b[0], c[1] + y * b[1], c[2] + y * b[2])
725  return geom.normalize(self.yrot[root](v, 1.0, 0.0))
726 
727  def _intersect(self, pixels, poly, root, box):
728  dx = box[1] - box[0]
729  dy = box[3] - box[2]
730  if dx == 1 and dy == 1:
731  i = self.id(root, box[0], box[2])
732  if poly.intersects(self.getGeometry(i)):
733  pixels.append(i)
734  return
735  if dx >= dy:
736  # Split x range and recurse
737  xsplit = box[0] + dx / 2
738  p = poly.clip(self.xplane[root][xsplit][1])
739  if p != None:
740  self._intersect(pixels, p, root, (box[0], xsplit, box[2], box[3]))
741  p = poly.clip(self.xplane[root][xsplit][2])
742  if p != None:
743  self._intersect(pixels, p, root, (xsplit, box[1], box[2], box[3]))
744  else:
745  # Split y range and recurse
746  ysplit = box[2] + dy / 2
747  p = poly.clip(self.yplane[root][ysplit][1])
748  if p != None:
749  self._intersect(pixels, p, root, (box[0], box[1], box[2], ysplit))
750  p = poly.clip(self.yplane[root][ysplit][2])
751  if p != None:
752  self._intersect(pixels, p, root, (box[0], box[1], ysplit, box[3]))
753 
754 
755 def createQuadSpherePixelization(policy=None):
756  policyFile = pexPolicy.DefaultPolicyFile(
757  "skypix", "QuadSpherePixelizationDictionary.paf", "policy")
758  defaults = pexPolicy.Policy.createPolicy(
759  policyFile, policyFile.getRepositoryPath())
760  if policy is None:
761  policy = pexPolicy.Policy()
762  policy.mergeDefaults(defaults)
763  # Obtain resolution and padding from policy
764  resolution = policy.get('resolutionPix')
765  padding = math.radians(policy.get('paddingArcsec') / 3600.0)
766  return QuadSpherePixelization(resolution, padding)
767 
768 def imageToPolygon(wcs, widthPix, heightPix, padRad=0.0):
769  """Computes and returns a spherical convex polygon approximation to the
770  region of the unit sphere covered by an image specified with a WCS and
771  a width/height (pixels). The polygon is computed by connecting the
772  world coordinates of the 4 image corners with great circles, and can
773  optionally be padded by padRad radians.
774  """
775  # Compute image corners
776  cd = wcs.getCDMatrix()
777  xpad = math.degrees(padRad) / math.sqrt(cd[0,0]**2 + cd[0,1]**2)
778  ypad = math.degrees(padRad) / math.sqrt(cd[1,0]**2 + cd[1,1]**2)
779  xmin, ymin = -0.5 - xpad, -0.5 - ypad
780  xmax, ymax = widthPix + xpad - 0.5, heightPix + ypad - 0.5
781  # Produce a lsst.afw.coord.coordLib.Coord object for each vertex
782  coords = [wcs.pixelToSky(xmin, ymin), wcs.pixelToSky(xmax, ymin),
783  wcs.pixelToSky(xmax, ymax), wcs.pixelToSky(xmin, ymax)]
784  # Map these to python tuples containing cartesian unit vectors
785  verts = []
786  for c in coords:
787  verts.append(tuple(c.getVector()))
788  # and turn them into a spherical convex polygon
789  convex, cc = geom.convex(verts)
790  if not convex:
791  raise RuntimeError('Image corners do not form a convex polygon: ' + cc)
792  elif not cc:
793  verts.reverse()
794  return geom.SphericalConvexPolygon(verts)
795 
iter
int iter
Definition: GaussianCentroid.cc:73
lsst::afw::geom::all
bool all(CoordinateExpr< N > const &expr)
Return true if all elements are true.
Definition: CoordinateExpr.h:80
lsst::pex::policy::Policy
a container for holding hierarchical configuration data in memory.
Definition: Policy.h:169
lsst::pex::policy::DefaultPolicyFile
a representation of a default Policy file that is stored as a file in the installation directory of a...
Definition: DefaultPolicyFile.h:69
Generated on Thu Sep 24 2015 02:29:24 for LSSTApplications by   doxygen 1.8.5