LSST Applications g070148d5b3+33e5256705,g0d53e28543+25c8b88941,g0da5cf3356+2dd1178308,g1081da9e2a+62d12e78cb,g17e5ecfddb+7e422d6136,g1c76d35bf8+ede3a706f7,g295839609d+225697d880,g2e2c1a68ba+cc1f6f037e,g2ffcdf413f+853cd4dcde,g38293774b4+62d12e78cb,g3b44f30a73+d953f1ac34,g48ccf36440+885b902d19,g4b2f1765b6+7dedbde6d2,g5320a0a9f6+0c5d6105b6,g56b687f8c9+ede3a706f7,g5c4744a4d9+ef6ac23297,g5ffd174ac0+0c5d6105b6,g6075d09f38+66af417445,g667d525e37+2ced63db88,g670421136f+2ced63db88,g71f27ac40c+2ced63db88,g774830318a+463cbe8d1f,g7876bc68e5+1d137996f1,g7985c39107+62d12e78cb,g7fdac2220c+0fd8241c05,g96f01af41f+368e6903a7,g9ca82378b8+2ced63db88,g9d27549199+ef6ac23297,gabe93b2c52+e3573e3735,gb065e2a02a+3dfbe639da,gbc3249ced9+0c5d6105b6,gbec6a3398f+0c5d6105b6,gc9534b9d65+35b9f25267,gd01420fc67+0c5d6105b6,geee7ff78d7+a14128c129,gf63283c776+ede3a706f7,gfed783d017+0c5d6105b6,w.2022.47
LSST Data Management Base Package
Loading...
Searching...
No Matches
matchPessimisticB.py
Go to the documentation of this file.
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__ = ["MatchPessimisticBTask", "MatchPessimisticBConfig",
23 "MatchTolerancePessimistic"]
24
25import numpy as np
26from scipy.spatial import cKDTree
27
28import lsst.pex.config as pexConfig
29import lsst.pipe.base as pipeBase
30import lsst.geom as geom
31import lsst.afw.table as afwTable
32from lsst.utils.timer import timeMethod
33
34from .matchOptimisticBTask import MatchTolerance
35
36from .pessimistic_pattern_matcher_b_3D import PessimisticPatternMatcherB
37
38
40 """Stores match tolerances for use in AstrometryTask and later
41 iterations of the matcher.
42
43 MatchPessimisticBTask relies on several state variables to be
44 preserved over different iterations in the
45 AstrometryTask.matchAndFitWcs loop of AstrometryTask.
46
47 Parameters
48 ----------
49 maxMatchDist : `lsst.geom.Angle`
50 Maximum distance to consider a match from the previous match/fit
51 iteration.
52 autoMaxMatchDist : `lsst.geom.Angle`
53 Automated estimation of the maxMatchDist from the sky statistics of the
54 source and reference catalogs.
55 maxShift : `lsst.geom.Angle`
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
59 data.
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.
65 """
66
67 def __init__(self, maxMatchDist=None, autoMaxMatchDist=None,
68 maxShift=None, lastMatchedPattern=None,
69 failedPatternList=None, PPMbObj=None):
70 self.maxMatchDistmaxMatchDist = maxMatchDist
71 self.autoMaxMatchDist = autoMaxMatchDist
72 self.maxShift = maxShift
73 self.lastMatchedPattern = lastMatchedPattern
74 self.PPMbObj = PPMbObj
75 if failedPatternList is None:
77 else:
78 self.failedPatternList = failedPatternList
79
80
81class MatchPessimisticBConfig(pexConfig.Config):
82 """Configuration for MatchPessimisticBTask
83 """
84 numBrightStars = pexConfig.RangeField(
85 doc="Number of bright stars to use. Sets the max number of patterns "
86 "that can be tested.",
87 dtype=int,
88 default=200,
89 min=2,
90 )
91 minMatchedPairs = pexConfig.RangeField(
92 doc="Minimum number of matched pairs; see also minFracMatchedPairs.",
93 dtype=int,
94 default=30,
95 min=2,
96 )
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 "
101 "minMatchedPairs.",
102 dtype=float,
103 default=0.3,
104 min=0,
105 max=1,
106 )
107 matcherIterations = pexConfig.RangeField(
108 doc="Number of softening iterations in matcher.",
109 dtype=int,
110 default=5,
111 min=1,
112 )
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.",
117 dtype=int,
118 default=250,
119 max=4000,
120 )
121 maxRotationDeg = pexConfig.RangeField(
122 doc="Rotation angle allowed between sources and position reference "
123 "objects (degrees).",
124 dtype=float,
125 default=1.0,
126 max=6.0,
127 )
128 numPointsForShape = pexConfig.Field(
129 doc="Number of points to define a shape for matching.",
130 dtype=int,
131 default=6,
132 )
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.",
138 dtype=int,
139 default=6,
140 )
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.",
146 dtype=float,
147 default=1.0,
148 min=0.0,
149 max=6.0,
150 )
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 "
156 "numBrightStars.",
157 dtype=int,
158 default=3,
159 )
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.",
167 dtype=int,
168 default=1000,
169 )
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.",
173 dtype=int,
174 default=2**16,
175 min=0,
176 max=2**16 + 1,
177 )
178
179 def validate(self):
180 pexConfig.Config.validate(self)
182 raise ValueError("numPointsForShapeAttempt must be greater than "
183 "or equal to numPointsForShape.")
184 if self.numPointsForShape > self.numBrightStars:
185 raise ValueError("numBrightStars must be greater than "
186 "numPointsForShape.")
187
188
189# The following block adds links to this task from the Task Documentation page.
190# \addtogroup LSST_task_documentation
191# \{
192# \page measAstrom_MatchPessimisticBTask
193# \ref MatchPessimisticBTask "MatchPessimisticBTask"
194# Match sources to reference objects
195# \}
196
197
198class MatchPessimisticBTask(pipeBase.Task):
199 """Match sources to reference objects.
200 """
201
202 ConfigClass = MatchPessimisticBConfig
203 _DefaultName = "matchObjectsToSources"
204
205 def __init__(self, **kwargs):
206 pipeBase.Task.__init__(self, **kwargs)
207
208 @timeMethod
209 def matchObjectsToSources(self, refCat, sourceCat, wcs, sourceFluxField, refFluxField,
210 match_tolerance=None):
211 """Match sources to position reference stars
212
214 catalog of reference objects that overlap the exposure; reads
215 fields for:
216
217 - coord
218 - the specified flux field
219
220 sourceCat : `lsst.afw.table.SourceCatalog`
221 Catalog of sources found on an exposure. This should already be
222 down-selected to "good"/"usable" sources in the calling Task.
224 estimated WCS
225 sourceFluxField: `str`
226 field of sourceCat to use for flux
227 refFluxField : `str`
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
234 fitting for the wcs.
235
236 Returns
237 -------
238 result : `lsst.pipe.base.Struct`
239 Result struct with components:
240
241 - ``matches`` : source to reference matches found (`list` of
243 - ``usableSourceCat`` : a catalog of sources potentially usable for
244 matching and WCS fitting (`lsst.afw.table.SourceCatalog`).
245 - ``match_tolerance`` : a MatchTolerance object containing the
246 resulting state variables from the match
248 """
249 import lsstDebug
250 debug = lsstDebug.Info(__name__)
251
252 # If we get an empty tolerance struct create the variables we need for
253 # this matcher.
254 if match_tolerance is None:
255 match_tolerance = MatchTolerancePessimistic()
256
257 # Make a name alias here for consistency with older code, and to make
258 # it clear that this is a good/usable (cleaned) source catalog.
259 goodSourceCat = sourceCat
260
261 numUsableSources = len(goodSourceCat)
262
263 if len(goodSourceCat) == 0:
264 raise pipeBase.TaskError("No sources are good")
265
266 minMatchedPairs = min(self.config.minMatchedPairs,
267 int(self.config.minFracMatchedPairs
268 * min([len(refCat), len(goodSourceCat)])))
269
270 if len(refCat) > self.config.maxRefObjects:
271 self.log.warning(
272 "WARNING: Reference catalog larger that maximum allowed. "
273 "Trimming to %i", self.config.maxRefObjects)
274 trimmedRefCat = self._filterRefCat(refCat, refFluxField)
275 else:
276 trimmedRefCat = refCat
277
278 doMatchReturn = self._doMatch(
279 refCat=trimmedRefCat,
280 sourceCat=goodSourceCat,
281 wcs=wcs,
282 refFluxField=refFluxField,
283 numUsableSources=numUsableSources,
284 minMatchedPairs=minMatchedPairs,
285 match_tolerance=match_tolerance,
286 sourceFluxField=sourceFluxField,
287 verbose=debug.verbose,
288 )
289 matches = doMatchReturn.matches
290 match_tolerance = doMatchReturn.match_tolerance
291
292 if len(matches) == 0:
293 raise RuntimeError("Unable to match sources")
294
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")
298
299 return pipeBase.Struct(
300 matches=matches,
301 usableSourceCat=goodSourceCat,
302 match_tolerance=match_tolerance,
303 )
304
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.
308
309 No trimming is done if len(refCat) > config.maxRefObjects.
310
311 Parameters
312 ----------
314 Catalog of reference objects to trim.
315 refFluxField : `str`
316 field of refCat to use for flux
317 Returns
318 -------
320 Catalog trimmed to the number set in the task config from the
321 brightest flux down.
322 """
323 # Find the flux cut that gives us the desired number of objects.
324 if len(refCat) <= self.config.maxRefObjects:
325 return refCat
326 fluxArray = refCat.get(refFluxField)
327 sortedFluxArray = fluxArray[fluxArray.argsort()]
328 minFlux = sortedFluxArray[-(self.config.maxRefObjects + 1)]
329
330 selected = (refCat.get(refFluxField) > minFlux)
331
332 outCat = afwTable.SimpleCatalog(refCat.schema)
333 outCat.reserve(self.config.maxRefObjects)
334 outCat.extend(refCat[selected])
335
336 return outCat
337
338 @timeMethod
339 def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources,
340 minMatchedPairs, match_tolerance, sourceFluxField, verbose):
341 """Implementation of matching sources to position reference objects
342
343 Unlike matchObjectsToSources, this method does not check if the sources
344 are suitable.
345
346 Parameters
347 ----------
349 catalog of position reference objects that overlap an exposure
350 sourceCat : `lsst.afw.table.SourceCatalog`
351 catalog of sources found on the exposure
353 estimated WCS of exposure
354 refFluxField : `str`
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.
366 verbose : `bool`
367 Set true to print diagnostic information to std::cout
368
369 Returns
370 -------
371 result :
372 Results struct with components:
373
374 - ``matches`` : a list the matches found
375 (`list` of `lsst.afw.table.ReferenceMatch`).
376 - ``match_tolerance`` : MatchTolerance containing updated values from
377 this fit iteration (`lsst.meas.astrom.MatchTolerancePessimistic`)
378 """
379
380 # Load the source and reference catalog as spherical points
381 # in numpy array. We do this rather than relying on internal
382 # lsst C objects for simplicity and because we require
383 # objects contiguous in memory. We need to do these slightly
384 # differently for the reference and source cats as they are
385 # different catalog objects with different fields.
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, :] = \
393 self._latlong_flux_to_xyz_mag(theta, phi, flux)
394
395 if match_tolerance.PPMbObj is None or \
396 match_tolerance.autoMaxMatchDist is None:
397 # The reference catalog is fixed per AstrometryTask so we only
398 # create the data needed if this is the first step in the match
399 # fit cycle.
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, :] = \
406 self._latlong_flux_to_xyz_mag(theta, phi, flux)
407 # Create our matcher object.
408 match_tolerance.PPMbObj = PessimisticPatternMatcherB(
409 ref_array[:, :3], self.log)
410 self.log.debug("Computing source statistics...")
411 maxMatchDistArcSecSrc = self._get_pair_pattern_statistics(
412 src_array)
413 self.log.debug("Computing reference statistics...")
414 maxMatchDistArcSecRef = self._get_pair_pattern_statistics(
415 ref_array)
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)
423
424 # Set configurable defaults when we encounter None type or set
425 # state based on previous run of AstrometryTask._matchAndFitWcs.
426 if match_tolerance.maxShift is None:
427 maxShiftArcseconds = (self.config.maxOffsetPix
428 * wcs.getPixelScale().asArcseconds())
429 else:
430 # We don't want to clamp down too hard on the allowed shift so
431 # we test that the smallest we ever allow is the pixel scale.
432 maxShiftArcseconds = np.max(
433 (match_tolerance.maxShift.asArcseconds(),
434 self.config.minMatchDistPixels
435 * wcs.getPixelScale().asArcseconds()))
436
437 # If our tolerances are not set from a previous run, estimate a
438 # starting tolerance guess from the statistics of patterns we can
439 # create on both the source and reference catalog. We use the smaller
440 # of the two.
441 if match_tolerance.maxMatchDist is None:
442 match_tolerance.maxMatchDist = match_tolerance.autoMaxMatchDist
443 else:
444 maxMatchDistArcSec = np.max(
445 (self.config.minMatchDistPixels
446 * wcs.getPixelScale().asArcseconds(),
447 np.min((match_tolerance.maxMatchDist.asArcseconds(),
448 match_tolerance.autoMaxMatchDist.asArcseconds()))))
449
450 # Make sure the data we are considering is dense enough to require
451 # the consensus mode of the matcher. If not default to Optimistic
452 # pattern matcher behavior. We enforce pessimistic mode if the
453 # reference cat is sufficiently large, avoiding false positives.
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:
461 numConsensus = 1
462
463 self.log.debug("Current tol maxDist: %.4f arcsec",
464 maxMatchDistArcSec)
465 self.log.debug("Current shift: %.4f arcsec",
466 maxShiftArcseconds)
467
468 match_found = False
469 # Start the iteration over our tolerances.
470 for soften_dist in range(self.config.matcherIterations):
471 if soften_dist == 0 and \
472 match_tolerance.lastMatchedPattern is not None:
473 # If we are on the first, most stringent tolerance,
474 # and have already found a match, the matcher should behave
475 # like an optimistic pattern matcher. Exiting at the first
476 # match.
477 run_n_consent = 1
478 else:
479 # If we fail or this is the first match attempt, set the
480 # pattern consensus to the specified config value.
481 run_n_consent = numConsensus
482 # We double the match dist tolerance each round and add 1 to the
483 # to the number of candidate spokes to check.
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)
496 )
497
498 if soften_dist == 0 and \
499 len(matcher_struct.match_ids) == 0 and \
500 match_tolerance.lastMatchedPattern is not None:
501 # If we found a pattern on a previous match-fit iteration and
502 # can't find an optimistic match on our first try with the
503 # tolerances as found in the previous match-fit,
504 # the match we found in the last iteration was likely bad. We
505 # append the bad match's index to the a list of
506 # patterns/matches to skip on subsequent iterations.
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:
513 # Match found, save a bit a state regarding this pattern
514 # in the match tolerance class object and exit.
515 match_tolerance.maxShift = \
516 matcher_struct.shift * geom.arcseconds
517 match_tolerance.lastMatchedPattern = \
518 matcher_struct.pattern_idx
519 match_found = True
520 break
521
522 # If we didn't find a match, exit early.
523 if not match_found:
524 return pipeBase.Struct(
525 matches=[],
526 match_tolerance=match_tolerance,
527 )
528
529 # The matcher returns all the nearest neighbors that agree between
530 # the reference and source catalog. For the current astrometric solver
531 # we need to remove as many false positives as possible before sending
532 # the matches off to the solver. The low value of 100 and high value of
533 # 2 are the low number of sigma and high respectively. The exact values
534 # were found after testing on data of various reference/source
535 # densities and astrometric distortion quality, specifically the
536 # visits: HSC (3358), DECam (406285, 410827),
537 # CFHT (793169, 896070, 980526).
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()))
542
543 # A match has been found, return our list of matches and
544 # return.
545 matches = []
546 for match_id_pair, dist_arcsec in zip(matcher_struct.match_ids,
547 distances_arcsec):
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])]
552 # We compute the true distance along and sphere. This isn't
553 # used in the WCS fitter however it is used in the unittest
554 # to confirm the matches computed.
555 match.distance = match.first.getCoord().separation(
556 match.second.getCoord()).asArcseconds()
557 matches.append(match)
558
559 return pipeBase.Struct(
560 matches=matches,
561 match_tolerance=match_tolerance,
562 )
563
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.
567
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.
571
572 Parameters
573 ----------
574 theta : `float`
575 Angle from the north pole (z axis) of the sphere
576 phi : `float`
577 Rotation around the sphere
578
579 Return
580 ------
581 output_array : `numpy.ndarray`, (N, 4)
582 Spherical unit vector x, y, z with flux.
583 """
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)
588 if flux > 0:
589 output_array[3] = -2.5 * np.log10(flux)
590 else:
591 # Set flux to a very faint mag if its for some reason it
592 # does not exist
593 output_array[3] = 99.
594
595 return output_array
596
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.
600
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
605 spoke length space.
606
607 Parameters
608 ----------
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.
612
613 Returns
614 -------
615 dist_tol : `float`
616 Suggested max match tolerance distance calculated from comparisons
617 between pinwheel patterns used in optimistic/pessimistic pattern
618 matcher.
619 """
620
621 self.log.debug("Starting automated tolerance calculation...")
622
623 # Create an empty array of all the patterns we possibly make
624 # sorting from brightest to faintest.
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])
629
630 # Sort our input array.
631 tmp_sort_array = cat_array[flux_args_array]
632
633 # Start making patterns.
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)
643
644 # When we store the length of each spoke in our pattern we
645 # sort from shortest to longest so we have a defined space
646 # to compare them in.
647 pattern_array[start_idx, :] = pattern_array[
648 start_idx, np.argsort(pattern_array[start_idx, :])]
649
650 # Create a searchable tree object of the patterns and find
651 # for any given pattern the closest pattern in the sorted
652 # spoke length space.
653 dist_tree = cKDTree(
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()
659
660 # We use the two closest patterns to set our tolerance.
661 dist_idx = 0
662 dist_tol = (np.degrees(dist_nearest_array[dist_idx]) * 3600.
663 / (self.config.numPointsForShape - 1.))
664
665 self.log.debug("Automated tolerance")
666 self.log.debug("\tdistance/match tol: %.4f [arcsec]", dist_tol)
667
668 return dist_tol
int min
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
Custom catalog class for record/table subclasses that are guaranteed to have an ID,...
Definition: SortedCatalog.h:42
A class representing an angle.
Definition: Angle.h:128
def matchObjectsToSources(self, refCat, sourceCat, wcs, sourceFluxField, refFluxField, match_tolerance=None)
def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources, minMatchedPairs, match_tolerance, sourceFluxField, verbose)
def __init__(self, maxMatchDist=None, autoMaxMatchDist=None, maxShift=None, lastMatchedPattern=None, failedPatternList=None, PPMbObj=None)
Lightweight representation of a geometric match between two records.
Definition: Match.h:67