LSSTApplications  16.0-10-g0ee56ad+5,16.0-11-ga33d1f2+5,16.0-12-g3ef5c14+3,16.0-12-g71e5ef5+18,16.0-12-gbdf3636+3,16.0-13-g118c103+3,16.0-13-g8f68b0a+3,16.0-15-gbf5c1cb+4,16.0-16-gfd17674+3,16.0-17-g7c01f5c+3,16.0-18-g0a50484+1,16.0-20-ga20f992+8,16.0-21-g0e05fd4+6,16.0-21-g15e2d33+4,16.0-22-g62d8060+4,16.0-22-g847a80f+4,16.0-25-gf00d9b8+1,16.0-28-g3990c221+4,16.0-3-gf928089+3,16.0-32-g88a4f23+5,16.0-34-gd7987ad+3,16.0-37-gc7333cb+2,16.0-4-g10fc685+2,16.0-4-g18f3627+26,16.0-4-g5f3a788+26,16.0-5-gaf5c3d7+4,16.0-5-gcc1f4bb+1,16.0-6-g3b92700+4,16.0-6-g4412fcd+3,16.0-6-g7235603+4,16.0-69-g2562ce1b+2,16.0-8-g14ebd58+4,16.0-8-g2df868b+1,16.0-8-g4cec79c+6,16.0-8-gadf6c7a+1,16.0-8-gfc7ad86,16.0-82-g59ec2a54a+1,16.0-9-g5400cdc+2,16.0-9-ge6233d7+5,master-g2880f2d8cf+3,v17.0.rc1
LSSTDataManagementBasePackage
imagePsfMatch.py
Go to the documentation of this file.
1 # LSST Data Management System
2 # Copyright 2008-2016 LSST Corporation.
3 #
4 # This product includes software developed by the
5 # LSST Project (http://www.lsst.org/).
6 #
7 # This program is free software: you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation, either version 3 of the License, or
10 # (at your option) any later version.
11 #
12 # This program is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
16 #
17 # You should have received a copy of the LSST License Statement and
18 # the GNU General Public License along with this program. If not,
19 # see <http://www.lsstcorp.org/LegalNotices/>.
20 #
21 
22 __all__ = ["ImagePsfMatchConfig", "ImagePsfMatchTask", "subtractAlgorithmRegistry"]
23 
24 import numpy as np
25 
26 import lsst.daf.base as dafBase
27 import lsst.pex.config as pexConfig
28 import lsst.afw.image as afwImage
29 import lsst.afw.math as afwMath
30 import lsst.afw.geom as afwGeom
31 import lsst.afw.table as afwTable
32 import lsst.pipe.base as pipeBase
33 from lsst.meas.algorithms import SourceDetectionTask, SubtractBackgroundTask
34 from lsst.meas.base import SingleFrameMeasurementTask
35 from .makeKernelBasisList import makeKernelBasisList
36 from .psfMatch import PsfMatchTask, PsfMatchConfigDF, PsfMatchConfigAL
37 from . import utils as dituils
38 from . import diffimLib
39 from . import diffimTools
40 import lsst.afw.display.ds9 as ds9
41 
42 sigma2fwhm = 2. * np.sqrt(2. * np.log(2.))
43 
44 
45 class ImagePsfMatchConfig(pexConfig.Config):
46  """Configuration for image-to-image Psf matching"""
47  kernel = pexConfig.ConfigChoiceField(
48  doc="kernel type",
49  typemap=dict(
50  AL=PsfMatchConfigAL,
51  DF=PsfMatchConfigDF
52  ),
53  default="AL",
54  )
55  selectDetection = pexConfig.ConfigurableField(
56  target=SourceDetectionTask,
57  doc="Initial detections used to feed stars to kernel fitting",
58  )
59  selectMeasurement = pexConfig.ConfigurableField(
60  target=SingleFrameMeasurementTask,
61  doc="Initial measurements used to feed stars to kernel fitting",
62  )
63 
64  def setDefaults(self):
65  # High sigma detections only
66  self.selectDetection.reEstimateBackground = False
67  self.selectDetection.thresholdValue = 10.0
68 
69  # Minimal set of measurments for star selection
70  self.selectMeasurement.algorithms.names.clear()
71  self.selectMeasurement.algorithms.names = ('base_SdssCentroid', 'base_PsfFlux', 'base_PixelFlags',
72  'base_SdssShape', 'base_GaussianFlux', 'base_SkyCoord')
73  self.selectMeasurement.slots.modelFlux = None
74  self.selectMeasurement.slots.apFlux = None
75  self.selectMeasurement.slots.calibFlux = None
76 
77 
79  """Psf-match two MaskedImages or Exposures using the sources in the images
80 
81  Parameters
82  ----------
83  args :
84  arguments to be passed to lsst.ip.diffim.PsfMatchTask.__init__
85  kwargs :
86  keyword arguments to be passed to lsst.ip.diffim.PsfMatchTask.__init__
87 
88  Notes
89  -----
90  Upon initialization, the kernel configuration is defined by self.config.kernel.active.
91  The task creates an lsst.afw.math.Warper from the subConfig self.config.kernel.active.warpingConfig.
92  A schema for the selection and measurement of candidate lsst.ip.diffim.KernelCandidates is
93  defined, and used to initize subTasks selectDetection (for candidate detection) and selectMeasurement
94  (for candidate measurement).
95 
96  Description
97 
98  Build a Psf-matching kernel using two input images, either as MaskedImages (in which case they need
99  to be astrometrically aligned) or Exposures (in which case astrometric alignment will happen by
100  default but may be turned off). This requires a list of input Sources which may be provided
101  by the calling Task; if not, the Task will perform a coarse source detection
102  and selection for this purpose. Sources are vetted for signal-to-noise and masked pixels
103  (in both the template and science image), and substamps around each acceptable
104  source are extracted and used to create an instance of KernelCandidate.
105  Each KernelCandidate is then placed within a lsst.afw.math.SpatialCellSet, which is used by an ensemble of
106  lsst.afw.math.CandidateVisitor instances to build the Psf-matching kernel. These visitors include, in
107  the order that they are called: BuildSingleKernelVisitor, KernelSumVisitor, BuildSpatialKernelVisitor,
108  and AssessSpatialKernelVisitor.
109 
110  Sigma clipping of KernelCandidates is performed as follows:
111 
112  - BuildSingleKernelVisitor, using the substamp diffim residuals from the per-source kernel fit,
113  if PsfMatchConfig.singleKernelClipping is True
114  - KernelSumVisitor, using the mean and standard deviation of the kernel sum from all candidates,
115  if PsfMatchConfig.kernelSumClipping is True
116  - AssessSpatialKernelVisitor, using the substamp diffim ressiduals from the spatial kernel fit,
117  if PsfMatchConfig.spatialKernelClipping is True
118 
119  The actual solving for the kernel (and differential background model) happens in
120  lsst.ip.diffim.PsfMatchTask._solve. This involves a loop over the SpatialCellSet that first builds the
121  per-candidate matching kernel for the requested number of KernelCandidates per cell
122  (PsfMatchConfig.nStarPerCell). The quality of this initial per-candidate difference image is examined,
123  using moments of the pixel residuals in the difference image normalized by the square root of the variance
124  (i.e. sigma); ideally this should follow a normal (0, 1) distribution,
125  but the rejection thresholds are set
126  by the config (PsfMatchConfig.candidateResidualMeanMax and PsfMatchConfig.candidateResidualStdMax).
127  All candidates that pass this initial build are then examined en masse to find the
128  mean/stdev of the kernel sums across all candidates.
129  Objects that are significantly above or below the mean,
130  typically due to variability or sources that are saturated in one image but not the other,
131  are also rejected.This threshold is defined by PsfMatchConfig.maxKsumSigma.
132  Finally, a spatial model is built using all currently-acceptable candidates,
133  and the spatial model used to derive a second set of (spatial) residuals
134  which are again used to reject bad candidates, using the same thresholds as above.
135 
136  Invoking the Task
137 
138  There is no run() method for this Task. Instead there are 4 methods that
139  may be used to invoke the Psf-matching. These are
140  lsst.ip.diffim.imagePsfMatch.ImagePsfMatchTask.matchMaskedImages matchMaskedImages
141  lsst.ip.diffim.imagePsfMatch.ImagePsfMatchTask.subtractMaskedImages subtractMaskedImages,
142  lsst.ip.diffim.imagePsfMatch.ImagePsfMatchTask.matchExposures matchExposures, and
143  lsst.ip.diffim.imagePsfMatch.ImagePsfMatchTask.subtractExposures subtractExposures.
144 
145  The methods that operate on lsst.afw.image.MaskedImage require that the images already be astrometrically
146  aligned, and are the same shape. The methods that operate on lsst.afw.image.Exposure allow for the
147  input images to be misregistered and potentially be different sizes; by default a
148  lsst.afw.math.LanczosWarpingKernel is used to perform the astrometric alignment. The methods
149  that "match" images return a Psf-matched image, while the methods that "subtract" images
150  return a Psf-matched and template subtracted image.
151 
152  See each method's returned lsst.pipe.base.Struct for more details.
153 
154  Debug variables
155 
156  The lsst.pipe.base.cmdLineTask.CmdLineTask command line task interface supports a
157  flag -d/--debug to import debug.py from your PYTHONPATH. The relevant contents of debug.py
158  for this Task include:
159 
160  .. code-block:: py
161 
162  import sys
163  import lsstDebug
164  def DebugInfo(name):
165  di = lsstDebug.getInfo(name)
166  if name == "lsst.ip.diffim.psfMatch":
167  di.display = True # enable debug output
168  di.maskTransparency = 80 # ds9 mask transparency
169  di.displayCandidates = True # show all the candidates and residuals
170  di.displayKernelBasis = False # show kernel basis functions
171  di.displayKernelMosaic = True # show kernel realized across the image
172  di.plotKernelSpatialModel = False # show coefficients of spatial model
173  di.showBadCandidates = True # show the bad candidates (red) along with good (green)
174  elif name == "lsst.ip.diffim.imagePsfMatch":
175  di.display = True # enable debug output
176  di.maskTransparency = 30 # ds9 mask transparency
177  di.displayTemplate = True # show full (remapped) template
178  di.displaySciIm = True # show science image to match to
179  di.displaySpatialCells = True # show spatial cells
180  di.displayDiffIm = True # show difference image
181  di.showBadCandidates = True # show the bad candidates (red) along with good (green)
182  elif name == "lsst.ip.diffim.diaCatalogSourceSelector":
183  di.display = False # enable debug output
184  di.maskTransparency = 30 # ds9 mask transparency
185  di.displayExposure = True # show exposure with candidates indicated
186  di.pauseAtEnd = False # pause when done
187  return di
188  lsstDebug.Info = DebugInfo
189  lsstDebug.frame = 1
190 
191  Note that if you want addional logging info, you may add to your scripts:
192 
193  .. code-block:: py
194 
195  import lsst.log.utils as logUtils
196  logUtils.traceSetAt("ip.diffim", 4)
197 
198  Examples
199  --------
200  A complete example of using ImagePsfMatchTask
201 
202  This code is imagePsfMatchTask.py in the examples directory, and can be run as e.g.
203 
204  ..code-block:: none
205 
206  examples/imagePsfMatchTask.py --debug
207  examples/imagePsfMatchTask.py --debug --mode="matchExposures"
208  examples/imagePsfMatchTask.py --debug --template /path/to/templateExp.fits
209  --science /path/to/scienceExp.fits
210 
211  Create a subclass of ImagePsfMatchTask that allows us to either match exposures, or subtract exposures:
212 
213  ..code-block:: none
214 
215  class MyImagePsfMatchTask(ImagePsfMatchTask):
216 
217  def __init__(self, args, kwargs):
218  ImagePsfMatchTask.__init__(self, args, kwargs)
219 
220  def run(self, templateExp, scienceExp, mode):
221  if mode == "matchExposures":
222  return self.matchExposures(templateExp, scienceExp)
223  elif mode == "subtractExposures":
224  return self.subtractExposures(templateExp, scienceExp)
225 
226  And allow the user the freedom to either run the script in default mode,
227  or point to their own images on disk.
228  Note that these images must be readable as an lsst.afw.image.Exposure.
229 
230  We have enabled some minor display debugging in this script via the --debug option. However, if you
231  have an lsstDebug debug.py in your PYTHONPATH you will get additional debugging displays. The following
232  block checks for this script:
233 
234  .. code-block:: py
235 
236  if args.debug:
237  try:
238  import debug
239  # Since I am displaying 2 images here, set the starting frame number for the LSST debug LSST
240  debug.lsstDebug.frame = 3
241  except ImportError as e:
242  print(e, file=sys.stderr)
243 
244  Finally, we call a run method that we define below.
245  First set up a Config and modify some of the parameters.
246  E.g. use an "Alard-Lupton" sum-of-Gaussian basis,
247  fit for a differential background, and use low order spatial
248  variation in the kernel and background:
249 
250  .. code-block:: py
251 
252  def run(args):
253  #
254  # Create the Config and use sum of gaussian basis
255  #
256  config = ImagePsfMatchTask.ConfigClass()
257  config.kernel.name = "AL"
258  config.kernel.active.fitForBackground = True
259  config.kernel.active.spatialKernelOrder = 1
260  config.kernel.active.spatialBgOrder = 0
261 
262  Make sure the images (if any) that were sent to the script exist on disk and are readable. If no images
263  are sent, make some fake data up for the sake of this example script (have a look at the code if you want
264  more details on generateFakeImages):
265 
266  .. code-block:: py
267 
268  # Run the requested method of the Task
269  if args.template is not None and args.science is not None:
270  if not os.path.isfile(args.template):
271  raise Exception("Template image %s does not exist" % (args.template))
272  if not os.path.isfile(args.science):
273  raise Exception("Science image %s does not exist" % (args.science))
274  try:
275  templateExp = afwImage.ExposureF(args.template)
276  except Exception as e:
277  raise Exception("Cannot read template image %s" % (args.template))
278  try:
279  scienceExp = afwImage.ExposureF(args.science)
280  except Exception as e:
281  raise Exception("Cannot read science image %s" % (args.science))
282  else:
283  templateExp, scienceExp = generateFakeImages()
284  config.kernel.active.sizeCellX = 128
285  config.kernel.active.sizeCellY = 128
286 
287  Create and run the Task:
288 
289  .. code-block:: py
290 
291  # Create the Task
292  psfMatchTask = MyImagePsfMatchTask(config=config)
293  # Run the Task
294  result = psfMatchTask.run(templateExp, scienceExp, args.mode)
295 
296  And finally provide some optional debugging displays:
297 
298  .. code-block:: py
299  if args.debug:
300  # See if the LSST debug has incremented the frame number; if not start with frame 3
301  try:
302  frame = debug.lsstDebug.frame + 1
303  except Exception:
304  frame = 3
305  ds9.mtv(result.matchedExposure, frame=frame, title="Example script: Matched Template Image")
306  if "subtractedExposure" in result.getDict():
307  ds9.mtv(result.subtractedExposure, frame=frame+1, title="Example script: Subtracted Image")
308 
309  """
310 
311  ConfigClass = ImagePsfMatchConfig
312 
313  def __init__(self, *args, **kwargs):
314  """Create the ImagePsfMatchTask
315 
316  """
317  PsfMatchTask.__init__(self, *args, **kwargs)
318  self.kConfig = self.config.kernel.active
319  self._warper = afwMath.Warper.fromConfig(self.kConfig.warpingConfig)
320  # the background subtraction task uses a config from an unusual location,
321  # so cannot easily be constructed with makeSubtask
322  self.background = SubtractBackgroundTask(config=self.kConfig.afwBackgroundConfig, name="background",
323  parentTask=self)
324  self.selectSchema = afwTable.SourceTable.makeMinimalSchema()
326  self.makeSubtask("selectDetection", schema=self.selectSchema)
327  self.makeSubtask("selectMeasurement", schema=self.selectSchema, algMetadata=self.selectAlgMetadata)
328 
329  def getFwhmPix(self, psf):
330  """Return the FWHM in pixels of a Psf"""
331  sigPix = psf.computeShape().getDeterminantRadius()
332  return sigPix * sigma2fwhm
333 
334  @pipeBase.timeMethod
335  def matchExposures(self, templateExposure, scienceExposure,
336  templateFwhmPix=None, scienceFwhmPix=None,
337  candidateList=None, doWarping=True, convolveTemplate=True):
338  """Warp and PSF-match an exposure to the reference
339 
340  Do the following, in order:
341 
342  - Warp templateExposure to match scienceExposure,
343  if doWarping True and their WCSs do not already match
344  - Determine a PSF matching kernel and differential background model
345  that matches templateExposure to scienceExposure
346  - Convolve templateExposure by PSF matching kernel
347 
348  Parameters
349  ----------
350  templateExposure : `lsst.afw.image.Exposure`
351  Exposure to warp and PSF-match to the reference masked image
352  scienceExposure : `lsst.afw.image.Exposure`
353  Exposure whose WCS and PSF are to be matched to
354  templateFwhmPix :`float`
355  FWHM (in pixels) of the Psf in the template image (image to convolve)
356  scienceFwhmPix : `float`
357  FWHM (in pixels) of the Psf in the science image
358  candidateList : `list`, optional
359  a list of footprints/maskedImages for kernel candidates;
360  if None then source detection is run.
361 
362  - Currently supported: list of Footprints or measAlg.PsfCandidateF
363 
364  doWarping : `bool`
365  what to do if templateExposure's and scienceExposure's WCSs do not match:
366 
367  - if True then warp templateExposure to match scienceExposure
368  - if False then raise an Exception
369 
370  convolveTemplate : `bool`
371  convolve the template image or the science image:
372 
373  - if True, templateExposure is warped if doWarping, templateExposure is convolved
374  - if False, templateExposure is warped if doWarping, scienceExposure is convolved
375 
376  Returns
377  -------
378  results : `Struct`
379  a `lsst.pipe.base.Struct` containing these fields:
380 
381  - ``matchedImage`` : the PSF-matched exposure =
382  warped templateExposure convolved by psfMatchingKernel. This has:
383 
384  - the same parent bbox, Wcs and Calib as scienceExposure
385  - the same filter as templateExposure
386  - no Psf (because the PSF-matching process does not compute one)
387 
388  - ``psfMatchingKernel`` : the PSF matching kernel
389  - ``backgroundModel`` : differential background model
390  - ``kernelCellSet`` : SpatialCellSet used to solve for the PSF matching kernel
391 
392  Raises
393  ------
394  RuntimeError
395  if doWarping is False and templateExposure's and scienceExposure's
396  WCSs do not match
397  """
398  if not self._validateWcs(templateExposure, scienceExposure):
399  if doWarping:
400  self.log.info("Astrometrically registering template to science image")
401  templatePsf = templateExposure.getPsf()
402  templateExposure = self._warper.warpExposure(scienceExposure.getWcs(),
403  templateExposure,
404  destBBox=scienceExposure.getBBox())
405  templateExposure.setPsf(templatePsf)
406  else:
407  self.log.error("ERROR: Input images not registered")
408  raise RuntimeError("Input images not registered")
409 
410  if templateFwhmPix is None:
411  if not templateExposure.hasPsf():
412  self.log.warn("No estimate of Psf FWHM for template image")
413  else:
414  templateFwhmPix = self.getFwhmPix(templateExposure.getPsf())
415  self.log.info("templateFwhmPix: {}".format(templateFwhmPix))
416 
417  if scienceFwhmPix is None:
418  if not scienceExposure.hasPsf():
419  self.log.warn("No estimate of Psf FWHM for science image")
420  else:
421  scienceFwhmPix = self.getFwhmPix(scienceExposure.getPsf())
422  self.log.info("scienceFwhmPix: {}".format(scienceFwhmPix))
423 
424  kernelSize = makeKernelBasisList(self.kConfig, templateFwhmPix, scienceFwhmPix)[0].getWidth()
425  candidateList = self.makeCandidateList(templateExposure, scienceExposure, kernelSize, candidateList)
426 
427  if convolveTemplate:
428  results = self.matchMaskedImages(
429  templateExposure.getMaskedImage(), scienceExposure.getMaskedImage(), candidateList,
430  templateFwhmPix=templateFwhmPix, scienceFwhmPix=scienceFwhmPix)
431  else:
432  results = self.matchMaskedImages(
433  scienceExposure.getMaskedImage(), templateExposure.getMaskedImage(), candidateList,
434  templateFwhmPix=scienceFwhmPix, scienceFwhmPix=templateFwhmPix)
435 
436  psfMatchedExposure = afwImage.makeExposure(results.matchedImage, scienceExposure.getWcs())
437  psfMatchedExposure.setFilter(templateExposure.getFilter())
438  psfMatchedExposure.setCalib(scienceExposure.getCalib())
439  results.warpedExposure = templateExposure
440  results.matchedExposure = psfMatchedExposure
441  return results
442 
443  @pipeBase.timeMethod
444  def matchMaskedImages(self, templateMaskedImage, scienceMaskedImage, candidateList,
445  templateFwhmPix=None, scienceFwhmPix=None):
446  """PSF-match a MaskedImage (templateMaskedImage) to a reference MaskedImage (scienceMaskedImage)
447 
448  Do the following, in order:
449 
450  - Determine a PSF matching kernel and differential background model
451  that matches templateMaskedImage to scienceMaskedImage
452  - Convolve templateMaskedImage by the PSF matching kernel
453 
454  Parameters
455  ----------
456  templateMaskedImage : `lsst.afw.image.MaskedImage`
457  masked image to PSF-match to the reference masked image;
458  must be warped to match the reference masked image
459  scienceMaskedImage : `lsst.afw.image.MaskedImage`
460  maskedImage whose PSF is to be matched to
461  templateFwhmPix : `float`
462  FWHM (in pixels) of the Psf in the template image (image to convolve)
463  scienceFwhmPix : `float`
464  FWHM (in pixels) of the Psf in the science image
465  candidateList : `list`, optional
466  a list of footprints/maskedImages for kernel candidates;
467  if None then source detection is run.
468 
469  - Currently supported: list of Footprints or measAlg.PsfCandidateF
470 
471  Returns
472  -------
473  result : `callable`
474  a `lsst.pipe.base.Struct` containing these fields:
475 
476  - psfMatchedMaskedImage: the PSF-matched masked image =
477  templateMaskedImage convolved with psfMatchingKernel.
478  This has the same xy0, dimensions and wcs as scienceMaskedImage.
479  - psfMatchingKernel: the PSF matching kernel
480  - backgroundModel: differential background model
481  - kernelCellSet: SpatialCellSet used to solve for the PSF matching kernel
482 
483  Raises
484  ------
485  RuntimeError
486  if input images have different dimensions
487  """
488 
489  import lsstDebug
490  display = lsstDebug.Info(__name__).display
491  displayTemplate = lsstDebug.Info(__name__).displayTemplate
492  displaySciIm = lsstDebug.Info(__name__).displaySciIm
493  displaySpatialCells = lsstDebug.Info(__name__).displaySpatialCells
494  maskTransparency = lsstDebug.Info(__name__).maskTransparency
495  if not maskTransparency:
496  maskTransparency = 0
497  if display:
498  ds9.setMaskTransparency(maskTransparency)
499 
500  if not candidateList:
501  raise RuntimeError("Candidate list must be populated by makeCandidateList")
502 
503  if not self._validateSize(templateMaskedImage, scienceMaskedImage):
504  self.log.error("ERROR: Input images different size")
505  raise RuntimeError("Input images different size")
506 
507  if display and displayTemplate:
508  ds9.mtv(templateMaskedImage, frame=lsstDebug.frame, title="Image to convolve")
509  lsstDebug.frame += 1
510 
511  if display and displaySciIm:
512  ds9.mtv(scienceMaskedImage, frame=lsstDebug.frame, title="Image to not convolve")
513  lsstDebug.frame += 1
514 
515  kernelCellSet = self._buildCellSet(templateMaskedImage,
516  scienceMaskedImage,
517  candidateList)
518 
519  if display and displaySpatialCells:
520  dituils.showKernelSpatialCells(scienceMaskedImage, kernelCellSet,
521  symb="o", ctype=ds9.CYAN, ctypeUnused=ds9.YELLOW, ctypeBad=ds9.RED,
522  size=4, frame=lsstDebug.frame, title="Image to not convolve")
523  lsstDebug.frame += 1
524 
525  if templateFwhmPix and scienceFwhmPix:
526  self.log.info("Matching Psf FWHM %.2f -> %.2f pix", templateFwhmPix, scienceFwhmPix)
527 
528  if self.kConfig.useBicForKernelBasis:
529  tmpKernelCellSet = self._buildCellSet(templateMaskedImage,
530  scienceMaskedImage,
531  candidateList)
532  nbe = diffimTools.NbasisEvaluator(self.kConfig, templateFwhmPix, scienceFwhmPix)
533  bicDegrees = nbe(tmpKernelCellSet, self.log)
534  basisList = makeKernelBasisList(self.kConfig, templateFwhmPix, scienceFwhmPix,
535  alardDegGauss=bicDegrees[0], metadata=self.metadata)
536  del tmpKernelCellSet
537  else:
538  basisList = makeKernelBasisList(self.kConfig, templateFwhmPix, scienceFwhmPix,
539  metadata=self.metadata)
540 
541  spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)
542 
543  psfMatchedMaskedImage = afwImage.MaskedImageF(templateMaskedImage.getBBox())
544  doNormalize = False
545  afwMath.convolve(psfMatchedMaskedImage, templateMaskedImage, psfMatchingKernel, doNormalize)
546  return pipeBase.Struct(
547  matchedImage=psfMatchedMaskedImage,
548  psfMatchingKernel=psfMatchingKernel,
549  backgroundModel=backgroundModel,
550  kernelCellSet=kernelCellSet,
551  )
552 
553  @pipeBase.timeMethod
554  def subtractExposures(self, templateExposure, scienceExposure,
555  templateFwhmPix=None, scienceFwhmPix=None,
556  candidateList=None, doWarping=True, convolveTemplate=True):
557  """Register, Psf-match and subtract two Exposures
558 
559  Do the following, in order:
560 
561  - Warp templateExposure to match scienceExposure, if their WCSs do not already match
562  - Determine a PSF matching kernel and differential background model
563  that matches templateExposure to scienceExposure
564  - PSF-match templateExposure to scienceExposure
565  - Compute subtracted exposure (see return values for equation).
566 
567  Parameters
568  ----------
569  templateExposure : `lsst.afw.image.Exposure`
570  exposure to PSF-match to scienceExposure
571  scienceExposure : `lsst.afw.image.Exposure`
572  reference Exposure
573  templateFwhmPix : `float`
574  FWHM (in pixels) of the Psf in the template image (image to convolve)
575  scienceFwhmPix : `float`
576  FWHM (in pixels) of the Psf in the science image
577  candidateList : `list`, optional
578  a list of footprints/maskedImages for kernel candidates;
579  if None then source detection is run.
580 
581  - Currently supported: list of Footprints or measAlg.PsfCandidateF
582 
583  doWarping : `bool`
584  what to do if templateExposure's and scienceExposure's WCSs do not match:
585 
586  - if True then warp templateExposure to match scienceExposure
587  - if False then raise an Exception
588 
589  convolveTemplate : `bool`
590  convolve the template image or the science image
591 
592  - if True, templateExposure is warped if doWarping, templateExposure is convolved
593  - if False, templateExposure is warped if doWarping, scienceExposure is convolved
594 
595  Returns
596  -------
597  result: `Struct`
598  a `lsst.pipe.base.Struct` containing these fields:
599 
600  - ``subtractedExposure`` : subtracted Exposure
601  scienceExposure - (matchedImage + backgroundModel)
602  - ``matchedImage`` : templateExposure after warping to match templateExposure (if doWarping true),
603  and convolving with psfMatchingKernel
604  - ``psfMatchingKernel`` : PSF matching kernel
605  - ``backgroundModel`` : differential background model
606  - ``kernelCellSet`` : SpatialCellSet used to determine PSF matching kernel
607  """
608 
609  results = self.matchExposures(
610  templateExposure=templateExposure,
611  scienceExposure=scienceExposure,
612  templateFwhmPix=templateFwhmPix,
613  scienceFwhmPix=scienceFwhmPix,
614  candidateList=candidateList,
615  doWarping=doWarping,
616  convolveTemplate=convolveTemplate
617  )
618 
619  subtractedExposure = afwImage.ExposureF(scienceExposure, True)
620  if convolveTemplate:
621  subtractedMaskedImage = subtractedExposure.getMaskedImage()
622  subtractedMaskedImage -= results.matchedExposure.getMaskedImage()
623  subtractedMaskedImage -= results.backgroundModel
624  else:
625  subtractedExposure.setMaskedImage(results.warpedExposure.getMaskedImage())
626  subtractedMaskedImage = subtractedExposure.getMaskedImage()
627  subtractedMaskedImage -= results.matchedExposure.getMaskedImage()
628  subtractedMaskedImage -= results.backgroundModel
629 
630  # Preserve polarity of differences
631  subtractedMaskedImage *= -1
632 
633  # Place back on native photometric scale
634  subtractedMaskedImage /= results.psfMatchingKernel.computeImage(
635  afwImage.ImageD(results.psfMatchingKernel.getDimensions()), False)
636 
637  import lsstDebug
638  display = lsstDebug.Info(__name__).display
639  displayDiffIm = lsstDebug.Info(__name__).displayDiffIm
640  maskTransparency = lsstDebug.Info(__name__).maskTransparency
641  if not maskTransparency:
642  maskTransparency = 0
643  if display:
644  ds9.setMaskTransparency(maskTransparency)
645  if display and displayDiffIm:
646  ds9.mtv(templateExposure, frame=lsstDebug.frame, title="Template")
647  lsstDebug.frame += 1
648  ds9.mtv(results.matchedExposure, frame=lsstDebug.frame, title="Matched template")
649  lsstDebug.frame += 1
650  ds9.mtv(scienceExposure, frame=lsstDebug.frame, title="Science Image")
651  lsstDebug.frame += 1
652  ds9.mtv(subtractedExposure, frame=lsstDebug.frame, title="Difference Image")
653  lsstDebug.frame += 1
654 
655  results.subtractedExposure = subtractedExposure
656  return results
657 
658  @pipeBase.timeMethod
659  def subtractMaskedImages(self, templateMaskedImage, scienceMaskedImage, candidateList,
660  templateFwhmPix=None, scienceFwhmPix=None):
661  """Psf-match and subtract two MaskedImages
662 
663  Do the following, in order:
664 
665  - PSF-match templateMaskedImage to scienceMaskedImage
666  - Determine the differential background
667  - Return the difference: scienceMaskedImage
668  ((warped templateMaskedImage convolved with psfMatchingKernel) + backgroundModel)
669 
670  Parameters
671  ----------
672  templateMaskedImage : `lsst.afw.image.MaskedImage`
673  MaskedImage to PSF-match to scienceMaskedImage
674  scienceMaskedImage : `lsst.afw.image.MaskedImage`
675  reference MaskedImage
676  templateFwhmPix : `float`
677  FWHM (in pixels) of the Psf in the template image (image to convolve)
678  scienceFwhmPix : `float`
679  FWHM (in pixels) of the Psf in the science image
680  candidateList : `list`, optional
681  a list of footprints/maskedImages for kernel candidates;
682  if None then source detection is run.
683 
684  - Currently supported: list of Footprints or measAlg.PsfCandidateF
685 
686  Returns
687  -------
688  results : `Struct`
689  a `lsst.pipe.base.Struct` containing these fields:
690 
691  - ``subtractedMaskedImage`` : scienceMaskedImage - (matchedImage + backgroundModel)
692  - ``matchedImage`` : templateMaskedImage convolved with psfMatchingKernel
693  - `psfMatchingKernel`` : PSF matching kernel
694  - ``backgroundModel`` : differential background model
695  - ``kernelCellSet`` : SpatialCellSet used to determine PSF matching kernel
696 
697  """
698  if not candidateList:
699  raise RuntimeError("Candidate list must be populated by makeCandidateList")
700 
701  results = self.matchMaskedImages(
702  templateMaskedImage=templateMaskedImage,
703  scienceMaskedImage=scienceMaskedImage,
704  candidateList=candidateList,
705  templateFwhmPix=templateFwhmPix,
706  scienceFwhmPix=scienceFwhmPix,
707  )
708 
709  subtractedMaskedImage = afwImage.MaskedImageF(scienceMaskedImage, True)
710  subtractedMaskedImage -= results.matchedImage
711  subtractedMaskedImage -= results.backgroundModel
712  results.subtractedMaskedImage = subtractedMaskedImage
713 
714  import lsstDebug
715  display = lsstDebug.Info(__name__).display
716  displayDiffIm = lsstDebug.Info(__name__).displayDiffIm
717  maskTransparency = lsstDebug.Info(__name__).maskTransparency
718  if not maskTransparency:
719  maskTransparency = 0
720  if display:
721  ds9.setMaskTransparency(maskTransparency)
722  if display and displayDiffIm:
723  ds9.mtv(subtractedMaskedImage, frame=lsstDebug.frame)
724  lsstDebug.frame += 1
725 
726  return results
727 
728  def getSelectSources(self, exposure, sigma=None, doSmooth=True, idFactory=None):
729  """Get sources to use for Psf-matching
730 
731  This method runs detection and measurement on an exposure.
732  The returned set of sources will be used as candidates for
733  Psf-matching.
734 
735  Parameters
736  ----------
737  exposure : `lsst.afw.image.Exposure`
738  Exposure on which to run detection/measurement
739  sigma : `float`
740  Detection threshold
741  doSmooth : `bool`
742  Whether or not to smooth the Exposure with Psf before detection
743  idFactory :
744  Factory for the generation of Source ids
745 
746  Returns
747  -------
748  selectSources :
749  source catalog containing candidates for the Psf-matching
750  """
751 
752  if idFactory:
753  table = afwTable.SourceTable.make(self.selectSchema, idFactory)
754  else:
755  table = afwTable.SourceTable.make(self.selectSchema)
756  mi = exposure.getMaskedImage()
757 
758  imArr = mi.getImage().getArray()
759  maskArr = mi.getMask().getArray()
760  miArr = np.ma.masked_array(imArr, mask=maskArr)
761  try:
762  bkgd = self.background.fitBackground(mi).getImageF()
763  except Exception:
764  self.log.warn("Failed to get background model. Falling back to median background estimation")
765  bkgd = np.ma.extras.median(miArr)
766 
767  # Take off background for detection
768  mi -= bkgd
769  try:
770  table.setMetadata(self.selectAlgMetadata)
771  detRet = self.selectDetection.makeSourceCatalog(
772  table=table,
773  exposure=exposure,
774  sigma=sigma,
775  doSmooth=doSmooth
776  )
777  selectSources = detRet.sources
778  self.selectMeasurement.run(measCat=selectSources, exposure=exposure)
779  finally:
780  # Put back on the background in case it is needed down stream
781  mi += bkgd
782  del bkgd
783  return selectSources
784 
785  def makeCandidateList(self, templateExposure, scienceExposure, kernelSize, candidateList=None):
786  """Make a list of acceptable KernelCandidates
787 
788  Accept or generate a list of candidate sources for
789  Psf-matching, and examine the Mask planes in both of the
790  images for indications of bad pixels
791 
792  Parameters
793  ----------
794  templateExposure : `lsst.afw.image.Exposure`
795  Exposure that will be convolved
796  scienceExposure : `lsst.afw.image.Exposure`
797  Exposure that will be matched-to
798  kernelSize : `float`
799  Dimensions of the Psf-matching Kernel, used to grow detection footprints
800  candidateList : `list`, optional
801  List of Sources to examine. Elements must be of type afw.table.Source
802  or a type that wraps a Source and has a getSource() method, such as
803  meas.algorithms.PsfCandidateF.
804 
805  Returns
806  -------
807  candidateList : `list` of `dict`
808  a list of dicts having a "source" and "footprint"
809  field for the Sources deemed to be appropriate for Psf
810  matching
811  """
812  if candidateList is None:
813  candidateList = self.getSelectSources(scienceExposure)
814 
815  if len(candidateList) < 1:
816  raise RuntimeError("No candidates in candidateList")
817 
818  listTypes = set(type(x) for x in candidateList)
819  if len(listTypes) > 1:
820  raise RuntimeError("Candidate list contains mixed types: %s" % [l for l in listTypes])
821 
822  if not isinstance(candidateList[0], afwTable.SourceRecord):
823  try:
824  candidateList[0].getSource()
825  except Exception as e:
826  raise RuntimeError("Candidate List is of type: %s. " % (type(candidateList[0])) +
827  "Can only make candidate list from list of afwTable.SourceRecords, " +
828  "measAlg.PsfCandidateF or other type with a getSource() method: %s" % (e))
829  candidateList = [c.getSource() for c in candidateList]
830 
831  candidateList = diffimTools.sourceToFootprintList(candidateList,
832  templateExposure, scienceExposure,
833  kernelSize,
834  self.kConfig.detectionConfig,
835  self.log)
836  if len(candidateList) == 0:
837  raise RuntimeError("Cannot find any objects suitable for KernelCandidacy")
838 
839  return candidateList
840 
841  def _adaptCellSize(self, candidateList):
842  """NOT IMPLEMENTED YET"""
843  return self.kConfig.sizeCellX, self.kConfig.sizeCellY
844 
845  def _buildCellSet(self, templateMaskedImage, scienceMaskedImage, candidateList):
846  """Build a SpatialCellSet for use with the solve method
847 
848  Parameters
849  ----------
850  templateMaskedImage : `lsst.afw.image.MaskedImage`
851  MaskedImage to PSF-matched to scienceMaskedImage
852  scienceMaskedImage : `lsst.afw.image.MaskedImage`
853  reference MaskedImage
854  candidateList : `list`
855  a list of footprints/maskedImages for kernel candidates;
856 
857  - Currently supported: list of Footprints or measAlg.PsfCandidateF
858 
859  Returns
860  -------
861  kernelCellSet :
862  a SpatialCellSet for use with self._solve
863  """
864  if not candidateList:
865  raise RuntimeError("Candidate list must be populated by makeCandidateList")
866 
867  sizeCellX, sizeCellY = self._adaptCellSize(candidateList)
868 
869  # Object to store the KernelCandidates for spatial modeling
870  kernelCellSet = afwMath.SpatialCellSet(templateMaskedImage.getBBox(),
871  sizeCellX, sizeCellY)
872 
873  policy = pexConfig.makePolicy(self.kConfig)
874  # Place candidates within the spatial grid
875  for cand in candidateList:
876  bbox = cand['footprint'].getBBox()
877 
878  tmi = afwImage.MaskedImageF(templateMaskedImage, bbox)
879  smi = afwImage.MaskedImageF(scienceMaskedImage, bbox)
880  cand = diffimLib.makeKernelCandidate(cand['source'], tmi, smi, policy)
881 
882  self.log.debug("Candidate %d at %f, %f", cand.getId(), cand.getXCenter(), cand.getYCenter())
883  kernelCellSet.insertCandidate(cand)
884 
885  return kernelCellSet
886 
887  def _validateSize(self, templateMaskedImage, scienceMaskedImage):
888  """Return True if two image-like objects are the same size
889  """
890  return templateMaskedImage.getDimensions() == scienceMaskedImage.getDimensions()
891 
892  def _validateWcs(self, templateExposure, scienceExposure):
893  """Return True if the WCS of the two Exposures have the same origin and extent
894  """
895  templateWcs = templateExposure.getWcs()
896  scienceWcs = scienceExposure.getWcs()
897  templateBBox = templateExposure.getBBox()
898  scienceBBox = scienceExposure.getBBox()
899 
900  # LLC
901  templateOrigin = templateWcs.pixelToSky(afwGeom.Point2D(templateBBox.getBegin()))
902  scienceOrigin = scienceWcs.pixelToSky(afwGeom.Point2D(scienceBBox.getBegin()))
903 
904  # URC
905  templateLimit = templateWcs.pixelToSky(afwGeom.Point2D(templateBBox.getEnd()))
906  scienceLimit = scienceWcs.pixelToSky(afwGeom.Point2D(scienceBBox.getEnd()))
907 
908  self.log.info("Template Wcs : %f,%f -> %f,%f",
909  templateOrigin[0], templateOrigin[1],
910  templateLimit[0], templateLimit[1])
911  self.log.info("Science Wcs : %f,%f -> %f,%f",
912  scienceOrigin[0], scienceOrigin[1],
913  scienceLimit[0], scienceLimit[1])
914 
915  templateBBox = afwGeom.Box2D(templateOrigin.getPosition(afwGeom.degrees),
916  templateLimit.getPosition(afwGeom.degrees))
917  scienceBBox = afwGeom.Box2D(scienceOrigin.getPosition(afwGeom.degrees),
918  scienceLimit.getPosition(afwGeom.degrees))
919  if not (templateBBox.overlaps(scienceBBox)):
920  raise RuntimeError("Input images do not overlap at all")
921 
922  if ((templateOrigin != scienceOrigin) or
923  (templateLimit != scienceLimit) or
924  (templateExposure.getDimensions() != scienceExposure.getDimensions())):
925  return False
926  return True
927 
928 
929 subtractAlgorithmRegistry = pexConfig.makeRegistry(
930  doc="A registry of subtraction algorithms for use as a subtask in imageDifference",
931 )
932 
933 subtractAlgorithmRegistry.register('al', ImagePsfMatchTask)
A floating-point coordinate rectangle geometry.
Definition: Box.h:294
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
def makeKernelBasisList(config, targetFwhmPix=None, referenceFwhmPix=None, basisDegGauss=None, metadata=None)
def getSelectSources(self, exposure, sigma=None, doSmooth=True, idFactory=None)
def _validateSize(self, templateMaskedImage, scienceMaskedImage)
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations...
daf::base::PropertySet * set
Definition: fits.cc:832
def _buildCellSet(self, templateMaskedImage, scienceMaskedImage, candidateList)
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
Definition: Exposure.h:446
A collection of SpatialCells covering an entire image.
Definition: SpatialCell.h:387
def matchMaskedImages(self, templateMaskedImage, scienceMaskedImage, candidateList, templateFwhmPix=None, scienceFwhmPix=None)
table::Key< int > type
Definition: Detector.cc:164
def _solve(self, kernelCellSet, basisList, returnOnExcept=False)
Definition: psfMatch.py:869
def matchExposures(self, templateExposure, scienceExposure, templateFwhmPix=None, scienceFwhmPix=None, candidateList=None, doWarping=True, convolveTemplate=True)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:168
def makeCandidateList(self, templateExposure, scienceExposure, kernelSize, candidateList=None)
def subtractMaskedImages(self, templateMaskedImage, scienceMaskedImage, candidateList, templateFwhmPix=None, scienceFwhmPix=None)
Record class that contains measurements made on a single exposure.
Definition: Source.h:82
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, bool doNormalize, bool doCopyEdge=false)
Old, deprecated version of convolve.
def subtractExposures(self, templateExposure, scienceExposure, templateFwhmPix=None, scienceFwhmPix=None, candidateList=None, doWarping=True, convolveTemplate=True)
def _validateWcs(self, templateExposure, scienceExposure)