LSST Applications g04c3c9f7ca+2075667efa,g1e125bf412+5f448d5fcf,g2079a07aa2+3e9fd84d81,g2305ad1205+b635cf1488,g2bbee38e9b+6c6beb4891,g337abbeb29+6c6beb4891,g33d1c0ed96+6c6beb4891,g3a166c0a6a+6c6beb4891,g50ff169b8f+96c6868917,g52b1c1532d+585e252eca,g591dd9f2cf+42f171e1e6,g5c3423f6d4+d536b04327,g607f77f49a+d536b04327,g6f43f06aed+ca1339dc19,g858d7b2824+d536b04327,g8ee334c5b4+d7f9608c2f,g9963eaa53e+b3dc1655d3,g998f4353bf+d536b04327,g99cad8db69+8ef2408349,g9ddcbc5298+9a081db1e4,ga1e77700b3+2cbb763275,gadfd92a7e4+aec2f3b930,gae0086650b+585e252eca,gb0e22166c9+0e73c8378f,gb3b7280ab2+cb5fdb229e,gbb8dafda3b+a327199e22,gc120e1dc64+88074880ea,gc28159a63d+6c6beb4891,gcdd4ae20e8+bd241b2308,gcde1bda545+903e937d91,gcf0d15dbbd+bd241b2308,gdaeeff99f8+f9a426f77a,gddc38dedce+585e252eca,ge79ae78c31+6c6beb4891,gfbcc870c63+b310236976,w.2024.23
LSST Data Management Base Package
Loading...
Searching...
No Matches
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"
31
33
34namespace lsst {
35namespace afw {
36
37template std::shared_ptr<math::SeparableKernel> table::io::PersistableFacade<
38 math::SeparableKernel>::dynamicCast(std::shared_ptr<table::io::Persistable> const&);
39
40namespace 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) {
53}
54
55SeparableKernel::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) {
69}
70
71SeparableKernel::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
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
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
185void 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
198double 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
268namespace {
272void _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
319int SeparableKernel::getCacheSize() const { return _kernelColCache.size(); };
320} // namespace math
321} // namespace afw
322} // namespace lsst
#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
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 class to represent a 2-dimensional array of pixels.
Definition Image.h:51
A Function taking one argument.
Definition Function.h:202
virtual std::shared_ptr< Function1< ReturnT > > clone() const =0
Return a pointer to a deep copy of this function.
std::string toString(std::string const &prefix="") const override
Return a string representation of the function.
Definition Function.h:240
void setParameter(unsigned int ind, double newValue)
Set one function parameter without range checking.
Definition Function.h:143
std::vector< double > const & getParameters() const noexcept
Return all function parameters.
Definition Function.h:129
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 reset(T... args)
T size(T... args)
table::Key< int > cacheSize