LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
TanWcsFormatter.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 
38 #ifndef __GNUC__
39 # define __attribute__(x) /*NOTHING*/
40 #endif
41 static char const* SVNid __attribute__((unused)) = "$Id$";
42 
43 //#include "boost/serialization/shared_ptr.hpp"
44 #include <boost/archive/binary_iarchive.hpp>
45 #include <boost/archive/binary_oarchive.hpp>
46 #include <boost/archive/text_iarchive.hpp>
47 #include <boost/archive/text_oarchive.hpp>
48 
49 #include "wcslib/wcs.h"
50 
51 #include "lsst/daf/base.h"
52 #include "lsst/daf/persistence.h"
54 #include "lsst/pex/exceptions.h"
55 #include "lsst/pex/logging/Trace.h"
59 #include "lsst/afw/image/TanWcs.h"
60 
61 #define EXEC_TRACE 20
62 static void execTrace(std::string s, int level = EXEC_TRACE) {
63  lsst::pex::logging::Trace("afw.TanWcsFormatter", level, s);
64 }
65 
66 
67 namespace afwForm = lsst::afw::formatters;
68 namespace afwImg = lsst::afw::image;
69 namespace dafBase = lsst::daf::base;
70 namespace dafPersist = lsst::daf::persistence;
71 namespace pexPolicy = lsst::pex::policy;
72 namespace pexExcept = lsst::pex::exceptions;
73 
74 
76  "TanWcs", typeid(afwImg::TanWcs), createInstance);
77 
80  dafPersist::Formatter(typeid(this)) {
81 }
82 
84 }
85 
87  dafBase::Persistable const* persistable,
90  execTrace("TamWcsFormatter write start");
91  afwImg::TanWcs const* ip = dynamic_cast<afwImg::TanWcs const*>(persistable);
92  if (ip == 0) {
93  throw LSST_EXCEPT(pexExcept::RuntimeError, "Persisting non-TanWcs");
94  }
95  if (typeid(*storage) == typeid(dafPersist::BoostStorage)) {
96  execTrace("TanWcsFormatter write BoostStorage");
97  dafPersist::BoostStorage* boost = dynamic_cast<dafPersist::BoostStorage*>(storage.get());
98  boost->getOArchive() & *ip;
99  execTrace("TanWcsFormatter write end");
100  return;
101  }
102  throw LSST_EXCEPT(pexExcept::RuntimeError, "Unrecognized Storage for TanWcs");
103 }
104 
106  dafPersist::Storage::Ptr storage,
107  dafBase::PropertySet::Ptr additionalData) {
108  execTrace("TanWcsFormatter read start");
109  if (typeid(*storage) == typeid(dafPersist::BoostStorage)) {
110  afwImg::TanWcs* ip = new afwImg::TanWcs;
111  execTrace("TanWcsFormatter read BoostStorage");
112  dafPersist::BoostStorage* boost = dynamic_cast<dafPersist::BoostStorage*>(storage.get());
113  boost->getIArchive() & *ip;
114  execTrace("TanWcsFormatter read end");
115  return ip;
116  }
117  else if (typeid(*storage) == typeid(dafPersist::FitsStorage)) {
118  execTrace("TanWcsFormatter read FitsStorage");
119  dafPersist::FitsStorage* fits = dynamic_cast<dafPersist::FitsStorage*>(storage.get());
120  int hdu = additionalData->get<int>("hdu", 0);
122  afwImg::readMetadata(fits->getPath(), hdu);
123  afwImg::TanWcs* ip = new afwImg::TanWcs(md);
124  execTrace("TanWcsFormatter read end");
125  return ip;
126  }
127  throw LSST_EXCEPT(pexExcept::RuntimeError, "Unrecognized Storage for TanWcs");
128 }
129 
134  throw LSST_EXCEPT(pexExcept::RuntimeError, "Unexpected call to update for TanWcs");
135 }
136 
137 
139 template <class Archive>
140 void serializeEigenArray(Archive& ar, Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& m) {
141  int rows = m.rows();
142  int cols = m.cols();
143  ar & rows & cols;
144  if (Archive::is_loading::value) {
145  m = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>(rows, cols);
146  }
147  for (int j = 0; j < m.cols(); ++j) {
148  for (int i = 0; i < m.rows(); ++i) {
149  ar & m(i,j);
150  }
151  }
152 }
153 
154 
155 static void encodeSipHeader(lsst::daf::base::PropertySet::Ptr wcsProps,
156  std::string const& which,
157  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> const& m) {
158  int order = m.rows();
159  if (m.cols() != order) {
160  throw LSST_EXCEPT(pexExcept::DomainError,
161  "sip" + which + " matrix is not square");
162  }
163  if (order > 0) {
164  order -= 1; // match SIP convention
165  wcsProps->add(which + "_ORDER", static_cast<int>(order));
166  for (int i = 0; i <= order; ++i) {
167  for (int j = 0; j <= order; ++j) {
168  double val = m(i, j);
169  if (val != 0.0) {
170  wcsProps->add((boost::format("%1%_%2%_%3%")
171  % which % i % j).str(), val);
172  }
173  }
174  }
175  }
176 }
177 
180  // Only generates properties for the first wcsInfo.
182 
183  if (wcs._wcsInfo == NULL) { // nothing to add
184  return wcsProps;
185  }
186 
187  wcsProps->add("NAXIS", wcs._wcsInfo[0].naxis, "number of data axes");
188  wcsProps->add("EQUINOX", wcs._wcsInfo[0].equinox, "Equinox of coordinates");
189  wcsProps->add("RADESYS", std::string(wcs._wcsInfo[0].radesys), "Coordinate system for equinox");
190  wcsProps->add("CRPIX1", wcs._wcsInfo[0].crpix[0], "WCS Coordinate reference pixel");
191  wcsProps->add("CRPIX2", wcs._wcsInfo[0].crpix[1], "WCS Coordinate reference pixel");
192  wcsProps->add("CD1_1", wcs._wcsInfo[0].cd[0], "WCS Coordinate scale matrix");
193  wcsProps->add("CD1_2", wcs._wcsInfo[0].cd[1], "WCS Coordinate scale matrix");
194  wcsProps->add("CD2_1", wcs._wcsInfo[0].cd[2], "WCS Coordinate scale matrix");
195  wcsProps->add("CD2_2", wcs._wcsInfo[0].cd[3], "WCS Coordinate scale matrix");
196  wcsProps->add("CRVAL1", wcs._wcsInfo[0].crval[0], "WCS Ref value (RA in decimal degrees)");
197  wcsProps->add("CRVAL2", wcs._wcsInfo[0].crval[1], "WCS Ref value (DEC in decimal degrees)");
198  wcsProps->add("CUNIT1", std::string(wcs._wcsInfo[0].cunit[0]));
199  wcsProps->add("CUNIT2", std::string(wcs._wcsInfo[0].cunit[1]));
200 
201  //Hack. Because wcslib4.3 gets confused when it's passed RA---TAN-SIP,
202  //we set the value of ctypes to just RA---TAN, regardless of whether
203  //the SIP types are present. But when we persist to a file, we need to
204  //check whether the SIP polynomials were actually there and correct
205  //ctypes if necessary. Bad things will happen if someone tries to
206  //use a system other than RA---TAN and DEC--TAN
207  std::string ctype1(wcs._wcsInfo[0].ctype[0]);
208  std::string ctype2(wcs._wcsInfo[0].ctype[1]);
209 
210  if (wcs._hasDistortion) {
211  if (ctype1.rfind("-SIP") == std::string::npos) {
212  ctype1 += "-SIP";
213  }
214  if (ctype2.rfind("-SIP") == std::string::npos) {
215  ctype2 += "-SIP";
216  }
217  encodeSipHeader(wcsProps, "A", wcs._sipA);
218  encodeSipHeader(wcsProps, "B", wcs._sipB);
219  encodeSipHeader(wcsProps, "AP", wcs._sipAp);
220  encodeSipHeader(wcsProps, "BP", wcs._sipBp);
221  }
222  wcsProps->add("CTYPE1", ctype1, "WCS Coordinate type");
223  wcsProps->add("CTYPE2", ctype2, "WCS Coordinate type");
224 
225  return wcsProps;
226 }
227 
228 template <class Archive>
230  Archive& ar, int const, dafBase::Persistable* persistable) {
231  execTrace("TanWcsFormatter delegateSerialize start");
232  afwImg::TanWcs* ip = dynamic_cast<afwImg::TanWcs*>(persistable);
233  if (ip == 0) {
234  throw LSST_EXCEPT(pexExcept::RuntimeError, "Serializing non-TanWcs");
235  }
236 
237  // Serialize most fields normally
238  ar & ip->_nWcsInfo & ip->_relax;
239  ar & ip->_wcsfixCtrl & ip->_wcshdrCtrl & ip->_nReject;
240  ar & ip->_coordSystem;
241 
242  ar & ip->_hasDistortion;
243 
244  if(ip->_hasDistortion) {
245  serializeEigenArray(ar, ip->_sipA);
246  serializeEigenArray(ar, ip->_sipAp);
247  serializeEigenArray(ar, ip->_sipB);
248  serializeEigenArray(ar, ip->_sipBp);
249  }
250 
251  // If we are loading, create the array of Wcs parameter structs
252  if (Archive::is_loading::value) {
253  ip->_wcsInfo =
254  reinterpret_cast<wcsprm*>(malloc(ip->_nWcsInfo * sizeof(wcsprm)));
255  }
256 
257  for (int i = 0; i < ip->_nWcsInfo; ++i) {
258  // If we are loading, initialize the struct first
259  if (Archive::is_loading::value) {
260  ip->_wcsInfo[i].flag = -1;
261  wcsini(1, 2, &(ip->_wcsInfo[i]));
262  }
263 
264  // Serialize only critical Wcs parameters
265  ar & ip->_wcsInfo[i].naxis;
266  ar & ip->_wcsInfo[i].equinox;
267  ar & ip->_wcsInfo[i].radesys;
268  ar & ip->_wcsInfo[i].crpix[0];
269  ar & ip->_wcsInfo[i].crpix[1];
270  ar & ip->_wcsInfo[i].cd[0];
271  ar & ip->_wcsInfo[i].cd[1];
272  ar & ip->_wcsInfo[i].cd[2];
273  ar & ip->_wcsInfo[i].cd[3];
274  ar & ip->_wcsInfo[i].crval[0];
275  ar & ip->_wcsInfo[i].crval[1];
276  ar & ip->_wcsInfo[i].cunit[0];
277  ar & ip->_wcsInfo[i].cunit[1];
278  ar & ip->_wcsInfo[i].ctype[0];
279  ar & ip->_wcsInfo[i].ctype[1];
280  ar & ip->_wcsInfo[i].altlin;
281 
282  // If we are loading, compute intermediate values given those above
283  if (Archive::is_loading::value) {
284  ip->_wcsInfo[i].flag = 0;
285  wcsset(&(ip->_wcsInfo[i]));
286  }
287  }
288  execTrace("TanWcsFormatter delegateSerialize end");
289 }
290 
292  boost::archive::text_oarchive & , int, dafBase::Persistable*);
294  boost::archive::text_iarchive & , int, dafBase::Persistable*);
296  boost::archive::binary_oarchive & , int, dafBase::Persistable*);
298  boost::archive::binary_iarchive & , int, dafBase::Persistable*);
299 
301  pexPolicy::Policy::Ptr policy) {
303 }
304 
int _wcsfixCtrl
Do potentially unsafe translations of non-standard unit strings? 0/1 = no/yes.
Definition: Wcs.h:394
Interface for TanWcsFormatter class.
static void delegateSerialize(Archive &ar, int const version, lsst::daf::base::Persistable *persistable)
boost::shared_ptr< Formatter > Ptr
Definition: Formatter.h:81
virtual void write(lsst::daf::base::Persistable const *persistable, lsst::daf::persistence::Storage::Ptr storage, lsst::daf::base::PropertySet::Ptr additionalData)
Eigen::MatrixXd _sipBp
Definition: TanWcs.h:205
coord::CoordSystem _coordSystem
Definition: Wcs.h:397
Class for storing ordered metadata with comments.
Definition: PropertyList.h:81
daf_persistence package header file
virtual boost::archive::text_iarchive & getIArchive(void)
boost::shared_ptr< PropertySet > Ptr
Definition: PropertySet.h:90
virtual std::string const & getPath(void)
Definition: FitsStorage.cc:109
definition of the Trace messaging facilities
virtual void update(lsst::daf::base::Persistable *persistable, lsst::daf::persistence::Storage::Ptr storage, lsst::daf::base::PropertySet::Ptr additionalData)
tbl::Key< int > wcs
#define __attribute__(x)
Eigen::MatrixXd _sipB
Definition: TanWcs.h:205
boost::shared_ptr< Policy > Ptr
Definition: Policy.h:172
limited backward compatibility to the DC2 run-time trace facilities
Definition: Trace.h:93
Construct a static instance of this helper class to register a Formatter subclass in the FormatterReg...
Definition: Formatter.h:138
Class for FITS file storage.
Definition: FitsStorage.h:52
table::Key< std::string > ctype2
Definition: Wcs.cc:1039
Interface for MaskedImageFormatter class.
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
struct wcsprm * _wcsInfo
Definition: Wcs.h:391
boost::shared_ptr< PropertyList > Ptr
Definition: PropertyList.h:84
Implementation of the WCS standard for the special case of the Gnomonic (tangent plane) projection...
Definition: TanWcs.h:68
Interface for PropertySetFormatter class.
virtual lsst::daf::base::Persistable * read(lsst::daf::persistence::Storage::Ptr storage, lsst::daf::base::PropertySet::Ptr additionalData)
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
virtual boost::archive::text_oarchive & getOArchive(void)
int _relax
Degree of permissiveness for wcspih (0 for strict); see wcshdr.h for details.
Definition: Wcs.h:393
tuple m
Definition: lsstimport.py:48
boost::shared_ptr< daf::base::PropertySet > readMetadata(std::string const &fileName, int hdu=0, bool strip=false)
Return the metadata (header entries) from a FITS file.
static lsst::daf::persistence::FormatterRegistration registration
boost::shared_ptr< Storage > Ptr
Definition: Storage.h:62
Eigen::MatrixXd _sipAp
Definition: TanWcs.h:205
Class for boost::serialization storage.
Definition: BoostStorage.h:58
Base class for all persistable classes.
Definition: Persistable.h:74
Eigen::MatrixXd _sipA
Definition: TanWcs.h:205
static lsst::daf::base::PropertyList::Ptr generatePropertySet(lsst::afw::image::TanWcs const &wcs)
#define EXEC_TRACE
static lsst::daf::persistence::Formatter::Ptr createInstance(lsst::pex::policy::Policy::Ptr policy)
int _wcshdrCtrl
Controls messages to stderr from wcshdr (0 for none); see wcshdr.h for details.
Definition: Wcs.h:395
ImageT val
Definition: CR.cc:154
TanWcsFormatter(lsst::pex::policy::Policy::Ptr policy)
table::Key< std::string > ctype1
Definition: Wcs.cc:1038
Class implementing persistence and retrieval for TanWcs objects.
Include files required for standard LSST Exception handling.
Interface for ImageFormatter class.
void serializeEigenArray(Archive &ar, Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > &m)
Provide a function to serialise an Eigen::Matrix so we can persist the SIP matrices.