LSST Applications g063fba187b+66a50001ff,g0f08755f38+1a22dc2551,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g217e2c1bcf+12e87a5bd8,g246886dfd9+466c7b6c06,g28da252d5a+858b171e04,g2bbee38e9b+404b60ec9b,g2bc492864f+404b60ec9b,g3156d2b45e+6e55a43351,g347aa1857d+404b60ec9b,g35bb328faa+a8ce1bb630,g3a166c0a6a+404b60ec9b,g3e281a1b8c+c5dd892a6c,g414038480c+6b9177ef31,g41af890bb2+70bea58702,g599934f4f4+b8c5400ca5,g781aacb6e4+a8ce1bb630,g7af13505b9+b5b9cefdb8,g80478fca09+c2997882f3,g82479be7b0+8974e6af0f,g858d7b2824+1a22dc2551,g89c8672015+f4add4ffd5,g8f1c07a47a+de51c9b0a5,g9125e01d80+a8ce1bb630,ga5288a1d22+b66f8cf76b,gb58c049af0+d64f4d3760,gc28159a63d+404b60ec9b,gcab2d0539d+66cf1de5d4,gcf0d15dbbd+12cb7e2563,gda6a2b7d83+12cb7e2563,gdaeeff99f8+1711a396fd,ge79ae78c31+404b60ec9b,gef2f8181fd+414189b318,gf0baf85859+c1f95f4921,gf0c06eb49c+1a22dc2551,gfa517265be+1a22dc2551,gfa999e8aa5+17cd334064,v28.0.0.rc2
LSST Data Management Base Package
Loading...
Searching...
No Matches
TransmissionCurve.cc
Go to the documentation of this file.
1/*
2 * LSST Data Management System
3 * Copyright 2017 LSST/AURA.
4 *
5 * This product includes software developed by the
6 * LSST Project (http://www.lsst.org/).
7 *
8 * This program is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the LSST License Statement and
19 * the GNU General Public License along with this program. If not,
20 * see <http://www.lsstcorp.org/LegalNotices/>.
21 */
22
23#include <algorithm>
24#include <memory>
25
26#include "ndarray.h"
27
28#include "gsl/gsl_interp.h"
29#include "gsl/gsl_interp2d.h"
30#include "gsl/gsl_errno.h"
31
36
37using namespace std::string_literals;
38
39#define LSST_CHECK_GSL(type, status) \
40 if (status) throw LSST_EXCEPT(type, "GSL error: "s + ::gsl_strerror(status))
41
42namespace lsst {
43namespace afw {
44
45template std::shared_ptr<image::TransmissionCurve> table::io::PersistableFacade<
47
48namespace image {
49
50namespace {
51
52/*
53 * The TransmissionCurve implementation returned by TransmissionCurve::makeIdentity.
54 *
55 * This is zero-state singleton whose throughput is always exactly 1.
56 */
57class IdentityTransmissionCurve : public TransmissionCurve {
58public:
59 static constexpr char const* NAME = "IdentityTransmissionCurve";
60
62 static std::shared_ptr<IdentityTransmissionCurve> instance(new IdentityTransmissionCurve());
63 return instance;
64 }
65
66 std::pair<double, double> getWavelengthBounds() const override {
67 constexpr double inf = std::numeric_limits<double>::infinity();
68 return std::make_pair(-inf, inf);
69 }
70
71 std::pair<double, double> getThroughputAtBounds() const override { return std::make_pair(1.0, 1.0); }
72
73 void sampleAt(lsst::geom::Point2D const&, ndarray::Array<double const, 1, 1> const& wavelengths,
74 ndarray::Array<double, 1, 1> const& out) const override {
75 out.deep() = 1.0;
76 }
77
78 bool isPersistable() const noexcept override { return true; }
79
80protected:
81 // transforming an IdentityTransmissionCurve is a no-op
84 return shared_from_this();
85 }
86
87 // multiplying an IdentityTransmissionCurve always yields the other operand
89 std::shared_ptr<TransmissionCurve const> other) const override {
90 return other;
91 }
92
93 std::string getPersistenceName() const override { return NAME; }
94
95 void write(OutputArchiveHandle& handle) const override { handle.saveEmpty(); }
96
97 class Factory : public table::io::PersistableFactory {
98 public:
99 std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
100 CatalogVector const& catalogs) const override {
101 LSST_ARCHIVE_ASSERT(catalogs.size() == 0u);
102 return get();
103 }
104
105 Factory() : table::io::PersistableFactory(NAME) {}
106 };
107
108 static Factory registration;
109
110private:
111 IdentityTransmissionCurve() = default;
112};
113
114IdentityTransmissionCurve::Factory IdentityTransmissionCurve::registration;
115
116/*
117 * InterpolatedTransmissionCurve: implements makeSpatiallyConstant and makeRadial.
118 *
119 * InterpolatedTransmissionCurve is templated on an implementation class
120 * that does most of the work, letting it handle the boilerplate that is
121 * common to both spatially-constant and radially-varying curves.
122 *
123 * We use two tricks to avoid repetition that bear some explanation:
124 *
125 * - Even though the two implementations have different state, we can use
126 * the same PersistenceHelper template for both by using a std::vector
127 * to hold the keys for the arrays. Impl2d has one extra array to save,
128 * so its vector is one element longer.
129 *
130 * - The throughput array held by Impl2d is conceptually a 2-d array, but we
131 * actually store a flattened view into it. That's okay because both GSL
132 * and the persistence layer only care about the pointer and the size. That
133 * requires going through some hoops in makeRadial, but from then on it's
134 * clean sailing.
135 */
136
137template <typename T>
138using GslPtr = std::unique_ptr<T, void (*)(T*)>;
139
140template <typename T>
141GslPtr<T> makeGslPtr(T* p, void (*free)(T*)) {
142 if (p == nullptr) {
143 throw std::bad_alloc();
144 }
145 return GslPtr<T>(p, free);
146}
147
148using ArrayKeyVector = std::vector<table::Key<table::Array<double>>>;
149
150/*
151 * Implementation object as a template parameter for the instantiation of
152 * InterpolatedTransmissionCurve used by makeSpatiallyConstant.
153 */
154class Impl1d {
155public:
156 static constexpr bool isSpatiallyConstant = true;
157 static constexpr char const* NAME = "SpatiallyConstantTransmissionCurve";
158
159 // Initialize the GSL interpolator and take ownership of the given arrays.
160 // Array size consistency checks are done by makeSpatiallyConstant,
161 // so we don't need to worry about them here.
162 Impl1d(ndarray::Array<double, 1, 1> const& throughput, ndarray::Array<double, 1, 1> const& wavelengths)
163 : _throughput(throughput),
164 _wavelengths(wavelengths),
165 _interp(makeGslPtr(::gsl_interp_alloc(::gsl_interp_linear, _wavelengths.size()),
166 &::gsl_interp_free)) {
167 int status = ::gsl_interp_init(_interp.get(), _wavelengths.getData(), _throughput.getData(),
168 _wavelengths.size());
169 LSST_CHECK_GSL(pex::exceptions::LogicError, status);
170 }
171
172 std::pair<double, double> getWavelengthBounds() const {
173 return std::make_pair(_wavelengths.front(), _wavelengths.back());
174 }
175
176 static void setupPersistence(table::Schema& schema, ArrayKeyVector& keys) {
177 keys.push_back(schema.addField<table::Array<double>>("throughput", "array of known throughput values",
178 "", 0));
179 keys.push_back(schema.addField<table::Array<double>>(
180 "wavelengths", "array of known wavelength values", "angstrom", 0));
181 }
182
183 void persist(table::BaseRecord& record, ArrayKeyVector const& keys) const {
184 record.set(keys[0], _throughput);
185 record.set(keys[1], _wavelengths);
186 }
187
188 static Impl1d unpersist(table::BaseRecord& record, ArrayKeyVector const& keys) {
189 return Impl1d(record[keys[0]], record[keys[1]]);
190 }
191
192 // A helper object constructed every time InterpolatedTransmissionCurve::sampleAt
193 // is called, and then invoked at every iteration of the loop therein.
194 struct Functor {
195 Functor(Impl1d const&, lsst::geom::Point2D const&)
196 : _accel(makeGslPtr(::gsl_interp_accel_alloc(), &gsl_interp_accel_free)) {}
197
198 double operator()(Impl1d const& parent, double wavelength) {
199 double result = 0.0;
200 int status = ::gsl_interp_eval_e(parent._interp.get(), parent._wavelengths.getData(),
201 parent._throughput.getData(), wavelength, _accel.get(), &result);
202 LSST_CHECK_GSL(pex::exceptions::RuntimeError, status);
203 return result;
204 }
205
206 private:
207 GslPtr<::gsl_interp_accel> _accel;
208 };
209
210private:
211 ndarray::Array<double, 1, 1> _throughput;
212 ndarray::Array<double, 1, 1> _wavelengths;
213 GslPtr<::gsl_interp> _interp;
214};
215
216/*
217 * Implementation object as a template parameter for the instantiation of
218 * InterpolatedTransmissionCurve used by makeRadial.
219 */
220class Impl2d {
221public:
222 static constexpr bool isSpatiallyConstant = false;
223 static constexpr char const* NAME = "RadialTransmissionCurve";
224
225 Impl2d(ndarray::Array<double, 1, 1> const& throughput, ndarray::Array<double, 1, 1> const& wavelengths,
226 ndarray::Array<double, 1, 1> const& radii)
227 : _throughput(throughput),
228 _wavelengths(wavelengths),
229 _radii(radii),
230 _interp(makeGslPtr(
231 ::gsl_interp2d_alloc(::gsl_interp2d_bilinear, _wavelengths.size(), _radii.size()),
232 &::gsl_interp2d_free)) {
233 int status = ::gsl_interp2d_init(_interp.get(), _wavelengths.getData(), _radii.getData(),
234 _throughput.getData(), _wavelengths.size(), _radii.size());
235 LSST_CHECK_GSL(pex::exceptions::LogicError, status);
236 }
237
238 std::pair<double, double> getWavelengthBounds() const {
239 return std::make_pair(_wavelengths.front(), _wavelengths.back());
240 }
241
242 static void setupPersistence(table::Schema& schema, ArrayKeyVector& keys) {
243 keys.push_back(schema.addField<table::Array<double>>(
244 "throughput",
245 "flattenned 2-d array of known throughput values, with radius dimension fastest", "", 0));
246 keys.push_back(schema.addField<table::Array<double>>(
247 "wavelengths", "array of known wavelength values", "angstrom", 0));
248 keys.push_back(schema.addField<table::Array<double>>("radii", "array of known radius values", "", 0));
249 }
250
251 void persist(table::BaseRecord& record, ArrayKeyVector const& keys) const {
252 record.set(keys[0], _throughput);
253 record.set(keys[1], _wavelengths);
254 record.set(keys[2], _radii);
255 }
256
257 static Impl2d unpersist(table::BaseRecord& record, ArrayKeyVector const& keys) {
258 return Impl2d(record[keys[0]], record[keys[1]], record[keys[2]]);
259 }
260
261 // A helper object constructed every time InterpolatedTransmissionCurve::sampleAt
262 // is called, and then invoked at every iteration of the loop therein.
263 struct Functor {
264 Functor(Impl2d const& parent, lsst::geom::Point2D const& point)
265 : _radius(point.asEigen().norm()),
266 _radiusAccel(makeGslPtr(::gsl_interp_accel_alloc(), &gsl_interp_accel_free)),
267 _wavelengthAccel(makeGslPtr(::gsl_interp_accel_alloc(), &gsl_interp_accel_free)) {
268 _radius = std::max(_radius, parent._radii.front());
269 _radius = std::min(_radius, parent._radii.back());
270 }
271
272 double operator()(Impl2d const& parent, double wavelength) {
273 double result = 0.0;
274 int status =
275 ::gsl_interp2d_eval_e(parent._interp.get(), parent._wavelengths.getData(),
276 parent._radii.getData(), parent._throughput.getData(), wavelength,
277 _radius, _radiusAccel.get(), _wavelengthAccel.get(), &result);
278 LSST_CHECK_GSL(pex::exceptions::RuntimeError, status);
279 return result;
280 }
281
282 private:
283 double _radius;
284 GslPtr<::gsl_interp_accel> _radiusAccel;
285 GslPtr<::gsl_interp_accel> _wavelengthAccel;
286 };
287
288private:
289 ndarray::Array<double, 1, 1> _throughput;
290 ndarray::Array<double, 1, 1> _wavelengths;
291 ndarray::Array<double, 1, 1> _radii;
292 GslPtr<::gsl_interp2d> _interp;
293};
294
295template <typename Impl>
296class InterpolatedTransmissionCurve : public TransmissionCurve {
297public:
298 InterpolatedTransmissionCurve(Impl impl, std::pair<double, double> throughputAtBounds)
299 : _atBounds(throughputAtBounds), _impl(std::move(impl)) {
300 if (!std::isfinite(_atBounds.first) || !std::isfinite(_atBounds.second)) {
301 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError,
302 "Throughput values at bounds must be finite");
303 }
304 }
305
306 std::pair<double, double> getWavelengthBounds() const override { return _impl.getWavelengthBounds(); }
307
308 std::pair<double, double> getThroughputAtBounds() const override { return _atBounds; }
309
310 void sampleAt(lsst::geom::Point2D const& point, ndarray::Array<double const, 1, 1> const& wavelengths,
311 ndarray::Array<double, 1, 1> const& out) const override {
312 LSST_THROW_IF_NE(wavelengths.getSize<0>(), out.getSize<0>(), pex::exceptions::LengthError,
313 "Length of wavelength array (%d) does not match size of output array (%d)");
314 typename Impl::Functor functor(_impl, point);
315 auto bounds = _impl.getWavelengthBounds();
316 auto wlIter = wavelengths.begin();
317 for (auto outIter = out.begin(); outIter != out.end(); ++outIter, ++wlIter) {
318 double& y = *outIter;
319 if (*wlIter < bounds.first) {
320 y = _atBounds.first;
321 } else if (*wlIter > bounds.second) {
322 y = _atBounds.second;
323 } else {
324 y = functor(_impl, *wlIter);
325 }
326 }
327 }
328
329 bool isPersistable() const noexcept override { return true; }
330
331protected:
334 if (_impl.isSpatiallyConstant) {
335 return shared_from_this();
336 } else {
338 }
339 }
340
341 std::string getPersistenceName() const override { return Impl::NAME; }
342
343 struct PersistenceHelper {
344 table::Schema schema;
345 table::Key<double> throughputAtMin;
346 table::Key<double> throughputAtMax;
347 ArrayKeyVector arrays;
348
349 static PersistenceHelper const& get() {
350 static PersistenceHelper const instance;
351 return instance;
352 }
353
354 private:
355 PersistenceHelper()
356 : schema(),
358 schema.addField<double>("throughputAtMin", "throughput below minimum wavelength")),
360 schema.addField<double>("throughputAtMax", "throughput above minimum wavelength")) {
361 Impl::setupPersistence(schema, arrays);
362 }
363 };
364
365 void write(OutputArchiveHandle& handle) const override {
366 auto const& keys = PersistenceHelper::get();
367 auto catalog = handle.makeCatalog(keys.schema);
368 auto record = catalog.addNew();
369 record->set(keys.throughputAtMin, _atBounds.first);
370 record->set(keys.throughputAtMax, _atBounds.second);
371 _impl.persist(*record, keys.arrays);
372 handle.saveCatalog(catalog);
373 }
374
375 class Factory : public table::io::PersistableFactory {
376 public:
377 std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
378 CatalogVector const& catalogs) const override {
379 auto const& keys = PersistenceHelper::get();
380 LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
381 LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
382 auto& record = catalogs.front().front();
383 return std::make_shared<InterpolatedTransmissionCurve>(
384 Impl::unpersist(record, keys.arrays),
385 std::make_pair(record.get(keys.throughputAtMin), record.get(keys.throughputAtMax)));
386 }
387
388 Factory() : table::io::PersistableFactory(Impl::NAME) {}
389 };
390
391 static Factory registration;
392
393private:
395 Impl _impl;
396};
397
398template <typename Impl>
399typename InterpolatedTransmissionCurve<Impl>::Factory InterpolatedTransmissionCurve<Impl>::registration;
400
401template class InterpolatedTransmissionCurve<Impl1d>;
402template class InterpolatedTransmissionCurve<Impl2d>;
403
404/*
405 * ProductTransmissionCurve: default for TransmissionCurve::multipliedBy().
406 *
407 * This is a straightforward lazy-evaluation object. Its only state is the
408 * two operands it delegates to.
409 */
410class ProductTransmissionCurve : public TransmissionCurve {
411public:
412 static constexpr char const* NAME = "ProductTransmissionCurve";
413
414 ProductTransmissionCurve(std::shared_ptr<TransmissionCurve const> a,
416 : _a(std::move(a)), _b(std::move(b)) {}
417
418 std::pair<double, double> getWavelengthBounds() const override {
419 auto aWavelengthBounds = _a->getWavelengthBounds();
420 auto bWavelengthBounds = _b->getWavelengthBounds();
421 auto aThroughputAtBounds = _a->getThroughputAtBounds();
422 auto bThroughputAtBounds = _b->getThroughputAtBounds();
423
424 auto determineWavelengthBound = [](double aWavelength, double bWavelength, double aThroughput,
425 double bThroughput, auto isFirstOuter) -> double {
426 // Use the outermost wavelength bound only if its throughput
427 // values are not being multiplied by zeros from the operand with
428 // the innermost wavelength bound.
429 if (isFirstOuter(aWavelength, bWavelength)) {
430 return (bThroughput == 0.0) ? bWavelength : aWavelength;
431 } else {
432 return (aThroughput == 0.0) ? aWavelength : bWavelength;
433 }
434 };
435
436 return std::make_pair(determineWavelengthBound(aWavelengthBounds.first, bWavelengthBounds.first,
437 aThroughputAtBounds.first, bThroughputAtBounds.first,
439 determineWavelengthBound(aWavelengthBounds.second, bWavelengthBounds.second,
440 aThroughputAtBounds.second, bThroughputAtBounds.second,
442 }
443
444 std::pair<double, double> getThroughputAtBounds() const override {
445 auto aAtBounds = _a->getThroughputAtBounds();
446 auto bAtBounds = _b->getThroughputAtBounds();
447 return std::make_pair(aAtBounds.first * bAtBounds.first, aAtBounds.second * bAtBounds.second);
448 }
449
450 void sampleAt(lsst::geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths,
451 ndarray::Array<double, 1, 1> const& out) const override {
452 _a->sampleAt(position, wavelengths, out);
453 ndarray::Array<double, 1, 1> tmp = ndarray::allocate(wavelengths.getSize<0>());
454 _b->sampleAt(position, wavelengths, tmp);
455 out.deep() *= tmp;
456 }
457
458 bool isPersistable() const noexcept override { return _a->isPersistable() && _b->isPersistable(); }
459
460protected:
461 std::string getPersistenceName() const override { return NAME; }
462
463 struct PersistenceHelper {
464 table::Schema schema;
465 table::Key<int> a;
466 table::Key<int> b;
467
468 static PersistenceHelper const& get() {
469 static PersistenceHelper const instance;
470 return instance;
471 }
472
473 private:
474 PersistenceHelper()
475 : schema(),
476 a(schema.addField<int>("a", "archive ID of first operand")),
477 b(schema.addField<int>("b", "archive ID of second operand")) {}
478 };
479
480 class Factory : public table::io::PersistableFactory {
481 public:
482 std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
483 CatalogVector const& catalogs) const override {
484 auto const& keys = PersistenceHelper::get();
485 LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
486 LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
487 auto const& record = catalogs.front().front();
488 return std::make_shared<ProductTransmissionCurve>(
489 archive.get<TransmissionCurve>(record.get(keys.a)),
490 archive.get<TransmissionCurve>(record.get(keys.b)));
491 }
492
493 Factory() : table::io::PersistableFactory(NAME) {}
494 };
495
496 void write(OutputArchiveHandle& handle) const override {
497 auto const& keys = PersistenceHelper::get();
498 auto catalog = handle.makeCatalog(keys.schema);
499 auto record = catalog.addNew();
500 record->set(keys.a, handle.put(_a));
501 record->set(keys.b, handle.put(_b));
502 handle.saveCatalog(catalog);
503 }
504
505 static Factory registration;
506
507private:
510};
511
512ProductTransmissionCurve::Factory ProductTransmissionCurve::registration;
513
514/*
515 * TransformedTransmissionCurve: default for TransmissionCurve::transform.
516 *
517 * This is a another straightforward lazy-evaluation object. Its only state
518 * is the two operands it delegates to.
519 */
520class TransformedTransmissionCurve : public TransmissionCurve {
521public:
522 static constexpr char const* NAME = "TransformedTransmissionCurve";
523
524 TransformedTransmissionCurve(std::shared_ptr<TransmissionCurve const> nested,
526 : _nested(std::move(nested)), _transform(std::move(transform)) {}
527
528 std::pair<double, double> getWavelengthBounds() const override { return _nested->getWavelengthBounds(); }
529
530 std::pair<double, double> getThroughputAtBounds() const override {
531 return _nested->getThroughputAtBounds();
532 }
533
534 void sampleAt(lsst::geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths,
535 ndarray::Array<double, 1, 1> const& out) const override {
536 return _nested->sampleAt(_transform->applyInverse(position), wavelengths, out);
537 }
538
539 bool isPersistable() const noexcept override {
540 return _nested->isPersistable() && _transform->isPersistable();
541 }
542
543protected:
544 // transforming a TransformedTransmissionCurve composes the transforms
547 return std::make_shared<TransformedTransmissionCurve>(_nested, transform->then(*_transform));
548 }
549
550 std::string getPersistenceName() const override { return NAME; }
551
552 struct PersistenceHelper {
553 table::Schema schema;
554 table::Key<int> nested;
555 table::Key<int> transform;
556
557 static PersistenceHelper const& get() {
558 static PersistenceHelper const instance;
559 return instance;
560 }
561
562 private:
563 PersistenceHelper()
564 : schema(),
565 nested(schema.addField<int>("nested", "archive ID of the nested TransmissionCurve")),
566 transform(schema.addField<int>("transform", "archive ID of the coordinate transform")) {}
567 };
568
569 void write(OutputArchiveHandle& handle) const override {
570 auto const& keys = PersistenceHelper::get();
571 auto catalog = handle.makeCatalog(keys.schema);
572 auto record = catalog.addNew();
573 record->set(keys.nested, handle.put(_nested));
574 record->set(keys.transform, handle.put(_transform));
575 handle.saveCatalog(catalog);
576 }
577
578 class Factory : public table::io::PersistableFactory {
579 public:
580 std::shared_ptr<table::io::Persistable> read(InputArchive const& archive,
581 CatalogVector const& catalogs) const override {
582 auto const& keys = PersistenceHelper::get();
583 LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
584 LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
585 auto const& record = catalogs.front().front();
586 return std::make_shared<TransformedTransmissionCurve>(
587 archive.get<TransmissionCurve>(record.get(keys.nested)),
588 archive.get<geom::TransformPoint2ToPoint2>(record.get(keys.transform)));
589 }
590
591 Factory() : table::io::PersistableFactory(NAME) {}
592 };
593
594 static Factory registration;
595
596private:
599};
600
601TransformedTransmissionCurve::Factory TransformedTransmissionCurve::registration;
602
603} // namespace
604
605/*
606 * TransmissionCurve itself
607 */
608
610 return IdentityTransmissionCurve::get();
611}
612
614 ndarray::Array<double const, 1> const& throughput, ndarray::Array<double const, 1> const& wavelengths,
615 double throughputAtMin, double throughputAtMax) {
616 ::gsl_set_error_handler_off();
617 LSST_THROW_IF_NE(wavelengths.getSize<0>(), throughput.getSize<0>(), pex::exceptions::LengthError,
618 "Length of wavelength array (%d) does not match size of throughput array (%d)");
619 return std::make_shared<InterpolatedTransmissionCurve<Impl1d>>(
620 Impl1d(ndarray::copy(throughput), ndarray::copy(wavelengths)),
622}
623
625 ndarray::Array<double const, 2> const& throughput, ndarray::Array<double const, 1> const& wavelengths,
626 ndarray::Array<double const, 1> const& radii, double throughputAtMin, double throughputAtMax) {
627 ::gsl_set_error_handler_off();
629 wavelengths.getSize<0>(), throughput.getSize<0>(), pex::exceptions::LengthError,
630 "Length of wavelength array (%d) does not match first dimension of of throughput array (%d)");
632 radii.getSize<0>(), throughput.getSize<1>(), pex::exceptions::LengthError,
633 "Length of radii array (%d) does not match second dimension of of throughput array (%d)");
634 // GSL wants a column major array (Array<T,2,-2>). But ndarray can only flatten row-major arrays
635 // (Array<T,2,2>). So we allocate a row-major array, assign the caller's throughput array to a
636 // transposed view of it, and then flatten the row-major array.
637 ndarray::Array<double, 2, 2> throughputTransposed = ndarray::allocate(throughput.getShape().reverse());
638 throughputTransposed.transpose() = throughput;
639 ndarray::Array<double, 1, 1> throughputFlat = ndarray::flatten<1>(throughputTransposed);
640 return std::make_shared<InterpolatedTransmissionCurve<Impl2d>>(
641 Impl2d(throughputFlat, ndarray::copy(wavelengths), ndarray::copy(radii)),
643}
644
646 TransmissionCurve const& other) const {
647 auto a = shared_from_this();
648 auto b = other.shared_from_this();
649 auto result = a->_multipliedByImpl(b);
650 if (result == nullptr) {
651 result = b->_multipliedByImpl(a);
652 }
653 if (result == nullptr) {
654 result = std::make_shared<ProductTransmissionCurve>(std::move(a), std::move(b));
655 }
656 return result;
657}
658
663
664ndarray::Array<double, 1, 1> TransmissionCurve::sampleAt(
665 lsst::geom::Point2D const& position, ndarray::Array<double const, 1, 1> const& wavelengths) const {
666 ndarray::Array<double, 1, 1> out = ndarray::allocate(wavelengths.getSize<0>());
667 sampleAt(position, wavelengths, out);
668 return out;
669}
670
675
680
681std::string TransmissionCurve::getPythonModule() const { return "lsst.afw.image"; }
682
683} // namespace image
684} // namespace afw
685} // namespace lsst
py::object result
Definition _schema.cc:429
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
int y
Definition SpanSet.cc:48
table::Key< int > transform
#define LSST_CHECK_GSL(type, status)
table::Key< int > nested
ArrayKeyVector arrays
table::Key< double > throughputAtMax
table::Key< double > throughputAtMin
table::Key< int > b
table::Key< int > a
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition Persistable.h:48
table::Schema schema
Definition python.h:134
#define LSST_THROW_IF_NE(N1, N2, EXC_CLASS, MSG)
Check whether the given values are equal, and throw an LSST Exception if they are not.
Definition asserts.h:38
A spatially-varying transmission curve as a function of wavelength.
static std::shared_ptr< TransmissionCurve const > makeRadial(ndarray::Array< double const, 2 > const &throughput, ndarray::Array< double const, 1 > const &wavelengths, ndarray::Array< double const, 1 > const &radii, double throughputAtMin=0.0, double throughputAtMax=0.0)
Create a new TransmissionCurve with throughput varying as function of radius.
virtual void sampleAt(lsst::geom::Point2D const &position, ndarray::Array< double const, 1, 1 > const &wavelengths, ndarray::Array< double, 1, 1 > const &out) const =0
Evaluate the throughput at a position into a provided output array.
std::shared_ptr< TransmissionCurve const > multipliedBy(TransmissionCurve const &other) const
Return a new TransmissionCurve that simply multiplies the values of two others.
virtual std::shared_ptr< TransmissionCurve const > _transformedByImpl(std::shared_ptr< geom::TransformPoint2ToPoint2 > transform) const
Polymorphic implementation for transformedBy().
std::string getPythonModule() const override
Return the fully-qualified Python module that should be imported to guarantee that its factory is reg...
std::shared_ptr< TransmissionCurve const > transformedBy(std::shared_ptr< geom::TransformPoint2ToPoint2 > transform) const
Return a view of a TransmissionCurve in a different coordinate system.
static std::shared_ptr< TransmissionCurve const > makeSpatiallyConstant(ndarray::Array< double const, 1 > const &throughput, ndarray::Array< double const, 1 > const &wavelengths, double throughputAtMin=0.0, double throughputAtMax=0.0)
Create a new TransmissionCurve with spatially-constant throughput.
virtual std::shared_ptr< TransmissionCurve const > _multipliedByImpl(std::shared_ptr< TransmissionCurve const > other) const
One-way polymorphic implementation for multipliedBy().
static std::shared_ptr< TransmissionCurve const > makeIdentity()
Create a new TranmissionCurve that has unit thoughput at all wavelengths everywhere.
EigenVector const & asEigen() const noexcept(IS_ELEMENT_NOTHROW_COPYABLE)
Return a fixed-size Eigen representation of the coordinate object.
Reports attempts to exceed implementation-defined length limits for some classes.
Definition Runtime.h:76
T infinity(T... args)
T isfinite(T... args)
T make_pair(T... args)
T max(T... args)
T min(T... args)
T move(T... args)
STL namespace.
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override