LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
makeKernelBasisList.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 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 import diffimLib
23 import lsst.pex.logging as pexLog
24 import numpy as np
25 sigma2fwhm = 2. * np.sqrt(2. * np.log(2.))
26 
27 def makeKernelBasisList(config, targetFwhmPix=None, referenceFwhmPix=None,
28  basisDegGauss=None, metadata=None):
29  """Generate the appropriate Kernel basis based on the Config"""
30  if config.kernelBasisSet == "alard-lupton":
31  return generateAlardLuptonBasisList(config, targetFwhmPix=targetFwhmPix,
32  referenceFwhmPix=referenceFwhmPix,
33  basisDegGauss=basisDegGauss,
34  metadata=metadata)
35  elif config.kernelBasisSet == "delta-function":
36  kernelSize = config.kernelSize
37  return diffimLib.makeDeltaFunctionBasisList(kernelSize, kernelSize)
38  else:
39  raise ValueError("Cannot generate %s basis set" % (config.kernelBasisSet))
40 
41 def generateAlardLuptonBasisList(config, targetFwhmPix=None, referenceFwhmPix=None,
42  basisDegGauss=None, metadata=None):
43  """Generate an Alard-Lupton kernel basis based upon the Config and
44  the input FWHM of the science and template images"""
45 
46  if config.kernelBasisSet != "alard-lupton":
47  raise RuntimeError("Cannot generate %s basis within generateAlardLuptonBasisList" % config.kernelBasisSet)
48 
49  kernelSize = config.kernelSize
50  fwhmScaling = config.kernelSizeFwhmScaling
51  basisNGauss = config.alardNGauss
52  basisSigmaGauss = config.alardSigGauss
53  basisGaussBeta = config.alardGaussBeta
54  basisMinSigma = config.alardMinSig
55  if basisDegGauss is None:
56  basisDegGauss = config.alardDegGauss
57 
58  if len(basisDegGauss) != basisNGauss:
59  raise ValueError("len(basisDegGauss) != basisNGauss : %d vs %d" % (len(basisDegGauss), basisNGauss))
60  if len(basisSigmaGauss) != basisNGauss:
61  raise ValueError("len(basisSigmaGauss) != basisNGauss : %d vs %d" % (len(basisSigmaGauss), basisNGauss))
62  if (kernelSize % 2) != 1:
63  raise ValueError("Only odd-sized Alard-Lupton bases allowed")
64 
65  if (targetFwhmPix is None) or (referenceFwhmPix is None) or (not config.scaleByFwhm):
66  if metadata is not None:
67  metadata.add("ALBasisNGauss", basisNGauss)
68  metadata.add("ALBasisDegGauss", basisDegGauss)
69  metadata.add("ALBasisSigGauss", basisSigmaGauss)
70  metadata.add("ALKernelSize", kernelSize)
71 
72  return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
73 
74  targetSigma = targetFwhmPix / sigma2fwhm
75  referenceSigma = referenceFwhmPix / sigma2fwhm
76  pexLog.Trace("lsst.ip.diffim.generateAlardLuptonBasisList", 2,
77  "Generating matching bases for sigma %.2f pix -> %.2f pix" % (targetSigma, referenceSigma))
78 
79  # Modify the size of Alard Lupton kernels based upon the images FWHM
80  #
81  # Note the operation is : template.x.kernel = science
82  #
83  # Assuming the template and science image Psfs are Gaussians with
84  # the Fwhm above, Fwhm_T **2 + Fwhm_K **2 = Fwhm_S **2
85  #
86  if targetSigma == referenceSigma:
87  # Leave defaults as-is
88  pass
89  elif referenceSigma > targetSigma:
90  # Normal convolution
91 
92  # First Gaussian has the sigma that comes from the convolution
93  # of two Gaussians : Sig_S**2 = Sig_T**2 + Sig_K**2
94  #
95  # If it's larger than basisMinSigma * basisGaussBeta, make it the
96  # second kernel. Else make it the smallest kernel. Unless
97  # only 1 kernel is asked for.
98  kernelSigma = np.sqrt(referenceSigma**2 - targetSigma**2)
99  if kernelSigma < basisMinSigma:
100  kernelSigma = basisMinSigma
101 
102  basisSigmaGauss = []
103  if basisNGauss == 1:
104  basisSigmaGauss.append(kernelSigma)
105  nAppended = 1
106  else:
107  if (kernelSigma/basisGaussBeta) > basisMinSigma:
108  basisSigmaGauss.append(kernelSigma/basisGaussBeta)
109  basisSigmaGauss.append(kernelSigma)
110  nAppended = 2
111  else:
112  basisSigmaGauss.append(kernelSigma)
113  nAppended = 1
114 
115  # Any other Gaussians above basisNGauss=1 come from a scaling
116  # relationship: Sig_i+1 / Sig_i = basisGaussBeta
117  for i in range(nAppended,basisNGauss):
118  basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
119 
120  kernelSize = int(fwhmScaling * basisSigmaGauss[-1])
121  kernelSize += 0 if kernelSize%2 else 1 # Make sure it's odd
122  kernelSize = min(config.kernelSizeMax, max(kernelSize, config.kernelSizeMin))
123 
124  else:
125  # Deconvolution; Define the progression of Gaussians using a
126  # method to derive a deconvolution sum-of-Gaussians from it's
127  # convolution counterpart. Only use 3 since the algorithm
128  # assumes 3 components.
129  #
130  # http://iopscience.iop.org/0266-5611/26/8/085002 Equation 40
131 
132  # Use specializations for deconvolution
133  basisNGauss = config.alardNGaussDeconv
134  basisMinSigma = config.alardMinSigDeconv
135 
136  kernelSigma = np.sqrt(targetSigma**2 - referenceSigma**2)
137  if kernelSigma < basisMinSigma:
138  kernelSigma = basisMinSigma
139 
140  basisSigmaGauss = []
141  if (kernelSigma/basisGaussBeta) > basisMinSigma:
142  basisSigmaGauss.append(kernelSigma/basisGaussBeta)
143  basisSigmaGauss.append(kernelSigma)
144  nAppended = 2
145  else:
146  basisSigmaGauss.append(kernelSigma)
147  nAppended = 1
148 
149  for i in range(nAppended,basisNGauss):
150  basisSigmaGauss.append(basisSigmaGauss[-1]*basisGaussBeta)
151 
152  kernelSize = int(fwhmScaling * basisSigmaGauss[-1])
153  kernelSize += 0 if kernelSize%2 else 1 # Make sure it's odd
154  kernelSize = min(config.kernelSizeMax, max(kernelSize, config.kernelSizeMin))
155 
156  # Now build a deconvolution set from these sigmas
157  sig0 = basisSigmaGauss[0]
158  sig1 = basisSigmaGauss[1]
159  sig2 = basisSigmaGauss[2]
160  basisSigmaGauss = []
161  for n in range(1, 3):
162  for j in range(n):
163  sigma2jn = (n - j) * sig1**2
164  sigma2jn += j * sig2**2
165  sigma2jn -= (n + 1) * sig0**2
166  sigmajn = np.sqrt(sigma2jn)
167  basisSigmaGauss.append(sigmajn)
168 
169  basisSigmaGauss.sort()
170  basisNGauss = len(basisSigmaGauss)
171  basisDegGauss = [config.alardDegGaussDeconv for x in basisSigmaGauss]
172 
173  if metadata is not None:
174  metadata.add("ALBasisNGauss", basisNGauss)
175  metadata.add("ALBasisDegGauss", basisDegGauss)
176  metadata.add("ALBasisSigGauss", basisSigmaGauss)
177  metadata.add("ALKernelSize", kernelSize)
178 
179  return diffimLib.makeAlardLuptonBasisList(kernelSize//2, basisNGauss, basisSigmaGauss, basisDegGauss)
180 
limited backward compatibility to the DC2 run-time trace facilities
Definition: Trace.h:93