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
imageDifference.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 #
3 # LSST Data Management System
4 # Copyright 2012 LSST Corporation.
5 #
6 # This product includes software developed by the
7 # LSST Project (http://www.lsst.org/).
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the LSST License Statement and
20 # the GNU General Public License along with this program. If not,
21 # see <http://www.lsstcorp.org/LegalNotices/>.
22 #
23 import math
24 import random
25 import numpy
26 
27 import lsst.pex.config as pexConfig
28 import lsst.pipe.base as pipeBase
29 import lsst.daf.base as dafBase
30 import lsst.afw.geom as afwGeom
31 import lsst.afw.image as afwImage
32 import lsst.afw.math as afwMath
33 import lsst.afw.table as afwTable
34 import lsst.meas.astrom as measAstrom
35 from lsst.pipe.tasks.registerImage import RegisterTask
36 from lsst.meas.algorithms import SourceDetectionTask, SourceMeasurementTask, \
37  starSelectorRegistry, PsfAttributes, SingleGaussianPsf
38 from lsst.ip.diffim import ImagePsfMatchTask, DipoleMeasurementTask, DipoleAnalysis, \
39  SourceFlagChecker, KernelCandidateF, cast_KernelCandidateF, makeKernelBasisList, \
40  KernelCandidateQa, DiaCatalogSourceSelector, DiaCatalogSourceSelectorConfig
41 import lsst.ip.diffim.diffimTools as diffimTools
42 
43 FwhmPerSigma = 2 * math.sqrt(2 * math.log(2))
44 IqrToSigma = 0.741
45 
46 class ImageDifferenceConfig(pexConfig.Config):
47  """Config for ImageDifferenceTask
48  """
49  doAddCalexpBackground = pexConfig.Field(dtype=bool, default=True,
50  doc="Add background to calexp before processing it. Useful as ipDiffim does background matching.")
51  doUseRegister = pexConfig.Field(dtype=bool, default=True,
52  doc="Use image-to-image registration to align template with science image")
53  doDebugRegister = pexConfig.Field(dtype=bool, default=False,
54  doc="Writing debugging data for doUseRegister")
55  doSelectSources = pexConfig.Field(dtype=bool, default=True,
56  doc="Select stars to use for kernel fitting")
57  doSelectDcrCatalog = pexConfig.Field(dtype=bool, default=False,
58  doc="Select stars of extreme color as part of the control sample")
59  doSelectVariableCatalog = pexConfig.Field(dtype=bool, default=False,
60  doc="Select stars that are variable to be part of the control sample")
61  doSubtract = pexConfig.Field(dtype=bool, default=True, doc="Compute subtracted exposure?")
62  doPreConvolve = pexConfig.Field(dtype=bool, default=True,
63  doc="Convolve science image by its PSF before PSF-matching?")
64  useGaussianForPreConvolution = pexConfig.Field(dtype=bool, default=True,
65  doc="Use a simple gaussian PSF model for pre-convolution (else use fit PSF)? "
66  "Ignored if doPreConvolve false.")
67  doDetection = pexConfig.Field(dtype=bool, default=True, doc="Detect sources?")
68  doMerge = pexConfig.Field(dtype=bool, default=True,
69  doc="Merge positive and negative diaSources with grow radius set by growFootprint")
70  doMatchSources = pexConfig.Field(dtype=bool, default=True,
71  doc="Match diaSources with input calexp sources and ref catalog sources")
72  doMeasurement = pexConfig.Field(dtype=bool, default=True, doc="Measure diaSources?")
73  doWriteSubtractedExp = pexConfig.Field(dtype=bool, default=True, doc="Write difference exposure?")
74  doWriteMatchedExp = pexConfig.Field(dtype=bool, default=False,
75  doc="Write warped and PSF-matched template coadd exposure?")
76  doWriteSources = pexConfig.Field(dtype=bool, default=True, doc="Write sources?")
77  doAddMetrics = pexConfig.Field(dtype=bool, default=True,
78  doc="Add columns to the source table to hold analysis metrics?")
79 
80  coaddName = pexConfig.Field(
81  doc="coadd name: typically one of deep or goodSeeing",
82  dtype=str,
83  default="deep",
84  )
85  convolveTemplate = pexConfig.Field(
86  doc="Which image gets convolved (default = template)",
87  dtype=bool,
88  default=True
89  )
90  astrometer = pexConfig.ConfigurableField(
91  target = measAstrom.AstrometryTask,
92  doc = "astrometry task; used to match sources to reference objects, but not to fit a WCS",
93  )
94  sourceSelector = starSelectorRegistry.makeField("Source selection algorithm", default="diacatalog")
95  subtract = pexConfig.ConfigurableField(
96  target=ImagePsfMatchTask,
97  doc="Warp and PSF match template to exposure, then subtract",
98  )
99  detection = pexConfig.ConfigurableField(
100  target=SourceDetectionTask,
101  doc="Low-threshold detection for final measurement",
102  )
103  dipoleMeasurement = pexConfig.ConfigurableField(
104  target=DipoleMeasurementTask,
105  doc="Final source measurement on low-threshold detections; dipole fitting enabled.",
106  )
107  measurement = pexConfig.ConfigurableField(
108  target=SourceMeasurementTask,
109  doc="Final source measurement on low-threshold detections; dipole fitting NOT enabled",
110  )
111  controlStepSize = pexConfig.Field(
112  doc="What step size (every Nth one) to select a control sample from the kernelSources",
113  dtype=int,
114  default=5
115  )
116  controlRandomSeed = pexConfig.Field(
117  doc="Random seed for shuffing the control sample",
118  dtype=int,
119  default=10
120  )
121  register = pexConfig.ConfigurableField(
122  target=RegisterTask,
123  doc="Task to enable image-to-image image registration (warping)",
124  )
125 
126  templateBorderSize = pexConfig.Field(dtype=int, default=10,
127  doc="Number of pixels to grow the requested template image to account for warping")
128 
129  templateSipOrder = pexConfig.Field(dtype=int, default=2,
130  doc="Sip Order for fitting the Template Wcs (default is too high, overfitting)")
131 
132  growFootprint = pexConfig.Field(dtype=int, default=2,
133  doc="Grow positive and negative footprints by this amount before merging")
134 
135  diaSourceMatchRadius = pexConfig.Field(dtype=float, default=0.5,
136  doc="Match radius (in arcseconds) for DiaSource to Source association")
137 
138  maxDiaSourcesToMeasure = pexConfig.Field(dtype=int, default=200,
139  doc="Do not measure more than this many sources with dipoleMeasurement, use measurement")
140 
141  def setDefaults(self):
142  # Set default source selector and configure defaults for that one and some common alternatives
143  self.sourceSelector.name = "diacatalog"
144  self.sourceSelector["secondMoment"].clumpNSigma = 2.0
145  # defaults are OK for catalog and diacatalog
146 
147  # DiaSource Detection
148  self.detection.thresholdPolarity = "both"
149  self.detection.reEstimateBackground = False
150  self.detection.thresholdType = "pixel_stdev"
151 
152  # Add filtered flux measurement, the correct measurement for pre-convolved images.
153  # Enable all measurements, regardless of doPreConvolved, as it makes data harvesting easier.
154  # To change that you must modify algorithms.names in the task's applyOverrides method,
155  # after the user has set doPreConvolved.
156  self.measurement.algorithms.names.add("flux.peakLikelihood")
157  self.dipoleMeasurement.algorithms.names.add("flux.peakLikelihood")
158 
159  # For shuffling the control sample
160  random.seed(self.controlRandomSeed)
161 
162  def validate(self):
163  pexConfig.Config.validate(self)
164  if self.doMeasurement and not self.doDetection:
165  raise ValueError("Cannot run source measurement without source detection.")
166  if self.doMerge and not self.doDetection:
167  raise ValueError("Cannot run source merging without source detection.")
168 
169 class ImageDifferenceTask(pipeBase.CmdLineTask):
170  """Subtract an image from a template coadd and measure the result
171  """
172  ConfigClass = ImageDifferenceConfig
173  _DefaultName = "imageDifference"
174 
175  def __init__(self, **kwargs):
176  pipeBase.CmdLineTask.__init__(self, **kwargs)
177  self.makeSubtask("subtract")
178 
179  if self.config.doUseRegister:
180  self.makeSubtask("register")
181 
182  if self.config.doSelectSources:
183  self.sourceSelector = self.config.sourceSelector.apply()
184  self.astrometer = self.makeSubtask("astrometer")
185  self.schema = afwTable.SourceTable.makeMinimalSchema()
187  if self.config.doDetection:
188  self.makeSubtask("detection", schema=self.schema)
189  if self.config.doMeasurement:
190  self.makeSubtask("dipoleMeasurement", schema=self.schema, algMetadata=self.algMetadata)
191  self.makeSubtask("measurement", schema=afwTable.SourceTable.makeMinimalSchema(),
192  algMetadata=self.algMetadata)
193  self.schema.addField(self.dipoleMeasurement._ClassificationFlag, "F",
194  "probability of being a dipole")
195 
196  if self.config.doMatchSources:
197  self.schema.addField("refMatchId", "L", "unique id of reference catalog match")
198  self.schema.addField("srcMatchId", "L", "unique id of source match")
199 
200  @pipeBase.timeMethod
201  def run(self, sensorRef):
202  """Subtract an image from a template coadd and measure the result
203 
204  Steps include:
205  - warp template coadd to match WCS of image
206  - PSF match image to warped template
207  - subtract image from PSF-matched, warped template
208  - persist difference image
209  - detect sources
210  - measure sources
211 
212  @param sensorRef: sensor-level butler data reference, used for the following data products:
213  Input only:
214  - calexp
215  - psf
216  - ccdExposureId
217  - ccdExposureId_bits
218  - self.config.coaddName + "Coadd_skyMap"
219  - self.config.coaddName + "Coadd"
220  Input or output, depending on config:
221  - self.config.coaddName + "Diff_subtractedExp"
222  Output, depending on config:
223  - self.config.coaddName + "Diff_matchedExp"
224  - self.config.coaddName + "Diff_src"
225 
226  @return pipe_base Struct containing these fields:
227  - subtractedExposure: exposure after subtracting template;
228  the unpersisted version if subtraction not run but detection run
229  None if neither subtraction nor detection run (i.e. nothing useful done)
230  - subtractRes: results of subtraction task; None if subtraction not run
231  - sources: detected and possibly measured sources; None if detection not run
232  """
233  self.log.info("Processing %s" % (sensorRef.dataId))
234 
235  # initialize outputs and some intermediate products
236  subtractedExposure = None
237  subtractRes = None
238  selectSources = None
239  kernelSources = None
240  controlSources = None
241  diaSources = None
242 
243  # We make one IdFactory that will be used by both icSrc and src datasets;
244  # I don't know if this is the way we ultimately want to do things, but at least
245  # this ensures the source IDs are fully unique.
246  expBits = sensorRef.get("ccdExposureId_bits")
247  expId = long(sensorRef.get("ccdExposureId"))
248  idFactory = afwTable.IdFactory.makeSource(expId, 64 - expBits)
249 
250  # Retrieve the science image we wish to analyze
251  exposure = sensorRef.get("calexp", immediate=True)
252  if self.config.doAddCalexpBackground:
253  mi = exposure.getMaskedImage()
254  mi += sensorRef.get("calexpBackground").getImage()
255  if not exposure.hasPsf():
256  raise pipeBase.TaskError("Exposure has no psf")
257  sciencePsf = exposure.getPsf()
258 
259  subtractedExposureName = self.config.coaddName + "Diff_differenceExp"
260  templateExposure = None # Stitched coadd exposure
261  templateSources = None # Sources on the template image
262  if self.config.doSubtract:
263  templateExposure, templateSources = self.getTemplate(exposure, sensorRef)
264 
265  # compute scienceSigmaOrig: sigma of PSF of science image before pre-convolution
266  ctr = afwGeom.Box2D(exposure.getBBox()).getCenter()
267  psfAttr = PsfAttributes(sciencePsf, afwGeom.Point2I(ctr))
268  scienceSigmaOrig = psfAttr.computeGaussianWidth(psfAttr.ADAPTIVE_MOMENT)
269 
270  # sigma of PSF of template image before warping
271  ctr = afwGeom.Box2D(templateExposure.getBBox()).getCenter()
272  psfAttr = PsfAttributes(templateExposure.getPsf(), afwGeom.Point2I(ctr))
273  templateSigma = psfAttr.computeGaussianWidth(psfAttr.ADAPTIVE_MOMENT)
274 
275  # if requested, convolve the science exposure with its PSF
276  # (properly, this should be a cross-correlation, but our code does not yet support that)
277  # compute scienceSigmaPost: sigma of science exposure with pre-convolution, if done,
278  # else sigma of original science exposure
279  if self.config.doPreConvolve:
280  convControl = afwMath.ConvolutionControl()
281  # cannot convolve in place, so make a new MI to receive convolved image
282  srcMI = exposure.getMaskedImage()
283  destMI = srcMI.Factory(srcMI.getDimensions())
284  srcPsf = sciencePsf
285  if self.config.useGaussianForPreConvolution:
286  # convolve with a simplified PSF model: a double Gaussian
287  kWidth, kHeight = sciencePsf.getLocalKernel().getDimensions()
288  preConvPsf = SingleGaussianPsf(kWidth, kHeight, scienceSigmaOrig)
289  else:
290  # convolve with science exposure's PSF model
291  preConvPsf = srcPsf
292  afwMath.convolve(destMI, srcMI, preConvPsf.getLocalKernel(), convControl)
293  exposure.setMaskedImage(destMI)
294  scienceSigmaPost = scienceSigmaOrig * math.sqrt(2)
295  else:
296  scienceSigmaPost = scienceSigmaOrig
297 
298  # If requested, find sources in the image
299  if self.config.doSelectSources:
300  if not sensorRef.datasetExists("src"):
301  self.log.warn("Src product does not exist; running detection, measurement, selection")
302  # Run own detection and measurement; necessary in nightly processing
303  selectSources = self.subtract.getSelectSources(
304  exposure,
305  sigma = scienceSigmaPost,
306  doSmooth = not self.doPreConvolve,
307  idFactory = idFactory,
308  )
309  else:
310  self.log.info("Source selection via src product")
311  # Sources already exist; for data release processing
312  selectSources = sensorRef.get("src")
313 
314  # Number of basis functions
315  nparam = len(makeKernelBasisList(self.subtract.config.kernel.active,
316  referenceFwhmPix=scienceSigmaPost * FwhmPerSigma,
317  targetFwhmPix=templateSigma * FwhmPerSigma))
318 
319  if self.config.doAddMetrics:
320  # Modify the schema of all Sources
321  kcQa = KernelCandidateQa(nparam)
322  selectSources = kcQa.addToSchema(selectSources)
323 
324  astromRet = self.astrometer.loadAndMatch(exposure=exposure, sourceCat=selectSources)
325  matches = astromRet.matches
326  kernelSources = self.sourceSelector.selectSources(exposure, selectSources, matches=matches)
327  random.shuffle(kernelSources, random.random)
328  controlSources = kernelSources[::self.config.controlStepSize]
329  kernelSources = [k for i,k in enumerate(kernelSources) if i % self.config.controlStepSize]
330 
331  if self.config.doSelectDcrCatalog:
332  redSelector = DiaCatalogSourceSelector(
333  DiaCatalogSourceSelectorConfig(grMin=self.sourceSelector.config.grMax, grMax=99.999))
334  redSources = redSelector.selectSources(exposure, selectSources, matches=matches)
335  controlSources.extend(redSources)
336 
337  blueSelector = DiaCatalogSourceSelector(
338  DiaCatalogSourceSelectorConfig(grMin=-99.999, grMax=self.sourceSelector.config.grMin))
339  blueSources = blueSelector.selectSources(exposure, selectSources, matches=matches)
340  controlSources.extend(blueSources)
341 
342  if self.config.doSelectVariableCatalog:
343  varSelector = DiaCatalogSourceSelector(
344  DiaCatalogSourceSelectorConfig(includeVariable=True))
345  varSources = varSelector.selectSources(exposure, selectSources, matches=matches)
346  controlSources.extend(varSources)
347 
348  self.log.info("Selected %d / %d sources for Psf matching (%d for control sample)"
349  % (len(kernelSources), len(selectSources), len(controlSources)))
350  allresids = {}
351  if self.config.doUseRegister:
352  self.log.info("Registering images")
353 
354  if templateSources is None:
355  # Run detection on the template, which is
356  # temporarily background-subtracted
357  templateSources = self.subtract.getSelectSources(
358  templateExposure,
359  sigma=templateSigma,
360  doSmooth=True,
361  idFactory=idFactory
362  )
363 
364  # Third step: we need to fit the relative astrometry.
365  #
366  wcsResults = self.fitAstrometry(templateSources, templateExposure, selectSources)
367  warpedExp = self.register.warpExposure(templateExposure, wcsResults.wcs,
368  exposure.getWcs(), exposure.getBBox())
369  templateExposure = warpedExp
370 
371  # Create debugging outputs on the astrometric
372  # residuals as a function of position. Persistence
373  # not yet implemented; expected on (I believe) #2636.
374  if self.config.doDebugRegister:
375  # Grab matches to reference catalog
376  srcToMatch = {x.second.getId() : x.first for x in matches}
377 
378  refCoordKey = wcsResults.matches[0].first.getTable().getCoordKey()
379  inCentroidKey = wcsResults.matches[0].second.getTable().getCentroidKey()
380  sids = [m.first.getId() for m in wcsResults.matches]
381  positions = [m.first.get(refCoordKey) for m in wcsResults.matches]
382  residuals = [m.first.get(refCoordKey).getOffsetFrom(wcsResults.wcs.pixelToSky(
383  m.second.get(inCentroidKey))) for m in wcsResults.matches]
384  allresids = dict(zip(sids, zip(positions, residuals)))
385 
386  cresiduals = [m.first.get(refCoordKey).getTangentPlaneOffset(
387  wcsResults.wcs.pixelToSky(
388  m.second.get(inCentroidKey))) for m in wcsResults.matches]
389  colors = numpy.array([-2.5*numpy.log10(srcToMatch[x].get("g"))
390  + 2.5*numpy.log10(srcToMatch[x].get("r"))
391  for x in sids if x in srcToMatch.keys()])
392  dlong = numpy.array([r[0].asArcseconds() for s,r in zip(sids, cresiduals)
393  if s in srcToMatch.keys()])
394  dlat = numpy.array([r[1].asArcseconds() for s,r in zip(sids, cresiduals)
395  if s in srcToMatch.keys()])
396  idx1 = numpy.where(colors<self.sourceSelector.config.grMin)
397  idx2 = numpy.where((colors>=self.sourceSelector.config.grMin)&
398  (colors<=self.sourceSelector.config.grMax))
399  idx3 = numpy.where(colors>self.sourceSelector.config.grMax)
400  rms1Long = IqrToSigma*(numpy.percentile(dlong[idx1],75)-numpy.percentile(dlong[idx1],25))
401  rms1Lat = IqrToSigma*(numpy.percentile(dlat[idx1],75)-numpy.percentile(dlat[idx1],25))
402  rms2Long = IqrToSigma*(numpy.percentile(dlong[idx2],75)-numpy.percentile(dlong[idx2],25))
403  rms2Lat = IqrToSigma*(numpy.percentile(dlat[idx2],75)-numpy.percentile(dlat[idx2],25))
404  rms3Long = IqrToSigma*(numpy.percentile(dlong[idx3],75)-numpy.percentile(dlong[idx3],25))
405  rms3Lat = IqrToSigma*(numpy.percentile(dlat[idx3],75)-numpy.percentile(dlat[idx3],25))
406  self.log.info("Blue star offsets'': %.3f %.3f, %.3f %.3f" % (numpy.median(dlong[idx1]),
407  rms1Long,
408  numpy.median(dlat[idx1]),
409  rms1Lat))
410  self.log.info("Green star offsets'': %.3f %.3f, %.3f %.3f" % (numpy.median(dlong[idx2]),
411  rms2Long,
412  numpy.median(dlat[idx2]),
413  rms2Lat))
414  self.log.info("Red star offsets'': %.3f %.3f, %.3f %.3f" % (numpy.median(dlong[idx3]),
415  rms3Long,
416  numpy.median(dlat[idx3]),
417  rms3Lat))
418 
419  self.metadata.add("RegisterBlueLongOffsetMedian", numpy.median(dlong[idx1]))
420  self.metadata.add("RegisterGreenLongOffsetMedian", numpy.median(dlong[idx2]))
421  self.metadata.add("RegisterRedLongOffsetMedian", numpy.median(dlong[idx3]))
422  self.metadata.add("RegisterBlueLongOffsetStd", rms1Long)
423  self.metadata.add("RegisterGreenLongOffsetStd", rms2Long)
424  self.metadata.add("RegisterRedLongOffsetStd", rms3Long)
425 
426  self.metadata.add("RegisterBlueLatOffsetMedian", numpy.median(dlat[idx1]))
427  self.metadata.add("RegisterGreenLatOffsetMedian", numpy.median(dlat[idx2]))
428  self.metadata.add("RegisterRedLatOffsetMedian", numpy.median(dlat[idx3]))
429  self.metadata.add("RegisterBlueLatOffsetStd", rms1Lat)
430  self.metadata.add("RegisterGreenLatOffsetStd", rms2Lat)
431  self.metadata.add("RegisterRedLatOffsetStd", rms3Lat)
432 
433  # warp template exposure to match exposure,
434  # PSF match template exposure to exposure,
435  # then return the difference
436 
437  #Return warped template... Construct sourceKernelCand list after subtract
438  self.log.info("Subtracting images")
439  subtractRes = self.subtract.subtractExposures(
440  templateExposure=templateExposure,
441  scienceExposure=exposure,
442  scienceFwhmPix=scienceSigmaPost * FwhmPerSigma,
443  templateFwhmPix=templateSigma * FwhmPerSigma,
444  candidateList=kernelSources,
445  convolveTemplate=self.config.convolveTemplate,
446  doWarping=not self.config.doUseRegister
447  )
448  subtractedExposure = subtractRes.subtractedExposure
449 
450  if self.config.doWriteMatchedExp:
451  sensorRef.put(subtractRes.matchedExposure, self.config.coaddName + "Diff_matchedExp")
452 
453  if self.config.doDetection:
454  self.log.info("Running diaSource detection")
455  if subtractedExposure is None:
456  subtractedExposure = sensorRef.get(subtractedExposureName)
457 
458  # Get Psf from the appropriate input image if it doesn't exist
459  if not subtractedExposure.hasPsf():
460  if self.config.convolveTemplate:
461  subtractedExposure.setPsf(exposure.getPsf())
462  else:
463  if templateExposure is None:
464  templateExposure, templateSources = self.getTemplate(exposure, sensorRef)
465  subtractedExposure.setPsf(templateExposure.getPsf())
466 
467  # Erase existing detection mask planes
468  mask = subtractedExposure.getMaskedImage().getMask()
469  mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
470 
471  table = afwTable.SourceTable.make(self.schema, idFactory)
472  table.setMetadata(self.algMetadata)
473  results = self.detection.makeSourceCatalog(
474  table=table,
475  exposure=subtractedExposure,
476  doSmooth=not self.config.doPreConvolve
477  )
478 
479  if self.config.doMerge:
480  fpSet = results.fpSets.positive
481  fpSet.merge(results.fpSets.negative, self.config.growFootprint,
482  self.config.growFootprint, False)
483  diaSources = afwTable.SourceCatalog(table)
484  fpSet.makeSources(diaSources)
485  self.log.info("Merging detections into %d sources" % (len(diaSources)))
486  else:
487  diaSources = results.sources
488 
489  if self.config.doMeasurement:
490  if len(diaSources) < self.config.maxDiaSourcesToMeasure:
491  self.log.info("Running diaSource dipole measurement")
492  self.dipoleMeasurement.run(subtractedExposure, diaSources)
493  else:
494  self.log.info("Running diaSource measurement")
495  self.measurement.run(subtractedExposure, diaSources)
496 
497  # Match with the calexp sources if possible
498  if self.config.doMatchSources:
499  if sensorRef.datasetExists("src"):
500  # Create key,val pair where key=diaSourceId and val=sourceId
501  matchRadAsec = self.config.diaSourceMatchRadius
502  matchRadPixel = matchRadAsec / exposure.getWcs().pixelScale().asArcseconds()
503 
504  srcMatches = afwTable.matchXy(sensorRef.get("src"), diaSources, matchRadPixel, True)
505  srcMatchDict = dict([(srcMatch.second.getId(), srcMatch.first.getId()) for
506  srcMatch in srcMatches])
507  self.log.info("Matched %d / %d diaSources to sources" % (len(srcMatchDict),
508  len(diaSources)))
509  else:
510  self.log.warn("Src product does not exist; cannot match with diaSources")
511  srcMatchDict = {}
512 
513  # Create key,val pair where key=diaSourceId and val=refId
514  refAstromConfig = measAstrom.AstrometryConfig()
515  refAstromConfig.matcher.maxMatchDistArcSec = matchRadAsec
516  refAstrometer = measAstrom.AstrometryTask(refAstromConfig)
517  astromRet = refAstrometer.run(exposure=exposure, sourceCat=diaSources)
518  refMatches = astromRet.matches
519  if refMatches is None:
520  self.log.warn("No diaSource matches with reference catalog")
521  refMatchDict = {}
522  else:
523  self.log.info("Matched %d / %d diaSources to reference catalog" % (len(refMatches),
524  len(diaSources)))
525  refMatchDict = dict([(refMatch.second.getId(), refMatch.first.getId()) for \
526  refMatch in refMatches])
527 
528  # Assign source Ids
529  for diaSource in diaSources:
530  sid = diaSource.getId()
531  if srcMatchDict.has_key(sid):
532  diaSource.set("srcMatchId", srcMatchDict[sid])
533  if refMatchDict.has_key(sid):
534  diaSource.set("refMatchId", refMatchDict[sid])
535 
536  if diaSources is not None and self.config.doWriteSources:
537  sensorRef.put(diaSources, self.config.coaddName + "Diff_diaSrc")
538 
539  if self.config.doAddMetrics and self.config.doSelectSources:
540  self.log.info("Evaluating metrics and control sample")
541 
542  kernelCandList = []
543  for cell in subtractRes.kernelCellSet.getCellList():
544  for cand in cell.begin(False): # include bad candidates
545  kernelCandList.append(cast_KernelCandidateF(cand))
546 
547  # Get basis list to build control sample kernels
548  basisList = afwMath.cast_LinearCombinationKernel(
549  kernelCandList[0].getKernel(KernelCandidateF.ORIG)).getKernelList()
550 
551  controlCandList = \
552  diffimTools.sourceTableToCandidateList(controlSources,
553  subtractRes.warpedExposure, exposure,
554  self.config.subtract.kernel.active,
555  self.config.subtract.kernel.active.detectionConfig,
556  self.log, doBuild=True, basisList=basisList)
557 
558  kcQa.apply(kernelCandList, subtractRes.psfMatchingKernel, subtractRes.backgroundModel,
559  dof=nparam)
560  kcQa.apply(controlCandList, subtractRes.psfMatchingKernel, subtractRes.backgroundModel)
561 
562  if self.config.doDetection:
563  kcQa.aggregate(selectSources, self.metadata, allresids, diaSources)
564  else:
565  kcQa.aggregate(selectSources, self.metadata, allresids)
566 
567  sensorRef.put(selectSources, self.config.coaddName + "Diff_kernelSrc")
568 
569  if self.config.doWriteSubtractedExp:
570  sensorRef.put(subtractedExposure, subtractedExposureName)
571 
572  self.runDebug(exposure, subtractRes, selectSources, kernelSources, diaSources)
573  return pipeBase.Struct(
574  subtractedExposure=subtractedExposure,
575  subtractRes=subtractRes,
576  sources=diaSources,
577  )
578 
579  def fitAstrometry(self, templateSources, templateExposure, selectSources):
580  """Fit the relative astrometry between templateSources and selectSources
581 
582  @todo remove this method. It originally fit a new WCS to the template before calling register.run
583  because our TAN-SIP fitter behaved badly for points far from CRPIX, but that's been fixed.
584  It remains because a subtask overrides it.
585  """
586  results = self.register.run(templateSources, templateSources.getWcs(),
587  templateExposure.getBBox(), selectSources)
588  return results
589 
590  def runDebug(self, exposure, subtractRes, selectSources, kernelSources, diaSources):
591  import lsstDebug
592  import lsst.afw.display.ds9 as ds9
593  display = lsstDebug.Info(__name__).display
594  showSubtracted = lsstDebug.Info(__name__).showSubtracted
595  showPixelResiduals = lsstDebug.Info(__name__).showPixelResiduals
596  showDiaSources = lsstDebug.Info(__name__).showDiaSources
597  showDipoles = lsstDebug.Info(__name__).showDipoles
598  maskTransparency = lsstDebug.Info(__name__).maskTransparency
599  if not maskTransparency:
600  maskTransparency = 0
601  ds9.setMaskTransparency(maskTransparency)
602 
603  if display and showSubtracted:
604  ds9.mtv(subtractRes.subtractedExposure, frame=lsstDebug.frame, title="Subtracted image")
605  mi = subtractRes.subtractedExposure.getMaskedImage()
606  x0, y0 = mi.getX0(), mi.getY0()
607  with ds9.Buffering():
608  for s in diaSources:
609  x, y = s.getX() - x0, s.getY() - y0
610  ctype = "red" if s.get("flags.negative") else "yellow"
611  if (s.get("flags.pixel.interpolated.center") or s.get("flags.pixel.saturated.center") or
612  s.get("flags.pixel.cr.center")):
613  ptype = "x"
614  elif (s.get("flags.pixel.interpolated.any") or s.get("flags.pixel.saturated.any") or
615  s.get("flags.pixel.cr.any")):
616  ptype = "+"
617  else:
618  ptype = "o"
619  ds9.dot(ptype, x, y, size=4, frame=lsstDebug.frame, ctype=ctype)
620  lsstDebug.frame += 1
621 
622  if display and showPixelResiduals and selectSources:
623  import lsst.ip.diffim.utils as diUtils
624  nonKernelSources = []
625  for source in selectSources:
626  if not source in kernelSources:
627  nonKernelSources.append(source)
628 
629  diUtils.plotPixelResiduals(exposure,
630  subtractRes.warpedExposure,
631  subtractRes.subtractedExposure,
632  subtractRes.kernelCellSet,
633  subtractRes.psfMatchingKernel,
634  subtractRes.backgroundModel,
635  nonKernelSources,
636  self.subtract.config.kernel.active.detectionConfig,
637  origVariance=False)
638  diUtils.plotPixelResiduals(exposure,
639  subtractRes.warpedExposure,
640  subtractRes.subtractedExposure,
641  subtractRes.kernelCellSet,
642  subtractRes.psfMatchingKernel,
643  subtractRes.backgroundModel,
644  nonKernelSources,
645  self.subtract.config.kernel.active.detectionConfig,
646  origVariance=True)
647  if display and showDiaSources:
648  import lsst.ip.diffim.utils as diUtils
649  flagChecker = SourceFlagChecker(diaSources)
650  isFlagged = [flagChecker(x) for x in diaSources]
651  isDipole = [x.get("classification.dipole") for x in diaSources]
652  diUtils.showDiaSources(diaSources, subtractRes.subtractedExposure, isFlagged, isDipole,
653  frame=lsstDebug.frame)
654  lsstDebug.frame += 1
655 
656  if display and showDipoles:
657  DipoleAnalysis().displayDipoles(subtractRes.subtractedExposure, diaSources,
658  frame=lsstDebug.frame)
659  lsstDebug.frame += 1
660 
661  def getTemplate(self, exposure, sensorRef):
662  """Return a template coadd exposure that overlaps the exposure
663 
664  @param[in] exposure: exposure
665  @param[in] sensorRef: a Butler data reference that can be used to obtain coadd data
666 
667  @return coaddExposure: a template coadd exposure assembled out of patches
668  """
669  skyMap = sensorRef.get(datasetType=self.config.coaddName + "Coadd_skyMap")
670  expWcs = exposure.getWcs()
671  expBoxD = afwGeom.Box2D(exposure.getBBox())
672  expBoxD.grow(self.config.templateBorderSize)
673  ctrSkyPos = expWcs.pixelToSky(expBoxD.getCenter())
674  tractInfo = skyMap.findTract(ctrSkyPos)
675  self.log.info("Using skyMap tract %s" % (tractInfo.getId(),))
676  skyCorners = [expWcs.pixelToSky(pixPos) for pixPos in expBoxD.getCorners()]
677  patchList = tractInfo.findPatchList(skyCorners)
678 
679  if not patchList:
680  raise RuntimeError("No suitable tract found")
681  self.log.info("Assembling %s coadd patches" % (len(patchList),))
682 
683  # compute coadd bbox
684  coaddWcs = tractInfo.getWcs()
685  coaddBBox = afwGeom.Box2D()
686  for skyPos in skyCorners:
687  coaddBBox.include(coaddWcs.skyToPixel(skyPos))
688  coaddBBox = afwGeom.Box2I(coaddBBox)
689  self.log.info("exposure dimensions=%s; coadd dimensions=%s" % \
690  (exposure.getDimensions(), coaddBBox.getDimensions()))
691 
692  # assemble coadd exposure from subregions of patches
693  coaddExposure = afwImage.ExposureF(coaddBBox, coaddWcs)
694  coaddExposure.getMaskedImage().set(numpy.nan, afwImage.MaskU.getPlaneBitMask("NO_DATA"), numpy.nan)
695  nPatchesFound = 0
696  coaddFilter = None
697  coaddPsf = None
698  coaddSources = None
699  for patchInfo in patchList:
700  patchSubBBox = patchInfo.getOuterBBox()
701  patchSubBBox.clip(coaddBBox)
702  patchArgDict = dict(
703  datasetType=self.config.coaddName + "Coadd_sub",
704  bbox=patchSubBBox,
705  tract=tractInfo.getId(),
706  patch="%s,%s" % (patchInfo.getIndex()[0], patchInfo.getIndex()[1]),
707  )
708  if patchSubBBox.isEmpty():
709  self.log.info("skip tract=%(tract)s; no overlapping pixels" % patchArgDict)
710  continue
711  if not sensorRef.datasetExists(**patchArgDict):
712  self.log.warn("%(datasetType)s, tract=%(tract)s, patch=%(patch)s does not exist" \
713  % patchArgDict)
714  continue
715 
716  nPatchesFound += 1
717  self.log.info("Reading patch %s" % patchArgDict)
718  coaddPatch = sensorRef.get(patchArgDict, immediate=True)
719  coaddView = afwImage.MaskedImageF(coaddExposure.getMaskedImage(), coaddPatch.getBBox())
720  coaddView <<= coaddPatch.getMaskedImage()
721  if coaddFilter is None:
722  coaddFilter = coaddPatch.getFilter()
723 
724  # Retrieve the PSF for this coadd tract, if not already retrieved
725  if coaddPsf is None and coaddPatch.hasPsf():
726  coaddPsf = coaddPatch.getPsf()
727 
728  if nPatchesFound == 0:
729  raise RuntimeError("No patches found!")
730 
731  if coaddPsf is None:
732  raise RuntimeError("No coadd Psf found!")
733 
734  coaddExposure.setPsf(coaddPsf)
735  coaddExposure.setFilter(coaddFilter)
736  return coaddExposure, coaddSources
737 
738  def _getConfigName(self):
739  """Return the name of the config dataset
740  """
741  return "%sDiff_config" % (self.config.coaddName,)
742 
743  def _getMetadataName(self):
744  """Return the name of the metadata dataset
745  """
746  return "%sDiff_metadata" % (self.config.coaddName,)
747 
748  def getSchemaCatalogs(self):
749  """Return a dict of empty catalogs for each catalog dataset produced by this task."""
750  diaSrc = afwTable.SourceCatalog(self.schema)
751  diaSrc.getTable().setMetadata(self.algMetadata)
752  return {self.config.coaddName + "Diff_diaSrc": diaSrc}
753 
754  @classmethod
756  """Create an argument parser
757  """
758  parser = pipeBase.ArgumentParser(name=cls._DefaultName)
759  parser.add_id_argument("--id", "calexp", help="data ID, e.g. --id visit=12345 ccd=1,2")
760  return parser
761 
763  winter2013TemplateId = pexConfig.Field(dtype=int, default=88868666,
764  doc="88868666 for sparse data; 22222200 (g) and 11111100 (i) for dense data")
765  winter2013WcsShift = pexConfig.Field(dtype=float, default=0.0,
766  doc="Shift stars going into RegisterTask by this amount")
767  winter2013WcsRms = pexConfig.Field(dtype=float, default=0.0,
768  doc="Perturb stars going into RegisterTask by this amount")
769 
771  ConfigClass = Winter2013ImageDifferenceConfig
772  _DefaultName = "winter2013ImageDifference"
773 
774  def __init__(self, **kwargs):
775  ImageDifferenceTask.__init__(self, **kwargs)
776 
777  def getTemplate(self, exposure, sensorRef):
778  """Return a template coadd exposure that overlaps the exposure
779 
780  @param[in] exposure: exposure
781  @param[in] sensorRef: a Butler data reference that can be used to obtain coadd data
782 
783  @return coaddExposure: a template coadd exposure assembled out of patches
784  """
785  # Using a deep simulated calexp instead of a coadd
786  self.log.warn("USING WINTER2013 : DEEP CALEXP AS TEMPLATE")
787  templateId = type(sensorRef.dataId)(sensorRef.dataId)
788  templateId["visit"] = self.config.winter2013TemplateId
789  template = sensorRef.getButler().get(datasetType="calexp", dataId=templateId)
790  if self.config.doAddCalexpBackground:
791  templateBg = sensorRef.getButler().get(datasetType="calexpBackground", dataId=templateId)
792  mi = template.getMaskedImage()
793  mi += templateBg
794  if not template.hasPsf():
795  raise pipeBase.TaskError("Template has no psf")
796  templateSources = sensorRef.getButler().get(datasetType="src", dataId=templateId)
797  return template, templateSources
798 
799 
800  def fitAstrometry(self, templateSources, templateExposure, selectSources):
801  """Fit the relative astrometry between templateSources and selectSources"""
802  if self.config.winter2013WcsShift > 0.0:
803  offset = afwGeom.Extent2D(self.config.winter2013WcsShift,
804  self.config.winter2013WcsShift)
805  cKey = templateSources[0].getTable().getCentroidKey()
806  for source in templateSources:
807  centroid = source.get(cKey)
808  source.set(cKey, centroid+offset)
809  elif self.config.winter2013WcsRms > 0.0:
810  cKey = templateSources[0].getTable().getCentroidKey()
811  for source in templateSources:
812  offset = afwGeom.Extent2D(self.config.winter2013WcsRms*numpy.random.normal(),
813  self.config.winter2013WcsRms*numpy.random.normal())
814  centroid = source.get(cKey)
815  source.set(cKey, centroid+offset)
816 
817  results = self.register.run(templateSources, templateExposure.getWcs(),
818  templateExposure.getBBox(), selectSources)
819  return results
SourceMatchVector matchXy(SourceCatalog const &cat, double radius, bool symmetric=true)
Class for storing ordered metadata with comments.
Definition: PropertyList.h:81
Parameters to control convolution.
Definition: ConvolveImage.h:58
An integer coordinate rectangle.
Definition: Box.h:53
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
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, bool doNormalize, bool doCopyEdge=false)
Old, deprecated version of convolve.
A floating-point coordinate rectangle geometry.
Definition: Box.h:271