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