SourceXtractorPlusPlus  0.11
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 template<typename T>
47  : m_filename(filename), m_manager(manager), m_hdu_number(hdu_number) {
48  int status = 0;
49  int bitpix, naxis;
50  long naxes[2] = {1,1};
51 
52  m_fits_file = m_manager->getFitsFile(filename);
53  auto fptr = m_fits_file->getFitsFilePtr();
54 
55  if (m_hdu_number <= 0) {
56  if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) {
57  throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename;
58  }
59  }
60  else {
61  switchHdu(fptr, m_hdu_number);
62  }
63 
64  fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status);
65  if (status != 0 || naxis != 2) {
66  throw Elements::Exception() << "Can't find 2D image in FITS file: " << filename << "[" << m_hdu_number << "]";
67  }
68 
69  m_width = naxes[0];
70  m_height = naxes[1];
71 }
72 
73 
74 template<typename T>
76  const std::shared_ptr<CoordinateSystem> coord_system,
78  : m_filename(filename), m_manager(manager), m_hdu_number(1) {
79  m_width = width;
80  m_height = height;
81 
82  {
83  // CREATE NEW FITS FILE
84  int status = 0;
85  fitsfile* fptr = nullptr;
86  fits_create_file(&fptr, ("!"+filename).c_str(), &status);
87  if (status != 0) {
88  throw Elements::Exception() << "Can't create or overwrite FITS file: " << filename;
89  }
90  assert(fptr != nullptr);
91 
92  long naxes[2] = {width, height};
93  fits_create_img(fptr, getImageType(), 2, naxes, &status);
94 
95  if (coord_system) {
96  auto headers = coord_system->getFitsHeaders();
97  for (auto& h : headers) {
98  std::ostringstream padded_key, serializer;
99  padded_key << std::setw(8) << std::left << h.first;
100 
101  serializer << padded_key.str() << "= " << std::left << std::setw(70) << h.second;
102  auto str = serializer.str();
103 
104  fits_update_card(fptr, padded_key.str().c_str(), str.c_str(), &status);
105  if (status != 0) {
106  char err_txt[31];
107  fits_get_errstatus(status, err_txt);
108  throw Elements::Exception() << "Couldn't write the WCS headers (" << err_txt << "): " << str;
109  }
110  }
111  }
112 
113  std::vector<T> buffer(width);
114  for (int i = 0; i<height; i++) {
115  long first_pixel[2] = {1, i+1};
116  fits_write_pix(fptr, getDataType(), first_pixel, width, &buffer[0], &status);
117  }
118  fits_close_file(fptr, &status);
119 
120  if (status != 0) {
121  throw Elements::Exception() << "Couldn't allocate space for new FITS file: " << filename;
122  }
123  }
124 
125  // Reopens the newly created file through theFitsFileManager
126  m_fits_file = m_manager->getFitsFile(filename, true);
127  auto fptr = m_fits_file->getFitsFilePtr();
128  switchHdu(fptr, m_hdu_number);
129 }
130 
131 
132 template<typename T>
134  auto fptr = m_fits_file->getFitsFilePtr();
135  switchHdu(fptr, m_hdu_number);
136 
137  auto tile = std::make_shared<ImageTile<T>>((const_cast<FitsImageSource *>(this))->shared_from_this(), x, y, width,
138  height);
139 
140  long first_pixel[2] = {x + 1, y + 1};
141  long last_pixel[2] = {x + width, y + height};
142  long increment[2] = {1, 1};
143  int status = 0;
144 
145  auto image = tile->getImage();
146  fits_read_subset(fptr, getDataType(), first_pixel, last_pixel, increment,
147  nullptr, &image->getData()[0], nullptr, &status);
148  if (status != 0) {
149  throw Elements::Exception() << "Error reading image tile from FITS file.";
150  }
151 
152  return tile;
153 }
154 
155 
156 template<typename T>
158  auto fptr = m_fits_file->getFitsFilePtr();
159  switchHdu(fptr, m_hdu_number);
160 
161  auto image = tile.getImage();
162 
163  int x = tile.getPosX();
164  int y = tile.getPosY();
165  int width = image->getWidth();
166  int height = image->getHeight();
167 
168  long first_pixel[2] = {x + 1, y + 1};
169  long last_pixel[2] = {x + width, y + height};
170  int status = 0;
171 
172  fits_write_subset(fptr, getDataType(), first_pixel, last_pixel, &image->getData()[0], &status);
173  if (status != 0) {
174  throw Elements::Exception() << "Error saving image tile to FITS file.";
175  }
176 }
177 
178 
179 template<typename T>
180 void FitsImageSource<T>::switchHdu(fitsfile *fptr, int hdu_number) const {
181  int status = 0;
182  int hdu_type = 0;
183 
184  fits_movabs_hdu(fptr, hdu_number, &hdu_type, &status);
185 
186  if (status != 0) {
187  throw Elements::Exception() << "Could not switch to HDU # " << hdu_number << " in file " << m_filename;
188  }
189  if (hdu_type != IMAGE_HDU) {
190  throw Elements::Exception() << "Trying to access non-image HDU in file " << m_filename;
191  }
192 }
193 
194 
195 
196 template<typename T>
198  number_of_records = 0;
199  std::string records;
200 
201  auto& headers = getMetadata();
202  for (auto record : headers) {
203  auto key = record.first;
204 
205  std::string record_string(key);
206  if (record_string.size() > 8) {
207  throw Elements::Exception() << "FITS keyword longer than 8 characters";
208  } else if (record_string.size() < 8) {
209  record_string += std::string(8 - record_string.size(), ' ');
210  }
211 
212  if (headers.at(key).m_value.type() == typeid(std::string)) {
213  record_string += "= '" + VariantCast<std::string>(headers.at(key).m_value) + "'";
214  }
215  else {
216  record_string += "= " + VariantCast<std::string>(headers.at(key).m_value);
217  }
218 
219  if (record_string.size() > 80) {
220  throw Elements::Exception() << "FITS record longer than 80 characters";
221  }
222 
223 
224  if (record_string.size() < 80) {
225  record_string += std::string(80 - record_string.size(), ' ');
226  }
227 
228  records += record_string;
229  number_of_records++;
230  }
231 
232  std::string record_string("END");
233  record_string += std::string(80 - record_string.size(), ' ');
234  records += record_string;
235 
236  std::unique_ptr<std::vector<char>> buffer(new std::vector<char>(records.begin(), records.end()));
237  buffer->emplace_back(0);
238  return buffer;
239 }
240 
241 template <>
242 int FitsImageSource<double>::getDataType() const { return TDOUBLE; }
243 
244 template <>
245 int FitsImageSource<float>::getDataType() const { return TFLOAT; }
246 
247 template <>
248 int FitsImageSource<unsigned int>::getDataType() const { return TUINT; }
249 
250 template <>
251 int FitsImageSource<int>::getDataType() const { return TINT; }
252 
253 //FIXME what if compiled on 32bit system?
254 template <>
255 int FitsImageSource<long>::getDataType() const { return TLONGLONG; }
256 
257 template <>
258 int FitsImageSource<long long>::getDataType() const { return TLONGLONG; }
259 
260 template <>
261 int FitsImageSource<double>::getImageType() const { return DOUBLE_IMG; }
262 
263 template <>
264 int FitsImageSource<float>::getImageType() const { return FLOAT_IMG; }
265 
266 template <>
267 int FitsImageSource<unsigned int>::getImageType() const { return LONG_IMG; }
268 
269 template <>
270 int FitsImageSource<int>::getImageType() const { return LONG_IMG; }
271 
272 //FIXME what if compiled on 32bit system?
273 template <>
274 int FitsImageSource<long>::getImageType() const { return LONGLONG_IMG; }
275 
276 template <>
277 int FitsImageSource<long long>::getImageType() const { return LONGLONG_IMG; }
278 
279 //FIXME add missing types
280 
281 
283 template class FitsImageSource<SeFloat>;
284 template class FitsImageSource<int64_t>;
285 template class FitsImageSource<unsigned int>;
286 
287 }
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)
T end(T...args)
std::shared_ptr< FitsFile > m_fits_file
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.
STL class.
T begin(T...args)
void switchHdu(fitsfile *fptr, int hdu_number) const
void saveTile(ImageTile< T > &tile) override
FitsImageSource(const std::string &filename, int hdu_number=0, std::shared_ptr< FitsFileManager > manager=FitsFileManager::getInstance())
std::shared_ptr< VectorImage< T > > & getImage()
Definition: ImageTile.h:87