SourceXtractorPlusPlus  0.13
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  : m_filename(filename), m_manager(manager), m_hdu_number(hdu_number) {
49  int status = 0;
50  int bitpix, naxis;
51  long naxes[2] = {1,1};
52 
53  m_fits_file = m_manager->getFitsFile(filename);
54  auto fptr = m_fits_file->getFitsFilePtr();
55 
56  if (m_hdu_number <= 0) {
57  if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) {
58  throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename;
59  }
60  }
61  else {
62  switchHdu(fptr, m_hdu_number);
63  }
64 
65  fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status);
66  if (status != 0 || naxis != 2) {
67  throw Elements::Exception() << "Can't find 2D image in FITS file: " << filename << "[" << m_hdu_number << "]";
68  }
69 
70  m_width = naxes[0];
71  m_height = naxes[1];
72 
73  if (image_type < 0) {
74  switch (bitpix) {
75  case FLOAT_IMG:
77  break;
78  case DOUBLE_IMG:
80  break;
81  case LONG_IMG:
83  break;
84  case ULONG_IMG:
86  break;
87  case LONGLONG_IMG:
89  break;
90  default:
91  throw Elements::Exception() << "Unsupported FITS image type: " << bitpix;
92  }
93  } else {
94  m_image_type = image_type;
95  }
96 }
97 
98 
100  const std::shared_ptr<CoordinateSystem> coord_system, bool append,
102  : m_filename(filename), m_manager(manager), m_width(width), m_height(height), m_image_type(image_type) {
103  {
104 
105  m_manager->closeAndPurgeFile(filename);
106 
107  int status = 0;
108  fitsfile* fptr = nullptr;
109 
110  if (append) {
111  fits_open_image(&fptr, filename.c_str(), READWRITE, &status);
112  if (status != 0) {
113  status = 0;
114  fits_create_file(&fptr, filename.c_str(), &status);
115  if (status != 0) {
116  throw Elements::Exception() << "Can't open or create FITS file to add HDU: " << filename;
117  }
118  }
119  } else {
120  // Overwrite file
121  fits_create_file(&fptr, ("!"+filename).c_str(), &status);
122  if (status != 0) {
123  throw Elements::Exception() << "Can't create or overwrite FITS file: " << filename;
124  }
125  }
126  assert(fptr != nullptr);
127 
128  long naxes[2] = {width, height};
129  fits_create_img(fptr, getImageType(), 2, naxes, &status);
130 
131  if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) {
132  throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename;
133  }
134 
135  int hdutype = 0;
136  fits_movabs_hdu(fptr, m_hdu_number, &hdutype, &status);
137 
138  if (coord_system) {
139  auto headers = coord_system->getFitsHeaders();
140  for (auto& h : headers) {
141  std::ostringstream padded_key, serializer;
142  padded_key << std::setw(8) << std::left << h.first;
143 
144  serializer << padded_key.str() << "= " << std::left << std::setw(70) << h.second;
145  auto str = serializer.str();
146 
147  fits_update_card(fptr, padded_key.str().c_str(), str.c_str(), &status);
148  if (status != 0) {
149  char err_txt[31];
150  fits_get_errstatus(status, err_txt);
151  throw Elements::Exception() << "Couldn't write the WCS headers (" << err_txt << "): " << str;
152  }
153  }
154  }
155 
156  std::vector<char> buffer(width * ImageTile::getTypeSize(image_type));
157  for (int i = 0; i<height; i++) {
158  long first_pixel[2] = {1, i+1};
159  fits_write_pix(fptr, getDataType(), first_pixel, width, &buffer[0], &status);
160  }
161  fits_close_file(fptr, &status);
162 
163  if (status != 0) {
164  throw Elements::Exception() << "Couldn't allocate space for new FITS file: " << filename;
165  }
166  }
167 
168  // Reopens the newly created file through theFitsFileManager
169  m_fits_file = m_manager->getFitsFile(filename, true);
170  auto fptr = m_fits_file->getFitsFilePtr();
171  switchHdu(fptr, m_hdu_number);
172 }
173 
175  auto fptr = m_fits_file->getFitsFilePtr();
176  switchHdu(fptr, m_hdu_number);
177 
178  auto tile = ImageTile::create(m_image_type, x, y, width, height,
179  std::const_pointer_cast<ImageSource>(shared_from_this()));
180 
181  long first_pixel[2] = {x + 1, y + 1};
182  long last_pixel[2] = {x + width, y + height};
183  long increment[2] = {1, 1};
184  int status = 0;
185 
186  fits_read_subset(fptr, getDataType(), first_pixel, last_pixel, increment,
187  nullptr, tile->getDataPtr(), nullptr, &status);
188  if (status != 0) {
189  throw Elements::Exception() << "Error reading image tile from FITS file.";
190  }
191 
192  return tile;
193 }
194 
196  auto fptr = m_fits_file->getFitsFilePtr();
197  switchHdu(fptr, m_hdu_number);
198 
199  int x = tile.getPosX();
200  int y = tile.getPosY();
201  int width = tile.getWidth();
202  int height = tile.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, tile.getDataPtr(), &status);
209  if (status != 0) {
210  throw Elements::Exception() << "Error saving image tile to FITS file.";
211  }
212 }
213 
214 void FitsImageSource::switchHdu(fitsfile *fptr, int hdu_number) const {
215  int status = 0;
216  int hdu_type = 0;
217 
218  fits_movabs_hdu(fptr, hdu_number, &hdu_type, &status);
219 
220  if (status != 0) {
221  throw Elements::Exception() << "Could not switch to HDU # " << hdu_number << " in file " << m_filename;
222  }
223  if (hdu_type != IMAGE_HDU) {
224  throw Elements::Exception() << "Trying to access non-image HDU in file " << m_filename;
225  }
226 }
227 
228 
229 
231  number_of_records = 0;
232  std::string records;
233 
234  auto& headers = getMetadata();
235  for (auto record : headers) {
236  auto key = record.first;
237 
238  std::string record_string(key);
239  if (record_string.size() > 8) {
240  throw Elements::Exception() << "FITS keyword longer than 8 characters";
241  } else if (record_string.size() < 8) {
242  record_string += std::string(8 - record_string.size(), ' ');
243  }
244 
245  if (headers.at(key).m_value.type() == typeid(std::string)) {
246  record_string += "= '" + VariantCast<std::string>(headers.at(key).m_value) + "'";
247  }
248  else {
249  record_string += "= " + VariantCast<std::string>(headers.at(key).m_value);
250  }
251 
252  if (record_string.size() > 80) {
253  throw Elements::Exception() << "FITS record longer than 80 characters";
254  }
255 
256 
257  if (record_string.size() < 80) {
258  record_string += std::string(80 - record_string.size(), ' ');
259  }
260 
261  records += record_string;
262  number_of_records++;
263  }
264 
265  std::string record_string("END");
266  record_string += std::string(80 - record_string.size(), ' ');
267  records += record_string;
268 
269  std::unique_ptr<std::vector<char>> buffer(new std::vector<char>(records.begin(), records.end()));
270  buffer->emplace_back(0);
271  return buffer;
272 }
273 
275  auto fptr = m_fits_file->getFitsFilePtr();
276  switchHdu(fptr, m_hdu_number);
277 
278  std::ostringstream padded_key, serializer;
279  padded_key << std::setw(8) << std::left << key;
280 
281  serializer << padded_key.str() << "= " << std::left << std::setw(70) << VariantCast<std::string>(value.m_value);
282  auto str = serializer.str();
283 
284  int status = 0;
285  fits_update_card(fptr, padded_key.str().c_str(), str.c_str(), &status);
286 
287  if (status != 0) {
288  char err_txt[31];
289  fits_get_errstatus(status, err_txt);
290  throw Elements::Exception() << "Couldn't write the metadata (" << err_txt << "): " << str;
291  }
292 
293  // update the metadata
294  m_fits_file->getHDUHeaders(m_hdu_number)[key] = value;
295 }
296 
298  switch (m_image_type) {
299  default:
301  return TFLOAT;
303  return TDOUBLE;
304  case ImageTile::IntImage:
305  return TINT;
307  return TUINT;
309  return TLONGLONG;
310  }
311 }
312 
314  switch (m_image_type) {
315  default:
317  return FLOAT_IMG;
319  return DOUBLE_IMG;
320  case ImageTile::IntImage:
321  return LONG_IMG;
323  return ULONG_IMG;
325  return LONGLONG_IMG;
326  }
327 }
328 
329 //FIXME add missing types
330 
331 }
std::unique_ptr< std::vector< char > > getFitsHeaders(int &number_of_records) const
const std::map< std::string, MetadataEntry > getMetadata() const override
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
std::shared_ptr< DependentParameter< std::shared_ptr< EngineParameter > > > y
std::shared_ptr< FitsFile > m_fits_file
T setw(T...args)
int getPosX() const
Definition: ImageTile.h:56
STL class.
std::shared_ptr< FitsFileManager > m_manager
string filename
Definition: conf.py:63
static size_t getTypeSize(ImageType image_type)
Definition: ImageTile.h:129
void setMetadata(std::string key, MetadataEntry value) override
int getWidth() const
Definition: ImageTile.h:66
void switchHdu(fitsfile *fptr, int hdu_number) const
T size(T...args)
STL class.
STL class.
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:96
T begin(T...args)
int getHeight() const
Definition: ImageTile.h:70
T c_str(T...args)
void saveTile(ImageTile &tile) override
int getPosY() const
Definition: ImageTile.h:60
virtual void * getDataPtr()=0
FitsImageSource(const std::string &filename, int hdu_number=0, ImageTile::ImageType image_type=ImageTile::AutoType, std::shared_ptr< FitsFileManager > manager=FitsFileManager::getInstance())