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 Define the progression of Gaussians using a
143 method to derive a deconvolution sum-of-Gaussians from it's
144 convolution counterpart. [1]_ Only use 3 since the algorithm
145 assumes 3 components.
146
147 **Metadata fields**
148
149 ALBasisNGauss : `int`
150 The number of base Gaussians in the AL basis functions.
151 ALBasisDegGauss : `list` of `int`
152 Polynomial order of spatial modification of the base Gaussian functions.
153 ALBasisSigGauss : `list` of `float`
154 Sigmas in pixels of the base Gaussians.
155 ALKernelSize : `int`
156 Kernel stamp size is (ALKernelSize pix, ALKernelSize pix).
157 ALBasisMode : `str`, either of ``config``, ``convolution``, ``deconvolution``
158 Indicates whether the config file values, the convolution or deconvolution algorithm
159 was used to determine the base Gaussian sigmas and the kernel stamp size.
160
161 References
162 ----------
163 .. [1] Ulmer, W.: Inverse problem of linear combinations of Gaussian convolution kernels
164 (deconvolution) and some applications to proton/photon dosimetry and image
165 processing. http://iopscience.iop.org/0266-5611/26/8/085002 Equation 40
166 """
167
168 if config.kernelBasisSet != "alard-lupton":
169 raise RuntimeError("Cannot generate %s basis within generateAlardLuptonBasisList" %
170 config.kernelBasisSet)
171
172 kernelSize = config.kernelSize
173 fwhmScaling = config.kernelSizeFwhmScaling
174 basisNGauss = config.alardNGauss
175 basisGaussBeta = config.alardGaussBeta
176 basisMinSigma = config.alardMinSig
177 if basisDegGauss is None:
178 basisDegGauss = config.alardDegGauss
179 if basisSigmaGauss is None:
180 basisSigmaGauss = config.alardSigGauss
181
182 if len(basisDegGauss) != basisNGauss:
183 raise ValueError("len(basisDegGauss) != basisNGauss : %d vs %d" % (len(basisDegGauss), basisNGauss))
184 if len(basisSigmaGauss) != basisNGauss:
185 raise ValueError("len(basisSigmaGauss) != basisNGauss : %d vs %d" %
186 (len(basisSigmaGauss), basisNGauss))
187 if (kernelSize % 2) != 1:
188 raise ValueError("Only odd-sized Alard-Lupton bases allowed")
189
190 logger = getLogger("lsst.ip.diffim.generateAlardLuptonBasisList")
191 if (targetFwhmPix is None) or (referenceFwhmPix is None) or (not config.scaleByFwhm):
192 logger.info("PSF sigmas are not available or scaling by fwhm disabled, "
193 "falling back to config values")
194 if metadata is not None:
195 metadata.add("ALBasisNGauss", basisNGauss)
196 metadata.add("ALBasisDegGauss", basisDegGauss)
197 metadata.add("ALBasisSigGauss", basisSigmaGauss)
198 metadata.add("ALKernelSize", kernelSize)
199 metadata.add("ALBasisMode", "config")
200
201 return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
202
203 targetSigma = targetFwhmPix / sigma2fwhm
204 referenceSigma = referenceFwhmPix / sigma2fwhm
205 logger.debug("Generating matching bases for sigma %.2f pix -> %.2f pix", targetSigma, referenceSigma)
206
207
208
209
210
211
212
213
214 if targetSigma == referenceSigma:
215
216 logger.debug("Target and reference psf fwhms are equal, falling back to config values")
217 basisMode = "config"
218 elif referenceSigma > targetSigma:
219
220
221
222
223
224
225
226
227 logger.debug("Reference psf fwhm is the greater, normal convolution mode")
228 basisMode = "convolution"
229 kernelSigma = np.sqrt(referenceSigma**2 - targetSigma**2)
230 if kernelSigma < basisMinSigma:
231 kernelSigma = basisMinSigma
232
233 basisSigmaGauss = []
234 if basisNGauss == 1:
235 basisSigmaGauss.append(kernelSigma)
236 nAppended = 1
237 else:
238 if (kernelSigma/basisGaussBeta) > basisMinSigma:
239 basisSigmaGauss.append(kernelSigma/basisGaussBeta)
240 basisSigmaGauss.append(kernelSigma)
241 nAppended = 2
242 else:
243 basisSigmaGauss.append(kernelSigma)
244 nAppended = 1
245
246
247
248 for i in range(nAppended, basisNGauss):
249 basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
250
251 kernelSize = int(fwhmScaling * basisSigmaGauss[-1])
252 kernelSize += 0 if kernelSize%2 else 1
253 kernelSize =
min(config.kernelSizeMax,
max(kernelSize, config.kernelSizeMin))
254
255 else:
256
257
258
259
260
261
262
263
264 logger.debug("Target psf fwhm is the greater, deconvolution mode")
265 basisMode = "deconvolution"
266 basisNGauss = config.alardNGaussDeconv
267 basisMinSigma = config.alardMinSigDeconv
268
269 kernelSigma = np.sqrt(targetSigma**2 - referenceSigma**2)
270 if kernelSigma < basisMinSigma:
271 kernelSigma = basisMinSigma
272
273 basisSigmaGauss = []
274 if (kernelSigma/basisGaussBeta) > basisMinSigma:
275 basisSigmaGauss.append(kernelSigma/basisGaussBeta)
276 basisSigmaGauss.append(kernelSigma)
277 nAppended = 2
278 else:
279 basisSigmaGauss.append(kernelSigma)
280 nAppended = 1
281
282 for i in range(nAppended, basisNGauss):
283 basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
284
285 kernelSize = int(fwhmScaling * basisSigmaGauss[-1])
286 kernelSize += 0 if kernelSize%2 else 1
287 kernelSize =
min(config.kernelSizeMax,
max(kernelSize, config.kernelSizeMin))
288
289
290 sig0 = basisSigmaGauss[0]
291 sig1 = basisSigmaGauss[1]
292 sig2 = basisSigmaGauss[2]
293 basisSigmaGauss = []
294 for n in range(1, 3):
295 for j in range(n):
296 sigma2jn = (n - j)*sig1**2
297 sigma2jn += j * sig2**2
298 sigma2jn -= (n + 1)*sig0**2
299 sigmajn = np.sqrt(sigma2jn)
300 basisSigmaGauss.append(sigmajn)
301
302 basisSigmaGauss.sort()
303 basisNGauss = len(basisSigmaGauss)
304 basisDegGauss = [config.alardDegGaussDeconv for x in basisSigmaGauss]
305
306 if metadata is not None:
307 metadata.add("ALBasisNGauss", basisNGauss)
308 metadata.add("ALBasisDegGauss", basisDegGauss)
309 metadata.add("ALBasisSigGauss", basisSigmaGauss)
310 metadata.add("ALKernelSize", kernelSize)
311 metadata.add("ALBasisMode", basisMode)
312
313 logger.debug("basisSigmaGauss: %s basisDegGauss: %s",
314 ','.join(['{:.1f}'.format(v) for v in basisSigmaGauss]),
315 ','.join(['{:d}'.format(v) for v in basisDegGauss]))
316
317 return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)