SourceXtractorPlusPlus  0.10
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 <fitsio.h>
25 
26 namespace wcslib {
27 
28 #include <wcslib/wcs.h>
29 #include <wcslib/wcshdr.h>
30 #include <wcslib/wcsfix.h>
31 #include <wcslib/wcsprintf.h>
32 #include <wcslib/getwcstab.h>
33 
34 }
35 
39 #include <boost/algorithm/string/trim.hpp>
40 #include <mutex>
41 
42 namespace SourceXtractor {
43 
44 static auto logger = Elements::Logging::getLogger("WCS");
45 
46 using namespace wcslib;
47 
48 decltype(&lincpy) safe_lincpy = &lincpy;
49 
55 static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst) {
56  static std::mutex lincpy_mutex;
57  std::lock_guard<std::mutex> lock(lincpy_mutex);
58 
59  return lincpy(alloc, linsrc, lindst);
60 }
61 
62 
63 WCS::WCS(const std::string& fits_file_path, int hdu_number) : m_wcs(nullptr, nullptr) {
64  fitsfile *fptr = NULL;
65  int status = 0;
66  fits_open_file(&fptr, fits_file_path.c_str(), READONLY, &status);
67 
68  int hdu_type;
69  fits_movabs_hdu(fptr, hdu_number, &hdu_type, &status);
70 
71  if (status != 0 || hdu_type != IMAGE_HDU) {
72  throw Elements::Exception() << "Can't read WCS information from " << fits_file_path << " HDU " << hdu_number;
73  }
74 
75  int nkeyrec;
76  char* header;
77  fits_hdr2str(fptr, 1, NULL, 0, &header, &nkeyrec, &status);
78 
79  if (hdu_type == IMAGE_HDU) {
80  int nreject = 0, nwcs = 0;
81  wcsprm* wcs;
82  wcspih(header, nkeyrec, WCSHDR_all, 0, &nreject, &nwcs, &wcs);
83  wcsset(wcs);
84 
85  m_wcs = decltype(m_wcs)(wcs, [nwcs](wcsprm* wcs) {
86  int nwcs_copy = nwcs;
87  wcsfree(wcs);
88  wcsvfree(&nwcs_copy, &wcs);
89  });
90  }
91 
92  free(header);
93  fits_close_file(fptr, &status);
94 
95  int wcsver[3];
96  wcslib_version(wcsver);
97  if (wcsver[0] < 5 || (wcsver[0] == 5 && wcsver[1] < 18)) {
98  logger.info() << "wcslib " << wcsver[0] << "." << wcsver[1]
99  << " is not fully thread safe, using wrapped lincpy call!";
101  }
102 }
103 
104 WCS::~WCS() {
105 }
107 WorldCoordinate WCS::imageToWorld(ImageCoordinate image_coordinate) const {
108  // wcsprm is in/out, since its member lin is modified by wcsp2s
109  wcslib::wcsprm wcs_copy = *m_wcs;
110  wcs_copy.lin.flag = -1;
111  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
112  linset(&wcs_copy.lin);
114  // +1 as fits standard coordinates start at 1
115  double pc_array[2] {image_coordinate.m_x + 1, image_coordinate.m_y + 1};
117  double ic_array[2] {0, 0};
118  double wc_array[2] {0, 0};
119  double phi, theta;
121  int status = 0;
122  wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
123  linfree(&wcs_copy.lin);
125  return WorldCoordinate(wc_array[0], wc_array[1]);
126 }
127 
128 ImageCoordinate WCS::worldToImage(WorldCoordinate world_coordinate) const {
129  // wcsprm is in/out, since its member lin is modified by wcss2p
130  wcslib::wcsprm wcs_copy = *m_wcs;
131  wcs_copy.lin.flag = -1;
132  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
133  linset(&wcs_copy.lin);
135  double pc_array[2] {0, 0};
136  double ic_array[2] {0, 0};
137  double wc_array[2] {world_coordinate.m_alpha, world_coordinate.m_delta};
138  double phi, theta;
139 
140  int status = 0;
141  wcss2p(&wcs_copy, 1, 1, wc_array, &phi, &theta, ic_array, pc_array, &status);
142  linfree(&wcs_copy.lin);
144  return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1
145 }
147 std::map<std::string, std::string> WCS::getFitsHeaders() const {
148  int nkeyrec;
149  char *raw_header;
151  if (wcshdo(WCSHDO_none, m_wcs.get(), &nkeyrec, &raw_header) != 0) {
152  throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system";
153  }
156  for (int i = 0; i < nkeyrec; ++i) {
157  char *hptr = &raw_header[80 * i];
158  std::string key(hptr, hptr + 8);
159  boost::trim(key);
160  std::string value(hptr + 9, hptr + 72);
161  boost::trim(value);
162  if (!key.empty()) {
163  headers.emplace(std::make_pair(key, value));
164  }
165  }
166 
167  free(raw_header);
168  return headers;
169 }
170 
171 }
static Elements::Logging logger
Definition: WCS.h:29
int wcspih(char *header, int nkeyrec, int relax, int ctrl, int *nreject, int *nwcs, struct wcsprm **wcs)
decltype(&lincpy) safe_lincpy
Definition: WCS.cpp:48
WCS(const std::string &fits_file_path, int hdu_number=1)
Definition: WCS.cpp:63
static auto logger
Definition: WCS.cpp:44
int wcsvfree(int *nwcs, struct wcsprm **wcs)
STL class.
int wcsfree(struct wcsprm *wcs)
int wcshdo(int relax, struct wcsprm *wcs, int *nkeyrec, char **header)
T make_pair(T... args)
const char * wcslib_version(int vers[3])
#define WCSHDR_all
Definition: WCS.cpp:1156
int wcsp2s(struct wcsprm *wcs, int ncoord, int nelem, const double pixcrd[], double imgcrd[], double phi[], double theta[], double world[], int stat[])
int linset(struct linprm *lin)
static int wrapped_lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst)
Definition: WCS.cpp:55
int wcss2p(struct wcsprm *wcs, int ncoord, int nelem, const double world[], double phi[], double theta[], double imgcrd[], double pixcrd[], int stat[])
T c_str(T... args)
T emplace(T... args)
int wcsset(struct wcsprm *wcs)
int linfree(struct linprm *lin)
static Logging getLogger(const std::string &name="")
struct linprm lin
Definition: WCS.cpp:1692
int lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst)