Skip to content

Commit

Permalink
Add wcsNearlyEqualOverBBox free function
Browse files Browse the repository at this point in the history 
Both the new wcsNearlyEqualOverBBox and the existing
assertWcsNearlyEqualOverBBox are thin wrappers
around a new private function _compareWcsOverBBox.
Also add a unit test for wcsNearlyEqualOverBBox.
  • Loading branch information
@r-owen
r-owen committed Aug 5, 2015
1 parent 9c467a9 commit 55d8bb9
Show file tree
Hide file tree
Showing 2 changed files with 114 additions and 27 deletions.
111 changes: 91 additions & 20 deletions python/lsst/afw/image/basicUtils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@

"""Application Framework image-related classes including Image, Mask and MaskedImage
"""
import itertools
import math

import numpy
Expand All @@ -32,7 +33,7 @@
from . import imageLib

__all__ = ["makeImageFromArray", "makeMaskFromArray", "makeMaskedImageFromArrays",
"assertWcsNearlyEqualOverBBox"]
"wcsNearlyEqualOverBBox", "assertWcsNearlyEqualOverBBox"]

suffixes = {str(numpy.uint16): "U", str(numpy.int32): "I", str(numpy.float32): "F", str(numpy.float64): "D"}

Expand All @@ -58,13 +59,10 @@ def makeMaskedImageFromArrays(image, mask=None, variance=None):
cls = getattr(imageLib, "MaskedImage%s" % (suffixes[str(image.dtype.type)],))
return cls(makeImageFromArray(image), makeMaskFromArray(mask), makeImageFromArray(variance))

@lsst.utils.tests.inTestCase
def assertWcsNearlyEqualOverBBox(testCase, wcs0, wcs1, bbox, maxDiffSky=0.01*afwGeom.arcseconds, maxDiffPix=0.01,
nx=5, ny=5, msg="WCSs differ"):
"""Compare pixelToSky and skyToPixel for two WCS over a rectangular grid of pixel positions
def _compareWcsOverBBox(wcs0, wcs1, bbox, maxDiffSky=0.01*afwGeom.arcseconds,
maxDiffPix=0.01, nx=5, ny=5, doShortCircuit=True):
"""!Compare two WCS over a rectangular grid of pixel positions
@param[in] testCase unittest.TestCase instance the test is part of;
an object supporting one method: fail(self, msgStr)
@param[in] wcs0 WCS 0 (an lsst.afw.image.Wcs)
@param[in] wcs1 WCS 1 (an lsst.afw.image.Wcs)
@param[in] bbox boundaries of pixel grid over which to compare the WCSs (an lsst.afw.geom.Box2I or Box2D)
Expand All @@ -73,32 +71,45 @@ def assertWcsNearlyEqualOverBBox(testCase, wcs0, wcs1, bbox, maxDiffSky=0.01*afw
@param[in] maxDiffPix maximum separation between pixel positions computed using Wcs.skyToPixel
@param[in] nx number of points in x for the grid of pixel positions
@param[in] ny number of points in y for the grid of pixel positions
@param[in] msg exception message prefix; details of the error are appended after ": "
@param[in] doShortCircuit if True then stop at the first error, else test all values in the grid
and return information about the worst violations found
@throw AssertionError if the two WCSs do not match sufficiently closely
@return return an empty string if the WCS are sufficiently close; else return a string describing
the largest error measured in pixel coordinates (if sky to pixel error was excessive) and sky coordinates
(if pixel to sky error was excessive). If doShortCircuit is true then the reported error is likely to be
much less than the maximum error across the whole pixel grid.
"""
if nx < 1 or ny < 1:
raise RuntimeError("nx = %s and ny = %s must both be positive" % (nx, ny))
if maxDiffSky <= 0*afwGeom.arcseconds:
raise RuntimeError("maxDiffSky = %s must be positive" % (maxDiffSky,))
if maxDiffPix <= 0:
raise RuntimeError("maxDiffPix = %s must be positive" % (maxDiffPix,))

bboxd = afwGeom.Box2D(bbox)
xList = numpy.linspace(bboxd.getMinX(), bboxd.getMaxX(), nx)
yList = numpy.linspace(bboxd.getMinY(), bboxd.getMaxY(), ny)
measDiffSky = (afwGeom.Angle(0), "?") # (sky diff, pix pos)
measDiffPix = (0, "?") # (pix diff, sky pos)
for x in xList:
for y in yList:
fromPixPos = afwGeom.Point2D(x, y)
sky0 = wcs0.pixelToSky(fromPixPos)
sky1 = wcs1.pixelToSky(fromPixPos)
diffSky = sky0.angularSeparation(sky1)
for x, y in itertools.product(xList, yList):
fromPixPos = afwGeom.Point2D(x, y)
sky0 = wcs0.pixelToSky(fromPixPos)
sky1 = wcs1.pixelToSky(fromPixPos)
diffSky = sky0.angularSeparation(sky1)
if diffSky > maxDiffSky:
if diffSky > measDiffSky[0]:
measDiffSky = (diffSky, fromPixPos)
if doShortCircuit:
break

toPixPos0 = wcs0.skyToPixel(sky0)
toPixPos1 = wcs1.skyToPixel(sky0)
diffPix = math.hypot(*(toPixPos0 - toPixPos1))
toPixPos0 = wcs0.skyToPixel(sky0)
toPixPos1 = wcs1.skyToPixel(sky0)
diffPix = math.hypot(*(toPixPos0 - toPixPos1))
if diffPix > maxDiffPix:
if diffPix > measDiffPix[0]:
measDiffPix = (diffPix, sky0)
if doShortCircuit:
break

msgList = []
if measDiffSky[0] > maxDiffSky:
Expand All @@ -107,5 +118,65 @@ def assertWcsNearlyEqualOverBBox(testCase, wcs0, wcs1, bbox, maxDiffSky=0.01*afw
if measDiffPix[0] > maxDiffPix:
msgList.append("%s max measured pix error > %s max allowed pix error at sky pos=%s" %
(measDiffPix[0], maxDiffPix, measDiffPix[1]))
if msgList:
testCase.fail("%s: %s" % (msg, "; ".join(msgList)))

return "; ".join(msgList)

def wcsNearlyEqualOverBBox(wcs0, wcs1, bbox, maxDiffSky=0.01*afwGeom.arcseconds,
maxDiffPix=0.01, nx=5, ny=5):
"""!Return True if two WCS are nearly equal over a grid of pixel positions, else False
@param[in] wcs0 WCS 0 (an lsst.afw.image.Wcs)
@param[in] wcs1 WCS 1 (an lsst.afw.image.Wcs)
@param[in] bbox boundaries of pixel grid over which to compare the WCSs (an lsst.afw.geom.Box2I or Box2D)
@param[in] maxDiffSky maximum separation between sky positions computed using Wcs.pixelToSky
(an lsst.afw.geom.Angle)
@param[in] maxDiffPix maximum separation between pixel positions computed using Wcs.skyToPixel
@param[in] nx number of points in x for the grid of pixel positions
@param[in] ny number of points in y for the grid of pixel positions
@param[in] doShortCircuit if True then stop at the first error, else test all values in the grid
and return information about the worst violations found
"""
return not bool(_compareWcsOverBBox(
wcs0 = wcs0,
wcs1 = wcs1,
bbox = bbox,
maxDiffSky = maxDiffSky,
maxDiffPix = maxDiffPix,
nx = nx,
ny = ny,
doShortCircuit = True,
))

@lsst.utils.tests.inTestCase
def assertWcsNearlyEqualOverBBox(testCase, wcs0, wcs1, bbox, maxDiffSky=0.01*afwGeom.arcseconds, maxDiffPix=0.01,
nx=5, ny=5, msg="WCSs differ"):
"""!Compare pixelToSky and skyToPixel for two WCS over a rectangular grid of pixel positions
If the WCS are too divergent, call testCase.fail; the message describes the largest error measured
in pixel coordinates (if sky to pixel error was excessive) and sky coordinates (if pixel to sky error
was excessive) across the entire pixel grid.
@param[in] testCase unittest.TestCase instance the test is part of;
an object supporting one method: fail(self, msgStr)
@param[in] wcs0 WCS 0 (an lsst.afw.image.Wcs)
@param[in] wcs1 WCS 1 (an lsst.afw.image.Wcs)
@param[in] bbox boundaries of pixel grid over which to compare the WCSs (an lsst.afw.geom.Box2I or Box2D)
@param[in] maxDiffSky maximum separation between sky positions computed using Wcs.pixelToSky
(an lsst.afw.geom.Angle)
@param[in] maxDiffPix maximum separation between pixel positions computed using Wcs.skyToPixel
@param[in] nx number of points in x for the grid of pixel positions
@param[in] ny number of points in y for the grid of pixel positions
@param[in] msg exception message prefix; details of the error are appended after ": "
"""
errMsg = _compareWcsOverBBox(
wcs0 = wcs0,
wcs1 = wcs1,
bbox = bbox,
maxDiffSky = maxDiffSky,
maxDiffPix = maxDiffPix,
nx = nx,
ny = ny,
doShortCircuit = False,
)
if errMsg:
testCase.fail("%s: %s" % (msg, errMsg))
30 changes: 23 additions & 7 deletions tests/testTestMethods.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
import lsst.afw.coord as afwCoord
import lsst.afw.geom as afwGeom
import lsst.afw.image as afwImage
from lsst.afw.image.basicUtils import _compareWcsOverBBox

class TestTestUtils(utilsTests.TestCase):
"""Test test methods added to lsst.utils.tests.TestCase
Expand Down Expand Up @@ -161,7 +162,7 @@ def testAssertPairsNearlyEqual(self):
pair0, pair1, maxDiff=radialDiff-1e-7)

def testAssertWcssNearlyEqualOverBBox(self):
"""Test assertWcsNearlyEqualOverBBox"""
"""Test assertWcsNearlyEqualOverBBox and wcsNearlyEqualOverBBox"""
bbox = afwGeom.Box2I(afwGeom.Point2I(0, 0), afwGeom.Extent2I(3001, 3001))
ctrPix = afwGeom.Point2I(1500, 1500)
metadata = dafBase.PropertySet()
Expand All @@ -185,15 +186,30 @@ def testAssertWcssNearlyEqualOverBBox(self):

self.assertWcsNearlyEqualOverBBox(wcs0, wcs0, bbox,
maxDiffSky=1e-7*afwGeom.arcseconds, maxDiffPix=1e-7)
self.assertTrue(afwImage.wcsNearlyEqualOverBBox(wcs0, wcs0, bbox,
maxDiffSky=1e-7*afwGeom.arcseconds, maxDiffPix=1e-7))

self.assertWcsNearlyEqualOverBBox(wcs0, wcs1, bbox,
maxDiffSky=0.04*afwGeom.arcseconds, maxDiffPix=0.02)

self.assertRaises(AssertionError, self.assertWcsNearlyEqualOverBBox,
wcs0, wcs1, bbox, maxDiffSky=0.001*afwGeom.arcseconds, maxDiffPix=0.02)

self.assertRaises(AssertionError, self.assertWcsNearlyEqualOverBBox,
wcs0, wcs1, bbox, maxDiffSky=0.04*afwGeom.arcseconds, maxDiffPix=0.001)
self.assertTrue(afwImage.wcsNearlyEqualOverBBox(wcs0, wcs1, bbox,
maxDiffSky=0.04*afwGeom.arcseconds, maxDiffPix=0.02))

self.assertRaises(AssertionError, self.assertWcsNearlyEqualOverBBox, wcs0, wcs1, bbox,
maxDiffSky=0.001*afwGeom.arcseconds, maxDiffPix=0.02)
self.assertFalse(afwImage.wcsNearlyEqualOverBBox(wcs0, wcs1, bbox,
maxDiffSky=0.001*afwGeom.arcseconds, maxDiffPix=0.02))

self.assertRaises(AssertionError, self.assertWcsNearlyEqualOverBBox, wcs0, wcs1, bbox,
maxDiffSky=0.04*afwGeom.arcseconds, maxDiffPix=0.001)
self.assertFalse(afwImage.wcsNearlyEqualOverBBox(wcs0, wcs1, bbox,
maxDiffSky=0.04*afwGeom.arcseconds, maxDiffPix=0.001))

# check that doShortCircuit works in the private implementation
errStr1 = _compareWcsOverBBox(wcs0, wcs1, bbox,
maxDiffSky=0.001*afwGeom.arcseconds, maxDiffPix=0.001, doShortCircuit=False)
errStr2 = _compareWcsOverBBox(wcs0, wcs1, bbox,
maxDiffSky=0.001*afwGeom.arcseconds, maxDiffPix=0.001, doShortCircuit=True)
self.assertNotEqual(errStr1, errStr2)

#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

Expand Down

0 comments on commit 55d8bb9

Please sign in to comment.