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 
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);
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 
105 }
106 
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);
113 
114  // +1 as fits standard coordinates start at 1
115  double pc_array[2] {image_coordinate.m_x + 1, image_coordinate.m_y + 1};
116 
117  double ic_array[2] {0, 0};
118  double wc_array[2] {0, 0};
119  double phi, theta;
120 
121  int status = 0;
122  wcsp2s(&wcs_copy, 1, 1, pc_array, ic_array, &phi, &theta, wc_array, &status);
123  linfree(&wcs_copy.lin);
124 
125  return WorldCoordinate(wc_array[0], wc_array[1]);
126 }
127 
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);
134 
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);
143 
144  return ImageCoordinate(pc_array[0] - 1, pc_array[1] - 1); // -1 as fits standard coordinates start at 1
145 }
146 
148  int nkeyrec;
149  char *raw_header;
150 
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  }
154 
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 }
std::lock
T lock(T... args)
SourceXtractor::WCS::getFitsHeaders
std::map< std::string, std::string > getFitsHeaders() const override
Definition: WCS.cpp:147
std::string
STL class.
SourceXtractor::WCS::m_wcs
std::unique_ptr< wcslib::wcsprm, std::function< void(wcslib::wcsprm *)> > m_wcs
Definition: WCS.h:46
std::map::emplace
T emplace(T... args)
std::unique_ptr::get
T get(T... args)
std::lock_guard
STL class.
SourceXtractor::safe_lincpy
decltype(&lincpy) safe_lincpy
Definition: WCS.cpp:48
SourceXtractor::WorldCoordinate::m_alpha
double m_alpha
Definition: CoordinateSystem.h:34
SourceXtractor::WorldCoordinate::m_delta
double m_delta
Definition: CoordinateSystem.h:34
SourceXtractor::WorldCoordinate
Definition: CoordinateSystem.h:33
SourceXtractor
Definition: Aperture.h:30
WCS.h
Exception.h
std::string::c_str
T c_str(T... args)
SourceXtractor::ImageCoordinate::m_y
double m_y
Definition: CoordinateSystem.h:43
Elements::Exception
std::map< std::string, std::string >
SourceXtractor::ImageCoordinate
Definition: CoordinateSystem.h:42
SourceXtractor::WCS::WCS
WCS(const std::string &fits_file_path, int hdu_number=1)
Definition: WCS.cpp:63
SourceXtractor::logger
static Elements::Logging logger
Definition: PluginManager.cpp:45
Elements::Logging::info
void info(const std::string &logMessage)
Elements::Logging::getLogger
static Logging getLogger(const std::string &name="")
SourceXtractor::WCS::imageToWorld
WorldCoordinate imageToWorld(ImageCoordinate image_coordinate) const override
Definition: WCS.cpp:107
wcslib
Definition: WCS.h:29
SourceXtractor::ImageCoordinate::m_x
double m_x
Definition: CoordinateSystem.h:43
std::free
T free(T... args)
SourceXtractor::WCS::~WCS
virtual ~WCS()
Definition: WCS.cpp:104
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:55
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:128