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