SourceXtractorPlusPlus  0.13
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  throw Elements::Exception() << "WCS exception: " << wcs_errmsg[ret_code];
83  }
84 }
85 
89 static std::set<std::string> wcsExtractKeywords(const char *header, int number_of_records) {
90  std::set<std::string> result;
91  for (const char *p = header; *p != '\0' && number_of_records; --number_of_records, p += 80) {
92  size_t len = strcspn(p, "= ");
93  result.insert(std::string(p, len));
94  }
95  return result;
96 }
97 
101 static void wcsCheckHeaders(const wcsprm *wcs, const char *headers_str, int number_of_records) {
102  auto headers = wcsExtractKeywords(headers_str, number_of_records);
103 
104  // SIP present, but not specified in CTYPE
105  // See https://github.com/astrorama/SourceXtractorPlusPlus/issues/254#issuecomment-765235869
106  if (wcs->lin.dispre) {
107  bool sip_used = false, sip_specified = false;
108  for (int i = 0; i < wcs->naxis; ++i) {
109  sip_used |= (strncmp(wcs->lin.dispre->dtype[i], "SIP", 3) == 0);
110  size_t ctype_len = strlen(wcs->ctype[i]);
111  sip_specified |= (strncmp(wcs->ctype[i] + ctype_len - 4, "-SIP", 4) == 0);
112  }
113  if (sip_used && !sip_specified) {
114  logger.warn() << "SIP coefficients present, but CTYPE has not the '-SIP' suffix";
115  logger.warn() << "SIP distortion will be applied, but this may not be desired";
116  logger.warn() << "To suppress this warning, explicitly add the '-SIP' suffix to the CTYPE,";
117  logger.warn() << "or remove the SIP distortion coefficients";
118  }
119  }
120 }
121 
125 static void wcsReportWarnings(const char *err_buffer) {
126  if (err_buffer[0]) {
127  logger.warn() << "WCS generated some errors in strict mode. This may be OK.";
128  logger.warn() << "Will run in relaxed mode.";
129  const char *eol;
130  do {
131  eol = strchrnul(err_buffer, '\n');
132  logger.warn() << std::string(err_buffer, eol - err_buffer);
133  err_buffer = eol + 1;
134  } while (*eol);
135  }
136 }
137 
143 static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst) {
144  static std::mutex lincpy_mutex;
145  std::lock_guard<std::mutex> lock(lincpy_mutex);
146 
147  return lincpy(alloc, linsrc, lindst);
148 }
149 
150 
151 WCS::WCS(const FitsImageSource& fits_image_source) : m_wcs(nullptr, nullptr) {
152  int number_of_records = 0;
153  auto fits_headers = fits_image_source.getFitsHeaders(number_of_records);
154 
155  init(&(*fits_headers)[0], number_of_records);
156 }
157 
158 WCS::WCS(const WCS& original) : m_wcs(nullptr, nullptr) {
159 
160  //FIXME Horrible hack: I couldn't figure out how to properly do a deep copy wcsprm so instead
161  // of making a copy, I use the ascii headers output from the original to recreate a new one
162 
163  int number_of_records;
164  char *raw_header;
165 
166  if (wcshdo(WCSHDO_none, original.m_wcs.get(), &number_of_records, &raw_header) != 0) {
167  throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system when copying WCS";
168  }
169 
170  init(raw_header, number_of_records);
171 
172  free(raw_header);
173 }
174 
175 
176 void WCS::init(char* headers, int number_of_records) {
177  wcserr_enable(1);
178 
179  int nreject = 0, nwcs = 0, nreject_strict = 0;
180  wcsprm* wcs;
181 
182  // Write warnings to a buffer
183  wcsprintf_set(nullptr);
184 
185  // Do a first pass, in strict mode, and ignore the result.
186  // Log the reported errors as warnings
187  int ret = wcspih(headers, number_of_records, WCSHDR_strict, 2, &nreject_strict, &nwcs, &wcs);
189  wcsReportWarnings(wcsprintf_buf());
190 
191  // Do a second pass, in relaxed mode. We use the result.
192  ret = wcspih(headers, number_of_records, WCSHDR_all, 0, &nreject, &nwcs, &wcs);
194  wcsset(wcs);
195 
196  // There are some things worth reporting about which WCS will not necessarily complain
197  wcsCheckHeaders(wcs, headers, number_of_records);
198 
199  m_wcs = decltype(m_wcs)(wcs, [nwcs](wcsprm* wcs) {
200  int nwcs_copy = nwcs;
201  wcsfree(wcs);
202  wcsvfree(&nwcs_copy, &wcs);
203  });
204 
205  int wcsver[3];
206  wcslib_version(wcsver);
207  if (wcsver[0] < 5 || (wcsver[0] == 5 && wcsver[1] < 18)) {
208  logger.info() << "wcslib " << wcsver[0] << "." << wcsver[1]
209  << " is not fully thread safe, using wrapped lincpy call!";
211  }
212 }
213 
214 
216 }
217 
219  // wcsprm is in/out, since its member lin is modified by wcsp2s
220  wcsprm wcs_copy = *m_wcs;
221  wcs_copy.lin.flag = -1;
222  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
223  linset(&wcs_copy.lin);
224 
225  // +1 as fits standard coordinates start at 1
226  double pc_array[2] {image_coordinate.m_x + 1, image_coordinate.m_y + 1};
227 
228  double ic_array[2] {0, 0};
229  double wc_array[2] {0, 0};
230  double phi, theta;
231 
232  int status = 0;
233  wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
234  int ret_val = wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
235  linfree(&wcs_copy.lin);
236  wcsRaiseOnTransformError(&wcs_copy, ret_val);
237 
238  return WorldCoordinate(wc_array[0], wc_array[1]);
239 }
240 
242  // wcsprm is in/out, since its member lin is modified by wcss2p
243  wcsprm wcs_copy = *m_wcs;
244  wcs_copy.lin.flag = -1;
245  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
246  linset(&wcs_copy.lin);
247 
248  double pc_array[2] {0, 0};
249  double ic_array[2] {0, 0};
250  double wc_array[2] {world_coordinate.m_alpha, world_coordinate.m_delta};
251  double phi, theta;
252 
253  int status = 0;
254  int ret_val = wcss2p(&wcs_copy, 1, 1, wc_array, &phi, &theta, ic_array, pc_array, &status);
255  linfree(&wcs_copy.lin);
256  wcsRaiseOnTransformError(&wcs_copy, ret_val);
257 
258  return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1
259 }
260 
262  int nkeyrec;
263  char *raw_header;
264 
265  if (wcshdo(WCSHDO_none, m_wcs.get(), &nkeyrec, &raw_header) != 0) {
266  throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system";
267  }
268 
270  for (int i = 0; i < nkeyrec; ++i) {
271  char *hptr = &raw_header[80 * i];
272  std::string key(hptr, hptr + 8);
273  boost::trim(key);
274  std::string value(hptr + 9, hptr + 72);
275  boost::trim(value);
276  if (!key.empty()) {
277  headers.emplace(std::make_pair(key, value));
278  }
279  }
280 
281  free(raw_header);
282  return headers;
283 }
284 
286  m_wcs->crpix[0] -= pc.m_x;
287  m_wcs->crpix[1] -= pc.m_y;
288 }
289 
290 
291 }
SourceXtractor::FitsImageSource
Definition: FitsImageSource.h:42
std::lock
T lock(T... args)
SourceXtractor::PixelCoordinate
A pixel coordinate made of two integers m_x and m_y.
Definition: PixelCoordinate.h:37
SourceXtractor::WCS::getFitsHeaders
std::map< std::string, std::string > getFitsHeaders() const override
Definition: WCS.cpp:261
std::strlen
T strlen(T... args)
SourceXtractor::WCS::init
void init(char *headers, int number_of_records)
Definition: WCS.cpp:176
SourceXtractor::WCS::WCS
WCS(const FitsImageSource &fits_image_source)
Definition: WCS.cpp:151
std::string
STL class.
SourceXtractor::wcsCheckHeaders
static void wcsCheckHeaders(const wcsprm *wcs, const char *headers_str, int number_of_records)
Definition: WCS.cpp:101
SourceXtractor::wcsLogErr
static void wcsLogErr(wcserr *err)
Definition: WCS.cpp:65
std::map::emplace
T emplace(T... args)
std::unique_ptr::get
T get(T... args)
std::lock_guard
STL class.
SourceXtractor::wcsExtractKeywords
static std::set< std::string > wcsExtractKeywords(const char *header, int number_of_records)
Definition: WCS.cpp:89
SourceXtractor::safe_lincpy
decltype(&lincpy) safe_lincpy
Definition: WCS.cpp:46
SourceXtractor::WorldCoordinate::m_alpha
double m_alpha
Definition: CoordinateSystem.h:34
SourceXtractor::WorldCoordinate::m_delta
double m_delta
Definition: CoordinateSystem.h:34
SourceXtractor::FitsImageSource::getFitsHeaders
std::unique_ptr< std::vector< char > > getFitsHeaders(int &number_of_records) const
Definition: FitsImageSource.cpp:230
pc
constexpr double pc
SourceXtractor::WorldCoordinate
Definition: CoordinateSystem.h:33
SourceXtractor
Definition: Aperture.h:30
SourceXtractor::WCS::m_wcs
std::unique_ptr< wcsprm, std::function< void(wcsprm *)> > m_wcs
Definition: WCS.h:54
WCS.h
Exception.h
SourceXtractor::ImageCoordinate::m_y
double m_y
Definition: CoordinateSystem.h:43
Elements::Exception
SourceXtractor::wcsRaiseOnTransformError
static void wcsRaiseOnTransformError(wcsprm *wcs, int ret_code)
Definition: WCS.cpp:75
std::map< std::string, std::string >
SourceXtractor::wcsReportWarnings
static void wcsReportWarnings(const char *err_buffer)
Definition: WCS.cpp:125
SourceXtractor::ImageCoordinate
Definition: CoordinateSystem.h:42
SourceXtractor::logger
static auto logger
Definition: WCS.cpp:44
Elements::Logging::getLogger
static Logging getLogger(const std::string &name="")
SourceXtractor::WCS::imageToWorld
WorldCoordinate imageToWorld(ImageCoordinate image_coordinate) const override
Definition: WCS.cpp:218
std::strcspn
T strcspn(T... args)
std::strncmp
T strncmp(T... args)
SourceXtractor::ImageCoordinate::m_x
double m_x
Definition: CoordinateSystem.h:43
SourceXtractor::WCS
Definition: WCS.h:37
std::set::insert
T insert(T... args)
std::free
T free(T... args)
SourceXtractor::WCS::~WCS
virtual ~WCS()
Definition: WCS.cpp:215
std::string::empty
T empty(T... args)
SourceXtractor::wrapped_lincpy
static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst)
Definition: WCS.cpp:143
std::mutex
STL class.
Logging.h
std::make_pair
T make_pair(T... args)
SourceXtractor::WCS::worldToImage
ImageCoordinate worldToImage(WorldCoordinate world_coordinate) const override
Definition: WCS.cpp:241
SourceXtractor::WCS::addOffset
void addOffset(PixelCoordinate pc)
Definition: WCS.cpp:285
std::set
STL class.
SourceXtractor::wcsRaiseOnParseError
static void wcsRaiseOnParseError(int ret_code)
Definition: WCS.cpp:51