22__all__ = [
"MatchPessimisticBTask",
"MatchPessimisticBConfig",
23 "MatchTolerancePessimistic"]
26from scipy.spatial
import cKDTree
32from lsst.utils.timer
import timeMethod
34from .matchOptimisticBTask
import MatchTolerance
36from .pessimistic_pattern_matcher_b_3D
import PessimisticPatternMatcherB
40 """Stores match tolerances for use in AstrometryTask and later
41 iterations of the matcher.
43 MatchPessimisticBTask relies on several state variables to be
44 preserved over different iterations in the
45 AstrometryTask.matchAndFitWcs loop of AstrometryTask.
50 Maximum distance to consider a match
from the previous match/fit
53 Automated estimation of the maxMatchDist
from the sky statistics of the
54 source
and reference catalogs.
56 Maximum shift found
in the previous match/fit cycle.
57 lastMatchedPattern : `int`
58 Index of the last source pattern that was matched into the reference
60 failedPatternList : `list` of `int`
61 Previous matches were found to be false positives.
63 Initialized Pessimistic pattern matcher object. Storing this prevents
64 the need
for recalculation of the searchable distances
in the PPMB.
67 def __init__(self, maxMatchDist=None, autoMaxMatchDist=None,
68 maxShift=None, lastMatchedPattern=None,
69 failedPatternList=None, PPMbObj=None):
75 if failedPatternList
is None:
82 """Configuration for MatchPessimisticBTask
84 numBrightStars = pexConfig.RangeField(
85 doc="Number of bright stars to use. Sets the max number of patterns "
86 "that can be tested.",
91 minMatchedPairs = pexConfig.RangeField(
92 doc=
"Minimum number of matched pairs; see also minFracMatchedPairs.",
97 minFracMatchedPairs = pexConfig.RangeField(
98 doc=
"Minimum number of matched pairs as a fraction of the smaller of "
99 "the number of reference stars or the number of good sources; "
100 "the actual minimum is the smaller of this value or "
107 matcherIterations = pexConfig.RangeField(
108 doc=
"Number of softening iterations in matcher.",
113 maxOffsetPix = pexConfig.RangeField(
114 doc=
"Maximum allowed shift of WCS, due to matching (pixel). "
115 "When changing this value, the "
116 "LoadReferenceObjectsConfig.pixelMargin should also be updated.",
121 maxRotationDeg = pexConfig.RangeField(
122 doc=
"Rotation angle allowed between sources and position reference "
123 "objects (degrees).",
128 numPointsForShape = pexConfig.Field(
129 doc=
"Number of points to define a shape for matching.",
133 numPointsForShapeAttempt = pexConfig.Field(
134 doc=
"Number of points to try for creating a shape. This value should "
135 "be greater than or equal to numPointsForShape. Besides "
136 "loosening the signal to noise cut in the 'matcher' SourceSelector, "
137 "increasing this number will solve CCDs where no match was found.",
141 minMatchDistPixels = pexConfig.RangeField(
142 doc=
"Distance in units of pixels to always consider a source-"
143 "reference pair a match. This prevents the astrometric fitter "
144 "from over-fitting and removing stars that should be matched and "
145 "allows for inclusion of new matches as the wcs improves.",
151 numPatternConsensus = pexConfig.Field(
152 doc=
"Number of implied shift/rotations from patterns that must agree "
153 "before it a given shift/rotation is accepted. This is only used "
154 "after the first softening iteration fails and if both the "
155 "number of reference and source objects is greater than "
160 numRefRequireConsensus = pexConfig.Field(
161 doc=
"If the available reference objects exceeds this number, "
162 "consensus/pessimistic mode will enforced regardless of the "
163 "number of available sources. Below this optimistic mode ("
164 "exit at first match rather than requiring numPatternConsensus to "
165 "be matched) can be used. If more sources are required to match, "
166 "decrease the signal to noise cut in the sourceSelector.",
170 maxRefObjects = pexConfig.RangeField(
171 doc=
"Maximum number of reference objects to use for the matcher. The "
172 "absolute maximum allowed for is 2 ** 16 for memory reasons.",
180 pexConfig.Config.validate(self)
182 raise ValueError(
"numPointsForShapeAttempt must be greater than "
183 "or equal to numPointsForShape.")
185 raise ValueError(
"numBrightStars must be greater than "
186 "numPointsForShape.")
199 """Match sources to reference objects.
202 ConfigClass = MatchPessimisticBConfig
203 _DefaultName = "matchObjectsToSources"
206 pipeBase.Task.__init__(self, **kwargs)
210 match_tolerance=None):
211 """Match sources to position reference stars
214 catalog of reference objects that overlap the exposure; reads
218 - the specified flux field
221 Catalog of sources found on an exposure. This should already be
222 down-selected to
"good"/
"usable" sources
in the calling Task.
225 sourceFluxField: `str`
226 field of sourceCat to use
for flux
228 field of refCat to use
for flux
230 is a MatchTolerance
class object or `
None`. This this
class is used
231 to communicate state between AstrometryTask
and MatcherTask.
232 AstrometryTask will also set the MatchTolerance
class variable
233 maxMatchDist based on the scatter AstrometryTask has found after
238 result : `lsst.pipe.base.Struct`
239 Result struct
with components:
241 - ``matches`` : source to reference matches found (`list` of
243 - ``usableSourceCat`` : a catalog of sources potentially usable
for
245 - ``match_tolerance`` : a MatchTolerance object containing the
246 resulting state variables
from the match
254 if match_tolerance
is None:
259 goodSourceCat = sourceCat
261 numUsableSources = len(goodSourceCat)
263 if len(goodSourceCat) == 0:
264 raise pipeBase.TaskError(
"No sources are good")
266 minMatchedPairs =
min(self.config.minMatchedPairs,
267 int(self.config.minFracMatchedPairs
268 *
min([len(refCat), len(goodSourceCat)])))
270 if len(refCat) > self.config.maxRefObjects:
272 "WARNING: Reference catalog larger that maximum allowed. "
273 "Trimming to %i", self.config.maxRefObjects)
276 trimmedRefCat = refCat
279 refCat=trimmedRefCat,
280 sourceCat=goodSourceCat,
282 refFluxField=refFluxField,
283 numUsableSources=numUsableSources,
284 minMatchedPairs=minMatchedPairs,
285 match_tolerance=match_tolerance,
286 sourceFluxField=sourceFluxField,
287 verbose=debug.verbose,
289 matches = doMatchReturn.matches
290 match_tolerance = doMatchReturn.match_tolerance
292 if len(matches) == 0:
293 raise RuntimeError(
"Unable to match sources")
295 self.log.info(
"Matched %d sources", len(matches))
296 if len(matches) < minMatchedPairs:
297 self.log.warning(
"Number of matches is smaller than request")
299 return pipeBase.Struct(
301 usableSourceCat=goodSourceCat,
302 match_tolerance=match_tolerance,
305 def _filterRefCat(self, refCat, refFluxField):
306 """Sub-select a number of reference objects starting from the brightest
307 and maxing out at the number specified by maxRefObjects
in the config.
309 No trimming
is done
if len(refCat) > config.maxRefObjects.
314 Catalog of reference objects to trim.
316 field of refCat to use
for flux
320 Catalog trimmed to the number set
in the task config
from the
324 if len(refCat) <= self.config.maxRefObjects:
326 fluxArray = refCat.get(refFluxField)
327 sortedFluxArray = fluxArray[fluxArray.argsort()]
328 minFlux = sortedFluxArray[-(self.config.maxRefObjects + 1)]
330 selected = (refCat.get(refFluxField) > minFlux)
333 outCat.reserve(self.config.maxRefObjects)
334 outCat.extend(refCat[selected])
339 def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources,
340 minMatchedPairs, match_tolerance, sourceFluxField, verbose):
341 """Implementation of matching sources to position reference objects
343 Unlike matchObjectsToSources, this method does not check
if the sources
349 catalog of position reference objects that overlap an exposure
351 catalog of sources found on the exposure
353 estimated WCS of exposure
355 field of refCat to use
for flux
356 numUsableSources : `int`
357 number of usable sources (sources
with known centroid that are
not
358 near the edge, but may be saturated)
359 minMatchedPairs : `int`
360 minimum number of matches
362 a MatchTolerance object containing variables specifying matcher
363 tolerances
and state
from possible previous runs.
364 sourceFluxField : `str`
365 Name of the flux field
in the source catalog.
367 Set true to
print diagnostic information to std::cout
372 Results struct
with components:
374 - ``matches`` : a list the matches found
376 - ``match_tolerance`` : MatchTolerance containing updated values
from
386 src_array = np.empty((len(sourceCat), 4), dtype=np.float64)
387 for src_idx, srcObj
in enumerate(sourceCat):
388 coord = wcs.pixelToSky(srcObj.getCentroid())
389 theta = np.pi / 2 - coord.getLatitude().asRadians()
390 phi = coord.getLongitude().asRadians()
391 flux = srcObj[sourceFluxField]
392 src_array[src_idx, :] = \
395 if match_tolerance.PPMbObj
is None or \
396 match_tolerance.autoMaxMatchDist
is None:
400 ref_array = np.empty((len(refCat), 4), dtype=np.float64)
401 for ref_idx, refObj
in enumerate(refCat):
402 theta = np.pi / 2 - refObj.getDec().asRadians()
403 phi = refObj.getRa().asRadians()
404 flux = refObj[refFluxField]
405 ref_array[ref_idx, :] = \
409 ref_array[:, :3], self.log)
410 self.log.debug(
"Computing source statistics...")
413 self.log.debug(
"Computing reference statistics...")
416 maxMatchDistArcSec = np.max((
417 self.config.minMatchDistPixels
418 * wcs.getPixelScale().asArcseconds(),
419 np.min((maxMatchDistArcSecSrc,
420 maxMatchDistArcSecRef))))
421 match_tolerance.autoMaxMatchDist =
geom.Angle(
422 maxMatchDistArcSec, geom.arcseconds)
426 if match_tolerance.maxShift
is None:
427 maxShiftArcseconds = (self.config.maxOffsetPix
428 * wcs.getPixelScale().asArcseconds())
432 maxShiftArcseconds = np.max(
433 (match_tolerance.maxShift.asArcseconds(),
434 self.config.minMatchDistPixels
435 * wcs.getPixelScale().asArcseconds()))
441 if match_tolerance.maxMatchDist
is None:
442 match_tolerance.maxMatchDist = match_tolerance.autoMaxMatchDist
444 maxMatchDistArcSec = np.max(
445 (self.config.minMatchDistPixels
446 * wcs.getPixelScale().asArcseconds(),
447 np.min((match_tolerance.maxMatchDist.asArcseconds(),
448 match_tolerance.autoMaxMatchDist.asArcseconds()))))
454 numConsensus = self.config.numPatternConsensus
455 if len(refCat) < self.config.numRefRequireConsensus:
456 minObjectsForConsensus = \
457 self.config.numBrightStars + \
458 self.config.numPointsForShapeAttempt
459 if len(refCat) < minObjectsForConsensus
or \
460 len(sourceCat) < minObjectsForConsensus:
463 self.log.debug(
"Current tol maxDist: %.4f arcsec",
465 self.log.debug(
"Current shift: %.4f arcsec",
470 for soften_dist
in range(self.config.matcherIterations):
471 if soften_dist == 0
and \
472 match_tolerance.lastMatchedPattern
is not None:
481 run_n_consent = numConsensus
484 matcher_struct = match_tolerance.PPMbObj.match(
485 source_array=src_array,
486 n_check=self.config.numPointsForShapeAttempt,
487 n_match=self.config.numPointsForShape,
488 n_agree=run_n_consent,
489 max_n_patterns=self.config.numBrightStars,
490 max_shift=maxShiftArcseconds,
491 max_rotation=self.config.maxRotationDeg,
492 max_dist=maxMatchDistArcSec * 2. ** soften_dist,
493 min_matches=minMatchedPairs,
494 pattern_skip_array=np.array(
495 match_tolerance.failedPatternList)
498 if soften_dist == 0
and \
499 len(matcher_struct.match_ids) == 0
and \
500 match_tolerance.lastMatchedPattern
is not None:
507 match_tolerance.failedPatternList.append(
508 match_tolerance.lastMatchedPattern)
509 match_tolerance.lastMatchedPattern =
None
510 maxShiftArcseconds = \
511 self.config.maxOffsetPix * wcs.getPixelScale().asArcseconds()
512 elif len(matcher_struct.match_ids) > 0:
515 match_tolerance.maxShift = \
516 matcher_struct.shift * geom.arcseconds
517 match_tolerance.lastMatchedPattern = \
518 matcher_struct.pattern_idx
524 return pipeBase.Struct(
526 match_tolerance=match_tolerance,
538 distances_arcsec = np.degrees(matcher_struct.distances_rad) * 3600
539 dist_cut_arcsec = np.max(
540 (np.degrees(matcher_struct.max_dist_rad) * 3600,
541 self.config.minMatchDistPixels * wcs.getPixelScale().asArcseconds()))
546 for match_id_pair, dist_arcsec
in zip(matcher_struct.match_ids,
548 if dist_arcsec < dist_cut_arcsec:
550 match.first = refCat[int(match_id_pair[1])]
551 match.second = sourceCat[int(match_id_pair[0])]
555 match.distance = match.first.getCoord().separation(
556 match.second.getCoord()).asArcseconds()
557 matches.append(match)
559 return pipeBase.Struct(
561 match_tolerance=match_tolerance,
564 def _latlong_flux_to_xyz_mag(self, theta, phi, flux):
565 """Convert angles theta and phi and a flux into unit sphere
566 x, y, z, and a relative magnitude.
568 Takes
in a afw catalog object
and converts the catalog object RA, DECs
569 to points on the unit sphere. Also converts the flux into a simple,
570 non-zero-pointed magnitude
for relative sorting.
575 Angle
from the north pole (z axis) of the sphere
577 Rotation around the sphere
581 output_array : `numpy.ndarray`, (N, 4)
582 Spherical unit vector x, y, z
with flux.
584 output_array = np.empty(4, dtype=np.float64)
585 output_array[0] = np.sin(theta)*np.cos(phi)
586 output_array[1] = np.sin(theta)*np.sin(phi)
587 output_array[2] = np.cos(theta)
589 output_array[3] = -2.5 * np.log10(flux)
593 output_array[3] = 99.
597 def _get_pair_pattern_statistics(self, cat_array):
598 """ Compute the tolerances for the matcher automatically by comparing
599 pinwheel patterns as we would
in the matcher.
601 We test how similar the patterns we can create
from a given set of
602 objects by computing the spoke lengths
for each pattern
and sorting
603 them
from smallest to largest. The match tolerance
is the average
604 distance per spoke between the closest two patterns
in the sorted
609 cat_array : `numpy.ndarray`, (N, 3)
610 array of 3 vectors representing the x, y, z position of catalog
611 objects on the unit sphere.
616 Suggested max match tolerance distance calculated
from comparisons
617 between pinwheel patterns used
in optimistic/pessimistic pattern
621 self.log.debug("Starting automated tolerance calculation...")
625 pattern_array = np.empty(
626 (cat_array.shape[0] - self.config.numPointsForShape,
627 self.config.numPointsForShape - 1))
628 flux_args_array = np.argsort(cat_array[:, -1])
631 tmp_sort_array = cat_array[flux_args_array]
634 for start_idx
in range(cat_array.shape[0]
635 - self.config.numPointsForShape):
636 pattern_points = tmp_sort_array[start_idx:start_idx
637 + self.config.numPointsForShape, :-1]
638 pattern_delta = pattern_points[1:, :] - pattern_points[0, :]
639 pattern_array[start_idx, :] = np.sqrt(
640 pattern_delta[:, 0] ** 2
641 + pattern_delta[:, 1] ** 2
642 + pattern_delta[:, 2] ** 2)
647 pattern_array[start_idx, :] = pattern_array[
648 start_idx, np.argsort(pattern_array[start_idx, :])]
654 pattern_array[:, :(self.config.numPointsForShape - 1)])
655 dist_nearest_array, ids = dist_tree.query(
656 pattern_array[:, :(self.config.numPointsForShape - 1)], k=2)
657 dist_nearest_array = dist_nearest_array[:, 1]
658 dist_nearest_array.sort()
662 dist_tol = (np.degrees(dist_nearest_array[dist_idx]) * 3600.
663 / (self.config.numPointsForShape - 1.))
665 self.log.debug(
"Automated tolerance")
666 self.log.debug(
"\tdistance/match tol: %.4f [arcsec]", dist_tol)
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
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.