LSST Applications 24.1.6,g063fba187b+56b85ce14a,g0f08755f38+df8a265115,g12f32b3c4e+891a09f10d,g1524ad2192+7a5d7b3fbd,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g28da252d5a+07cb1400be,g2bbee38e9b+ae03bbfc84,g2bc492864f+ae03bbfc84,g3156d2b45e+6e55a43351,g347aa1857d+ae03bbfc84,g35bb328faa+a8ce1bb630,g3a166c0a6a+ae03bbfc84,g3e281a1b8c+c5dd892a6c,g414038480c+6b9177ef31,g41af890bb2+8f257c4c0b,g781aacb6e4+a8ce1bb630,g7af13505b9+7137b3b17d,g80478fca09+6df6903293,g82479be7b0+091ce1d07f,g858d7b2824+df8a265115,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,g9726552aa6+414189b318,ga5288a1d22+4a2bca08d7,gacef1a1666+c9a8ff65f4,gb58c049af0+d64f4d3760,gbcfae0f0a0+de1d42d831,gc28159a63d+ae03bbfc84,gcf0d15dbbd+72117bf34e,gda6a2b7d83+72117bf34e,gdaeeff99f8+1711a396fd,ge500cccec5+c8c9c9af63,ge79ae78c31+ae03bbfc84,gf0baf85859+c1f95f4921,gfa517265be+df8a265115,gfa999e8aa5+17cd334064,gfb92a5be7c+df8a265115
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
Definition __init__.py:1
T size(T... args)