Inheritance diagram for lsst.pipe.tasks.imageDifference.ImageDifferenceFromTemplateConfig:
Detailed Description
doAddCalexpBackground = pexConfig.Field(dtype=bool, default=False,
doc="Add background to calexp before processing it. "
"Useful as ipDiffim does background matching.")
doUseRegister = pexConfig.Field(dtype=bool, default=False,
doc="Re-compute astrometry on the template. "
"Use image-to-image registration to align template with "
"science image (AL only).")
doDebugRegister = pexConfig.Field(dtype=bool, default=False,
doc="Writing debugging data for doUseRegister")
doSelectSources = pexConfig.Field(dtype=bool, default=False,
doc="Select stars to use for kernel fitting (AL only)")
doSelectDcrCatalog = pexConfig.Field(dtype=bool, default=False,
doc="Select stars of extreme color as part "
"of the control sample (AL only)")
doSelectVariableCatalog = pexConfig.Field(dtype=bool, default=False,
doc="Select stars that are variable to be part "
"of the control sample (AL only)")
doSubtract = pexConfig.Field(dtype=bool, default=True, doc="Compute subtracted exposure?")
doPreConvolve = pexConfig.Field(dtype=bool, default=False,
doc="Not in use. Superseded by useScoreImageDetection.",
deprecated="This option superseded by useScoreImageDetection."
" Will be removed after v22.")
useScoreImageDetection = pexConfig.Field(
dtype=bool, default=False, doc="Calculate the pre-convolved AL likelihood or "
"the Zogy score image. Use it for source detection (if doDetection=True).")
doWriteScoreExp = pexConfig.Field(
dtype=bool, default=False, doc="Write AL likelihood or Zogy score exposure?")
doScaleTemplateVariance = pexConfig.Field(dtype=bool, default=False,
doc="Scale variance of the template before PSF matching")
doScaleDiffimVariance = pexConfig.Field(dtype=bool, default=True,
doc="Scale variance of the diffim before PSF matching. "
"You may do either this or template variance scaling, "
"or neither. (Doing both is a waste of CPU.)")
useGaussianForPreConvolution = pexConfig.Field(
dtype=bool, default=False, doc="Use a simple gaussian PSF model for pre-convolution "
"(oherwise use exposure PSF)? (AL and if useScoreImageDetection=True only)")
doDetection = pexConfig.Field(dtype=bool, default=True, doc="Detect sources?")
doDecorrelation = pexConfig.Field(dtype=bool, default=True,
doc="Perform diffim decorrelation to undo pixel correlation due to A&L "
"kernel convolution (AL only)? If True, also update the diffim PSF.")
doMerge = pexConfig.Field(dtype=bool, default=True,
doc="Merge positive and negative diaSources with grow radius "
"set by growFootprint")
doMatchSources = pexConfig.Field(dtype=bool, default=False,
doc="Match diaSources with input calexp sources and ref catalog sources")
doMeasurement = pexConfig.Field(dtype=bool, default=True, doc="Measure diaSources?")
doDipoleFitting = pexConfig.Field(dtype=bool, default=True, doc="Measure dipoles using new algorithm?")
doForcedMeasurement = pexConfig.Field(
dtype=bool,
default=True,
doc="Force photometer diaSource locations on PVI?")
doWriteSubtractedExp = pexConfig.Field(
dtype=bool, default=True, doc="Write difference exposure (AL and Zogy) ?")
doWriteWarpedExp = pexConfig.Field(
dtype=bool, default=False, doc="Write WCS, warped template coadd exposure?")
doWriteMatchedExp = pexConfig.Field(dtype=bool, default=False,
doc="Write warped and PSF-matched template coadd exposure?")
doWriteSources = pexConfig.Field(dtype=bool, default=True, doc="Write sources?")
doAddMetrics = pexConfig.Field(dtype=bool, default=False,
doc="Add columns to the source table to hold analysis metrics?")
coaddName = pexConfig.Field(
doc="coadd name: typically one of deep, goodSeeing, or dcr",
dtype=str,
default="deep",
)
convolveTemplate = pexConfig.Field(
doc="Which image gets convolved (default = template)",
dtype=bool,
default=True
)
refObjLoader = pexConfig.ConfigurableField(
target=LoadIndexedReferenceObjectsTask,
doc="reference object loader",
)
astrometer = pexConfig.ConfigurableField(
target=AstrometryTask,
doc="astrometry task; used to match sources to reference objects, but not to fit a WCS",
)
sourceSelector = pexConfig.ConfigurableField(
target=ObjectSizeStarSelectorTask,
doc="Source selection algorithm",
)
subtract = subtractAlgorithmRegistry.makeField("Subtraction Algorithm", default="al")
decorrelate = pexConfig.ConfigurableField(
target=DecorrelateALKernelSpatialTask,
doc="Decorrelate effects of A&L kernel convolution on image difference, only if doSubtract is True. "
"If this option is enabled, then detection.thresholdValue should be set to 5.0 (rather than the "
"default of 5.5).",
)
# Old style ImageMapper grid. ZogyTask has its own grid option
doSpatiallyVarying = pexConfig.Field(
dtype=bool,
default=False,
doc="Perform A&L decorrelation on a grid across the "
"image in order to allow for spatial variations. Zogy does not use this option."
)
detection = pexConfig.ConfigurableField(
target=SourceDetectionTask,
doc="Low-threshold detection for final measurement",
)
measurement = pexConfig.ConfigurableField(
target=DipoleFitTask,
doc="Enable updated dipole fitting method",
)
doApCorr = lsst.pex.config.Field(
dtype=bool,
default=True,
doc="Run subtask to apply aperture corrections"
)
applyApCorr = lsst.pex.config.ConfigurableField(
target=ApplyApCorrTask,
doc="Subtask to apply aperture corrections"
)
forcedMeasurement = pexConfig.ConfigurableField(
target=ForcedMeasurementTask,
doc="Subtask to force photometer PVI at diaSource location.",
)
getTemplate = pexConfig.ConfigurableField(
target=GetCoaddAsTemplateTask,
doc="Subtask to retrieve template exposure and sources",
)
scaleVariance = pexConfig.ConfigurableField(
target=ScaleVarianceTask,
doc="Subtask to rescale the variance of the template "
"to the statistically expected level"
)
controlStepSize = pexConfig.Field(
doc="What step size (every Nth one) to select a control sample from the kernelSources",
dtype=int,
default=5
)
controlRandomSeed = pexConfig.Field(
doc="Random seed for shuffing the control sample",
dtype=int,
default=10
)
register = pexConfig.ConfigurableField(
target=RegisterTask,
doc="Task to enable image-to-image image registration (warping)",
)
kernelSourcesFromRef = pexConfig.Field(
doc="Select sources to measure kernel from reference catalog if True, template if false",
dtype=bool,
default=False
)
templateSipOrder = pexConfig.Field(
dtype=int, default=2,
doc="Sip Order for fitting the Template Wcs (default is too high, overfitting)"
)
growFootprint = pexConfig.Field(
dtype=int, default=2,
doc="Grow positive and negative footprints by this amount before merging"
)
diaSourceMatchRadius = pexConfig.Field(
dtype=float, default=0.5,
doc="Match radius (in arcseconds) for DiaSource to Source association"
)
requiredTemplateFraction = pexConfig.Field(
dtype=float, default=0.1,
doc="Do not attempt to run task if template covers less than this fraction of pixels."
"Setting to 0 will always attempt image subtraction"
)
doSkySources = pexConfig.Field(
dtype=bool,
default=False,
doc="Generate sky sources?",
)
skySources = pexConfig.ConfigurableField(
target=SkyObjectsTask,
doc="Generate sky sources",
)
def setDefaults(self):
# defaults are OK for catalog and diacatalog
self.subtract['al'].kernel.name = "AL"
self.subtract['al'].kernel.active.fitForBackground = True
self.subtract['al'].kernel.active.spatialKernelOrder = 1
self.subtract['al'].kernel.active.spatialBgOrder = 2
# DiaSource Detection
self.detection.thresholdPolarity = "both"
self.detection.thresholdValue = 5.0
self.detection.reEstimateBackground = False
self.detection.thresholdType = "pixel_stdev"
# Add filtered flux measurement, the correct measurement for pre-convolved images.
# Enable all measurements, regardless of doPreConvolve, as it makes data harvesting easier.
# To change that you must modify algorithms.names in the task's applyOverrides method,
# after the user has set doPreConvolve.
self.measurement.algorithms.names.add('base_PeakLikelihoodFlux')
self.measurement.plugins.names |= ['ext_trailedSources_Naive',
'base_LocalPhotoCalib',
'base_LocalWcs']
self.forcedMeasurement.plugins = ["base_TransformedCentroid", "base_PsfFlux"]
self.forcedMeasurement.copyColumns = {
"id": "objectId", "parent": "parentObjectId", "coord_ra": "coord_ra", "coord_dec": "coord_dec"}
self.forcedMeasurement.slots.centroid = "base_TransformedCentroid"
self.forcedMeasurement.slots.shape = None
# For shuffling the control sample
random.seed(self.controlRandomSeed)
def validate(self):
pexConfig.Config.validate(self)
if not self.doSubtract and not self.doDetection:
raise ValueError("Either doSubtract or doDetection must be enabled.")
if self.doMeasurement and not self.doDetection:
raise ValueError("Cannot run source measurement without source detection.")
if self.doMerge and not self.doDetection:
raise ValueError("Cannot run source merging without source detection.")
if self.doSkySources and not self.doDetection:
raise ValueError("Cannot run sky source creation without source detection.")
if self.doUseRegister and not self.doSelectSources:
raise ValueError("doUseRegister=True and doSelectSources=False. "
"Cannot run RegisterTask without selecting sources.")
if hasattr(self.getTemplate, "coaddName"):
if self.getTemplate.coaddName != self.coaddName:
raise ValueError("Mis-matched coaddName and getTemplate.coaddName in the config.")
if self.doScaleDiffimVariance and self.doScaleTemplateVariance:
raise ValueError("Scaling the diffim variance and scaling the template variance "
"are both set. Please choose one or the other.")
# We cannot allow inconsistencies that would lead to None or not available output products
if self.subtract.name == 'zogy':
if self.doWriteMatchedExp:
raise ValueError("doWriteMatchedExp=True Matched exposure is not "
"calculated in zogy subtraction.")
if self.doAddMetrics:
raise ValueError("doAddMetrics=True Kernel metrics does not exist in zogy subtraction.")
if self.doDecorrelation:
raise ValueError(
"doDecorrelation=True The decorrelation afterburner does not exist in zogy subtraction.")
if self.doSelectSources:
raise ValueError(
"doSelectSources=True Selecting sources for PSF matching is not a zogy option.")
if self.useGaussianForPreConvolution:
raise ValueError(
"useGaussianForPreConvolution=True This is an AL subtraction only option.")
else:
# AL only consistency checks
if self.useScoreImageDetection and not self.convolveTemplate:
raise ValueError(
"convolveTemplate=False and useScoreImageDetection=True "
"Pre-convolution and matching of the science image is not a supported operation.")
if self.doWriteSubtractedExp and self.useScoreImageDetection:
raise ValueError(
"doWriteSubtractedExp=True and useScoreImageDetection=True "
"Regular difference image is not calculated. "
"AL subtraction calculates either the regular difference image or the score image.")
if self.doWriteScoreExp and not self.useScoreImageDetection:
raise ValueError(
"doWriteScoreExp=True and useScoreImageDetection=False "
"Score image is not calculated. "
"AL subtraction calculates either the regular difference image or the score image.")
if self.doAddMetrics and not self.doSubtract:
raise ValueError("Subtraction must be enabled for kernel metrics calculation.")
if self.useGaussianForPreConvolution and not self.useScoreImageDetection:
raise ValueError(
"useGaussianForPreConvolution=True and useScoreImageDetection=False "
"Gaussian PSF approximation exists only for AL subtraction w/ pre-convolution.")
class ImageDifferenceTaskRunner(pipeBase.ButlerInitializedTaskRunner):
@staticmethod
def getTargetList(parsedCmd, **kwargs):
return pipeBase.TaskRunner.getTargetList(parsedCmd, templateIdList=parsedCmd.templateId.idList,
**kwargs)
class ImageDifferenceTask(pipeBase.CmdLineTask, pipeBase.PipelineTask):
ConfigClass = ImageDifferenceConfig RunnerClass = ImageDifferenceTaskRunner _DefaultName = "imageDifference" def __init__(self, butler=None, **kwargs):
super().__init__(**kwargs)
self.makeSubtask("getTemplate")
self.makeSubtask("subtract")
if self.config.subtract.name == 'al' and self.config.doDecorrelation:
self.makeSubtask("decorrelate")
if self.config.doScaleTemplateVariance or self.config.doScaleDiffimVariance:
self.makeSubtask("scaleVariance")
if self.config.doUseRegister:
self.makeSubtask("register")
self.schema = afwTable.SourceTable.makeMinimalSchema()
if self.config.doSelectSources:
self.makeSubtask("sourceSelector")
if self.config.kernelSourcesFromRef:
self.makeSubtask('refObjLoader', butler=butler)
self.makeSubtask("astrometer", refObjLoader=self.refObjLoader)
self.algMetadata = dafBase.PropertyList()
if self.config.doDetection:
self.makeSubtask("detection", schema=self.schema)
if self.config.doMeasurement:
self.makeSubtask("measurement", schema=self.schema,
algMetadata=self.algMetadata)
if self.config.doApCorr:
self.makeSubtask("applyApCorr", schema=self.measurement.schema)
if self.config.doForcedMeasurement:
self.schema.addField(
"ip_diffim_forced_PsfFlux_instFlux", "D",
"Forced PSF flux measured on the direct image.",
units="count")
self.schema.addField(
"ip_diffim_forced_PsfFlux_instFluxErr", "D",
"Forced PSF flux error measured on the direct image.",
units="count")
self.schema.addField(
"ip_diffim_forced_PsfFlux_area", "F",
"Forced PSF flux effective area of PSF.",
units="pixel")
self.schema.addField(
"ip_diffim_forced_PsfFlux_flag", "Flag",
"Forced PSF flux general failure flag.")
self.schema.addField(
"ip_diffim_forced_PsfFlux_flag_noGoodPixels", "Flag",
"Forced PSF flux not enough non-rejected pixels in data to attempt the fit.")
self.schema.addField(
"ip_diffim_forced_PsfFlux_flag_edge", "Flag",
"Forced PSF flux object was too close to the edge of the image to use the full PSF model.")
self.makeSubtask("forcedMeasurement", refSchema=self.schema)
if self.config.doMatchSources:
self.schema.addField("refMatchId", "L", "unique id of reference catalog match")
self.schema.addField("srcMatchId", "L", "unique id of source match")
if self.config.doSkySources:
self.makeSubtask("skySources")
self.skySourceKey = self.schema.addField("sky_source", type="Flag", doc="Sky objects.")
# initialize InitOutputs
self.outputSchema = afwTable.SourceCatalog(self.schema)
self.outputSchema.getTable().setMetadata(self.algMetadata)
@staticmethod
def makeIdFactory(expId, expBits):
return ExposureIdInfo(expId, expBits).makeSourceIdFactory()
@lsst.utils.inheritDoc(pipeBase.PipelineTask)
def runQuantum(self, butlerQC: pipeBase.ButlerQuantumContext,
inputRefs: pipeBase.InputQuantizedConnection,
outputRefs: pipeBase.OutputQuantizedConnection):
inputs = butlerQC.get(inputRefs)
self.log.info("Processing %s", butlerQC.quantum.dataId)
expId, expBits = butlerQC.quantum.dataId.pack("visit_detector",
returnMaxBits=True)
idFactory = self.makeIdFactory(expId=expId, expBits=expBits)
if self.config.coaddName == 'dcr':
templateExposures = inputRefs.dcrCoadds
else:
templateExposures = inputRefs.coaddExposures
templateStruct = self.getTemplate.runQuantum(
inputs['exposure'], butlerQC, inputRefs.skyMap, templateExposures
)
self.checkTemplateIsSufficient(templateStruct.exposure)
outputs = self.run(exposure=inputs['exposure'],
templateExposure=templateStruct.exposure,
idFactory=idFactory)
# Consistency with runDataref gen2 handling
if outputs.diaSources is None:
del outputs.diaSources
butlerQC.put(outputs, outputRefs)
@timeMethod
def runDataRef(self, sensorRef, templateIdList=None):
subtractedExposureName = self.config.coaddName + "Diff_differenceExp"
subtractedExposure = None
selectSources = None
calexpBackgroundExposure = None
self.log.info("Processing %s", sensorRef.dataId)
# We make one IdFactory that will be used by both icSrc and src datasets;
# I don't know if this is the way we ultimately want to do things, but at least
# this ensures the source IDs are fully unique.
idFactory = self.makeIdFactory(expId=int(sensorRef.get("ccdExposureId")),
expBits=sensorRef.get("ccdExposureId_bits"))
if self.config.doAddCalexpBackground:
calexpBackgroundExposure = sensorRef.get("calexpBackground")
# Retrieve the science image we wish to analyze
exposure = sensorRef.get("calexp", immediate=True)
# Retrieve the template image
template = self.getTemplate.runDataRef(exposure, sensorRef, templateIdList=templateIdList)
if sensorRef.datasetExists("src"):
self.log.info("Source selection via src product")
# Sources already exist; for data release processing
selectSources = sensorRef.get("src")
if not self.config.doSubtract and self.config.doDetection:
# If we don't do subtraction, we need the subtracted exposure from the repo
subtractedExposure = sensorRef.get(subtractedExposureName)
# Both doSubtract and doDetection cannot be False
results = self.run(exposure=exposure,
selectSources=selectSources,
templateExposure=template.exposure,
templateSources=template.sources,
idFactory=idFactory,
calexpBackgroundExposure=calexpBackgroundExposure,
subtractedExposure=subtractedExposure)
if self.config.doWriteSources and results.diaSources is not None:
sensorRef.put(results.diaSources, self.config.coaddName + "Diff_diaSrc")
if self.config.doWriteWarpedExp:
sensorRef.put(results.warpedExposure, self.config.coaddName + "Diff_warpedExp")
if self.config.doWriteMatchedExp:
sensorRef.put(results.matchedExposure, self.config.coaddName + "Diff_matchedExp")
if self.config.doAddMetrics and self.config.doSelectSources:
sensorRef.put(results.selectSources, self.config.coaddName + "Diff_kernelSrc")
if self.config.doWriteSubtractedExp:
sensorRef.put(results.subtractedExposure, subtractedExposureName)
if self.config.doWriteScoreExp:
sensorRef.put(results.scoreExposure, self.config.coaddName + "Diff_scoreExp")
return results
@timeMethod
def run(self, exposure=None, selectSources=None, templateExposure=None, templateSources=None,
idFactory=None, calexpBackgroundExposure=None, subtractedExposure=None):
subtractRes = None
controlSources = None
subtractedExposure = None
scoreExposure = None
diaSources = None
kernelSources = None
# We'll clone exposure if modified but will still need the original
exposureOrig = exposure
if self.config.doAddCalexpBackground:
mi = exposure.getMaskedImage()
mi += calexpBackgroundExposure.getImage()
if not exposure.hasPsf():
raise pipeBase.TaskError("Exposure has no psf")
sciencePsf = exposure.getPsf()
if self.config.doSubtract:
if self.config.doScaleTemplateVariance:
self.log.info("Rescaling template variance")
templateVarFactor = self.scaleVariance.run(
templateExposure.getMaskedImage())
self.log.info("Template variance scaling factor: %.2f", templateVarFactor)
self.metadata.add("scaleTemplateVarianceFactor", templateVarFactor)
self.metadata.add("psfMatchingAlgorithm", self.config.subtract.name)
if self.config.subtract.name == 'zogy':
subtractRes = self.subtract.run(exposure, templateExposure, doWarping=True)
scoreExposure = subtractRes.scoreExp
subtractedExposure = subtractRes.diffExp
subtractRes.subtractedExposure = subtractedExposure
subtractRes.matchedExposure = None
elif self.config.subtract.name == 'al':
# compute scienceSigmaOrig: sigma of PSF of science image before pre-convolution
# Just need a rough estimate; average positions are fine
sciAvgPos = sciencePsf.getAveragePosition()
scienceSigmaOrig = sciencePsf.computeShape(sciAvgPos).getDeterminantRadius()
templatePsf = templateExposure.getPsf()
templateAvgPos = templatePsf.getAveragePosition()
templateSigma = templatePsf.computeShape(templateAvgPos).getDeterminantRadius()
# if requested, convolve the science exposure with its PSF
# (properly, this should be a cross-correlation, but our code does not yet support that)
# compute scienceSigmaPost: sigma of science exposure with pre-convolution, if done,
# else sigma of original science exposure
# TODO: DM-22762 This functional block should be moved into its own method
preConvPsf = None
if self.config.useScoreImageDetection:
self.log.warning("AL likelihood image: pre-convolution of PSF is not implemented.")
convControl = afwMath.ConvolutionControl()
# cannot convolve in place, so need a new image anyway
srcMI = exposure.maskedImage
exposure = exposure.clone() # New deep copy
srcPsf = sciencePsf
if self.config.useGaussianForPreConvolution:
self.log.info(
"AL likelihood image: Using Gaussian (sigma=%.2f) PSF estimation "
"for science image pre-convolution", scienceSigmaOrig)
# convolve with a simplified PSF model: a double Gaussian
kWidth, kHeight = sciencePsf.getLocalKernel().getDimensions()
preConvPsf = SingleGaussianPsf(kWidth, kHeight, scienceSigmaOrig)
else:
# convolve with science exposure's PSF model
self.log.info(
"AL likelihood image: Using the science image PSF for pre-convolution.")
preConvPsf = srcPsf
afwMath.convolve(exposure.maskedImage, srcMI, preConvPsf.getLocalKernel(), convControl)
scienceSigmaPost = scienceSigmaOrig*math.sqrt(2)
else:
scienceSigmaPost = scienceSigmaOrig
# If requested, find and select sources from the image
# else, AL subtraction will do its own source detection
# TODO: DM-22762 This functional block should be moved into its own method
if self.config.doSelectSources:
if selectSources is None:
self.log.warning("Src product does not exist; running detection, measurement,"
" selection")
# Run own detection and measurement; necessary in nightly processing
selectSources = self.subtract.getSelectSources(
exposure,
sigma=scienceSigmaPost,
doSmooth=not self.config.useScoreImageDetection,
idFactory=idFactory,
)
if self.config.doAddMetrics:
# Number of basis functions
nparam = len(makeKernelBasisList(self.subtract.config.kernel.active,
referenceFwhmPix=scienceSigmaPost*FwhmPerSigma,
targetFwhmPix=templateSigma*FwhmPerSigma))
# Modify the schema of all Sources
# DEPRECATED: This is a data dependent (nparam) output product schema
# outside the task constructor.
# NOTE: The pre-determination of nparam at this point
# may be incorrect as the template psf is warped later in
# ImagePsfMatchTask.matchExposures()
kcQa = KernelCandidateQa(nparam)
selectSources = kcQa.addToSchema(selectSources)
if self.config.kernelSourcesFromRef:
# match exposure sources to reference catalog
astromRet = self.astrometer.loadAndMatch(exposure=exposure, sourceCat=selectSources)
matches = astromRet.matches
elif templateSources:
# match exposure sources to template sources
mc = afwTable.MatchControl()
mc.findOnlyClosest = False
matches = afwTable.matchRaDec(templateSources, selectSources, 1.0*geom.arcseconds,
mc)
else:
raise RuntimeError("doSelectSources=True and kernelSourcesFromRef=False,"
"but template sources not available. Cannot match science "
"sources with template sources. Run process* on data from "
"which templates are built.")
kernelSources = self.sourceSelector.run(selectSources, exposure=exposure,
matches=matches).sourceCat
random.shuffle(kernelSources, random.random)
controlSources = kernelSources[::self.config.controlStepSize]
kernelSources = [k for i, k in enumerate(kernelSources)
if i % self.config.controlStepSize]
if self.config.doSelectDcrCatalog:
redSelector = DiaCatalogSourceSelectorTask(
DiaCatalogSourceSelectorConfig(grMin=self.sourceSelector.config.grMax,
grMax=99.999))
redSources = redSelector.selectStars(exposure, selectSources, matches=matches).starCat
controlSources.extend(redSources)
blueSelector = DiaCatalogSourceSelectorTask(
DiaCatalogSourceSelectorConfig(grMin=-99.999,
grMax=self.sourceSelector.config.grMin))
blueSources = blueSelector.selectStars(exposure, selectSources,
matches=matches).starCat
controlSources.extend(blueSources)
if self.config.doSelectVariableCatalog:
varSelector = DiaCatalogSourceSelectorTask(
DiaCatalogSourceSelectorConfig(includeVariable=True))
varSources = varSelector.selectStars(exposure, selectSources, matches=matches).starCat
controlSources.extend(varSources)
self.log.info("Selected %d / %d sources for Psf matching (%d for control sample)",
len(kernelSources), len(selectSources), len(controlSources))
allresids = {}
# TODO: DM-22762 This functional block should be moved into its own method
if self.config.doUseRegister:
self.log.info("Registering images")
if templateSources is None:
# Run detection on the template, which is
# temporarily background-subtracted
# sigma of PSF of template image before warping
templateSources = self.subtract.getSelectSources(
templateExposure,
sigma=templateSigma,
doSmooth=True,
idFactory=idFactory
)
# Third step: we need to fit the relative astrometry.
#
wcsResults = self.fitAstrometry(templateSources, templateExposure, selectSources)
warpedExp = self.register.warpExposure(templateExposure, wcsResults.wcs,
exposure.getWcs(), exposure.getBBox())
templateExposure = warpedExp
# Create debugging outputs on the astrometric
# residuals as a function of position. Persistence
# not yet implemented; expected on (I believe) #2636.
if self.config.doDebugRegister:
# Grab matches to reference catalog
srcToMatch = {x.second.getId(): x.first for x in matches}
refCoordKey = wcsResults.matches[0].first.getTable().getCoordKey()
inCentroidKey = wcsResults.matches[0].second.getTable().getCentroidSlot().getMeasKey()
sids = [m.first.getId() for m in wcsResults.matches]
positions = [m.first.get(refCoordKey) for m in wcsResults.matches]
residuals = [m.first.get(refCoordKey).getOffsetFrom(wcsResults.wcs.pixelToSky(
m.second.get(inCentroidKey))) for m in wcsResults.matches]
allresids = dict(zip(sids, zip(positions, residuals)))
cresiduals = [m.first.get(refCoordKey).getTangentPlaneOffset(
wcsResults.wcs.pixelToSky(
m.second.get(inCentroidKey))) for m in wcsResults.matches]
colors = numpy.array([-2.5*numpy.log10(srcToMatch[x].get("g"))
+ 2.5*numpy.log10(srcToMatch[x].get("r"))
for x in sids if x in srcToMatch.keys()])
dlong = numpy.array([r[0].asArcseconds() for s, r in zip(sids, cresiduals)
if s in srcToMatch.keys()])
dlat = numpy.array([r[1].asArcseconds() for s, r in zip(sids, cresiduals)
if s in srcToMatch.keys()])
idx1 = numpy.where(colors < self.sourceSelector.config.grMin)
idx2 = numpy.where((colors >= self.sourceSelector.config.grMin)
& (colors <= self.sourceSelector.config.grMax))
idx3 = numpy.where(colors > self.sourceSelector.config.grMax)
rms1Long = IqrToSigma*(
(numpy.percentile(dlong[idx1], 75) - numpy.percentile(dlong[idx1], 25)))
rms1Lat = IqrToSigma*(numpy.percentile(dlat[idx1], 75)
- numpy.percentile(dlat[idx1], 25))
rms2Long = IqrToSigma*(
(numpy.percentile(dlong[idx2], 75) - numpy.percentile(dlong[idx2], 25)))
rms2Lat = IqrToSigma*(numpy.percentile(dlat[idx2], 75)
- numpy.percentile(dlat[idx2], 25))
rms3Long = IqrToSigma*(
(numpy.percentile(dlong[idx3], 75) - numpy.percentile(dlong[idx3], 25)))
rms3Lat = IqrToSigma*(numpy.percentile(dlat[idx3], 75)
- numpy.percentile(dlat[idx3], 25))
self.log.info("Blue star offsets'': %.3f %.3f, %.3f %.3f",
numpy.median(dlong[idx1]), rms1Long,
numpy.median(dlat[idx1]), rms1Lat)
self.log.info("Green star offsets'': %.3f %.3f, %.3f %.3f",
numpy.median(dlong[idx2]), rms2Long,
numpy.median(dlat[idx2]), rms2Lat)
self.log.info("Red star offsets'': %.3f %.3f, %.3f %.3f",
numpy.median(dlong[idx3]), rms3Long,
numpy.median(dlat[idx3]), rms3Lat)
self.metadata.add("RegisterBlueLongOffsetMedian", numpy.median(dlong[idx1]))
self.metadata.add("RegisterGreenLongOffsetMedian", numpy.median(dlong[idx2]))
self.metadata.add("RegisterRedLongOffsetMedian", numpy.median(dlong[idx3]))
self.metadata.add("RegisterBlueLongOffsetStd", rms1Long)
self.metadata.add("RegisterGreenLongOffsetStd", rms2Long)
self.metadata.add("RegisterRedLongOffsetStd", rms3Long)
self.metadata.add("RegisterBlueLatOffsetMedian", numpy.median(dlat[idx1]))
self.metadata.add("RegisterGreenLatOffsetMedian", numpy.median(dlat[idx2]))
self.metadata.add("RegisterRedLatOffsetMedian", numpy.median(dlat[idx3]))
self.metadata.add("RegisterBlueLatOffsetStd", rms1Lat)
self.metadata.add("RegisterGreenLatOffsetStd", rms2Lat)
self.metadata.add("RegisterRedLatOffsetStd", rms3Lat)
# warp template exposure to match exposure,
# PSF match template exposure to exposure,
# then return the difference
# Return warped template... Construct sourceKernelCand list after subtract
self.log.info("Subtracting images")
subtractRes = self.subtract.subtractExposures(
templateExposure=templateExposure,
scienceExposure=exposure,
candidateList=kernelSources,
convolveTemplate=self.config.convolveTemplate,
doWarping=not self.config.doUseRegister
)
if self.config.useScoreImageDetection:
scoreExposure = subtractRes.subtractedExposure
else:
subtractedExposure = subtractRes.subtractedExposure
if self.config.doDetection:
self.log.info("Computing diffim PSF")
# Get Psf from the appropriate input image if it doesn't exist
if subtractedExposure is not None and not subtractedExposure.hasPsf():
if self.config.convolveTemplate:
subtractedExposure.setPsf(exposure.getPsf())
else:
subtractedExposure.setPsf(templateExposure.getPsf())
# If doSubtract is False, then subtractedExposure was fetched from disk (above),
# thus it may have already been decorrelated. Thus, we do not decorrelate if
# doSubtract is False.
# NOTE: At this point doSubtract == True
if self.config.doDecorrelation and self.config.doSubtract:
preConvKernel = None
if self.config.useGaussianForPreConvolution:
preConvKernel = preConvPsf.getLocalKernel()
if self.config.useScoreImageDetection:
scoreExposure = self.decorrelate.run(exposureOrig, subtractRes.warpedExposure,
scoreExposure,
subtractRes.psfMatchingKernel,
spatiallyVarying=self.config.doSpatiallyVarying,
preConvKernel=preConvKernel,
templateMatched=True,
preConvMode=True).correctedExposure
# Note that the subtracted exposure is always decorrelated,
# even if the score image is used for detection
subtractedExposure = self.decorrelate.run(exposureOrig, subtractRes.warpedExposure,
subtractedExposure,
subtractRes.psfMatchingKernel,
spatiallyVarying=self.config.doSpatiallyVarying,
preConvKernel=None,
templateMatched=self.config.convolveTemplate,
preConvMode=False).correctedExposure
# END (if subtractAlgorithm == 'AL')
# END (if self.config.doSubtract)
if self.config.doDetection:
self.log.info("Running diaSource detection")
# subtractedExposure - reserved for task return value
# in zogy, it is always the proper difference image
# in AL, it may be (yet) pre-convolved and/or decorrelated
#
# detectionExposure - controls which exposure to use for detection
# in-place modifications will appear in task return
if self.config.useScoreImageDetection:
# zogy with score image detection enabled
self.log.info("Detection, diffim rescaling and measurements are "
"on AL likelihood or Zogy score image.")
detectionExposure = scoreExposure
else:
# AL or zogy with no score image detection
detectionExposure = subtractedExposure
# Rescale difference image variance plane
if self.config.doScaleDiffimVariance:
self.log.info("Rescaling diffim variance")
diffimVarFactor = self.scaleVariance.run(detectionExposure.getMaskedImage())
self.log.info("Diffim variance scaling factor: %.2f", diffimVarFactor)
self.metadata.add("scaleDiffimVarianceFactor", diffimVarFactor)
# Erase existing detection mask planes
mask = detectionExposure.getMaskedImage().getMask()
mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE"))
table = afwTable.SourceTable.make(self.schema, idFactory)
table.setMetadata(self.algMetadata)
results = self.detection.run(
table=table,
exposure=detectionExposure,
doSmooth=not self.config.useScoreImageDetection
)
if self.config.doMerge:
fpSet = results.fpSets.positive
fpSet.merge(results.fpSets.negative, self.config.growFootprint,
self.config.growFootprint, False)
diaSources = afwTable.SourceCatalog(table)
fpSet.makeSources(diaSources)
self.log.info("Merging detections into %d sources", len(diaSources))
else:
diaSources = results.sources
# Inject skySources before measurement.
if self.config.doSkySources:
skySourceFootprints = self.skySources.run(
mask=detectionExposure.mask,
seed=detectionExposure.info.id)
if skySourceFootprints:
for foot in skySourceFootprints:
s = diaSources.addNew()
s.setFootprint(foot)
s.set(self.skySourceKey, True)
if self.config.doMeasurement:
newDipoleFitting = self.config.doDipoleFitting
self.log.info("Running diaSource measurement: newDipoleFitting=%r", newDipoleFitting)
if not newDipoleFitting:
# Just fit dipole in diffim
self.measurement.run(diaSources, detectionExposure)
else:
# Use (matched) template and science image (if avail.) to constrain dipole fitting
if self.config.doSubtract and 'matchedExposure' in subtractRes.getDict():
self.measurement.run(diaSources, detectionExposure, exposure,
subtractRes.matchedExposure)
else:
self.measurement.run(diaSources, detectionExposure, exposure)
if self.config.doApCorr:
self.applyApCorr.run(
catalog=diaSources,
apCorrMap=detectionExposure.getInfo().getApCorrMap()
)
if self.config.doForcedMeasurement:
# Run forced psf photometry on the PVI at the diaSource locations.
# Copy the measured flux and error into the diaSource.
forcedSources = self.forcedMeasurement.generateMeasCat(
exposure, diaSources, detectionExposure.getWcs())
self.forcedMeasurement.run(forcedSources, exposure, diaSources, detectionExposure.getWcs())
mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFlux")[0],
"ip_diffim_forced_PsfFlux_instFlux", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFluxErr")[0],
"ip_diffim_forced_PsfFlux_instFluxErr", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_area")[0],
"ip_diffim_forced_PsfFlux_area", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag")[0],
"ip_diffim_forced_PsfFlux_flag", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_noGoodPixels")[0],
"ip_diffim_forced_PsfFlux_flag_noGoodPixels", True)
mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_edge")[0],
"ip_diffim_forced_PsfFlux_flag_edge", True)
for diaSource, forcedSource in zip(diaSources, forcedSources):
diaSource.assign(forcedSource, mapper)
# Match with the calexp sources if possible
if self.config.doMatchSources:
if selectSources is not None:
# Create key,val pair where key=diaSourceId and val=sourceId
matchRadAsec = self.config.diaSourceMatchRadius
matchRadPixel = matchRadAsec/exposure.getWcs().getPixelScale().asArcseconds()
srcMatches = afwTable.matchXy(selectSources, diaSources, matchRadPixel)
srcMatchDict = dict([(srcMatch.second.getId(), srcMatch.first.getId()) for
srcMatch in srcMatches])
self.log.info("Matched %d / %d diaSources to sources",
len(srcMatchDict), len(diaSources))
else:
self.log.warning("Src product does not exist; cannot match with diaSources")
srcMatchDict = {}
# Create key,val pair where key=diaSourceId and val=refId
refAstromConfig = AstrometryConfig()
refAstromConfig.matcher.maxMatchDistArcSec = matchRadAsec
refAstrometer = AstrometryTask(refAstromConfig)
astromRet = refAstrometer.run(exposure=exposure, sourceCat=diaSources)
refMatches = astromRet.matches
if refMatches is None:
self.log.warning("No diaSource matches with reference catalog")
refMatchDict = {}
else:
self.log.info("Matched %d / %d diaSources to reference catalog",
len(refMatches), len(diaSources))
refMatchDict = dict([(refMatch.second.getId(), refMatch.first.getId()) for
refMatch in refMatches])
# Assign source Ids
for diaSource in diaSources:
sid = diaSource.getId()
if sid in srcMatchDict:
diaSource.set("srcMatchId", srcMatchDict[sid])
if sid in refMatchDict:
diaSource.set("refMatchId", refMatchDict[sid])
if self.config.doAddMetrics and self.config.doSelectSources:
self.log.info("Evaluating metrics and control sample")
kernelCandList = []
for cell in subtractRes.kernelCellSet.getCellList():
for cand in cell.begin(False): # include bad candidates
kernelCandList.append(cand)
# Get basis list to build control sample kernels
basisList = kernelCandList[0].getKernel(KernelCandidateF.ORIG).getKernelList()
nparam = len(kernelCandList[0].getKernel(KernelCandidateF.ORIG).getKernelParameters())
controlCandList = (
diffimTools.sourceTableToCandidateList(controlSources,
subtractRes.warpedExposure, exposure,
self.config.subtract.kernel.active,
self.config.subtract.kernel.active.detectionConfig,
self.log, doBuild=True, basisList=basisList))
KernelCandidateQa.apply(kernelCandList, subtractRes.psfMatchingKernel,
subtractRes.backgroundModel, dof=nparam)
KernelCandidateQa.apply(controlCandList, subtractRes.psfMatchingKernel,
subtractRes.backgroundModel)
if self.config.doDetection:
KernelCandidateQa.aggregate(selectSources, self.metadata, allresids, diaSources)
else:
KernelCandidateQa.aggregate(selectSources, self.metadata, allresids)
self.runDebug(exposure, subtractRes, selectSources, kernelSources, diaSources)
return pipeBase.Struct(
subtractedExposure=subtractedExposure,
scoreExposure=scoreExposure,
warpedExposure=subtractRes.warpedExposure,
matchedExposure=subtractRes.matchedExposure,
subtractRes=subtractRes,
diaSources=diaSources,
selectSources=selectSources
)
def fitAstrometry(self, templateSources, templateExposure, selectSources):
Definition at line 1339 of file imageDifference.py.
The documentation for this class was generated from the following file:
- /j/snowflake/release/lsstsw/stack/lsst-scipipe-0.7.0/Linux64/pipe_tasks/21.0.0-172-gfb10e10a+18fedfabac/python/lsst/pipe/tasks/imageDifference.py
