SourceXtractorPlusPlus  0.10
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 
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  if (fits_get_hdu_num(fptr, &m_hdu_number) < 0) {
73  throw Elements::Exception() << "Can't get the active HDU from the FITS file: " << filename;
74  }
75  }
76  else {
77  switchHdu(fptr, m_hdu_number);
78  }
79 
80  fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status);
81  if (status != 0 || naxis != 2) {
82  throw Elements::Exception() << "Can't find 2D image in FITS file: " << filename << "[" << m_hdu_number << "]";
83  }
84 
85  m_width = naxes[0];
86  m_height = naxes[1];
87 
88  m_header = loadFitsHeader(fptr);
89  loadHeadFile();
90 }
91 
92 
93 template<typename T>
95  const std::shared_ptr<CoordinateSystem> coord_system,
97  : m_filename(filename), m_manager(manager), m_hdu_number(1) {
98  m_width = width;
99  m_height = height;
100 
101  {
102  // CREATE NEW FITS FILE
103  int status = 0;
104  fitsfile* fptr = nullptr;
105  fits_create_file(&fptr, ("!"+filename).c_str(), &status);
106  if (status != 0) {
107  throw Elements::Exception() << "Can't create or overwrite FITS file: " << filename;
108  }
109  assert(fptr != nullptr);
110 
111  long naxes[2] = {width, height};
112  fits_create_img(fptr, getImageType(), 2, naxes, &status);
113 
114  if (coord_system) {
115  auto headers = coord_system->getFitsHeaders();
116  for (auto& h : headers) {
117  std::ostringstream padded_key, serializer;
118  padded_key << std::setw(8) << std::left << h.first;
119 
120  serializer << padded_key.str() << "= " << std::left << std::setw(70) << h.second;
121  auto str = serializer.str();
122 
123  fits_update_card(fptr, padded_key.str().c_str(), str.c_str(), &status);
124  if (status != 0) {
125  char err_txt[31];
126  fits_get_errstatus(status, err_txt);
127  throw Elements::Exception() << "Couldn't write the WCS headers (" << err_txt << "): " << str;
128  }
129  }
130  }
131 
132  std::vector<T> buffer(width);
133  for (int i = 0; i<height; i++) {
134  long first_pixel[2] = {1, i+1};
135  fits_write_pix(fptr, getDataType(), first_pixel, width, &buffer[0], &status);
136  }
137  fits_close_file(fptr, &status);
138 
139  if (status != 0) {
140  throw Elements::Exception() << "Couldn't allocate space for new FITS file: " << filename;
141  }
142  }
143 
144  auto fptr = m_manager->getFitsFile(filename, true);
145  switchHdu(fptr, m_hdu_number);
146 }
147 
148 
149 template<typename T>
150 std::shared_ptr<ImageTile<T>> FitsImageSource<T>::getImageTile(int x, int y, int width, int height) const {
151  auto fptr = m_manager->getFitsFile(m_filename);
152  switchHdu(fptr, m_hdu_number);
153 
154  auto tile = std::make_shared<ImageTile<T>>((const_cast<FitsImageSource *>(this))->shared_from_this(), x, y, width,
155  height);
156 
157  long first_pixel[2] = {x + 1, y + 1};
158  long last_pixel[2] = {x + width, y + height};
159  long increment[2] = {1, 1};
160  int status = 0;
161 
162  auto image = tile->getImage();
163  fits_read_subset(fptr, getDataType(), first_pixel, last_pixel, increment,
164  nullptr, &image->getData()[0], nullptr, &status);
165  if (status != 0) {
166  throw Elements::Exception() << "Error reading image tile from FITS file.";
167  }
168 
169  return tile;
170 }
171 
172 
173 template<typename T>
175  auto fptr = m_manager->getFitsFile(m_filename, true);
176  switchHdu(fptr, m_hdu_number);
177 
178  auto image = tile.getImage();
179 
180  int x = tile.getPosX();
181  int y = tile.getPosY();
182  int width = image->getWidth();
183  int height = image->getHeight();
184 
185  long first_pixel[2] = {x + 1, y + 1};
186  long last_pixel[2] = {x + width, y + height};
187  int status = 0;
188 
189  fits_write_subset(fptr, getDataType(), first_pixel, last_pixel, &image->getData()[0], &status);
190  if (status != 0) {
191  throw Elements::Exception() << "Error saving image tile to FITS file.";
192  }
193 }
194 
195 
196 template<typename T>
197 void FitsImageSource<T>::switchHdu(fitsfile *fptr, int hdu_number) const {
198  int status = 0;
199  int hdu_type = 0;
200 
201  fits_movabs_hdu(fptr, hdu_number, &hdu_type, &status);
202 
203  if (status != 0) {
204  throw Elements::Exception() << "Could not switch to HDU # " << hdu_number << " in file " << m_filename;
205  }
206  if (hdu_type != IMAGE_HDU) {
207  throw Elements::Exception() << "Trying to access non-image HDU in file " << m_filename;
208  }
209 }
210 
211 
212 template<typename T>
214  auto filename = boost::filesystem::path(m_filename);
215  auto base_name = filename.stem();
216  base_name.replace_extension(".head");
217  auto head_filename = filename.parent_path() / base_name;
218 
219  int current_hdu = 1;
220 
221  if (boost::filesystem::exists(head_filename)) {
222  std::ifstream file;
223 
224  // open the file and check
225  file.open(head_filename.native());
226  if (!file.good() || !file.is_open()) {
227  throw Elements::Exception() << "Cannot load ascii header file: " << head_filename;
228  }
229 
230  while (file.good()) {
231  std::string line;
232  std::getline(file, line);
233 
234  line = boost::regex_replace(line, boost::regex("\\s*#.*"), std::string(""));
235  line = boost::regex_replace(line, boost::regex("\\s*$"), std::string(""));
236 
237  if (line.size() == 0) {
238  continue;
239  }
240 
241  if (boost::to_upper_copy(line) == "END") {
242  current_hdu++;
243  }
244  else if (current_hdu == m_hdu_number) {
245  boost::smatch sub_matches;
246  if (boost::regex_match(line, sub_matches, boost::regex("(.+)=(.+)")) && sub_matches.size() == 3) {
247  auto keyword = boost::to_upper_copy(sub_matches[1].str());
248  m_header[keyword] = sub_matches[2];
249  }
250  }
251  }
252  }
253 }
254 
255 
256 template <>
257 int FitsImageSource<double>::getDataType() const { return TDOUBLE; }
258 
259 template <>
260 int FitsImageSource<float>::getDataType() const { return TFLOAT; }
261 
262 template <>
263 int FitsImageSource<unsigned int>::getDataType() const { return TUINT; }
264 
265 template <>
266 int FitsImageSource<int>::getDataType() const { return TINT; }
267 
268 //FIXME what if compiled on 32bit system?
269 template <>
270 int FitsImageSource<long>::getDataType() const { return TLONGLONG; }
271 
272 template <>
273 int FitsImageSource<long long>::getDataType() const { return TLONGLONG; }
274 
275 template <>
276 int FitsImageSource<double>::getImageType() const { return DOUBLE_IMG; }
277 
278 template <>
279 int FitsImageSource<float>::getImageType() const { return FLOAT_IMG; }
280 
281 template <>
282 int FitsImageSource<unsigned int>::getImageType() const { return LONG_IMG; }
283 
284 template <>
285 int FitsImageSource<int>::getImageType() const { return LONG_IMG; }
286 
287 //FIXME what if compiled on 32bit system?
288 template <>
289 int FitsImageSource<long>::getImageType() const { return LONGLONG_IMG; }
290 
291 template <>
292 int FitsImageSource<long long>::getImageType() const { return LONGLONG_IMG; }
293 
294 //FIXME add missing types
295 
296 
298 template class FitsImageSource<SeFloat>;
299 template class FitsImageSource<int64_t>;
300 template class FitsImageSource<unsigned int>;
301 
302 }
STL class.
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.
void switchHdu(fitsfile *fptr, int hdu_number) const
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 str(T... args)
T size(T... args)
STL class.
T strncmp(T... args)
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)