LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
LinearCombinationKernel.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 #include <numeric>
25 #include <sstream>
26 #include <stdexcept>
27 #include <typeinfo>
28 
29 #include "boost/format.hpp"
30 
31 #include "lsst/pex/exceptions.h"
32 #include "lsst/geom.h"
34 #include "lsst/afw/math/Kernel.h"
37 
39 
40 namespace lsst {
41 namespace afw {
42 
43 template std::shared_ptr<math::LinearCombinationKernel> table::io::PersistableFacade<
44  math::LinearCombinationKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
45 
46 namespace math {
47 
49  : Kernel(),
50  _kernelList(),
51  _kernelImagePtrList(),
52  _kernelSumList(),
53  _kernelParams(),
54  _isDeltaFunctionBasis(false) {}
55 
57  std::vector<double> const &kernelParameters)
58  : Kernel(kernelList[0]->getWidth(), kernelList[0]->getHeight(), kernelList.size()),
59  _kernelList(),
60  _kernelImagePtrList(),
61  _kernelSumList(),
62  _kernelParams(kernelParameters),
63  _isDeltaFunctionBasis(false) {
64  if (kernelList.size() != kernelParameters.size()) {
66  os << "kernelList.size() = " << kernelList.size() << " != " << kernelParameters.size() << " = "
67  << "kernelParameters.size()";
69  }
70  checkKernelList(kernelList);
71  _setKernelList(kernelList);
72 }
73 
75  Kernel::SpatialFunction const &spatialFunction)
76  : Kernel(kernelList[0]->getWidth(), kernelList[0]->getHeight(), kernelList.size(), spatialFunction),
77  _kernelList(),
78  _kernelImagePtrList(),
79  _kernelSumList(),
80  _kernelParams(std::vector<double>(kernelList.size())),
81  _isDeltaFunctionBasis(false) {
82  checkKernelList(kernelList);
83  _setKernelList(kernelList);
84 }
85 
87  KernelList const &kernelList, std::vector<Kernel::SpatialFunctionPtr> const &spatialFunctionList)
88  : Kernel(kernelList[0]->getWidth(), kernelList[0]->getHeight(), spatialFunctionList),
89  _kernelList(),
90  _kernelImagePtrList(),
91  _kernelSumList(),
92  _kernelParams(std::vector<double>(kernelList.size())),
93  _isDeltaFunctionBasis(false) {
94  if (kernelList.size() != spatialFunctionList.size()) {
96  os << "kernelList.size() = " << kernelList.size() << " != " << spatialFunctionList.size() << " = "
97  << "spatialFunctionList.size()";
99  }
100  checkKernelList(kernelList);
101  _setKernelList(kernelList);
102 }
103 
106  if (this->isSpatiallyVarying()) {
107  retPtr.reset(new LinearCombinationKernel(this->_kernelList, this->_spatialFunctionList));
108  } else {
109  retPtr.reset(new LinearCombinationKernel(this->_kernelList, this->_kernelParams));
110  }
111  retPtr->setCtr(this->getCtr());
112  return retPtr;
113 }
114 
116  KernelList kernelList;
117  kernelList.reserve(getKernelList().size());
118  for (const std::shared_ptr<Kernel> &kIter : getKernelList()) {
119  kernelList.push_back(kIter->resized(width, height));
120  }
121 
123  if (isSpatiallyVarying()) {
124  retPtr = std::make_shared<LinearCombinationKernel>(kernelList, _spatialFunctionList);
125  } else {
126  retPtr = std::make_shared<LinearCombinationKernel>(kernelList, _kernelParams);
127  }
128 
129  return retPtr;
130 }
131 
133  if (kernelList.size() < 1) {
134  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "kernelList has no elements");
135  }
136 
137  lsst::geom::Extent2I const dim0 = kernelList[0]->getDimensions();
138  lsst::geom::Point2I const ctr0 = kernelList[0]->getCtr();
139 
140  for (unsigned int ii = 0; ii < kernelList.size(); ++ii) {
141  if (kernelList[ii]->getDimensions() != dim0) {
143  (boost::format("kernel %d has different size than kernel 0") % ii).str());
144  }
145  if (kernelList[ii]->getCtr() != ctr0) {
147  (boost::format("kernel %d has different center than kernel 0") % ii).str());
148  }
149  if (kernelList[ii]->isSpatiallyVarying()) {
151  (boost::format("kernel %d is spatially varying") % ii).str());
152  }
153  }
154 }
155 
156 KernelList const &LinearCombinationKernel::getKernelList() const { return _kernelList; }
157 
159 
161 
163  if (!this->isSpatiallyVarying()) {
164  return std::shared_ptr<Kernel>();
165  }
166  Kernel::SpatialFunctionPtr const firstSpFuncPtr = this->_spatialFunctionList[0];
167  if (!firstSpFuncPtr->isLinearCombination()) {
168  return std::shared_ptr<Kernel>();
169  }
170 
171  typedef image::Image<Kernel::Pixel> KernelImage;
172  typedef std::shared_ptr<KernelImage> KernelImagePtr;
173  typedef std::vector<KernelImagePtr> KernelImageList;
174 
175  // create kernel images for new refactored basis kernels
176  int const nSpatialParameters = this->getNSpatialParameters();
177  KernelImageList newKernelImagePtrList;
178  newKernelImagePtrList.reserve(nSpatialParameters);
179  for (int i = 0; i < nSpatialParameters; ++i) {
180  KernelImagePtr kernelImagePtr(new KernelImage(this->getDimensions()));
181  newKernelImagePtrList.push_back(kernelImagePtr);
182  }
183  KernelImage kernelImage(this->getDimensions());
185  this->_spatialFunctionList.begin();
186  KernelList::const_iterator kIter = _kernelList.begin();
187  KernelList::const_iterator const kEnd = _kernelList.end();
188  auto &firstSpFunc = *firstSpFuncPtr;
189  auto &firstType = typeid(firstSpFunc); // noncopyable object of static storage duration
190  for (; kIter != kEnd; ++kIter, ++spFuncPtrIter) {
191  auto &spFunc = **spFuncPtrIter;
192  if (typeid(spFunc) != firstType) {
193  return std::shared_ptr<Kernel>();
194  }
195 
196  (**kIter).computeImage(kernelImage, false);
197  for (int i = 0; i < nSpatialParameters; ++i) {
198  double spParam = (*spFuncPtrIter)->getParameter(i);
199  newKernelImagePtrList[i]->scaledPlus(spParam, kernelImage);
200  }
201  }
202 
203  // create new kernel; the basis kernels are fixed kernels computed above
204  // and the corresponding spatial model is the same function as the original kernel,
205  // but with all coefficients zero except coeff_i = 1.0
206  KernelList newKernelList;
207  newKernelList.reserve(nSpatialParameters);
208  KernelImageList::iterator newKImPtrIter = newKernelImagePtrList.begin();
209  KernelImageList::iterator const newKImPtrEnd = newKernelImagePtrList.end();
210  for (; newKImPtrIter != newKImPtrEnd; ++newKImPtrIter) {
211  newKernelList.push_back(std::shared_ptr<Kernel>(new FixedKernel(**newKImPtrIter)));
212  }
213  std::vector<SpatialFunctionPtr> newSpFunctionPtrList;
214  for (int i = 0; i < nSpatialParameters; ++i) {
215  std::vector<double> newSpParameters(nSpatialParameters, 0.0);
216  newSpParameters[i] = 1.0;
217  SpatialFunctionPtr newSpFunctionPtr = firstSpFuncPtr->clone();
218  newSpFunctionPtr->setParameters(newSpParameters);
219  newSpFunctionPtrList.push_back(newSpFunctionPtr);
220  }
222  new LinearCombinationKernel(newKernelList, newSpFunctionPtrList));
223  refactoredKernel->setCtr(this->getCtr());
224  return refactoredKernel;
225 }
226 
229  os << prefix << "LinearCombinationKernel:" << std::endl;
230  os << prefix << "..Kernels:" << std::endl;
231  for (KernelList::const_iterator i = _kernelList.begin(); i != _kernelList.end(); ++i) {
232  os << (*i)->toString(prefix + "\t");
233  }
234  os << "..parameters: [ ";
235  for (std::vector<double>::const_iterator i = _kernelParams.begin(); i != _kernelParams.end(); ++i) {
236  if (i != _kernelParams.begin()) os << ", ";
237  os << *i;
238  }
239  os << " ]" << std::endl;
240  os << Kernel::toString(prefix + "\t");
241  return os.str();
242 }
243 
244 //
245 // Protected Member Functions
246 //
248  image = 0.0;
249  double imSum = 0.0;
250  std::vector<std::shared_ptr<image::Image<Pixel>>>::const_iterator kImPtrIter =
251  _kernelImagePtrList.begin();
252  std::vector<double>::const_iterator kSumIter = _kernelSumList.begin();
253  std::vector<double>::const_iterator kParIter = _kernelParams.begin();
254  for (; kImPtrIter != _kernelImagePtrList.end(); ++kImPtrIter, ++kSumIter, ++kParIter) {
255  image.scaledPlus(*kParIter, **kImPtrIter);
256  imSum += (*kSumIter) * (*kParIter);
257  }
258 
259  if (doNormalize) {
260  if (imSum == 0) {
261  throw LSST_EXCEPT(pexExcept::OverflowError, "Cannot normalize; kernel sum is 0");
262  }
263  image /= imSum;
264  imSum = 1;
265  }
266 
267  return imSum;
268 }
269 
270 void LinearCombinationKernel::setKernelParameter(unsigned int ind, double value) const {
271  this->_kernelParams[ind] = value;
272 }
273 
274 //
275 // Private Member Functions
276 //
277 void LinearCombinationKernel::_setKernelList(KernelList const &kernelList) {
278  _kernelSumList.clear();
279  _kernelImagePtrList.clear();
280  _kernelList.clear();
281  _isDeltaFunctionBasis = true;
282  for (KernelList::const_iterator kIter = kernelList.begin(), kEnd = kernelList.end(); kIter != kEnd;
283  ++kIter) {
284  std::shared_ptr<Kernel> basisKernelPtr = (*kIter)->clone();
285  if (dynamic_cast<DeltaFunctionKernel const *>(&(*basisKernelPtr)) == 0) {
286  _isDeltaFunctionBasis = false;
287  }
288  _kernelList.push_back(basisKernelPtr);
290  _kernelSumList.push_back(basisKernelPtr->computeImage(*kernelImagePtr, false));
291  _kernelImagePtrList.push_back(kernelImagePtr);
292  }
293 }
294 
295 // ------ Persistence ---------------------------------------------------------------------------------------
296 
297 namespace {
298 
299 struct LinearCombinationKernelPersistenceHelper : public Kernel::PersistenceHelper {
300  table::Key<table::Array<double>> amplitudes;
301  table::Key<table::Array<int>> components;
302 
303  LinearCombinationKernelPersistenceHelper(int nComponents, bool isSpatiallyVarying)
304  : Kernel::PersistenceHelper(isSpatiallyVarying ? nComponents : 0),
305  components(schema.addField<table::Array<int>>("components", "archive IDs of component kernel",
306  nComponents)) {
307  if (!isSpatiallyVarying) {
308  amplitudes = schema.addField<table::Array<double>>("amplitudes", "amplitudes component kernel",
309  nComponents);
310  }
311  }
312 
313  LinearCombinationKernelPersistenceHelper(table::Schema const &schema_)
314  : Kernel::PersistenceHelper(schema_), components(schema["components"]) {
315  if (!spatialFunctions.isValid()) {
316  amplitudes = schema["amplitudes"];
317  LSST_ARCHIVE_ASSERT(amplitudes.getSize() == components.getSize());
318  } else {
319  LSST_ARCHIVE_ASSERT(spatialFunctions.getSize() == components.getSize());
320  }
321  }
322 };
323 
324 } // namespace
325 
327 public:
329  CatalogVector const &catalogs) const override {
330  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
331  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
332  LinearCombinationKernelPersistenceHelper const keys(catalogs.front().getSchema());
333  afw::table::BaseRecord const &record = catalogs.front().front();
334  lsst::geom::Extent2I dimensions(record.get(keys.dimensions));
335  std::vector<std::shared_ptr<Kernel>> componentList(keys.components.getSize());
336  for (std::size_t i = 0; i < componentList.size(); ++i) {
337  componentList[i] = archive.get<Kernel>(record[keys.components[i]]);
338  }
340  if (keys.spatialFunctions.isValid()) {
341  std::vector<SpatialFunctionPtr> spatialFunctionList = keys.readSpatialFunctions(archive, record);
342  result.reset(new LinearCombinationKernel(componentList, spatialFunctionList));
343  } else {
344  std::vector<double> kernelParameters(keys.amplitudes.getSize());
345  for (std::size_t i = 0; i < kernelParameters.size(); ++i) {
346  kernelParameters[i] = record[keys.amplitudes[i]];
347  }
348  result.reset(new LinearCombinationKernel(componentList, kernelParameters));
349  }
350  LSST_ARCHIVE_ASSERT(result->getDimensions() == dimensions);
351  result->setCtr(record.get(keys.center));
352  return result;
353  }
354 
355  explicit Factory(std::string const &name) : afw::table::io::PersistableFactory(name) {}
356 };
357 
358 namespace {
359 
360 std::string getLinearCombinationKernelPersistenceName() { return "LinearCombinationKernel"; }
361 
362 LinearCombinationKernel::Factory registration(getLinearCombinationKernelPersistenceName());
363 
364 } // namespace
365 
367  return getLinearCombinationKernelPersistenceName();
368 }
369 
371  bool isVarying = isSpatiallyVarying();
372  LinearCombinationKernelPersistenceHelper const keys(getNBasisKernels(), isVarying);
373  std::shared_ptr<afw::table::BaseRecord> record = keys.write(handle, *this);
374  if (isVarying) {
375  for (int n = 0; n < keys.components.getSize(); ++n) {
376  record->set(keys.components[n], handle.put(_kernelList[n]));
377  record->set(keys.spatialFunctions[n], handle.put(_spatialFunctionList[n]));
378  }
379  } else {
380  for (int n = 0; n < keys.components.getSize(); ++n) {
381  record->set(keys.components[n], handle.put(_kernelList[n]));
382  record->set(keys.amplitudes[n], _kernelParams[n]);
383  }
384  }
385 }
386 } // namespace math
387 } // namespace afw
388 } // namespace lsst
py::object result
Definition: _schema.cc:430
table::Key< std::string > name
Definition: Amplifier.cc:116
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
table::Key< table::Array< int > > components
table::Key< table::Array< double > > amplitudes
std::ostream * os
Definition: Schema.cc:746
std::string prefix
Definition: SchemaMapper.cc:79
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
table::Schema schema
Definition: python.h:134
T begin(T... args)
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
A kernel that has only one non-zero pixel (of value 1)
Definition: Kernel.h:644
A kernel created from an Image.
Definition: Kernel.h:472
A Function taking two arguments.
Definition: Function.h:259
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:111
lsst::geom::Extent2I const getDimensions() const
Return the Kernel's dimensions (width, height)
Definition: Kernel.h:213
std::vector< SpatialFunctionPtr > _spatialFunctionList
Definition: Kernel.h:450
lsst::geom::Point2I getCtr() const
Return index of kernel's center.
Definition: Kernel.h:235
int getNSpatialParameters() const
Return the number of spatial parameters (0 if not spatially varying)
Definition: Kernel.h:252
virtual std::string toString(std::string const &prefix="") const
Return a string representation of the kernel.
Definition: Kernel.cc:195
bool isSpatiallyVarying() const
Return true iff the kernel is spatially varying (has a spatial function)
Definition: Kernel.h:334
std::shared_ptr< afw::table::io::Persistable > read(InputArchive const &archive, CatalogVector const &catalogs) const override
Construct a new object from the given InputArchive and vector of catalogs.
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
std::shared_ptr< Kernel > resized(int width, int height) const override
Return a pointer to a clone with specified kernel dimensions.
std::vector< double > getKernelSumList() const
Get the sum of the pixels of each fixed basis kernel.
std::shared_ptr< Kernel > refactor() const
Refactor the kernel as a linear combination of N bases where N is the number of parameters for the sp...
virtual KernelList const & getKernelList() const
Get the fixed basis kernels.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
LinearCombinationKernel()
Construct an empty LinearCombinationKernel of size 0x0.
std::shared_ptr< Kernel > clone() const override
Return a pointer to a deep copy of this kernel.
void setKernelParameter(unsigned int ind, double value) const override
Set one kernel parameter.
void checkKernelList(const KernelList &kernelList) const
Check that all kernels have the same size and center and that none are spatially varying.
double doComputeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize) const override
Low-level version of computeImage.
int getNBasisKernels() const
Get the number of basis kernels.
Definition: Kernel.h:769
std::vector< double > getKernelParameters() const override
Return the current kernel parameters.
std::string toString(std::string const &prefix="") const override
Return a string representation of the kernel.
Base class for all records.
Definition: BaseRecord.h:31
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
A vector of catalogs used by Persistable.
Definition: CatalogVector.h:29
A multi-catalog archive object used to load table::io::Persistable objects.
Definition: InputArchive.h:31
std::shared_ptr< Persistable > get(int id) const
Load the Persistable with the given ID and return it.
An object passed to Persistable::write to allow it to persist itself.
int put(Persistable const *obj, bool permissive=false)
Save an object to the archive and return a unique ID that can be used to retrieve it from an InputArc...
A base class for factory classes used to reconstruct objects from records.
Definition: Persistable.h:228
PersistableFactory(std::string const &name)
Constructor for the factory.
Definition: Persistable.cc:74
Reports invalid arguments.
Definition: Runtime.h:66
Reports when the result of an arithmetic operation is too large for the destination type.
Definition: Runtime.h:124
T clear(T... args)
T end(T... args)
T endl(T... args)
T front(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
A base class for image defects.
STL namespace.
T push_back(T... args)
T reserve(T... args)
T reset(T... args)
T size(T... args)