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
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 <algorithm>
25 #include <iterator>
26 #include <sstream>
27 
28 #include "lsst/pex/exceptions.h"
29 #include "lsst/afw/math/Kernel.h"
30 
31 namespace pexExcept = lsst::pex::exceptions;
32 namespace afwImage = lsst::afw::image;
33 namespace afwMath = lsst::afw::math;
34 
36 :
37  Kernel(),
38  _kernelColFunctionPtr(),
39  _kernelRowFunctionPtr(),
40  _localColList(0), _localRowList(0),
41  _kernelX(0), _kernelY(0),
42  _kernelRowCache(0), _kernelColCache(0)
43 {
44  _setKernelXY();
45 }
46 
48  int width,
49  int height,
50  KernelFunction const& kernelColFunction,
51  KernelFunction const& kernelRowFunction,
52  Kernel::SpatialFunction const& spatialFunction
53 ) :
54  Kernel(width, height, kernelColFunction.getNParameters() + kernelRowFunction.getNParameters(),
55  spatialFunction),
56  _kernelColFunctionPtr(kernelColFunction.clone()),
57  _kernelRowFunctionPtr(kernelRowFunction.clone()),
58  _localColList(width), _localRowList(height),
59  _kernelX(width), _kernelY(height),
60  _kernelRowCache(0), _kernelColCache(0)
61 {
62  _setKernelXY();
63 }
64 
66  int width,
67  int height,
68  KernelFunction const& kernelColFunction,
69  KernelFunction const& kernelRowFunction,
70  std::vector<Kernel::SpatialFunctionPtr> const& spatialFunctionList
71 ) :
72  Kernel(width, height, spatialFunctionList),
73  _kernelColFunctionPtr(kernelColFunction.clone()),
74  _kernelRowFunctionPtr(kernelRowFunction.clone()),
75  _localColList(width), _localRowList(height),
76  _kernelX(width), _kernelY(height),
77  _kernelRowCache(0), _kernelColCache(0)
78 {
79  if (kernelColFunction.getNParameters() + kernelRowFunction.getNParameters()
80  != spatialFunctionList.size()) {
81  std::ostringstream os;
82  os << "kernelColFunction.getNParameters() + kernelRowFunction.getNParameters() = "
83  << kernelColFunction.getNParameters() << " + " << kernelRowFunction.getNParameters()
84  << " != " << spatialFunctionList.size() << " = " << "spatialFunctionList.size()";
85  throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
86  }
87 
88  _setKernelXY();
89 }
90 
91 PTR(afwMath::Kernel) afwMath::SeparableKernel::clone() const {
92  PTR(afwMath::Kernel) retPtr;
93  if (this->isSpatiallyVarying()) {
94  retPtr.reset(new afwMath::SeparableKernel(this->getWidth(), this->getHeight(),
95  *(this->_kernelColFunctionPtr), *(this->_kernelRowFunctionPtr), this->_spatialFunctionList));
96  } else {
97  retPtr.reset(new afwMath::SeparableKernel(this->getWidth(), this->getHeight(),
98  *(this->_kernelColFunctionPtr), *(this->_kernelRowFunctionPtr)));
99  }
100  retPtr->setCtr(this->getCtr());
101  retPtr->computeCache(this->getCacheSize());
102  return retPtr;
103 }
104 
106  std::vector<Pixel> &colList,
107  std::vector<Pixel> &rowList,
108  bool doNormalize,
109  double x,
110  double y
111 ) const {
112  if (static_cast<int>(colList.size()) != this->getWidth()
113  || static_cast<int>(rowList.size()) != this->getHeight()) {
114  std::ostringstream os;
115  os << "colList.size(), rowList.size() = ("
116  << colList.size() << ", " << rowList.size()
117  << ") != ("<< this->getWidth() << ", " << this->getHeight()
118  << ") = " << "kernel dimensions";
119  throw LSST_EXCEPT(pexExcept::InvalidParameterError, os.str());
120  }
121  if (this->isSpatiallyVarying()) {
122  this->setKernelParametersFromSpatialModel(x, y);
123  }
124 
125  return basicComputeVectors(colList, rowList, doNormalize);
126 }
127 
129 ) const {
130  return _kernelColFunctionPtr->clone();
131 }
132 
134 ) const {
135  return _kernelRowFunctionPtr->clone();
136 }
137 
138 std::string afwMath::SeparableKernel::toString(std::string const& prefix) const {
139  std::ostringstream os;
140  os << prefix << "SeparableKernel:" << std::endl;
141  os << prefix << "..x (width) function: "
142  << (_kernelColFunctionPtr ? _kernelColFunctionPtr->toString() : "None") << std::endl;
143  os << prefix << "..y (rows) function: "
144  << (_kernelRowFunctionPtr ? _kernelRowFunctionPtr->toString() : "None") << std::endl;
145  os << Kernel::toString(prefix + "\t");
146  return os.str();
147 }
148 
150  std::vector<double> allParams = _kernelColFunctionPtr->getParameters();
151  std::vector<double> yParams = _kernelRowFunctionPtr->getParameters();
152  std::copy(yParams.begin(), yParams.end(), std::back_inserter(allParams));
153  return allParams;
154 }
155 
156 //
157 // Protected Member Functions
158 //
159 
162  bool doNormalize
163 ) const {
164  double imSum = basicComputeVectors(_localColList, _localRowList, doNormalize);
165 
166  for (int y = 0; y != image.getHeight(); ++y) {
168  for (std::vector<Pixel>::iterator colIter = _localColList.begin();
169  colIter != _localColList.end(); ++colIter, ++imPtr) {
170  *imPtr = (*colIter)*_localRowList[y];
171  }
172  }
173 
174  return imSum;
175 }
176 
177 void afwMath::SeparableKernel::setKernelParameter(unsigned int ind, double value) const {
178  unsigned int const nColParams = _kernelColFunctionPtr->getNParameters();
179  if (ind < nColParams) {
180  _kernelColFunctionPtr->setParameter(ind, value);
181  } else {
182  _kernelRowFunctionPtr->setParameter(ind - nColParams, value);
183  }
184 }
185 
186 //
187 // Private Member Functions
188 //
189 
191  std::vector<Pixel> &colList,
192  std::vector<Pixel> &rowList,
193  bool doNormalize
194 ) const {
195  double colSum = 0.0;
196  if (_kernelColCache.empty()) {
197  for (unsigned int i = 0; i != colList.size(); ++i) {
198  double colFuncValue = (*_kernelColFunctionPtr)(_kernelX[i]);
199  colList[i] = colFuncValue;
200  colSum += colFuncValue;
201  }
202  } else {
203  int const cacheSize = _kernelColCache.size();
204 
205  int const indx = this->getKernelParameter(0)*cacheSize;
206 
207  std::vector<double> &cachedValues = _kernelColCache.at(indx);
208  for (unsigned int i = 0; i != colList.size(); ++i) {
209  double colFuncValue = cachedValues[i];
210  colList[i] = colFuncValue;
211  colSum += colFuncValue;
212  }
213  }
214 
215  double rowSum = 0.0;
216  if (_kernelRowCache.empty()) {
217  for (unsigned int i = 0; i != rowList.size(); ++i) {
218  double rowFuncValue = (*_kernelRowFunctionPtr)(_kernelY[i]);
219  rowList[i] = rowFuncValue;
220  rowSum += rowFuncValue;
221  }
222  } else {
223  int const cacheSize = _kernelRowCache.size();
224 
225  int const indx = this->getKernelParameter(1)*cacheSize;
226 
227  std::vector<double> &cachedValues = _kernelRowCache.at(indx);
228  for (unsigned int i = 0; i != rowList.size(); ++i) {
229  double rowFuncValue = cachedValues[i];
230  rowList[i] = rowFuncValue;
231  rowSum += rowFuncValue;
232 
233 #if 0
234  if (indx == cacheSize/2) {
235  if (::fabs(rowFuncValue - (*_kernelRowFunctionPtr)(_kernelX[i])) > 1e-2) {
236  std::cout << indx << " " << i << " "
237  << rowFuncValue << " "
238  << (*_kernelRowFunctionPtr)(_kernelX[i])
239  << std::endl;
240  }
241  }
242 #endif
243  }
244  }
245 
246  double imSum = colSum * rowSum;
247  if (doNormalize) {
248  if ((colSum == 0) || (rowSum == 0)) {
249  throw LSST_EXCEPT(pexExcept::OverflowError, "Cannot normalize; kernel sum is 0");
250  }
251  for (std::vector<Pixel>::iterator colIter = colList.begin(); colIter != colList.end(); ++colIter) {
252  *colIter /= colSum;
253  }
254 
255  for (std::vector<Pixel>::iterator rowIter = rowList.begin(); rowIter != rowList.end(); ++rowIter) {
256  *rowIter /= rowSum;
257  }
258  imSum = 1.0;
259  }
260  return imSum;
261 }
262 
263 /************************************************************************************************************/
264 namespace {
268  void _computeCache(int const cacheSize,
269  std::vector<double> const& x,
271  std::vector<std::vector<double> > *kernelCache)
272  {
273  if (cacheSize <= 0) {
274  kernelCache->erase(kernelCache->begin(), kernelCache->end());
275  return;
276  }
277 
278  if (kernelCache[0].size() != x.size()) { // invalid
279  kernelCache->erase(kernelCache->begin(), kernelCache->end());
280  }
281 
282  int const old_cacheSize = kernelCache->size();
283 
284  if (cacheSize == old_cacheSize) {
285  return; // nothing to do
286  }
287 
288  if (cacheSize < old_cacheSize) {
289  kernelCache->erase(kernelCache->begin() + cacheSize, kernelCache->end());
290  } else {
291  kernelCache->resize(cacheSize);
292  for (int i = old_cacheSize; i != cacheSize; ++i) {
293  (*kernelCache)[i].resize(x.size());
294  }
295  }
296  //
297  // Actually fill the cache
298  //
299  for (int i = 0; i != cacheSize; ++i) {
300  func->setParameter(0, (i + 0.5)/static_cast<double>(cacheSize));
301  for (unsigned int j = 0; j != x.size(); ++j) {
302  (*kernelCache)[i][j] = (*func)(x[j]);
303  }
304  }
305  }
306 }
307 
309  int const cacheSize
310 ) {
312 
313  func = getKernelColFunction();
314  _computeCache(cacheSize, _kernelY, func, &_kernelColCache);
315 
316  func = getKernelRowFunction();
317  _computeCache(cacheSize, _kernelX, func, &_kernelRowCache);
318 }
319 
321  return _kernelColCache.size();
322 };
int y
virtual int getCacheSize() const
Get the current cache size (0 if none)
Declare the Kernel class and subclasses.
#define PTR(...)
Definition: base.h:41
x_iterator row_begin(int y) const
Definition: Image.h:319
virtual std::string toString(std::string const &prefix="") const
Return a string representation of the kernel.
Definition: Kernel.cc:224
virtual std::vector< double > getKernelParameters() const
Return the current kernel parameters.
SeparableKernel()
Construct an empty spatially invariant SeparableKernel of size 0x0.
SelectEigenView< T >::Type copy(Eigen::EigenBase< T > const &other)
Copy an arbitrary Eigen expression into a new EigenView.
Definition: eigen.h:390
A kernel described by a pair of functions: func(x, y) = colFunc(x) * rowFunc(y)
Definition: Kernel.h:986
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) ...
tbl::Key< int > cacheSize
Definition: CoaddPsf.cc:318
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
unsigned int getNParameters() const
Return the number of function parameters.
Definition: Function.h:125
virtual void computeCache(int const cacheSize)
Compute a cache of Kernel values, if desired.
virtual void setKernelParameter(unsigned int ind, double value) const
Set one kernel parameter.
A Function taking one argument.
Definition: Function.h:229
int x
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
KernelFunctionPtr getKernelColFunction() const
Get a deep copy of the col kernel function.
KernelFunctionPtr getKernelRowFunction() const
Get a deep copy of the row kernel function.
virtual void _setKernelXY()
Definition: Kernel.h:1138
virtual double doComputeImage(lsst::afw::image::Image< Pixel > &image, bool doNormalize) const
Low-level version of computeImage.
virtual std::string toString(std::string const &prefix="") const
Return a string representation of the kernel.
double basicComputeVectors(std::vector< Pixel > &colList, std::vector< Pixel > &rowList, bool doNormalize) const
Compute the column and row arrays in place, where kernel(col, row) = colList(col) * rowList(row) ...
int getHeight() const
Return the number of rows in the image.
Definition: Image.h:239
boost::shared_ptr< KernelFunction > KernelFunctionPtr
Definition: Kernel.h:991
Kernels are used for convolution with MaskedImages and (eventually) Images.
Definition: Kernel.h:134
Include files required for standard LSST Exception handling.