LSSTApplications  17.0+118,17.0+12,17.0+68,18.0.0+32,18.0.0+6,18.0.0+73,18.0.0-4-g68ffd23+2,18.1.0-1-g0001055+10,18.1.0-1-g03d53ef+3,18.1.0-1-g1349e88+48,18.1.0-1-g2505f39+38,18.1.0-1-g5315e5e+2,18.1.0-1-g5e4b7ea+12,18.1.0-1-g7e8fceb+2,18.1.0-1-g85f8cd4+41,18.1.0-1-g8ff0b9f+1,18.1.0-1-gd55f500+29,18.1.0-13-gfe4edf0b+4,18.1.0-14-g259bd21+13,18.1.0-16-gf53fc9f+2,18.1.0-2-g5f9922c+17,18.1.0-2-gd3b74e5+6,18.1.0-2-gfbf3545+25,18.1.0-2-gfefb8b5+37,18.1.0-24-gd64d31bc+2,18.1.0-24-ged780bc+2,18.1.0-27-g8f4a40ebd,18.1.0-3-g52aa583+20,18.1.0-3-g8ea57af+2,18.1.0-3-g8f4a2b1+35,18.1.0-3-gb69f684+34,18.1.0-5-g1dd662b,18.1.0-5-g6dbcb01+34,18.1.0-6-gae77429+1,18.1.0-7-g9d75d83+2,18.1.0-7-gae09a6d+22,18.1.0-8-gc69d46e+20,18.1.0-9-g1af92ce+2,18.1.0-9-gee19f03+6,w.2019.44
LSSTDataManagementBasePackage
Transform.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2008-2017 LSST Corporation.
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 <exception>
24 #include <memory>
25 #include <ostream>
26 #include <sstream>
27 #include <vector>
28 
29 #include "astshim.h"
32 #include "lsst/afw/geom/Endpoint.h"
34 #include "lsst/afw/geom/SkyWcs.h"
39 
40 namespace lsst {
41 namespace afw {
42 namespace geom {
43 
44 template <class FromEndpoint, class ToEndpoint>
46  : _fromEndpoint(mapping.getNIn()),
47  _mapping(simplify ? mapping.simplified() : mapping.copy()),
48  _toEndpoint(mapping.getNOut()) {}
49 
50 template <typename FromEndpoint, typename ToEndpoint>
52  : _fromEndpoint(frameSet.getNIn()), _mapping(), _toEndpoint(frameSet.getNOut()) {
53  auto frameSetCopy = frameSet.copy();
54  // Normalize the base and current frame in a way that affects its behavior as a mapping.
55  // To do this one must set the current frame to the frame to be normalized
56  // and normalize the frame set as a frame (i.e. normalize the frame "in situ").
57  // The obvious alternative of normalizing a shallow copy of the frame does not work;
58  // the frame is altered but not the associated mapping!
59 
60  // Normalize the current frame by normalizing the frameset as a frame
61  _toEndpoint.normalizeFrame(frameSetCopy);
62 
63  // Normalize the base frame by temporarily making it the current frame,
64  // normalizing the frameset as a frame, then making it the base frame again
65  const int baseIndex = frameSetCopy->getBase();
66  const int currentIndex = frameSetCopy->getCurrent();
67  frameSetCopy->setCurrent(baseIndex);
68  _fromEndpoint.normalizeFrame(frameSetCopy);
69  frameSetCopy->setBase(baseIndex);
70  frameSetCopy->setCurrent(currentIndex);
71  _mapping = simplify ? frameSetCopy->getMapping()->simplified() : frameSetCopy->getMapping();
72 }
73 
74 template <typename FromEndpoint, typename ToEndpoint>
76  : _fromEndpoint(mapping->getNIn()), _mapping(mapping), _toEndpoint(mapping->getNOut()) {}
77 
78 template <class FromEndpoint, class ToEndpoint>
80  typename FromEndpoint::Point const &point) const {
81  auto const rawFromData = _fromEndpoint.dataFromPoint(point);
82  auto rawToData = _mapping->applyForward(rawFromData);
83  return _toEndpoint.pointFromData(rawToData);
84 }
85 
86 template <class FromEndpoint, class ToEndpoint>
88  typename FromEndpoint::Array const &array) const {
89  auto const rawFromData = _fromEndpoint.dataFromArray(array);
90  auto rawToData = _mapping->applyForward(rawFromData);
91  return _toEndpoint.arrayFromData(rawToData);
92 }
93 
94 template <class FromEndpoint, class ToEndpoint>
96  typename ToEndpoint::Point const &point) const {
97  auto const rawFromData = _toEndpoint.dataFromPoint(point);
98  auto rawToData = _mapping->applyInverse(rawFromData);
99  return _fromEndpoint.pointFromData(rawToData);
100 }
101 
102 template <class FromEndpoint, class ToEndpoint>
104  typename ToEndpoint::Array const &array) const {
105  auto const rawFromData = _toEndpoint.dataFromArray(array);
106  auto rawToData = _mapping->applyInverse(rawFromData);
107  return _fromEndpoint.arrayFromData(rawToData);
108 }
109 
110 template <class FromEndpoint, class ToEndpoint>
112  auto inverse = std::dynamic_pointer_cast<ast::Mapping>(_mapping->inverted());
113  if (!inverse) {
114  // don't throw std::bad_cast because it doesn't let you provide debugging info
115  std::ostringstream buffer;
116  buffer << "Mapping.inverted() does not return a Mapping. Called from: " << _mapping;
118  }
119  return std::make_shared<Transform<ToEndpoint, FromEndpoint>>(*inverse);
120 }
121 
122 template <class FromEndpoint, class ToEndpoint>
124  int const nIn = _fromEndpoint.getNAxes();
125  int const nOut = _toEndpoint.getNAxes();
126  std::vector<double> const point = _fromEndpoint.dataFromPoint(x);
127 
128  Eigen::MatrixXd jacobian(nOut, nIn);
129  for (int i = 0; i < nOut; ++i) {
130  for (int j = 0; j < nIn; ++j) {
131  jacobian(i, j) = _mapping->rate(point, i + 1, j + 1);
132  }
133  }
134  return jacobian;
135 }
136 
137 template <class FromEndpoint, class ToEndpoint>
140  os << "Transform" << FromEndpoint::getClassPrefix() << "To" << ToEndpoint::getClassPrefix();
141  return os.str();
142 }
143 
144 template <class FromEndpoint, class ToEndpoint>
146  std::istream &is) {
147  return detail::readStream<Transform<FromEndpoint, ToEndpoint>>(is);
148 }
149 
150 template <class FromEndpoint, class ToEndpoint>
152  std::string &str) {
153  std::istringstream is(str);
155 }
156 
157 template <class FromEndpoint, class ToEndpoint>
159  detail::writeStream<Transform<FromEndpoint, ToEndpoint>>(*this, os);
160 }
161 
162 template <class FromEndpoint, class ToEndpoint>
165  writeStream(os);
166  return os.str();
167 }
168 
169 template <class FromEndpoint, class ToEndpoint>
170 template <class NextToEndpoint>
172  Transform<ToEndpoint, NextToEndpoint> const &next, bool simplify) const {
173  if (_toEndpoint.getNAxes() == next.getFromEndpoint().getNAxes()) {
174  auto nextMapping = next.getMapping();
175  auto combinedMapping = getMapping()->then(*next.getMapping());
176  if (simplify) {
177  return std::make_shared<Transform<FromEndpoint, NextToEndpoint>>(*combinedMapping.simplified());
178  } else {
179  return std::make_shared<Transform<FromEndpoint, NextToEndpoint>>(combinedMapping);
180  }
181  } else {
182  auto message = "Cannot match " + std::to_string(_toEndpoint.getNAxes()) + "-D to-endpoint to " +
183  std::to_string(next.getFromEndpoint().getNAxes()) + "-D from-endpoint.";
185  }
186 }
187 
188 template <class FromEndpoint, class ToEndpoint>
189 std::ostream &operator<<(std::ostream &os, Transform<FromEndpoint, ToEndpoint> const &transform) {
190  os << "Transform<" << transform.getFromEndpoint() << ", " << transform.getToEndpoint() << ">";
191  return os;
192 };
193 
194 namespace {
195 
196 class TransformPersistenceHelper {
197 public:
198  table::Schema schema;
199  table::Key<table::Array<std::uint8_t>> bytes;
200 
201  static TransformPersistenceHelper const &get() {
202  static TransformPersistenceHelper instance;
203  return instance;
204  }
205 
206  // No copying
207  TransformPersistenceHelper(TransformPersistenceHelper const &) = delete;
208  TransformPersistenceHelper &operator=(TransformPersistenceHelper const &) = delete;
209 
210  // No moving
211  TransformPersistenceHelper(TransformPersistenceHelper &&) = delete;
212  TransformPersistenceHelper &operator=(TransformPersistenceHelper &&) = delete;
213 
214 private:
215  TransformPersistenceHelper()
216  : schema(),
217  bytes(schema.addField<table::Array<std::uint8_t>>(
218  "bytes", "a bytestring containing the output of Transform.writeString", "")) {}
219 };
220 
221 template <typename FromEndpoint, typename ToEndpoint>
222 class TransformFactory : public table::io::PersistableFactory {
223 public:
224  explicit TransformFactory(std::string const &name) : table::io::PersistableFactory(name) {}
225 
226  std::shared_ptr<table::io::Persistable> read(InputArchive const &archive,
227  CatalogVector const &catalogs) const override {
228  auto const &keys = TransformPersistenceHelper::get();
229  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
230  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
231  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
232  auto const &record = catalogs.front().front();
233  std::string stringRep = formatters::bytesToString(record.get(keys.bytes));
234  return Transform<FromEndpoint, ToEndpoint>::readString(stringRep);
235  }
236 };
237 
238 } // namespace
239 
240 template <class FromEndpoint, class ToEndpoint>
242  auto const &keys = TransformPersistenceHelper::get();
243  table::BaseCatalog cat = handle.makeCatalog(keys.schema);
245  record->set(keys.bytes, formatters::stringToBytes(writeString()));
246  handle.saveCatalog(cat);
247 }
248 
249 #define INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, NextToEndpoint) \
250  template std::shared_ptr<Transform<FromEndpoint, NextToEndpoint>> \
251  Transform<FromEndpoint, ToEndpoint>::then<NextToEndpoint>( \
252  Transform<ToEndpoint, NextToEndpoint> const &next, bool) const;
253 
254 #define INSTANTIATE_TRANSFORM(FromEndpoint, ToEndpoint) \
255  } /* namespace geom */ \
256  template std::shared_ptr<geom::Transform<geom::FromEndpoint, geom::ToEndpoint>> \
257  table::io::PersistableFacade<geom::Transform<geom::FromEndpoint, geom::ToEndpoint>>::dynamicCast( \
258  std::shared_ptr<table::io::Persistable> const &); \
259  namespace geom { \
260  template class Transform<FromEndpoint, ToEndpoint>; \
261  template std::ostream &operator<<<FromEndpoint, ToEndpoint>( \
262  std::ostream &os, Transform<FromEndpoint, ToEndpoint> const &transform); \
263  namespace { \
264  TransformFactory<FromEndpoint, ToEndpoint> registration##FromEndpoint##ToEndpoint( \
265  Transform<FromEndpoint, ToEndpoint>::getShortClassName()); \
266  } /* namespace */ \
267  INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, GenericEndpoint) \
268  INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, Point2Endpoint) \
269  INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, SpherePointEndpoint)
270 
271 // explicit instantiations
272 INSTANTIATE_TRANSFORM(GenericEndpoint, GenericEndpoint);
273 INSTANTIATE_TRANSFORM(GenericEndpoint, Point2Endpoint);
274 INSTANTIATE_TRANSFORM(GenericEndpoint, SpherePointEndpoint);
275 INSTANTIATE_TRANSFORM(Point2Endpoint, GenericEndpoint);
276 INSTANTIATE_TRANSFORM(Point2Endpoint, Point2Endpoint);
277 INSTANTIATE_TRANSFORM(Point2Endpoint, SpherePointEndpoint);
278 INSTANTIATE_TRANSFORM(SpherePointEndpoint, GenericEndpoint);
279 INSTANTIATE_TRANSFORM(SpherePointEndpoint, Point2Endpoint);
280 INSTANTIATE_TRANSFORM(SpherePointEndpoint, SpherePointEndpoint);
281 
282 } // namespace geom
283 } // namespace afw
284 } // namespace lsst
table::Key< table::Array< std::uint8_t > > bytes
Definition: Transform.cc:199
def write(self, patchRef, catalog)
Write the output.
FromEndpoint getFromEndpoint() const
Get the "from" endpoint.
Definition: Transform.h:130
ToPoint applyForward(FromPoint const &point) const
Transform one point in the forward direction ("from" to "to")
Definition: Transform.cc:79
#define INSTANTIATE_TRANSFORM(FromEndpoint, ToEndpoint)
Definition: Transform.cc:254
An object passed to Persistable::write to allow it to persist itself.
std::shared_ptr< const ast::Mapping > getMapping() const
Get the contained mapping.
Definition: Transform.h:135
std::shared_ptr< Transform > readStream(std::istream &is)
Deserialize a Transform from an input stream.
T to_string(T... args)
Transform(Transform const &)=default
void writeStream(Transform const &transform, std::ostream &os)
Serialize a Transform to an output stream.
An abstract base class for objects which transform one set of coordinates to another.
Definition: Mapping.h:59
std::string bytesToString(ndarray::Array< std::uint8_t const, 1, 1 > const &bytes)
Decode a std::string from a vector of uint8 returned by stringToBytes.
Definition: Utils.cc:173
STL class.
STL class.
PolynomialFunction1d simplified(ScaledPolynomialFunction1d const &f)
Calculate the standard polynomial function that is equivalent to a scaled standard polynomial functio...
std::shared_ptr< FrameSet > copy() const
Return a deep copy of this object.
Definition: FrameSet.h:147
A base class for image defects.
FromPoint applyInverse(ToPoint const &point) const
Transform one point in the inverse direction ("to" to "from")
Definition: Transform.cc:95
T str(T... args)
T dynamic_pointer_cast(T... args)
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
double x
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
table::Schema schema
Definition: Transform.cc:198
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
Reports invalid arguments.
Definition: Runtime.h:66
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
table::Key< int > transform
ndarray::Array< std::uint8_t, 1, 1 > stringToBytes(std::string const &str)
Encode a std::string as a vector of uint8.
Definition: Utils.cc:162
std::ostream * os
Definition: Schema.cc:746
STL class.
A FrameSet consists of a set of one or more Frames (which describe coordinate systems), connected together by Mappings (which describe how the coordinate systems are inter-related).
Definition: FrameSet.h:99
Transform LSST spatial data, such as lsst::geom::Point2D and lsst::geom::SpherePoint, using an AST mapping.
Definition: Transform.h:67
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:485