SourceXtractorPlusPlus  0.8
Please provide a description of the project.
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FitsImageSource.cpp
Go to the documentation of this file.
1 
17 /*
18  * FitsImageSource.cpp
19  *
20  * Created on: Mar 21, 2018
21  * Author: mschefer
22  */
23 
26 #include <iomanip>
27 #include <fstream>
28 #include <boost/filesystem/path.hpp>
29 #include <boost/filesystem/operations.hpp>
30 #include <boost/regex.hpp>
31 #include <boost/algorithm/string/case_conv.hpp>
32 #include <boost/algorithm/string/trim.hpp>
33 
34 
35 namespace SourceXtractor {
36 
37 
40  char record[81];
41  int keynum = 1, status = 0;
42 
43  fits_read_record(fptr, keynum, record, &status);
44  while (status == 0 && strncmp(record, "END", 3) != 0) {
45  static boost::regex regex("(.+)=([^\\/]*)(.*)");
46  std::string record_str(record);
47 
48  boost::smatch sub_matches;
49  if (boost::regex_match(record_str, sub_matches, regex)) {
50  auto keyword = boost::to_upper_copy(sub_matches[1].str());
51  boost::trim(keyword);
52  auto i = headers.emplace(keyword, sub_matches[2]);
53  boost::trim(i.first->second);
54  }
55  fits_read_record(fptr, ++keynum, record, &status);
56  }
57 
58  return headers;
59 }
60 
61 template<typename T>
64  : m_filename(filename), m_manager(manager), m_hdu_number(hdu_number) {
65  int status = 0;
66  int bitpix, naxis;
67  long naxes[2] = {1,1};
68 
69  auto fptr = m_manager->getFitsFile(filename);
70 
71  if (m_hdu_number <= 0) {
72  int hdunum, hdutype;
73  fits_get_num_hdus(fptr, &hdunum, &status);
74  if (status != 0) {
75  throw Elements::Exception() << "Can't get the number of HDUs in the FITS file: " << filename;
76  }
77 
78  m_hdu_number = 0;
79  do {
80  ++m_hdu_number;
81  switchHdu(fptr, m_hdu_number);
82  fits_get_hdu_type(fptr, &hdutype, &status);
83  if (status == 0 && hdutype == IMAGE_HDU) {
84  fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status);
85  }
86  if (status != 0) {
87  throw Elements::Exception() << "Can't get the type of the HDU " << m_hdu_number << " in the FITS file: "
88  << filename;
89  }
90  } while (m_hdu_number <=hdunum && !(hdutype == IMAGE_HDU && naxis == 2));
91  if (hdutype != IMAGE_HDU) {
92  throw Elements::Exception() << "Can't find 2D image in FITS file: " << filename;
93  }
94  }
95  else {
96  switchHdu(fptr, m_hdu_number);
97  }
98 
99  fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status);
100  if (status != 0 || naxis != 2) {
101  throw Elements::Exception() << "Can't find 2D image in FITS file: " << filename << "[" << m_hdu_number << "]";
102  }
103 
104  m_width = naxes[0];
105  m_height = naxes[1];
106 
107  m_header = loadFitsHeader(fptr);
108  loadHeadFile();
109 }
110 
111 
112 template<typename T>
114  const std::shared_ptr<CoordinateSystem> coord_system,
116  : m_filename(filename), m_manager(manager), m_hdu_number(1) {
117  m_width = width;
118  m_height = height;
119 
120  {
121  // CREATE NEW FITS FILE
122  int status = 0;
123  fitsfile* fptr = nullptr;
124  fits_create_file(&fptr, ("!"+filename).c_str(), &status);
125  if (status != 0) {
126  throw Elements::Exception() << "Can't create or overwrite FITS file: " << filename;
127  }
128  assert(fptr != nullptr);
129 
130  long naxes[2] = {width, height};
131  fits_create_img(fptr, getImageType(), 2, naxes, &status);
132 
133  if (coord_system) {
134  auto headers = coord_system->getFitsHeaders();
135  for (auto& h : headers) {
136  std::ostringstream padded_key, serializer;
137  padded_key << std::setw(8) << std::left << h.first;
138 
139  serializer << padded_key.str() << "= " << std::left << std::setw(70) << h.second;
140  auto str = serializer.str();
141 
142  fits_update_card(fptr, padded_key.str().c_str(), str.c_str(), &status);
143  if (status != 0) {
144  char err_txt[31];
145  fits_get_errstatus(status, err_txt);
146  throw Elements::Exception() << "Couldn't write the WCS headers (" << err_txt << "): " << str;
147  }
148  }
149  }
150 
151  std::vector<T> buffer(width);
152  for (int i = 0; i<height; i++) {
153  long first_pixel[2] = {1, i+1};
154  fits_write_pix(fptr, getDataType(), first_pixel, width, &buffer[0], &status);
155  }
156  fits_close_file(fptr, &status);
157 
158  if (status != 0) {
159  throw Elements::Exception() << "Couldn't allocate space for new FITS file: " << filename;
160  }
161  }
162 
163  auto fptr = m_manager->getFitsFile(filename, true);
164  switchHdu(fptr, m_hdu_number);
165 }
166 
167 
168 template<typename T>
169 std::shared_ptr<ImageTile<T>> FitsImageSource<T>::getImageTile(int x, int y, int width, int height) const {
170  auto fptr = m_manager->getFitsFile(m_filename);
171  switchHdu(fptr, m_hdu_number);
172 
173  auto tile = std::make_shared<ImageTile<T>>((const_cast<FitsImageSource *>(this))->shared_from_this(), x, y, width,
174  height);
175 
176  long first_pixel[2] = {x + 1, y + 1};
177  long last_pixel[2] = {x + width, y + height};
178  long increment[2] = {1, 1};
179  int status = 0;
180 
181  auto image = tile->getImage();
182  fits_read_subset(fptr, getDataType(), first_pixel, last_pixel, increment,
183  nullptr, &image->getData()[0], nullptr, &status);
184  if (status != 0) {
185  throw Elements::Exception() << "Error reading image tile from FITS file.";
186  }
187 
188  return tile;
189 }
190 
191 
192 template<typename T>
194  auto fptr = m_manager->getFitsFile(m_filename, true);
195  switchHdu(fptr, m_hdu_number);
196 
197  auto image = tile.getImage();
198 
199  int x = tile.getPosX();
200  int y = tile.getPosY();
201  int width = image->getWidth();
202  int height = image->getHeight();
203 
204  long first_pixel[2] = {x + 1, y + 1};
205  long last_pixel[2] = {x + width, y + height};
206  int status = 0;
207 
208  fits_write_subset(fptr, getDataType(), first_pixel, last_pixel, &image->getData()[0], &status);
209  if (status != 0) {
210  throw Elements::Exception() << "Error saving image tile to FITS file.";
211  }
212 }
213 
214 
215 template<typename T>
216 void FitsImageSource<T>::switchHdu(fitsfile *fptr, int hdu_number) const {
217  int status = 0;
218  int hdu_type = 0;
219 
220  fits_movabs_hdu(fptr, hdu_number, &hdu_type, &status);
221 
222  if (status != 0) {
223  throw Elements::Exception() << "Could not switch to HDU # " << hdu_number << " in file " << m_filename;
224  }
225  if (hdu_type != IMAGE_HDU) {
226  throw Elements::Exception() << "Trying to access non-image HDU in file " << m_filename;
227  }
228 }
229 
230 
231 template<typename T>
233  auto filename = boost::filesystem::path(m_filename);
234  auto base_name = filename.stem();
235  base_name.replace_extension(".head");
236  auto head_filename = filename.parent_path() / base_name;
237 
238  int current_hdu = 1;
239 
240  if (boost::filesystem::exists(head_filename)) {
241  std::ifstream file;
242 
243  // open the file and check
244  file.open(head_filename.native());
245  if (!file.good() || !file.is_open()) {
246  throw Elements::Exception() << "Cannot load ascii header file: " << head_filename;
247  }
248 
249  while (file.good()) {
250  std::string line;
251  std::getline(file, line);
252 
253  line = boost::regex_replace(line, boost::regex("\\s*#.*"), std::string(""));
254  line = boost::regex_replace(line, boost::regex("\\s*$"), std::string(""));
255 
256  if (line.size() == 0) {
257  continue;
258  }
259 
260  if (boost::to_upper_copy(line) == "END") {
261  current_hdu++;
262  }
263  else if (current_hdu == m_hdu_number) {
264  boost::smatch sub_matches;
265  if (boost::regex_match(line, sub_matches, boost::regex("(.+)=(.+)")) && sub_matches.size() == 3) {
266  auto keyword = boost::to_upper_copy(sub_matches[1].str());
267  m_header[keyword] = sub_matches[2];
268  }
269  }
270  }
271  }
272 }
273 
274 
275 template <>
276 int FitsImageSource<double>::getDataType() const { return TDOUBLE; }
277 
278 template <>
279 int FitsImageSource<float>::getDataType() const { return TFLOAT; }
280 
281 template <>
282 int FitsImageSource<unsigned int>::getDataType() const { return TUINT; }
283 
284 template <>
285 int FitsImageSource<int>::getDataType() const { return TINT; }
286 
287 //FIXME what if compiled on 32bit system?
288 template <>
289 int FitsImageSource<long>::getDataType() const { return TLONGLONG; }
290 
291 template <>
292 int FitsImageSource<long long>::getDataType() const { return TLONGLONG; }
293 
294 template <>
295 int FitsImageSource<double>::getImageType() const { return DOUBLE_IMG; }
296 
297 template <>
298 int FitsImageSource<float>::getImageType() const { return FLOAT_IMG; }
299 
300 template <>
301 int FitsImageSource<unsigned int>::getImageType() const { return LONG_IMG; }
302 
303 template <>
304 int FitsImageSource<int>::getImageType() const { return LONG_IMG; }
305 
306 //FIXME what if compiled on 32bit system?
307 template <>
308 int FitsImageSource<long>::getImageType() const { return LONGLONG_IMG; }
309 
310 template <>
311 int FitsImageSource<long long>::getImageType() const { return LONGLONG_IMG; }
312 
313 //FIXME add missing types
314 
315 
317 template class FitsImageSource<SeFloat>;
318 template class FitsImageSource<int64_t>;
319 template class FitsImageSource<unsigned int>;
320 
321 }
T open(T...args)
std::shared_ptr< VectorImage< T > > getImage()
Definition: ImageTile.h:87
T good(T...args)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
T getline(T...args)
T left(T...args)
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
T setw(T...args)
STL class.
std::shared_ptr< FitsFileManager > m_manager
string filename
Definition: conf.py:63
std::shared_ptr< ImageTile< T > > getImageTile(int x, int y, int width, int height) const override
T size(T...args)
STL class.
T strncmp(T...args)
void switchHdu(fitsfile *fptr, int hdu_number) const
void saveTile(ImageTile< T > &tile) override
T emplace(T...args)
std::map< std::string, std::string > m_header
T is_open(T...args)
STL class.
FitsImageSource(const std::string &filename, int hdu_number=0, std::shared_ptr< FitsFileManager > manager=FitsFileManager::getInstance())
static std::map< std::string, std::string > loadFitsHeader(fitsfile *fptr)