SourceXtractorPlusPlus  0.15
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 
24 #include <iomanip>
25 #include <fstream>
26 #include <numeric>
27 #include <string>
28 
29 #include <boost/filesystem/path.hpp>
30 #include <boost/filesystem/operations.hpp>
31 #include <boost/regex.hpp>
32 #include <boost/algorithm/string/case_conv.hpp>
33 #include <boost/algorithm/string/trim.hpp>
34 
36 
38 #include "SEUtils/VariantCast.h"
39 
41 
42 namespace SourceXtractor {
43 
44 namespace {
45 }
46 
48  ImageTile::ImageType image_type,
50  : m_filename(filename), m_handler(manager->getFileHandler(filename)), m_hdu_number(hdu_number) {
51  int status = 0;
52  int bitpix, naxis;
53  long naxes[2] = {1, 1};
54 
55  auto acc = m_handler->getAccessor<FitsFile>();
56  auto fptr = acc->m_fd.getFitsFilePtr();
57 
58  if (m_hdu_number <= 0) {
59  if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) {
60  throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename;
61  }
62  }
63  else {
64  switchHdu(fptr, m_hdu_number);
65  }
66 
67  fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status);
68  if (status != 0 || naxis != 2) {
69  throw Elements::Exception() << "Can't find 2D image in FITS file: " << filename << "[" << m_hdu_number << "]";
70  }
71 
72  m_width = naxes[0];
73  m_height = naxes[1];
74 
75  if (image_type < 0) {
76  switch (bitpix) {
77  case FLOAT_IMG:
79  break;
80  case DOUBLE_IMG:
82  break;
83  case LONG_IMG:
85  break;
86  case ULONG_IMG:
88  break;
89  case LONGLONG_IMG:
91  break;
92  default:
93  throw Elements::Exception() << "Unsupported FITS image type: " << bitpix;
94  }
95  }
96  else {
97  m_image_type = image_type;
98  }
99 }
100 
101 
103  ImageTile::ImageType image_type,
104  const std::shared_ptr<CoordinateSystem> coord_system, bool append,
105  bool empty_primary,
107  : m_filename(filename)
108  , m_handler(manager->getFileHandler(filename))
109  , m_width(width)
110  , m_height(height)
111  , m_image_type(image_type) {
112 
113  int status = 0;
114  fitsfile* fptr = nullptr;
115 
116  if (!append) {
117  // delete file if it exists already
118  boost::filesystem::remove(m_filename);
119  }
120 
121  {
122  auto acc = m_handler->getAccessor<FitsFile>(FileHandler::kWrite);
123  fptr = acc->m_fd.getFitsFilePtr();
124 
125  assert(fptr != nullptr);
126 
127  if (empty_primary && acc->m_fd.getImageHdus().size() == 0) {
128  fits_create_img(fptr, FLOAT_IMG, 0, nullptr, &status);
129  if (status != 0) {
130  throw Elements::Exception() << "Can't create empty hdu: " << filename << " status: " << status;
131  }
132  }
133 
134  long naxes[2] = {width, height};
135  fits_create_img(fptr, getImageType(), 2, naxes, &status);
136 
137  if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) {
138  throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename << " status: " << status;
139  }
140 
141  int hdutype = 0;
142  fits_movabs_hdu(fptr, m_hdu_number, &hdutype, &status);
143 
144  if (coord_system) {
145  auto headers = coord_system->getFitsHeaders();
146  for (auto& h : headers) {
147  std::ostringstream padded_key, serializer;
148  padded_key << std::setw(8) << std::left << h.first;
149 
150  serializer << padded_key.str() << "= " << std::left << std::setw(70) << h.second;
151  auto str = serializer.str();
152 
153  fits_update_card(fptr, padded_key.str().c_str(), str.c_str(), &status);
154  if (status != 0) {
155  char err_txt[31];
156  fits_get_errstatus(status, err_txt);
157  throw Elements::Exception() << "Couldn't write the WCS headers (" << err_txt << "): " << str << " status: " << status;
158  }
159  }
160  }
161 
162  std::vector<char> buffer(width * ImageTile::getTypeSize(image_type));
163  for (int i = 0; i < height; i++) {
164  long first_pixel[2] = {1, i + 1};
165  fits_write_pix(fptr, getDataType(), first_pixel, width, &buffer[0], &status);
166  }
167  fits_close_file(fptr, &status);
168 
169  if (status != 0) {
170  throw Elements::Exception() << "Couldn't allocate space for new FITS file: " << filename << " status: " << status;
171  }
172 
173  acc->m_fd.refresh(); // make sure changes to the file structure are taken into account
174 
175  }
176 
177  // work around for canonical name bug:
178  // our file's canonical name might be wrong if it didn't exist, so we need to make sure we get the correct handler
179  // after we created the file
180 
181  m_handler = nullptr;
182  m_handler = manager->getFileHandler(filename);
183 }
184 
186  auto acc = m_handler->getAccessor<FitsFile>();
187  auto fptr = acc->m_fd.getFitsFilePtr();
188  switchHdu(fptr, m_hdu_number);
189 
190  auto tile = ImageTile::create(m_image_type, x, y, width, height,
191  std::const_pointer_cast<ImageSource>(shared_from_this()));
192 
193  long first_pixel[2] = {x + 1, y + 1};
194  long last_pixel[2] = {x + width, y + height};
195  long increment[2] = {1, 1};
196  int status = 0;
197 
198  fits_read_subset(fptr, getDataType(), first_pixel, last_pixel, increment,
199  nullptr, tile->getDataPtr(), nullptr, &status);
200  if (status != 0) {
201  throw Elements::Exception() << "Error reading image tile from FITS file.";
202  }
203 
204  return tile;
205 }
206 
208  auto acc = m_handler->getAccessor<FitsFile>(FileHandler::kWrite);
209  auto fptr = acc->m_fd.getFitsFilePtr();
210  switchHdu(fptr, m_hdu_number);
211 
212  int x = tile.getPosX();
213  int y = tile.getPosY();
214  int width = tile.getWidth();
215  int height = tile.getHeight();
216 
217  long first_pixel[2] = {x + 1, y + 1};
218  long last_pixel[2] = {x + width, y + height};
219  int status = 0;
220 
221  fits_write_subset(fptr, getDataType(), first_pixel, last_pixel, tile.getDataPtr(), &status);
222  if (status != 0) {
223  throw Elements::Exception() << "Error saving image tile to FITS file.";
224  }
225  fits_flush_buffer(fptr, 0, &status);
226 }
227 
228 void FitsImageSource::switchHdu(fitsfile *fptr, int hdu_number) const {
229  int status = 0;
230  int hdu_type = 0;
231 
232  fits_movabs_hdu(fptr, hdu_number, &hdu_type, &status);
233 
234  if (status != 0) {
235  throw Elements::Exception() << "Could not switch to HDU # " << hdu_number << " in file "
236  << m_filename;
237  }
238  if (hdu_type != IMAGE_HDU) {
239  throw Elements::Exception() << "Trying to access non-image HDU in file " << m_filename;
240  }
241 }
242 
243 
245  number_of_records = 0;
246  std::string records;
247 
248  auto& headers = getMetadata();
249  for (auto record : headers) {
250  auto key = record.first;
251 
252  std::string record_string(key);
253  if (record_string.size() > 8) {
254  throw Elements::Exception() << "FITS keyword longer than 8 characters";
255  }
256  else if (record_string.size() < 8) {
257  record_string += std::string(8 - record_string.size(), ' ');
258  }
259 
260  if (headers.at(key).m_value.type() == typeid(std::string)) {
261  record_string += "= '" + VariantCast<std::string>(headers.at(key).m_value) + "'";
262  }
263  else {
264  record_string += "= " + VariantCast<std::string>(headers.at(key).m_value);
265  }
266 
267  if (record_string.size() > 80) {
268  throw Elements::Exception() << "FITS record longer than 80 characters";
269  }
270 
271 
272  if (record_string.size() < 80) {
273  record_string += std::string(80 - record_string.size(), ' ');
274  }
275 
276  records += record_string;
277  number_of_records++;
278  }
279 
280  std::string record_string("END");
281  record_string += std::string(80 - record_string.size(), ' ');
282  records += record_string;
283 
284  std::unique_ptr<std::vector<char>> buffer(new std::vector<char>(records.begin(), records.end()));
285  buffer->emplace_back(0);
286  return buffer;
287 }
288 
290  auto acc = m_handler->getAccessor<FitsFile>();
291  return acc->m_fd.getHDUHeaders(m_hdu_number);
292 }
293 
295  auto acc = m_handler->getAccessor<FitsFile>();
296  auto fptr = acc->m_fd.getFitsFilePtr();
297  switchHdu(fptr, m_hdu_number);
298 
299  std::ostringstream padded_key, serializer;
300  padded_key << std::setw(8) << std::left << key;
301 
302  serializer << padded_key.str() << "= " << std::left << std::setw(70)
303  << VariantCast<std::string>(value.m_value);
304  auto str = serializer.str();
305 
306  int status = 0;
307  fits_update_card(fptr, padded_key.str().c_str(), str.c_str(), &status);
308 
309  if (status != 0) {
310  char err_txt[31];
311  fits_get_errstatus(status, err_txt);
312  throw Elements::Exception() << "Couldn't write the metadata (" << err_txt << "): " << str;
313  }
314 
315  // update the metadata
316  acc->m_fd.getHDUHeaders(m_hdu_number)[key] = value;
317 }
318 
320  switch (m_image_type) {
321  default:
323  return TFLOAT;
325  return TDOUBLE;
326  case ImageTile::IntImage:
327  return TINT;
329  return TUINT;
331  return TLONGLONG;
332  }
333 }
334 
336  switch (m_image_type) {
337  default:
339  return FLOAT_IMG;
341  return DOUBLE_IMG;
342  case ImageTile::IntImage:
343  return LONG_IMG;
345  return ULONG_IMG;
347  return LONGLONG_IMG;
348  }
349 }
350 
351 //FIXME add missing types
352 
353 }
std::unique_ptr< std::vector< char > > getFitsHeaders(int &number_of_records) const
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
T left(T...args)
ImageTile::ImageType m_image_type
T end(T...args)
std::shared_ptr< ImageTile > getImageTile(int x, int y, int width, int height) const override
STL class.
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
T setw(T...args)
int getPosX() const
Definition: ImageTile.h:56
STL class.
FitsImageSource(const std::string &filename, int hdu_number=0, ImageTile::ImageType image_type=ImageTile::AutoType, std::shared_ptr< FileManager > manager=FileManager::getDefault())
string filename
Definition: conf.py:63
static size_t getTypeSize(ImageType image_type)
Definition: ImageTile.h:119
void setMetadata(std::string key, MetadataEntry value) override
represents access to a whole FITS file and handles loading and caching FITS headers ...
Definition: FitsFile.h:43
int getWidth() const
Definition: ImageTile.h:66
void switchHdu(fitsfile *fptr, int hdu_number) const
T size(T...args)
STL class.
STL class.
std::map< std::string, MetadataEntry > & getHDUHeaders(int hdu)
Definition: FitsFile.cpp:110
static std::shared_ptr< ImageTile > create(ImageType image_type, int x, int y, int width, int height, std::shared_ptr< ImageSource > source=nullptr)
Definition: ImageTile.cpp:24
T begin(T...args)
int getHeight() const
Definition: ImageTile.h:70
void saveTile(ImageTile &tile) override
fitsfile * getFitsFilePtr()
Definition: FitsFile.cpp:102
int getPosY() const
Definition: ImageTile.h:60
const std::map< std::string, MetadataEntry > getMetadata() const override
std::shared_ptr< FileHandler > m_handler
virtual void * getDataPtr()=0