SourceXtractorPlusPlus  0.13
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
WCS.cpp
Go to the documentation of this file.
1 
17 /*
18  * WCS.cpp
19  *
20  * Created on: Nov 17, 2016
21  * Author: mschefer
22  */
23 
24 #include <mutex>
25 
26 #include <boost/algorithm/string/trim.hpp>
27 
28 #include <fitsio.h>
29 
30 #include <wcslib/wcs.h>
31 #include <wcslib/wcshdr.h>
32 #include <wcslib/wcsfix.h>
33 #include <wcslib/wcsprintf.h>
34 #include <wcslib/getwcstab.h>
35 #if 0
36 #include <wcslib/dis.h>
37 #endif
38 
40 #include "ElementsKernel/Logging.h"
41 
43 
44 namespace SourceXtractor {
45 
46 static auto logger = Elements::Logging::getLogger("WCS");
47 
48 decltype(&lincpy) safe_lincpy = &lincpy;
49 
53 static void wcsRaiseOnParseError(int ret_code) {
54  switch (ret_code) {
55  case WCSHDRERR_SUCCESS:
56  break;
57  case WCSHDRERR_MEMORY:
58  throw Elements::Exception() << "wcslib failed to allocate memory";
59  case WCSHDRERR_PARSER:
60  throw Elements::Exception() << "wcslib failed to parse the FITS headers";
61  default:
62  throw Elements::Exception() << "Unexpected error when parsing the FITS WCS header: "
63  << ret_code;
64  }
65 }
66 
67 static void wcsLogErr(wcserr *err) {
68  if (err) {
69  logger.error() << err->file << ":" << err->line_no << " " << err->function;
70  logger.error() << err->msg;
71  }
72 }
73 
77 static void wcsRaiseOnTransformError(wcsprm *wcs, int ret_code) {
78  if (ret_code != WCSERR_SUCCESS) {
79  wcsLogErr(wcs->err);
80  wcsLogErr(wcs->lin.err);
81 #if 0
82  if (wcs->lin.dispre) {
83  wcsLogErr(wcs->lin.dispre->err);
84  }
85  throw Elements::Exception() << "WCS exception: " << wcs_errmsg[ret_code];
86 #endif
87  }
88 }
89 
93 static std::set<std::string> wcsExtractKeywords(const char *header, int number_of_records) {
94  std::set<std::string> result;
95  for (const char *p = header; *p != '\0' && number_of_records; --number_of_records, p += 80) {
96  size_t len = strcspn(p, "= ");
97  result.insert(std::string(p, len));
98  }
99  return result;
100 }
101 
105 static void wcsCheckHeaders(const wcsprm *wcs, const char *headers_str, int number_of_records) {
106  auto headers = wcsExtractKeywords(headers_str, number_of_records);
107 
108 #if 0
109  // SIP present, but not specified in CTYPE
110  // See https://github.com/astrorama/SourceXtractorPlusPlus/issues/254#issuecomment-765235869
111  if (wcs->lin.dispre) {
112  bool sip_used = false, sip_specified = false;
113  for (int i = 0; i < wcs->naxis; ++i) {
114  sip_used |= (strncmp(wcs->lin.dispre->dtype[i], "SIP", 3) == 0);
115  size_t ctype_len = strlen(wcs->ctype[i]);
116  sip_specified |= (strncmp(wcs->ctype[i] + ctype_len - 4, "-SIP", 4) == 0);
117  }
118  if (sip_used && !sip_specified) {
119  logger.warn() << "SIP coefficients present, but CTYPE has not the '-SIP' suffix";
120  logger.warn() << "SIP distortion will be applied, but this may not be desired";
121  logger.warn() << "To suppress this warning, explicitly add the '-SIP' suffix to the CTYPE,";
122  logger.warn() << "or remove the SIP distortion coefficients";
123  }
124  }
125 #endif
126 }
127 
131 static void wcsReportWarnings(const char *err_buffer) {
132  if (err_buffer[0]) {
133  logger.warn() << "WCS generated some errors in strict mode. This may be OK.";
134  logger.warn() << "Will run in relaxed mode.";
135  const char *eol;
136  do {
137  eol = strchrnul(err_buffer, '\n');
138  logger.warn() << std::string(err_buffer, eol - err_buffer);
139  err_buffer = eol + 1;
140  } while (*eol);
141  }
142 }
143 
149 static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst) {
150  static std::mutex lincpy_mutex;
151  std::lock_guard<std::mutex> lock(lincpy_mutex);
152 
153  return lincpy(alloc, linsrc, lindst);
154 }
155 
156 
157 WCS::WCS(const FitsImageSource& fits_image_source) : m_wcs(nullptr, nullptr) {
158  int number_of_records = 0;
159  auto fits_headers = fits_image_source.getFitsHeaders(number_of_records);
160 
161  init(&(*fits_headers)[0], number_of_records);
162 }
163 
164 WCS::WCS(const WCS& original) : m_wcs(nullptr, nullptr) {
165 
166  //FIXME Horrible hack: I couldn't figure out how to properly do a deep copy wcsprm so instead
167  // of making a copy, I use the ascii headers output from the original to recreate a new one
168 
169  int number_of_records;
170  char *raw_header;
171 
172  if (wcshdo(WCSHDO_none, original.m_wcs.get(), &number_of_records, &raw_header) != 0) {
173  throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system when copying WCS";
174  }
175 
176  init(raw_header, number_of_records);
177 
178  free(raw_header);
179 }
180 
181 
182 void WCS::init(char* headers, int number_of_records) {
183  wcserr_enable(1);
184 
185  int nreject = 0, nwcs = 0, nreject_strict = 0;
186  wcsprm* wcs;
187 
188  // Write warnings to a buffer
189  wcsprintf_set(nullptr);
190 
191 #if 0
192  // Do a first pass, in strict mode, and ignore the result.
193  // Log the reported errors as warnings
194  int ret = wcspih(headers, number_of_records, WCSHDR_strict, 2, &nreject_strict, &nwcs, &wcs);
196  wcsReportWarnings(wcsprintf_buf());
197 #endif
198 
199  // Do a second pass, in relaxed mode. We use the result.
200  int ret = wcspih(headers, number_of_records, WCSHDR_all, 0, &nreject, &nwcs, &wcs);
202  wcsset(wcs);
203 
204  // There are some things worth reporting about which WCS will not necessarily complain
205  wcsCheckHeaders(wcs, headers, number_of_records);
206 
207  m_wcs = decltype(m_wcs)(wcs, [nwcs](wcsprm* wcs) {
208  int nwcs_copy = nwcs;
209  wcsfree(wcs);
210  wcsvfree(&nwcs_copy, &wcs);
211  });
212 
213 #if 0
214  int wcsver[3];
215  wcslib_version(wcsver);
216  if (wcsver[0] < 5 || (wcsver[0] == 5 && wcsver[1] < 18)) {
217  logger.info() << "wcslib " << wcsver[0] << "." << wcsver[1]
218  << " is not fully thread safe, using wrapped lincpy call!";
220  }
221 #else
223 #endif
224 }
225 
226 
228 }
229 
231  // wcsprm is in/out, since its member lin is modified by wcsp2s
232  wcsprm wcs_copy = *m_wcs;
233  wcs_copy.lin.flag = -1;
234  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
235  linset(&wcs_copy.lin);
236 
237  // +1 as fits standard coordinates start at 1
238  double pc_array[2] {image_coordinate.m_x + 1, image_coordinate.m_y + 1};
239 
240  double ic_array[2] {0, 0};
241  double wc_array[2] {0, 0};
242  double phi, theta;
243 
244  int status = 0;
245  wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
246  int ret_val = wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
247  linfree(&wcs_copy.lin);
248  wcsRaiseOnTransformError(&wcs_copy, ret_val);
249 
250  return WorldCoordinate(wc_array[0], wc_array[1]);
251 }
252 
254  // wcsprm is in/out, since its member lin is modified by wcss2p
255  wcsprm wcs_copy = *m_wcs;
256  wcs_copy.lin.flag = -1;
257  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
258  linset(&wcs_copy.lin);
259 
260  double pc_array[2] {0, 0};
261  double ic_array[2] {0, 0};
262  double wc_array[2] {world_coordinate.m_alpha, world_coordinate.m_delta};
263  double phi, theta;
264 
265  int status = 0;
266  int ret_val = wcss2p(&wcs_copy, 1, 1, wc_array, &phi, &theta, ic_array, pc_array, &status);
267  linfree(&wcs_copy.lin);
268  wcsRaiseOnTransformError(&wcs_copy, ret_val);
269 
270  return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1
271 }
272 
274  int nkeyrec;
275  char *raw_header;
276 
277  if (wcshdo(WCSHDO_none, m_wcs.get(), &nkeyrec, &raw_header) != 0) {
278  throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system";
279  }
280 
282  for (int i = 0; i < nkeyrec; ++i) {
283  char *hptr = &raw_header[80 * i];
284  std::string key(hptr, hptr + 8);
285  boost::trim(key);
286  std::string value(hptr + 9, hptr + 72);
287  boost::trim(value);
288  if (!key.empty()) {
289  headers.emplace(std::make_pair(key, value));
290  }
291  }
292 
293  free(raw_header);
294  return headers;
295 }
296 
298  m_wcs->crpix[0] -= pc.m_x;
299  m_wcs->crpix[1] -= pc.m_y;
300 }
301 
302 
303 }
WCS(const FitsImageSource &fits_image_source)
Definition: WCS.cpp:157
T empty(T...args)
static auto logger
Definition: WCS.cpp:46
std::unique_ptr< std::vector< char > > getFitsHeaders(int &number_of_records) const
static void wcsCheckHeaders(const wcsprm *wcs, const char *headers_str, int number_of_records)
Definition: WCS.cpp:105
decltype(&lincpy) safe_lincpy
Definition: WCS.cpp:48
static void wcsRaiseOnParseError(int ret_code)
Definition: WCS.cpp:53
T free(T...args)
void init(char *headers, int number_of_records)
Definition: WCS.cpp:182
STL class.
WorldCoordinate imageToWorld(ImageCoordinate image_coordinate) const override
Definition: WCS.cpp:230
static std::set< std::string > wcsExtractKeywords(const char *header, int number_of_records)
Definition: WCS.cpp:93
static void wcsRaiseOnTransformError(wcsprm *wcs, int ret_code)
Definition: WCS.cpp:77
T strlen(T...args)
T strcspn(T...args)
T lock(T...args)
T make_pair(T...args)
std::map< std::string, std::string > getFitsHeaders() const override
Definition: WCS.cpp:273
A pixel coordinate made of two integers m_x and m_y.
std::unique_ptr< wcsprm, std::function< void(wcsprm *)> > m_wcs
Definition: WCS.h:54
virtual ~WCS()
Definition: WCS.cpp:227
static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst)
Definition: WCS.cpp:149
void addOffset(PixelCoordinate pc)
Definition: WCS.cpp:297
T insert(T...args)
STL class.
static void wcsReportWarnings(const char *err_buffer)
Definition: WCS.cpp:131
T strncmp(T...args)
static void wcsLogErr(wcserr *err)
Definition: WCS.cpp:67
T emplace(T...args)
ImageCoordinate worldToImage(WorldCoordinate world_coordinate) const override
Definition: WCS.cpp:253
static Logging getLogger(const std::string &name="")