LSSTApplications  18.1.0
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"
38 #include "lsst/afw/table/io/Persistable.cc"
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  schema.getCitizen().markPersistent();
220  }
221 };
222 
223 template <typename FromEndpoint, typename ToEndpoint>
224 class TransformFactory : public table::io::PersistableFactory {
225 public:
226  explicit TransformFactory(std::string const &name) : table::io::PersistableFactory(name) {}
227 
228  std::shared_ptr<table::io::Persistable> read(InputArchive const &archive,
229  CatalogVector const &catalogs) const override {
230  auto const &keys = TransformPersistenceHelper::get();
231  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
232  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
233  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
234  auto const &record = catalogs.front().front();
235  std::string stringRep = formatters::bytesToString(record.get(keys.bytes));
236  return Transform<FromEndpoint, ToEndpoint>::readString(stringRep);
237  }
238 };
239 
240 } // namespace
241 
242 template <class FromEndpoint, class ToEndpoint>
244  auto const &keys = TransformPersistenceHelper::get();
245  table::BaseCatalog cat = handle.makeCatalog(keys.schema);
247  record->set(keys.bytes, formatters::stringToBytes(writeString()));
248  handle.saveCatalog(cat);
249 }
250 
251 #define INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, NextToEndpoint) \
252  template std::shared_ptr<Transform<FromEndpoint, NextToEndpoint>> \
253  Transform<FromEndpoint, ToEndpoint>::then<NextToEndpoint>( \
254  Transform<ToEndpoint, NextToEndpoint> const &next, bool) const;
255 
256 #define INSTANTIATE_TRANSFORM(FromEndpoint, ToEndpoint) \
257  } /* namespace geom */ \
258  template std::shared_ptr<geom::Transform<geom::FromEndpoint, geom::ToEndpoint>> \
259  table::io::PersistableFacade<geom::Transform<geom::FromEndpoint, geom::ToEndpoint>>::dynamicCast( \
260  std::shared_ptr<table::io::Persistable> const &); \
261  namespace geom { \
262  template class Transform<FromEndpoint, ToEndpoint>; \
263  template std::ostream &operator<<<FromEndpoint, ToEndpoint>( \
264  std::ostream &os, Transform<FromEndpoint, ToEndpoint> const &transform); \
265  namespace { \
266  TransformFactory<FromEndpoint, ToEndpoint> registration##FromEndpoint##ToEndpoint( \
267  Transform<FromEndpoint, ToEndpoint>::getShortClassName()); \
268  } /* namespace */ \
269  INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, GenericEndpoint) \
270  INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, Point2Endpoint) \
271  INSTANTIATE_OVERLOADS(FromEndpoint, ToEndpoint, SpherePointEndpoint)
272 
273 // explicit instantiations
274 INSTANTIATE_TRANSFORM(GenericEndpoint, GenericEndpoint);
275 INSTANTIATE_TRANSFORM(GenericEndpoint, Point2Endpoint);
276 INSTANTIATE_TRANSFORM(GenericEndpoint, SpherePointEndpoint);
277 INSTANTIATE_TRANSFORM(Point2Endpoint, GenericEndpoint);
278 INSTANTIATE_TRANSFORM(Point2Endpoint, Point2Endpoint);
279 INSTANTIATE_TRANSFORM(Point2Endpoint, SpherePointEndpoint);
280 INSTANTIATE_TRANSFORM(SpherePointEndpoint, GenericEndpoint);
281 INSTANTIATE_TRANSFORM(SpherePointEndpoint, Point2Endpoint);
282 INSTANTIATE_TRANSFORM(SpherePointEndpoint, SpherePointEndpoint);
283 
284 } // namespace geom
285 } // namespace afw
286 } // namespace lsst
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
An object passed to Persistable::write to allow it to persist itself.
table::Key< table::Array< std::uint8_t > > bytes
Definition: Transform.cc:199
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
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
table::Schema schema
Definition: Transform.cc:198
#define INSTANTIATE_TRANSFORM(FromEndpoint, ToEndpoint)
Definition: Transform.cc:256
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
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
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
std::ostream * os
Definition: Schema.cc:746
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