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