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
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"
33 #include "lsst/afw/math/Kernel.h"
35 #include "lsst/afw/geom.h"
36 
37 namespace pexExcept = lsst::pex::exceptions;
38 namespace afwMath = lsst::afw::math;
39 namespace afwImage = lsst::afw::image;
40 namespace afwGeom = lsst::afw::geom;
41 
43 :
44  Kernel(),
45  _kernelList(),
46  _kernelImagePtrList(),
47  _kernelSumList(),
48  _kernelParams(),
49  _isDeltaFunctionBasis(false)
50 { }
51 
53  KernelList const &kernelList,
54  std::vector<double> const &kernelParameters
55 ) :
56  Kernel(kernelList[0]->getWidth(), kernelList[0]->getHeight(), kernelList.size()),
57  _kernelList(),
58  _kernelImagePtrList(),
59  _kernelSumList(),
60  _kernelParams(kernelParameters),
61  _isDeltaFunctionBasis(false)
62 {
63  if (kernelList.size() != kernelParameters.size()) {
64  std::ostringstream os;
65  os << "kernelList.size() = " << kernelList.size()
66  << " != " << kernelParameters.size() << " = " << "kernelParameters.size()";
67  throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
68  }
69  checkKernelList(kernelList);
70  _setKernelList(kernelList);
71 }
72 
74  KernelList const &kernelList,
75  Kernel::SpatialFunction const &spatialFunction
76 ) :
77  Kernel(kernelList[0]->getWidth(), kernelList[0]->getHeight(), kernelList.size(), spatialFunction),
78  _kernelList(),
79  _kernelImagePtrList(),
80  _kernelSumList(),
81  _kernelParams(std::vector<double>(kernelList.size())),
82  _isDeltaFunctionBasis(false)
83 {
84  checkKernelList(kernelList);
85  _setKernelList(kernelList);
86 }
87 
89  KernelList const &kernelList,
90  std::vector<Kernel::SpatialFunctionPtr> const &spatialFunctionList
91 ) :
92  Kernel(kernelList[0]->getWidth(), kernelList[0]->getHeight(), spatialFunctionList),
93  _kernelList(),
94  _kernelImagePtrList(),
95  _kernelSumList(),
96  _kernelParams(std::vector<double>(kernelList.size())),
97  _isDeltaFunctionBasis(false)
98 {
99  if (kernelList.size() != spatialFunctionList.size()) {
100  std::ostringstream os;
101  os << "kernelList.size() = " << kernelList.size()
102  << " != " << spatialFunctionList.size() << " = " << "spatialFunctionList.size()";
103  throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
104  }
105  checkKernelList(kernelList);
106  _setKernelList(kernelList);
107 }
108 
109 PTR(afwMath::Kernel) afwMath::LinearCombinationKernel::clone() const {
110  PTR(Kernel) retPtr;
111  if (this->isSpatiallyVarying()) {
112  retPtr.reset(new afwMath::LinearCombinationKernel(this->_kernelList, this->_spatialFunctionList));
113  } else {
114  retPtr.reset(new afwMath::LinearCombinationKernel(this->_kernelList, this->_kernelParams));
115  }
116  retPtr->setCtr(this->getCtr());
117  return retPtr;
118 }
119 
121  if (kernelList.size() < 1) {
122  throw LSST_EXCEPT(pexExcept::InvalidParameterError, "kernelList has no elements");
123  }
124 
125  afwGeom::Extent2I const dim0 = kernelList[0]->getDimensions();
126  afwGeom::Point2I const ctr0 = kernelList[0]->getCtr();
127 
128  for (unsigned int ii = 0; ii < kernelList.size(); ++ii) {
129  if (kernelList[ii]->getDimensions() != dim0) {
130  throw LSST_EXCEPT(pexExcept::InvalidParameterError,
131  (boost::format("kernel %d has different size than kernel 0") % ii).str());
132  }
133  if (kernelList[ii]->getCtr() != ctr0) {
134  throw LSST_EXCEPT(pexExcept::InvalidParameterError,
135  (boost::format("kernel %d has different center than kernel 0") % ii).str());
136  }
137  if (kernelList[ii]->isSpatiallyVarying()) {
138  throw LSST_EXCEPT(pexExcept::InvalidParameterError,
139  (boost::format("kernel %d is spatially varying") % ii).str());
140  }
141  }
142 }
143 
145  return _kernelList;
146 }
147 
149  return _kernelSumList;
150 }
151 
153  return _kernelParams;
154 }
155 
156 PTR(afwMath::Kernel) afwMath::LinearCombinationKernel::refactor() const {
157  if (!this->isSpatiallyVarying()) {
158  return PTR(Kernel)();
159  }
160  Kernel::SpatialFunctionPtr const firstSpFuncPtr = this->_spatialFunctionList[0];
161  if (!firstSpFuncPtr->isLinearCombination()) {
162  return PTR(Kernel)();
163  }
164 
165  typedef lsst::afw::image::Image<Kernel::Pixel> KernelImage;
166  typedef boost::shared_ptr<KernelImage> KernelImagePtr;
167  typedef std::vector<KernelImagePtr> KernelImageList;
168 
169  // create kernel images for new refactored basis kernels
170  int const nSpatialParameters = this->getNSpatialParameters();
171  KernelImageList newKernelImagePtrList;
172  newKernelImagePtrList.reserve(nSpatialParameters);
173  for (int i = 0; i < nSpatialParameters; ++i) {
174  KernelImagePtr kernelImagePtr(new KernelImage(this->getDimensions()));
175  newKernelImagePtrList.push_back(kernelImagePtr);
176  }
177  KernelImage kernelImage(this->getDimensions());
178  std::vector<Kernel::SpatialFunctionPtr>::const_iterator spFuncPtrIter =
179  this->_spatialFunctionList.begin();
180  afwMath::KernelList::const_iterator kIter = _kernelList.begin();
181  afwMath::KernelList::const_iterator const kEnd = _kernelList.end();
182  for ( ; kIter != kEnd; ++kIter, ++spFuncPtrIter) {
183  if (typeid(**spFuncPtrIter) != typeid(*firstSpFuncPtr)) {
184  return PTR(Kernel)();
185  }
186 
187  (**kIter).computeImage(kernelImage, false);
188  for (int i = 0; i < nSpatialParameters; ++i) {
189  double spParam = (*spFuncPtrIter)->getParameter(i);
190  newKernelImagePtrList[i]->scaledPlus(spParam, kernelImage);
191  }
192  }
193 
194  // create new kernel; the basis kernels are fixed kernels computed above
195  // and the corresponding spatial model is the same function as the original kernel,
196  // but with all coefficients zero except coeff_i = 1.0
197  afwMath::KernelList newKernelList;
198  newKernelList.reserve(nSpatialParameters);
199  KernelImageList::iterator newKImPtrIter = newKernelImagePtrList.begin();
200  KernelImageList::iterator const newKImPtrEnd = newKernelImagePtrList.end();
201  for ( ; newKImPtrIter != newKImPtrEnd; ++newKImPtrIter) {
202  newKernelList.push_back(PTR(Kernel)(new afwMath::FixedKernel(**newKImPtrIter)));
203  }
204  std::vector<SpatialFunctionPtr> newSpFunctionPtrList;
205  for (int i = 0; i < nSpatialParameters; ++i) {
206  std::vector<double> newSpParameters(nSpatialParameters, 0.0);
207  newSpParameters[i] = 1.0;
208  SpatialFunctionPtr newSpFunctionPtr = firstSpFuncPtr->clone();
209  newSpFunctionPtr->setParameters(newSpParameters);
210  newSpFunctionPtrList.push_back(newSpFunctionPtr);
211  }
212  PTR(LinearCombinationKernel) refactoredKernel(
213  new LinearCombinationKernel(newKernelList, newSpFunctionPtrList));
214  refactoredKernel->setCtr(this->getCtr());
215  return refactoredKernel;
216 }
217 
218 std::string afwMath::LinearCombinationKernel::toString(std::string const& prefix) const {
219  std::ostringstream os;
220  os << prefix << "LinearCombinationKernel:" << std::endl;
221  os << prefix << "..Kernels:" << std::endl;
222  for (KernelList::const_iterator i = _kernelList.begin(); i != _kernelList.end(); ++i) {
223  os << (*i)->toString(prefix + "\t");
224  }
225  os << "..parameters: [ ";
226  for (std::vector<double>::const_iterator i = _kernelParams.begin(); i != _kernelParams.end(); ++i) {
227  if (i != _kernelParams.begin()) os << ", ";
228  os << *i;
229  }
230  os << " ]" << std::endl;
231  os << Kernel::toString(prefix + "\t");
232  return os.str();
233 }
234 
235 //
236 // Protected Member Functions
237 //
240  bool doNormalize
241 ) const {
242  image = 0.0;
243  double imSum = 0.0;
244  std::vector<PTR(afwImage::Image<Pixel>)>::const_iterator kImPtrIter = _kernelImagePtrList.begin();
245  std::vector<double>::const_iterator kSumIter = _kernelSumList.begin();
246  std::vector<double>::const_iterator kParIter = _kernelParams.begin();
247  for ( ; kImPtrIter != _kernelImagePtrList.end(); ++kImPtrIter, ++kSumIter, ++kParIter) {
248  image.scaledPlus(*kParIter, **kImPtrIter);
249  imSum += (*kSumIter) * (*kParIter);
250  }
251 
252  if (doNormalize) {
253  if (imSum == 0) {
254  throw LSST_EXCEPT(pexExcept::OverflowError, "Cannot normalize; kernel sum is 0");
255  }
256  image /= imSum;
257  imSum = 1;
258  }
259 
260  return imSum;
261 }
262 
263 void afwMath::LinearCombinationKernel::setKernelParameter(unsigned int ind, double value) const {
264  this->_kernelParams[ind] = value;
265 }
266 
267 //
268 // Private Member Functions
269 //
271  _kernelSumList.clear();
272  _kernelImagePtrList.clear();
273  _kernelList.clear();
274  _isDeltaFunctionBasis = true;
275  for (KernelList::const_iterator kIter = kernelList.begin(), kEnd = kernelList.end();
276  kIter != kEnd; ++kIter) {
277  PTR(Kernel) basisKernelPtr = (*kIter)->clone();
278  if (dynamic_cast<afwMath::DeltaFunctionKernel const *>(&(*basisKernelPtr)) == 0) {
279  _isDeltaFunctionBasis = false;
280  }
281  _kernelList.push_back(basisKernelPtr);
282  PTR(afwImage::Image<Pixel>) kernelImagePtr(new afwImage::Image<Pixel>(this->getDimensions()));
283  _kernelSumList.push_back(basisKernelPtr->computeImage(*kernelImagePtr, false));
284  _kernelImagePtrList.push_back(kernelImagePtr);
285  }
286 }
287 
288 // ------ Persistence ---------------------------------------------------------------------------------------
289 
290 namespace lsst { namespace afw { namespace math {
291 
292 namespace {
293 
294 struct LinearCombinationKernelPersistenceHelper : public Kernel::PersistenceHelper {
295  table::Key< table::Array<double> > amplitudes;
296  table::Key< table::Array<int> > components;
297 
298  LinearCombinationKernelPersistenceHelper(int nComponents, bool isSpatiallyVarying) :
299  Kernel::PersistenceHelper(isSpatiallyVarying ? nComponents : 0),
300  components(
301  schema.addField< table::Array<int> >("components", "archive IDs of component kernel",
302  nComponents)
303  )
304  {
305  if (!isSpatiallyVarying) {
306  amplitudes = schema.addField< table::Array<double> >("amplitudes", "amplitudes component kernel",
307  nComponents);
308  }
309  }
310 
311  LinearCombinationKernelPersistenceHelper(table::Schema const & schema_) :
312  Kernel::PersistenceHelper(schema_), components(schema["components"])
313  {
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 
324 } // anonymous
325 
327 public:
328 
330  read(InputArchive const & archive, CatalogVector const & catalogs) const {
331  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
332  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
333  LinearCombinationKernelPersistenceHelper const keys(catalogs.front().getSchema());
334  afw::table::BaseRecord const & record = catalogs.front().front();
335  geom::Extent2I dimensions(record.get(keys.dimensions));
336  std::vector<PTR(Kernel)> componentList(keys.components.getSize());
337  for (std::size_t i = 0; i < componentList.size(); ++i) {
338  componentList[i] = archive.get<Kernel>(record[keys.components[i]]);
339  }
341  if (keys.spatialFunctions.isValid()) {
342  std::vector<SpatialFunctionPtr> spatialFunctionList = keys.readSpatialFunctions(archive, record);
343  result.reset(new LinearCombinationKernel(componentList, spatialFunctionList));
344  } else {
345  std::vector<double> kernelParameters(keys.amplitudes.getSize());
346  for (std::size_t i = 0; i < kernelParameters.size(); ++i) {
347  kernelParameters[i] = record[keys.amplitudes[i]];
348  }
349  result.reset(new LinearCombinationKernel(componentList, kernelParameters));
350  }
351  LSST_ARCHIVE_ASSERT(result->getDimensions() == dimensions);
352  result->setCtr(record.get(keys.center));
353  return result;
354  }
355 
356  explicit Factory(std::string const & name) : afw::table::io::PersistableFactory(name) {}
357 };
358 
359 namespace {
360 
361 std::string getLinearCombinationKernelPersistenceName() { return "LinearCombinationKernel"; }
362 
363 LinearCombinationKernel::Factory registration(getLinearCombinationKernelPersistenceName());
364 
365 } // anonymous
366 
368  return getLinearCombinationKernelPersistenceName();
369 }
370 
372  bool isVarying = isSpatiallyVarying();
373  LinearCombinationKernelPersistenceHelper const keys(getNBasisKernels(), isVarying);
374  PTR(afw::table::BaseRecord) record = keys.write(handle, *this);
375  if (isVarying) {
376  for (int n = 0; n < keys.components.getSize(); ++n) {
377  record->set(keys.components[n], handle.put(_kernelList[n]));
378  record->set(keys.spatialFunctions[n], handle.put(_spatialFunctionList[n]));
379  }
380  } else {
381  for (int n = 0; n < keys.components.getSize(); ++n) {
382  record->set(keys.components[n], handle.put(_kernelList[n]));
383  record->set(keys.amplitudes[n], _kernelParams[n]);
384  }
385  }
386 }
387 
388 }}} // namespace lsst::afw::math
An include file to include the header files for lsst::afw::geom.
Declare the Kernel class and subclasses.
boost::shared_ptr< lsst::afw::math::Function2< double > > SpatialFunctionPtr
Definition: Kernel.h:143
#define PTR(...)
Definition: base.h:41
table::Key< std::string > name
Definition: ApCorrMap.cc:71
virtual KernelList const & getKernelList() const
Get the fixed basis kernels.
An object passed to Persistable::write to allow it to persist itself.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:47
int put(Persistable const *obj, bool permissive=false)
Save a nested Persistable to the same archive.
virtual std::string toString(std::string const &prefix="") const
Return a string representation of the kernel.
Definition: Kernel.cc:224
void checkKernelList(const KernelList &kernelList) const
Check that all kernels have the same size and center and that none are spatially varying.
A base class for factory classes used to reconstruct objects from records.
Definition: Persistable.h:231
LinearCombinationKernel()
Construct an empty LinearCombinationKernel of size 0x0.
virtual std::string getPersistenceName() const
Return the unique name used to persist this object and look up its factory.
Define a collection of useful Functions.
std::vector< double > getKernelSumList() const
Get the sum of the pixels of each fixed basis kernel.
void _setKernelList(KernelList const &kernelList)
Set _kernelList by cloning each input kernel and update the kernel image cache.
iterator begin() const
Definition: Image.cc:325
std::map< Citizen const *, CitizenInfo > table
Definition: Citizen.h:93
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
A base class for objects that can be persisted via afw::table::io Archive classes.
Definition: Persistable.h:74
table::Key< table::Array< double > > amplitudes
afw::table::PointKey< int > dimensions
A kernel that is a linear combination of fixed basis kernels.
Definition: Kernel.h:814
virtual void setKernelParameter(unsigned int ind, double value) const
Set one kernel parameter.
virtual std::vector< double > getKernelParameters() const
Return the current kernel parameters.
tbl::Schema schema
virtual double doComputeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize) const
Low-level version of computeImage.
table::Key< table::Array< int > > components
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
A vector of catalogs used by Persistable.
Definition: CatalogVector.h:26
Base class for all records.
Definition: BaseRecord.h:27
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
std::vector< boost::shared_ptr< Kernel > > KernelList
Definition: Kernel.h:542
A multi-catalog archive object used to load table::io::Persistable objects.
Definition: InputArchive.h:28
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:134
virtual std::string toString(std::string const &prefix="") const
Return a string representation of the kernel.
A class to represent a 2-dimensional array of pixels.
Definition: PSF.h:43
A kernel created from an Image.
Definition: Kernel.h:551
Include files required for standard LSST Exception handling.
void scaledPlus(double const c, Image< PixelT >const &rhs)
Add Image c*rhs to lhs.
Definition: Image.cc:677