Coverage for tests / test_matchPessimisticB.py: 16%
153 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 09:11 +0000
1# This file is part of meas_astrom.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22import math
23import os
24import unittest
26import numpy as np
28import lsst.geom
29import lsst.afw.geom as afwGeom
30import lsst.afw.table as afwTable
31import lsst.utils.tests
32from lsst.meas.algorithms import convertReferenceCatalog
33from lsst.meas.astrom.sip import genDistortedImage
34import lsst.meas.astrom as measAstrom
37class TestMatchPessimisticB(unittest.TestCase):
39 def setUp(self):
40 np.random.seed(12345)
42 self.config = measAstrom.MatchPessimisticBTask.ConfigClass()
43 # Value below is to assure all matches are selected. The
44 # original test is set for a 3 arcsecond max match distance
45 # using matchOptimisticB.
46 self.config.minMatchDistPixels = 2.0
47 self.MatchPessimisticB = measAstrom.MatchPessimisticBTask(
48 config=self.config)
50 self.wcs = afwGeom.makeSkyWcs(crpix=lsst.geom.Point2D(791.4, 559.7),
51 crval=lsst.geom.SpherePoint(36.930640, -4.939560, lsst.geom.degrees),
52 cdMatrix=afwGeom.makeCdMatrix(scale=5.17e-5*lsst.geom.degrees))
53 self.distortedWcs = self.wcs
55 self.filename = os.path.join(os.path.dirname(__file__), "cat.xy.fits")
56 self.tolArcsec = .4
57 self.tolPixel = .1
59 # 3 of the objects are removed by the source selector and are used in
60 # matching hence the 183 number vs the total of 186. This is also why
61 # these three objects are missing in the testReferenceFilter test.
62 self.expectedMatches = 183
64 def testLinearXDistort(self):
65 self.singleTestInstance(self.filename, genDistortedImage.linearXDistort)
67 def testLinearYDistort(self):
68 self.singleTestInstance(self.filename, genDistortedImage.linearYDistort)
70 def testQuadraticDistort(self):
71 self.singleTestInstance(self.filename, genDistortedImage.quadraticDistort)
73 def testLargeDistortion(self):
74 # This transform is about as extreme as I can get:
75 # using 0.0005 in the last value appears to produce numerical issues.
77 # It produces a maximum deviation of 459 pixels, which should be
78 # sufficient.
79 pixelsToTanPixels = afwGeom.makeRadialTransform([0.0, 1.1, 0.0004])
80 self.distortedWcs = afwGeom.makeModifiedWcs(pixelTransform=pixelsToTanPixels,
81 wcs=self.wcs,
82 modifyActualPixels=False)
84 def applyDistortion(src):
85 out = src.table.copyRecord(src)
86 out.set(out.table.getCentroidSlot().getMeasKey(),
87 pixelsToTanPixels.applyInverse(src.getCentroid()))
88 return out
90 self.singleTestInstance(self.filename, applyDistortion)
92 def singleTestInstance(self, filename, distortFunc, doPlot=False):
93 sourceCat = self.loadSourceCatalog(self.filename)
94 refCat = self.computePosRefCatalog(sourceCat)
96 # Apply source selector to sourceCat, using the astrometry config defaults
97 tempConfig = measAstrom.AstrometryTask.ConfigClass()
98 # This field isn't in the old test catalog.
99 tempConfig.sourceSelector["science"].flags.bad.remove("base_PixelFlags_flag_nodata")
100 tempConfig.sourceSelector["science"].doCentroidErrorLimit = False
101 tempSolver = measAstrom.AstrometryTask(config=tempConfig, refObjLoader=None)
102 sourceSelection = tempSolver.sourceSelector.run(sourceCat)
104 distortedCat = genDistortedImage.distortList(sourceSelection.sourceCat, distortFunc)
106 if doPlot:
107 import matplotlib.pyplot as plt
109 undistorted = [self.wcs.skyToPixel(self.distortedWcs.pixelToSky(ss.getCentroid()))
110 for ss in distortedCat]
111 refs = [self.wcs.skyToPixel(ss.getCoord()) for ss in refCat]
113 def plot(catalog, symbol):
114 plt.plot([ss.getX() for ss in catalog],
115 [ss.getY() for ss in catalog], symbol)
117 plot(distortedCat, 'b+') # Distorted positions: blue +
118 plot(undistorted, 'g+') # Undistorted positions: green +
119 plot(refs, 'rx') # Reference catalog: red x
120 # The green + should overlap with the red x, because that's how
121 # MatchPessimisticB does it.
123 plt.show()
125 sourceCat = distortedCat
127 matchRes = self.MatchPessimisticB.matchObjectsToSources(
128 refCat=refCat,
129 sourceCat=sourceCat,
130 wcs=self.distortedWcs,
131 sourceFluxField='slot_ApFlux_instFlux',
132 refFluxField="r_flux",
133 )
134 matches = matchRes.matches
135 if doPlot:
136 measAstrom.plotAstrometry(matches=matches, refCat=refCat,
137 sourceCat=sourceCat)
138 self.assertEqual(len(matches), self.expectedMatches)
140 refCoordKey = afwTable.CoordKey(refCat.schema["coord"])
141 srcCoordKey = afwTable.CoordKey(sourceCat.schema["coord"])
142 refCentroidKey = afwTable.Point2DKey(refCat.getSchema()["centroid"])
143 maxDistErr = 0*lsst.geom.radians
145 for refObj, source, distRad in matches:
146 sourceCoord = source.get(srcCoordKey)
147 refCoord = refObj.get(refCoordKey)
148 predDist = sourceCoord.separation(refCoord)
149 distErr = abs(predDist - distRad*lsst.geom.radians)
150 maxDistErr = max(distErr, maxDistErr)
152 if refObj.getId() != source.getId():
153 refCentroid = refObj.get(refCentroidKey)
154 sourceCentroid = source.getCentroid()
155 radius = math.hypot(*(refCentroid - sourceCentroid))
156 self.fail(
157 "ID mismatch: %s at %s != %s at %s; error = %0.1f pix" %
158 (refObj.getId(), refCentroid, source.getId(),
159 sourceCentroid, radius))
161 self.assertLess(maxDistErr.asArcseconds(), 1e-7)
163 def testPassingMatcherState(self):
164 """Test that results of the matcher can be propagated to to in
165 subsequent iterations.
166 """
167 sourceCat = self.loadSourceCatalog(self.filename)
168 refCat = self.computePosRefCatalog(sourceCat)
170 # Apply source selector to sourceCat, using the astrometry config defaults
171 tempConfig = measAstrom.AstrometryTask.ConfigClass()
172 # This field isn't in the old test catalog.
173 tempConfig.sourceSelector["science"].flags.bad.remove("base_PixelFlags_flag_nodata")
174 tempConfig.sourceSelector["science"].doCentroidErrorLimit = False
175 tempSolver = measAstrom.AstrometryTask(config=tempConfig, refObjLoader=None)
176 sourceSelection = tempSolver.sourceSelector.run(sourceCat)
178 distortedCat = genDistortedImage.distortList(sourceSelection.sourceCat,
179 genDistortedImage.linearXDistort)
181 sourceCat = distortedCat
183 matchRes = self.MatchPessimisticB.matchObjectsToSources(
184 refCat=refCat,
185 sourceCat=sourceCat,
186 wcs=self.distortedWcs,
187 sourceFluxField='slot_ApFlux_instFlux',
188 refFluxField="r_flux",
189 )
191 maxShift = matchRes.matchTolerance.maxShift * 300
192 # Force the matcher to use a different pattern thatn the previous
193 # "iteration".
194 matchTol = measAstrom.MatchTolerancePessimistic(
195 maxMatchDist=matchRes.matchTolerance.maxMatchDist,
196 autoMaxMatchDist=matchRes.matchTolerance.autoMaxMatchDist,
197 maxShift=maxShift,
198 lastMatchedPattern=0,
199 failedPatternList=[0],
200 PPMbObj=matchRes.matchTolerance.PPMbObj,
201 )
203 matchRes = self.MatchPessimisticB.matchObjectsToSources(
204 refCat=refCat,
205 sourceCat=sourceCat,
206 wcs=self.distortedWcs,
207 sourceFluxField='slot_ApFlux_instFlux',
208 refFluxField="r_flux",
209 matchTolerance=matchTol,
210 )
212 self.assertEqual(len(matchRes.matches), self.expectedMatches)
213 self.assertLess(matchRes.matchTolerance.maxShift, maxShift)
214 self.assertEqual(matchRes.matchTolerance.lastMatchedPattern, 1)
215 self.assertIsNotNone(matchRes.matchTolerance.maxMatchDist)
216 self.assertIsNotNone(matchRes.matchTolerance.autoMaxMatchDist)
217 self.assertIsNotNone(matchRes.matchTolerance.lastMatchedPattern)
218 self.assertIsNotNone(matchRes.matchTolerance.failedPatternList)
219 self.assertIsNotNone(matchRes.matchTolerance.PPMbObj)
221 def testReferenceFilter(self):
222 """Test sub-selecting reference objects by flux."""
223 sourceCat = self.loadSourceCatalog(self.filename)
224 refCat = self.computePosRefCatalog(sourceCat)
226 # Apply source selector to sourceCat, using the astrometry config defaults
227 tempConfig = measAstrom.AstrometryTask.ConfigClass()
228 # This field isn't in the old test catalog.
229 tempConfig.sourceSelector["science"].flags.bad.remove("base_PixelFlags_flag_nodata")
230 tempConfig.sourceSelector["science"].doCentroidErrorLimit = False
231 tempSolver = measAstrom.AstrometryTask(config=tempConfig, refObjLoader=None)
232 sourceSelection = tempSolver.sourceSelector.run(sourceCat)
234 distortedCat = genDistortedImage.distortList(sourceSelection.sourceCat,
235 genDistortedImage.linearXDistort)
237 matchPessConfig = measAstrom.MatchPessimisticBTask.ConfigClass()
238 matchPessConfig.maxRefObjects = 150
239 matchPessConfig.minMatchDistPixels = 5.0
241 matchPess = measAstrom.MatchPessimisticBTask(config=matchPessConfig)
242 trimmedRefCat = matchPess._filterRefCat(refCat, 'r_flux')
243 self.assertEqual(len(trimmedRefCat), matchPessConfig.maxRefObjects)
245 matchRes = matchPess.matchObjectsToSources(
246 refCat=refCat,
247 sourceCat=distortedCat,
248 wcs=self.distortedWcs,
249 sourceFluxField='slot_ApFlux_instFlux',
250 refFluxField="r_flux",
251 )
253 self.assertEqual(len(matchRes.matches), matchPessConfig.maxRefObjects - 3)
255 def computePosRefCatalog(self, sourceCat):
256 """Generate a position reference catalog from a source catalog
257 """
258 minimalPosRefSchema = convertReferenceCatalog._makeSchema(filterNameList=["r"], addCentroid=True)
259 refCat = afwTable.SimpleCatalog(minimalPosRefSchema)
260 refCat.reserve(len(sourceCat))
261 for source in sourceCat:
262 refObj = refCat.addNew()
263 refObj.setCoord(source.getCoord())
264 refObj.set("centroid_x", source.getX())
265 refObj.set("centroid_y", source.getY())
266 refObj.set("hasCentroid", True)
267 refObj.set("r_flux", np.random.uniform(1, 10000))
268 refObj.set("r_fluxErr", source.get("slot_ApFlux_instFluxErr"))
269 refObj.setId(source.getId())
270 return refCat
272 def loadSourceCatalog(self, filename):
273 """Load a list of xy points from a file, set coord, and return a
274 SourceSet of points
276 """
277 sourceCat = afwTable.SourceCatalog.readFits(filename)
278 aliasMap = sourceCat.schema.getAliasMap()
279 aliasMap.set("slot_ApFlux", "base_PsfFlux")
280 instFluxKey = sourceCat.schema["slot_ApFlux_instFlux"].asKey()
281 instFluxErrKey = sourceCat.schema["slot_ApFlux_instFluxErr"].asKey()
283 # Source x,y positions are ~ (500,1500) x (500,1500)
284 centroidKey = sourceCat.table.getCentroidSlot().getMeasKey()
285 for src in sourceCat:
286 adjCentroid = src.get(centroidKey) - lsst.geom.Extent2D(500, 500)
287 src.set(centroidKey, adjCentroid)
288 src.set(instFluxKey, 1000)
289 src.set(instFluxErrKey, 1)
291 # Set catalog coord
292 for src in sourceCat:
293 src.updateCoord(self.wcs)
294 return sourceCat
297class MemoryTester(lsst.utils.tests.MemoryTestCase):
298 pass
301def setup_module(module):
302 lsst.utils.tests.init()
305if __name__ == "__main__": 305 ↛ 307line 305 didn't jump to line 307 because the condition on line 305 was never true
307 lsst.utils.tests.init()
308 unittest.main()