diff --git a/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py b/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py index d6e9f69f..efd3ba86 100644 --- a/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py +++ b/python/lsst/ip/diffim/computeSpatiallySampledMetrics.py @@ -27,7 +27,7 @@ import lsst.pipe.base as pipeBase import lsst.pex.config as pexConfig -from lsst.ip.diffim.utils import getPsfFwhm, angleMean, evaluateMaskFraction +from lsst.ip.diffim.utils import getPsfFwhm, angleMean, evaluateMaskFraction, getKernelCenterDisplacement from lsst.meas.algorithms import SkyObjectsTask from lsst.pex.exceptions import InvalidParameterError, RangeError from lsst.utils.timer import timeMethod @@ -72,6 +72,12 @@ class SpatiallySampledMetricsConnections(pipeBase.PipelineTaskConnections, storageClass="SourceCatalog", name="{fakesType}{coaddName}Diff_candidateDiaSrc", ) + psfMatchingKernel = pipeBase.connectionTypes.Input( + doc="Kernel used to PSF match the science and template images.", + dimensions=("instrument", "visit", "detector"), + storageClass="MatchingKernel", + name="{fakesType}{coaddName}Diff_psfMatchKernel", + ) spatiallySampledMetrics = pipeBase.connectionTypes.Output( doc="Summary metrics computed at randomized locations.", dimensions=("instrument", "visit", "detector"), @@ -164,9 +170,32 @@ def __init__(self, **kwargs): "%s_mask_fraction"%maskPlane.lower(), "F", "Fraction of pixels with %s mask"%maskPlane ) + self.schema.addField( + "psfMatchingKernel_sum", "F", + "PSF matching kernel sum at location.") + self.schema.addField( + "psfMatchingKernel_dx", "F", + "PSF matching kernel centroid offset in x at location.", + units="pixel") + self.schema.addField( + "psfMatchingKernel_dy", "F", + "PSF matching kernel centroid offset in y at location.", + units="pixel") + self.schema.addField( + "psfMatchingKernel_length", "F", + "PSF matching kernel centroid offset module.", + units="arcsecond") + self.schema.addField( + "psfMatchingKernel_position_angle", "F", + "PSF matching kernel centroid offset position angle.", + units="radian") + self.schema.addField( + "psfMatchingKernel_direction", "F", + "PSF matching kernel centroid offset direction in detector plane.", + units="radian") @timeMethod - def run(self, science, matchedTemplate, template, difference, diaSources): + def run(self, science, matchedTemplate, template, difference, diaSources, psfMatchingKernel): """Calculate difference image metrics on specific locations across the images Parameters @@ -183,6 +212,8 @@ def run(self, science, matchedTemplate, template, difference, diaSources): Result of subtracting template from the science image. diaSources : `lsst.afw.table.SourceCatalog` The catalog of detected sources. + psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel` + The PSF matching kernel of the subtraction to evaluate. Returns ------- @@ -207,12 +238,13 @@ def run(self, science, matchedTemplate, template, difference, diaSources): for src in spatiallySampledMetrics: self._evaluateLocalMetric(src, science, matchedTemplate, template, difference, diaSources, - metricsMaskPlanes=metricsMaskPlanes) + metricsMaskPlanes=metricsMaskPlanes, + psfMatchingKernel=psfMatchingKernel) return pipeBase.Struct(spatiallySampledMetrics=spatiallySampledMetrics.asAstropy()) def _evaluateLocalMetric(self, src, science, matchedTemplate, template, difference, diaSources, - metricsMaskPlanes): + metricsMaskPlanes, psfMatchingKernel): """Calculate image quality metrics at spatially sampled locations. Parameters @@ -229,6 +261,8 @@ def _evaluateLocalMetric(self, src, science, matchedTemplate, template, differen Result of subtracting template from the science image. metricsMaskPlanes : `list` of `str` Mask planes to calculate metrics from. + psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel` + The PSF matching kernel of the subtraction to evaluate. """ bbox = src.getFootprint().getBBox() pix = bbox.getCenter() @@ -271,3 +305,24 @@ def _evaluateLocalMetric(self, src, science, matchedTemplate, template, differen src.set("%s_mask_fraction"%maskPlane.lower(), evaluateMaskFraction(difference.mask[bbox], maskPlane) ) + + krnlSum, dx, dy, direction, length = getKernelCenterDisplacement( + psfMatchingKernel, src.get('x'), src.get('y')) + + point1 = lsst.geom.SpherePoint( + src.get('coord_ra'), src.get('coord_dec'), + lsst.geom.radians) + point2 = science.wcs.pixelToSky(src.get('x') + dx, src.get('y') + dy) + bearing = point1.bearingTo(point2) + pa_ref_angle = lsst.geom.Angle(np.pi/2, lsst.geom.radians) + pa = pa_ref_angle - bearing + # Wrap around to get Delta_RA from -pi to +pi + pa = pa.wrapCtr() + position_angle = pa.asRadians() + + src.set('psfMatchingKernel_sum', krnlSum) + src.set('psfMatchingKernel_dx', dx) + src.set('psfMatchingKernel_dy', dy) + src.set('psfMatchingKernel_length', length*pixScale.asArcseconds()) + src.set('psfMatchingKernel_position_angle', position_angle) # in E of N position angle + src.set('psfMatchingKernel_direction', direction) # direction offset in detector diff --git a/python/lsst/ip/diffim/utils.py b/python/lsst/ip/diffim/utils.py index 851f7961..9d66349e 100644 --- a/python/lsst/ip/diffim/utils.py +++ b/python/lsst/ip/diffim/utils.py @@ -22,7 +22,7 @@ """Support utilities for Measuring sources""" # Export DipoleTestImage to expose fake image generating funcs -__all__ = ["DipoleTestImage", "evaluateMeanPsfFwhm", "getPsfFwhm"] +__all__ = ["DipoleTestImage", "evaluateMeanPsfFwhm", "getPsfFwhm", "getKernelCenterDisplacement"] import itertools import numpy as np @@ -853,6 +853,65 @@ def detectDipoleSources(self, doMerge=True, diffim=None, detectSigma=5.5, grow=3 return detectTask, schema +def getKernelCenterDisplacement(kernel, x, y, image=None): + """Calculate the PSF matching kernel peak offset from the nominal + position. + + Parameters + ---------- + kernel : `~lsst.afw.math.LinearCombinationKernel` + The PSF matching kernel to evaluate. + x : `float` + The x position on the detector to evaluate the kernel + y : `float` + The y position on the detector to evaluate the kernel + image : `~lsst.afw.image.ImageD` + The image to use as base for computing kernel pixel values + + Returns + ------- + kernel_sum : `float` + The sum of the kernel on the desired location + dx : `float` + The displacement of the kernel averaged peak, with respect to the + center of the extraction of the kernel + dy : `float` + The displacement of the kernel averaged peak, with respect to the + center of the extraction of the kernel + pos_angle: `float` + The position angle in detector coordinates of the displacement + length : `float` + The displacement module of the kernel centroid in pixel units + """ + + if image is None: + image = afwImage.ImageD(kernel.getDimensions()) + + # obtain the kernel image + hsize = kernel.getWidth()//2 + kernel_sum = kernel.computeImage(image, doNormalize=False, x=x, y=y) + + data = image.array + h, w = data.shape + xx = np.arange(w) + yy = np.arange(h) + + # create sum vectors and estimate weighted average + vx = data.sum(axis=0) + vx /= vx.sum() + dx = np.dot(vx, xx) - hsize + + vy = data.sum(axis=1) + vy /= vy.sum() + dy = np.dot(vy, yy) - hsize + + # obtain position angle and norm of displacement + pos_angle = np.arctan2(dy, dx) + length = np.sqrt(dx**2 + dy**2) + + return kernel_sum, dx, dy, pos_angle, length + + def getPsfFwhm(psf, average=True, position=None): """Directly calculate the horizontal and vertical widths of a PSF at half its maximum value.