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
modelPsfMatch.py
Go to the documentation of this file.
1 # LSST Data Management System
2 # Copyright 2008, 2009, 2010 LSST Corporation.
3 #
4 # This product includes software developed by the
5 # LSST Project (http://www.lsst.org/).
6 #
7 # This program is free software: you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation, either version 3 of the License, or
10 # (at your option) any later version.
11 #
12 # This program is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
16 #
17 # You should have received a copy of the LSST License Statement and
18 # the GNU General Public License along with this program. If not,
19 # see <http://www.lsstcorp.org/LegalNotices/>.
20 #
21 import numpy as num
22 import diffimLib
23 import lsst.afw.geom as afwGeom
24 import lsst.afw.image as afwImage
25 import lsst.afw.math as afwMath
26 import lsst.pex.logging as pexLog
27 import lsst.pex.config as pexConfig
28 import lsst.meas.algorithms as measAlg
29 import lsst.pipe.base as pipeBase
30 from makeKernelBasisList import makeKernelBasisList
31 from psfMatch import PsfMatchTask, PsfMatchConfigAL
32 from . import utils as diUtils
33 import lsst.afw.display.ds9 as ds9
34 
35 sigma2fwhm = 2. * num.sqrt(2. * num.log(2.))
36 
37 class ModelPsfMatchConfig(pexConfig.Config):
38  """!Configuration for model-to-model Psf matching"""
39 
40  kernel = pexConfig.ConfigChoiceField(
41  doc = "kernel type",
42  typemap = dict(
43  AL = PsfMatchConfigAL,
44  ),
45  default = "AL",
46  )
47 
48  def setDefaults(self):
49  # No sigma clipping
50  self.kernel.active.singleKernelClipping = False
51  self.kernel.active.kernelSumClipping = False
52  self.kernel.active.spatialKernelClipping = False
53  self.kernel.active.checkConditionNumber = False
54 
55  # Variance is ill defined
56  self.kernel.active.constantVarianceWeighting = True
57 
58  # Psfs are typically small; reduce the kernel size
59  self.kernel.active.kernelSizeMin = 11
60  self.kernel.active.kernelSize = 11
61 
62 ## \addtogroup LSST_task_documentation
63 ## \{
64 ## \page ModelPsfMatchTask
65 ## \ref ModelPsfMatchTask_ "ModelPsfMatchTask"
66 ## \copybrief ModelPsfMatchTask
67 ## \}
68 
69 class ModelPsfMatchTask(PsfMatchTask):
70  """!
71 \anchor ModelPsfMatchTask_
72 
73 \brief Matching of two model Psfs, and application of the Psf-matching kernel to an input Exposure
74 
75 \section ip_diffim_modelpsfmatch_Contents Contents
76 
77  - \ref ip_diffim_modelpsfmatch_Purpose
78  - \ref ip_diffim_modelpsfmatch_Initialize
79  - \ref ip_diffim_modelpsfmatch_IO
80  - \ref ip_diffim_modelpsfmatch_Config
81  - \ref ip_diffim_modelpsfmatch_Metadata
82  - \ref ip_diffim_modelpsfmatch_Debug
83  - \ref ip_diffim_modelpsfmatch_Example
84 
85 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
86 
87 \section ip_diffim_modelpsfmatch_Purpose Description
88 
89 This Task differs from ImagePsfMatchTask in that it matches two Psf _models_, by realizing
90 them in an Exposure-sized SpatialCellSet and then inserting each Psf-image pair into KernelCandidates.
91 Because none of the pairs of sources that are to be matched should be invalid, all sigma clipping is
92 turned off in ModelPsfMatchConfig. And because there is no tracked _variance_ in the Psf images, the
93 debugging and logging QA info should be interpreted with caution.
94 
95 One item of note is that the sizes of Psf models are fixed (e.g. its defined as a 21x21 matrix). When the
96 Psf-matching kernel is being solved for, the Psf "image" is convolved with each kernel basis function,
97 leading to a loss of information around the borders. This pixel loss will be problematic for the numerical
98 stability of the kernel solution if the size of the convolution kernel (set by ModelPsfMatchConfig.kernelSize)
99 is much bigger than: psfSize//2. Thus the sizes of Psf-model matching kernels are typically smaller
100 than their image-matching counterparts. If the size of the kernel is too small, the convolved stars will
101 look "boxy"; if the kernel is too large, the kernel solution will be "noisy". This is a trade-off that
102 needs careful attention for a given dataset.
103 
104 The primary use case for this Task is in matching an Exposure to a constant-across-the-sky Psf model for the
105 purposes of image coaddition. It is important to note that in the code, the "template" Psf is the Psf
106 that the science image gets matched to. In this sense the order of template and science image are
107 reversed, compared to ImagePsfMatchTask, which operates on the template image.
108 
109 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
110 
111 \section ip_diffim_modelpsfmatch_Initialize Task initialization
112 
113 \copydoc \_\_init\_\_
114 
115 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
116 
117 \section ip_diffim_modelpsfmatch_IO Invoking the Task
118 
119 \copydoc run
120 
121 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
122 
123 \section ip_diffim_modelpsfmatch_Config Configuration parameters
124 
125 See \ref ModelPsfMatchConfig
126 
127 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
128 
129 \section ip_diffim_modelpsfmatch_Metadata Quantities set in Metadata
130 
131 See \ref ip_diffim_psfmatch_Metadata "PsfMatchTask"
132 
133 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
134 
135 \section ip_diffim_modelpsfmatch_Debug Debug variables
136 
137 The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a
138 flag \c -d/--debug to import \b debug.py from your \c PYTHONPATH. The relevant contents of debug.py
139 for this Task include:
140 
141 \code{.py}
142  import sys
143  import lsstDebug
144  def DebugInfo(name):
145  di = lsstDebug.getInfo(name)
146  if name == "lsst.ip.diffim.psfMatch":
147  di.display = True # global
148  di.maskTransparency = 80 # ds9 mask transparency
149  di.displayCandidates = True # show all the candidates and residuals
150  di.displayKernelBasis = False # show kernel basis functions
151  di.displayKernelMosaic = True # show kernel realized across the image
152  di.plotKernelSpatialModel = False # show coefficients of spatial model
153  di.showBadCandidates = True # show the bad candidates (red) along with good (green)
154  elif name == "lsst.ip.diffim.modelPsfMatch":
155  di.display = True # global
156  di.maskTransparency = 30 # ds9 mask transparency
157  di.displaySpatialCells = True # show spatial cells before the fit
158  return di
159  lsstDebug.Info = DebugInfo
160  lsstDebug.frame = 1
161 \endcode
162 
163 Note that if you want addional logging info, you may add to your scripts:
164 \code{.py}
165 import lsst.pex.logging as pexLog
166 pexLog.Trace_setVerbosity('lsst.ip.diffim', 5)
167 \endcode
168 
169 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
170 
171 \section ip_diffim_modelpsfmatch_Example A complete example of using ModelPsfMatchTask
172 
173 This code is modelPsfMatchTask.py in the examples directory, and can be run as \em e.g.
174 \code
175 examples/modelPsfMatchTask.py
176 examples/modelPsfMatchTask.py --debug
177 examples/modelPsfMatchTask.py --debug --template /path/to/templateExp.fits --science /path/to/scienceExp.fits
178 \endcode
179 
180 \dontinclude modelPsfMatchTask.py
181 Create a subclass of ModelPsfMatchTask that accepts two exposures. Note that the "template" exposure
182 contains the Psf that will get matched to, and the "science" exposure is the one that will be convolved:
183 \skip MyModelPsfMatchTask
184 \until return
185 
186 And allow the user the freedom to either run the script in default mode, or point to their own images on disk.
187 Note that these images must be readable as an lsst.afw.image.Exposure:
188 \skip main
189 \until parse_args
190 
191 We have enabled some minor display debugging in this script via the --debug option. However, if you
192 have an lsstDebug debug.py in your PYTHONPATH you will get additional debugging displays. The following
193 block checks for this script:
194 \skip args.debug
195 \until sys.stderr
196 
197 \dontinclude modelPsfMatchTask.py
198 Finally, we call a run method that we define below. First set up a Config and modify some of the parameters.
199 In particular we don't want to "grow" the sizes of the kernel or KernelCandidates, since we are operating with
200 fixed--size images (i.e. the size of the input Psf models).
201 \skip run(args)
202 \until False
203 
204 Make sure the images (if any) that were sent to the script exist on disk and are readable. If no images
205 are sent, make some fake data up for the sake of this example script (have a look at the code if you want
206 more details on generateFakeData):
207 \skip requested
208 \until sizeCellY
209 
210 Display the two images if --debug:
211 \skip args.debug
212 \until Science
213 
214 Create and run the Task:
215 \skip Create
216 \until result
217 
218 And finally provide optional debugging display of the Psf-matched (via the Psf models) science image:
219 \skip args.debug
220 \until result.psfMatchedExposure
221 
222 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
223 
224  """
225  ConfigClass = ModelPsfMatchConfig
226 
227  def __init__(self, *args, **kwargs):
228  """!Create a ModelPsfMatchTask
229 
230  \param *args arguments to be passed to lsst.ip.diffim.PsfMatchTask.__init__
231  \param **kwargs keyword arguments to be passed to lsst.ip.diffim.PsfMatchTask.__init__
232 
233  Upon initialization, the kernel configuration is defined by self.config.kernel.active. This Task
234  does have a run() method, which is the default way to call the Task.
235  """
236  PsfMatchTask.__init__(self, *args, **kwargs)
237  self.kConfig = self.config.kernel.active
238 
239  @pipeBase.timeMethod
240  def run(self, exposure, referencePsfModel, kernelSum=1.0):
241  """!Psf-match an exposure to a model Psf
242 
243  @param exposure: Exposure to Psf-match to the reference Psf model;
244  it must return a valid PSF model via exposure.getPsf()
245  @param referencePsfModel: The Psf model to match to (an lsst.afw.detection.Psf)
246  @param kernelSum: A multipicative factor to apply to the kernel sum (default=1.0)
247 
248  @return
249  - psfMatchedExposure: the Psf-matched Exposure. This has the same parent bbox, Wcs, Calib and
250  Filter as the input Exposure but no Psf. In theory the Psf should equal referencePsfModel but
251  the match is likely not exact.
252  - psfMatchingKernel: the spatially varying Psf-matching kernel
253  - kernelCellSet: SpatialCellSet used to solve for the Psf-matching kernel
254 
255  Raise a RuntimeError if the Exposure does not contain a Psf model
256  """
257  if not exposure.hasPsf():
258  raise RuntimeError("exposure does not contain a Psf model")
259 
260  maskedImage = exposure.getMaskedImage()
261 
262  self.log.log(pexLog.Log.INFO, "compute Psf-matching kernel")
263  kernelCellSet = self._buildCellSet(exposure, referencePsfModel)
264  width, height = referencePsfModel.getLocalKernel().getDimensions()
265  psfAttr1 = measAlg.PsfAttributes(exposure.getPsf(), width//2, height//2)
266  psfAttr2 = measAlg.PsfAttributes(referencePsfModel, width//2, height//2)
267  s1 = psfAttr1.computeGaussianWidth(psfAttr1.ADAPTIVE_MOMENT) # gaussian sigma in pixels
268  s2 = psfAttr2.computeGaussianWidth(psfAttr2.ADAPTIVE_MOMENT) # gaussian sigma in pixels
269  fwhm1 = s1 * sigma2fwhm # science Psf
270  fwhm2 = s2 * sigma2fwhm # template Psf
271 
272  basisList = makeKernelBasisList(self.kConfig, fwhm1, fwhm2, metadata = self.metadata)
273  spatialSolution, psfMatchingKernel, backgroundModel = self._solve(kernelCellSet, basisList)
274 
275  if psfMatchingKernel.isSpatiallyVarying():
276  sParameters = num.array(psfMatchingKernel.getSpatialParameters())
277  sParameters[0][0] = kernelSum
278  psfMatchingKernel.setSpatialParameters(sParameters)
279  else:
280  kParameters = num.array(psfMatchingKernel.getKernelParameters())
281  kParameters[0] = kernelSum
282  psfMatchingKernel.setKernelParameters(kParameters)
283 
284  self.log.log(pexLog.Log.INFO, "Psf-match science exposure to reference")
285  psfMatchedExposure = afwImage.ExposureF(exposure.getBBox(), exposure.getWcs())
286  psfMatchedExposure.setFilter(exposure.getFilter())
287  psfMatchedExposure.setCalib(exposure.getCalib())
288  psfMatchedMaskedImage = psfMatchedExposure.getMaskedImage()
289 
290  # Normalize the psf-matching kernel while convolving since its magnitude is meaningless
291  # when PSF-matching one model to another.
292  doNormalize = True
293  afwMath.convolve(psfMatchedMaskedImage, maskedImage, psfMatchingKernel, doNormalize)
294 
295  self.log.log(pexLog.Log.INFO, "done")
296  return pipeBase.Struct(psfMatchedExposure=psfMatchedExposure,
297  psfMatchingKernel=psfMatchingKernel,
298  kernelCellSet=kernelCellSet,
299  metadata=self.metadata)
300 
301  def _diagnostic(self, kernelCellSet, spatialSolution, spatialKernel, spatialBg):
302  """!Print diagnostic information on spatial kernel and background fit
303 
304  The debugging diagnostics are not really useful here, since the images we are matching have
305  no variance. Thus override the _diagnostic method to generate no logging information"""
306  return
307 
308  def _buildCellSet(self, exposure, referencePsfModel):
309  """!Build a SpatialCellSet for use with the solve method
310 
311  @param exposure: The science exposure that will be convolved; must contain a Psf
312  @param referencePsfModel: Psf model to match to
313 
314  @return kernelCellSet: a SpatialCellSet to be used by self._solve
315 
316  Raise a RuntimeError if the reference Psf model and science Psf model have different dimensions
317  """
318  scienceBBox = exposure.getBBox()
319  sciencePsfModel = exposure.getPsf()
320  # The Psf base class does not support getKernel() in general, as there are some Psf
321  # classes for which this is not meaningful.
322  # Many Psfs we use in practice are KernelPsfs, and this algorithm will work fine for them,
323  # but someday it should probably be modified to support arbitrary Psfs.
324  referencePsfModel = measAlg.KernelPsf.swigConvert(referencePsfModel)
325  sciencePsfModel = measAlg.KernelPsf.swigConvert(sciencePsfModel)
326  if referencePsfModel is None or sciencePsfModel is None:
327  raise RuntimeError("ERROR: Psf matching is only implemented for KernelPsfs")
328  if (referencePsfModel.getKernel().getDimensions() != sciencePsfModel.getKernel().getDimensions()):
329  pexLog.Trace(self.log.getName(), 1,
330  "ERROR: Dimensions of reference Psf and science Psf different; exiting")
331  raise RuntimeError, "ERROR: Dimensions of reference Psf and science Psf different; exiting"
332 
333  psfWidth, psfHeight = referencePsfModel.getKernel().getDimensions()
334  maxKernelSize = min(psfWidth, psfHeight) - 1
335  if maxKernelSize % 2 == 0:
336  maxKernelSize -= 1
337  if self.kConfig.kernelSize > maxKernelSize:
338  raise ValueError, "Kernel size (%d) too big to match Psfs of size %d; reduce to at least %d" % (
339  self.kConfig.kernelSize, psfWidth, maxKernelSize)
340 
341  # Infer spatial order of Psf model!
342  #
343  # Infer from the number of spatial parameters.
344  # (O + 1) * (O + 2) / 2 = N
345  # O^2 + 3 * O + 2 * (1 - N) = 0
346  #
347  # Roots are [-3 +/- sqrt(9 - 8 * (1 - N))] / 2
348  #
349  nParameters = sciencePsfModel.getKernel().getNSpatialParameters()
350  root = num.sqrt(9 - 8 * (1 - nParameters))
351  if (root != root // 1): # We know its an integer solution
352  pexLog.Trace(self.log.getName(), 3, "Problem inferring spatial order of image's Psf")
353  else:
354  order = (root - 3) / 2
355  if (order != order // 1):
356  pexLog.Trace(self.log.getName(), 3, "Problem inferring spatial order of image's Psf")
357  else:
358  pexLog.Trace(self.log.getName(), 2, "Spatial order of Psf = %d; matching kernel order = %d" % (
359  order, self.kConfig.spatialKernelOrder))
360 
361  regionSizeX, regionSizeY = scienceBBox.getDimensions()
362  scienceX0, scienceY0 = scienceBBox.getMin()
363 
364  sizeCellX = self.kConfig.sizeCellX
365  sizeCellY = self.kConfig.sizeCellY
366 
367  kernelCellSet = afwMath.SpatialCellSet(
368  afwGeom.Box2I(afwGeom.Point2I(scienceX0, scienceY0),
369  afwGeom.Extent2I(regionSizeX, regionSizeY)),
370  sizeCellX, sizeCellY
371  )
372 
373  nCellX = regionSizeX // sizeCellX
374  nCellY = regionSizeY // sizeCellY
375  dimenR = referencePsfModel.getKernel().getDimensions()
376  dimenS = sciencePsfModel.getKernel().getDimensions()
377 
378  policy = pexConfig.makePolicy(self.kConfig)
379  for row in range(nCellY):
380  # place at center of cell
381  posY = sizeCellY * row + sizeCellY // 2 + scienceY0
382 
383  for col in range(nCellX):
384  # place at center of cell
385  posX = sizeCellX * col + sizeCellX // 2 + scienceX0
386 
387  pexLog.Trace(self.log.getName(), 5, "Creating Psf candidate at %.1f %.1f" % (posX, posY))
388 
389  # reference kernel image, at location of science subimage
390  kernelImageR = referencePsfModel.computeImage(afwGeom.Point2D(posX, posY)).convertF()
391  kernelMaskR = afwImage.MaskU(dimenR)
392  kernelMaskR.set(0)
393  kernelVarR = afwImage.ImageF(dimenR)
394  kernelVarR.set(1.0)
395  referenceMI = afwImage.MaskedImageF(kernelImageR, kernelMaskR, kernelVarR)
396 
397  # kernel image we are going to convolve
398  kernelImageS = sciencePsfModel.computeImage(afwGeom.Point2D(posX, posY)).convertF()
399  kernelMaskS = afwImage.MaskU(dimenS)
400  kernelMaskS.set(0)
401  kernelVarS = afwImage.ImageF(dimenS)
402  kernelVarS.set(1.0)
403  scienceMI = afwImage.MaskedImageF(kernelImageS, kernelMaskS, kernelVarS)
404 
405  # The image to convolve is the science image, to the reference Psf.
406  kc = diffimLib.makeKernelCandidate(posX, posY, scienceMI, referenceMI, policy)
407  kernelCellSet.insertCandidate(kc)
408 
409  import lsstDebug
410  display = lsstDebug.Info(__name__).display
411  displaySpatialCells = lsstDebug.Info(__name__).displaySpatialCells
412  maskTransparency = lsstDebug.Info(__name__).maskTransparency
413  if not maskTransparency:
414  maskTransparency = 0
415  if display:
416  ds9.setMaskTransparency(maskTransparency)
417  if display and displaySpatialCells:
418  diUtils.showKernelSpatialCells(exposure.getMaskedImage(), kernelCellSet,
419  symb="o", ctype=ds9.CYAN, ctypeUnused=ds9.YELLOW, ctypeBad=ds9.RED,
420  size=4, frame=lsstDebug.frame, title="Image to be convolved")
421  lsstDebug.frame += 1
422  return kernelCellSet
def _diagnostic
Print diagnostic information on spatial kernel and background fit.
def __init__
Create a ModelPsfMatchTask.
limited backward compatibility to the DC2 run-time trace facilities
Definition: Trace.h:93
An integer coordinate rectangle.
Definition: Box.h:53
A collection of SpatialCells covering an entire image.
Definition: SpatialCell.h:378
Matching of two model Psfs, and application of the Psf-matching kernel to an input Exposure...
Configuration for model-to-model Psf matching.
def _buildCellSet
Build a SpatialCellSet for use with the solve method.
def run
Psf-match an exposure to a model Psf.
void convolve(OutImageT &convolvedImage, InImageT const &inImage, KernelT const &kernel, bool doNormalize, bool doCopyEdge=false)
Old, deprecated version of convolve.