LSST Applications g0b6bd0c080+a72a5dd7e6,g1182afd7b4+2a019aa3bb,g17e5ecfddb+2b8207f7de,g1d67935e3f+06cf436103,g38293774b4+ac198e9f13,g396055baef+6a2097e274,g3b44f30a73+6611e0205b,g480783c3b1+98f8679e14,g48ccf36440+89c08d0516,g4b93dc025c+98f8679e14,g5c4744a4d9+a302e8c7f0,g613e996a0d+e1c447f2e0,g6c8d09e9e7+25247a063c,g7271f0639c+98f8679e14,g7a9cd813b8+124095ede6,g9d27549199+a302e8c7f0,ga1cf026fa3+ac198e9f13,ga32aa97882+7403ac30ac,ga786bb30fb+7a139211af,gaa63f70f4e+9994eb9896,gabf319e997+ade567573c,gba47b54d5d+94dc90c3ea,gbec6a3398f+06cf436103,gc6308e37c7+07dd123edb,gc655b1545f+ade567573c,gcc9029db3c+ab229f5caf,gd01420fc67+06cf436103,gd877ba84e5+06cf436103,gdb4cecd868+6f279b5b48,ge2d134c3d5+cc4dbb2e3f,ge448b5faa6+86d1ceac1d,gecc7e12556+98f8679e14,gf3ee170dca+25247a063c,gf4ac96e456+ade567573c,gf9f5ea5b4d+ac198e9f13,gff490e6085+8c2580be5c,w.2022.27
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"
37
39
40namespace lsst {
41namespace afw {
42
43template std::shared_ptr<math::LinearCombinationKernel> table::io::PersistableFacade<
44 math::LinearCombinationKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const &);
45
46namespace 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
156KernelList const &LinearCombinationKernel::getKernelList() const { return _kernelList; }
157
159
161
163 if (!this->isSpatiallyVarying()) {
165 }
166 Kernel::SpatialFunctionPtr const firstSpFuncPtr = this->_spatialFunctionList[0];
167 if (!firstSpFuncPtr->isLinearCombination()) {
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) {
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
270void LinearCombinationKernel::setKernelParameter(unsigned int ind, double value) const {
271 this->_kernelParams[ind] = value;
272}
273
274//
275// Private Member Functions
276//
277void 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
296namespace {
297
298struct 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
326public:
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
357namespace {
358
359std::string getLinearCombinationKernelPersistenceName() { return "LinearCombinationKernel"; }
360
361LinearCombinationKernel::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)