LSST Applications g04e9c324dd+8c5ae1fdc5,g134cb467dc+1b3060144d,g18429d2f64+f642bf4753,g199a45376c+0ba108daf9,g1fd858c14a+2dcf163641,g262e1987ae+7b8c96d2ca,g29ae962dfc+3bd6ecb08a,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+53e1a9e7c5,g4595892280+fef73a337f,g47891489e3+2efcf17695,g4d44eb3520+642b70b07e,g53246c7159+8c5ae1fdc5,g67b6fd64d1+2efcf17695,g67fd3c3899+b70e05ef52,g74acd417e5+317eb4c7d4,g786e29fd12+668abc6043,g87389fa792+8856018cbb,g89139ef638+2efcf17695,g8d7436a09f+3be3c13596,g8ea07a8fe4+9f5ccc88ac,g90f42f885a+a4e7b16d9b,g97be763408+ad77d7208f,g9dd6db0277+b70e05ef52,ga681d05dcb+a3f46e7fff,gabf8522325+735880ea63,gac2eed3f23+2efcf17695,gb89ab40317+2efcf17695,gbf99507273+8c5ae1fdc5,gd8ff7fe66e+b70e05ef52,gdab6d2f7ff+317eb4c7d4,gdc713202bf+b70e05ef52,gdfd2d52018+b10e285e0f,ge365c994fd+310e8507c4,ge410e46f29+2efcf17695,geaed405ab2+562b3308c0,gffca2db377+8c5ae1fdc5,w.2025.35
LSST Data Management Base Package
Loading...
Searching...
No Matches
makeKernelBasisList.py
Go to the documentation of this file.
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
25from . import diffimLib
26import numpy as np
27
28from lsst.utils.logging import getLogger
29
30sigma2fwhm = 2. * np.sqrt(2. * np.log(2.))
31
32
33def makeKernelBasisList(config, targetFwhmPix=None, referenceFwhmPix=None,
34 basisDegGauss=None, basisSigmaGauss=None, metadata=None):
35 """Generate the delta function or Alard-Lupton kernel bases depending on the Config.
36 Wrapper to call either `lsst.ip.diffim.makeDeltaFunctionBasisList` or
37 `lsst.ip.diffim.generateAlardLuptonBasisList`.
38
39 Parameters
40 ----------
41 config : `lsst.ip.diffim.PsfMatchConfigAL`
42 Configuration object.
43 targetFwhmPix : `float`, optional
44 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
45 Not used for delta function basis sets.
46 referenceFwhmPix : `float`, optional
47 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
48 Not used for delta function basis sets.
49 basisDegGauss : `list` of `int`, optional
50 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
51 Not used for delta function basis sets.
52 basisSigmaGauss : `list` of `int`, optional
53 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
54 Not used for delta function basis sets.
55 metadata : `lsst.daf.base.PropertySet`, optional
56 Passed on to `lsst.ip.diffim.generateAlardLuptonBasisList`.
57 Not used for delta function basis sets.
58
59 Returns
60 -------
61 basisList: `list` of `lsst.afw.math.kernel.FixedKernel`
62 List of basis kernels.
63
64 Notes
65 -----
66 See `lsst.ip.diffim.generateAlardLuptonBasisList` and
67 `lsst.ip.diffim.makeDeltaFunctionBasisList` for more information.
68
69 Raises
70 ------
71 ValueError
72 If ``config.kernelBasisSet`` has an invalid value (not "alard-lupton" or "delta-function").
73 """
74 if config.kernelBasisSet == "alard-lupton":
75 return generateAlardLuptonBasisList(config, targetFwhmPix=targetFwhmPix,
76 referenceFwhmPix=referenceFwhmPix,
77 basisDegGauss=basisDegGauss,
78 basisSigmaGauss=basisSigmaGauss,
79 metadata=metadata)
80 elif config.kernelBasisSet == "delta-function":
81 kernelSize = config.kernelSize
82 return diffimLib.makeDeltaFunctionBasisList(kernelSize, kernelSize)
83 else:
84 raise ValueError("Cannot generate %s basis set" % (config.kernelBasisSet))
85
86
87def generateAlardLuptonBasisList(config, targetFwhmPix=None, referenceFwhmPix=None,
88 basisDegGauss=None, basisSigmaGauss=None, metadata=None):
89 """Generate an Alard-Lupton kernel basis list based upon the Config and
90 the input FWHM of the science and template images.
91
92 Parameters
93 ----------
94 config : `lsst.ip.diffim.PsfMatchConfigAL`
95 Configuration object for the Alard-Lupton algorithm.
96 targetFwhmPix : `float`, optional
97 Fwhm width (pixel) of the template exposure characteristic psf.
98 This is the _target_ that will be matched to the science exposure.
99 referenceFwhmPix : `float`, optional
100 Fwhm width (pixel) of the science exposure characteristic psf.
101 basisDegGauss : `list` of `int`, optional
102 Polynomial degree of each Gaussian (sigma) basis. If None, defaults to `config.alardDegGauss`.
103 basisSigmaGauss : `list` of `int`, optional
104 Sigmas of each Gaussian basis. If None, defaults to `config.alardSigGauss`.
105 metadata : `lsst.daf.base.PropertySet`, optional
106 If specified, object to collect metadata fields about the kernel basis list.
107
108 Returns
109 -------
110 basisList : `list` of `lsst.afw.math.kernel.FixedKernel`
111 List of basis kernels. For each degree value ``n`` in ``config.basisDegGauss`` (n+2)(n+1)/2 kernels
112 are generated and appended to the list in the order of the polynomial parameter number.
113 See `lsst.afw.math.polynomialFunction2D` documentation for more details.
114
115 Raises
116 ------
117 RuntimeError
118 - if ``config.kernelBasisSet`` is not equal to "alard-lupton"
119 ValueError
120 - if ``config.kernelSize`` is even
121 - if the number of Gaussians and the number of given
122 sigma values are not equal or
123 - if the number of Gaussians and the number of given
124 polynomial degree values are not equal
125
126 Notes
127 -----
128 The polynomial functions (``f``) are always evaluated in the -1.0, +1.0 range in both x, y directions,
129 edge to edge, with ``f(0,0)`` evaluated at the kernel center pixel, ``f(-1.0,-1.0)`` at the kernel
130 ``(0,0)`` pixel. They are not scaled by the sigmas of the Gaussians.
131
132 Base Gaussian widths (sigmas in pixels) of the kernels are determined as:
133 - If not all fwhm parameters are provided or ``config.scaleByFwhm==False``
134 then ``basisSigmaGauss`` is used. If ``basisSigmaGauss`` is not
135 provided, then ``config.alardSigGauss`` is used. In both cases, the
136 length of sigmas must be equal to ``config.alardNGauss``.
137 - If ``targetFwhmPix<referenceFwhmPix`` (normal convolution):
138 First sigma ``Sig_K`` is determined to satisfy: ``Sig_reference**2 = Sig_target**2 + Sig_K**2``.
139 If it's larger than ``config.alardMinSig * config.alardGaussBeta``, make it the
140 second kernel. Else make it the smallest kernel, unless only 1 kernel is asked for.
141 - If ``referenceFwhmPix < targetFwhmPix`` (deconvolution):
142 Compute the basis list with a significantly smaller ``config.alardMinSig``.
143
144 **Metadata fields**
145
146 ALBasisNGauss : `int`
147 The number of base Gaussians in the AL basis functions.
148 ALBasisDegGauss : `list` of `int`
149 Polynomial order of spatial modification of the base Gaussian functions.
150 ALBasisSigGauss : `list` of `float`
151 Sigmas in pixels of the base Gaussians.
152 ALKernelSize : `int`
153 Kernel stamp size is (ALKernelSize pix, ALKernelSize pix).
154 ALBasisMode : `str`, either of ``config``, ``convolution``, ``deconvolution``
155 Indicates whether the config file values, the convolution or deconvolution algorithm
156 was used to determine the base Gaussian sigmas and the kernel stamp size.
157
158 References
159 ----------
160 .. [1] Ulmer, W.: Inverse problem of linear combinations of Gaussian convolution kernels
161 (deconvolution) and some applications to proton/photon dosimetry and image
162 processing. http://iopscience.iop.org/0266-5611/26/8/085002 Equation 40
163 """
164
165 if config.kernelBasisSet != "alard-lupton":
166 raise RuntimeError("Cannot generate %s basis within generateAlardLuptonBasisList" %
167 config.kernelBasisSet)
168
169 kernelSize = config.kernelSize
170 fwhmScaling = config.kernelSizeFwhmScaling
171 basisNGauss = config.alardNGauss
172 basisGaussBeta = config.alardGaussBeta
173 basisMinSigma = config.alardMinSig
174 if basisDegGauss is None:
175 basisDegGauss = config.alardDegGauss
176 if basisSigmaGauss is None:
177 basisSigmaGauss = config.alardSigGauss
178
179 if len(basisDegGauss) != basisNGauss:
180 raise ValueError("len(basisDegGauss) != basisNGauss : %d vs %d" % (len(basisDegGauss), basisNGauss))
181 if len(basisSigmaGauss) != basisNGauss:
182 raise ValueError("len(basisSigmaGauss) != basisNGauss : %d vs %d" %
183 (len(basisSigmaGauss), basisNGauss))
184 if (kernelSize % 2) != 1:
185 raise ValueError("Only odd-sized Alard-Lupton bases allowed")
186
187 logger = getLogger("lsst.ip.diffim.generateAlardLuptonBasisList")
188 if (targetFwhmPix is None) or (referenceFwhmPix is None) or (not config.scaleByFwhm):
189 logger.info("PSF sigmas are not available or scaling by fwhm disabled, "
190 "falling back to config values")
191 if metadata is not None:
192 metadata.add("ALBasisNGauss", basisNGauss)
193 metadata.add("ALBasisDegGauss", basisDegGauss)
194 metadata.add("ALBasisSigGauss", basisSigmaGauss)
195 metadata.add("ALKernelSize", kernelSize)
196 metadata.add("ALBasisMode", "config")
197
198 return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
199
200 targetSigma = targetFwhmPix / sigma2fwhm
201 referenceSigma = referenceFwhmPix / sigma2fwhm
202 logger.debug("Generating matching bases for sigma %.2f pix -> %.2f pix", targetSigma, referenceSigma)
203
204 # Modify the size of Alard Lupton kernels based upon the images FWHM
205 #
206 # Note the operation is : template.x.kernel = science
207 #
208 # Assuming the template and science image Psfs are Gaussians with
209 # the Fwhm above, Fwhm_T **2 + Fwhm_K **2 = Fwhm_S **2
210 #
211 if targetSigma == referenceSigma:
212 # Leave defaults as-is
213 logger.debug("Target and reference psf fwhms are equal, falling back to config values")
214 basisMode = "config"
215 elif referenceSigma > targetSigma:
216 # Normal convolution
217
218 # First Gaussian has the sigma that comes from the convolution
219 # of two Gaussians : Sig_S**2 = Sig_T**2 + Sig_K**2
220 #
221 # If it's larger than basisMinSigma * basisGaussBeta, make it the
222 # second kernel. Else make it the smallest kernel. Unless
223 # only 1 kernel is asked for.
224 logger.debug("Reference psf fwhm is the greater, normal convolution mode")
225 basisMode = "convolution"
226
227 basisSigmaGauss = _calculateBasisSigmas(referenceSigma,
228 targetSigma,
229 basisMinSigma,
230 basisGaussBeta,
231 basisNGauss,
232 )
233
234 else:
235 # Swap order of reference and target sigmas,
236 # and potentially allow smaller basis modes.
237 logger.debug("Target psf fwhm is the greater, deconvolution mode")
238 basisMode = "deconvolution"
239 basisMinSigma = config.alardMinSigDeconv
240 basisSigmaGauss = _calculateBasisSigmas(targetSigma,
241 referenceSigma,
242 basisMinSigma,
243 basisGaussBeta,
244 basisNGauss,
245 )
246
247 basisNGauss = len(basisSigmaGauss)
248
249 kernelSize = int(fwhmScaling * basisSigmaGauss[-1])
250 kernelSize += 0 if kernelSize%2 else 1 # Make sure it's odd
251 kernelSize = min(config.kernelSizeMax, max(kernelSize, config.kernelSizeMin))
252
253 if metadata is not None:
254 metadata.add("ALBasisNGauss", basisNGauss)
255 metadata.add("ALBasisDegGauss", basisDegGauss)
256 metadata.add("ALBasisSigGauss", basisSigmaGauss)
257 metadata.add("ALKernelSize", kernelSize)
258 metadata.add("ALBasisMode", basisMode)
259
260 logger.debug("basisSigmaGauss: %s basisDegGauss: %s",
261 ','.join(['{:.1f}'.format(v) for v in basisSigmaGauss]),
262 ','.join(['{:d}'.format(v) for v in basisDegGauss]))
263
264 return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
265
266
267def _calculateBasisSigmas(referenceSigma, targetSigma, basisMinSigma, basisGaussBeta, basisNGauss):
268 """Calculate a series of Gaussian sigmas to use for fitting the matching
269 kernel between the reference and target images.
270
271 Parameters
272 ----------
273 referenceSigma : `float`
274 Sigma (pixel) of the reference exposure characteristic psf.
275 targetSigma : `float`
276 Sigma (pixel) of the science exposure characteristic psf.
277 basisMinSigma : `float`
278 Minimum sigma (pixels) for base Gaussians.
279 basisGaussBeta : `float`
280 Scaling multiplier of base Gaussian sigmas for automated sigma
281 determination
282 basisNGauss : `int`
283 The number of base Gaussians in the AL basis functions.
284
285 Returns
286 -------
287 basisSigmaGauss : `list` of `float`
288 Sigmas (widths) of the basis modes, in pixels.
289 """
290 kernelSigma = np.sqrt(referenceSigma**2 - targetSigma**2)
291 if kernelSigma < basisMinSigma:
292 kernelSigma = basisMinSigma
293
294 # If more than one gaussian is requested and kernelSigma is more than a
295 # factor of basisGaussBeta larger than the minimum sigma, use
296 # kernelSigma/basisGaussBeta as the size of the first gaussian.
297 if ((kernelSigma/basisGaussBeta) > basisMinSigma) & (basisNGauss > 1):
298 basisSigmaGauss = [kernelSigma/basisGaussBeta, ]
299 else:
300 basisSigmaGauss = [kernelSigma, ]
301
302 for i in range(1, basisNGauss):
303 basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
304 return basisSigmaGauss
_calculateBasisSigmas(referenceSigma, targetSigma, basisMinSigma, basisGaussBeta, basisNGauss)
generateAlardLuptonBasisList(config, targetFwhmPix=None, referenceFwhmPix=None, basisDegGauss=None, basisSigmaGauss=None, metadata=None)