Coverage for python / lsst / meas / astrom / matchOptimisticBTask.py: 25%

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 

22__all__ = ["MatchOptimisticBTask", "MatchOptimisticBConfig", 

23 "MatchTolerance"] 

24 

25import math 

26 

27import lsst.pex.config as pexConfig 

28import lsst.pipe.base as pipeBase 

29from lsst.utils.timer import timeMethod 

30 

31from .setMatchDistance import setMatchDistance 

32from .matchOptimisticB import matchOptimisticB, MatchOptimisticBControl 

33 

34 

35class MatchTolerance: 

36 """Stores match tolerances for use in `lsst.meas.astrom.AstrometryTask` and 

37 later iterations of the matcher. 

38 

39 MatchOptimsiticBTask relies on a maximum distance for matching 

40 set by either the default in MatchOptimisticBConfig or the 2 sigma 

41 scatter found after AstrometryTask has fit for a wcs. 

42 

43 Parameters 

44 ---------- 

45 maxMatchDist : `lsst.geom.Angle` 

46 Current maximum distance to consider a match. 

47 """ 

48 

49 def __init__(self, maxMatchDist=None): 

50 self.maxMatchDist = maxMatchDist 

51 

52 

53class MatchOptimisticBConfig(pexConfig.Config): 

54 """Configuration for MatchOptimisticBTask 

55 """ 

56 maxMatchDistArcSec = pexConfig.RangeField( 

57 doc="Maximum separation between reference objects and sources " 

58 "beyond which they will not be considered a match (arcsec)", 

59 dtype=float, 

60 default=2.0, 

61 min=0, 

62 ) 

63 numBrightStars = pexConfig.RangeField( 

64 doc="Maximum number of bright stars to use in fit.", 

65 dtype=int, 

66 default=150, 

67 min=2, 

68 ) 

69 minMatchedPairs = pexConfig.RangeField( 

70 doc="Minimum number of matched pairs; see also minFracMatchedPairs", 

71 dtype=int, 

72 default=30, 

73 min=2, 

74 ) 

75 minFracMatchedPairs = pexConfig.RangeField( 

76 doc="Minimum number of matched pairs as a fraction of the smaller of " 

77 "the number of reference stars or the number of good sources; " 

78 "the actual minimum is the smaller of this value or minMatchedPairs", 

79 dtype=float, 

80 default=0.3, 

81 min=0, 

82 max=1, 

83 ) 

84 maxOffsetPix = pexConfig.RangeField( 

85 doc="Maximum allowed shift of WCS, due to matching (pixel). " 

86 "When changing this value, the LoadReferenceObjectsConfig.pixelMargin should also be updated.", 

87 dtype=int, 

88 default=250, 

89 max=4000, 

90 ) 

91 maxRotationDeg = pexConfig.RangeField( 

92 doc="Rotation angle allowed between sources and position reference objects (degrees)", 

93 dtype=float, 

94 default=1.0, 

95 max=6.0, 

96 ) 

97 allowedNonperpDeg = pexConfig.RangeField( 

98 doc="Allowed non-perpendicularity of x and y (degree)", 

99 dtype=float, 

100 default=0.2, 

101 max=45.0, 

102 ) 

103 numPointsForShape = pexConfig.Field( 

104 doc="number of points to define a shape for matching", 

105 dtype=int, 

106 default=6, 

107 ) 

108 maxDeterminant = pexConfig.Field( 

109 doc="maximum determinant of linear transformation matrix for a usable solution", 

110 dtype=float, 

111 default=0.02, 

112 ) 

113 

114 

115# The following block adds links to this task from the Task Documentation page. 

116# \addtogroup LSST_task_documentation 

117# \{ 

118# \page measAstrom_matchOptimisticBTask 

119# \ref MatchOptimisticBTask "MatchOptimisticBTask" 

120# Match sources to reference objects 

121# \} 

122 

123 

124class MatchOptimisticBTask(pipeBase.Task): 

125 """Match sources to reference objects using the Optimistic Pattern Matcher 

126 B algorithm of Tabur 2007. 

127 """ 

128 ConfigClass = MatchOptimisticBConfig 

129 _DefaultName = "matchObjectsToSources" 

130 

131 def __init__(self, **kwargs): 

132 pipeBase.Task.__init__(self, **kwargs) 

133 

134 def filterStars(self, refCat): 

135 """Extra filtering pass; subclass if desired. 

136 

137 Parameters 

138 ---------- 

139 refCat : `lsst.afw.table.SimpleCatalog` 

140 Catalog of reference objects. 

141 

142 Returns 

143 ------- 

144 trimmedRefCat : `lsst.afw.table.SimpleCatalog` 

145 Reference catalog with some filtering applied. Currently no 

146 filtering is applied. 

147 """ 

148 return refCat 

149 

150 @timeMethod 

151 def matchObjectsToSources(self, refCat, sourceCat, wcs, sourceFluxField, refFluxField, 

152 matchTolerance=None, bbox=None): 

153 """Match sources to position reference stars. 

154 

155 Parameters 

156 ---------- 

157 refCat : `lsst.afw.table.SimpleCatalog` 

158 Reference catalog to match. 

159 sourceCat : `lsst.afw.table.SourceCatalog` 

160 Catalog of sources found on an exposure. This should already be 

161 down-selected to "good"/"usable" sources in the calling Task. 

162 wcs : `lsst.afw.geom.SkyWcs` 

163 Current WCS of the exposure containing the sources. 

164 sourceFluxField : `str` 

165 Field of the sourceCat to use for flux 

166 refFluxField : `str` 

167 Field of the refCat to use for flux 

168 matchTolerance : `lsst.meas.astrom.MatchTolerance` 

169 Object containing information from previous 

170 `lsst.meas.astrom.AstrometryTask` match/fit cycles for use in 

171 matching. If `None` is config defaults. 

172 bbox : `lsst.geom.Box2I`, optional 

173 Bounding box of the exposure for evaluating the local pixelScale 

174 (defaults to the Sky Origin of the ``wcs`` provided if ``bbox`` 

175 is `None`). 

176 

177 Returns 

178 ------- 

179 matchResult : `lsst.pipe.base.Struct` 

180 Result struct with components 

181 

182 - ``matches`` : List of matches with distance below the maximum match 

183 distance (`list` of `lsst.afw.table.ReferenceMatch`). 

184 - ``useableSourceCat`` : Catalog of sources matched and suited for 

185 WCS fitting (`lsst.afw.table.SourceCatalog`). 

186 - ``matchTolerance`` : MatchTolerance object updated from this 

187 match iteration (`lsst.meas.astrom.MatchTolerance`). 

188 """ 

189 import lsstDebug 

190 debug = lsstDebug.Info(__name__) 

191 

192 preNumObj = len(refCat) 

193 refCat = self.filterStars(refCat) 

194 numRefObj = len(refCat) 

195 

196 if self.log: 

197 self.log.info("filterStars purged %d reference stars, leaving %d stars", 

198 preNumObj - numRefObj, numRefObj) 

199 

200 if matchTolerance is None: 

201 matchTolerance = MatchTolerance() 

202 

203 # Make a name alias here for consistency with older code, and to make 

204 # it clear that this is a good/usable (cleaned) source catalog. 

205 usableSourceCat = sourceCat 

206 

207 numUsableSources = len(usableSourceCat) 

208 

209 if len(usableSourceCat) == 0: 

210 raise pipeBase.TaskError("No sources are usable") 

211 

212 minMatchedPairs = min(self.config.minMatchedPairs, 

213 int(self.config.minFracMatchedPairs * min([len(refCat), len(usableSourceCat)]))) 

214 

215 # match usable (possibly saturated) sources and then purge saturated sources from the match list 

216 usableMatches = self._doMatch( 

217 refCat=refCat, 

218 sourceCat=usableSourceCat, 

219 wcs=wcs, 

220 refFluxField=refFluxField, 

221 numUsableSources=numUsableSources, 

222 minMatchedPairs=minMatchedPairs, 

223 maxMatchDist=matchTolerance.maxMatchDist, 

224 sourceFluxField=sourceFluxField, 

225 verbose=debug.verbose, 

226 bbox=bbox, 

227 ) 

228 

229 # cull non-good sources 

230 matches = [] 

231 self._getIsGoodKeys(usableSourceCat.schema) 

232 for match in usableMatches: 

233 if self._isGoodTest(match.second): 

234 # Append the isGood match. 

235 matches.append(match) 

236 

237 self.log.debug("Found %d usable matches, of which %d had good sources", 

238 len(usableMatches), len(matches)) 

239 

240 if len(matches) == 0: 

241 raise RuntimeError("Unable to match sources") 

242 

243 self.log.info("Matched %d sources", len(matches)) 

244 if len(matches) < minMatchedPairs: 

245 self.log.warning("Number of matches is smaller than request") 

246 

247 return pipeBase.Struct( 

248 matches=matches, 

249 usableSourceCat=usableSourceCat, 

250 matchTolerance=matchTolerance, 

251 ) 

252 

253 def _getIsGoodKeys(self, schema): 

254 """Retrieve the keys needed for the isGoodTest from the source catalog 

255 schema. 

256 

257 Parameters 

258 ---------- 

259 schema : `lsst.afw.table.Schema` 

260 Source schema to retrieve `lsst.afw.table.Key` s from. 

261 """ 

262 self.edgeKey = schema["base_PixelFlags_flag_edge"].asKey() 

263 self.interpolatedCenterKey = schema["base_PixelFlags_flag_interpolatedCenter"].asKey() 

264 self.saturatedKey = schema["base_PixelFlags_flag_saturated"].asKey() 

265 

266 def _isGoodTest(self, source): 

267 """Test that an object is good for use in the WCS fitter. 

268 

269 This is a hard coded version of the isGood flag from the old SourceInfo 

270 class that used to be part of this class. 

271 

272 Parameters 

273 ---------- 

274 source : `lsst.afw.table.SourceRecord` 

275 Source to test. 

276 

277 Returns 

278 ------- 

279 isGood : `bool` 

280 Source passes CCD edge and saturated tests. 

281 """ 

282 return (not source.get(self.edgeKey) 

283 and not source.get(self.interpolatedCenterKey) 

284 and not source.get(self.saturatedKey)) 

285 

286 @timeMethod 

287 def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources, minMatchedPairs, 

288 maxMatchDist, sourceFluxField, verbose, bbox=None): 

289 """Implementation of matching sources to position reference stars. 

290 

291 Unlike matchObjectsToSources, this method does not check if the sources 

292 are suitable. 

293 

294 Parameters 

295 ---------- 

296 refCat : `lsst.afw.table.SimpleCatalog` 

297 Catalog of reference objects. 

298 sourceCat : `lsst.afw.table.SourceCatalog` 

299 Catalog of detected sources. 

300 wcs : `lsst.afw.geom.SkyWcs` 

301 Current best WCS of the image. 

302 refFluxFioeld : `str` 

303 Name of flux field in refCat to use. 

304 numUsableSources : `int` 

305 Total number of source usable for matching. 

306 mintMatchPairs : `int` 

307 Minimum number of objects to match between the refCat and sourceCat 

308 to consider a valid match. 

309 maxMatchDist : `lsst.geom.Angle` 

310 Maximum separation to considering a reference and a source a match. 

311 sourceFluxField : `str` 

312 Name of source catalog flux field. 

313 verbose : `bool` 

314 Print diagnostic information std::cout 

315 bbox : `lsst.geom.Box2I`, optional 

316 Bounding box of the exposure for evaluating the local pixelScale 

317 (defaults to the Sky Origin of the ``wcs`` provided if ``bbox`` 

318 is None). 

319 

320 Returns 

321 ------- 

322 matches : `list` of `lsst.afw.table.ReferenceMatch` 

323 """ 

324 numSources = len(sourceCat) 

325 posRefBegInd = numUsableSources - numSources 

326 if maxMatchDist is None: 

327 maxMatchDistArcSec = self.config.maxMatchDistArcSec 

328 else: 

329 maxMatchDistArcSec = min(maxMatchDist.asArcseconds(), self.config.maxMatchDistArcSec) 

330 

331 if bbox is not None: 

332 pixelScale = wcs.getPixelScale(bbox.getCenter()).asArcseconds() 

333 else: 

334 pixelScale = wcs.getPixelScale().asArcseconds() 

335 

336 configMatchDistPix = maxMatchDistArcSec/pixelScale 

337 

338 matchControl = MatchOptimisticBControl() 

339 matchControl.refFluxField = refFluxField 

340 matchControl.sourceFluxField = sourceFluxField 

341 matchControl.numBrightStars = self.config.numBrightStars 

342 matchControl.minMatchedPairs = self.config.minMatchedPairs 

343 matchControl.maxOffsetPix = self.config.maxOffsetPix 

344 matchControl.numPointsForShape = self.config.numPointsForShape 

345 matchControl.maxDeterminant = self.config.maxDeterminant 

346 

347 for maxRotInd in range(4): 

348 matchControl.maxRotationDeg = self.config.maxRotationDeg * math.pow(2.0, 0.5*maxRotInd) 

349 for matchRadInd in range(3): 

350 matchControl.matchingAllowancePix = configMatchDistPix * math.pow(1.25, matchRadInd) 

351 

352 for angleDiffInd in range(3): 

353 matchControl.allowedNonperpDeg = self.config.allowedNonperpDeg*(angleDiffInd+1) 

354 matches = matchOptimisticB( 

355 refCat, 

356 sourceCat, 

357 matchControl, 

358 wcs, 

359 posRefBegInd, 

360 verbose, 

361 ) 

362 if matches is not None and len(matches) > 0: 

363 setMatchDistance(matches) 

364 return matches 

365 return matches