SourceXtractorPlusPlus  0.10
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
37 #include "ElementsKernel/Logging.h"
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);
95 }
96 
98 }
99 
101  // wcsprm is in/out, since its member lin is modified by wcsp2s
102  wcslib::wcsprm wcs_copy = *m_wcs;
103  wcs_copy.lin.flag = -1;
104  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
105  linset(&wcs_copy.lin);
106 
107  // +1 as fits standard coordinates start at 1
108  double pc_array[2] {image_coordinate.m_x + 1, image_coordinate.m_y + 1};
109 
110  double ic_array[2] {0, 0};
111  double wc_array[2] {0, 0};
112  double phi, theta;
113 
114  int status = 0;
115  wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
116  linfree(&wcs_copy.lin);
117 
118  return WorldCoordinate(wc_array[0], wc_array[1]);
119 }
120 
122  // wcsprm is in/out, since its member lin is modified by wcss2p
123  wcslib::wcsprm wcs_copy = *m_wcs;
124  wcs_copy.lin.flag = -1;
125  safe_lincpy(true, &m_wcs->lin, &wcs_copy.lin);
126  linset(&wcs_copy.lin);
128  double pc_array[2] {0, 0};
129  double ic_array[2] {0, 0};
130  double wc_array[2] {world_coordinate.m_alpha, world_coordinate.m_delta};
131  double phi, theta;
132 
133  int status = 0;
134  wcss2p(&wcs_copy, 1, 1, wc_array, &phi, &theta, ic_array, pc_array, &status);
135  linfree(&wcs_copy.lin);
136 
137  return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1
138 }
139 
141  int nkeyrec;
142  char *raw_header;
143 
144  if (wcshdo(WCSHDO_none, m_wcs.get(), &nkeyrec, &raw_header) != 0) {
145  throw Elements::Exception() << "Failed to get the FITS headers for the WCS coordinate system";
146  }
147 
149  for (int i = 0; i < nkeyrec; ++i) {
150  char *hptr = &raw_header[80 * i];
151  std::string key(hptr, hptr + 8);
152  boost::trim(key);
153  std::string value(hptr + 9, hptr + 72);
154  boost::trim(value);
155  if (!key.empty()) {
156  headers.emplace(std::make_pair(key, value));
157  }
158  }
160  free(raw_header);
161  return headers;
162 }
163 
164 }
T empty(T...args)
static Elements::Logging logger
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
T free(T...args)
WCS(const std::string &fits_file_path, int hdu_number=1)
Definition: WCS.cpp:63
int wcsvfree(int *nwcs, struct wcsprm **wcs)
STL class.
WorldCoordinate imageToWorld(ImageCoordinate image_coordinate) const override
Definition: WCS.cpp:100
int wcsfree(struct wcsprm *wcs)
int wcshdo(int relax, struct wcsprm *wcs, int *nkeyrec, char **header)
T lock(T...args)
T make_pair(T...args)
std::map< std::string, std::string > getFitsHeaders() const override
Definition: WCS.cpp:140
#define WCSHDR_all
Definition: WCS.cpp:1064
virtual ~WCS()
Definition: WCS.cpp:97
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)
ImageCoordinate worldToImage(WorldCoordinate world_coordinate) const override
Definition: WCS.cpp:121
std::unique_ptr< wcslib::wcsprm, std::function< void(wcslib::wcsprm *)> > m_wcs
Definition: WCS.h:46
static Logging getLogger(const std::string &name="")
struct linprm lin
Definition: WCS.cpp:1539
int lincpy(int alloc, const struct linprm *linsrc, struct linprm *lindst)