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
SeparableKernel.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 <vector>
25 #include <iostream>
26 
27 #include "lsst/pex/exceptions.h"
28 #include "lsst/afw/math/Kernel.h"
31 
33 
34 namespace lsst {
35 namespace afw {
36 
37 template std::shared_ptr<math::SeparableKernel> table::io::PersistableFacade<
38  math::SeparableKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const&);
39 
40 namespace math {
41 
43  : Kernel(),
44  _kernelColFunctionPtr(),
45  _kernelRowFunctionPtr(),
46  _localColList(0),
47  _localRowList(0),
48  _kernelX(0),
49  _kernelY(0),
50  _kernelRowCache(0),
51  _kernelColCache(0) {
52  _setKernelXY();
53 }
54 
55 SeparableKernel::SeparableKernel(int width, int height, KernelFunction const& kernelColFunction,
56  KernelFunction const& kernelRowFunction,
57  Kernel::SpatialFunction const& spatialFunction)
58  : Kernel(width, height, kernelColFunction.getNParameters() + kernelRowFunction.getNParameters(),
59  spatialFunction),
60  _kernelColFunctionPtr(kernelColFunction.clone()),
61  _kernelRowFunctionPtr(kernelRowFunction.clone()),
62  _localColList(width),
63  _localRowList(height),
64  _kernelX(width),
65  _kernelY(height),
66  _kernelRowCache(0),
67  _kernelColCache(0) {
68  _setKernelXY();
69 }
70 
71 SeparableKernel::SeparableKernel(int width, int height, KernelFunction const& kernelColFunction,
72  KernelFunction const& kernelRowFunction,
73  std::vector<Kernel::SpatialFunctionPtr> const& spatialFunctionList)
74  : Kernel(width, height, spatialFunctionList),
75  _kernelColFunctionPtr(kernelColFunction.clone()),
76  _kernelRowFunctionPtr(kernelRowFunction.clone()),
77  _localColList(width),
78  _localRowList(height),
79  _kernelX(width),
80  _kernelY(height),
81  _kernelRowCache(0),
82  _kernelColCache(0) {
83  if (kernelColFunction.getNParameters() + kernelRowFunction.getNParameters() !=
84  spatialFunctionList.size()) {
86  os << "kernelColFunction.getNParameters() + kernelRowFunction.getNParameters() = "
87  << kernelColFunction.getNParameters() << " + " << kernelRowFunction.getNParameters()
88  << " != " << spatialFunctionList.size() << " = "
89  << "spatialFunctionList.size()";
91  }
92 
93  _setKernelXY();
94 }
95 
98  if (this->isSpatiallyVarying()) {
99  retPtr.reset(new SeparableKernel(this->getWidth(), this->getHeight(), *(this->_kernelColFunctionPtr),
100  *(this->_kernelRowFunctionPtr), this->_spatialFunctionList));
101  } else {
102  retPtr.reset(new SeparableKernel(this->getWidth(), this->getHeight(), *(this->_kernelColFunctionPtr),
103  *(this->_kernelRowFunctionPtr)));
104  }
105  retPtr->setCtr(this->getCtr());
106  retPtr->computeCache(this->getCacheSize());
107  return retPtr;
108 }
109 
110 std::shared_ptr<Kernel> SeparableKernel::resized(int width, int height) const {
112  if (isSpatiallyVarying()) {
113  retPtr = std::make_shared<SeparableKernel>(width, height, *_kernelColFunctionPtr,
114  *_kernelRowFunctionPtr, _spatialFunctionList);
115  } else {
116  retPtr = std::make_shared<SeparableKernel>(width, height, *_kernelColFunctionPtr,
117  *_kernelRowFunctionPtr);
118  }
119  return retPtr;
120 }
121 
123  bool doNormalize, double x, double y) const {
124  if (static_cast<int>(colList.size()) != this->getWidth() ||
125  static_cast<int>(rowList.size()) != this->getHeight()) {
127  os << "colList.size(), rowList.size() = (" << colList.size() << ", " << rowList.size() << ") != ("
128  << this->getWidth() << ", " << this->getHeight() << ") = "
129  << "kernel dimensions";
131  }
132  if (this->isSpatiallyVarying()) {
134  }
135 
136  return basicComputeVectors(colList, rowList, doNormalize);
137 }
138 
140  return _kernelColFunctionPtr->clone();
141 }
142 
144  return _kernelRowFunctionPtr->clone();
145 }
146 
149  os << prefix << "SeparableKernel:" << std::endl;
150  os << prefix
151  << "..x (width) function: " << (_kernelColFunctionPtr ? _kernelColFunctionPtr->toString() : "None")
152  << std::endl;
153  os << prefix
154  << "..y (rows) function: " << (_kernelRowFunctionPtr ? _kernelRowFunctionPtr->toString() : "None")
155  << std::endl;
156  os << Kernel::toString(prefix + "\t");
157  return os.str();
158 }
159 
161  std::vector<double> allParams = _kernelColFunctionPtr->getParameters();
162  std::vector<double> yParams = _kernelRowFunctionPtr->getParameters();
163  std::copy(yParams.begin(), yParams.end(), std::back_inserter(allParams));
164  return allParams;
165 }
166 
167 //
168 // Protected Member Functions
169 //
170 
172  double imSum = basicComputeVectors(_localColList, _localRowList, doNormalize);
173 
174  for (int y = 0; y != image.getHeight(); ++y) {
175  image::Image<Pixel>::x_iterator imPtr = image.row_begin(y);
176  for (std::vector<Pixel>::iterator colIter = _localColList.begin(); colIter != _localColList.end();
177  ++colIter, ++imPtr) {
178  *imPtr = (*colIter) * _localRowList[y];
179  }
180  }
181 
182  return imSum;
183 }
184 
185 void SeparableKernel::setKernelParameter(unsigned int ind, double value) const {
186  unsigned int const nColParams = _kernelColFunctionPtr->getNParameters();
187  if (ind < nColParams) {
188  _kernelColFunctionPtr->setParameter(ind, value);
189  } else {
190  _kernelRowFunctionPtr->setParameter(ind - nColParams, value);
191  }
192 }
193 
194 //
195 // Private Member Functions
196 //
197 
198 double SeparableKernel::basicComputeVectors(std::vector<Pixel>& colList, std::vector<Pixel>& rowList,
199  bool doNormalize) const {
200  double colSum = 0.0;
201  if (_kernelColCache.empty()) {
202  for (unsigned int i = 0; i != colList.size(); ++i) {
203  double colFuncValue = (*_kernelColFunctionPtr)(_kernelX[i]);
204  colList[i] = colFuncValue;
205  colSum += colFuncValue;
206  }
207  } else {
208  int const cacheSize = _kernelColCache.size();
209 
210  int const indx = this->getKernelParameter(0) * cacheSize;
211 
212  std::vector<double>& cachedValues = _kernelColCache.at(indx);
213  for (unsigned int i = 0; i != colList.size(); ++i) {
214  double colFuncValue = cachedValues[i];
215  colList[i] = colFuncValue;
216  colSum += colFuncValue;
217  }
218  }
219 
220  double rowSum = 0.0;
221  if (_kernelRowCache.empty()) {
222  for (unsigned int i = 0; i != rowList.size(); ++i) {
223  double rowFuncValue = (*_kernelRowFunctionPtr)(_kernelY[i]);
224  rowList[i] = rowFuncValue;
225  rowSum += rowFuncValue;
226  }
227  } else {
228  int const cacheSize = _kernelRowCache.size();
229 
230  int const indx = this->getKernelParameter(1) * cacheSize;
231 
232  std::vector<double>& cachedValues = _kernelRowCache.at(indx);
233  for (unsigned int i = 0; i != rowList.size(); ++i) {
234  double rowFuncValue = cachedValues[i];
235  rowList[i] = rowFuncValue;
236  rowSum += rowFuncValue;
237 
238 #if 0
239  if (indx == cacheSize/2) {
240  if (::fabs(rowFuncValue - (*_kernelRowFunctionPtr)(_kernelX[i])) > 1e-2) {
241  std::cout << indx << " " << i << " "
242  << rowFuncValue << " "
243  << (*_kernelRowFunctionPtr)(_kernelX[i])
244  << std::endl;
245  }
246  }
247 #endif
248  }
249  }
250 
251  double imSum = colSum * rowSum;
252  if (doNormalize) {
253  if ((colSum == 0) || (rowSum == 0)) {
254  throw LSST_EXCEPT(pexExcept::OverflowError, "Cannot normalize; kernel sum is 0");
255  }
256  for (double & colIter : colList) {
257  colIter /= colSum;
258  }
259 
260  for (double & rowIter : rowList) {
261  rowIter /= rowSum;
262  }
263  imSum = 1.0;
264  }
265  return imSum;
266 }
267 
268 namespace {
272 void _computeCache(int const cacheSize, std::vector<double> const& x,
274  if (cacheSize <= 0) {
275  kernelCache->erase(kernelCache->begin(), kernelCache->end());
276  return;
277  }
278 
279  if (kernelCache[0].size() != x.size()) { // invalid
280  kernelCache->erase(kernelCache->begin(), kernelCache->end());
281  }
282 
283  int const old_cacheSize = kernelCache->size();
284 
285  if (cacheSize == old_cacheSize) {
286  return; // nothing to do
287  }
288 
289  if (cacheSize < old_cacheSize) {
290  kernelCache->erase(kernelCache->begin() + cacheSize, kernelCache->end());
291  } else {
292  kernelCache->resize(cacheSize);
293  for (int i = old_cacheSize; i != cacheSize; ++i) {
294  (*kernelCache)[i].resize(x.size());
295  }
296  }
297  //
298  // Actually fill the cache
299  //
300  for (int i = 0; i != cacheSize; ++i) {
301  func->setParameter(0, (i + 0.5) / static_cast<double>(cacheSize));
302  for (unsigned int j = 0; j != x.size(); ++j) {
303  (*kernelCache)[i][j] = (*func)(x[j]);
304  }
305  }
306 }
307 } // namespace
308 
311 
312  func = getKernelColFunction();
313  _computeCache(cacheSize, _kernelY, func, &_kernelColCache);
314 
315  func = getKernelRowFunction();
316  _computeCache(cacheSize, _kernelX, func, &_kernelRowCache);
317 }
318 
319 int SeparableKernel::getCacheSize() const { return _kernelColCache.size(); };
320 } // namespace math
321 } // namespace afw
322 } // namespace lsst
double x
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
std::ostream * os
Definition: Schema.cc:557
std::string prefix
Definition: SchemaMapper.cc:72
int y
Definition: SpanSet.cc:48
T at(T... args)
T back_inserter(T... args)
T begin(T... args)
typename _view_t::x_iterator x_iterator
An iterator for traversing the pixels in a row.
Definition: ImageBase.h:133
A Function taking one argument.
Definition: Function.h:202
std::string toString(std::string const &prefix="") const override
Return a string representation of the function.
Definition: Function.h:240
virtual std::shared_ptr< Function1< ReturnT > > clone() const =0
Return a pointer to a deep copy of this function.
A Function taking two arguments.
Definition: Function.h:259
std::vector< double > const & getParameters() const noexcept
Return all function parameters.
Definition: Function.h:129
void setParameter(unsigned int ind, double newValue)
Set one function parameter without range checking.
Definition: Function.h:143
unsigned int getNParameters() const noexcept
Return the number of function parameters.
Definition: Function.h:112
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:110
std::vector< SpatialFunctionPtr > _spatialFunctionList
Definition: Kernel.h:449
int getHeight() const
Return the Kernel's height.
Definition: Kernel.h:229
lsst::geom::Point2I getCtr() const
Return index of kernel's center.
Definition: Kernel.h:234
virtual std::string toString(std::string const &prefix="") const
Return a string representation of the kernel.
Definition: Kernel.cc:185
int getWidth() const
Return the Kernel's width.
Definition: Kernel.h:224
bool isSpatiallyVarying() const
Return true iff the kernel is spatially varying (has a spatial function)
Definition: Kernel.h:333
void setKernelParametersFromSpatialModel(double x, double y) const
Set the kernel parameters from the spatial model (if any).
Definition: Kernel.cc:217
KernelFunctionPtr getKernelRowFunction() const
Get a deep copy of the row kernel function.
std::string toString(std::string const &prefix="") const override
Return a string representation of the kernel.
std::shared_ptr< Kernel > clone() const override
Return a pointer to a deep copy of this kernel.
std::vector< double > getKernelParameters() const override
Return the current kernel parameters.
std::shared_ptr< KernelFunction > KernelFunctionPtr
Definition: Kernel.h:863
double getKernelParameter(unsigned int i) const override
Return a particular Kernel Parameter (no bounds checking).
Definition: Kernel.h:933
int getCacheSize() const override
Get the current cache size (0 if none)
double doComputeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize) const override
Low-level version of computeImage.
KernelFunctionPtr getKernelColFunction() const
Get a deep copy of the col kernel function.
double computeVectors(std::vector< Pixel > &colList, std::vector< Pixel > &rowList, bool doNormalize, double x=0.0, double y=0.0) const
Compute the column and row arrays in place, where kernel(col, row) = colList(col) * rowList(row)
std::shared_ptr< Kernel > resized(int width, int height) const override
Return a pointer to a clone with specified kernel dimensions.
virtual void _setKernelXY() override
Definition: Kernel.h:1006
void setKernelParameter(unsigned int ind, double value) const override
Set one kernel parameter.
SeparableKernel()
Construct an empty spatially invariant SeparableKernel of size 0x0.
void computeCache(int const cacheSize) override
Compute a cache of Kernel values, if desired.
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 copy(T... args)
T empty(T... args)
T end(T... args)
T endl(T... args)
T fabs(T... args)
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
A base class for image defects.
T reset(T... args)
T size(T... args)
table::Key< int > cacheSize