LSSTApplications  10.0+286,10.0+36,10.0+46,10.0-2-g4f67435,10.1+152,10.1+37,11.0,11.0+1,11.0-1-g47edd16,11.0-1-g60db491,11.0-1-g7418c06,11.0-2-g04d2804,11.0-2-g68503cd,11.0-2-g818369d,11.0-2-gb8b8ce7
LSSTDataManagementBasePackage
Public Member Functions | Public Attributes | Static Public Attributes | Private Member Functions | Static Private Attributes | List of all members
lsst.pipe.tasks.imageDifference.ImageDifferenceTask Class Reference
Inheritance diagram for lsst.pipe.tasks.imageDifference.ImageDifferenceTask:
lsst.pipe.tasks.imageDifference.Winter2013ImageDifferenceTask

Public Member Functions

def __init__
 
def run
 
def fitAstrometry
 
def runDebug
 
def getTemplate
 
def getSchemaCatalogs
 

Public Attributes

 sourceSelector
 
 astrometer
 
 schema
 
 algMetadata
 

Static Public Attributes

 ConfigClass = ImageDifferenceConfig
 

Private Member Functions

def _getConfigName
 
def _getMetadataName
 
def _makeArgumentParser
 

Static Private Attributes

string _DefaultName = "imageDifference"
 

Detailed Description

Subtract an image from a template coadd and measure the result

Definition at line 169 of file imageDifference.py.

Constructor & Destructor Documentation

def lsst.pipe.tasks.imageDifference.ImageDifferenceTask.__init__ (   self,
  kwargs 
)

Definition at line 175 of file imageDifference.py.

175 
176  def __init__(self, **kwargs):
177  pipeBase.CmdLineTask.__init__(self, **kwargs)
178  self.makeSubtask("subtract")
179 
180  if self.config.doUseRegister:
181  self.makeSubtask("register")
182 
183  if self.config.doSelectSources:
184  self.sourceSelector = self.config.sourceSelector.apply()
185  self.astrometer = self.makeSubtask("astrometer")
186  self.schema = afwTable.SourceTable.makeMinimalSchema()
187  self.algMetadata = dafBase.PropertyList()
188  if self.config.doDetection:
189  self.makeSubtask("detection", schema=self.schema)
190  if self.config.doMeasurement:
191  self.makeSubtask("dipoleMeasurement", schema=self.schema, algMetadata=self.algMetadata)
192  self.makeSubtask("measurement", schema=afwTable.SourceTable.makeMinimalSchema(),
193  algMetadata=self.algMetadata)
194  self.schema.addField(self.dipoleMeasurement._ClassificationFlag, "F",
195  "probability of being a dipole")
196 
197  if self.config.doMatchSources:
198  self.schema.addField("refMatchId", "L", "unique id of reference catalog match")
199  self.schema.addField("srcMatchId", "L", "unique id of source match")
lsst::daf::base::PropertyList
Class for storing ordered metadata with comments.
Definition: PropertyList.h:81

Member Function Documentation

def lsst.pipe.tasks.imageDifference.ImageDifferenceTask._getConfigName (   self)
private 
Return the name of the config dataset

Definition at line 738 of file imageDifference.py.

738 
739  def _getConfigName(self):
740  """Return the name of the config dataset
741  """
742  return "%sDiff_config" % (self.config.coaddName,)
def lsst.pipe.tasks.imageDifference.ImageDifferenceTask._getMetadataName (   self)
private 
Return the name of the metadata dataset

Definition at line 743 of file imageDifference.py.

743 
744  def _getMetadataName(self):
745  """Return the name of the metadata dataset
746  """
747  return "%sDiff_metadata" % (self.config.coaddName,)
def lsst.pipe.tasks.imageDifference.ImageDifferenceTask._makeArgumentParser (   cls)
private 
Create an argument parser

Definition at line 755 of file imageDifference.py.

755 
756  def _makeArgumentParser(cls):
757  """Create an argument parser
758  """
759  parser = pipeBase.ArgumentParser(name=cls._DefaultName)
760  parser.add_id_argument("--id", "calexp", help="data ID, e.g. --id visit=12345 ccd=1,2")
761  return parser
def lsst.pipe.tasks.imageDifference.ImageDifferenceTask.fitAstrometry (   self,
  templateSources,
  templateExposure,
  selectSources 
) 
Fit the relative astrometry between templateSources and selectSources

@todo remove this method. It originally fit a new WCS to the template before calling register.run
because our TAN-SIP fitter behaved badly for points far from CRPIX, but that's been fixed.
It remains because a subtask overrides it.

Definition at line 579 of file imageDifference.py.

579 
580  def fitAstrometry(self, templateSources, templateExposure, selectSources):
581  """Fit the relative astrometry between templateSources and selectSources
582 
583  @todo remove this method. It originally fit a new WCS to the template before calling register.run
584  because our TAN-SIP fitter behaved badly for points far from CRPIX, but that's been fixed.
585  It remains because a subtask overrides it.
586  """
587  results = self.register.run(templateSources, templateSources.getWcs(),
588  templateExposure.getBBox(), selectSources)
589  return results
def lsst.pipe.tasks.imageDifference.ImageDifferenceTask.getSchemaCatalogs (   self) 
Return a dict of empty catalogs for each catalog dataset produced by this task.

Definition at line 748 of file imageDifference.py.

748 
749  def getSchemaCatalogs(self):
750  """Return a dict of empty catalogs for each catalog dataset produced by this task."""
751  diaSrc = afwTable.SourceCatalog(self.schema)
752  diaSrc.getTable().setMetadata(self.algMetadata)
753  return {self.config.coaddName + "Diff_diaSrc": diaSrc}
lsst::afw::table::SortedCatalogT
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
def lsst.pipe.tasks.imageDifference.ImageDifferenceTask.getTemplate (   self,
  exposure,
  sensorRef 
) 
Return a template coadd exposure that overlaps the exposure

@param[in] exposure: exposure
@param[in] sensorRef: a Butler data reference that can be used to obtain coadd data

@return coaddExposure: a template coadd exposure assembled out of patches

Definition at line 661 of file imageDifference.py.

661 
662  def getTemplate(self, exposure, sensorRef):
663  """Return a template coadd exposure that overlaps the exposure
664 
665  @param[in] exposure: exposure
666  @param[in] sensorRef: a Butler data reference that can be used to obtain coadd data
667 
668  @return coaddExposure: a template coadd exposure assembled out of patches
669  """
670  skyMap = sensorRef.get(datasetType=self.config.coaddName + "Coadd_skyMap")
671  expWcs = exposure.getWcs()
672  expBoxD = afwGeom.Box2D(exposure.getBBox())
673  expBoxD.grow(self.config.templateBorderSize)
674  ctrSkyPos = expWcs.pixelToSky(expBoxD.getCenter())
675  tractInfo = skyMap.findTract(ctrSkyPos)
676  self.log.info("Using skyMap tract %s" % (tractInfo.getId(),))
677  skyCorners = [expWcs.pixelToSky(pixPos) for pixPos in expBoxD.getCorners()]
678  patchList = tractInfo.findPatchList(skyCorners)
679 
680  if not patchList:
681  raise RuntimeError("No suitable tract found")
682  self.log.info("Assembling %s coadd patches" % (len(patchList),))
683 
684  # compute coadd bbox
685  coaddWcs = tractInfo.getWcs()
686  coaddBBox = afwGeom.Box2D()
687  for skyPos in skyCorners:
688  coaddBBox.include(coaddWcs.skyToPixel(skyPos))
689  coaddBBox = afwGeom.Box2I(coaddBBox)
690  self.log.info("exposure dimensions=%s; coadd dimensions=%s" % \
691  (exposure.getDimensions(), coaddBBox.getDimensions()))
692 
693  # assemble coadd exposure from subregions of patches
694  coaddExposure = afwImage.ExposureF(coaddBBox, coaddWcs)
695  coaddExposure.getMaskedImage().set(numpy.nan, afwImage.MaskU.getPlaneBitMask("NO_DATA"), numpy.nan)
696  nPatchesFound = 0
697  coaddFilter = None
698  coaddPsf = None
699  coaddSources = None
700  for patchInfo in patchList:
701  patchSubBBox = patchInfo.getOuterBBox()
702  patchSubBBox.clip(coaddBBox)
703  patchArgDict = dict(
704  datasetType=self.config.coaddName + "Coadd_sub",
705  bbox=patchSubBBox,
706  tract=tractInfo.getId(),
707  patch="%s,%s" % (patchInfo.getIndex()[0], patchInfo.getIndex()[1]),
708  )
709  if patchSubBBox.isEmpty():
710  self.log.info("skip tract=%(tract)s; no overlapping pixels" % patchArgDict)
711  continue
712  if not sensorRef.datasetExists(**patchArgDict):
713  self.log.warn("%(datasetType)s, tract=%(tract)s, patch=%(patch)s does not exist" \
714  % patchArgDict)
715  continue
716 
717  nPatchesFound += 1
718  self.log.info("Reading patch %s" % patchArgDict)
719  coaddPatch = sensorRef.get(patchArgDict, immediate=True)
720  coaddView = afwImage.MaskedImageF(coaddExposure.getMaskedImage(), coaddPatch.getBBox())
721  coaddView <<= coaddPatch.getMaskedImage()
722  if coaddFilter is None:
723  coaddFilter = coaddPatch.getFilter()
724 
725  # Retrieve the PSF for this coadd tract, if not already retrieved
726  if coaddPsf is None and coaddPatch.hasPsf():
727  coaddPsf = coaddPatch.getPsf()
728 
729  if nPatchesFound == 0:
730  raise RuntimeError("No patches found!")
731 
732  if coaddPsf is None:
733  raise RuntimeError("No coadd Psf found!")
734 
735  coaddExposure.setPsf(coaddPsf)
736  coaddExposure.setFilter(coaddFilter)
737  return coaddExposure, coaddSources
lsst::afw::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:53
lsst::afw::geom::Box2D
A floating-point coordinate rectangle geometry.
Definition: Box.h:271
def lsst.pipe.tasks.imageDifference.ImageDifferenceTask.run (   self,
  sensorRef 
) 
Subtract an image from a template coadd and measure the result

Steps include:
- warp template coadd to match WCS of image
- PSF match image to warped template
- subtract image from PSF-matched, warped template
- persist difference image
- detect sources
- measure sources

@param sensorRef: sensor-level butler data reference, used for the following data products:
Input only:
- calexp
- psf
- ccdExposureId
- ccdExposureId_bits
- self.config.coaddName + "Coadd_skyMap"
- self.config.coaddName + "Coadd"
Input or output, depending on config:
- self.config.coaddName + "Diff_subtractedExp"
Output, depending on config:
- self.config.coaddName + "Diff_matchedExp"
- self.config.coaddName + "Diff_src"

@return pipe_base Struct containing these fields:
- subtractedExposure: exposure after subtracting template;
    the unpersisted version if subtraction not run but detection run
    None if neither subtraction nor detection run (i.e. nothing useful done)
- subtractRes: results of subtraction task; None if subtraction not run
- sources: detected and possibly measured sources; None if detection not run

Definition at line 201 of file imageDifference.py.

201 
202  def run(self, sensorRef):
203  """Subtract an image from a template coadd and measure the result
204 
205  Steps include:
206  - warp template coadd to match WCS of image
207  - PSF match image to warped template
208  - subtract image from PSF-matched, warped template
209  - persist difference image
210  - detect sources
211  - measure sources
212 
213  @param sensorRef: sensor-level butler data reference, used for the following data products:
214  Input only:
215  - calexp
216  - psf
217  - ccdExposureId
218  - ccdExposureId_bits
219  - self.config.coaddName + "Coadd_skyMap"
220  - self.config.coaddName + "Coadd"
221  Input or output, depending on config:
222  - self.config.coaddName + "Diff_subtractedExp"
223  Output, depending on config:
224  - self.config.coaddName + "Diff_matchedExp"
225  - self.config.coaddName + "Diff_src"
226 
227  @return pipe_base Struct containing these fields:
228  - subtractedExposure: exposure after subtracting template;
229  the unpersisted version if subtraction not run but detection run
230  None if neither subtraction nor detection run (i.e. nothing useful done)
231  - subtractRes: results of subtraction task; None if subtraction not run
232  - sources: detected and possibly measured sources; None if detection not run
233  """
234  self.log.info("Processing %s" % (sensorRef.dataId))
235 
236  # initialize outputs and some intermediate products
237  subtractedExposure = None
238  subtractRes = None
239  selectSources = None
240  kernelSources = None
241  controlSources = None
242  diaSources = None
243 
244  # We make one IdFactory that will be used by both icSrc and src datasets;
245  # I don't know if this is the way we ultimately want to do things, but at least
246  # this ensures the source IDs are fully unique.
247  expBits = sensorRef.get("ccdExposureId_bits")
248  expId = long(sensorRef.get("ccdExposureId"))
249  idFactory = afwTable.IdFactory.makeSource(expId, 64 - expBits)
250 
251  # Retrieve the science image we wish to analyze
252  exposure = sensorRef.get("calexp", immediate=True)
253  if self.config.doAddCalexpBackground:
254  mi = exposure.getMaskedImage()
255  mi += sensorRef.get("calexpBackground").getImage()
256  if not exposure.hasPsf():
257  raise pipeBase.TaskError("Exposure has no psf")
258  sciencePsf = exposure.getPsf()
259 
260  subtractedExposureName = self.config.coaddName + "Diff_differenceExp"
261  templateExposure = None # Stitched coadd exposure
262  templateSources = None # Sources on the template image
263  if self.config.doSubtract:
264  templateExposure, templateSources = self.getTemplate(exposure, sensorRef)
265 
266  # compute scienceSigmaOrig: sigma of PSF of science image before pre-convolution
267  ctr = afwGeom.Box2D(exposure.getBBox()).getCenter()
268  psfAttr = PsfAttributes(sciencePsf, afwGeom.Point2I(ctr))
269  scienceSigmaOrig = psfAttr.computeGaussianWidth(psfAttr.ADAPTIVE_MOMENT)
270 
271  # sigma of PSF of template image before warping
272  ctr = afwGeom.Box2D(templateExposure.getBBox()).getCenter()
273  psfAttr = PsfAttributes(templateExposure.getPsf(), afwGeom.Point2I(ctr))
274  templateSigma = psfAttr.computeGaussianWidth(psfAttr.ADAPTIVE_MOMENT)
275 
276  # if requested, convolve the science exposure with its PSF
277  # (properly, this should be a cross-correlation, but our code does not yet support that)
278  # compute scienceSigmaPost: sigma of science exposure with pre-convolution, if done,
279  # else sigma of original science exposure
280  if self.config.doPreConvolve:
281  convControl = afwMath.ConvolutionControl()
282  # cannot convolve in place, so make a new MI to receive convolved image
283  srcMI = exposure.getMaskedImage()
284  destMI = srcMI.Factory(srcMI.getDimensions())
285  srcPsf = sciencePsf
286  if self.config.useGaussianForPreConvolution:
287  # convolve with a simplified PSF model: a double Gaussian
288  kWidth, kHeight = sciencePsf.getLocalKernel().getDimensions()
289  preConvPsf = SingleGaussianPsf(kWidth, kHeight, scienceSigmaOrig)
290  else:
291  # convolve with science exposure's PSF model
292  preConvPsf = srcPsf
293  afwMath.convolve(destMI, srcMI, preConvPsf.getLocalKernel(), convControl)
294  exposure.setMaskedImage(destMI)
295  scienceSigmaPost = scienceSigmaOrig * math.sqrt(2)
296  else:
297  scienceSigmaPost = scienceSigmaOrig
298 
299  # If requested, find sources in the image
300  if self.config.doSelectSources:
301  if not sensorRef.datasetExists("src"):
302  self.log.warn("Src product does not exist; running detection, measurement, selection")
303  # Run own detection and measurement; necessary in nightly processing
304  selectSources = self.subtract.getSelectSources(
305  exposure,
306  sigma = scienceSigmaPost,
307  doSmooth = not self.doPreConvolve,
308  idFactory = idFactory,
309  )
310  else:
311  self.log.info("Source selection via src product")
312  # Sources already exist; for data release processing
313  selectSources = sensorRef.get("src")
314 
315  # Number of basis functions
316  nparam = len(makeKernelBasisList(self.subtract.config.kernel.active,
317  referenceFwhmPix=scienceSigmaPost * FwhmPerSigma,
318  targetFwhmPix=templateSigma * FwhmPerSigma))
319 
320  if self.config.doAddMetrics:
321  # Modify the schema of all Sources
322  kcQa = KernelCandidateQa(nparam)
323  selectSources = kcQa.addToSchema(selectSources)
324 
325  astromRet = self.astrometer.loadAndMatch(exposure=exposure, sourceCat=selectSources)
326  matches = astromRet.matches
327  kernelSources = self.sourceSelector.selectSources(exposure, selectSources, matches=matches)
328  random.shuffle(kernelSources, random.random)
329  controlSources = kernelSources[::self.config.controlStepSize]
330  kernelSources = [k for i,k in enumerate(kernelSources) if i % self.config.controlStepSize]
331 
332  if self.config.doSelectDcrCatalog:
333  redSelector = DiaCatalogSourceSelector(
334  DiaCatalogSourceSelectorConfig(grMin=self.sourceSelector.config.grMax, grMax=99.999))
335  redSources = redSelector.selectSources(exposure, selectSources, matches=matches)
336  controlSources.extend(redSources)
337 
338  blueSelector = DiaCatalogSourceSelector(
339  DiaCatalogSourceSelectorConfig(grMin=-99.999, grMax=self.sourceSelector.config.grMin))
340  blueSources = blueSelector.selectSources(exposure, selectSources, matches=matches)
341  controlSources.extend(blueSources)
342 
343  if self.config.doSelectVariableCatalog:
344  varSelector = DiaCatalogSourceSelector(
345  DiaCatalogSourceSelectorConfig(includeVariable=True))
346  varSources = varSelector.selectSources(exposure, selectSources, matches=matches)
347  controlSources.extend(varSources)
348 
349  self.log.info("Selected %d / %d sources for Psf matching (%d for control sample)"
350  % (len(kernelSources), len(selectSources), len(controlSources)))
351  allresids = {}
352  if self.config.doUseRegister:
353  self.log.info("Registering images")
354 
355  if templateSources is None:
356  # Run detection on the template, which is
357  # temporarily background-subtracted
358  templateSources = self.subtract.getSelectSources(
359  templateExposure,
360  sigma=templateSigma,
361  doSmooth=True,
362  idFactory=idFactory
363  )
364 
365  # Third step: we need to fit the relative astrometry.
366  #
367  wcsResults = self.fitAstrometry(templateSources, templateExposure, selectSources)
368  warpedExp = self.register.warpExposure(templateExposure, wcsResults.wcs,
369  exposure.getWcs(), exposure.getBBox())
370  templateExposure = warpedExp
371 
372  # Create debugging outputs on the astrometric
373  # residuals as a function of position. Persistence
374  # not yet implemented; expected on (I believe) #2636.
375  if self.config.doDebugRegister:
376  # Grab matches to reference catalog
377  srcToMatch = {x.second.getId() : x.first for x in matches}
378 
379  refCoordKey = wcsResults.matches[0].first.getTable().getCoordKey()
380  inCentroidKey = wcsResults.matches[0].second.getTable().getCentroidKey()
381  sids = [m.first.getId() for m in wcsResults.matches]
382  positions = [m.first.get(refCoordKey) for m in wcsResults.matches]
383  residuals = [m.first.get(refCoordKey).getOffsetFrom(wcsResults.wcs.pixelToSky(
384  m.second.get(inCentroidKey))) for m in wcsResults.matches]
385  allresids = dict(zip(sids, zip(positions, residuals)))
386 
387  cresiduals = [m.first.get(refCoordKey).getTangentPlaneOffset(
388  wcsResults.wcs.pixelToSky(
389  m.second.get(inCentroidKey))) for m in wcsResults.matches]
390  colors = numpy.array([-2.5*numpy.log10(srcToMatch[x].get("g"))
391  + 2.5*numpy.log10(srcToMatch[x].get("r"))
392  for x in sids if x in srcToMatch.keys()])
393  dlong = numpy.array([r[0].asArcseconds() for s,r in zip(sids, cresiduals)
394  if s in srcToMatch.keys()])
395  dlat = numpy.array([r[1].asArcseconds() for s,r in zip(sids, cresiduals)
396  if s in srcToMatch.keys()])
397  idx1 = numpy.where(colors<self.sourceSelector.config.grMin)
398  idx2 = numpy.where((colors>=self.sourceSelector.config.grMin)&
399  (colors<=self.sourceSelector.config.grMax))
400  idx3 = numpy.where(colors>self.sourceSelector.config.grMax)
401  rms1Long = IqrToSigma*(numpy.percentile(dlong[idx1],75)-numpy.percentile(dlong[idx1],25))
402  rms1Lat = IqrToSigma*(numpy.percentile(dlat[idx1],75)-numpy.percentile(dlat[idx1],25))
403  rms2Long = IqrToSigma*(numpy.percentile(dlong[idx2],75)-numpy.percentile(dlong[idx2],25))
404  rms2Lat = IqrToSigma*(numpy.percentile(dlat[idx2],75)-numpy.percentile(dlat[idx2],25))
405  rms3Long = IqrToSigma*(numpy.percentile(dlong[idx3],75)-numpy.percentile(dlong[idx3],25))
406  rms3Lat = IqrToSigma*(numpy.percentile(dlat[idx3],75)-numpy.percentile(dlat[idx3],25))
407  self.log.info("Blue star offsets'': %.3f %.3f, %.3f %.3f" % (numpy.median(dlong[idx1]),
408  rms1Long,
409  numpy.median(dlat[idx1]),
410  rms1Lat))
411  self.log.info("Green star offsets'': %.3f %.3f, %.3f %.3f" % (numpy.median(dlong[idx2]),
412  rms2Long,
413  numpy.median(dlat[idx2]),
414  rms2Lat))
415  self.log.info("Red star offsets'': %.3f %.3f, %.3f %.3f" % (numpy.median(dlong[idx3]),
416  rms3Long,
417  numpy.median(dlat[idx3]),
418  rms3Lat))
419 
420  self.metadata.add("RegisterBlueLongOffsetMedian", numpy.median(dlong[idx1]))
421  self.metadata.add("RegisterGreenLongOffsetMedian", numpy.median(dlong[idx2]))
422  self.metadata.add("RegisterRedLongOffsetMedian", numpy.median(dlong[idx3]))
423  self.metadata.add("RegisterBlueLongOffsetStd", rms1Long)
424  self.metadata.add("RegisterGreenLongOffsetStd", rms2Long)
425  self.metadata.add("RegisterRedLongOffsetStd", rms3Long)
426 
427  self.metadata.add("RegisterBlueLatOffsetMedian", numpy.median(dlat[idx1]))
428  self.metadata.add("RegisterGreenLatOffsetMedian", numpy.median(dlat[idx2]))
429  self.metadata.add("RegisterRedLatOffsetMedian", numpy.median(dlat[idx3]))
430  self.metadata.add("RegisterBlueLatOffsetStd", rms1Lat)
431  self.metadata.add("RegisterGreenLatOffsetStd", rms2Lat)
432  self.metadata.add("RegisterRedLatOffsetStd", rms3Lat)
433 
434  # warp template exposure to match exposure,
435  # PSF match template exposure to exposure,
436  # then return the difference
437 
438  #Return warped template... Construct sourceKernelCand list after subtract
439  self.log.info("Subtracting images")
440  subtractRes = self.subtract.subtractExposures(
441  templateExposure=templateExposure,
442  scienceExposure=exposure,
443  scienceFwhmPix=scienceSigmaPost * FwhmPerSigma,
444  templateFwhmPix=templateSigma * FwhmPerSigma,
445  candidateList=kernelSources,
446  convolveTemplate=self.config.convolveTemplate,
447  doWarping=not self.config.doUseRegister
448  )
449  subtractedExposure = subtractRes.subtractedExposure
450 
451  if self.config.doWriteMatchedExp:
452  sensorRef.put(subtractRes.matchedExposure, self.config.coaddName + "Diff_matchedExp")
453 
454  if self.config.doDetection:
455  self.log.info("Running diaSource detection")
456  if subtractedExposure is None:
457  subtractedExposure = sensorRef.get(subtractedExposureName)
458 
459  # Get Psf from the appropriate input image if it doesn't exist
460  if not subtractedExposure.hasPsf():
461  if self.config.convolveTemplate:
462  subtractedExposure.setPsf(exposure.getPsf())
463  else:
464  if templateExposure is None:
465  templateExposure, templateSources = self.getTemplate(exposure, sensorRef)
466  subtractedExposure.setPsf(templateExposure.getPsf())
467 
468  # Erase existing detection mask planes
469  mask = subtractedExposure.getMaskedImage().getMask()
470  mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
471 
472  table = afwTable.SourceTable.make(self.schema, idFactory)
473  table.setMetadata(self.algMetadata)
474  results = self.detection.makeSourceCatalog(
475  table=table,
476  exposure=subtractedExposure,
477  doSmooth=not self.config.doPreConvolve
478  )
479 
480  if self.config.doMerge:
481  fpSet = results.fpSets.positive
482  fpSet.merge(results.fpSets.negative, self.config.growFootprint,
483  self.config.growFootprint, False)
484  diaSources = afwTable.SourceCatalog(table)
485  fpSet.makeSources(diaSources)
486  self.log.info("Merging detections into %d sources" % (len(diaSources)))
487  else:
488  diaSources = results.sources
489 
490  if self.config.doMeasurement:
491  if len(diaSources) < self.config.maxDiaSourcesToMeasure:
492  self.log.info("Running diaSource dipole measurement")
493  self.dipoleMeasurement.run(subtractedExposure, diaSources)
494  else:
495  self.log.info("Running diaSource measurement")
496  self.measurement.run(subtractedExposure, diaSources)
497 
498  # Match with the calexp sources if possible
499  if self.config.doMatchSources:
500  if sensorRef.datasetExists("src"):
501  # Create key,val pair where key=diaSourceId and val=sourceId
502  matchRadAsec = self.config.diaSourceMatchRadius
503  matchRadPixel = matchRadAsec / exposure.getWcs().pixelScale().asArcseconds()
504 
505  srcMatches = afwTable.matchXy(sensorRef.get("src"), diaSources, matchRadPixel, True)
506  srcMatchDict = dict([(srcMatch.second.getId(), srcMatch.first.getId()) for
507  srcMatch in srcMatches])
508  self.log.info("Matched %d / %d diaSources to sources" % (len(srcMatchDict),
509  len(diaSources)))
510  else:
511  self.log.warn("Src product does not exist; cannot match with diaSources")
512  srcMatchDict = {}
513 
514  # Create key,val pair where key=diaSourceId and val=refId
515  refAstromConfig = measAstrom.AstrometryConfig()
516  refAstromConfig.matcher.maxMatchDistArcSec = matchRadAsec
517  refAstrometer = measAstrom.AstrometryTask(refAstromConfig)
518  astromRet = refAstrometer.run(exposure=exposure, sourceCat=diaSources)
519  refMatches = astromRet.matches
520  if refMatches is None:
521  self.log.warn("No diaSource matches with reference catalog")
522  refMatchDict = {}
523  else:
524  self.log.info("Matched %d / %d diaSources to reference catalog" % (len(refMatches),
525  len(diaSources)))
526  refMatchDict = dict([(refMatch.second.getId(), refMatch.first.getId()) for \
527  refMatch in refMatches])
528 
529  # Assign source Ids
530  for diaSource in diaSources:
531  sid = diaSource.getId()
532  if srcMatchDict.has_key(sid):
533  diaSource.set("srcMatchId", srcMatchDict[sid])
534  if refMatchDict.has_key(sid):
535  diaSource.set("refMatchId", refMatchDict[sid])
536 
537  if diaSources is not None and self.config.doWriteSources:
538  sensorRef.put(diaSources, self.config.coaddName + "Diff_diaSrc")
539 
540  if self.config.doAddMetrics and self.config.doSelectSources:
541  self.log.info("Evaluating metrics and control sample")
542 
543  kernelCandList = []
544  for cell in subtractRes.kernelCellSet.getCellList():
545  for cand in cell.begin(False): # include bad candidates
546  kernelCandList.append(cast_KernelCandidateF(cand))
547 
548  # Get basis list to build control sample kernels
549  basisList = afwMath.cast_LinearCombinationKernel(
550  kernelCandList[0].getKernel(KernelCandidateF.ORIG)).getKernelList()
551 
552  controlCandList = \
553  diffimTools.sourceTableToCandidateList(controlSources,
554  subtractRes.warpedExposure, exposure,
555  self.config.subtract.kernel.active,
556  self.config.subtract.kernel.active.detectionConfig,
557  self.log, doBuild=True, basisList=basisList)
558 
559  kcQa.apply(kernelCandList, subtractRes.psfMatchingKernel, subtractRes.backgroundModel,
560  dof=nparam)
561  kcQa.apply(controlCandList, subtractRes.psfMatchingKernel, subtractRes.backgroundModel)
562 
563  if self.config.doDetection:
564  kcQa.aggregate(selectSources, self.metadata, allresids, diaSources)
565  else:
566  kcQa.aggregate(selectSources, self.metadata, allresids)
567 
568  sensorRef.put(selectSources, self.config.coaddName + "Diff_kernelSrc")
569 
570  if self.config.doWriteSubtractedExp:
571  sensorRef.put(subtractedExposure, subtractedExposureName)
572 
573  self.runDebug(exposure, subtractRes, selectSources, kernelSources, diaSources)
574  return pipeBase.Struct(
575  subtractedExposure=subtractedExposure,
576  subtractRes=subtractRes,
577  sources=diaSources,
578  )
lsst::afw::table::matchXy
SourceMatchVector matchXy(SourceCatalog const &cat, double radius, bool symmetric=true)
lsst::afw::math::ConvolutionControl
Parameters to control convolution.
Definition: ConvolveImage.h:58
lsst::afw::table::SortedCatalogT
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
lsst::afw::math::convolve
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, bool doNormalize, bool doCopyEdge=false)
Old, deprecated version of convolve.
Definition: ConvolveImage.cc:347
lsst::afw::geom::Box2D
A floating-point coordinate rectangle geometry.
Definition: Box.h:271
def lsst.pipe.tasks.imageDifference.ImageDifferenceTask.runDebug (   self,
  exposure,
  subtractRes,
  selectSources,
  kernelSources,
  diaSources 
)

Definition at line 590 of file imageDifference.py.

590 
591  def runDebug(self, exposure, subtractRes, selectSources, kernelSources, diaSources):
592  import lsstDebug
593  import lsst.afw.display.ds9 as ds9
594  display = lsstDebug.Info(__name__).display
595  showSubtracted = lsstDebug.Info(__name__).showSubtracted
596  showPixelResiduals = lsstDebug.Info(__name__).showPixelResiduals
597  showDiaSources = lsstDebug.Info(__name__).showDiaSources
598  showDipoles = lsstDebug.Info(__name__).showDipoles
599  maskTransparency = lsstDebug.Info(__name__).maskTransparency
600  if not maskTransparency:
601  maskTransparency = 0
602  ds9.setMaskTransparency(maskTransparency)
603 
604  if display and showSubtracted:
605  ds9.mtv(subtractRes.subtractedExposure, frame=lsstDebug.frame, title="Subtracted image")
606  mi = subtractRes.subtractedExposure.getMaskedImage()
607  x0, y0 = mi.getX0(), mi.getY0()
608  with ds9.Buffering():
609  for s in diaSources:
610  x, y = s.getX() - x0, s.getY() - y0
611  ctype = "red" if s.get("flags.negative") else "yellow"
612  if (s.get("flags.pixel.interpolated.center") or s.get("flags.pixel.saturated.center") or
613  s.get("flags.pixel.cr.center")):
614  ptype = "x"
615  elif (s.get("flags.pixel.interpolated.any") or s.get("flags.pixel.saturated.any") or
616  s.get("flags.pixel.cr.any")):
617  ptype = "+"
618  else:
619  ptype = "o"
620  ds9.dot(ptype, x, y, size=4, frame=lsstDebug.frame, ctype=ctype)
621  lsstDebug.frame += 1
622 
623  if display and showPixelResiduals and selectSources:
624  import lsst.ip.diffim.utils as diUtils
625  nonKernelSources = []
626  for source in selectSources:
627  if not source in kernelSources:
628  nonKernelSources.append(source)
629 
630  diUtils.plotPixelResiduals(exposure,
631  subtractRes.warpedExposure,
632  subtractRes.subtractedExposure,
633  subtractRes.kernelCellSet,
634  subtractRes.psfMatchingKernel,
635  subtractRes.backgroundModel,
636  nonKernelSources,
637  self.subtract.config.kernel.active.detectionConfig,
638  origVariance=False)
639  diUtils.plotPixelResiduals(exposure,
640  subtractRes.warpedExposure,
641  subtractRes.subtractedExposure,
642  subtractRes.kernelCellSet,
643  subtractRes.psfMatchingKernel,
644  subtractRes.backgroundModel,
645  nonKernelSources,
646  self.subtract.config.kernel.active.detectionConfig,
647  origVariance=True)
648  if display and showDiaSources:
649  import lsst.ip.diffim.utils as diUtils
650  flagChecker = SourceFlagChecker(diaSources)
651  isFlagged = [flagChecker(x) for x in diaSources]
652  isDipole = [x.get("classification.dipole") for x in diaSources]
653  diUtils.showDiaSources(diaSources, subtractRes.subtractedExposure, isFlagged, isDipole,
654  frame=lsstDebug.frame)
655  lsstDebug.frame += 1
656 
657  if display and showDipoles:
658  DipoleAnalysis().displayDipoles(subtractRes.subtractedExposure, diaSources,
659  frame=lsstDebug.frame)
660  lsstDebug.frame += 1

Member Data Documentation

string lsst.pipe.tasks.imageDifference.ImageDifferenceTask._DefaultName = "imageDifference"
staticprivate

Definition at line 173 of file imageDifference.py.

lsst.pipe.tasks.imageDifference.ImageDifferenceTask.algMetadata

Definition at line 186 of file imageDifference.py.

lsst.pipe.tasks.imageDifference.ImageDifferenceTask.astrometer

Definition at line 184 of file imageDifference.py.

lsst.pipe.tasks.imageDifference.ImageDifferenceTask.ConfigClass = ImageDifferenceConfig
static

Definition at line 172 of file imageDifference.py.

lsst.pipe.tasks.imageDifference.ImageDifferenceTask.schema

Definition at line 185 of file imageDifference.py.

lsst.pipe.tasks.imageDifference.ImageDifferenceTask.sourceSelector

Definition at line 183 of file imageDifference.py.

The documentation for this class was generated from the following file:
Generated on Thu Sep 24 2015 02:29:42 for LSSTApplications by   doxygen 1.8.5