SourceXtractorPlusPlus  0.16
Please provide a description of the project.
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 
168  if (status != 0) {
169  throw Elements::Exception() << "Couldn't allocate space for new FITS file: " << filename << " status: " << status;
170  }
171 
172  acc->m_fd.refresh(); // make sure changes to the file structure are taken into account
173 
174  }
175 
176  // work around for canonical name bug:
177  // our file's canonical name might be wrong if it didn't exist, so we need to make sure we get the correct handler
178  // after we created the file
179 
180  m_handler = nullptr;
181  m_handler = manager->getFileHandler(filename);
182 }
183 
184 std::shared_ptr<ImageTile> FitsImageSource::getImageTile(int x, int y, int width, int height) const {
185  auto acc = m_handler->getAccessor<FitsFile>();
186  auto fptr = acc->m_fd.getFitsFilePtr();
187  switchHdu(fptr, m_hdu_number);
188 
189  auto tile = ImageTile::create(m_image_type, x, y, width, height,
190  std::const_pointer_cast<ImageSource>(shared_from_this()));
191 
192  long first_pixel[2] = {x + 1, y + 1};
193  long last_pixel[2] = {x + width, y + height};
194  long increment[2] = {1, 1};
195  int status = 0;
196 
197  fits_read_subset(fptr, getDataType(), first_pixel, last_pixel, increment,
198  nullptr, tile->getDataPtr(), nullptr, &status);
199  if (status != 0) {
200  throw Elements::Exception() << "Error reading image tile from FITS file.";
201  }
202 
203  return tile;
204 }
205 
207  auto acc = m_handler->getAccessor<FitsFile>(FileHandler::kWrite);
208  auto fptr = acc->m_fd.getFitsFilePtr();
209  switchHdu(fptr, m_hdu_number);
210 
211  int x = tile.getPosX();
212  int y = tile.getPosY();
213  int width = tile.getWidth();
214  int height = tile.getHeight();
215 
216  long first_pixel[2] = {x + 1, y + 1};
217  long last_pixel[2] = {x + width, y + height};
218  int status = 0;
219 
220  fits_write_subset(fptr, getDataType(), first_pixel, last_pixel, tile.getDataPtr(), &status);
221  if (status != 0) {
222  throw Elements::Exception() << "Error saving image tile to FITS file.";
223  }
224  fits_flush_buffer(fptr, 0, &status);
225 }
226 
227 void FitsImageSource::switchHdu(fitsfile *fptr, int hdu_number) const {
228  int status = 0;
229  int hdu_type = 0;
230 
231  fits_movabs_hdu(fptr, hdu_number, &hdu_type, &status);
232 
233  if (status != 0) {
234  throw Elements::Exception() << "Could not switch to HDU # " << hdu_number << " in file "
235  << m_filename;
236  }
237  if (hdu_type != IMAGE_HDU) {
238  throw Elements::Exception() << "Trying to access non-image HDU in file " << m_filename;
239  }
240 }
241 
242 
244  number_of_records = 0;
245  std::string records;
246 
247  auto& headers = getMetadata();
248  for (auto record : headers) {
249  auto key = record.first;
250 
251  std::string record_string(key);
252  if (record_string.size() > 8) {
253  throw Elements::Exception() << "FITS keyword longer than 8 characters";
254  }
255  else if (record_string.size() < 8) {
256  record_string += std::string(8 - record_string.size(), ' ');
257  }
258 
259  if (headers.at(key).m_value.type() == typeid(std::string)) {
260  record_string += "= '" + VariantCast<std::string>(headers.at(key).m_value) + "'";
261  }
262  else {
263  record_string += "= " + VariantCast<std::string>(headers.at(key).m_value);
264  }
265 
266  if (record_string.size() > 80) {
267  throw Elements::Exception() << "FITS record longer than 80 characters";
268  }
269 
270 
271  if (record_string.size() < 80) {
272  record_string += std::string(80 - record_string.size(), ' ');
273  }
274 
275  records += record_string;
276  number_of_records++;
277  }
278 
279  std::string record_string("END");
280  record_string += std::string(80 - record_string.size(), ' ');
281  records += record_string;
282 
283  std::unique_ptr<std::vector<char>> buffer(new std::vector<char>(records.begin(), records.end()));
284  buffer->emplace_back(0);
285  return buffer;
286 }
287 
289  auto acc = m_handler->getAccessor<FitsFile>();
290  return acc->m_fd.getHDUHeaders(m_hdu_number);
291 }
292 
294  auto acc = m_handler->getAccessor<FitsFile>();
295  auto fptr = acc->m_fd.getFitsFilePtr();
296  switchHdu(fptr, m_hdu_number);
297 
298  std::ostringstream padded_key, serializer;
299  padded_key << std::setw(8) << std::left << key;
300 
301  serializer << padded_key.str() << "= " << std::left << std::setw(70)
302  << VariantCast<std::string>(value.m_value);
303  auto str = serializer.str();
304 
305  int status = 0;
306  fits_update_card(fptr, padded_key.str().c_str(), str.c_str(), &status);
307 
308  if (status != 0) {
309  char err_txt[31];
310  fits_get_errstatus(status, err_txt);
311  throw Elements::Exception() << "Couldn't write the metadata (" << err_txt << "): " << str;
312  }
313 
314  // update the metadata
315  acc->m_fd.getHDUHeaders(m_hdu_number)[key] = value;
316 }
317 
319  switch (m_image_type) {
320  default:
322  return TFLOAT;
324  return TDOUBLE;
325  case ImageTile::IntImage:
326  return TINT;
328  return TUINT;
330  return TLONGLONG;
331  }
332 }
333 
335  switch (m_image_type) {
336  default:
338  return FLOAT_IMG;
340  return DOUBLE_IMG;
341  case ImageTile::IntImage:
342  return LONG_IMG;
344  return ULONG_IMG;
346  return LONGLONG_IMG;
347  }
348 }
349 
350 //FIXME add missing types
351 
352 }
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > x
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
T begin(T... args)
represents access to a whole FITS file and handles loading and caching FITS headers
Definition: FitsFile.h:43
std::map< std::string, MetadataEntry > & getHDUHeaders(int hdu)
Definition: FitsFile.cpp:110
fitsfile * getFitsFilePtr()
Definition: FitsFile.cpp:102
const std::map< std::string, MetadataEntry > getMetadata() const override
std::shared_ptr< ImageTile > getImageTile(int x, int y, int width, int height) const override
ImageTile::ImageType m_image_type
void saveTile(ImageTile &tile) override
FitsImageSource(const std::string &filename, int hdu_number=0, ImageTile::ImageType image_type=ImageTile::AutoType, std::shared_ptr< FileManager > manager=FileManager::getDefault())
std::unique_ptr< std::vector< char > > getFitsHeaders(int &number_of_records) const
void setMetadata(std::string key, MetadataEntry value) override
std::shared_ptr< FileHandler > m_handler
void switchHdu(fitsfile *fptr, int hdu_number) const
static size_t getTypeSize(ImageType image_type)
Definition: ImageTile.h:117
virtual void * getDataPtr()=0
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 end(T... args)
T left(T... args)
string filename
Definition: conf.py:63
T setw(T... args)
T size(T... args)
T str(T... args)