Coverage for python / lsst / ip / diffim / makeKernelBasisList.py: 7%

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 

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