LSSTApplications  11.0-13-gbb96280,12.1+18,12.1+7,12.1-1-g14f38d3+72,12.1-1-g16c0db7+5,12.1-1-g5961e7a+84,12.1-1-ge22e12b+23,12.1-11-g06625e2+4,12.1-11-g0d7f63b+4,12.1-19-gd507bfc,12.1-2-g7dda0ab+38,12.1-2-gc0bc6ab+81,12.1-21-g6ffe579+2,12.1-21-gbdb6c2a+4,12.1-24-g941c398+5,12.1-3-g57f6835+7,12.1-3-gf0736f3,12.1-37-g3ddd237,12.1-4-gf46015e+5,12.1-5-g06c326c+20,12.1-5-g648ee80+3,12.1-5-gc2189d7+4,12.1-6-ga608fc0+1,12.1-7-g3349e2a+5,12.1-7-gfd75620+9,12.1-9-g577b946+5,12.1-9-gc4df26a+10
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 1721 of file geometry.py.

1721 
1722 def _feasible1D(points, y):
1723  xmin = NEG_INF
1724  xmax = INF
1725  for v in points:
1726  if abs(v[0]) <= MIN_FLOAT:
1727  if y * v[1] <= 0.0:
1728  return False
1729  # inequality is trivially true, skip it
1730  else:
1731  xlim = - y * v[1] / v[0]
1732  if v[0] > 0.0:
1733  xmin = max(xmin, xlim)
1734  else:
1735  xmax = min(xmax, xlim)
1736  if xmax <= xmin:
1737  return False
1738  return True
1739 
def lsst.geom.geometry._feasible2D (   points,
  z 
)
private

Definition at line 1649 of file geometry.py.

1649 
1650 def _feasible2D(points, z):
1651  I1 = []
1652  I2 = []
1653  xmin = NEG_INF
1654  xmax = INF
1655  # transform each constraint of the form x*v[0] + y*v[1] + z*v[2] > 0
1656  # into y op a*x + b or x op c, where op is < or >
1657  for v in points:
1658  if abs(v[1]) <= MIN_FLOAT:
1659  if abs(v[0]) <= MIN_FLOAT:
1660  if z * v[2] <= 0.0:
1661  # inequalities trivially lack a solution
1662  return False
1663  # current inequality is trivially true, skip it
1664  else:
1665  xlim = - z * v[2] / v[0]
1666  if v[0] > 0.0:
1667  xmin = max(xmin, xlim)
1668  else:
1669  xmax = min(xmax, xlim)
1670  if xmax <= xmin:
1671  # inequalities trivially lack a solution
1672  return False
1673  else:
1674  # finite since |z|, |v[i]| <= 1 and 1/MIN_FLOAT is finite
1675  coeffs = (v[0] / v[1], - z * v[2] / v[1])
1676  if v[1] > 0.0:
1677  I1.append(coeffs)
1678  else:
1679  I2.append(coeffs)
1680  # At this point (xmin, xmax) is non-empty - if either I1 or I2 is empty
1681  # then a solution trivially exists
1682  if len(I1) == 0 or len(I2) == 0:
1683  return True
1684  # Check for a feasible solution to the inequalities I1 of the form
1685  # form y > a*x + b, I2 of the form y < a*x + b, x > xmin and x < xmax
1686  numConstraints = 0
1687  while True:
1688  intersections = _prune(I1, xmin, xmax, operator.gt)
1689  intersections.extend(_prune(I2, xmin, xmax, operator.lt))
1690  if len(intersections) == 0:
1691  # I1 and I2 each contain exactly one constraint
1692  a1, b1 = I1[0]
1693  a2, b2 = I2[0]
1694  if a1 == a2:
1695  xi = INF
1696  else:
1697  xi = (b2 - b1) / (a1 - a2)
1698  if isinf(xi):
1699  return b1 < b2
1700  return (xi > xmin or a1 < a2) and (xi < xmax or a1 > a2)
1701  elif numConstraints == len(I1) + len(I2):
1702  # No constraints were pruned during search interval refinement,
1703  # and g was not conclusively less than h. Conservatively return
1704  # False to avoid an infinite loop.
1705  return False
1706  numConstraints = len(I1) + len(I2)
1707  x, xerr = median(intersections)
1708  # If g(x) < h(x), x is a feasible solution. Otherwise, refine the
1709  # search interval by examining the one-sided derivates of g/h.
1710  gmin, gmax, sg, Sg = _gh(I1, x, xerr, operator.gt)
1711  hmin, hmax, sh, Sh = _gh(I2, x, xerr, operator.lt)
1712  if gmax <= hmin:
1713  return True
1714  elif sg > Sh:
1715  xmax = x + xerr
1716  elif Sg < sh:
1717  xmin = x - xerr
1718  else:
1719  return False
1720 
def lsst.geom.geometry._gh (   constraints,
  x,
  xerr,
  op 
)
private

Definition at line 1630 of file geometry.py.

1630 
1631 def _gh(constraints, x, xerr, op):
1632  a, b = constraints[0]
1633  amin, amax = a, a
1634  vmin, vmax = _vrange(x, xerr, a, b)
1635  for i in range(1, len(constraints)):
1636  a, b = constraints[i]
1637  vimin, vimax = _vrange(x, xerr, a, b)
1638  if vimax < vmin or vimin > vmax:
1639  if op(vimin, vmin):
1640  amin = a
1641  amax = a
1642  vmin = vimin
1643  vmax = vimax
1644  else:
1645  amin = min(amin, a)
1646  amax = max(amax, a)
1647  return vmin, vmax, amin, amax
1648 
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 1525 of file geometry.py.

1525 
1526 def _medianOfMedians(array, left, right):
1527  """Returns the "median of medians" for an array. This primitive is used
1528  for pivot selection in the median finding algorithm.
1529  """
1530  while True:
1531  if right - left <= 5:
1532  return _medianOfN(array, left, right - left)
1533  i = left
1534  j = left
1535  while i + 4 < right:
1536  k = _medianOfN(array, i, 5)
1537  tmp = array[j]
1538  array[j] = array[k]
1539  array[k] = tmp
1540  j += 1
1541  i += 5
1542  right = j
1543 
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 1502 of file geometry.py.

1502 
1503 def _medianOfN(array, i, n):
1504  """Finds the median of n consecutive elements in an array starting
1505  at index i (efficient for small n). The index of the median element
1506  is returned.
1507  """
1508  if n == 1:
1509  return i
1510  k = n >> 1
1511  e1 = i + k + 1
1512  e2 = i + n
1513  for j in range(i, e1):
1514  minIndex = j
1515  minValue = array[j]
1516  for s in range(j + 1, e2):
1517  if array[s] < minValue:
1518  minIndex = s
1519  minValue = array[s]
1520  tmp = array[j]
1521  array[j] = array[minIndex]
1522  array[minIndex] = tmp
1523  return i + k
1524 
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 1484 of file geometry.py.

1484 
1485 def _partition(array, left, right, i):
1486  """Partitions an array around the pivot value at index i,
1487  returning the new index of the pivot value.
1488  """
1489  pivot = array[i]
1490  array[i] = array[right - 1]
1491  j = left
1492  for k in range(left, right - 1):
1493  if array[k] < pivot:
1494  tmp = array[j]
1495  array[j] = array[k]
1496  array[k] = tmp
1497  j += 1
1498  array[right - 1] = array[j]
1499  array[j] = pivot
1500  return j
1501 
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 1573 of file geometry.py.

1573 
1574 def _prune(constraints, xmin, xmax, op):
1575  """Removes redundant constraints and reports intersection points
1576  of non-redundant pairs. The constraint list is modified in place.
1577  """
1578  intersections = []
1579  i = 0
1580  while i < len(constraints) - 1:
1581  a1, b1 = constraints[i]
1582  a2, b2 = constraints[i + 1]
1583  # if constraints are near parallel or intersect to the left or right
1584  # of the feasible x range, remove the higher/lower constraint
1585  da = a1 - a2
1586  if abs(da) < MIN_FLOAT / EPSILON:
1587  xi = INF
1588  else:
1589  xi = (b2 - b1) / da
1590  xierr = 2.0 * EPSILON * abs(xi)
1591  if isinf(xi):
1592  if op(b1, b2):
1593  constraints[i + 1] = constraints[-1]
1594  else:
1595  constraints[i] = constraints[-1]
1596  constraints.pop()
1597  else:
1598  if xi + xierr <= xmin:
1599  if op(a1, a2):
1600  constraints[i + 1] = constraints[-1]
1601  else:
1602  constraints[i] = constraints[-1]
1603  constraints.pop()
1604  elif xi - xierr >= xmax:
1605  if op(a1, a2):
1606  constraints[i] = constraints[-1]
1607  else:
1608  constraints[i + 1] = constraints[-1]
1609  constraints.pop()
1610  else:
1611  # save intersection for later
1612  intersections.append((xi, xierr))
1613  i += 2
1614  return intersections
1615 
def lsst.geom.geometry._vrange (   x,
  xerr,
  a,
  b 
)
private

Definition at line 1616 of file geometry.py.

1616 
1617 def _vrange(x, xerr, a, b):
1618  p = a * (x - xerr)
1619  v = p + b
1620  verr = EPSILON * abs(p) + EPSILON * abs(v)
1621  vmin = v - verr
1622  vmax = v + verr
1623  p = a * (x + xerr)
1624  v = p + b
1625  verr = EPSILON * abs(p) + EPSILON * abs(v)
1626  vmin = min(vmin, v - verr)
1627  vmax = max(vmax, v + verr)
1628  return vmin, vmax
1629 
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 222 of file geometry.py.

222 
223 def alpha(r, centerPhi, phi):
224  """Returns the longitude angle alpha of the intersections (alpha, phi),
225  (-alpha, phi) of the circle centered on (0, centerPhi) with radius r and
226  the plane z = sin(phi). If there is no intersection, None is returned.
227  """
228  if phi < centerPhi - r or phi > centerPhi + r:
229  return None
230  a = abs(centerPhi - phi)
231  if a <= r - (90.0 - phi) or a <= r - (90.0 + phi):
232  return None
233  p = math.radians(phi)
234  cp = math.radians(centerPhi)
235  x = math.cos(math.radians(r)) - math.sin(cp) * math.sin(cp)
236  u = math.cos(cp) * math.cos(p)
237  y = math.sqrt(abs(u * u - x * x))
238  return math.degrees(abs(math.atan2(y, x)))
239 
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 329 of file geometry.py.

329 
330 def between(p, n, v1, v2):
331  """Tests whether p lies on the shortest great circle arc from cartesian
332  unit vectors v1 to v2, assuming that p is a unit vector on the plane
333  defined by the origin, v1 and v2.
334  """
335  p1 = cross(n, v1)
336  p2 = cross(n, v2)
337  return dot(p1, p) >= 0.0 and dot(p2, p) <= 0.0
338 
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 260 of file geometry.py.

260 
261 def cartesianAngularSep(v1, v2):
262  """Returns the angular separation in degrees between points v1 and v2,
263  which must both be specified as cartesian 3-vectors.
264  """
265  cs = dot(v1, v2)
266  n = cross(v1, v2)
267  ss = math.sqrt(dot(n, n))
268  if cs == 0.0 and ss == 0.0:
269  return 0.0
270  return math.degrees(math.atan2(ss, cs))
271 
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 141 of file geometry.py.

141 
142 def cartesianUnitVector(*args):
143  """Returns a unit cartesian 3-vector corresponding to the input
144  coordinates, which can be spherical or 3D cartesian. The 2 (spherical)
145  or 3 (cartesian 3-vector) input coordinates can either be passed
146  individually or as a tuple/list, and can be of any type convertible
147  to a float.
148  """
149  if len(args) == 1:
150  args = args[0]
151  if len(args) == 2:
152  theta = math.radians(float(args[0]))
153  phi = math.radians(float(args[1]))
154  cosPhi = math.cos(phi)
155  return (math.cos(theta) * cosPhi,
156  math.sin(theta) * cosPhi,
157  math.sin(phi))
158  elif len(args) == 3:
159  x = float(args[0])
160  y = float(args[1])
161  z = float(args[2])
162  n = math.sqrt(x * x + y * y + z * z)
163  if n == 0.0:
164  raise RuntimeError('Cannot normalize a 3-vector with 0 magnitude')
165  return (x / n, y / n, z / n)
166  raise TypeError('Expecting 2 or 3 coordinates, or a tuple/list ' +
167  'containing 2 or 3 coordinates')
168 
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 316 of file geometry.py.

316 
317 def centroid(vertices):
318  """Computes the centroid of a set of vertices (which must be passed in
319  as a list of cartesian unit vectors) on the unit sphere.
320  """
321  x, y, z = 0.0, 0.0, 0.0
322  # Note: could use more accurate summation routines
323  for v in vertices:
324  x += v[0]
325  y += v[1]
326  z += v[2]
327  return normalize((x, y, z))
328 
def lsst.geom.geometry.clampPhi (   phi) 
Clamps the input latitude angle phi to [-90.0, 90.0] deg.

Definition at line 184 of file geometry.py.

184 
185 def clampPhi(phi):
186  """Clamps the input latitude angle phi to [-90.0, 90.0] deg.
187  """
188  if phi < -90.0:
189  return -90.0
190  elif phi >= 90.0:
191  return 90.0
192  return phi
193 
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 1782 of file geometry.py.

1782 
1783 def convex(vertices):
1784  """Tests whether an ordered list of vertices (which must be specified
1785  as cartesian unit vectors) form a spherical convex polygon and determines
1786  their winding order. Returns a 2-tuple as follows:
1787 
1788  (True, True): The vertex list forms a spherical convex polygon and is in
1789  counter-clockwise order when viewed from outside the unit
1790  sphere in a right handed coordinate system.
1791  (True, False): The vertex list forms a spherical convex polygon and is in
1792  in clockwise order.
1793  (False, msg): The vertex list does not form a spherical convex polygon -
1794  msg is a string describing why the test failed.
1795 
1796  The algorithm completes in O(n) time, where n is the number of
1797  input vertices.
1798  """
1799  if len(vertices) < 3:
1800  return (False, '3 or more vertices must be specified')
1801  if not hemispherical(vertices):
1802  return (False, 'vertices are not hemispherical')
1803  center = centroid(vertices)
1804  windingAngle = 0.0
1805  clockwise = False
1806  counterClockwise = False
1807  p1 = cross(center, vertices[-1])
1808  n2 = dot(p1, p1)
1809  if abs(n2) < CROSS_N2MIN:
1810  return (False, 'centroid of vertices is too close to a vertex')
1811  for i in range(len(vertices)):
1812  beg = vertices[i - 2]
1813  mid = vertices[i - 1]
1814  end = vertices[i]
1815  plane = cross(mid, end)
1816  n2 = dot(plane, plane)
1817  if dot(mid, end) >= COS_MAX or n2 < CROSS_N2MIN:
1818  return (False, 'vertex list contains [near-]duplicates')
1819  plane = invScale(plane, math.sqrt(n2))
1820  d = dot(plane, beg)
1821  if d > SIN_MIN:
1822  if clockwise:
1823  return (False, 'vertices wind around centroid in both ' +
1824  'clockwise and counter-clockwise order')
1825  counterClockwise = True
1826  elif d < -SIN_MIN:
1827  if counterClockwise:
1828  return (False, 'vertices wind around centroid in both ' +
1829  'clockwise and counter-clockwise order')
1830  clockwise = True
1831  # center must be inside polygon formed by vertices if they form a
1832  # convex polygon - an equivalent check is to test that polygon
1833  # vertices always wind around center in the same direction.
1834  d = dot(plane, center)
1835  if d < SIN_MIN and counterClockwise or d > -SIN_MIN and clockwise:
1836  return (False, 'centroid of vertices is not conclusively ' +
1837  'inside all edges')
1838  # sum up winding angle for edge (mid, end)
1839  p2 = cross(center, end)
1840  n2 = dot(p2, p2)
1841  if abs(n2) < CROSS_N2MIN:
1842  return (False, 'centroid of vertices is too close to a vertex')
1843  windingAngle += cartesianAngularSep(p1, p2)
1844  p1 = p2
1845  # For convex polygons, the closest multiple of 360 to
1846  # windingAngle is 1.
1847  m = math.floor(windingAngle / 360.0)
1848  if m == 0.0 and windingAngle > 180.0 or m == 1.0 and windingAngle < 540.0:
1849  return (True, counterClockwise)
1850  return (False, 'vertices do not completely wind around centroid, or ' +
1851  'wind around it multiple times')
1852 
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 1853 of file geometry.py.

1853 
1854 def convexHull(points):
1855  """Computes the convex hull (a spherical convex polygon) of an unordered
1856  list of points on the unit sphere, which must be passed in as cartesian
1857  unit vectors. The algorithm takes O(n log n) time, where n is the number
1858  of points.
1859  """
1860  if len(points) < 3:
1861  return None
1862  if not hemispherical(points):
1863  return None
1864  center = centroid(points)
1865  maxSep = 0.0
1866  extremum = 0
1867  # Vertex furthest from the center is on the hull
1868  for i in range(len(points)):
1869  sep = cartesianAngularSep(points[i], center)
1870  if sep > maxSep:
1871  extremum = i
1872  maxSep = sep
1873  anchor = points[extremum]
1874  refPlane = cross(center, anchor)
1875  n2 = dot(refPlane, refPlane)
1876  if n2 < CROSS_N2MIN:
1877  # vertex and center are too close
1878  return None
1879  refPlane = invScale(refPlane, math.sqrt(n2))
1880  # Order points by winding angle from the first (extreme) vertex
1881  av = [(0.0, anchor)]
1882  for i in range(0, len(points)):
1883  if i == extremum:
1884  continue
1885  v = points[i]
1886  plane = cross(center, v)
1887  n2 = dot(plane, plane)
1888  if n2 >= CROSS_N2MIN:
1889  plane = invScale(plane, math.sqrt(n2))
1890  p = cross(refPlane, plane)
1891  sa = math.sqrt(dot(p, p))
1892  if dot(p, center) < 0.0:
1893  sa = -sa
1894  angle = math.atan2(sa, dot(refPlane, plane))
1895  if angle < 0.0:
1896  angle += 2.0 * math.pi
1897  av.append((angle, v))
1898  if len(av) < 3:
1899  return None
1900  av.sort(key=lambda t: t[0]) # stable, so av[0] still contains anchor
1901  # Loop over vertices using a Graham scan adapted for spherical geometry.
1902  # Chan's algorithm would be an improvement, but seems likely to be slower
1903  # for small vertex counts (the expected case).
1904  verts = [anchor]
1905  edges = [None]
1906  edge = None
1907  i = 1
1908  while i < len(av):
1909  v = av[i][1]
1910  p = cross(anchor, v)
1911  n2 = dot(p, p)
1912  if dot(anchor, v) < COS_MAX and n2 >= CROSS_N2MIN:
1913  if edge is None:
1914  # Compute first edge
1915  edge = invScale(p, math.sqrt(n2))
1916  verts.append(v)
1917  edges.append(edge)
1918  anchor = v
1919  else:
1920  d = dot(v, edge)
1921  if d > SIN_MIN:
1922  # v is inside the edge defined by the last
1923  # 2 vertices on the hull
1924  edge = invScale(p, math.sqrt(n2))
1925  verts.append(v)
1926  edges.append(edge)
1927  anchor = v
1928  elif d < -SIN_MIN:
1929  # backtrack - the most recently added hull vertex
1930  # is not actually on the hull.
1931  verts.pop()
1932  edges.pop()
1933  anchor = verts[-1]
1934  edge = edges[-1]
1935  # reprocess v to decide whether to add it to the hull
1936  # or whether further backtracking is necessary
1937  continue
1938  # v is coplanar with edge, skip it
1939  i += 1
1940  # Handle backtracking necessary for last edge
1941  if len(verts) < 3:
1942  return None
1943  v = verts[0]
1944  while True:
1945  p = cross(anchor, v)
1946  n2 = dot(p, p)
1947  if dot(anchor, v) < COS_MAX and n2 >= CROSS_N2MIN:
1948  if dot(v, edge) > SIN_MIN:
1949  edges[0] = invScale(p, math.sqrt(n2))
1950  break
1951  verts.pop()
1952  edges.pop()
1953  anchor = verts[-1]
1954  edge = edges[-1]
1955  if len(verts) < 3:
1956  return None
1957  return SphericalConvexPolygon(verts, edges)
1958 
1959 
1960 # -- Generic unit sphere partitioning scheme base class ----
1961 
# 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 85 of file geometry.py.

85 
86 def cross(v1, v2):
87  """Returns the cross product of cartesian 3-vectors v1 and v2.
88  """
89  return (v1[1] * v2[2] - v1[2] * v2[1],
90  v1[2] * v2[0] - v1[0] * v2[2],
91  v1[0] * v2[1] - v1[1] * v2[0])
92 
def lsst.geom.geometry.dot (   v1,
  v2 
) 
Returns the dot product of cartesian 3-vectors v1 and v2.

Definition at line 79 of file geometry.py.

79 
80 def dot(v1, v2):
81  """Returns the dot product of cartesian 3-vectors v1 and v2.
82  """
83  return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]
84 
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 1740 of file geometry.py.

1740 
1741 def hemispherical(points):
1742  """Tests whether a set of points is hemispherical, i.e. whether a plane
1743  exists such that all points are strictly on one side of the plane. The
1744  algorithm used is Megiddo's algorithm for linear programming in R2 and
1745  has run-time O(n), where n is the number of points. Points must be passed
1746  in as a list of cartesian unit vectors.
1747  """
1748  # Check whether the set of linear equations
1749  # x * v[0] + y * v[1] + z * v[2] > 0.0 (for v in points)
1750  # has a solution (x, y, z). If (x,y,z) is a solution (is feasible),
1751  # so is C*(x,y,z), C > 0. Therefore we can fix z to 1, -1 and
1752  # perform 2D feasibility tests.
1753  if _feasible2D(points, 1.0):
1754  return True
1755  if _feasible2D(points, -1.0):
1756  return True
1757  # At this point a feasible solution must have z = 0. Fix y to 1, -1 and
1758  # perform 1D feasibility tests.
1759  if _feasible1D(points, 1.0):
1760  return True
1761  if _feasible1D(points, -1.0):
1762  return True
1763  # At this point a feasible solution must have y = z = 0. If all v[0]
1764  # are non-zero and have the same sign, then there is a feasible solution.
1765  havePos = False
1766  haveNeg = False
1767  for v in points:
1768  if v[0] > 0.0:
1769  if haveNeg:
1770  return False
1771  havePos = True
1772  elif v[0] < 0.0:
1773  if havePos:
1774  return False
1775  haveNeg = True
1776  else:
1777  return False
1778  return True
1779 
1780 
1781 # -- 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 93 of file geometry.py.

93 
94 def invScale(v, s):
95  """Returns a copy of the cartesian 3-vector v scaled by 1 / s.
96  """
97  return (v[0] / s, v[1] / s, v[2] / s)
98 
def lsst.geom.geometry.isinf (   x)

Definition at line 73 of file geometry.py.

73 
74  def isinf(x):
75  return x == INF or x == NEG_INF
76 
77 # -- Utility methods ----
78 
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 240 of file geometry.py.

240 
241 def maxAlpha(r, centerPhi):
242  """Computes alpha, the extent in longitude angle [-alpha, alpha] of the
243  circle with radius r and center (0, centerPhi) on the unit sphere.
244  Both r and centerPhi are assumed to be in units of degrees.
245  centerPhi is clamped to lie in the range [-90,90] and r must
246  lie in the range [0, 90].
247  """
248  assert r >= 0.0 and r <= 90.0
249  if r == 0.0:
250  return 0.0
251  centerPhi = clampPhi(centerPhi)
252  if abs(centerPhi) + r > 90.0 - POLE_EPSILON:
253  return 180.0
254  r = math.radians(r)
255  c = math.radians(centerPhi)
256  y = math.sin(r)
257  x = math.sqrt(abs(math.cos(c - r) * math.cos(c + r)))
258  return math.degrees(abs(math.atan(y / x)))
259 
def lsst.geom.geometry.median (   array) 
Finds the median element of the given array in linear time.
Examples:
imageStatistics.cc.

Definition at line 1544 of file geometry.py.

1544 
1545 def median(array):
1546  """Finds the median element of the given array in linear time.
1547  """
1548  left = 0
1549  right = len(array)
1550  if right == 0:
1551  return None
1552  k = right >> 1
1553  while True:
1554  i = _medianOfMedians(array, left, right)
1555  i = _partition(array, left, right, i)
1556  if k == i:
1557  return array[k]
1558  elif k < i:
1559  right = i
1560  else:
1561  left = i + 1
1562 
1563 
1564 # -- Testing whether a set of points is hemispherical ----
1565 #
1566 # This test is used by the convexity test and convex hull algorithm. It is
1567 # implemented using Megiddo's algorithm for linear programming in R2, see:
1568 #
1569 # Megiddo, N. 1982. Linear-time algorithms for linear programming in R3 and related problems.
1570 # In Proceedings of the 23rd Annual Symposium on Foundations of Computer Science (November 03 - 05, 1982).
1571 # SFCS. IEEE Computer Society, Washington, DC, 329-338. DOI=
1572 # 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 272 of file geometry.py.

272 
273 def minEdgeSep(p, n, v1, v2):
274  """Returns the minimum angular separation in degrees between p
275  and points on the great circle edge with plane normal n and
276  vertices v1, v2. All inputs must be unit cartesian 3-vectors.
277  """
278  p1 = cross(n, v1)
279  p2 = cross(n, v2)
280  if dot(p1, p) >= 0.0 and dot(p2, p) <= 0.0:
281  return abs(90.0 - cartesianAngularSep(p, n))
282  else:
283  return min(cartesianAngularSep(p, v1), cartesianAngularSep(p, v2))
284 
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 285 of file geometry.py.

285 
286 def minPhiEdgeSep(p, phi, minTheta, maxTheta):
287  """Returns the minimum angular separation in degrees between p
288  and points on the small circle edge with constant latitude angle
289  phi and vertices (minTheta, phi), (maxTheta, phi). p must be in
290  spherical coordinates.
291  """
292  if minTheta > maxTheta:
293  if p[0] >= minTheta or p[0] <= maxTheta:
294  return abs(p[1] - phi)
295  else:
296  if p[0] >= minTheta and p[0] <= maxTheta:
297  return abs(p[1] - phi)
298  return min(sphericalAngularSep(p, (minTheta, phi)),
299  sphericalAngularSep(p, (maxTheta, phi)))
300 
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 301 of file geometry.py.

301 
302 def minThetaEdgeSep(p, theta, minPhi, maxPhi):
303  """Returns the minimum angular separation in degrees between p
304  and points on the great circle edge with constant longitude angle
305  theta and vertices (theta, minPhi), (theta, maxPhi). p must be a
306  unit cartesian 3-vector.
307  """
308  v1 = cartesianUnitVector(theta, minPhi)
309  v2 = cartesianUnitVector(theta, maxPhi)
310  n = cross(v1, v2)
311  d2 = dot(n, n)
312  if d2 == 0.0:
313  return min(cartesianAngularSep(p, v1), cartesianAngularSep(p, v2))
314  return minEdgeSep(p, invScale(n, math.sqrt(d2)), v1, v2)
315 
def lsst.geom.geometry.normalize (   v) 
Returns a normalized copy of the cartesian 3-vector v.

Definition at line 99 of file geometry.py.

99 
100 def normalize(v):
101  """Returns a normalized copy of the cartesian 3-vector v.
102  """
103  n = math.sqrt(dot(v, v))
104  if n == 0.0:
105  raise RuntimeError('Cannot normalize a 3-vector with 0 magnitude')
106  return (v[0] / n, v[1] / n, v[2] / n)
107 
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 205 of file geometry.py.

205 
206 def northEast(v):
207  """Returns unit N,E basis vectors for a point v, which must be a
208  cartesian 3-vector.
209  """
210  north = (-v[0] * v[2],
211  -v[1] * v[2],
212  v[0] * v[0] + v[1] * v[1])
213  if north == (0.0, 0.0, 0.0):
214  # pick an arbitrary orthogonal basis with z = 0
215  north = (-1.0, 0.0, 0.0)
216  east = (0.0, 1.0, 0.0)
217  else:
218  north = normalize(north)
219  east = normalize(cross(north, v))
220  return north, east
221 
def lsst.geom.geometry.reduceTheta (   theta) 
Range reduces the given longitude angle to lie in the range
[0.0, 360.0).

Definition at line 194 of file geometry.py.

194 
195 def reduceTheta(theta):
196  """Range reduces the given longitude angle to lie in the range
197  [0.0, 360.0).
198  """
199  theta = math.fmod(theta, 360.0)
200  if theta < 0.0:
201  return theta + 360.0
202  else:
203  return theta
204 
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 339 of file geometry.py.

339 
340 def segments(phiMin, phiMax, width):
341  """Computes the number of segments to divide the given latitude angle
342  range [phiMin, phiMax] (degrees) into. Two points within the range
343  separated by at least one segment are guaranteed to have angular
344  separation of at least width degrees.
345  """
346  p = max(abs(phiMin), abs(phiMax))
347  if p > 90.0 - 1.0 * DEG_PER_ARCSEC:
348  return 1
349  if width >= 180.0:
350  return 1
351  elif width < 1.0 * DEG_PER_ARCSEC:
352  width = 1.0 * DEG_PER_ARCSEC
353  p = math.radians(p)
354  cw = math.cos(math.radians(width))
355  sp = math.sin(p)
356  cp = math.cos(p)
357  x = cw - sp * sp
358  u = cp * cp
359  y = math.sqrt(abs(u * u - x * x))
360  return int(math.floor((2.0 * math.pi) / abs(math.atan2(y, x))))
361 
362 
363 # -- Regions on the unit sphere ----
364 
# 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 169 of file geometry.py.

169 
170 def sphericalAngularSep(p1, p2):
171  """Returns the angular separation in degrees between points p1 and p2,
172  which must both be specified in spherical coordinates. The implementation
173  uses the halversine distance formula.
174  """
175  sdt = math.sin(math.radians(p1[0] - p2[0]) * 0.5)
176  sdp = math.sin(math.radians(p1[1] - p2[1]) * 0.5)
177  cc = math.cos(math.radians(p1[1])) * math.cos(math.radians(p2[1]))
178  s = math.sqrt(sdp * sdp + cc * sdt * sdt)
179  if s > 1.0:
180  return 180.0
181  else:
182  return 2.0 * math.degrees(math.asin(s))
183 
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 108 of file geometry.py.

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

Variable Documentation

float lsst.geom.geometry.ANGLE_EPSILON = 0.001

Definition at line 37 of file geometry.py.

float lsst.geom.geometry.ARCSEC_PER_DEG = 3600.0

Definition at line 33 of file geometry.py.

float lsst.geom.geometry.COS_MAX = 1.0

Definition at line 44 of file geometry.py.

int lsst.geom.geometry.CROSS_N2MIN = 2

Definition at line 48 of file geometry.py.

float lsst.geom.geometry.DEG_PER_ARCSEC = 1.0

Definition at line 34 of file geometry.py.

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

Definition at line 62 of file geometry.py.

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

Definition at line 55 of file geometry.py.

lsst.geom.geometry.isinf = math.isinf

Definition at line 71 of file geometry.py.

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

Definition at line 60 of file geometry.py.

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

Definition at line 61 of file geometry.py.

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

Definition at line 56 of file geometry.py.

float lsst.geom.geometry.POLE_EPSILON = 1.0

Definition at line 40 of file geometry.py.

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

Definition at line 52 of file geometry.py.

Generated on Mon Feb 27 2017 03:55:46 for LSSTApplications by   doxygen 1.8.5