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
measurePsf.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010, 2011 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 import lsst.afw.math as afwMath
23 import lsst.afw.display.ds9 as ds9
24 import lsst.meas.algorithms as measAlg
25 import lsst.meas.algorithms.utils as maUtils
26 import lsst.pex.config as pexConfig
27 import lsst.pipe.base as pipeBase
28 import lsst.afw.table as afwTable
29 
30 class MeasurePsfConfig(pexConfig.Config):
31  starSelector = measAlg.starSelectorRegistry.makeField("Star selection algorithm", default="secondMoment")
32  psfDeterminer = measAlg.psfDeterminerRegistry.makeField("PSF Determination algorithm", default="pca")
33 
34 ## \addtogroup LSST_task_documentation
35 ## \{
36 ## \page MeasurePsfTask
37 ## \ref MeasurePsfTask_ "MeasurePsfTask"
38 ## \copybrief MeasurePsfTask
39 ## \}
40 
41 class MeasurePsfTask(pipeBase.Task):
42  """!
43 \anchor MeasurePsfTask_
44 
45 \brief Measure the PSF
46 
47 \section pipe_tasks_measurePsf_Contents Contents
48 
49  - \ref pipe_tasks_measurePsf_Purpose
50  - \ref pipe_tasks_measurePsf_Initialize
51  - \ref pipe_tasks_measurePsf_IO
52  - \ref pipe_tasks_measurePsf_Config
53  - \ref pipe_tasks_measurePsf_Debug
54  - \ref pipe_tasks_measurePsf_Example
55 
56 \section pipe_tasks_measurePsf_Purpose Description
57 
58 A task that wraps two algorithms set via a pair of registries specified in the task's
59 \ref pipe_tasks_measurePsf_Config.
60 Both algorithms are classes with a constructor taking a pex.config.Config object (\em e.g.
61 lsst.meas.algorithms.objectSizeStarSelector.ObjectSizeStarSelector.__init__).
62 
63 The algorithms are:
64  - a star selector with API:
65 \code
66 selectStars(self, exposure, catalog, matches=None)
67 \endcode
68 which returns a list of lsst.meas.algorithms.PsfCandidate (\em e.g.
69 lsst.meas.algorithms.objectSizeStarSelector.ObjectSizeStarSelector.selectStars)
70 
71  - a psf estimator with API:
72 \code
73 determinePsf(exposure, psfCandidateList, metadata=None, flagKey=None)
74 \endcode
75 which returns an lsst.afw.detection.Psf and lsst.afw.math.SpatialCellSet (\em e.g.
76 lsst.meas.algorithms.pcaPsfDeterminer.PcaPsfDeterminer.determinePsf).
77 MeasurePsfTask calls determinePsf with \c flagKey set to
78 "calib.psf.used" if a schema is passed to its constructor (see \ref pipe_tasks_measurePsf_Initialize).
79 
80 See also lsst.meas.algorithms.starSelectorRegistry.starSelectorRegistry and
81 lsst.meas.algorithms.psfDeterminerRegistry.psfDeterminerRegistry.
82 
83 \warning
84 There is no establised set of configuration parameters for these algorithms, so once you start modifying
85 parameters (as we do in \ref pipe_tasks_measurePsf_Example) your code is no longer portable.
86 
87 \section pipe_tasks_measurePsf_Initialize Task initialisation
88 
89 \copydoc \_\_init\_\_
90 
91 \section pipe_tasks_measurePsf_IO Invoking the Task
92 
93 \copydoc run
94 
95 \section pipe_tasks_measurePsf_Config Configuration parameters
96 
97 See \ref MeasurePsfConfig.
98 
99 \warning
100 The star selector and psf determiner registries should be modified to return a class
101 which has a ConfigClass attribute and can be instantiated with a config. Until then, there's no
102 obvious way to get a registry algorithm's Config from another Config.
103 
104 \section pipe_tasks_measurePsf_Debug Debug variables
105 
106 The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a
107 flag \c -d to import \b debug.py from your \c PYTHONPATH; see \ref baseDebug for more about \b debug.py files.
108 
109 <DL>
110  <DT> \c display
111  <DD> If True, display debugging plots
112  <DT> displayExposure
113  <DD> display the Exposure + spatialCells
114  <DT> displayPsfCandidates
115  <DD> show mosaic of candidates
116  <DT> showBadCandidates
117  <DD> Include bad candidates
118  <DT> displayPsfMosaic
119  <DD> show mosaic of reconstructed PSF(xy)
120  <DT> displayResiduals
121  <DD> show residuals
122  <DT> normalizeResiduals
123  <DD> Normalise residuals by object amplitude
124 </DL>
125 
126 Additionally you can enable any debug outputs that your chosen star selector and psf determiner support.
127 
128 \section pipe_tasks_measurePsf_Example A complete example of using MeasurePsfTask
129 
130 This code is in \link measurePsfTask.py\endlink in the examples directory, and can be run as \em e.g.
131 \code
132 examples/measurePsfTask.py --ds9
133 \endcode
134 \dontinclude measurePsfTask.py
135 
136 The example also runs SourceDetectionTask and SourceMeasurementTask; see \ref meas_algorithms_measurement_Example for more explanation.
137 
138 Import the tasks (there are some other standard imports; read the file to see them all):
139 
140 \skip SourceDetectionTask
141 \until MeasurePsfTask
142 
143 We need to create the tasks before processing any data as the task constructor
144 can add an extra column to the schema, but first we need an almost-empty
145 Schema:
146 
147 \skipline makeMinimalSchema
148 
149 We can now call the constructors for the tasks we need to find and characterize candidate
150 PSF stars:
151 
152 \skip SourceDetectionTask.ConfigClass
153 \until measureTask
154 
155 Note that we've chosen a minimal set of measurement plugins: we need the
156 outputs of \c base_SdssCentroid, \c base_SdssShape and \c base_CircularApertureFlux as
157 inputs to the PSF measurement algorithm, while \c base_PixelFlags identifies
158 and flags bad sources (e.g. with pixels too close to the edge) so they can be
159 removed later.
160 
161 Now we can create and configure the task that we're interested in:
162 
163 \skip MeasurePsfTask
164 \until measurePsfTask
165 
166 We're now ready to process the data (we could loop over multiple exposures/catalogues using the same
167 task objects). First create the output table:
168 
169 \skipline afwTable
170 
171 And process the image:
172 
173 \skip sources =
174 \until result
175 
176 We can then unpack and use the results:
177 
178 \skip psf
179 \until cellSet
180 
181 If you specified \c --ds9 you can see the PSF candidates:
182 
183 \skip display
184 \until RED
185 
186 <HR>
187 
188 To investigate the \ref pipe_tasks_measurePsf_Debug, put something like
189 \code{.py}
190  import lsstDebug
191  def DebugInfo(name):
192  di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
193 
194  if name == "lsst.pipe.tasks.measurePsf" :
195  di.display = True
196  di.displayExposure = False # display the Exposure + spatialCells
197  di.displayPsfCandidates = True # show mosaic of candidates
198  di.displayPsfMosaic = True # show mosaic of reconstructed PSF(xy)
199  di.displayResiduals = True # show residuals
200  di.showBadCandidates = True # Include bad candidates
201  di.normalizeResiduals = False # Normalise residuals by object amplitude
202 
203  return di
204 
205  lsstDebug.Info = DebugInfo
206 \endcode
207 into your debug.py file and run measurePsfTask.py with the \c --debug flag.
208  """
209  ConfigClass = MeasurePsfConfig
210  _DefaultName = "measurePsf"
211 
212  def __init__(self, schema=None, **kwargs):
213  """!Create the detection task. Most arguments are simply passed onto pipe.base.Task.
214 
215  \param schema An lsst::afw::table::Schema used to create the output lsst.afw.table.SourceCatalog
216  \param **kwargs Keyword arguments passed to lsst.pipe.base.task.Task.__init__.
217 
218  If schema is not None, 'calib.psf.candidate' and 'calib.psf.used' fields will be added to
219  identify which stars were employed in the PSF estimation.
220 
221  \note This task can add fields to the schema, so any code calling this task must ensure that
222  these fields are indeed present in the input table.
223  """
224 
225  pipeBase.Task.__init__(self, **kwargs)
226  if schema is not None:
227  self.candidateKey = schema.addField(
228  "calib_psfCandidate", type="Flag",
229  doc=("Flag set if the source was a candidate for PSF determination, "
230  "as determined by the '%s' star selector.") % self.config.starSelector.name
231  )
232  self.usedKey = schema.addField(
233  "calib_psfUsed", type="Flag",
234  doc=("Flag set if the source was actually used for PSF determination, "
235  "as determined by the '%s' PSF determiner.") % self.config.psfDeterminer.name
236  )
237  else:
238  self.candidateKey = None
239  self.usedKey = None
240  self.starSelector = self.config.starSelector.apply()
241  self.psfDeterminer = self.config.psfDeterminer.apply()
242 
243  @pipeBase.timeMethod
244  def run(self, exposure, sources, matches=None):
245  """!Measure the PSF
246 
247  \param[in,out] exposure Exposure to process; measured PSF will be added.
248  \param[in,out] sources Measured sources on exposure; flag fields will be set marking
249  stars chosen by the star selector and the PSF determiner if a schema
250  was passed to the task constructor.
251 
252  \param[in] matches a list of lsst.afw.table.ReferenceMatch objects (\em i.e. of lsst.afw.table.Match
253  with \c first being of type lsst.afw.table.SimpleRecord and \c second
254  type lsst.afw.table.SourceRecord --- the reference object and detected
255  object respectively) as returned by \em e.g. the AstrometryTask.
256  Used by star selectors that choose to refer to an external catalog.
257 
258  \return a pipe.base.Struct with fields:
259  - psf: The measured PSF (also set in the input exposure)
260  - cellSet: an lsst.afw.math.SpatialCellSet containing the PSF candidates as returned by the psf determiner.
261  """
262  self.log.info("Measuring PSF")
263 
264  import lsstDebug
265  display = lsstDebug.Info(__name__).display
266  displayExposure = lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
267  displayPsfMosaic = lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)
268  displayPsfCandidates = lsstDebug.Info(__name__).displayPsfCandidates # show mosaic of candidates
269  displayResiduals = lsstDebug.Info(__name__).displayResiduals # show residuals
270  showBadCandidates = lsstDebug.Info(__name__).showBadCandidates # include bad candidates
271  normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals # normalise residuals by object peak
272 
273  #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
274  #
275  # Run star selector
276  #
277  psfCandidateList = self.starSelector.selectStars(exposure, sources, matches=matches)
278  if psfCandidateList and self.candidateKey is not None:
279  for cand in psfCandidateList:
280  source = cand.getSource()
281  source.set(self.candidateKey, True)
282 
283  self.log.info("PSF star selector found %d candidates" % len(psfCandidateList))
284 
285  if display:
286  frame = display
287  if displayExposure:
288  ds9.mtv(exposure, frame=frame, title="psf determination")
289 
290  #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
291  #
292  # Determine PSF
293  #
294  psf, cellSet = self.psfDeterminer.determinePsf(exposure, psfCandidateList, self.metadata,
295  flagKey=self.usedKey)
296  self.log.info("PSF determination using %d/%d stars." %
297  (self.metadata.get("numGoodStars"), self.metadata.get("numAvailStars")))
298 
299  exposure.setPsf(psf)
300 
301  if display:
302  frame = display
303  if displayExposure:
304  showPsfSpatialCells(exposure, cellSet, showBadCandidates, frame=frame)
305  frame += 1
306 
307  if displayPsfCandidates: # Show a mosaic of PSF candidates
308  plotPsfCandidates(cellSet, showBadCandidates, frame)
309  frame += 1
310 
311  if displayResiduals:
312  frame = plotResiduals(exposure, cellSet,
313  showBadCandidates=showBadCandidates,
314  normalizeResiduals=normalizeResiduals,
315  frame=frame)
316  if displayPsfMosaic:
317  maUtils.showPsfMosaic(exposure, psf, frame=frame, showFwhm=True)
318  ds9.scale(0, 1, "linear", frame=frame)
319  frame += 1
320 
321  return pipeBase.Struct(
322  psf = psf,
323  cellSet = cellSet,
324  )
325 
326 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
327 #
328 # Debug code
329 #
330 def showPsfSpatialCells(exposure, cellSet, showBadCandidates, frame=1):
331  maUtils.showPsfSpatialCells(exposure, cellSet,
332  symb="o", ctype=ds9.CYAN, ctypeUnused=ds9.YELLOW,
333  size=4, frame=frame)
334  for cell in cellSet.getCellList():
335  for cand in cell.begin(not showBadCandidates): # maybe include bad candidates
336  cand = measAlg.cast_PsfCandidateF(cand)
337  status = cand.getStatus()
338  ds9.dot('+', *cand.getSource().getCentroid(), frame=frame,
339  ctype=ds9.GREEN if status == afwMath.SpatialCellCandidate.GOOD else
340  ds9.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else ds9.RED)
341 
342 def plotPsfCandidates(cellSet, showBadCandidates=False, frame=1):
343  import lsst.afw.display.utils as displayUtils
344 
345  stamps = []
346  for cell in cellSet.getCellList():
347  for cand in cell.begin(not showBadCandidates): # maybe include bad candidates
348  cand = measAlg.cast_PsfCandidateF(cand)
349 
350  try:
351  im = cand.getMaskedImage()
352 
353  chi2 = cand.getChi2()
354  if chi2 < 1e100:
355  chi2 = "%.1f" % chi2
356  else:
357  chi2 = numpy.nan
358 
359  stamps.append((im, "%d%s" %
360  (maUtils.splitId(cand.getSource().getId(), True)["objId"], chi2),
361  cand.getStatus()))
362  except Exception, e:
363  continue
364 
365  mos = displayUtils.Mosaic()
366  for im, label, status in stamps:
367  im = type(im)(im, True)
368  try:
369  im /= afwMath.makeStatistics(im, afwMath.MAX).getValue()
370  except NotImplementedError:
371  pass
372 
373  mos.append(im, label,
374  ds9.GREEN if status == afwMath.SpatialCellCandidate.GOOD else
375  ds9.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else ds9.RED)
376 
377  if mos.images:
378  mos.makeMosaic(frame=frame, title="Psf Candidates")
379 
380 def plotResiduals(exposure, cellSet, showBadCandidates=False, normalizeResiduals=True, frame=2):
381  psf = exposure.getPsf()
382  while True:
383  try:
384  maUtils.showPsfCandidates(exposure, cellSet, psf=psf, frame=frame,
385  normalize=normalizeResiduals,
386  showBadCandidates=showBadCandidates)
387  frame += 1
388  maUtils.showPsfCandidates(exposure, cellSet, psf=psf, frame=frame,
389  normalize=normalizeResiduals,
390  showBadCandidates=showBadCandidates,
391  variance=True)
392  frame += 1
393  except Exception as e:
394  if not showBadCandidates:
395  showBadCandidates = True
396  continue
397  break
398 
399  return frame
def __init__
Create the detection task.
Definition: measurePsf.py:212
Statistics makeStatistics(afwImage::Mask< afwImage::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl)
Specialization to handle Masks.
Definition: Statistics.cc:1082