LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
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 <memory>
25 #include <vector>
26 #include <sstream>
27 #include <stdexcept>
28 #include <typeinfo>
29 
30 #include "boost/format.hpp"
31 
32 #include "lsst/pex/exceptions.h"
33 #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  using KernelImage = image::Image<Kernel::Pixel>;
172  using KernelImagePtr = std::shared_ptr<KernelImage>;
173  using KernelImageList = std::vector<KernelImagePtr>;
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 (auto const &i : _kernelList) {
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 (auto const &kIter : kernelList) {
283  std::shared_ptr<Kernel> basisKernelPtr = kIter->clone();
284  if (dynamic_cast<DeltaFunctionKernel const *>(&(*basisKernelPtr)) == nullptr) {
285  _isDeltaFunctionBasis = false;
286  }
287  _kernelList.push_back(basisKernelPtr);
289  _kernelSumList.push_back(basisKernelPtr->computeImage(*kernelImagePtr, false));
290  _kernelImagePtrList.push_back(kernelImagePtr);
291  }
292 }
293 
294 // ------ Persistence ---------------------------------------------------------------------------------------
295 
296 namespace {
297 
298 struct LinearCombinationKernelPersistenceHelper : public Kernel::PersistenceHelper {
299  table::Key<table::Array<double>> amplitudes;
300  table::Key<table::Array<int>> components;
301 
302  LinearCombinationKernelPersistenceHelper(int nComponents, bool isSpatiallyVarying)
303  : Kernel::PersistenceHelper(isSpatiallyVarying ? nComponents : 0),
304  components(schema.addField<table::Array<int>>("components", "archive IDs of component kernel",
305  nComponents)) {
306  if (!isSpatiallyVarying) {
307  amplitudes = schema.addField<table::Array<double>>("amplitudes", "amplitudes component kernel",
308  nComponents);
309  }
310  }
311 
312  LinearCombinationKernelPersistenceHelper(table::Schema const &schema_)
313  : Kernel::PersistenceHelper(schema_), components(schema["components"]) {
314  if (!spatialFunctions.isValid()) {
315  amplitudes = schema["amplitudes"];
316  LSST_ARCHIVE_ASSERT(amplitudes.getSize() == components.getSize());
317  } else {
318  LSST_ARCHIVE_ASSERT(spatialFunctions.getSize() == components.getSize());
319  }
320  }
321 };
322 
323 } // namespace
324 
326 public:
328  CatalogVector const &catalogs) const override {
329  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
330  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
331  LinearCombinationKernelPersistenceHelper const keys(catalogs.front().getSchema());
332  afw::table::BaseRecord const &record = catalogs.front().front();
333  lsst::geom::Extent2I dimensions(record.get(keys.dimensions));
334  std::vector<std::shared_ptr<Kernel>> componentList(keys.components.getSize());
335  for (std::size_t i = 0; i < componentList.size(); ++i) {
336  componentList[i] = archive.get<Kernel>(record[keys.components[i]]);
337  }
339  if (keys.spatialFunctions.isValid()) {
340  std::vector<SpatialFunctionPtr> spatialFunctionList = keys.readSpatialFunctions(archive, record);
341  result.reset(new LinearCombinationKernel(componentList, spatialFunctionList));
342  } else {
343  std::vector<double> kernelParameters(keys.amplitudes.getSize());
344  for (std::size_t i = 0; i < kernelParameters.size(); ++i) {
345  kernelParameters[i] = record[keys.amplitudes[i]];
346  }
347  result.reset(new LinearCombinationKernel(componentList, kernelParameters));
348  }
349  LSST_ARCHIVE_ASSERT(result->getDimensions() == dimensions);
350  result->setCtr(record.get(keys.center));
351  return result;
352  }
353 
354  explicit Factory(std::string const &name) : afw::table::io::PersistableFactory(name) {}
355 };
356 
357 namespace {
358 
359 std::string getLinearCombinationKernelPersistenceName() { return "LinearCombinationKernel"; }
360 
361 LinearCombinationKernel::Factory registration(getLinearCombinationKernelPersistenceName());
362 
363 } // namespace
364 
366  return getLinearCombinationKernelPersistenceName();
367 }
368 
370  bool isVarying = isSpatiallyVarying();
371  LinearCombinationKernelPersistenceHelper const keys(getNBasisKernels(), isVarying);
372  std::shared_ptr<afw::table::BaseRecord> record = keys.write(handle, *this);
373  if (isVarying) {
374  for (std::size_t n = 0; n < keys.components.getSize(); ++n) {
375  record->set(keys.components[n], handle.put(_kernelList[n]));
376  record->set(keys.spatialFunctions[n], handle.put(_spatialFunctionList[n]));
377  }
378  } else {
379  for (std::size_t n = 0; n < keys.components.getSize(); ++n) {
380  record->set(keys.components[n], handle.put(_kernelList[n]));
381  record->set(keys.amplitudes[n], _kernelParams[n]);
382  }
383  }
384 }
385 } // namespace math
386 } // namespace afw
387 } // namespace lsst
py::object result
Definition: _schema.cc:429
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:48
table::Key< table::Array< int > > components
table::Key< table::Array< double > > amplitudes
std::ostream * os
Definition: Schema.cc:557
std::string prefix
Definition: SchemaMapper.cc:72
#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:51
A kernel that has only one non-zero pixel (of value 1)
Definition: Kernel.h:643
A kernel created from an Image.
Definition: Kernel.h:471
A Function taking two arguments.
Definition: Function.h:259
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:110
lsst::geom::Extent2I const getDimensions() const
Return the Kernel's dimensions (width, height)
Definition: Kernel.h:212
std::vector< SpatialFunctionPtr > _spatialFunctionList
Definition: Kernel.h:449
lsst::geom::Point2I getCtr() const
Return index of kernel's center.
Definition: Kernel.h:234
int getNSpatialParameters() const
Return the number of spatial parameters (0 if not spatially varying)
Definition: Kernel.h:251
virtual std::string toString(std::string const &prefix="") const
Return a string representation of the kernel.
Definition: Kernel.cc:185
bool isSpatiallyVarying() const
Return true iff the kernel is spatially varying (has a spatial function)
Definition: Kernel.h:333
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:768
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)