3 from scipy.spatial
import cKDTree
9 from lsst.utils.timer
import timeMethod
11 from .matchOptimisticBTask
import MatchTolerance
13 from .pessimistic_pattern_matcher_b_3D
import PessimisticPatternMatcherB
15 __all__ = [
"MatchPessimisticBTask",
"MatchPessimisticBConfig",
16 "MatchTolerancePessimistic"]
20 """Stores match tolerances for use in AstrometryTask and later
21 iterations of the matcher.
23 MatchPessimisticBTask relies on several state variables to be
24 preserved over different iterations in the
25 AstrometryTask.matchAndFitWcs loop of AstrometryTask.
29 maxMatchDist : `lsst.geom.Angle`
30 Maximum distance to consider a match from the previous match/fit
32 autoMaxMatchDist : `lsst.geom.Angle`
33 Automated estimation of the maxMatchDist from the sky statistics of the
34 source and reference catalogs.
35 maxShift : `lsst.geom.Angle`
36 Maximum shift found in the previous match/fit cycle.
37 lastMatchedPattern : `int`
38 Index of the last source pattern that was matched into the reference
40 failedPatternList : `list` of `int`
41 Previous matches were found to be false positives.
42 PPMbObj : `lsst.meas.astrom.PessimisticPatternMatcherB`
43 Initialized Pessimistic pattern matcher object. Storing this prevents
44 the need for recalculation of the searchable distances in the PPMB.
47 def __init__(self, maxMatchDist=None, autoMaxMatchDist=None,
48 maxShift=None, lastMatchedPattern=None,
49 failedPatternList=None, PPMbObj=None):
55 if failedPatternList
is None:
62 """Configuration for MatchPessimisticBTask
64 numBrightStars = pexConfig.RangeField(
65 doc=
"Number of bright stars to use. Sets the max number of patterns "
66 "that can be tested.",
71 minMatchedPairs = pexConfig.RangeField(
72 doc=
"Minimum number of matched pairs; see also minFracMatchedPairs.",
77 minFracMatchedPairs = pexConfig.RangeField(
78 doc=
"Minimum number of matched pairs as a fraction of the smaller of "
79 "the number of reference stars or the number of good sources; "
80 "the actual minimum is the smaller of this value or "
87 matcherIterations = pexConfig.RangeField(
88 doc=
"Number of softening iterations in matcher.",
93 maxOffsetPix = pexConfig.RangeField(
94 doc=
"Maximum allowed shift of WCS, due to matching (pixel). "
95 "When changing this value, the "
96 "LoadReferenceObjectsConfig.pixelMargin should also be updated.",
101 maxRotationDeg = pexConfig.RangeField(
102 doc=
"Rotation angle allowed between sources and position reference "
103 "objects (degrees).",
108 numPointsForShape = pexConfig.Field(
109 doc=
"Number of points to define a shape for matching.",
113 numPointsForShapeAttempt = pexConfig.Field(
114 doc=
"Number of points to try for creating a shape. This value should "
115 "be greater than or equal to numPointsForShape. Besides "
116 "loosening the signal to noise cut in the 'matcher' SourceSelector, "
117 "increasing this number will solve CCDs where no match was found.",
121 minMatchDistPixels = pexConfig.RangeField(
122 doc=
"Distance in units of pixels to always consider a source-"
123 "reference pair a match. This prevents the astrometric fitter "
124 "from over-fitting and removing stars that should be matched and "
125 "allows for inclusion of new matches as the wcs improves.",
131 numPatternConsensus = pexConfig.Field(
132 doc=
"Number of implied shift/rotations from patterns that must agree "
133 "before it a given shift/rotation is accepted. This is only used "
134 "after the first softening iteration fails and if both the "
135 "number of reference and source objects is greater than "
140 numRefRequireConsensus = pexConfig.Field(
141 doc=
"If the available reference objects exceeds this number, "
142 "consensus/pessimistic mode will enforced regardless of the "
143 "number of available sources. Below this optimistic mode ("
144 "exit at first match rather than requiring numPatternConsensus to "
145 "be matched) can be used. If more sources are required to match, "
146 "decrease the signal to noise cut in the sourceSelector.",
150 maxRefObjects = pexConfig.RangeField(
151 doc=
"Maximum number of reference objects to use for the matcher. The "
152 "absolute maximum allowed for is 2 ** 16 for memory reasons.",
160 pexConfig.Config.validate(self)
162 raise ValueError(
"numPointsForShapeAttempt must be greater than "
163 "or equal to numPointsForShape.")
165 raise ValueError(
"numBrightStars must be greater than "
166 "numPointsForShape.")
179 """Match sources to reference objects.
182 ConfigClass = MatchPessimisticBConfig
183 _DefaultName =
"matchObjectsToSources"
186 pipeBase.Task.__init__(self, **kwargs)
190 match_tolerance=None):
191 """Match sources to position reference stars
193 refCat : `lsst.afw.table.SimpleCatalog`
194 catalog of reference objects that overlap the exposure; reads
198 - the specified flux field
200 sourceCat : `lsst.afw.table.SourceCatalog`
201 Catalog of sources found on an exposure. This should already be
202 down-selected to "good"/"usable" sources in the calling Task.
203 wcs : `lsst.afw.geom.SkyWcs`
205 sourceFluxField: `str`
206 field of sourceCat to use for flux
208 field of refCat to use for flux
209 match_tolerance : `lsst.meas.astrom.MatchTolerancePessimistic`
210 is a MatchTolerance class object or `None`. This this class is used
211 to communicate state between AstrometryTask and MatcherTask.
212 AstrometryTask will also set the MatchTolerance class variable
213 maxMatchDist based on the scatter AstrometryTask has found after
218 result : `lsst.pipe.base.Struct`
219 Result struct with components:
221 - ``matches`` : source to reference matches found (`list` of
222 `lsst.afw.table.ReferenceMatch`)
223 - ``usableSourceCat`` : a catalog of sources potentially usable for
224 matching and WCS fitting (`lsst.afw.table.SourceCatalog`).
225 - ``match_tolerance`` : a MatchTolerance object containing the
226 resulting state variables from the match
227 (`lsst.meas.astrom.MatchTolerancePessimistic`).
234 if match_tolerance
is None:
239 goodSourceCat = sourceCat
241 numUsableSources = len(goodSourceCat)
243 if len(goodSourceCat) == 0:
244 raise pipeBase.TaskError(
"No sources are good")
246 minMatchedPairs =
min(self.config.minMatchedPairs,
247 int(self.config.minFracMatchedPairs
248 *
min([len(refCat), len(goodSourceCat)])))
250 if len(refCat) > self.config.maxRefObjects:
252 "WARNING: Reference catalog larger that maximum allowed. "
253 "Trimming to %i", self.config.maxRefObjects)
254 trimmedRefCat = self.
_filterRefCat_filterRefCat(refCat, refFluxField)
256 trimmedRefCat = refCat
258 doMatchReturn = self.
_doMatch_doMatch(
259 refCat=trimmedRefCat,
260 sourceCat=goodSourceCat,
262 refFluxField=refFluxField,
263 numUsableSources=numUsableSources,
264 minMatchedPairs=minMatchedPairs,
265 match_tolerance=match_tolerance,
266 sourceFluxField=sourceFluxField,
267 verbose=debug.verbose,
269 matches = doMatchReturn.matches
270 match_tolerance = doMatchReturn.match_tolerance
272 if len(matches) == 0:
273 raise RuntimeError(
"Unable to match sources")
275 self.log.
info(
"Matched %d sources", len(matches))
276 if len(matches) < minMatchedPairs:
277 self.log.
warning(
"Number of matches is smaller than request")
279 return pipeBase.Struct(
281 usableSourceCat=goodSourceCat,
282 match_tolerance=match_tolerance,
285 def _filterRefCat(self, refCat, refFluxField):
286 """Sub-select a number of reference objects starting from the brightest
287 and maxing out at the number specified by maxRefObjects in the config.
289 No trimming is done if len(refCat) > config.maxRefObjects.
293 refCat : `lsst.afw.table.SimpleCatalog`
294 Catalog of reference objects to trim.
296 field of refCat to use for flux
299 outCat : `lsst.afw.table.SimpleCatalog`
300 Catalog trimmed to the number set in the task config from the
304 if len(refCat) <= self.config.maxRefObjects:
306 fluxArray = refCat.get(refFluxField)
307 sortedFluxArray = fluxArray[fluxArray.argsort()]
308 minFlux = sortedFluxArray[-(self.config.maxRefObjects + 1)]
310 selected = (refCat.get(refFluxField) > minFlux)
313 outCat.reserve(self.config.maxRefObjects)
314 outCat.extend(refCat[selected])
319 def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources,
320 minMatchedPairs, match_tolerance, sourceFluxField, verbose):
321 """Implementation of matching sources to position reference objects
323 Unlike matchObjectsToSources, this method does not check if the sources
328 refCat : `lsst.afw.table.SimpleCatalog`
329 catalog of position reference objects that overlap an exposure
330 sourceCat : `lsst.afw.table.SourceCatalog`
331 catalog of sources found on the exposure
332 wcs : `lsst.afw.geom.SkyWcs`
333 estimated WCS of exposure
335 field of refCat to use for flux
336 numUsableSources : `int`
337 number of usable sources (sources with known centroid that are not
338 near the edge, but may be saturated)
339 minMatchedPairs : `int`
340 minimum number of matches
341 match_tolerance : `lsst.meas.astrom.MatchTolerancePessimistic`
342 a MatchTolerance object containing variables specifying matcher
343 tolerances and state from possible previous runs.
344 sourceFluxField : `str`
345 Name of the flux field in the source catalog.
347 Set true to print diagnostic information to std::cout
352 Results struct with components:
354 - ``matches`` : a list the matches found
355 (`list` of `lsst.afw.table.ReferenceMatch`).
356 - ``match_tolerance`` : MatchTolerance containing updated values from
357 this fit iteration (`lsst.meas.astrom.MatchTolerancePessimistic`)
366 src_array = np.empty((len(sourceCat), 4), dtype=np.float64)
367 for src_idx, srcObj
in enumerate(sourceCat):
368 coord = wcs.pixelToSky(srcObj.getCentroid())
369 theta = np.pi / 2 - coord.getLatitude().asRadians()
370 phi = coord.getLongitude().asRadians()
371 flux = srcObj[sourceFluxField]
372 src_array[src_idx, :] = \
375 if match_tolerance.PPMbObj
is None or \
376 match_tolerance.autoMaxMatchDist
is None:
380 ref_array = np.empty((len(refCat), 4), dtype=np.float64)
381 for ref_idx, refObj
in enumerate(refCat):
382 theta = np.pi / 2 - refObj.getDec().asRadians()
383 phi = refObj.getRa().asRadians()
384 flux = refObj[refFluxField]
385 ref_array[ref_idx, :] = \
389 ref_array[:, :3], self.log)
390 self.log.
debug(
"Computing source statistics...")
393 self.log.
debug(
"Computing reference statistics...")
396 maxMatchDistArcSec = np.max((
397 self.config.minMatchDistPixels
398 * wcs.getPixelScale().asArcseconds(),
399 np.min((maxMatchDistArcSecSrc,
400 maxMatchDistArcSecRef))))
401 match_tolerance.autoMaxMatchDist =
geom.Angle(
402 maxMatchDistArcSec, geom.arcseconds)
406 if match_tolerance.maxShift
is None:
407 maxShiftArcseconds = (self.config.maxOffsetPix
408 * wcs.getPixelScale().asArcseconds())
412 maxShiftArcseconds = np.max(
413 (match_tolerance.maxShift.asArcseconds(),
414 self.config.minMatchDistPixels
415 * wcs.getPixelScale().asArcseconds()))
421 if match_tolerance.maxMatchDist
is None:
422 match_tolerance.maxMatchDist = match_tolerance.autoMaxMatchDist
424 maxMatchDistArcSec = np.max(
425 (self.config.minMatchDistPixels
426 * wcs.getPixelScale().asArcseconds(),
427 np.min((match_tolerance.maxMatchDist.asArcseconds(),
428 match_tolerance.autoMaxMatchDist.asArcseconds()))))
434 numConsensus = self.config.numPatternConsensus
435 if len(refCat) < self.config.numRefRequireConsensus:
436 minObjectsForConsensus = \
437 self.config.numBrightStars + \
438 self.config.numPointsForShapeAttempt
439 if len(refCat) < minObjectsForConsensus
or \
440 len(sourceCat) < minObjectsForConsensus:
443 self.log.
debug(
"Current tol maxDist: %.4f arcsec",
445 self.log.
debug(
"Current shift: %.4f arcsec",
450 for soften_dist
in range(self.config.matcherIterations):
451 if soften_dist == 0
and \
452 match_tolerance.lastMatchedPattern
is not None:
461 run_n_consent = numConsensus
464 matcher_struct = match_tolerance.PPMbObj.match(
465 source_array=src_array,
466 n_check=self.config.numPointsForShapeAttempt,
467 n_match=self.config.numPointsForShape,
468 n_agree=run_n_consent,
469 max_n_patterns=self.config.numBrightStars,
470 max_shift=maxShiftArcseconds,
471 max_rotation=self.config.maxRotationDeg,
472 max_dist=maxMatchDistArcSec * 2. ** soften_dist,
473 min_matches=minMatchedPairs,
474 pattern_skip_array=np.array(
475 match_tolerance.failedPatternList)
478 if soften_dist == 0
and \
479 len(matcher_struct.match_ids) == 0
and \
480 match_tolerance.lastMatchedPattern
is not None:
487 match_tolerance.failedPatternList.append(
488 match_tolerance.lastMatchedPattern)
489 match_tolerance.lastMatchedPattern =
None
490 maxShiftArcseconds = \
491 self.config.maxOffsetPix * wcs.getPixelScale().asArcseconds()
492 elif len(matcher_struct.match_ids) > 0:
495 match_tolerance.maxShift = \
496 matcher_struct.shift * geom.arcseconds
497 match_tolerance.lastMatchedPattern = \
498 matcher_struct.pattern_idx
504 return pipeBase.Struct(
506 match_tolerance=match_tolerance,
518 distances_arcsec = np.degrees(matcher_struct.distances_rad) * 3600
519 dist_cut_arcsec = np.max(
520 (np.degrees(matcher_struct.max_dist_rad) * 3600,
521 self.config.minMatchDistPixels * wcs.getPixelScale().asArcseconds()))
526 for match_id_pair, dist_arcsec
in zip(matcher_struct.match_ids,
528 if dist_arcsec < dist_cut_arcsec:
530 match.first = refCat[int(match_id_pair[1])]
531 match.second = sourceCat[int(match_id_pair[0])]
535 match.distance = match.first.getCoord().separation(
536 match.second.getCoord()).asArcseconds()
537 matches.append(match)
539 return pipeBase.Struct(
541 match_tolerance=match_tolerance,
544 def _latlong_flux_to_xyz_mag(self, theta, phi, flux):
545 """Convert angles theta and phi and a flux into unit sphere
546 x, y, z, and a relative magnitude.
548 Takes in a afw catalog object and converts the catalog object RA, DECs
549 to points on the unit sphere. Also converts the flux into a simple,
550 non-zero-pointed magnitude for relative sorting.
555 Angle from the north pole (z axis) of the sphere
557 Rotation around the sphere
561 output_array : `numpy.ndarray`, (N, 4)
562 Spherical unit vector x, y, z with flux.
564 output_array = np.empty(4, dtype=np.float64)
565 output_array[0] = np.sin(theta)*np.cos(phi)
566 output_array[1] = np.sin(theta)*np.sin(phi)
567 output_array[2] = np.cos(theta)
569 output_array[3] = -2.5 * np.log10(flux)
573 output_array[3] = 99.
577 def _get_pair_pattern_statistics(self, cat_array):
578 """ Compute the tolerances for the matcher automatically by comparing
579 pinwheel patterns as we would in the matcher.
581 We test how similar the patterns we can create from a given set of
582 objects by computing the spoke lengths for each pattern and sorting
583 them from smallest to largest. The match tolerance is the average
584 distance per spoke between the closest two patterns in the sorted
589 cat_array : `numpy.ndarray`, (N, 3)
590 array of 3 vectors representing the x, y, z position of catalog
591 objects on the unit sphere.
596 Suggested max match tolerance distance calculated from comparisons
597 between pinwheel patterns used in optimistic/pessimistic pattern
601 self.log.
debug(
"Starting automated tolerance calculation...")
605 pattern_array = np.empty(
606 (cat_array.shape[0] - self.config.numPointsForShape,
607 self.config.numPointsForShape - 1))
608 flux_args_array = np.argsort(cat_array[:, -1])
611 tmp_sort_array = cat_array[flux_args_array]
614 for start_idx
in range(cat_array.shape[0]
615 - self.config.numPointsForShape):
616 pattern_points = tmp_sort_array[start_idx:start_idx
617 + self.config.numPointsForShape, :-1]
618 pattern_delta = pattern_points[1:, :] - pattern_points[0, :]
619 pattern_array[start_idx, :] = np.sqrt(
620 pattern_delta[:, 0] ** 2
621 + pattern_delta[:, 1] ** 2
622 + pattern_delta[:, 2] ** 2)
627 pattern_array[start_idx, :] = pattern_array[
628 start_idx, np.argsort(pattern_array[start_idx, :])]
634 pattern_array[:, :(self.config.numPointsForShape - 1)])
635 dist_nearest_array, ids = dist_tree.query(
636 pattern_array[:, :(self.config.numPointsForShape - 1)], k=2)
637 dist_nearest_array = dist_nearest_array[:, 1]
638 dist_nearest_array.sort()
642 dist_tol = (np.degrees(dist_nearest_array[dist_idx]) * 3600.
643 / (self.config.numPointsForShape - 1.))
645 self.log.
debug(
"Automated tolerance")
646 self.log.
debug(
"\tdistance/match tol: %.4f [arcsec]", dist_tol)
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
A class representing an angle.
def _get_pair_pattern_statistics(self, cat_array)
def _latlong_flux_to_xyz_mag(self, theta, phi, flux)
def matchObjectsToSources(self, refCat, sourceCat, wcs, sourceFluxField, refFluxField, match_tolerance=None)
def _filterRefCat(self, refCat, refFluxField)
def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources, minMatchedPairs, match_tolerance, sourceFluxField, verbose)
def __init__(self, **kwargs)
def __init__(self, maxMatchDist=None, autoMaxMatchDist=None, maxShift=None, lastMatchedPattern=None, failedPatternList=None, PPMbObj=None)
Lightweight representation of a geometric match between two records.