LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
matchPessimisticB.py
Go to the documentation of this file.
1 
2 import numpy as np
3 from scipy.spatial import cKDTree
4 
5 import lsst.pex.config as pexConfig
6 import lsst.pipe.base as pipeBase
7 import lsst.geom as geom
8 import lsst.afw.table as afwTable
9 from lsst.utils.timer import timeMethod
10 
11 from .matchOptimisticBTask import MatchTolerance
12 
13 from .pessimistic_pattern_matcher_b_3D import PessimisticPatternMatcherB
14 
15 __all__ = ["MatchPessimisticBTask", "MatchPessimisticBConfig",
16  "MatchTolerancePessimistic"]
17 
18 
20  """Stores match tolerances for use in AstrometryTask and later
21  iterations of the matcher.
22 
23  MatchPessimisticBTask relies on several state variables to be
24  preserved over different iterations in the
25  AstrometryTask.matchAndFitWcs loop of AstrometryTask.
26 
27  Parameters
28  ----------
29  maxMatchDist : `lsst.geom.Angle`
30  Maximum distance to consider a match from the previous match/fit
31  iteration.
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
39  data.
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.
45  """
46 
47  def __init__(self, maxMatchDist=None, autoMaxMatchDist=None,
48  maxShift=None, lastMatchedPattern=None,
49  failedPatternList=None, PPMbObj=None):
50  self.maxMatchDistmaxMatchDistmaxMatchDist = maxMatchDist
51  self.autoMaxMatchDistautoMaxMatchDist = autoMaxMatchDist
52  self.maxShiftmaxShift = maxShift
53  self.lastMatchedPatternlastMatchedPattern = lastMatchedPattern
54  self.PPMbObjPPMbObj = PPMbObj
55  if failedPatternList is None:
56  self.failedPatternListfailedPatternList = []
57  else:
58  self.failedPatternListfailedPatternList = failedPatternList
59 
60 
61 class MatchPessimisticBConfig(pexConfig.Config):
62  """Configuration for MatchPessimisticBTask
63  """
64  numBrightStars = pexConfig.RangeField(
65  doc="Number of bright stars to use. Sets the max number of patterns "
66  "that can be tested.",
67  dtype=int,
68  default=200,
69  min=2,
70  )
71  minMatchedPairs = pexConfig.RangeField(
72  doc="Minimum number of matched pairs; see also minFracMatchedPairs.",
73  dtype=int,
74  default=30,
75  min=2,
76  )
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 "
81  "minMatchedPairs.",
82  dtype=float,
83  default=0.3,
84  min=0,
85  max=1,
86  )
87  matcherIterations = pexConfig.RangeField(
88  doc="Number of softening iterations in matcher.",
89  dtype=int,
90  default=5,
91  min=1,
92  )
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.",
97  dtype=int,
98  default=250,
99  max=4000,
100  )
101  maxRotationDeg = pexConfig.RangeField(
102  doc="Rotation angle allowed between sources and position reference "
103  "objects (degrees).",
104  dtype=float,
105  default=1.0,
106  max=6.0,
107  )
108  numPointsForShape = pexConfig.Field(
109  doc="Number of points to define a shape for matching.",
110  dtype=int,
111  default=6,
112  )
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.",
118  dtype=int,
119  default=6,
120  )
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.",
126  dtype=float,
127  default=1.0,
128  min=0.0,
129  max=6.0,
130  )
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 "
136  "numBrightStars.",
137  dtype=int,
138  default=3,
139  )
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.",
147  dtype=int,
148  default=1000,
149  )
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.",
153  dtype=int,
154  default=2**16,
155  min=0,
156  max=2**16 + 1,
157  )
158 
159  def validate(self):
160  pexConfig.Config.validate(self)
161  if self.numPointsForShapeAttemptnumPointsForShapeAttempt < self.numPointsForShapenumPointsForShape:
162  raise ValueError("numPointsForShapeAttempt must be greater than "
163  "or equal to numPointsForShape.")
164  if self.numPointsForShapenumPointsForShape > self.numBrightStarsnumBrightStars:
165  raise ValueError("numBrightStars must be greater than "
166  "numPointsForShape.")
167 
168 
169 # The following block adds links to this task from the Task Documentation page.
170 # \addtogroup LSST_task_documentation
171 # \{
172 # \page measAstrom_MatchPessimisticBTask
173 # \ref MatchPessimisticBTask "MatchPessimisticBTask"
174 # Match sources to reference objects
175 # \}
176 
177 
178 class MatchPessimisticBTask(pipeBase.Task):
179  """Match sources to reference objects.
180  """
181 
182  ConfigClass = MatchPessimisticBConfig
183  _DefaultName = "matchObjectsToSources"
184 
185  def __init__(self, **kwargs):
186  pipeBase.Task.__init__(self, **kwargs)
187 
188  @timeMethod
189  def matchObjectsToSources(self, refCat, sourceCat, wcs, sourceFluxField, refFluxField,
190  match_tolerance=None):
191  """Match sources to position reference stars
192 
193  refCat : `lsst.afw.table.SimpleCatalog`
194  catalog of reference objects that overlap the exposure; reads
195  fields for:
196 
197  - coord
198  - the specified flux field
199 
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`
204  estimated WCS
205  sourceFluxField: `str`
206  field of sourceCat to use for flux
207  refFluxField : `str`
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
214  fitting for the wcs.
215 
216  Returns
217  -------
218  result : `lsst.pipe.base.Struct`
219  Result struct with components:
220 
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`).
228  """
229  import lsstDebug
230  debug = lsstDebug.Info(__name__)
231 
232  # If we get an empty tolerance struct create the variables we need for
233  # this matcher.
234  if match_tolerance is None:
235  match_tolerance = MatchTolerancePessimistic()
236 
237  # Make a name alias here for consistency with older code, and to make
238  # it clear that this is a good/usable (cleaned) source catalog.
239  goodSourceCat = sourceCat
240 
241  numUsableSources = len(goodSourceCat)
242 
243  if len(goodSourceCat) == 0:
244  raise pipeBase.TaskError("No sources are good")
245 
246  minMatchedPairs = min(self.config.minMatchedPairs,
247  int(self.config.minFracMatchedPairs
248  * min([len(refCat), len(goodSourceCat)])))
249 
250  if len(refCat) > self.config.maxRefObjects:
251  self.log.warning(
252  "WARNING: Reference catalog larger that maximum allowed. "
253  "Trimming to %i", self.config.maxRefObjects)
254  trimmedRefCat = self._filterRefCat_filterRefCat(refCat, refFluxField)
255  else:
256  trimmedRefCat = refCat
257 
258  doMatchReturn = self._doMatch_doMatch(
259  refCat=trimmedRefCat,
260  sourceCat=goodSourceCat,
261  wcs=wcs,
262  refFluxField=refFluxField,
263  numUsableSources=numUsableSources,
264  minMatchedPairs=minMatchedPairs,
265  match_tolerance=match_tolerance,
266  sourceFluxField=sourceFluxField,
267  verbose=debug.verbose,
268  )
269  matches = doMatchReturn.matches
270  match_tolerance = doMatchReturn.match_tolerance
271 
272  if len(matches) == 0:
273  raise RuntimeError("Unable to match sources")
274 
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")
278 
279  return pipeBase.Struct(
280  matches=matches,
281  usableSourceCat=goodSourceCat,
282  match_tolerance=match_tolerance,
283  )
284 
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.
288 
289  No trimming is done if len(refCat) > config.maxRefObjects.
290 
291  Parameters
292  ----------
293  refCat : `lsst.afw.table.SimpleCatalog`
294  Catalog of reference objects to trim.
295  refFluxField : `str`
296  field of refCat to use for flux
297  Returns
298  -------
299  outCat : `lsst.afw.table.SimpleCatalog`
300  Catalog trimmed to the number set in the task config from the
301  brightest flux down.
302  """
303  # Find the flux cut that gives us the desired number of objects.
304  if len(refCat) <= self.config.maxRefObjects:
305  return refCat
306  fluxArray = refCat.get(refFluxField)
307  sortedFluxArray = fluxArray[fluxArray.argsort()]
308  minFlux = sortedFluxArray[-(self.config.maxRefObjects + 1)]
309 
310  selected = (refCat.get(refFluxField) > minFlux)
311 
312  outCat = afwTable.SimpleCatalog(refCat.schema)
313  outCat.reserve(self.config.maxRefObjects)
314  outCat.extend(refCat[selected])
315 
316  return outCat
317 
318  @timeMethod
319  def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources,
320  minMatchedPairs, match_tolerance, sourceFluxField, verbose):
321  """Implementation of matching sources to position reference objects
322 
323  Unlike matchObjectsToSources, this method does not check if the sources
324  are suitable.
325 
326  Parameters
327  ----------
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
334  refFluxField : `str`
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.
346  verbose : `bool`
347  Set true to print diagnostic information to std::cout
348 
349  Returns
350  -------
351  result :
352  Results struct with components:
353 
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`)
358  """
359 
360  # Load the source and reference catalog as spherical points
361  # in numpy array. We do this rather than relying on internal
362  # lsst C objects for simplicity and because we require
363  # objects contiguous in memory. We need to do these slightly
364  # differently for the reference and source cats as they are
365  # different catalog objects with different fields.
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, :] = \
373  self._latlong_flux_to_xyz_mag_latlong_flux_to_xyz_mag(theta, phi, flux)
374 
375  if match_tolerance.PPMbObj is None or \
376  match_tolerance.autoMaxMatchDist is None:
377  # The reference catalog is fixed per AstrometryTask so we only
378  # create the data needed if this is the first step in the match
379  # fit cycle.
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, :] = \
386  self._latlong_flux_to_xyz_mag_latlong_flux_to_xyz_mag(theta, phi, flux)
387  # Create our matcher object.
388  match_tolerance.PPMbObj = PessimisticPatternMatcherB(
389  ref_array[:, :3], self.log)
390  self.log.debug("Computing source statistics...")
391  maxMatchDistArcSecSrc = self._get_pair_pattern_statistics_get_pair_pattern_statistics(
392  src_array)
393  self.log.debug("Computing reference statistics...")
394  maxMatchDistArcSecRef = self._get_pair_pattern_statistics_get_pair_pattern_statistics(
395  ref_array)
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)
403 
404  # Set configurable defaults when we encounter None type or set
405  # state based on previous run of AstrometryTask._matchAndFitWcs.
406  if match_tolerance.maxShift is None:
407  maxShiftArcseconds = (self.config.maxOffsetPix
408  * wcs.getPixelScale().asArcseconds())
409  else:
410  # We don't want to clamp down too hard on the allowed shift so
411  # we test that the smallest we ever allow is the pixel scale.
412  maxShiftArcseconds = np.max(
413  (match_tolerance.maxShift.asArcseconds(),
414  self.config.minMatchDistPixels
415  * wcs.getPixelScale().asArcseconds()))
416 
417  # If our tolerances are not set from a previous run, estimate a
418  # starting tolerance guess from the statistics of patterns we can
419  # create on both the source and reference catalog. We use the smaller
420  # of the two.
421  if match_tolerance.maxMatchDist is None:
422  match_tolerance.maxMatchDist = match_tolerance.autoMaxMatchDist
423  else:
424  maxMatchDistArcSec = np.max(
425  (self.config.minMatchDistPixels
426  * wcs.getPixelScale().asArcseconds(),
427  np.min((match_tolerance.maxMatchDist.asArcseconds(),
428  match_tolerance.autoMaxMatchDist.asArcseconds()))))
429 
430  # Make sure the data we are considering is dense enough to require
431  # the consensus mode of the matcher. If not default to Optimistic
432  # pattern matcher behavior. We enforce pessimistic mode if the
433  # reference cat is sufficiently large, avoiding false positives.
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:
441  numConsensus = 1
442 
443  self.log.debug("Current tol maxDist: %.4f arcsec",
444  maxMatchDistArcSec)
445  self.log.debug("Current shift: %.4f arcsec",
446  maxShiftArcseconds)
447 
448  match_found = False
449  # Start the iteration over our tolerances.
450  for soften_dist in range(self.config.matcherIterations):
451  if soften_dist == 0 and \
452  match_tolerance.lastMatchedPattern is not None:
453  # If we are on the first, most stringent tolerance,
454  # and have already found a match, the matcher should behave
455  # like an optimistic pattern matcher. Exiting at the first
456  # match.
457  run_n_consent = 1
458  else:
459  # If we fail or this is the first match attempt, set the
460  # pattern consensus to the specified config value.
461  run_n_consent = numConsensus
462  # We double the match dist tolerance each round and add 1 to the
463  # to the number of candidate spokes to check.
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)
476  )
477 
478  if soften_dist == 0 and \
479  len(matcher_struct.match_ids) == 0 and \
480  match_tolerance.lastMatchedPattern is not None:
481  # If we found a pattern on a previous match-fit iteration and
482  # can't find an optimistic match on our first try with the
483  # tolerances as found in the previous match-fit,
484  # the match we found in the last iteration was likely bad. We
485  # append the bad match's index to the a list of
486  # patterns/matches to skip on subsequent iterations.
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:
493  # Match found, save a bit a state regarding this pattern
494  # in the match tolerance class object and exit.
495  match_tolerance.maxShift = \
496  matcher_struct.shift * geom.arcseconds
497  match_tolerance.lastMatchedPattern = \
498  matcher_struct.pattern_idx
499  match_found = True
500  break
501 
502  # If we didn't find a match, exit early.
503  if not match_found:
504  return pipeBase.Struct(
505  matches=[],
506  match_tolerance=match_tolerance,
507  )
508 
509  # The matcher returns all the nearest neighbors that agree between
510  # the reference and source catalog. For the current astrometric solver
511  # we need to remove as many false positives as possible before sending
512  # the matches off to the solver. The low value of 100 and high value of
513  # 2 are the low number of sigma and high respectively. The exact values
514  # were found after testing on data of various reference/source
515  # densities and astrometric distortion quality, specifically the
516  # visits: HSC (3358), DECam (406285, 410827),
517  # CFHT (793169, 896070, 980526).
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()))
522 
523  # A match has been found, return our list of matches and
524  # return.
525  matches = []
526  for match_id_pair, dist_arcsec in zip(matcher_struct.match_ids,
527  distances_arcsec):
528  if dist_arcsec < dist_cut_arcsec:
529  match = afwTable.ReferenceMatch()
530  match.first = refCat[int(match_id_pair[1])]
531  match.second = sourceCat[int(match_id_pair[0])]
532  # We compute the true distance along and sphere. This isn't
533  # used in the WCS fitter however it is used in the unittest
534  # to confirm the matches computed.
535  match.distance = match.first.getCoord().separation(
536  match.second.getCoord()).asArcseconds()
537  matches.append(match)
538 
539  return pipeBase.Struct(
540  matches=matches,
541  match_tolerance=match_tolerance,
542  )
543 
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.
547 
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.
551 
552  Parameters
553  ----------
554  theta : `float`
555  Angle from the north pole (z axis) of the sphere
556  phi : `float`
557  Rotation around the sphere
558 
559  Return
560  ------
561  output_array : `numpy.ndarray`, (N, 4)
562  Spherical unit vector x, y, z with flux.
563  """
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)
568  if flux > 0:
569  output_array[3] = -2.5 * np.log10(flux)
570  else:
571  # Set flux to a very faint mag if its for some reason it
572  # does not exist
573  output_array[3] = 99.
574 
575  return output_array
576 
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.
580 
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
585  spoke length space.
586 
587  Parameters
588  ----------
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.
592 
593  Returns
594  -------
595  dist_tol : `float`
596  Suggested max match tolerance distance calculated from comparisons
597  between pinwheel patterns used in optimistic/pessimistic pattern
598  matcher.
599  """
600 
601  self.log.debug("Starting automated tolerance calculation...")
602 
603  # Create an empty array of all the patterns we possibly make
604  # sorting from brightest to faintest.
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])
609 
610  # Sort our input array.
611  tmp_sort_array = cat_array[flux_args_array]
612 
613  # Start making patterns.
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)
623 
624  # When we store the length of each spoke in our pattern we
625  # sort from shortest to longest so we have a defined space
626  # to compare them in.
627  pattern_array[start_idx, :] = pattern_array[
628  start_idx, np.argsort(pattern_array[start_idx, :])]
629 
630  # Create a searchable tree object of the patterns and find
631  # for any given pattern the closest pattern in the sorted
632  # spoke length space.
633  dist_tree = cKDTree(
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()
639 
640  # We use the two closest patterns to set our tolerance.
641  dist_idx = 0
642  dist_tol = (np.degrees(dist_nearest_array[dist_idx]) * 3600.
643  / (self.config.numPointsForShape - 1.))
644 
645  self.log.debug("Automated tolerance")
646  self.log.debug("\tdistance/match tol: %.4f [arcsec]", dist_tol)
647 
648  return dist_tol
int min
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:127
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