LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
makeKernelBasisList.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008-2016 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 
23 __all__ = ["makeKernelBasisList", "generateAlardLuptonBasisList"]
24 
25 from . import diffimLib
26 from lsst.log import Log
27 import numpy as np
28 
29 sigma2fwhm = 2. * np.sqrt(2. * np.log(2.))
30 
31 
32 def makeKernelBasisList(config, targetFwhmPix=None, referenceFwhmPix=None,
33  basisDegGauss=None, basisSigmaGauss=None, metadata=None):
34  """Generate the delta function or Alard-Lupton kernel bases depending on the Config.
35  Wrapper to call either `lsst.ip.diffim.makeDeltaFunctionBasisList` or
36  `lsst.ip.diffim.generateAlardLuptonBasisList`.
37 
38  Parameters
39  ----------
40  config : `lsst.ip.diffim.PsfMatchConfigAL`
41  Configuration object.
42  targetFwhmPix : `float`, optional
43  Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
44  Not used for delta function basis sets.
45  referenceFwhmPix : `float`, optional
46  Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
47  Not used for delta function basis sets.
48  basisDegGauss : `list` of `int`, optional
49  Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
50  Not used for delta function basis sets.
51  basisSigmaGauss : `list` of `int`, optional
52  Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
53  Not used for delta function basis sets.
54  metadata : `lsst.daf.base.PropertySet`, optional
55  Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
56  Not used for delta function basis sets.
57 
58  Returns
59  -------
60  basisList: `list` of `lsst.afw.math.kernel.FixedKernel`
61  List of basis kernels.
62 
63  Notes
64  -----
65  See `lsst.ip.diffim.generateAlardLuptonBasisList` and
66  `lsst.ip.diffim.makeDeltaFunctionBasisList` for more information.
67 
68  Raises
69  ------
70  ValueError
71  If ``config.kernelBasisSet`` has an invalid value (not "alard-lupton" or "delta-function").
72  """
73  if config.kernelBasisSet == "alard-lupton":
74  return generateAlardLuptonBasisList(config, targetFwhmPix=targetFwhmPix,
75  referenceFwhmPix=referenceFwhmPix,
76  basisDegGauss=basisDegGauss,
77  basisSigmaGauss=basisSigmaGauss,
78  metadata=metadata)
79  elif config.kernelBasisSet == "delta-function":
80  kernelSize = config.kernelSize
81  return diffimLib.makeDeltaFunctionBasisList(kernelSize, kernelSize)
82  else:
83  raise ValueError("Cannot generate %s basis set" % (config.kernelBasisSet))
84 
85 
86 def generateAlardLuptonBasisList(config, targetFwhmPix=None, referenceFwhmPix=None,
87  basisDegGauss=None, basisSigmaGauss=None, metadata=None):
88  """Generate an Alard-Lupton kernel basis list based upon the Config and
89  the input FWHM of the science and template images.
90 
91  Parameters
92  ----------
93  config : `lsst.ip.diffim.PsfMatchConfigAL`
94  Configuration object for the Alard-Lupton algorithm.
95  targetFwhmPix : `float`, optional
96  Fwhm width (pixel) of the template exposure characteristic psf.
97  This is the _target_ that will be matched to the science exposure.
98  referenceFwhmPix : `float`, optional
99  Fwhm width (pixel) of the science exposure characteristic psf.
100  basisDegGauss : `list` of `int`, optional
101  Polynomial degree of each Gaussian (sigma) basis. If None, defaults to `config.alardDegGauss`.
102  basisSigmaGauss : `list` of `int`, optional
103  Sigmas of each Gaussian basis. If None, defaults to `config.alardSigGauss`.
104  metadata : `lsst.daf.base.PropertySet`, optional
105  If specified, object to collect metadata fields about the kernel basis list.
106 
107  Returns
108  -------
109  basisList : `list` of `lsst.afw.math.kernel.FixedKernel`
110  List of basis kernels. For each degree value ``n`` in ``config.basisDegGauss`` (n+2)(n+1)/2 kernels
111  are generated and appended to the list in the order of the polynomial parameter number.
112  See `lsst.afw.math.polynomialFunction2D` documentation for more details.
113 
114  Raises
115  ------
116  RuntimeError
117  - if ``config.kernelBasisSet`` is not equal to "alard-lupton"
118  ValueError
119  - if ``config.kernelSize`` is even
120  - if the number of Gaussians and the number of given
121  sigma values are not equal or
122  - if the number of Gaussians and the number of given
123  polynomial degree values are not equal
124 
125  Notes
126  -----
127  The polynomial functions (``f``) are always evaluated in the -1.0, +1.0 range in both x, y directions,
128  edge to edge, with ``f(0,0)`` evaluated at the kernel center pixel, ``f(-1.0,-1.0)`` at the kernel
129  ``(0,0)`` pixel. They are not scaled by the sigmas of the Gaussians.
130 
131  Base Gaussian widths (sigmas in pixels) of the kernels are determined as:
132  - If not all fwhm parameters are provided or ``config.scaleByFwhm==False``
133  then ``basisSigmaGauss`` is used. If ``basisSigmaGauss`` is not
134  provided, then ``config.alardSigGauss`` is used. In both cases, the
135  length of sigmas must be equal to ``config.alardNGauss``.
136  - If ``targetFwhmPix<referenceFwhmPix`` (normal convolution):
137  First sigma ``Sig_K`` is determined to satisfy: ``Sig_reference**2 = Sig_target**2 + Sig_K**2``.
138  If it's larger than ``config.alardMinSig * config.alardGaussBeta``, make it the
139  second kernel. Else make it the smallest kernel, unless only 1 kernel is asked for.
140  - If ``referenceFwhmPix < targetFwhmPix`` (deconvolution):
141  Define the progression of Gaussians using a
142  method to derive a deconvolution sum-of-Gaussians from it's
143  convolution counterpart. [1]_ Only use 3 since the algorithm
144  assumes 3 components.
145 
146  **Metadata fields**
147 
148  ALBasisNGauss : `int`
149  The number of base Gaussians in the AL basis functions.
150  ALBasisDegGauss : `list` of `int`
151  Polynomial order of spatial modification of the base Gaussian functions.
152  ALBasisSigGauss : `list` of `float`
153  Sigmas in pixels of the base Gaussians.
154  ALKernelSize : `int`
155  Kernel stamp size is (ALKernelSize pix, ALKernelSize pix).
156  ALBasisMode : `str`, either of ``config``, ``convolution``, ``deconvolution``
157  Indicates whether the config file values, the convolution or deconvolution algorithm
158  was used to determine the base Gaussian sigmas and the kernel stamp size.
159 
160  References
161  ----------
162  .. [1] Ulmer, W.: Inverse problem of linear combinations of Gaussian convolution kernels
163  (deconvolution) and some applications to proton/photon dosimetry and image
164  processing. http://iopscience.iop.org/0266-5611/26/8/085002 Equation 40
165  """
166 
167  if config.kernelBasisSet != "alard-lupton":
168  raise RuntimeError("Cannot generate %s basis within generateAlardLuptonBasisList" %
169  config.kernelBasisSet)
170 
171  kernelSize = config.kernelSize
172  fwhmScaling = config.kernelSizeFwhmScaling
173  basisNGauss = config.alardNGauss
174  basisGaussBeta = config.alardGaussBeta
175  basisMinSigma = config.alardMinSig
176  if basisDegGauss is None:
177  basisDegGauss = config.alardDegGauss
178  if basisSigmaGauss is None:
179  basisSigmaGauss = config.alardSigGauss
180 
181  if len(basisDegGauss) != basisNGauss:
182  raise ValueError("len(basisDegGauss) != basisNGauss : %d vs %d" % (len(basisDegGauss), basisNGauss))
183  if len(basisSigmaGauss) != basisNGauss:
184  raise ValueError("len(basisSigmaGauss) != basisNGauss : %d vs %d" %
185  (len(basisSigmaGauss), basisNGauss))
186  if (kernelSize % 2) != 1:
187  raise ValueError("Only odd-sized Alard-Lupton bases allowed")
188 
189  logger = Log.getLogger("lsst.ip.diffim.generateAlardLuptonBasisList")
190  if (targetFwhmPix is None) or (referenceFwhmPix is None) or (not config.scaleByFwhm):
191  logger.info("PSF sigmas are not available or scaling by fwhm disabled, "
192  "falling back to config values")
193  if metadata is not None:
194  metadata.add("ALBasisNGauss", basisNGauss)
195  metadata.add("ALBasisDegGauss", basisDegGauss)
196  metadata.add("ALBasisSigGauss", basisSigmaGauss)
197  metadata.add("ALKernelSize", kernelSize)
198  metadata.add("ALBasisMode", "config")
199 
200  return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
201 
202  targetSigma = targetFwhmPix / sigma2fwhm
203  referenceSigma = referenceFwhmPix / sigma2fwhm
204  logger.debug("Generating matching bases for sigma %.2f pix -> %.2f pix", targetSigma, referenceSigma)
205 
206  # Modify the size of Alard Lupton kernels based upon the images FWHM
207  #
208  # Note the operation is : template.x.kernel = science
209  #
210  # Assuming the template and science image Psfs are Gaussians with
211  # the Fwhm above, Fwhm_T **2 + Fwhm_K **2 = Fwhm_S **2
212  #
213  if targetSigma == referenceSigma:
214  # Leave defaults as-is
215  logger.debug("Target and reference psf fwhms are equal, falling back to config values")
216  basisMode = "config"
217  elif referenceSigma > targetSigma:
218  # Normal convolution
219 
220  # First Gaussian has the sigma that comes from the convolution
221  # of two Gaussians : Sig_S**2 = Sig_T**2 + Sig_K**2
222  #
223  # If it's larger than basisMinSigma * basisGaussBeta, make it the
224  # second kernel. Else make it the smallest kernel. Unless
225  # only 1 kernel is asked for.
226  logger.debug("Reference psf fwhm is the greater, normal convolution mode")
227  basisMode = "convolution"
228  kernelSigma = np.sqrt(referenceSigma**2 - targetSigma**2)
229  if kernelSigma < basisMinSigma:
230  kernelSigma = basisMinSigma
231 
232  basisSigmaGauss = []
233  if basisNGauss == 1:
234  basisSigmaGauss.append(kernelSigma)
235  nAppended = 1
236  else:
237  if (kernelSigma/basisGaussBeta) > basisMinSigma:
238  basisSigmaGauss.append(kernelSigma/basisGaussBeta)
239  basisSigmaGauss.append(kernelSigma)
240  nAppended = 2
241  else:
242  basisSigmaGauss.append(kernelSigma)
243  nAppended = 1
244 
245  # Any other Gaussians above basisNGauss=1 come from a scaling
246  # relationship: Sig_i+1 / Sig_i = basisGaussBeta
247  for i in range(nAppended, basisNGauss):
248  basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
249 
250  kernelSize = int(fwhmScaling * basisSigmaGauss[-1])
251  kernelSize += 0 if kernelSize%2 else 1 # Make sure it's odd
252  kernelSize = min(config.kernelSizeMax, max(kernelSize, config.kernelSizeMin))
253 
254  else:
255  # Deconvolution; Define the progression of Gaussians using a
256  # method to derive a deconvolution sum-of-Gaussians from it's
257  # convolution counterpart. Only use 3 since the algorithm
258  # assumes 3 components.
259  #
260  # http://iopscience.iop.org/0266-5611/26/8/085002 Equation 40
261 
262  # Use specializations for deconvolution
263  logger.debug("Target psf fwhm is the greater, deconvolution mode")
264  basisMode = "deconvolution"
265  basisNGauss = config.alardNGaussDeconv
266  basisMinSigma = config.alardMinSigDeconv
267 
268  kernelSigma = np.sqrt(targetSigma**2 - referenceSigma**2)
269  if kernelSigma < basisMinSigma:
270  kernelSigma = basisMinSigma
271 
272  basisSigmaGauss = []
273  if (kernelSigma/basisGaussBeta) > basisMinSigma:
274  basisSigmaGauss.append(kernelSigma/basisGaussBeta)
275  basisSigmaGauss.append(kernelSigma)
276  nAppended = 2
277  else:
278  basisSigmaGauss.append(kernelSigma)
279  nAppended = 1
280 
281  for i in range(nAppended, basisNGauss):
282  basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
283 
284  kernelSize = int(fwhmScaling * basisSigmaGauss[-1])
285  kernelSize += 0 if kernelSize%2 else 1 # Make sure it's odd
286  kernelSize = min(config.kernelSizeMax, max(kernelSize, config.kernelSizeMin))
287 
288  # Now build a deconvolution set from these sigmas
289  sig0 = basisSigmaGauss[0]
290  sig1 = basisSigmaGauss[1]
291  sig2 = basisSigmaGauss[2]
292  basisSigmaGauss = []
293  for n in range(1, 3):
294  for j in range(n):
295  sigma2jn = (n - j)*sig1**2
296  sigma2jn += j * sig2**2
297  sigma2jn -= (n + 1)*sig0**2
298  sigmajn = np.sqrt(sigma2jn)
299  basisSigmaGauss.append(sigmajn)
300 
301  basisSigmaGauss.sort()
302  basisNGauss = len(basisSigmaGauss)
303  basisDegGauss = [config.alardDegGaussDeconv for x in basisSigmaGauss]
304 
305  if metadata is not None:
306  metadata.add("ALBasisNGauss", basisNGauss)
307  metadata.add("ALBasisDegGauss", basisDegGauss)
308  metadata.add("ALBasisSigGauss", basisSigmaGauss)
309  metadata.add("ALKernelSize", kernelSize)
310  metadata.add("ALBasisMode", basisMode)
311 
312  logger.debug("basisSigmaGauss: %s basisDegGauss: %s",
313  ','.join(['{:.1f}'.format(v) for v in basisSigmaGauss]),
314  ','.join(['{:d}'.format(v) for v in basisDegGauss]))
315 
316  return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
int min
int max
def generateAlardLuptonBasisList(config, targetFwhmPix=None, referenceFwhmPix=None, basisDegGauss=None, basisSigmaGauss=None, metadata=None)
def makeKernelBasisList(config, targetFwhmPix=None, referenceFwhmPix=None, basisDegGauss=None, basisSigmaGauss=None, metadata=None)
Definition: Log.h:717
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174