LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
matchOptimisticB.py
Go to the documentation of this file.
1 from __future__ import absolute_import, division, print_function
2 import math
3 
4 import numpy as np
5 
6 import lsst.afw.table as afwTable
7 import lsst.pex.config as pexConfig
8 import lsst.pipe.base as pipeBase
9 from .setMatchDistance import setMatchDistance
10 from .astromLib import matchOptimisticB, MatchOptimisticBControl
11 
12 __all__ = ["MatchOptimisticBTask", "MatchOptimisticBConfig", "SourceInfo"]
13 
14 class MatchOptimisticBConfig(pexConfig.Config):
15  """Configuration for MatchOptimisticBTask
16  """
17  sourceFluxType = pexConfig.Field(
18  doc = "Type of source flux; typically one of Ap or Psf",
19  dtype = str,
20  default = "Ap",
21  )
22  maxMatchDistArcSec = pexConfig.RangeField(
23  doc = "Maximum separation between reference objects and sources "
24  "beyond which they will not be considered a match (arcsec)",
25  dtype = float,
26  default = 3,
27  min = 0,
28  )
29  numBrightStars = pexConfig.RangeField(
30  doc = "Number of bright stars to use",
31  dtype = int,
32  default = 50,
33  min = 2,
34  )
35  minMatchedPairs = pexConfig.RangeField(
36  doc = "Minimum number of matched pairs; see also minFracMatchedPairs",
37  dtype = int,
38  default = 30,
39  min = 2,
40  )
41  minFracMatchedPairs = pexConfig.RangeField(
42  doc = "Minimum number of matched pairs as a fraction of the smaller of "
43  "the number of reference stars or the number of good sources; "
44  "the actual minimum is the smaller of this value or minMatchedPairs",
45  dtype = float,
46  default = 0.3,
47  min = 0,
48  max = 1,
49  )
50  maxOffsetPix = pexConfig.RangeField(
51  doc = "Maximum allowed shift of WCS, due to matching (pixel)",
52  dtype = int,
53  default = 300,
54  max = 4000,
55  )
56  maxRotationDeg = pexConfig.RangeField(
57  doc = "Rotation angle allowed between sources and position reference objects (degrees)",
58  dtype = float,
59  default = 1.0,
60  max = 6.0,
61  )
62  allowedNonperpDeg = pexConfig.RangeField(
63  doc = "Allowed non-perpendicularity of x and y (degree)",
64  dtype = float,
65  default = 3.0,
66  max = 45.0,
67  )
68  numPointsForShape = pexConfig.Field(
69  doc = "number of points to define a shape for matching",
70  dtype = int,
71  default = 6,
72  )
73  maxDeterminant = pexConfig.Field(
74  doc = "maximum determinant of linear transformation matrix for a usable solution",
75  dtype = float,
76  default = 0.02,
77  )
78 
79 class SourceInfo(object):
80  """Provide usability tests and catalog keys for sources in a source catalog
81 
82  Fields set include:
83  - centroidKey key for centroid
84  - centroidFlagKey key for flag that is True if centroid is valid
85  - edgeKey key for field that is True if source is near an edge
86  - saturatedKey key for field that is True if source has any saturated pixels
87  - interpolatedCenterKey key for field that is True if center pixels have interpolated values;
88  interpolation is triggered by saturation, cosmic rays and bad pixels, and possibly other reasons
89  - fluxField name of flux field
90 
91  @throw RuntimeError if schema version unsupported or a needed field is not found
92  """
93  def __init__(self, schema, fluxType="Ap"):
94  """Construct a SourceInfo
95 
96  @param[in] schema source catalog schema
97  @param[in] fluxType flux type: typically one of "Ap" or "Psf"
98 
99  @throw RuntimeError if the flux field is not found
100  """
101  self.centroidKey = afwTable.Point2DKey(schema["slot_Centroid"])
102  self.centroidFlagKey = schema["slot_Centroid_flag"].asKey()
103  self.edgeKey = schema["base_PixelFlags_flag_edge"].asKey()
104  self.saturatedKey = schema["base_PixelFlags_flag_saturated"].asKey()
105  self.fluxField = "slot_%sFlux_flux" % (fluxType,)
106  self.interpolatedCenterKey = schema["base_PixelFlags_flag_interpolatedCenter"].asKey()
107  # extra keys that might be useful
108  self.parentKey = schema["parent"].asKey()
109 
110  if self.fluxField not in schema:
111  raise RuntimeError("Could not find flux field %s in source schema" % (self.fluxField,))
112 
113  def _isMultiple(self, source):
114  """Return True if source is likely multiple sources
115  """
116  if source.get(self.parentKey) != 0:
117  return True
118  footprint = source.getFootprint()
119  return footprint is not None and len(footprint.getPeaks()) > 1
120 
121  def hasCentroid(self, source):
122  """Return True if the source has a valid centroid
123  """
124  centroid = source.get(self.centroidKey)
125  return np.all(np.isfinite(centroid)) and not source.getCentroidFlag()
126 
127  def isUsable(self, source):
128  """Return True if the source is usable for matching, even if possibly saturated
129 
130  For a source to be usable it must:
131  - have a valid centroid
132  - be not too near the edge
133  """
134  return self.hasCentroid(source) and not source.get(self.edgeKey)
135 
136  def isGood(self, source):
137  """Return True if source is usable for matching (as per isUsable) and likely has a good centroid
138 
139  For a source to have a good centroid it should not be interpolated in the center;
140  this includes saturated sources so we don't have to test separately for that.
141  """
142  return self.isUsable(source) and not source.get(self.interpolatedCenterKey)
143 
144 
145 # The following block adds links to this task from the Task Documentation page.
146 ## \addtogroup LSST_task_documentation
147 ## \{
148 ## \page measAstrom_matchOptimisticBTask
149 ## \ref MatchOptimisticBTask "MatchOptimisticBTask"
150 ## Match sources to reference objects
151 ## \}
152 
153 class MatchOptimisticBTask(pipeBase.Task):
154  """!Match sources to reference objects
155 
156  @anchor MatchOptimisticBTask_
157 
158  @section meas_astrom_matchOptimisticB_Contents Contents
159 
160  - @ref meas_astrom_matchOptimisticB_Purpose
161  - @ref meas_astrom_matchOptimisticB_Initialize
162  - @ref meas_astrom_matchOptimisticB_IO
163  - @ref meas_astrom_matchOptimisticB_Config
164  - @ref meas_astrom_matchOptimisticB_Example
165  - @ref meas_astrom_matchOptimisticB_Debug
166 
167  @section meas_astrom_matchOptimisticB_Purpose Description
168 
169  Match sources to reference objects. This is often done as a preliminary step to fitting an astrometric
170  or photometric solution. For details about the matching algorithm see matchOptimisticB.h
171 
172  @section meas_astrom_matchOptimisticB_Initialize Task initialisation
173 
174  @copydoc \_\_init\_\_
175 
176  @section meas_astrom_matchOptimisticB_IO Invoking the Task
177 
178  @copydoc matchObjectsToSources
179 
180  @section meas_astrom_matchOptimisticB_Config Configuration parameters
181 
182  See @ref MatchOptimisticBConfig
183 
184  To modify the tests for usable sources and good sources, subclass SourceInfo and
185  set MatchOptimisticBTask.SourceInfoClass to your subclass.
186 
187  @section meas_astrom_matchOptimisticB_Example A complete example of using MatchOptimisticBTask
188 
189  MatchOptimisticBTask is a subtask of AstrometryTask, which is called by PhotoCalTask.
190  See \ref meas_photocal_photocal_Example.
191 
192  @section meas_astrom_matchOptimisticB_Debug Debug variables
193 
194  The @link lsst.pipe.base.cmdLineTask.CmdLineTask command line task@endlink interface supports a
195  flag @c -d to import @b debug.py from your @c PYTHONPATH; see @ref baseDebug for more about
196  @b debug.py files.
197 
198  The available variables in MatchOptimisticBTask are:
199  <DL>
200  <DT> @c verbose (bool)
201  <DD> If True then the matcher prints debug messages to stdout
202  </DL>
203 
204  To investigate the @ref meas_astrom_matchOptimisticB_Debug, put something like
205  @code{.py}
206  import lsstDebug
207  def DebugInfo(name):
208  debug = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
209  if name == "lsst.pipe.tasks.astrometry":
210  debug.verbose = True
211 
212  return debug
213 
214  lsstDebug.Info = DebugInfo
215  @endcode
216  into your debug.py file and run this task with the @c --debug flag.
217  """
218  ConfigClass = MatchOptimisticBConfig
219  _DefaultName = "matchObjectsToSources"
220  SourceInfoClass = SourceInfo
221 
222  def filterStars(self, refCat):
223  """Extra filtering pass; subclass if desired
224  """
225  return refCat
226 
227  @pipeBase.timeMethod
228  def matchObjectsToSources(self, refCat, sourceCat, wcs, refFluxField, maxMatchDistArcSec=None):
229  """!Match sources to position reference stars
230 
231  @param[in] refCat catalog of reference objects that overlap the exposure; reads fields for:
232  - coord
233  - the specified flux field
234  @param[in] sourceCat catalog of sources found on an exposure; reads fields for:
235  - centroid
236  - centroid flag
237  - edge flag
238  - saturated flag
239  - aperture flux, if found, else PSF flux
240  @param[in] wcs estimated WCS
241  @param[in] refFluxField field of refCat to use for flux
242  @param[in] maxMatchDistArcSec maximum distance between reference objects and sources (arcsec);
243  if specified then min(config.maxMatchDistArcSec, maxMatchDistArcSec) is used
244  if None then config.maxMatchDistArcSec is used
245  @return an lsst.pipe.base.Struct with fields:
246  - matches a list of matches, each instance of lsst.afw.table.ReferenceMatch
247  - usableSourcCat a catalog of sources potentially usable for matching.
248  For this fitter usable sources include unresolved sources not too near the edge.
249  It includes saturated sources, even those these are removed from the final match list,
250  because saturated sources may be used to determine the match list.
251  """
252  import lsstDebug
253  debug = lsstDebug.Info(__name__)
254 
255  preNumObj = len(refCat)
256  refCat = self.filterStars(refCat)
257  numRefObj = len(refCat)
258 
259  if self.log: self.log.info("filterStars purged %d reference stars, leaving %d stars" % \
260  (preNumObj - numRefObj, numRefObj))
261 
262  sourceInfo = self.SourceInfoClass(
263  schema = sourceCat.schema,
264  fluxType = self.config.sourceFluxType,
265  )
266 
267  # usableSourceCat: sources that are good but may be saturated
268  numSources = len(sourceCat)
269  usableSourceCat = afwTable.SourceCatalog(sourceCat.table)
270  usableSourceCat.extend(s for s in sourceCat if sourceInfo.isUsable(s))
271  numUsableSources = len(usableSourceCat)
272  self.log.info("Purged %d unusable sources, leaving %d usable sources" % \
273  (numSources - numUsableSources, numUsableSources))
274 
275  if len(usableSourceCat) == 0:
276  raise pipeBase.TaskError("No sources are usable")
277 
278  del sourceCat # avoid accidentally using sourceCat; use usableSourceCat or goodSourceCat from now on
279 
280  minMatchedPairs = min(self.config.minMatchedPairs,
281  int(self.config.minFracMatchedPairs * min([len(refCat), len(usableSourceCat)])))
282 
283  # match usable (possibly saturated) sources and then purge saturated sources from the match list
284  usableMatches = self._doMatch(
285  refCat = refCat,
286  sourceCat = usableSourceCat,
287  wcs = wcs,
288  refFluxField = refFluxField,
289  numUsableSources = numUsableSources,
290  minMatchedPairs = minMatchedPairs,
291  maxMatchDistArcSec = maxMatchDistArcSec,
292  sourceInfo = sourceInfo,
293  verbose = debug.verbose,
294  )
295 
296  # cull non-good sources
297  matches = []
298  for match in usableMatches:
299  if sourceInfo.isGood(match.second):
300  matches.append(match)
301 
302  self.log.logdebug("Found %d usable matches, of which %d had good sources" %
303  (len(usableMatches), len(matches)))
304 
305  if len(matches) == 0:
306  raise RuntimeError("Unable to match sources")
307 
308  self.log.info("Matched %d sources" % len(matches))
309  if len(matches) < minMatchedPairs:
310  self.log.warn("Number of matches is smaller than request")
311 
312  return pipeBase.Struct(
313  matches = matches,
314  usableSourceCat = usableSourceCat,
315  )
316 
317  @pipeBase.timeMethod
318  def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources, minMatchedPairs,
319  maxMatchDistArcSec, sourceInfo, verbose):
320  """!Implementation of matching sources to position reference stars
321 
322  Unlike matchObjectsToSources, this method does not check if the sources are suitable.
323 
324  @param[in] refCat catalog of position reference stars that overlap an exposure
325  @param[in] sourceCat catalog of sources found on the exposure
326  @param[in] wcs estimated WCS of exposure
327  @param[in] refFluxField field of refCat to use for flux
328  @param[in] numUsableSources number of usable sources (sources with known centroid
329  that are not near the edge, but may be saturated)
330  @param[in] minMatchedPairs minimum number of matches
331  @param[in] maxMatchDistArcSec maximum distance between reference objects and sources (arcsec);
332  if specified then min(config.maxMatchDistArcSec, maxMatchDistArcSec) is used
333  if None then config.maxMatchDistArcSec is used
334  @param[in] sourceInfo SourceInfo for the sourceCat
335  @param[in] verbose true to print diagnostic information to std::cout
336 
337  @return a list of matches, an instance of lsst.afw.table.ReferenceMatch
338  """
339  numSources = len(sourceCat)
340  posRefBegInd = numUsableSources - numSources
341  if maxMatchDistArcSec is None:
342  maxMatchDistArcSec = self.config.maxMatchDistArcSec
343  else:
344  maxMatchDistArcSec = min(maxMatchDistArcSec, self.config.maxMatchDistArcSec)
345  configMatchDistPix = maxMatchDistArcSec/wcs.pixelScale().asArcseconds()
346 
347  matchControl = MatchOptimisticBControl()
348  matchControl.refFluxField = refFluxField
349  matchControl.sourceFluxField = sourceInfo.fluxField
350  matchControl.numBrightStars = self.config.numBrightStars
351  matchControl.minMatchedPairs = self.config.minMatchedPairs
352  matchControl.maxOffsetPix = self.config.maxOffsetPix
353  matchControl.numPointsForShape = self.config.numPointsForShape
354  matchControl.maxDeterminant = self.config.maxDeterminant
355 
356  for maxRotInd in range(4):
357  matchControl.maxRotationDeg = self.config.maxRotationDeg * math.pow(2.0, 0.5*maxRotInd)
358  for matchRadInd in range(3):
359  matchControl.matchingAllowancePix = configMatchDistPix * math.pow(1.25, matchRadInd)
360 
361  for angleDiffInd in range(3):
362  matchControl.allowedNonperpDeg = self.config.allowedNonperpDeg*(angleDiffInd+1)
363  matches = matchOptimisticB(
364  refCat,
365  sourceCat,
366  matchControl,
367  wcs,
368  posRefBegInd,
369  verbose,
370  )
371  if matches is not None and len(matches) > 0:
372  setMatchDistance(matches)
373  return matches
374  return matches
def matchObjectsToSources
Match sources to position reference stars.
lsst::afw::table::ReferenceMatchVector matchOptimisticB(lsst::afw::table::SimpleCatalog const &posRefCat, lsst::afw::table::SourceCatalog const &sourceCat, MatchOptimisticBControl const &control, afw::image::Wcs const &wcs, int posRefBegInd=0, bool verbose=false)
def _doMatch
Implementation of matching sources to position reference stars.
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55