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.
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`.
106 If specified, object to collect metadata fields about the kernel basis list.
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.
118 -
if ``config.kernelBasisSet``
is not equal to
"alard-lupton"
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
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.
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.
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.
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.
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
168 if config.kernelBasisSet !=
"alard-lupton":
169 raise RuntimeError(
"Cannot generate %s basis within generateAlardLuptonBasisList" %
170 config.kernelBasisSet)
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
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")
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")
201 return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
203 targetSigma = targetFwhmPix / sigma2fwhm
204 referenceSigma = referenceFwhmPix / sigma2fwhm
205 logger.debug(
"Generating matching bases for sigma %.2f pix -> %.2f pix", targetSigma, referenceSigma)
214 if targetSigma == referenceSigma:
216 logger.debug(
"Target and reference psf fwhms are equal, falling back to config values")
218 elif referenceSigma > targetSigma:
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
235 basisSigmaGauss.append(kernelSigma)
238 if (kernelSigma/basisGaussBeta) > basisMinSigma:
239 basisSigmaGauss.append(kernelSigma/basisGaussBeta)
240 basisSigmaGauss.append(kernelSigma)
243 basisSigmaGauss.append(kernelSigma)
248 for i
in range(nAppended, basisNGauss):
249 basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
251 kernelSize =
int(fwhmScaling * basisSigmaGauss[-1])
252 kernelSize += 0
if kernelSize%2
else 1
253 kernelSize =
min(config.kernelSizeMax,
max(kernelSize, config.kernelSizeMin))
264 logger.debug(
"Target psf fwhm is the greater, deconvolution mode")
265 basisMode =
"deconvolution"
266 basisNGauss = config.alardNGaussDeconv
267 basisMinSigma = config.alardMinSigDeconv
269 kernelSigma = np.sqrt(targetSigma**2 - referenceSigma**2)
270 if kernelSigma < basisMinSigma:
271 kernelSigma = basisMinSigma
274 if (kernelSigma/basisGaussBeta) > basisMinSigma:
275 basisSigmaGauss.append(kernelSigma/basisGaussBeta)
276 basisSigmaGauss.append(kernelSigma)
279 basisSigmaGauss.append(kernelSigma)
282 for i
in range(nAppended, basisNGauss):
283 basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
285 kernelSize =
int(fwhmScaling * basisSigmaGauss[-1])
286 kernelSize += 0
if kernelSize%2
else 1
287 kernelSize =
min(config.kernelSizeMax,
max(kernelSize, config.kernelSizeMin))
290 sig0 = basisSigmaGauss[0]
291 sig1 = basisSigmaGauss[1]
292 sig2 = basisSigmaGauss[2]
294 for n
in range(1, 3):
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)
302 basisSigmaGauss.sort()
303 basisNGauss = len(basisSigmaGauss)
304 basisDegGauss = [config.alardDegGaussDeconv
for x
in basisSigmaGauss]
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)
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]))
317 return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
Class for storing generic metadata.
def getLogger(loggername)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)