LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
Mapping.cc
Go to the documentation of this file.
1/*
2 * LSST Data Management System
3 * Copyright 2017 AURA/LSST.
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 <https://www.lsstcorp.org/LegalNotices/>.
21 */
22#include <cmath>
23#include <limits>
24#include <memory>
25#include <sstream>
26#include <stdexcept>
27
28#include "astshim/base.h"
30#include "astshim/Frame.h"
31#include "astshim/Mapping.h"
32#include "astshim/ParallelMap.h"
33#include "astshim/SeriesMap.h"
34
35namespace ast {
36
37SeriesMap Mapping::then(Mapping const &next) const { return SeriesMap(*this, next); }
38
39ParallelMap Mapping::under(Mapping const &next) const { return ParallelMap(*this, next); }
40
42 auto rawCopy = reinterpret_cast<AstMapping *>(astCopy(getRawPtr()));
43 astInvert(rawCopy);
44 assertOK(reinterpret_cast<AstObject *>(rawCopy));
45 // use false because the pointer has already been copied
46 return Object::fromAstObject<Mapping>(reinterpret_cast<AstObject *>(rawCopy), false);
47}
48
49Array2D Mapping::linearApprox(PointD const &lbnd, PointD const &ubnd, double tol) const {
50 int const nIn = getNIn();
51 int const nOut = getNOut();
52 detail::assertEqual(lbnd.size(), "lbnd.size", static_cast<std::size_t>(nIn), "nIn");
53 detail::assertEqual(ubnd.size(), "ubnd.size", static_cast<std::size_t>(nIn), "nIn");
54 Array2D fit = ndarray::allocate(ndarray::makeVector(1 + nIn, nOut));
55 int isOK = astLinearApprox(getRawPtr(), lbnd.data(), ubnd.data(), tol, fit.getData());
56 assertOK();
57 if (!isOK) {
58 throw std::runtime_error("Mapping not sufficiently linear");
59 }
60 return fit;
61}
62
63template <typename Class>
65 if ((i < 0) || (i > 1)) {
67 os << "i =" << i << "; must be 0 or 1";
68 throw std::invalid_argument(os.str());
69 }
70 // Report pre-existing problems now so our later test for "not a compound object" is accurate
71 assertOK();
72
73 AstMapping *rawMap1;
74 AstMapping *rawMap2;
75 int series, invert1, invert2;
76 astDecompose(getRawPtr(), &rawMap1, &rawMap2, &series, &invert1, &invert2);
77 assertOK();
78
79 if (!rawMap2) {
80 // Not a compound object; free rawMap1 (rawMap2 is null, so no need to free it) and throw an exception
81 astAnnul(reinterpret_cast<AstObject *>(rawMap1));
83 os << "This " << getClassName() << " is not a compound object";
84 throw std::runtime_error(os.str());
85 }
86
87 // Make a deep copy of the returned object and free the shallow copies
88 AstMapping *retRawMap;
89 int invert;
90 if (i == 0) {
91 retRawMap = reinterpret_cast<AstMapping *>(astCopy(reinterpret_cast<AstObject *>(rawMap1)));
92 invert = invert1;
93 } else {
94 retRawMap = reinterpret_cast<AstMapping *>(astCopy(reinterpret_cast<AstObject *>(rawMap2)));
95 invert = invert2;
96 }
97 astAnnul(reinterpret_cast<AstObject *>(rawMap1));
98 astAnnul(reinterpret_cast<AstObject *>(rawMap2));
99 assertOK();
100
101 // If the mapping's internal invert flag does not match the value used when the CmpMap was made
102 // then invert the mapping. Note that it is not possible to create such objects in astshim
103 // but it is possible to read in objects created by other software.
104 if (invert != astGetI(retRawMap, "Invert")) {
105 astInvert(retRawMap);
106 assertOK();
107 }
108
109 return Object::fromAstObject<Class>(reinterpret_cast<AstObject *>(retRawMap), copy);
110}
111
112void Mapping::_tran(ConstArray2D const &from, bool doForward, Array2D const &to) const {
113 int const nFromAxes = doForward ? getNIn() : getNOut();
114 int const nToAxes = doForward ? getNOut() : getNIn();
115 detail::assertEqual(from.getSize<0>(), "from.size[0]", static_cast<std::size_t>(nFromAxes),
116 "from coords");
117 detail::assertEqual(to.getSize<0>(), "to.size[0]", static_cast<std::size_t>(nToAxes), "to coords");
118 detail::assertEqual(from.getSize<1>(), "from.size[1]", to.getSize<1>(), "to.size[1]");
119 int const nPts = from.getSize<1>();
120 // astTranN treats 0 points as an error and the call isn't needed anyway
121 if (nPts > 0) {
122 astTranN(getRawPtr(), nPts, nFromAxes, nPts, from.getData(), static_cast<int>(doForward), nToAxes, nPts,
123 to.getData());
124 }
125 assertOK();
127}
128
129void Mapping::_tranGrid(PointI const &lbnd, PointI const &ubnd, double tol, int maxpix, bool doForward,
130 Array2D const &to) const {
131 int const nFromAxes = doForward ? getNIn() : getNOut();
132 int const nToAxes = doForward ? getNOut() : getNIn();
133 detail::assertEqual(lbnd.size(), "lbnd.size", static_cast<std::size_t>(nFromAxes), "from coords");
134 detail::assertEqual(ubnd.size(), "ubnd.size", static_cast<std::size_t>(nFromAxes), "from coords");
135 detail::assertEqual(to.getSize<1>(), "to.size[1]", static_cast<std::size_t>(nToAxes), "to coords");
136 int const nPts = to.getSize<0>();
137 astTranGrid(getRawPtr(), nFromAxes, lbnd.data(), ubnd.data(), tol, maxpix, static_cast<int>(doForward),
138 nToAxes, nPts, to.getData());
139 assertOK();
141}
142
143// Explicit instantiations
144template std::shared_ptr<Frame> Mapping::decompose(int i, bool) const;
145template std::shared_ptr<Mapping> Mapping::decompose(int i, bool) const;
146
147} // namespace ast
std::ostream * os
Definition Schema.cc:557
An abstract base class for objects which transform one set of coordinates to another.
Definition Mapping.h:59
SeriesMap then(Mapping const &next) const
Return a series compound mapping this(first(input)) containing shallow copies of the original.
Definition Mapping.cc:37
std::shared_ptr< Class > decompose(int i, bool copy) const
Return a deep copy of one of the two component mappings.
Definition Mapping.cc:64
ParallelMap under(Mapping const &next) const
Return a parallel compound mapping containing shallow copies of the original.
Definition Mapping.cc:39
Array2D linearApprox(PointD const &lbnd, PointD const &ubnd, double tol) const
Compute a linear approximation to the forward transformation.
Definition Mapping.cc:49
int getNOut() const
Get NOut: the number of output axes.
Definition Mapping.h:82
int getNIn() const
Get NIn: the number of input axes.
Definition Mapping.h:77
std::shared_ptr< Mapping > inverted() const
Get an inverse mapping.
Definition Mapping.cc:41
std::string getClassName() const
Get Class: the name of the class (e.g.
Definition Object.h:139
AstObject const * getRawPtr() const
Get the raw AST pointer.
Definition Object.h:292
A parallel compound mapping where the first Mapping is used to transform the lower numbered coordinat...
Definition ParallelMap.h:53
A series compound mapping where the first Mapping is used to transform the coordinates of each point ...
Definition SeriesMap.h:53
T data(T... args)
void assertEqual(T1 val1, std::string const &descr1, T2 val2, std::string const &descr2)
Definition utils.h:48
void astBadToNan(std::vector< double > &p)
Replace AST__BAD with a quiet NaN in a vector.
Definition utils.h:59
AST wrapper classes and functions.
ndarray::Array< const double, 2, 2 > ConstArray2D
2D array of const double; typically used for lists of const points
Definition base.h:46
std::vector< int > PointI
Vector of ints; typically used for the bounds of Mapping.tranGridForward and inverse.
Definition base.h:50
ndarray::Array< double, 2, 2 > Array2D
2D array of const double; typically used for lists of const points
Definition base.h:42
void assertOK(AstObject *rawPtr1=nullptr, AstObject *rawPtr2=nullptr)
Throw std::runtime_error if AST's state is bad.
Definition base.cc:49
T size(T... args)