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
Classes | Functions | Variables
lsst.geom.geometry Namespace Reference

Classes

class  SphericalRegion
 
class  SphericalBox
 
class  SphericalCircle
 
class  SphericalEllipse
 
class  SphericalConvexPolygon
 
class  PartitionMap
 
class  _SubList
 
class  SphericalBoxPartitionMap
 

Functions

def isinf
 
def dot
 
def cross
 
def invScale
 
def normalize
 
def sphericalCoords
 
def cartesianUnitVector
 
def sphericalAngularSep
 
def clampPhi
 
def reduceTheta
 
def northEast
 
def alpha
 
def maxAlpha
 
def cartesianAngularSep
 
def minEdgeSep
 
def minPhiEdgeSep
 
def minThetaEdgeSep
 
def centroid
 
def between
 
def segments
 
def _partition
 
def _medianOfN
 
def _medianOfMedians
 
def median
 
def _prune
 
def _vrange
 
def _gh
 
def _feasible2D
 
def _feasible1D
 
def hemispherical
 
def convex
 
def convexHull
 

Variables

float ARCSEC_PER_DEG = 3600.0
 
float DEG_PER_ARCSEC = 1.0
 
float ANGLE_EPSILON = 0.001
 
float POLE_EPSILON = 1.0
 
float COS_MAX = 1.0
 
int CROSS_N2MIN = 2
 
tuple SIN_MIN = math.sqrt(CROSS_N2MIN)
 
tuple INF = float('inf')
 
tuple NEG_INF = float('-inf')
 
 MAX_FLOAT = float_info.max
 
 MIN_FLOAT = float_info.min
 
 EPSILON = float_info.epsilon
 
 isinf = math.isinf
 

Function Documentation

def lsst.geom.geometry._feasible1D (   points,
  y 
)
private

Definition at line 1672 of file geometry.py.

1672 
1673 def _feasible1D(points, y):
1674  xmin = NEG_INF
1675  xmax = INF
1676  for v in points:
1677  if abs(v[0]) <= MIN_FLOAT:
1678  if y * v[1] <= 0.0:
1679  return False
1680  # inequality is trivially true, skip it
1681  else:
1682  xlim = - y * v[1] / v[0]
1683  if v[0] > 0.0:
1684  xmin = max(xmin, xlim)
1685  else:
1686  xmax = min(xmax, xlim)
1687  if xmax <= xmin:
1688  return False
1689  return True
def lsst.geom.geometry._feasible2D (   points,
  z 
)
private

Definition at line 1601 of file geometry.py.

1601 
1602 def _feasible2D(points, z):
1603  I1 = []
1604  I2 = []
1605  xmin = NEG_INF
1606  xmax = INF
1607  # transform each constraint of the form x*v[0] + y*v[1] + z*v[2] > 0
1608  # into y op a*x + b or x op c, where op is < or >
1609  for v in points:
1610  if abs(v[1]) <= MIN_FLOAT:
1611  if abs(v[0]) <= MIN_FLOAT:
1612  if z * v[2] <= 0.0:
1613  # inequalities trivially lack a solution
1614  return False
1615  # current inequality is trivially true, skip it
1616  else:
1617  xlim = - z * v[2] / v[0]
1618  if v[0] > 0.0:
1619  xmin = max(xmin, xlim)
1620  else:
1621  xmax = min(xmax, xlim)
1622  if xmax <= xmin:
1623  # inequalities trivially lack a solution
1624  return False
1625  else:
1626  # finite since |z|, |v[i]| <= 1 and 1/MIN_FLOAT is finite
1627  coeffs = (v[0] / v[1], - z * v[2] / v[1])
1628  if v[1] > 0.0:
1629  I1.append(coeffs)
1630  else:
1631  I2.append(coeffs)
1632  # At this point (xmin, xmax) is non-empty - if either I1 or I2 is empty
1633  # then a solution trivially exists
1634  if len(I1) == 0 or len(I2) == 0:
1635  return True
1636  # Check for a feasible solution to the inequalities I1 of the form
1637  # form y > a*x + b, I2 of the form y < a*x + b, x > xmin and x < xmax
1638  numConstraints = 0
1639  while True:
1640  intersections = _prune(I1, xmin, xmax, operator.gt)
1641  intersections.extend(_prune(I2, xmin, xmax, operator.lt))
1642  if len(intersections) == 0:
1643  # I1 and I2 each contain exactly one constraint
1644  a1, b1 = I1[0]
1645  a2, b2 = I2[0]
1646  if a1 == a2:
1647  xi = INF
1648  else:
1649  xi = (b2 - b1) / (a1 - a2)
1650  if isinf(xi):
1651  return b1 < b2
1652  return (xi > xmin or a1 < a2) and (xi < xmax or a1 > a2)
1653  elif numConstraints == len(I1) + len(I2):
1654  # No constraints were pruned during search interval refinement,
1655  # and g was not conclusively less than h. Conservatively return
1656  # False to avoid an infinite loop.
1657  return False
1658  numConstraints = len(I1) + len(I2)
1659  x, xerr = median(intersections)
1660  # If g(x) < h(x), x is a feasible solution. Otherwise, refine the
1661  # search interval by examining the one-sided derivates of g/h.
1662  gmin, gmax, sg, Sg = _gh(I1, x, xerr, operator.gt)
1663  hmin, hmax, sh, Sh = _gh(I2, x, xerr, operator.lt)
1664  if gmax <= hmin:
1665  return True
1666  elif sg > Sh:
1667  xmax = x + xerr
1668  elif Sg < sh:
1669  xmin = x - xerr
1670  else:
1671  return False
def lsst.geom.geometry._gh (   constraints,
  x,
  xerr,
  op 
)
private

Definition at line 1583 of file geometry.py.

1583 
1584 def _gh(constraints, x, xerr, op):
1585  a, b = constraints[0]
1586  amin, amax = a, a
1587  vmin, vmax = _vrange(x, xerr, a, b)
1588  for i in xrange(1, len(constraints)):
1589  a, b = constraints[i]
1590  vimin, vimax = _vrange(x, xerr, a, b)
1591  if vimax < vmin or vimin > vmax:
1592  if op(vimin, vmin):
1593  amin = a
1594  amax = a
1595  vmin = vimin
1596  vmax = vimax
1597  else:
1598  amin = min(amin, a)
1599  amax = max(amax, a)
1600  return vmin, vmax, amin, amax
def lsst.geom.geometry._medianOfMedians (   array,
  left,
  right 
)
private 
Returns the "median of medians" for an array. This primitive is used
for pivot selection in the median finding algorithm.

Definition at line 1482 of file geometry.py.

1482 
1483 def _medianOfMedians(array, left, right):
1484  """Returns the "median of medians" for an array. This primitive is used
1485  for pivot selection in the median finding algorithm.
1486  """
1487  while True:
1488  if right - left <= 5:
1489  return _medianOfN(array, left, right - left)
1490  i = left
1491  j = left
1492  while i + 4 < right:
1493  k = _medianOfN(array, i, 5)
1494  tmp = array[j]
1495  array[j] = array[k]
1496  array[k] = tmp
1497  j += 1
1498  i += 5
1499  right = j
def lsst.geom.geometry._medianOfN (   array,
  i,
  n 
)
private 
Finds the median of n consecutive elements in an array starting
at index i (efficient for small n). The index of the median element
is returned.

Definition at line 1460 of file geometry.py.

1460 
1461 def _medianOfN(array, i, n):
1462  """Finds the median of n consecutive elements in an array starting
1463  at index i (efficient for small n). The index of the median element
1464  is returned.
1465  """
1466  if n == 1:
1467  return i
1468  k = n >> 1
1469  e1 = i + k + 1
1470  e2 = i + n
1471  for j in xrange(i, e1):
1472  minIndex = j
1473  minValue = array[j]
1474  for s in xrange(j + 1, e2):
1475  if array[s] < minValue:
1476  minIndex = s
1477  minValue = array[s]
1478  tmp = array[j]
1479  array[j] = array[minIndex]
1480  array[minIndex] = tmp
1481  return i + k
def lsst.geom.geometry._partition (   array,
  left,
  right,
  i 
)
private 
Partitions an array around the pivot value at index i,
returning the new index of the pivot value.

Definition at line 1443 of file geometry.py.

1443 
1444 def _partition(array, left, right, i):
1445  """Partitions an array around the pivot value at index i,
1446  returning the new index of the pivot value.
1447  """
1448  pivot = array[i]
1449  array[i] = array[right - 1]
1450  j = left
1451  for k in xrange(left, right - 1):
1452  if array[k] < pivot:
1453  tmp = array[j]
1454  array[j] = array[k]
1455  array[k] = tmp
1456  j += 1
1457  array[right - 1] = array[j]
1458  array[j] = pivot
1459  return j
def lsst.geom.geometry._prune (   constraints,
  xmin,
  xmax,
  op 
)
private 
Removes redundant constraints and reports intersection points
of non-redundant pairs. The constraint list is modified in place.

Definition at line 1528 of file geometry.py.

1528 
1529 def _prune(constraints, xmin, xmax, op):
1530  """Removes redundant constraints and reports intersection points
1531  of non-redundant pairs. The constraint list is modified in place.
1532  """
1533  intersections = []
1534  i = 0
1535  while i < len(constraints) - 1:
1536  a1, b1 = constraints[i]
1537  a2, b2 = constraints[i + 1]
1538  # if constraints are near parallel or intersect to the left or right
1539  # of the feasible x range, remove the higher/lower constraint
1540  da = a1 - a2
1541  if abs(da) < MIN_FLOAT / EPSILON:
1542  xi = INF
1543  else:
1544  xi = (b2 - b1) / da
1545  xierr = 2.0 * EPSILON * abs(xi)
1546  if isinf(xi):
1547  if op(b1, b2):
1548  constraints[i + 1] = constraints[-1]
1549  else:
1550  constraints[i] = constraints[-1]
1551  constraints.pop()
1552  else:
1553  if xi + xierr <= xmin:
1554  if op(a1, a2):
1555  constraints[i + 1] = constraints[-1]
1556  else:
1557  constraints[i] = constraints[-1]
1558  constraints.pop()
1559  elif xi - xierr >= xmax:
1560  if op(a1, a2):
1561  constraints[i] = constraints[-1]
1562  else:
1563  constraints[i + 1] = constraints[-1]
1564  constraints.pop()
1565  else:
1566  # save intersection for later
1567  intersections.append((xi, xierr))
1568  i += 2
1569  return intersections
def lsst.geom.geometry._vrange (   x,
  xerr,
  a,
  b 
)
private

Definition at line 1570 of file geometry.py.

1570 
1571 def _vrange(x, xerr, a, b):
1572  p = a * (x - xerr)
1573  v = p + b
1574  verr = EPSILON * abs(p) + EPSILON * abs(v)
1575  vmin = v - verr
1576  vmax = v + verr
1577  p = a * (x + xerr)
1578  v = p + b
1579  verr = EPSILON * abs(p) + EPSILON * abs(v)
1580  vmin = min(vmin, v - verr)
1581  vmax = max(vmax, v + verr)
1582  return vmin, vmax
def lsst.geom.geometry.alpha (   r,
  centerPhi,
  phi 
) 
Returns the longitude angle alpha of the intersections (alpha, phi),
(-alpha, phi) of the circle centered on (0, centerPhi) with radius r and
the plane z = sin(phi). If there is no intersection, None is returned.

Definition at line 206 of file geometry.py.

206 
207 def alpha(r, centerPhi, phi):
208  """Returns the longitude angle alpha of the intersections (alpha, phi),
209  (-alpha, phi) of the circle centered on (0, centerPhi) with radius r and
210  the plane z = sin(phi). If there is no intersection, None is returned.
211  """
212  if phi < centerPhi - r or phi > centerPhi + r:
213  return None
214  a = abs(centerPhi - phi)
215  if a <= r - (90.0 - phi) or a <= r - (90.0 + phi):
216  return None
217  p = math.radians(phi)
218  cp = math.radians(centerPhi)
219  x = math.cos(math.radians(r)) - math.sin(cp) * math.sin(cp)
220  u = math.cos(cp) * math.cos(p)
221  y = math.sqrt(abs(u * u - x * x))
222  return math.degrees(abs(math.atan2(y, x)))
def lsst.geom.geometry.between (   p,
  n,
  v1,
  v2 
) 
Tests whether p lies on the shortest great circle arc from cartesian
unit vectors v1 to v2, assuming that p is a unit vector on the plane
defined by the origin, v1 and v2.

Definition at line 306 of file geometry.py.

306 
307 def between(p, n, v1, v2):
308  """Tests whether p lies on the shortest great circle arc from cartesian
309  unit vectors v1 to v2, assuming that p is a unit vector on the plane
310  defined by the origin, v1 and v2.
311  """
312  p1 = cross(n, v1)
313  p2 = cross(n, v2)
314  return dot(p1, p) >= 0.0 and dot(p2, p) <= 0.0
def lsst.geom.geometry.cartesianAngularSep (   v1,
  v2 
) 
Returns the angular separation in degrees between points v1 and v2,
which must both be specified as cartesian 3-vectors.

Definition at line 242 of file geometry.py.

242 
243 def cartesianAngularSep(v1, v2):
244  """Returns the angular separation in degrees between points v1 and v2,
245  which must both be specified as cartesian 3-vectors.
246  """
247  cs = dot(v1, v2)
248  n = cross(v1, v2)
249  ss = math.sqrt(dot(n, n))
250  if cs == 0.0 and ss == 0.0:
251  return 0.0
252  return math.degrees(math.atan2(ss, cs))
def lsst.geom.geometry.cartesianUnitVector (   args) 
Returns a unit cartesian 3-vector corresponding to the input
coordinates, which can be spherical or 3D cartesian. The 2 (spherical)
or 3 (cartesian 3-vector) input coordinates can either be passed
individually or as a tuple/list, and can be of any type convertible
to a float.

Definition at line 130 of file geometry.py.

130 
131 def cartesianUnitVector(*args):
132  """Returns a unit cartesian 3-vector corresponding to the input
133  coordinates, which can be spherical or 3D cartesian. The 2 (spherical)
134  or 3 (cartesian 3-vector) input coordinates can either be passed
135  individually or as a tuple/list, and can be of any type convertible
136  to a float.
137  """
138  if len(args) == 1:
139  args = args[0]
140  if len(args) == 2:
141  theta = math.radians(float(args[0]))
142  phi = math.radians(float(args[1]))
143  cosPhi = math.cos(phi)
144  return (math.cos(theta) * cosPhi,
145  math.sin(theta) * cosPhi,
146  math.sin(phi))
147  elif len(args) == 3:
148  x = float(args[0])
149  y = float(args[1])
150  z = float(args[2])
151  n = math.sqrt(x * x + y * y + z * z)
152  if n == 0.0:
153  raise RuntimeError('Cannot normalize a 3-vector with 0 magnitude')
154  return (x / n, y / n, z / n)
155  raise TypeError('Expecting 2 or 3 coordinates, or a tuple/list ' +
156  'containing 2 or 3 coordinates')
def lsst.geom.geometry.centroid (   vertices) 
Computes the centroid of a set of vertices (which must be passed in
as a list of cartesian unit vectors) on the unit sphere.

Definition at line 294 of file geometry.py.

294 
295 def centroid(vertices):
296  """Computes the centroid of a set of vertices (which must be passed in
297  as a list of cartesian unit vectors) on the unit sphere.
298  """
299  x, y, z = 0.0, 0.0, 0.0
300  # Note: could use more accurate summation routines
301  for v in vertices:
302  x += v[0]
303  y += v[1]
304  z += v[2]
305  return normalize((x, y, z))
def lsst.geom.geometry.clampPhi (   phi) 
Clamps the input latitude angle phi to [-90.0, 90.0] deg.

Definition at line 171 of file geometry.py.

171 
172 def clampPhi(phi):
173  """Clamps the input latitude angle phi to [-90.0, 90.0] deg.
174  """
175  if phi < -90.0:
176  return -90.0
177  elif phi >= 90.0:
178  return 90.0
179  return phi
def lsst.geom.geometry.convex (   vertices) 
Tests whether an ordered list of vertices (which must be specified
as cartesian unit vectors) form a spherical convex polygon and determines
their winding order. Returns a 2-tuple as follows:

(True, True):   The vertex list forms a spherical convex polygon and is in
                counter-clockwise order when viewed from outside the unit
                sphere in a right handed coordinate system.
(True, False):  The vertex list forms a spherical convex polygon and is in
                in clockwise order.
(False, msg):   The vertex list does not form a spherical convex polygon -
                msg is a string describing why the test failed.

The algorithm completes in O(n) time, where n is the number of
input vertices.

Definition at line 1732 of file geometry.py.

1732 
1733 def convex(vertices):
1734  """Tests whether an ordered list of vertices (which must be specified
1735  as cartesian unit vectors) form a spherical convex polygon and determines
1736  their winding order. Returns a 2-tuple as follows:
1737 
1738  (True, True): The vertex list forms a spherical convex polygon and is in
1739  counter-clockwise order when viewed from outside the unit
1740  sphere in a right handed coordinate system.
1741  (True, False): The vertex list forms a spherical convex polygon and is in
1742  in clockwise order.
1743  (False, msg): The vertex list does not form a spherical convex polygon -
1744  msg is a string describing why the test failed.
1745 
1746  The algorithm completes in O(n) time, where n is the number of
1747  input vertices.
1748  """
1749  if len(vertices) < 3:
1750  return (False, '3 or more vertices must be specified')
1751  if not hemispherical(vertices):
1752  return (False, 'vertices are not hemispherical')
1753  center = centroid(vertices)
1754  windingAngle = 0.0
1755  clockwise = False
1756  counterClockwise = False
1757  p1 = cross(center, vertices[-1])
1758  n2 = dot(p1, p1)
1759  if abs(n2) < CROSS_N2MIN:
1760  return (False, 'centroid of vertices is too close to a vertex')
1761  for i in xrange(len(vertices)):
1762  beg = vertices[i - 2]
1763  mid = vertices[i - 1]
1764  end = vertices[i]
1765  plane = cross(mid, end)
1766  n2 = dot(plane, plane)
1767  if dot(mid, end) >= COS_MAX or n2 < CROSS_N2MIN:
1768  return (False, 'vertex list contains [near-]duplicates')
1769  plane = invScale(plane, math.sqrt(n2))
1770  d = dot(plane, beg)
1771  if d > SIN_MIN:
1772  if clockwise:
1773  return (False, 'vertices wind around centroid in both ' +
1774  'clockwise and counter-clockwise order')
1775  counterClockwise = True
1776  elif d < -SIN_MIN:
1777  if counterClockwise:
1778  return (False, 'vertices wind around centroid in both ' +
1779  'clockwise and counter-clockwise order')
1780  clockwise = True
1781  # center must be inside polygon formed by vertices if they form a
1782  # convex polygon - an equivalent check is to test that polygon
1783  # vertices always wind around center in the same direction.
1784  d = dot(plane, center)
1785  if d < SIN_MIN and counterClockwise or d > -SIN_MIN and clockwise:
1786  return (False, 'centroid of vertices is not conclusively ' +
1787  'inside all edges')
1788  # sum up winding angle for edge (mid, end)
1789  p2 = cross(center, end)
1790  n2 = dot(p2, p2)
1791  if abs(n2) < CROSS_N2MIN:
1792  return (False, 'centroid of vertices is too close to a vertex')
1793  windingAngle += cartesianAngularSep(p1, p2)
1794  p1 = p2
1795  # For convex polygons, the closest multiple of 360 to
1796  # windingAngle is 1.
1797  m = math.floor(windingAngle / 360.0)
1798  if m == 0.0 and windingAngle > 180.0 or m == 1.0 and windingAngle < 540.0:
1799  return (True, counterClockwise)
1800  return (False, 'vertices do not completely wind around centroid, or ' +
1801  'wind around it multiple times')
def lsst.geom.geometry.convexHull (   points) 
Computes the convex hull (a spherical convex polygon) of an unordered
list of points on the unit sphere, which must be passed in as cartesian
unit vectors. The algorithm takes O(n log n) time, where n is the number
of points.

Definition at line 1802 of file geometry.py.

1802 
1803 def convexHull(points):
1804  """Computes the convex hull (a spherical convex polygon) of an unordered
1805  list of points on the unit sphere, which must be passed in as cartesian
1806  unit vectors. The algorithm takes O(n log n) time, where n is the number
1807  of points.
1808  """
1809  if len(points) < 3:
1810  return None
1811  if not hemispherical(points):
1812  return None
1813  center = centroid(points)
1814  maxSep = 0.0
1815  extremum = 0
1816  # Vertex furthest from the center is on the hull
1817  for i in xrange(len(points)):
1818  sep = cartesianAngularSep(points[i], center)
1819  if sep > maxSep:
1820  extremum = i
1821  maxSep = sep
1822  anchor = points[extremum]
1823  refPlane = cross(center, anchor)
1824  n2 = dot(refPlane, refPlane)
1825  if n2 < CROSS_N2MIN:
1826  # vertex and center are too close
1827  return None
1828  refPlane = invScale(refPlane, math.sqrt(n2))
1829  # Order points by winding angle from the first (extreme) vertex
1830  av = [(0.0, anchor)]
1831  for i in xrange(0, len(points)):
1832  if i == extremum:
1833  continue
1834  v = points[i]
1835  plane = cross(center, v)
1836  n2 = dot(plane, plane)
1837  if n2 >= CROSS_N2MIN:
1838  plane = invScale(plane, math.sqrt(n2))
1839  p = cross(refPlane, plane)
1840  sa = math.sqrt(dot(p, p))
1841  if dot(p, center) < 0.0:
1842  sa = -sa
1843  angle = math.atan2(sa, dot(refPlane, plane))
1844  if angle < 0.0:
1845  angle += 2.0 * math.pi
1846  av.append((angle, v))
1847  if len(av) < 3:
1848  return None
1849  av.sort(key=lambda t: t[0]) # stable, so av[0] still contains anchor
1850  # Loop over vertices using a Graham scan adapted for spherical geometry.
1851  # Chan's algorithm would be an improvement, but seems likely to be slower
1852  # for small vertex counts (the expected case).
1853  verts = [anchor]
1854  edges = [None]
1855  edge = None
1856  i = 1
1857  while i < len(av):
1858  v = av[i][1]
1859  p = cross(anchor, v)
1860  n2 = dot(p, p)
1861  if dot(anchor, v) < COS_MAX and n2 >= CROSS_N2MIN:
1862  if edge == None:
1863  # Compute first edge
1864  edge = invScale(p, math.sqrt(n2))
1865  verts.append(v)
1866  edges.append(edge)
1867  anchor = v
1868  else:
1869  d = dot(v, edge)
1870  if d > SIN_MIN:
1871  # v is inside the edge defined by the last
1872  # 2 vertices on the hull
1873  edge = invScale(p, math.sqrt(n2))
1874  verts.append(v)
1875  edges.append(edge)
1876  anchor = v
1877  elif d < -SIN_MIN:
1878  # backtrack - the most recently added hull vertex
1879  # is not actually on the hull.
1880  verts.pop()
1881  edges.pop()
1882  anchor = verts[-1]
1883  edge = edges[-1]
1884  # reprocess v to decide whether to add it to the hull
1885  # or whether further backtracking is necessary
1886  continue
1887  # v is coplanar with edge, skip it
1888  i += 1
1889  # Handle backtracking necessary for last edge
1890  if len(verts) < 3:
1891  return None
1892  v = verts[0]
1893  while True:
1894  p = cross(anchor, v)
1895  n2 = dot(p, p)
1896  if dot(anchor, v) < COS_MAX and n2 >= CROSS_N2MIN:
1897  if dot(v, edge) > SIN_MIN:
1898  edges[0] = invScale(p, math.sqrt(n2))
1899  break
1900  verts.pop()
1901  edges.pop()
1902  anchor = verts[-1]
1903  edge = edges[-1]
1904  if len(verts) < 3:
1905  return None
1906  return SphericalConvexPolygon(verts, edges)
1907 
1908 
1909 # -- Generic unit sphere partitioning scheme base class ----
1910 
# TODO: make this an ABC once python 2.6 becomes the default
def lsst.geom.geometry.cross (   v1,
  v2 
) 
Returns the cross product of cartesian 3-vectors v1 and v2.

Definition at line 78 of file geometry.py.

78 
79 def cross(v1, v2):
80  """Returns the cross product of cartesian 3-vectors v1 and v2.
81  """
82  return (v1[1] * v2[2] - v1[2] * v2[1],
83  v1[2] * v2[0] - v1[0] * v2[2],
84  v1[0] * v2[1] - v1[1] * v2[0])
def lsst.geom.geometry.dot (   v1,
  v2 
) 
Returns the dot product of cartesian 3-vectors v1 and v2.

Definition at line 73 of file geometry.py.

73 
74 def dot(v1, v2):
75  """Returns the dot product of cartesian 3-vectors v1 and v2.
76  """
77  return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]
def lsst.geom.geometry.hemispherical (   points) 
Tests whether a set of points is hemispherical, i.e. whether a plane
exists such that all points are strictly on one side of the plane. The
algorithm used is Megiddo's algorithm for linear programming in R2 and
has run-time O(n), where n is the number of points. Points must be passed
in as a list of cartesian unit vectors.

Definition at line 1690 of file geometry.py.

1690 
1691 def hemispherical(points):
1692  """Tests whether a set of points is hemispherical, i.e. whether a plane
1693  exists such that all points are strictly on one side of the plane. The
1694  algorithm used is Megiddo's algorithm for linear programming in R2 and
1695  has run-time O(n), where n is the number of points. Points must be passed
1696  in as a list of cartesian unit vectors.
1697  """
1698  # Check whether the set of linear equations
1699  # x * v[0] + y * v[1] + z * v[2] > 0.0 (for v in points)
1700  # has a solution (x, y, z). If (x,y,z) is a solution (is feasible),
1701  # so is C*(x,y,z), C > 0. Therefore we can fix z to 1, -1 and
1702  # perform 2D feasibility tests.
1703  if _feasible2D(points, 1.0):
1704  return True
1705  if _feasible2D(points, -1.0):
1706  return True
1707  # At this point a feasible solution must have z = 0. Fix y to 1, -1 and
1708  # perform 1D feasibility tests.
1709  if _feasible1D(points, 1.0):
1710  return True
1711  if _feasible1D(points, -1.0):
1712  return True
1713  # At this point a feasible solution must have y = z = 0. If all v[0]
1714  # are non-zero and have the same sign, then there is a feasible solution.
1715  havePos = False
1716  haveNeg = False
1717  for v in points:
1718  if v[0] > 0.0:
1719  if haveNeg:
1720  return False
1721  havePos = True
1722  elif v[0] < 0.0:
1723  if havePos:
1724  return False
1725  haveNeg = True
1726  else:
1727  return False
1728  return True
1729 
1730 
1731 # -- Convexity test and convex hull algorithm ----
def lsst.geom.geometry.invScale (   v,
  s 
) 
Returns a copy of the cartesian 3-vector v scaled by 1 / s.

Definition at line 85 of file geometry.py.

85 
86 def invScale(v, s):
87  """Returns a copy of the cartesian 3-vector v scaled by 1 / s.
88  """
89  return (v[0] / s, v[1] / s, v[2] / s)
def lsst.geom.geometry.isinf (   x)

Definition at line 68 of file geometry.py.

68 
69  def isinf(x):
70  return x == INF or x == NEG_INF
71 
72 # -- Utility methods ----
def lsst.geom.geometry.maxAlpha (   r,
  centerPhi 
) 
Computes alpha, the extent in longitude angle [-alpha, alpha] of the
circle with radius r and center (0, centerPhi) on the unit sphere.
Both r and centerPhi are assumed to be in units of degrees.
centerPhi is clamped to lie in the range [-90,90] and r must
lie in the range [0, 90].

Definition at line 223 of file geometry.py.

223 
224 def maxAlpha(r, centerPhi):
225  """Computes alpha, the extent in longitude angle [-alpha, alpha] of the
226  circle with radius r and center (0, centerPhi) on the unit sphere.
227  Both r and centerPhi are assumed to be in units of degrees.
228  centerPhi is clamped to lie in the range [-90,90] and r must
229  lie in the range [0, 90].
230  """
231  assert r >= 0.0 and r <= 90.0
232  if r == 0.0:
233  return 0.0
234  centerPhi = clampPhi(centerPhi)
235  if abs(centerPhi) + r > 90.0 - POLE_EPSILON:
236  return 180.0
237  r = math.radians(r)
238  c = math.radians(centerPhi)
239  y = math.sin(r)
240  x = math.sqrt(abs(math.cos(c - r) * math.cos(c + r)))
241  return math.degrees(abs(math.atan(y / x)))
def lsst.geom.geometry.median (   array) 
Finds the median element of the given array in linear time.
Examples:
imageStatistics.cc.

Definition at line 1500 of file geometry.py.

1500 
1501 def median(array):
1502  """Finds the median element of the given array in linear time.
1503  """
1504  left = 0
1505  right = len(array)
1506  if right == 0:
1507  return None
1508  k = right >> 1
1509  while True:
1510  i = _medianOfMedians(array, left, right)
1511  i = _partition(array, left, right, i)
1512  if k == i:
1513  return array[k]
1514  elif k < i:
1515  right = i
1516  else:
1517  left = i + 1
1518 
1519 
1520 # -- Testing whether a set of points is hemispherical ----
1521 #
1522 # This test is used by the convexity test and convex hull algorithm. It is
1523 # implemented using Megiddo's algorithm for linear programming in R2, see:
1524 #
1525 # Megiddo, N. 1982. Linear-time algorithms for linear programming in R3 and related problems.
1526 # In Proceedings of the 23rd Annual Symposium on Foundations of Computer Science (November 03 - 05, 1982).
1527 # SFCS. IEEE Computer Society, Washington, DC, 329-338. DOI= http://dx.doi.org/10.1109/SFCS.1982.74
def lsst.geom.geometry.minEdgeSep (   p,
  n,
  v1,
  v2 
) 
Returns the minimum angular separation in degrees between p
and points on the great circle edge with plane normal n and
vertices v1, v2. All inputs must be unit cartesian 3-vectors.

Definition at line 253 of file geometry.py.

253 
254 def minEdgeSep(p, n, v1, v2):
255  """Returns the minimum angular separation in degrees between p
256  and points on the great circle edge with plane normal n and
257  vertices v1, v2. All inputs must be unit cartesian 3-vectors.
258  """
259  p1 = cross(n, v1)
260  p2 = cross(n, v2)
261  if dot(p1, p) >= 0.0 and dot(p2, p) <= 0.0:
262  return abs(90.0 - cartesianAngularSep(p, n))
263  else:
264  return min(cartesianAngularSep(p, v1), cartesianAngularSep(p, v2))
def lsst.geom.geometry.minPhiEdgeSep (   p,
  phi,
  minTheta,
  maxTheta 
) 
Returns the minimum angular separation in degrees between p
and points on the small circle edge with constant latitude angle
phi and vertices (minTheta, phi), (maxTheta, phi). p must be in
spherical coordinates.

Definition at line 265 of file geometry.py.

265 
266 def minPhiEdgeSep(p, phi, minTheta, maxTheta):
267  """Returns the minimum angular separation in degrees between p
268  and points on the small circle edge with constant latitude angle
269  phi and vertices (minTheta, phi), (maxTheta, phi). p must be in
270  spherical coordinates.
271  """
272  if minTheta > maxTheta:
273  if p[0] >= minTheta or p[0] <= maxTheta:
274  return abs(p[1] - phi)
275  else:
276  if p[0] >= minTheta and p[0] <= maxTheta:
277  return abs(p[1] - phi)
278  return min(sphericalAngularSep(p, (minTheta, phi)),
279  sphericalAngularSep(p, (maxTheta, phi)))
def lsst.geom.geometry.minThetaEdgeSep (   p,
  theta,
  minPhi,
  maxPhi 
) 
Returns the minimum angular separation in degrees between p
and points on the great circle edge with constant longitude angle
theta and vertices (theta, minPhi), (theta, maxPhi). p must be a
unit cartesian 3-vector.

Definition at line 280 of file geometry.py.

280 
281 def minThetaEdgeSep(p, theta, minPhi, maxPhi):
282  """Returns the minimum angular separation in degrees between p
283  and points on the great circle edge with constant longitude angle
284  theta and vertices (theta, minPhi), (theta, maxPhi). p must be a
285  unit cartesian 3-vector.
286  """
287  v1 = cartesianUnitVector(theta, minPhi)
288  v2 = cartesianUnitVector(theta, maxPhi)
289  n = cross(v1, v2)
290  d2 = dot(n, n)
291  if d2 == 0.0:
292  return min(cartesianAngularSep(p, v1), cartesianAngularSep(p, v2))
293  return minEdgeSep(p, invScale(n, math.sqrt(d2)), v1, v2)
def lsst.geom.geometry.normalize (   v) 
Returns a normalized copy of the cartesian 3-vector v.

Definition at line 90 of file geometry.py.

90 
91 def normalize(v):
92  """Returns a normalized copy of the cartesian 3-vector v.
93  """
94  n = math.sqrt(dot(v, v))
95  if n == 0.0:
96  raise RuntimeError('Cannot normalize a 3-vector with 0 magnitude')
97  return (v[0] / n, v[1] / n, v[2] / n)
def lsst.geom.geometry.northEast (   v) 
Returns unit N,E basis vectors for a point v, which must be a
cartesian 3-vector.

Definition at line 190 of file geometry.py.

190 
191 def northEast(v):
192  """Returns unit N,E basis vectors for a point v, which must be a
193  cartesian 3-vector.
194  """
195  north = (-v[0] * v[2],
196  -v[1] * v[2],
197  v[0] * v[0] + v[1] * v[1])
198  if north == (0.0, 0.0, 0.0):
199  # pick an arbitrary orthogonal basis with z = 0
200  north = (-1.0, 0.0, 0.0)
201  east = (0.0, 1.0, 0.0)
202  else:
203  north = normalize(north)
204  east = normalize(cross(north, v))
205  return north, east
def lsst.geom.geometry.reduceTheta (   theta) 
Range reduces the given longitude angle to lie in the range
[0.0, 360.0).

Definition at line 180 of file geometry.py.

180 
181 def reduceTheta(theta):
182  """Range reduces the given longitude angle to lie in the range
183  [0.0, 360.0).
184  """
185  theta = math.fmod(theta, 360.0)
186  if theta < 0.0:
187  return theta + 360.0
188  else:
189  return theta
def lsst.geom.geometry.segments (   phiMin,
  phiMax,
  width 
) 
Computes the number of segments to divide the given latitude angle
range [phiMin, phiMax] (degrees) into. Two points within the range
separated by at least one segment are guaranteed to have angular
separation of at least width degrees.

Definition at line 315 of file geometry.py.

315 
316 def segments(phiMin, phiMax, width):
317  """Computes the number of segments to divide the given latitude angle
318  range [phiMin, phiMax] (degrees) into. Two points within the range
319  separated by at least one segment are guaranteed to have angular
320  separation of at least width degrees.
321  """
322  p = max(abs(phiMin), abs(phiMax))
323  if p > 90.0 - 1.0 * DEG_PER_ARCSEC:
324  return 1
325  if width >= 180.0:
326  return 1
327  elif width < 1.0 * DEG_PER_ARCSEC:
328  width = 1.0 * DEG_PER_ARCSEC
329  p = math.radians(p)
330  cw = math.cos(math.radians(width))
331  sp = math.sin(p)
332  cp = math.cos(p)
333  x = cw - sp * sp
334  u = cp * cp
335  y = math.sqrt(abs(u * u - x * x))
336  return int(math.floor((2.0 * math.pi) / abs(math.atan2(y, x))))
337 
338 
339 # -- Regions on the unit sphere ----
340 
# TODO: make this an ABC once python 2.6 becomes the default
def lsst.geom.geometry.sphericalAngularSep (   p1,
  p2 
) 
Returns the angular separation in degrees between points p1 and p2,
which must both be specified in spherical coordinates. The implementation
uses the halversine distance formula.

Definition at line 157 of file geometry.py.

157 
158 def sphericalAngularSep(p1, p2):
159  """Returns the angular separation in degrees between points p1 and p2,
160  which must both be specified in spherical coordinates. The implementation
161  uses the halversine distance formula.
162  """
163  sdt = math.sin(math.radians(p1[0] - p2[0]) * 0.5)
164  sdp = math.sin(math.radians(p1[1] - p2[1]) * 0.5)
165  cc = math.cos(math.radians(p1[1])) * math.cos(math.radians(p2[1]))
166  s = math.sqrt(sdp * sdp + cc * sdt * sdt)
167  if s > 1.0:
168  return 180.0
169  else:
170  return 2.0 * math.degrees(math.asin(s))
def lsst.geom.geometry.sphericalCoords (   args) 
Returns spherical coordinates in degrees for the input coordinates,
which can be spherical or 3D cartesian. The 2 (spherical) or 3
(cartesian 3-vector) inputs can be passed either individually
or as a tuple/list, and can be of any type convertible to a float.

Definition at line 98 of file geometry.py.

98 
99 def sphericalCoords(*args):
100  """Returns spherical coordinates in degrees for the input coordinates,
101  which can be spherical or 3D cartesian. The 2 (spherical) or 3
102  (cartesian 3-vector) inputs can be passed either individually
103  or as a tuple/list, and can be of any type convertible to a float.
104  """
105  if len(args) == 1:
106  args = args[0]
107  if len(args) == 2:
108  t = (float(args[0]), float(args[1]))
109  if t[1] < -90.0 or t[1] > 90.0:
110  raise RuntimeError('Latitude angle is out of bounds')
111  return t
112  elif len(args) == 3:
113  x = float(args[0])
114  y = float(args[1])
115  z = float(args[2])
116  d2 = x * x + y * y
117  if d2 == 0.0:
118  theta = 0.0
119  else:
120  theta = math.degrees(math.atan2(y, x))
121  if theta < 0.0:
122  theta += 360.0
123  if z == 0.0:
124  phi = 0.0
125  else:
126  phi = math.degrees(math.atan2(z, math.sqrt(d2)))
127  return (theta, phi)
128  raise TypeError('Expecting 2 or 3 coordinates, or a tuple/list ' +
129  'containing 2 or 3 coordinates')

Variable Documentation

float lsst.geom.geometry.ANGLE_EPSILON = 0.001

Definition at line 32 of file geometry.py.

float lsst.geom.geometry.ARCSEC_PER_DEG = 3600.0

Definition at line 28 of file geometry.py.

float lsst.geom.geometry.COS_MAX = 1.0

Definition at line 39 of file geometry.py.

int lsst.geom.geometry.CROSS_N2MIN = 2

Definition at line 43 of file geometry.py.

float lsst.geom.geometry.DEG_PER_ARCSEC = 1.0

Definition at line 29 of file geometry.py.

float lsst.geom.geometry.EPSILON = float_info.epsilon

Definition at line 57 of file geometry.py.

tuple lsst.geom.geometry.INF = float('inf')

Definition at line 50 of file geometry.py.

lsst.geom.geometry.isinf = math.isinf

Definition at line 66 of file geometry.py.

float lsst.geom.geometry.MAX_FLOAT = float_info.max

Definition at line 55 of file geometry.py.

float lsst.geom.geometry.MIN_FLOAT = float_info.min

Definition at line 56 of file geometry.py.

tuple lsst.geom.geometry.NEG_INF = float('-inf')

Definition at line 51 of file geometry.py.

float lsst.geom.geometry.POLE_EPSILON = 1.0

Definition at line 35 of file geometry.py.

tuple lsst.geom.geometry.SIN_MIN = math.sqrt(CROSS_N2MIN)

Definition at line 47 of file geometry.py.

Generated on Thu Sep 24 2015 02:29:35 for LSSTApplications by   doxygen 1.8.5