83 basisDegGauss=
None, metadata=
None):
84 """Generate an Alard-Lupton kernel basis list based upon the Config and 85 the input FWHM of the science and template images. 89 config : `lsst.ip.diffim.PsfMatchConfigAL` 90 Configuration object for the Alard-Lupton algorithm. 91 targetFwhmPix : `float`, optional 92 Fwhm width (pixel) of the template exposure characteristic psf. 93 This is the _target_ that will be matched to the science exposure. 94 referenceFwhmPix : `float`, optional 95 Fwhm width (pixel) of the science exposure characteristic psf. 96 basisDegGauss : `list` of `int`, optional 97 Polynomial degree of each Gaussian (sigma) basis. If None, defaults to `config.alardDegGauss`. 98 metadata : `lsst.daf.base.PropertySet`, optional 99 If specified, object to collect metadata fields about the kernel basis list. 103 basisList : `list` of `lsst.afw.math.kernel.FixedKernel` 104 List of basis kernels. For each degree value ``n`` in ``config.basisDegGauss`` (n+2)(n+1)/2 kernels 105 are generated and appended to the list in the order of the polynomial parameter number. 106 See `lsst.afw.math.polynomialFunction2D` documentation for more details. 110 The polynomial functions (``f``) are always evaluated in the -1.0, +1.0 range in both x, y directions, 111 edge to edge, with ``f(0,0)`` evaluated at the kernel center pixel, ``f(-1.0,-1.0)`` at the kernel 112 ``(0,0)`` pixel. They are not scaled by the sigmas of the Gaussians. 114 Base Gaussian widths (sigmas in pixels) of the kernels are determined as: 115 - If not all fwhm parameters are provided or ``config.scaleByFwhm==False`` 116 then ``config.alardNGauss`` and ``config.alardSigGauss`` are used. 117 - If ``targetFwhmPix<referenceFwhmPix`` (normal convolution): 118 First sigma ``Sig_K`` is determined to satisfy: ``Sig_reference**2 = Sig_target**2 + Sig_K**2``. 119 If it's larger than ``config.alardMinSig * config.alardGaussBeta``, make it the 120 second kernel. Else make it the smallest kernel, unless only 1 kernel is asked for. 121 - If ``referenceFwhmPix < targetFwhmPix`` (deconvolution): 122 Define the progression of Gaussians using a 123 method to derive a deconvolution sum-of-Gaussians from it's 124 convolution counterpart. [1]_ Only use 3 since the algorithm 125 assumes 3 components. 130 .. [1] Ulmer, W.: Inverse problem of linear combinations of Gaussian convolution kernels 131 (deconvolution) and some applications to proton/photon dosimetry and image 132 processing. http://iopscience.iop.org/0266-5611/26/8/085002 Equation 40 137 - if ``config.kernelBasisSet`` is not equal to "alard-lupton" 139 - if ``config.kernelSize`` is even 140 - if the number of Gaussians and the number of given 141 sigma values are not equal or 142 - if the number of Gaussians and the number of given 143 polynomial degree values are not equal 146 if config.kernelBasisSet !=
"alard-lupton":
147 raise RuntimeError(
"Cannot generate %s basis within generateAlardLuptonBasisList" %
148 config.kernelBasisSet)
150 kernelSize = config.kernelSize
151 fwhmScaling = config.kernelSizeFwhmScaling
152 basisNGauss = config.alardNGauss
153 basisSigmaGauss = config.alardSigGauss
154 basisGaussBeta = config.alardGaussBeta
155 basisMinSigma = config.alardMinSig
156 if basisDegGauss
is None:
157 basisDegGauss = config.alardDegGauss
159 if len(basisDegGauss) != basisNGauss:
160 raise ValueError(
"len(basisDegGauss) != basisNGauss : %d vs %d" % (len(basisDegGauss), basisNGauss))
161 if len(basisSigmaGauss) != basisNGauss:
162 raise ValueError(
"len(basisSigmaGauss) != basisNGauss : %d vs %d" %
163 (len(basisSigmaGauss), basisNGauss))
164 if (kernelSize % 2) != 1:
165 raise ValueError(
"Only odd-sized Alard-Lupton bases allowed")
167 if (targetFwhmPix
is None)
or (referenceFwhmPix
is None)
or (
not config.scaleByFwhm):
168 if metadata
is not None:
169 metadata.add(
"ALBasisNGauss", basisNGauss)
170 metadata.add(
"ALBasisDegGauss", basisDegGauss)
171 metadata.add(
"ALBasisSigGauss", basisSigmaGauss)
172 metadata.add(
"ALKernelSize", kernelSize)
174 return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
176 targetSigma = targetFwhmPix / sigma2fwhm
177 referenceSigma = referenceFwhmPix / sigma2fwhm
178 logger = Log.getLogger(
"lsst.ip.diffim.generateAlardLuptonBasisList")
179 logger.debug(
"Generating matching bases for sigma %.2f pix -> %.2f pix", targetSigma, referenceSigma)
188 if targetSigma == referenceSigma:
191 elif referenceSigma > targetSigma:
200 kernelSigma = np.sqrt(referenceSigma**2 - targetSigma**2)
201 if kernelSigma < basisMinSigma:
202 kernelSigma = basisMinSigma
206 basisSigmaGauss.append(kernelSigma)
209 if (kernelSigma/basisGaussBeta) > basisMinSigma:
210 basisSigmaGauss.append(kernelSigma/basisGaussBeta)
211 basisSigmaGauss.append(kernelSigma)
214 basisSigmaGauss.append(kernelSigma)
219 for i
in range(nAppended, basisNGauss):
220 basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
222 kernelSize =
int(fwhmScaling * basisSigmaGauss[-1])
223 kernelSize += 0
if kernelSize%2
else 1
224 kernelSize =
min(config.kernelSizeMax,
max(kernelSize, config.kernelSizeMin))
235 basisNGauss = config.alardNGaussDeconv
236 basisMinSigma = config.alardMinSigDeconv
238 kernelSigma = np.sqrt(targetSigma**2 - referenceSigma**2)
239 if kernelSigma < basisMinSigma:
240 kernelSigma = basisMinSigma
243 if (kernelSigma/basisGaussBeta) > basisMinSigma:
244 basisSigmaGauss.append(kernelSigma/basisGaussBeta)
245 basisSigmaGauss.append(kernelSigma)
248 basisSigmaGauss.append(kernelSigma)
251 for i
in range(nAppended, basisNGauss):
252 basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
254 kernelSize =
int(fwhmScaling * basisSigmaGauss[-1])
255 kernelSize += 0
if kernelSize%2
else 1
256 kernelSize =
min(config.kernelSizeMax,
max(kernelSize, config.kernelSizeMin))
259 sig0 = basisSigmaGauss[0]
260 sig1 = basisSigmaGauss[1]
261 sig2 = basisSigmaGauss[2]
263 for n
in range(1, 3):
265 sigma2jn = (n - j)*sig1**2
266 sigma2jn += j * sig2**2
267 sigma2jn -= (n + 1)*sig0**2
268 sigmajn = np.sqrt(sigma2jn)
269 basisSigmaGauss.append(sigmajn)
271 basisSigmaGauss.sort()
272 basisNGauss = len(basisSigmaGauss)
273 basisDegGauss = [config.alardDegGaussDeconv
for x
in basisSigmaGauss]
275 if metadata
is not None:
276 metadata.add(
"ALBasisNGauss", basisNGauss)
277 metadata.add(
"ALBasisDegGauss", basisDegGauss)
278 metadata.add(
"ALBasisSigGauss", basisSigmaGauss)
279 metadata.add(
"ALKernelSize", kernelSize)
281 return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)