SourceXtractorPlusPlus  0.14
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 = strchr(err_buffer, '\n');
138  if (eol) {
139  logger.warn() << std::string(err_buffer, eol - err_buffer);
140  err_buffer = eol + 1;
141  }
142  else {
143  logger.warn() << std::string(err_buffer);
144  }
145  } while (eol);
146  }
147 }
148 
154 static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst) {
155  static std::mutex lincpy_mutex;
156  std::lock_guard<std::mutex> lock(lincpy_mutex);
157 
158  return lincpy(alloc, linsrc, lindst);
159 }
160 
161 
162 WCS::WCS(const FitsImageSource& fits_image_source) : m_wcs(nullptr, nullptr) {
163  int number_of_records = 0;
164  auto fits_headers = fits_image_source.getFitsHeaders(number_of_records);
165 
166  init(&(*fits_headers)[0], number_of_records);
167 }
168 
169 WCS::WCS(const WCS& original) : m_wcs(nullptr, nullptr) {
170 
171  //FIXME Horrible hack: I couldn't figure out how to properly do a deep copy wcsprm so instead
172  // of making a copy, I use the ascii headers output from the original to recreate a new one
173 
174  int number_of_records;
175  char *raw_header;
176 
177  if (wcshdo(WCSHDO_none, original.m_wcs.get(), &number_of_records, &raw_header) != 0) {
178  throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system when copying WCS";
179  }
180 
181  init(raw_header, number_of_records);
182 
183  free(raw_header);
184 }
185 
186 
187 void WCS::init(char* headers, int number_of_records) {
188  wcserr_enable(1);
189 
190  int nreject = 0, nwcs = 0, nreject_strict = 0;
191  wcsprm* wcs;
192 
193  // Write warnings to a buffer
194  wcsprintf_set(nullptr);
195 
196 #if 0
197  // Do a first pass, in strict mode, and ignore the result.
198  // Log the reported errors as warnings
199  int ret = wcspih(headers, number_of_records, WCSHDR_strict, 2, &nreject_strict, &nwcs, &wcs);
201  wcsReportWarnings(wcsprintf_buf());
202 #endif
203 
204  // Do a second pass, in relaxed mode. We use the result.
205  int ret = wcspih(headers, number_of_records, WCSHDR_all, 0, &nreject, &nwcs, &wcs);
207  wcsset(wcs);
208 
209  // There are some things worth reporting about which WCS will not necessarily complain
210  wcsCheckHeaders(wcs, headers, number_of_records);
211 
212  m_wcs = decltype(m_wcs)(wcs, [nwcs](wcsprm* wcs) {
213  int nwcs_copy = nwcs;
214  wcsfree(wcs);
215  wcsvfree(&nwcs_copy, &wcs);
216  });
217 
218 #if 0
219  int wcsver[3];
220  wcslib_version(wcsver);
221  if (wcsver[0] < 5 || (wcsver[0] == 5 && wcsver[1] < 18)) {
222  logger.info() << "wcslib " << wcsver[0] << "." << wcsver[1]
223  << " is not fully thread safe, using wrapped lincpy call!";
225  }
226 #else
228 #endif
229 }
230 
231 
233 }
234 
236  // wcsprm is in/out, since its member lin is modified by wcsp2s
237  wcsprm wcs_copy = *m_wcs;
238  wcs_copy.lin.flag = -1;
239  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
240  linset(&wcs_copy.lin);
241 
242  // +1 as fits standard coordinates start at 1
243  double pc_array[2] {image_coordinate.m_x + 1, image_coordinate.m_y + 1};
244 
245  double ic_array[2] {0, 0};
246  double wc_array[2] {0, 0};
247  double phi, theta;
248 
249  int status = 0;
250  wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
251  int ret_val = wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
252  linfree(&wcs_copy.lin);
253  wcsRaiseOnTransformError(&wcs_copy, ret_val);
254 
255  return WorldCoordinate(wc_array[0], wc_array[1]);
256 }
257 
259  // wcsprm is in/out, since its member lin is modified by wcss2p
260  wcsprm wcs_copy = *m_wcs;
261  wcs_copy.lin.flag = -1;
262  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
263  linset(&wcs_copy.lin);
264 
265  double pc_array[2] {0, 0};
266  double ic_array[2] {0, 0};
267  double wc_array[2] {world_coordinate.m_alpha, world_coordinate.m_delta};
268  double phi, theta;
269 
270  int status = 0;
271  int ret_val = wcss2p(&wcs_copy, 1, 1, wc_array, &phi, &theta, ic_array, pc_array, &status);
272  linfree(&wcs_copy.lin);
273  wcsRaiseOnTransformError(&wcs_copy, ret_val);
274 
275  return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1
276 }
277 
279  int nkeyrec;
280  char *raw_header;
281 
282  if (wcshdo(WCSHDO_none, m_wcs.get(), &nkeyrec, &raw_header) != 0) {
283  throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system";
284  }
285 
287  for (int i = 0; i < nkeyrec; ++i) {
288  char *hptr = &raw_header[80 * i];
289  std::string key(hptr, hptr + 8);
290  boost::trim(key);
291  std::string value(hptr + 9, hptr + 72);
292  boost::trim(value);
293  if (!key.empty()) {
294  headers.emplace(std::make_pair(key, value));
295  }
296  }
297 
298  free(raw_header);
299  return headers;
300 }
301 
303  m_wcs->crpix[0] -= pc.m_x;
304  m_wcs->crpix[1] -= pc.m_y;
305 }
306 
307 
308 }
WCS(const FitsImageSource &fits_image_source)
Definition: WCS.cpp:162
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:187
STL class.
WorldCoordinate imageToWorld(ImageCoordinate image_coordinate) const override
Definition: WCS.cpp:235
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:278
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:232
static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst)
Definition: WCS.cpp:154
void addOffset(PixelCoordinate pc)
Definition: WCS.cpp:302
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:258
static Logging getLogger(const std::string &name="")
T strchr(T...args)