Coverage for tests/test_matchPessimisticB.py: 16%

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/>. 

21 

22import math 

23import os 

24import unittest 

25 

26import numpy as np 

27 

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 

35 

36 

37class TestMatchPessimisticB(unittest.TestCase): 

38 

39 def setUp(self): 

40 np.random.seed(12345) 

41 

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) 

49 

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 

54 

55 self.filename = os.path.join(os.path.dirname(__file__), "cat.xy.fits") 

56 self.tolArcsec = .4 

57 self.tolPixel = .1 

58 

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 

63 

64 def testLinearXDistort(self): 

65 self.singleTestInstance(self.filename, genDistortedImage.linearXDistort) 

66 

67 def testLinearYDistort(self): 

68 self.singleTestInstance(self.filename, genDistortedImage.linearYDistort) 

69 

70 def testQuadraticDistort(self): 

71 self.singleTestInstance(self.filename, genDistortedImage.quadraticDistort) 

72 

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. 

76 

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) 

83 

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 

89 

90 self.singleTestInstance(self.filename, applyDistortion) 

91 

92 def singleTestInstance(self, filename, distortFunc, doPlot=False): 

93 sourceCat = self.loadSourceCatalog(self.filename) 

94 refCat = self.computePosRefCatalog(sourceCat) 

95 

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 tempSolver = measAstrom.AstrometryTask(config=tempConfig, refObjLoader=None) 

101 sourceSelection = tempSolver.sourceSelector.run(sourceCat) 

102 

103 distortedCat = genDistortedImage.distortList(sourceSelection.sourceCat, distortFunc) 

104 

105 if doPlot: 

106 import matplotlib.pyplot as plt 

107 

108 undistorted = [self.wcs.skyToPixel(self.distortedWcs.pixelToSky(ss.getCentroid())) 

109 for ss in distortedCat] 

110 refs = [self.wcs.skyToPixel(ss.getCoord()) for ss in refCat] 

111 

112 def plot(catalog, symbol): 

113 plt.plot([ss.getX() for ss in catalog], 

114 [ss.getY() for ss in catalog], symbol) 

115 

116 plot(distortedCat, 'b+') # Distorted positions: blue + 

117 plot(undistorted, 'g+') # Undistorted positions: green + 

118 plot(refs, 'rx') # Reference catalog: red x 

119 # The green + should overlap with the red x, because that's how 

120 # MatchPessimisticB does it. 

121 

122 plt.show() 

123 

124 sourceCat = distortedCat 

125 

126 matchRes = self.MatchPessimisticB.matchObjectsToSources( 

127 refCat=refCat, 

128 sourceCat=sourceCat, 

129 wcs=self.distortedWcs, 

130 sourceFluxField='slot_ApFlux_instFlux', 

131 refFluxField="r_flux", 

132 ) 

133 matches = matchRes.matches 

134 if doPlot: 

135 measAstrom.plotAstrometry(matches=matches, refCat=refCat, 

136 sourceCat=sourceCat) 

137 self.assertEqual(len(matches), self.expectedMatches) 

138 

139 refCoordKey = afwTable.CoordKey(refCat.schema["coord"]) 

140 srcCoordKey = afwTable.CoordKey(sourceCat.schema["coord"]) 

141 refCentroidKey = afwTable.Point2DKey(refCat.getSchema()["centroid"]) 

142 maxDistErr = 0*lsst.geom.radians 

143 

144 for refObj, source, distRad in matches: 

145 sourceCoord = source.get(srcCoordKey) 

146 refCoord = refObj.get(refCoordKey) 

147 predDist = sourceCoord.separation(refCoord) 

148 distErr = abs(predDist - distRad*lsst.geom.radians) 

149 maxDistErr = max(distErr, maxDistErr) 

150 

151 if refObj.getId() != source.getId(): 

152 refCentroid = refObj.get(refCentroidKey) 

153 sourceCentroid = source.getCentroid() 

154 radius = math.hypot(*(refCentroid - sourceCentroid)) 

155 self.fail( 

156 "ID mismatch: %s at %s != %s at %s; error = %0.1f pix" % 

157 (refObj.getId(), refCentroid, source.getId(), 

158 sourceCentroid, radius)) 

159 

160 self.assertLess(maxDistErr.asArcseconds(), 1e-7) 

161 

162 def testPassingMatcherState(self): 

163 """Test that results of the matcher can be propagated to to in 

164 subsequent iterations. 

165 """ 

166 sourceCat = self.loadSourceCatalog(self.filename) 

167 refCat = self.computePosRefCatalog(sourceCat) 

168 

169 # Apply source selector to sourceCat, using the astrometry config defaults 

170 tempConfig = measAstrom.AstrometryTask.ConfigClass() 

171 # This field isn't in the old test catalog. 

172 tempConfig.sourceSelector["science"].flags.bad.remove("base_PixelFlags_flag_nodata") 

173 tempSolver = measAstrom.AstrometryTask(config=tempConfig, refObjLoader=None) 

174 sourceSelection = tempSolver.sourceSelector.run(sourceCat) 

175 

176 distortedCat = genDistortedImage.distortList(sourceSelection.sourceCat, 

177 genDistortedImage.linearXDistort) 

178 

179 sourceCat = distortedCat 

180 

181 matchRes = self.MatchPessimisticB.matchObjectsToSources( 

182 refCat=refCat, 

183 sourceCat=sourceCat, 

184 wcs=self.distortedWcs, 

185 sourceFluxField='slot_ApFlux_instFlux', 

186 refFluxField="r_flux", 

187 ) 

188 

189 maxShift = matchRes.matchTolerance.maxShift * 300 

190 # Force the matcher to use a different pattern thatn the previous 

191 # "iteration". 

192 matchTol = measAstrom.MatchTolerancePessimistic( 

193 maxMatchDist=matchRes.matchTolerance.maxMatchDist, 

194 autoMaxMatchDist=matchRes.matchTolerance.autoMaxMatchDist, 

195 maxShift=maxShift, 

196 lastMatchedPattern=0, 

197 failedPatternList=[0], 

198 PPMbObj=matchRes.matchTolerance.PPMbObj, 

199 ) 

200 

201 matchRes = self.MatchPessimisticB.matchObjectsToSources( 

202 refCat=refCat, 

203 sourceCat=sourceCat, 

204 wcs=self.distortedWcs, 

205 sourceFluxField='slot_ApFlux_instFlux', 

206 refFluxField="r_flux", 

207 matchTolerance=matchTol, 

208 ) 

209 

210 self.assertEqual(len(matchRes.matches), self.expectedMatches) 

211 self.assertLess(matchRes.matchTolerance.maxShift, maxShift) 

212 self.assertEqual(matchRes.matchTolerance.lastMatchedPattern, 1) 

213 self.assertIsNotNone(matchRes.matchTolerance.maxMatchDist) 

214 self.assertIsNotNone(matchRes.matchTolerance.autoMaxMatchDist) 

215 self.assertIsNotNone(matchRes.matchTolerance.lastMatchedPattern) 

216 self.assertIsNotNone(matchRes.matchTolerance.failedPatternList) 

217 self.assertIsNotNone(matchRes.matchTolerance.PPMbObj) 

218 

219 def testReferenceFilter(self): 

220 """Test sub-selecting reference objects by flux.""" 

221 sourceCat = self.loadSourceCatalog(self.filename) 

222 refCat = self.computePosRefCatalog(sourceCat) 

223 

224 # Apply source selector to sourceCat, using the astrometry config defaults 

225 tempConfig = measAstrom.AstrometryTask.ConfigClass() 

226 # This field isn't in the old test catalog. 

227 tempConfig.sourceSelector["science"].flags.bad.remove("base_PixelFlags_flag_nodata") 

228 tempSolver = measAstrom.AstrometryTask(config=tempConfig, refObjLoader=None) 

229 sourceSelection = tempSolver.sourceSelector.run(sourceCat) 

230 

231 distortedCat = genDistortedImage.distortList(sourceSelection.sourceCat, 

232 genDistortedImage.linearXDistort) 

233 

234 matchPessConfig = measAstrom.MatchPessimisticBTask.ConfigClass() 

235 matchPessConfig.maxRefObjects = 150 

236 matchPessConfig.minMatchDistPixels = 5.0 

237 

238 matchPess = measAstrom.MatchPessimisticBTask(config=matchPessConfig) 

239 trimmedRefCat = matchPess._filterRefCat(refCat, 'r_flux') 

240 self.assertEqual(len(trimmedRefCat), matchPessConfig.maxRefObjects) 

241 

242 matchRes = matchPess.matchObjectsToSources( 

243 refCat=refCat, 

244 sourceCat=distortedCat, 

245 wcs=self.distortedWcs, 

246 sourceFluxField='slot_ApFlux_instFlux', 

247 refFluxField="r_flux", 

248 ) 

249 

250 self.assertEqual(len(matchRes.matches), matchPessConfig.maxRefObjects - 3) 

251 

252 def computePosRefCatalog(self, sourceCat): 

253 """Generate a position reference catalog from a source catalog 

254 """ 

255 minimalPosRefSchema = convertReferenceCatalog._makeSchema(filterNameList=["r"], addCentroid=True) 

256 refCat = afwTable.SimpleCatalog(minimalPosRefSchema) 

257 refCat.reserve(len(sourceCat)) 

258 for source in sourceCat: 

259 refObj = refCat.addNew() 

260 refObj.setCoord(source.getCoord()) 

261 refObj.set("centroid_x", source.getX()) 

262 refObj.set("centroid_y", source.getY()) 

263 refObj.set("hasCentroid", True) 

264 refObj.set("r_flux", np.random.uniform(1, 10000)) 

265 refObj.set("r_fluxErr", source.get("slot_ApFlux_instFluxErr")) 

266 refObj.setId(source.getId()) 

267 return refCat 

268 

269 def loadSourceCatalog(self, filename): 

270 """Load a list of xy points from a file, set coord, and return a 

271 SourceSet of points 

272 

273 """ 

274 sourceCat = afwTable.SourceCatalog.readFits(filename) 

275 aliasMap = sourceCat.schema.getAliasMap() 

276 aliasMap.set("slot_ApFlux", "base_PsfFlux") 

277 instFluxKey = sourceCat.schema["slot_ApFlux_instFlux"].asKey() 

278 instFluxErrKey = sourceCat.schema["slot_ApFlux_instFluxErr"].asKey() 

279 

280 # Source x,y positions are ~ (500,1500) x (500,1500) 

281 centroidKey = sourceCat.table.getCentroidSlot().getMeasKey() 

282 for src in sourceCat: 

283 adjCentroid = src.get(centroidKey) - lsst.geom.Extent2D(500, 500) 

284 src.set(centroidKey, adjCentroid) 

285 src.set(instFluxKey, 1000) 

286 src.set(instFluxErrKey, 1) 

287 

288 # Set catalog coord 

289 for src in sourceCat: 

290 src.updateCoord(self.wcs) 

291 return sourceCat 

292 

293 

294class MemoryTester(lsst.utils.tests.MemoryTestCase): 

295 pass 

296 

297 

298def setup_module(module): 

299 lsst.utils.tests.init() 

300 

301 

302if __name__ == "__main__": 302 ↛ 304line 302 didn't jump to line 304 because the condition on line 302 was never true

303 

304 lsst.utils.tests.init() 

305 unittest.main()