LSSTApplications  20.0.0
LSSTDataManagementBasePackage
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
schema
table::Schema schema
Definition: Amplifier.cc:115
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
std::string
STL class.
std::shared_ptr
STL class.
lsst::afw::math::Kernel::getNSpatialParameters
int getNSpatialParameters() const
Return the number of spatial parameters (0 if not spatially varying)
Definition: Kernel.h:274
lsst::afw::math::Kernel::getDimensions
lsst::geom::Extent2I const getDimensions() const
Return the Kernel's dimensions (width, height)
Definition: Kernel.h:213
lsst::afw::math::LinearCombinationKernel::LinearCombinationKernel
LinearCombinationKernel()
Construct an empty LinearCombinationKernel of size 0x0.
Definition: LinearCombinationKernel.cc:48
lsst::afw::math::LinearCombinationKernel::Factory::Factory
Factory(std::string const &name)
Definition: LinearCombinationKernel.cc:355
lsst::afw::math::Kernel::isSpatiallyVarying
bool isSpatiallyVarying() const
Return true iff the kernel is spatially varying (has a spatial function)
Definition: Kernel.h:380
lsst::afw::table::BaseRecord::get
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
lsst::afw::math::DeltaFunctionKernel
A kernel that has only one non-zero pixel (of value 1)
Definition: Kernel.h:690
lsst::afw::math::LinearCombinationKernel::getKernelParameters
std::vector< double > getKernelParameters() const override
Return the current kernel parameters.
Definition: LinearCombinationKernel.cc:160
std::vector::reserve
T reserve(T... args)
lsst::afw::table::io::OutputArchiveHandle
An object passed to Persistable::write to allow it to persist itself.
Definition: OutputArchive.h:118
lsst::afw::math::Function2::clone
virtual std::shared_ptr< Function2< ReturnT > > clone() const =0
Return a pointer to a deep copy of this function.
std::vector< std::shared_ptr< Kernel > >
std::vector::size
T size(T... args)
Persistable.cc
pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
lsst::afw::math::FixedKernel
A kernel created from an Image.
Definition: Kernel.h:518
FunctionLibrary.h
lsst::afw
Definition: imageAlgorithm.dox:1
lsst::afw::table::io::InputArchive
A multi-catalog archive object used to load table::io::Persistable objects.
Definition: InputArchive.h:31
astshim.keyMap.keyMapContinued.keys
def keys(self)
Definition: keyMapContinued.py:6
geom.h
lsst::afw::geom.transform.transformContinued.name
string name
Definition: transformContinued.py:32
std::vector::front
T front(T... args)
std::shared_ptr::reset
T reset(T... args)
lsst::afw::math::Kernel::getCtr
lsst::geom::Point2I getCtr() const
Return index of kernel's center.
Definition: Kernel.h:235
std::vector::clear
T clear(T... args)
lsst::afw::math::LinearCombinationKernel::Factory::read
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.
Definition: LinearCombinationKernel.cc:328
lsst::afw::table::io::PersistableFactory::PersistableFactory
PersistableFactory(std::string const &name)
Constructor for the factory.
Definition: Persistable.cc:74
std::vector::push_back
T push_back(T... args)
LSST_ARCHIVE_ASSERT
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
lsst::afw::math::Function::isLinearCombination
virtual bool isLinearCombination() const noexcept
Is the function a linear combination of its parameters?
Definition: Function.h:138
lsst::afw::math::LinearCombinationKernel::refactor
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...
Definition: LinearCombinationKernel.cc:162
lsst::afw::math::Kernel::toString
virtual std::string toString(std::string const &prefix="") const
Return a string representation of the kernel.
Definition: Kernel.cc:195
lsst::afw::math::LinearCombinationKernel::resized
std::shared_ptr< Kernel > resized(int width, int height) const override
Return a pointer to a clone with specified kernel dimensions.
Definition: LinearCombinationKernel.cc:115
lsst::afw::math::LinearCombinationKernel::toString
std::string toString(std::string const &prefix="") const override
Return a string representation of the kernel.
Definition: LinearCombinationKernel.cc:227
lsst::afw::math::LinearCombinationKernel::setKernelParameter
void setKernelParameter(unsigned int ind, double value) const override
Set one kernel parameter.
Definition: LinearCombinationKernel.cc:270
lsst::afw::math::LinearCombinationKernel::getNBasisKernels
int getNBasisKernels() const
Get the number of basis kernels.
Definition: Kernel.h:815
lsst::pex::exceptions::OverflowError
Reports when the result of an arithmetic operation is too large for the destination type.
Definition: Runtime.h:124
lsst::afw::table::BaseRecord
Base class for all records.
Definition: BaseRecord.h:31
lsst::afw::math::Function::setParameters
void setParameters(std::vector< double > const &params)
Set all function parameters.
Definition: Function.h:156
lsst::afw::table::io::CatalogVector
A vector of catalogs used by Persistable.
Definition: CatalogVector.h:29
lsst::afw::table::io::PersistableFactory
A base class for factory classes used to reconstruct objects from records.
Definition: Persistable.h:228
lsst::afw::math::LinearCombinationKernel::write
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
Definition: LinearCombinationKernel.cc:370
dimensions
afw::table::PointKey< int > dimensions
Definition: GaussianPsf.cc:49
result
py::object result
Definition: _schema.cc:429
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
lsst::afw::math::LinearCombinationKernel::Factory
Definition: LinearCombinationKernel.cc:326
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::ostringstream
STL class.
amplitudes
table::Key< table::Array< double > > amplitudes
Definition: LinearCombinationKernel.cc:300
lsst::afw::math::Kernel::_spatialFunctionList
std::vector< SpatialFunctionPtr > _spatialFunctionList
Definition: Kernel.h:496
lsst::afw::math::LinearCombinationKernel::clone
std::shared_ptr< Kernel > clone() const override
Return a pointer to a deep copy of this kernel.
Definition: LinearCombinationKernel.cc:104
lsst::afw::math::Function2< double >
os
std::ostream * os
Definition: Schema.cc:746
lsst::afw::math::LinearCombinationKernel::getKernelList
virtual KernelList const & getKernelList() const
Get the fixed basis kernels.
Definition: LinearCombinationKernel.cc:156
std::endl
T endl(T... args)
lsst::pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
KernelPersistenceHelper.h
lsst::afw::image.slicing.Factory
Factory
Definition: slicing.py:252
std::vector::begin
T begin(T... args)
std
STL namespace.
Kernel.h
lsst::geom::Point< int, 2 >
lsst::afw::math::LinearCombinationKernel::checkKernelList
void checkKernelList(const KernelList &kernelList) const
Check that all kernels have the same size and center and that none are spatially varying.
Definition: LinearCombinationKernel.cc:132
lsst::afw::math::LinearCombinationKernel::doComputeImage
double doComputeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize) const override
Low-level version of computeImage.
Definition: LinearCombinationKernel.cc:247
lsst::afw::math::LinearCombinationKernel::getKernelSumList
std::vector< double > getKernelSumList() const
Get the sum of the pixels of each fixed basis kernel.
Definition: LinearCombinationKernel.cc:158
lsst::pex::exceptions
Definition: Exception.h:37
std::size_t
lsst::afw::math::Kernel
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:111
lsst::afw::table::io::OutputArchiveHandle::put
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...
Definition: OutputArchive.cc:216
std::vector::end
T end(T... args)
lsst::afw::math::LinearCombinationKernel::getPersistenceName
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
Definition: LinearCombinationKernel.cc:366
lsst::afw::image::Image
A class to represent a 2-dimensional array of pixels.
Definition: Image.h:58
lsst::geom::Extent< int, 2 >
components
table::Key< table::Array< int > > components
Definition: LinearCombinationKernel.cc:301
prefix
std::string prefix
Definition: SchemaMapper.cc:79
exceptions.h
lsst::afw::table::io::InputArchive::get
std::shared_ptr< Persistable > get(int id) const
Load the Persistable with the given ID and return it.
Definition: InputArchive.cc:182