SourceXtractorPlusPlus  0.12
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 
37 #include "ElementsKernel/Logging.h"
38 
40 
41 namespace SourceXtractor {
42 
43 static auto logger = Elements::Logging::getLogger("WCS");
44 
45 decltype(&lincpy) safe_lincpy = &lincpy;
46 
52 static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst) {
53  static std::mutex lincpy_mutex;
54  std::lock_guard<std::mutex> lock(lincpy_mutex);
55 
56  return lincpy(alloc, linsrc, lindst);
57 }
58 
59 
60 WCS::WCS(const FitsImageSource<SeFloat>& fits_image_source) : m_wcs(nullptr, nullptr) {
61  int number_of_records = 0;
62  auto fits_headers = fits_image_source.getFitsHeaders(number_of_records);
63 
64  int nreject = 0, nwcs = 0;
65  wcsprm* wcs;
66  wcspih(&(*fits_headers)[0], number_of_records, WCSHDR_all, 0, &nreject, &nwcs, &wcs);
67  wcsset(wcs);
68 
69  m_wcs = decltype(m_wcs)(wcs, [nwcs](wcsprm* wcs) {
70  int nwcs_copy = nwcs;
71  wcsfree(wcs);
72  wcsvfree(&nwcs_copy, &wcs);
73  });
74 
75 
76  int wcsver[3];
77  wcslib_version(wcsver);
78  if (wcsver[0] < 5 || (wcsver[0] == 5 && wcsver[1] < 18)) {
79  logger.info() << "wcslib " << wcsver[0] << "." << wcsver[1]
80  << " is not fully thread safe, using wrapped lincpy call!";
82  }
83 }
84 
86 }
87 
89  // wcsprm is in/out, since its member lin is modified by wcsp2s
90  wcsprm wcs_copy = *m_wcs;
91  wcs_copy.lin.flag = -1;
92  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
93  linset(&wcs_copy.lin);
94 
95  // +1 as fits standard coordinates start at 1
96  double pc_array[2] {image_coordinate.m_x + 1, image_coordinate.m_y + 1};
97 
98  double ic_array[2] {0, 0};
99  double wc_array[2] {0, 0};
100  double phi, theta;
101 
102  int status = 0;
103  wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
104  int ret_val = wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
105  linfree(&wcs_copy.lin);
106  if (ret_val != 0) {
107  logger.error() << "wcslib's wcsp2s returned with error code: " << ret_val;
108  throw Elements::Exception() << "WCS exception";
109  }
110 
111  return WorldCoordinate(wc_array[0], wc_array[1]);
112 }
113 
115  // wcsprm is in/out, since its member lin is modified by wcss2p
116  wcsprm wcs_copy = *m_wcs;
117  wcs_copy.lin.flag = -1;
118  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
119  linset(&wcs_copy.lin);
120 
121  double pc_array[2] {0, 0};
122  double ic_array[2] {0, 0};
123  double wc_array[2] {world_coordinate.m_alpha, world_coordinate.m_delta};
124  double phi, theta;
125 
126  int status = 0;
127  int ret_val = wcss2p(&wcs_copy, 1, 1, wc_array, &phi, &theta, ic_array, pc_array, &status);
128  linfree(&wcs_copy.lin);
129  if (ret_val != 0) {
130  logger.error() << "wcslib's wcss2p returned with error code: " << ret_val;
131  throw Elements::Exception() << "WCS exception";
132  }
133 
134  return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1
135 }
136 
138  int nkeyrec;
139  char *raw_header;
140 
141  if (wcshdo(WCSHDO_none, m_wcs.get(), &nkeyrec, &raw_header) != 0) {
142  throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system";
143  }
144 
146  for (int i = 0; i < nkeyrec; ++i) {
147  char *hptr = &raw_header[80 * i];
148  std::string key(hptr, hptr + 8);
149  boost::trim(key);
150  std::string value(hptr + 9, hptr + 72);
151  boost::trim(value);
152  if (!key.empty()) {
153  headers.emplace(std::make_pair(key, value));
154  }
155  }
156 
157  free(raw_header);
158  return headers;
159 }
160 
161 }
T empty(T... args)
static Elements::Logging logger
WCS(const FitsImageSource< SeFloat > &fits_image_source)
Definition: WCS.cpp:60
decltype(&lincpy) safe_lincpy
Definition: WCS.cpp:45
std::unique_ptr< std::vector< char > > getFitsHeaders(int &number_of_records) const
void info(const std::string &logMessage)
T free(T... args)
STL class.
WorldCoordinate imageToWorld(ImageCoordinate image_coordinate) const override
Definition: WCS.cpp:88
T lock(T... args)
T make_pair(T... args)
std::map< std::string, std::string > getFitsHeaders() const override
Definition: WCS.cpp:137
std::unique_ptr< wcsprm, std::function< void(wcsprm *)> > m_wcs
Definition: WCS.h:48
virtual ~WCS()
Definition: WCS.cpp:85
static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst)
Definition: WCS.cpp:52
T get(T... args)
void error(const std::string &logMessage)
T emplace(T... args)
ImageCoordinate worldToImage(WorldCoordinate world_coordinate) const override
Definition: WCS.cpp:114
static Logging getLogger(const std::string &name="")