Coverage for python / lsst / meas / astrom / matchPessimisticB.py: 16%
180 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-14 23:59 +0000
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/>.
22__all__ = ["MatchPessimisticBTask", "MatchPessimisticBConfig",
23 "MatchTolerancePessimistic"]
25import numpy as np
26from scipy.spatial import cKDTree
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
34from . import exceptions
35from .matchOptimisticBTask import MatchTolerance
36from .pessimistic_pattern_matcher_b_3D import PessimisticPatternMatcherB
39class MatchTolerancePessimistic(MatchTolerance):
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.
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.
62 PPMbObj : `lsst.meas.astrom.PessimisticPatternMatcherB`
63 Initialized Pessimistic pattern matcher object. Storing this prevents
64 the need for recalculation of the searchable distances in the PPMB.
65 """
67 def __init__(self, maxMatchDist=None, autoMaxMatchDist=None,
68 maxShift=None, lastMatchedPattern=None,
69 failedPatternList=None, PPMbObj=None):
70 self.maxMatchDist = maxMatchDist
71 self.autoMaxMatchDist = autoMaxMatchDist
72 self.maxShift = maxShift
73 self.lastMatchedPattern = lastMatchedPattern
74 self.PPMbObj = PPMbObj
75 if failedPatternList is None:
76 self.failedPatternList = []
77 else:
78 self.failedPatternList = failedPatternList
81class MatchPessimisticBConfig(pexConfig.Config):
82 """Configuration for MatchPessimisticBTask
83 """
84 numBrightStars = pexConfig.RangeField(
85 doc="Maximum number of bright stars to use. Sets the max number of patterns "
86 "that can be tested.",
87 dtype=int,
88 default=150,
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**13,
175 min=0,
176 max=2**16 + 1,
177 )
179 def validate(self):
180 pexConfig.Config.validate(self)
181 if self.numPointsForShapeAttempt < self.numPointsForShape:
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.")
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# \}
198class MatchPessimisticBTask(pipeBase.Task):
199 """Match sources to reference objects.
200 """
202 ConfigClass = MatchPessimisticBConfig
203 _DefaultName = "matchPessimisticB"
205 def __init__(self, **kwargs):
206 pipeBase.Task.__init__(self, **kwargs)
208 @timeMethod
209 def matchObjectsToSources(self, refCat, sourceCat, wcs, sourceFluxField, refFluxField,
210 matchTolerance=None, bbox=None):
211 """Match sources to position reference stars
213 refCat : `lsst.afw.table.SimpleCatalog`
214 catalog of reference objects that overlap the exposure; reads
215 fields for:
217 - coord
218 - the specified flux field
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.
223 wcs : `lsst.afw.geom.SkyWcs`
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
229 matchTolerance : `lsst.meas.astrom.MatchTolerancePessimistic`
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 bbox : `lsst.geom.Box2I`, optional
236 Bounding box of the exposure for evaluating the local pixelScale
237 (defaults to the Sky Origin of the ``wcs`` provided if ``bbox``
238 is `None`).
240 Returns
241 -------
242 result : `lsst.pipe.base.Struct`
243 Result struct with components:
245 - ``matches`` : source to reference matches found (`list` of
246 `lsst.afw.table.ReferenceMatch`)
247 - ``usableSourceCat`` : a catalog of sources potentially usable for
248 matching and WCS fitting (`lsst.afw.table.SourceCatalog`).
249 - ``matchTolerance`` : a MatchTolerance object containing the
250 resulting state variables from the match
251 (`lsst.meas.astrom.MatchTolerancePessimistic`).
252 """
253 import lsstDebug
254 debug = lsstDebug.Info(__name__)
256 # If we get an empty tolerance struct create the variables we need for
257 # this matcher.
258 if matchTolerance is None:
259 matchTolerance = MatchTolerancePessimistic()
261 # Make a name alias here for consistency with older code, and to make
262 # it clear that this is a good/usable (cleaned) source catalog.
263 goodSourceCat = sourceCat
265 if (numUsableSources := len(goodSourceCat)) == 0:
266 raise exceptions.MatcherFailure("No sources are good")
268 minMatchedPairs = min(self.config.minMatchedPairs,
269 int(self.config.minFracMatchedPairs
270 * min([len(refCat), len(goodSourceCat)])))
272 if len(goodSourceCat) <= self.config.numPointsForShape:
273 msg = (f"Not enough catalog objects ({len(goodSourceCat)}) to make a "
274 f"shape for the matcher (need {self.config.numPointsForShape}).")
275 raise exceptions.MatcherFailure(msg)
276 if len(refCat) <= self.config.numPointsForShape:
277 msg = (f"Not enough refcat objects ({len(refCat)}) to make a "
278 f"shape for the matcher (need {self.config.numPointsForShape}).")
279 raise exceptions.MatcherFailure(msg)
281 if len(refCat) > self.config.maxRefObjects:
282 self.log.warning(
283 "Reference catalog larger than maximum allowed. "
284 "Trimming to %i", self.config.maxRefObjects)
285 trimmedRefCat = self._filterRefCat(refCat, refFluxField)
286 else:
287 trimmedRefCat = refCat
289 doMatchReturn = self._doMatch(
290 refCat=trimmedRefCat,
291 sourceCat=goodSourceCat,
292 wcs=wcs,
293 refFluxField=refFluxField,
294 numUsableSources=numUsableSources,
295 minMatchedPairs=minMatchedPairs,
296 matchTolerance=matchTolerance,
297 sourceFluxField=sourceFluxField,
298 verbose=debug.verbose,
299 bbox=bbox,
300 )
301 matches = doMatchReturn.matches
302 matchTolerance = doMatchReturn.matchTolerance
304 if (nMatches := len(matches)) == 0:
305 raise exceptions.MatcherFailure("No matches found")
307 self.log.info("Matched %d sources", nMatches)
308 if nMatches < minMatchedPairs:
309 self.log.warning("Number of matches (%s) is smaller than minimum requested (%s)",
310 nMatches, minMatchedPairs)
312 return pipeBase.Struct(
313 matches=matches,
314 usableSourceCat=goodSourceCat,
315 matchTolerance=matchTolerance,
316 )
318 def _filterRefCat(self, refCat, refFluxField):
319 """Sub-select a number of reference objects starting from the brightest
320 and maxing out at the number specified by maxRefObjects in the config.
322 No trimming is done if len(refCat) > config.maxRefObjects.
324 Parameters
325 ----------
326 refCat : `lsst.afw.table.SimpleCatalog`
327 Catalog of reference objects to trim.
328 refFluxField : `str`
329 field of refCat to use for flux
330 Returns
331 -------
332 outCat : `lsst.afw.table.SimpleCatalog`
333 Catalog trimmed to the number set in the task config from the
334 brightest flux down.
335 """
336 # Find the flux cut that gives us the desired number of objects.
337 if len(refCat) <= self.config.maxRefObjects:
338 return refCat
339 fluxArray = refCat.get(refFluxField)
340 sortedFluxArray = fluxArray[fluxArray.argsort()]
341 minFlux = sortedFluxArray[-(self.config.maxRefObjects + 1)]
343 selected = (refCat.get(refFluxField) > minFlux)
345 outCat = afwTable.SimpleCatalog(refCat.schema)
346 outCat.reserve(self.config.maxRefObjects)
347 outCat.extend(refCat[selected])
349 return outCat
351 @timeMethod
352 def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources,
353 minMatchedPairs, matchTolerance, sourceFluxField, verbose, bbox=None):
354 """Implementation of matching sources to position reference objects
356 Unlike matchObjectsToSources, this method does not check if the sources
357 are suitable.
359 Parameters
360 ----------
361 refCat : `lsst.afw.table.SimpleCatalog`
362 catalog of position reference objects that overlap an exposure
363 sourceCat : `lsst.afw.table.SourceCatalog`
364 catalog of sources found on the exposure
365 wcs : `lsst.afw.geom.SkyWcs`
366 estimated WCS of exposure
367 refFluxField : `str`
368 field of refCat to use for flux
369 numUsableSources : `int`
370 number of usable sources (sources with known centroid that are not
371 near the edge, but may be saturated)
372 minMatchedPairs : `int`
373 minimum number of matches
374 matchTolerance : `lsst.meas.astrom.MatchTolerancePessimistic`
375 a MatchTolerance object containing variables specifying matcher
376 tolerances and state from possible previous runs.
377 sourceFluxField : `str`
378 Name of the flux field in the source catalog.
379 verbose : `bool`
380 Set true to print diagnostic information to std::cout
381 bbox : `lsst.geom.Box2I`, optional
382 Bounding box of the exposure for evaluating the local pixelScale
383 (defaults to the Sky Origin of the ``wcs`` provided if ``bbox``
384 is `None`).
386 Returns
387 -------
388 result : `lsst.pipe.base.Struct`
389 Results struct with components:
391 - ``matches`` : a list the matches found
392 (`list` of `lsst.afw.table.ReferenceMatch`).
393 - ``matchTolerance`` : MatchTolerance containing updated values from
394 this fit iteration (`lsst.meas.astrom.MatchTolerancePessimistic`)
395 """
396 if bbox is not None:
397 pixelScale = wcs.getPixelScale(bbox.getCenter()).asArcseconds()
398 else:
399 pixelScale = wcs.getPixelScale().asArcseconds()
401 # Load the source and reference catalog as spherical points
402 # in numpy array. We do this rather than relying on internal
403 # lsst C objects for simplicity and because we require
404 # objects contiguous in memory. We need to do these slightly
405 # differently for the reference and source cats as they are
406 # different catalog objects with different fields.
407 src_array = np.empty((len(sourceCat), 4), dtype=np.float64)
408 for src_idx, srcObj in enumerate(sourceCat):
409 coord = wcs.pixelToSky(srcObj.getCentroid())
410 theta = np.pi / 2 - coord.getLatitude().asRadians()
411 phi = coord.getLongitude().asRadians()
412 flux = srcObj[sourceFluxField]
413 src_array[src_idx, :] = \
414 self._latlong_flux_to_xyz_mag(theta, phi, flux)
416 if matchTolerance.PPMbObj is None or \
417 matchTolerance.autoMaxMatchDist is None:
418 # The reference catalog is fixed per AstrometryTask so we only
419 # create the data needed if this is the first step in the match
420 # fit cycle.
421 ref_array = np.empty((len(refCat), 4), dtype=np.float64)
422 for ref_idx, refObj in enumerate(refCat):
423 theta = np.pi / 2 - refObj.getDec().asRadians()
424 phi = refObj.getRa().asRadians()
425 flux = refObj[refFluxField]
426 ref_array[ref_idx, :] = \
427 self._latlong_flux_to_xyz_mag(theta, phi, flux)
428 # Create our matcher object.
429 matchTolerance.PPMbObj = PessimisticPatternMatcherB(
430 ref_array[:, :3], self.log)
431 self.log.debug("Computing source statistics...")
432 maxMatchDistArcSecSrc = self._get_pair_pattern_statistics(
433 src_array)
434 self.log.debug("Computing reference statistics...")
435 maxMatchDistArcSecRef = self._get_pair_pattern_statistics(
436 ref_array)
437 maxMatchDistArcSec = np.max((
438 self.config.minMatchDistPixels * pixelScale,
439 np.min((maxMatchDistArcSecSrc,
440 maxMatchDistArcSecRef))))
441 matchTolerance.autoMaxMatchDist = geom.Angle(
442 maxMatchDistArcSec, geom.arcseconds)
444 # Set configurable defaults when we encounter None type or set
445 # state based on previous run of AstrometryTask._matchAndFitWcs.
446 if matchTolerance.maxShift is None:
447 maxShiftArcseconds = self.config.maxOffsetPix * pixelScale
448 else:
449 # We don't want to clamp down too hard on the allowed shift so
450 # we test that the smallest we ever allow is the pixel scale.
451 maxShiftArcseconds = np.max(
452 (matchTolerance.maxShift.asArcseconds(),
453 self.config.minMatchDistPixels * pixelScale))
455 # If our tolerances are not set from a previous run, estimate a
456 # starting tolerance guess from the statistics of patterns we can
457 # create on both the source and reference catalog. We use the smaller
458 # of the two.
459 if matchTolerance.maxMatchDist is None:
460 matchTolerance.maxMatchDist = matchTolerance.autoMaxMatchDist
461 else:
462 maxMatchDistArcSec = np.max(
463 (self.config.minMatchDistPixels * pixelScale,
464 np.min((matchTolerance.maxMatchDist.asArcseconds(),
465 matchTolerance.autoMaxMatchDist.asArcseconds()))))
467 # Make sure the data we are considering is dense enough to require
468 # the consensus mode of the matcher. If not default to Optimistic
469 # pattern matcher behavior. We enforce pessimistic mode if the
470 # reference cat is sufficiently large, avoiding false positives.
471 numConsensus = self.config.numPatternConsensus
472 if len(refCat) < self.config.numRefRequireConsensus:
473 minObjectsForConsensus = \
474 self.config.numBrightStars + \
475 self.config.numPointsForShapeAttempt
476 if len(refCat) < minObjectsForConsensus or \
477 len(sourceCat) < minObjectsForConsensus:
478 numConsensus = 1
480 self.log.debug("Current tol maxDist: %.4f arcsec",
481 maxMatchDistArcSec)
482 self.log.debug("Current shift: %.4f arcsec",
483 maxShiftArcseconds)
485 match_found = False
486 # Start the iteration over our tolerances.
487 for soften_dist in range(self.config.matcherIterations):
488 if soften_dist == 0 and \
489 matchTolerance.lastMatchedPattern is not None:
490 # If we are on the first, most stringent tolerance,
491 # and have already found a match, the matcher should behave
492 # like an optimistic pattern matcher. Exiting at the first
493 # match.
494 run_n_consent = 1
495 else:
496 # If we fail or this is the first match attempt, set the
497 # pattern consensus to the specified config value.
498 run_n_consent = numConsensus
499 # We double the match dist tolerance each round and add 1 to the
500 # to the number of candidate spokes to check.
501 matcher_struct = matchTolerance.PPMbObj.match(
502 source_array=src_array,
503 n_check=self.config.numPointsForShapeAttempt,
504 n_match=self.config.numPointsForShape,
505 n_agree=run_n_consent,
506 max_n_patterns=self.config.numBrightStars,
507 max_shift=maxShiftArcseconds,
508 max_rotation=self.config.maxRotationDeg,
509 max_dist=maxMatchDistArcSec * 2. ** soften_dist,
510 min_matches=minMatchedPairs,
511 pattern_skip_array=np.array(
512 matchTolerance.failedPatternList)
513 )
515 if soften_dist == 0 and \
516 len(matcher_struct.match_ids) == 0 and \
517 matchTolerance.lastMatchedPattern is not None:
518 # If we found a pattern on a previous match-fit iteration and
519 # can't find an optimistic match on our first try with the
520 # tolerances as found in the previous match-fit,
521 # the match we found in the last iteration was likely bad. We
522 # append the bad match's index to the a list of
523 # patterns/matches to skip on subsequent iterations.
524 matchTolerance.failedPatternList.append(
525 matchTolerance.lastMatchedPattern)
526 matchTolerance.lastMatchedPattern = None
527 maxShiftArcseconds = self.config.maxOffsetPix * pixelScale
528 elif len(matcher_struct.match_ids) > 0:
529 # Match found, save a bit a state regarding this pattern
530 # in the match tolerance class object and exit.
531 matchTolerance.maxShift = \
532 matcher_struct.shift * geom.arcseconds
533 matchTolerance.lastMatchedPattern = \
534 matcher_struct.pattern_idx
535 match_found = True
536 break
538 # If we didn't find a match, exit early.
539 if not match_found:
540 return pipeBase.Struct(
541 matches=[],
542 matchTolerance=matchTolerance,
543 )
545 # The matcher returns all the nearest neighbors that agree between
546 # the reference and source catalog. For the current astrometric solver
547 # we need to remove as many false positives as possible before sending
548 # the matches off to the solver. The low value of 100 and high value of
549 # 2 are the low number of sigma and high respectively. The exact values
550 # were found after testing on data of various reference/source
551 # densities and astrometric distortion quality, specifically the
552 # visits: HSC (3358), DECam (406285, 410827),
553 # CFHT (793169, 896070, 980526).
554 distances_arcsec = np.degrees(matcher_struct.distances_rad) * 3600
555 dist_cut_arcsec = np.max(
556 (np.degrees(matcher_struct.max_dist_rad) * 3600,
557 self.config.minMatchDistPixels * pixelScale))
559 # A match has been found, return our list of matches and
560 # return.
561 matches = []
562 for match_id_pair, dist_arcsec in zip(matcher_struct.match_ids,
563 distances_arcsec):
564 if dist_arcsec < dist_cut_arcsec:
565 match = afwTable.ReferenceMatch()
566 match.first = refCat[int(match_id_pair[1])]
567 match.second = sourceCat[int(match_id_pair[0])]
568 # We compute the true distance along and sphere. This isn't
569 # used in the WCS fitter however it is used in the unittest
570 # to confirm the matches computed.
571 match.distance = match.first.getCoord().separation(
572 match.second.getCoord()).asArcseconds()
573 matches.append(match)
575 return pipeBase.Struct(
576 matches=matches,
577 matchTolerance=matchTolerance,
578 )
580 def _latlong_flux_to_xyz_mag(self, theta, phi, flux):
581 """Convert angles theta and phi and a flux into unit sphere
582 x, y, z, and a relative magnitude.
584 Takes in a afw catalog object and converts the catalog object RA, DECs
585 to points on the unit sphere. Also converts the flux into a simple,
586 non-zero-pointed magnitude for relative sorting.
588 Parameters
589 ----------
590 theta : `float`
591 Angle from the north pole (z axis) of the sphere
592 phi : `float`
593 Rotation around the sphere
595 Return
596 ------
597 output_array : `numpy.ndarray`, (N, 4)
598 Spherical unit vector x, y, z with flux.
599 """
600 output_array = np.empty(4, dtype=np.float64)
601 output_array[0] = np.sin(theta)*np.cos(phi)
602 output_array[1] = np.sin(theta)*np.sin(phi)
603 output_array[2] = np.cos(theta)
604 if flux > 0:
605 output_array[3] = -2.5 * np.log10(flux)
606 else:
607 # Set flux to a very faint mag if its for some reason it
608 # does not exist
609 output_array[3] = 99.
611 return output_array
613 def _get_pair_pattern_statistics(self, cat_array):
614 """ Compute the tolerances for the matcher automatically by comparing
615 pinwheel patterns as we would in the matcher.
617 We test how similar the patterns we can create from a given set of
618 objects by computing the spoke lengths for each pattern and sorting
619 them from smallest to largest. The match tolerance is the average
620 distance per spoke between the closest two patterns in the sorted
621 spoke length space.
623 Parameters
624 ----------
625 cat_array : `numpy.ndarray`, (N, 3)
626 array of 3 vectors representing the x, y, z position of catalog
627 objects on the unit sphere.
629 Returns
630 -------
631 dist_tol : `float`
632 Suggested max match tolerance distance calculated from comparisons
633 between pinwheel patterns used in optimistic/pessimistic pattern
634 matcher.
635 """
637 self.log.debug("Starting automated tolerance calculation...")
639 # Create an empty array of all the patterns we possibly make
640 # sorting from brightest to faintest.
641 pattern_array = np.empty(
642 (cat_array.shape[0] - self.config.numPointsForShape,
643 self.config.numPointsForShape - 1))
644 flux_args_array = np.argsort(cat_array[:, -1])
646 # Sort our input array.
647 tmp_sort_array = cat_array[flux_args_array]
649 # Start making patterns.
650 for start_idx in range(cat_array.shape[0]
651 - self.config.numPointsForShape):
652 pattern_points = tmp_sort_array[start_idx:start_idx
653 + self.config.numPointsForShape, :-1]
654 pattern_delta = pattern_points[1:, :] - pattern_points[0, :]
655 pattern_array[start_idx, :] = np.sqrt(
656 pattern_delta[:, 0] ** 2
657 + pattern_delta[:, 1] ** 2
658 + pattern_delta[:, 2] ** 2)
660 # When we store the length of each spoke in our pattern we
661 # sort from shortest to longest so we have a defined space
662 # to compare them in.
663 pattern_array[start_idx, :] = pattern_array[
664 start_idx, np.argsort(pattern_array[start_idx, :])]
666 # Create a searchable tree object of the patterns and find
667 # for any given pattern the closest pattern in the sorted
668 # spoke length space.
669 dist_tree = cKDTree(
670 pattern_array[:, :(self.config.numPointsForShape - 1)])
671 dist_nearest_array, ids = dist_tree.query(
672 pattern_array[:, :(self.config.numPointsForShape - 1)], k=2)
673 dist_nearest_array = dist_nearest_array[:, 1]
674 dist_nearest_array.sort()
676 # We use the two closest patterns to set our tolerance.
677 dist_idx = 0
678 dist_tol = (np.degrees(dist_nearest_array[dist_idx]) * 3600.
679 / (self.config.numPointsForShape - 1.))
681 self.log.debug("Automated tolerance")
682 self.log.debug("\tdistance/match tol: %.4f [arcsec]", dist_tol)
684 return dist_tol