26 #include "ndarray/eigen.h"
55 -(2.0 *
bbox.getCenterY()) /
bbox.getHeight()));
63 _toChebyshevRange(makeChebyshevRangeTransform(
lsst::
geom::Box2D(
bbox))),
78 typedef detail::TrapezoidalPacker Packer;
81 void evaluateBasis1d(ndarray::Array<double, 1, 1>
const& t,
double x) {
82 int const n = t.getSize<0>();
89 for (
int i = 2; i < n; ++i) {
90 t[i] = 2.0 *
x * t[i - 1] - t[i - 2];
98 ndarray::Array<double, 2, 2> makeMatrix(ndarray::Array<double const, 1>
const&
x,
99 ndarray::Array<double const, 1>
const&
y,
101 Packer
const& packer,
Control const& ctrl) {
102 int const nPoints =
x.getSize<0>();
103 ndarray::Array<double, 1, 1> tx = ndarray::allocate(packer.nx);
104 ndarray::Array<double, 1, 1> ty = ndarray::allocate(packer.ny);
105 ndarray::Array<double, 2, 2> out = ndarray::allocate(nPoints, packer.size);
108 for (
int p = 0; p < nPoints; ++p) {
110 evaluateBasis1d(tx, sxy.getX());
111 evaluateBasis1d(ty, sxy.getY());
112 packer.pack(out[p], tx, ty);
123 Packer
const& packer,
Control const& ctrl) {
125 ndarray::Array<double, 2, 2> tx = ndarray::allocate(
bbox.getWidth(), packer.nx);
126 for (
int x =
bbox.getBeginX(), p = 0; p <
bbox.getWidth(); ++p, ++
x) {
133 ndarray::Array<double, 2, 2> out = ndarray::allocate(
bbox.getArea(), packer.size);
135 ndarray::Array<double, 1, 1> ty = ndarray::allocate(ctrl.orderY + 1);
136 for (
int y =
bbox.getBeginY(), i = 0; i <
bbox.getHeight(); ++i, ++
y) {
139 for (
int j = 0; j <
bbox.getWidth(); ++j, ++outIter) {
141 packer.pack(*outIter, tx[j], ty);
150 ndarray::Array<double const, 1>
const&
x,
151 ndarray::Array<double const, 1>
const&
y,
152 ndarray::Array<double const, 1>
const&
z,
158 Packer
const packer(ctrl);
160 ndarray::Array<double, 2, 2> matrix = makeMatrix(
x,
y,
result->_toChebyshevRange, packer, ctrl);
169 ndarray::Array<double const, 1>
const&
x,
170 ndarray::Array<double const, 1>
const&
y,
171 ndarray::Array<double const, 1>
const&
z,
172 ndarray::Array<double const, 1>
const&
w,
178 Packer
const packer(ctrl);
180 ndarray::Array<double, 2, 2> matrix = makeMatrix(
x,
y,
result->_toChebyshevRange, packer, ctrl);
183 ndarray::asEigenArray(matrix).colwise() *= ndarray::asEigenArray(
w);
184 ndarray::Array<double, 1, 1> wz = ndarray::copy(
z);
185 ndarray::asEigenArray(wz) *= ndarray::asEigenArray(
w);
193 template <
typename T>
201 Packer
const packer(ctrl);
202 ndarray::Array<double, 2, 2> matrix = makeMatrix(
bbox,
result->_toChebyshevRange, packer, ctrl);
204 ndarray::Array<double, 2, 2> imgCopy = ndarray::allocate(img.
getArray().getShape());
206 ndarray::Array<double const, 1, 1>
z = ndarray::flatten<1>(imgCopy);
220 (_coefficients.getSize<1>() - 1))
226 (_coefficients.getSize<0>() - 1))
233 ndarray::Array<double, 1, 1> packed = ndarray::allocate(packer.size);
241 return std::make_shared<ChebyshevBoundedField>(
bbox, _coefficients);
252 template <
typename CoeffGetter>
253 double evaluateFunction1d(CoeffGetter g,
double x,
int size) {
254 double b_kp2 = 0.0, b_kp1 = 0.0;
255 for (
int k = (size - 1); k > 0; --k) {
256 double b_k = g[k] + 2 *
x * b_kp1 - b_kp2;
260 return g[0] +
x * b_kp1 - b_kp2;
268 struct RecursionArrayImitator {
269 double operator[](
int i)
const {
273 RecursionArrayImitator(ndarray::Array<double const, 2, 2>
const& coefficients_,
double x_)
284 return evaluateFunction1d(RecursionArrayImitator(_coefficients, p.getX()), p.getY(),
285 _coefficients.getSize<0>());
294 return 2.0 / (1.0 - double(n * n));
300 for (ndarray::Size j = 0; j < _coefficients.getSize<0>(); j++) {
301 for (ndarray::Size i = 0; i < _coefficients.getSize<1>(); i++) {
305 return result * determinant;
314 struct PersistenceHelper {
320 PersistenceHelper(
int nx,
int ny)
322 orderX(
schema.addField<int>(
"order_x",
"maximum Chebyshev function order in x")),
325 "coefficients",
"Chebyshev function coefficients, ordered by y then x", nx * ny)) {}
327 PersistenceHelper(table::Schema
const& s)
331 class ChebyshevBoundedFieldFactory :
public table::io::PersistableFactory {
334 :
afw::table::io::PersistableFactory(
name) {}
337 CatalogVector
const& catalogs)
const override {
340 table::BaseRecord
const& record = catalogs.front().front();
341 PersistenceHelper
const keys(record.getSchema());
343 int nx = record.get(
keys.orderX) + 1;
344 int ny =
keys.coefficients.getSize() / nx;
346 ndarray::Array<double, 2, 2>
coefficients = ndarray::allocate(ny, nx);
352 std::string getChebyshevBoundedFieldPersistenceName() {
return "ChebyshevBoundedField"; }
354 ChebyshevBoundedFieldFactory registration(getChebyshevBoundedFieldPersistenceName());
359 return getChebyshevBoundedFieldPersistenceName();
365 PersistenceHelper
const keys(_coefficients.getSize<1>(), _coefficients.getSize<0>());
368 record->set(
keys.orderX, _coefficients.getSize<1>() - 1);
370 (*record)[
keys.coefficients].deep() = ndarray::flatten<1>(_coefficients);
382 if (!rhsCasted)
return false;
384 return (
getBBox() == rhsCasted->getBBox()) &&
385 (_coefficients.getShape() == rhsCasted->_coefficients.getShape()) &&
386 all(equal(_coefficients, rhsCasted->_coefficients));
389 std::string ChebyshevBoundedField::toString()
const {
391 os <<
"ChebyshevBoundedField (" << _coefficients.getShape() <<
" coefficients in y,x)";
399 #define INSTANTIATE(T) \
400 template std::shared_ptr<ChebyshevBoundedField> ChebyshevBoundedField::fit(image::Image<T> const& image, \