FitsFile.cpp

Go to the documentation of this file.

 
/*
 * FitsFile.cpp
 *
 *  Created on: Jun 9, 2020
 *      Author: mschefer
 */
 
#include <assert.h>
 
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
 
#include <boost/algorithm/string/case_conv.hpp>
#include <boost/algorithm/string/trim.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/path.hpp>
#include <boost/regex.hpp>
 
#include "ElementsKernel/Exception.h"
#include "SEFramework/FITS/FitsFile.h"
 
namespace SourceXtractor {
 
static typename MetadataEntry::value_t valueAutoCast(const std::string& value) {
  boost::regex float_regex("^[-+]?(\\d*\\.?\\d+|\\d+\\.?\\d*)([eE][-+]?\\d+)?$");
  boost::regex int_regex("^[-+]?\\d+$");
 
  try {
    if (value.empty()) {
      return value;
    } else if (boost::regex_match(value, int_regex)) {
      return static_cast<int64_t>(std::stoll(value));
    } else if (boost::regex_match(value, float_regex)) {
      return std::stod(value);
    } else if (value.size() == 1) {
      return value.at(0);
    }
  } catch (...) {
  }
 
  // Single quotes are used as escape code of another single quote, and
  // the string starts and ends with single quotes.
  // We used to use boost::io::quoted here, but it seems that starting with 1.73 it
  // does not work well when the escape code and the delimiter are the same
  std::string unquoted;
  bool        escape = false;
  unquoted.reserve(value.size());
  for (auto i = value.begin(); i != value.end(); ++i) {
    if (*i == '\'' && !escape) {
      escape = true;
      // skip this char
    } else {
      escape = false;
      unquoted.push_back(*i);
    }
  }
  return unquoted;
}
 
static void close_fits(fitsfile* ptr) {
  if (ptr != nullptr) {
    int status = 0;
    fits_close_file(ptr, &status);
  }
}
 
FitsFile::FitsFile(const boost::filesystem::path& path, bool writeable)
    : m_path(path), m_is_writeable(writeable), m_fits_ptr(nullptr, close_fits) {
 
  open();
  loadInfo();
}
 
FitsFile::~FitsFile() {}
 
fitsfile* FitsFile::getFitsFilePtr() {
  return m_fits_ptr.get();
}
 
const std::vector<int>& FitsFile::getImageHdus() const {
  return m_image_hdus;
}
 
std::map<std::string, MetadataEntry>& FitsFile::getHDUHeaders(int hdu) {
  return m_headers.at(hdu - 1);
}
 
void FitsFile::open() {
  int status = 0;
  fitsfile* ptr = nullptr;
 
  // Open
  fits_open_image(&ptr, m_path.native().c_str(), m_is_writeable ? READWRITE : READONLY, &status);
  m_fits_ptr.reset(ptr);
  if (status != 0) {
    if (m_is_writeable) {
      // Create file if it does not exists
      status = 0;
      fits_create_file(&ptr, m_path.native().c_str(), &status);
      m_fits_ptr.reset(ptr);
    }
    if (status != 0) {
      throw Elements::Exception() << "Can't open FITS file: " << m_path << " status: " << status;
    }
  }
  assert(ptr->Fptr->open_count == 1);
}
 
void FitsFile::refresh() {
 
  // After modifying the FITS file, We need to close and reopen the file before we can query
  // infos about it again without cfitsio crashing
 
  int status = 0;
  fitsfile* ptr = nullptr;
 
  m_fits_ptr.reset(nullptr);
 
  fits_open_image(&ptr, m_path.native().c_str(), m_is_writeable ? READWRITE : READONLY, &status);
  if (status != 0) {
    throw Elements::Exception() << "Can't close and reopen FITS file: " << m_path << " status: " << status;
  }
  assert(ptr->Fptr->open_count == 1);
  m_fits_ptr.reset(ptr);
 
  loadInfo();
}
 
void FitsFile::loadInfo() {
  int status = 0;
 
  fitsfile* ptr = m_fits_ptr.get();
 
  // save current HDU (if the file is opened with advanced cfitsio syntax it might be set already
  int original_hdu = 0;
  fits_get_hdu_num(ptr, &original_hdu);
 
  // Number of HDU
  m_image_hdus.clear();
  int number_of_hdus = 0;
  if (fits_get_num_hdus(ptr, &number_of_hdus, &status) < 0) {
    throw Elements::Exception() << "Can't get the number of HDUs in FITS file: " << m_path;
  }
  m_headers.clear();
  m_headers.resize(number_of_hdus);
 
  // loop over HDUs to determine which ones are images
  int hdu_type = 0;
  for (int hdu_number = 1; hdu_number <= number_of_hdus; ++hdu_number) {
    fits_movabs_hdu(ptr, hdu_number, &hdu_type, &status);
    if (status != 0) {
      throw Elements::Exception() << "Can't switch HDUs while opening: " << m_path;
    }
 
    if (hdu_type == IMAGE_HDU) {
      int  bitpix, naxis;
      long naxes[3] = {1, 1, 1};
 
      fits_get_img_param(ptr, 3, &bitpix, &naxis, naxes, &status);
      if (status == 0 && (naxis == 2 || naxis == 3)) {
        m_image_hdus.emplace_back(hdu_number);
      }
    }
  }
 
  // go back to saved HDU
  fits_movabs_hdu(ptr, original_hdu, &hdu_type, &status);
 
  // load all FITS headers
  loadFitsHeader();
 
  // load optional .head file to override headers
  loadHeadFile();
}
 
static std::map<std::string, MetadataEntry> loadHeadersFromFits(fitsfile* fptr) {
  std::map<std::string, MetadataEntry> headers;
  char record[81];
  int keynum = 1, status = 0;
 
  fits_read_record(fptr, keynum, record, &status);
  while (status == 0 && strncmp(record, "END", 3) != 0) {
    static boost::regex regex("([^=]{8})=([^\\/]*)(\\/(.*))?");
    std::string         record_str(record);
 
    boost::smatch sub_matches;
    if (boost::regex_match(record_str, sub_matches, regex)) {
      auto keyword = boost::to_upper_copy(sub_matches[1].str());
      auto value   = sub_matches[2].str();
      auto comment = sub_matches[4].str();
      boost::trim(keyword);
      boost::trim(value);
      boost::trim(comment);
      headers.emplace(keyword, MetadataEntry{valueAutoCast(value), {{"comment", comment}}});
    }
    fits_read_record(fptr, ++keynum, record, &status);
  }
 
  return headers;
}
 
void FitsFile::loadFitsHeader() {
  int status = 0;
 
  // save current HDU (if the file is opened with advanced cfitsio syntax it might be set already)
  int original_hdu = 0;
  fits_get_hdu_num(m_fits_ptr.get(), &original_hdu);
 
  int hdu_type = 0;
  for (unsigned int i = 0; i < m_headers.size(); i++) {
    fits_movabs_hdu(m_fits_ptr.get(), i + 1, &hdu_type, &status);  // +1 hdus start at 1
 
    m_headers[i] = loadHeadersFromFits(m_fits_ptr.get());
  }
 
  // go back to saved HDU
  fits_movabs_hdu(m_fits_ptr.get(), original_hdu, &hdu_type, &status);
}
 
void FitsFile::loadHeadFile() {
  auto base_name = m_path.stem();
  base_name.replace_extension(".head");
  auto head_filename = m_path.parent_path() / base_name;
 
  if (!boost::filesystem::exists(head_filename)) {
    return;
  }
 
  auto          hdu_iter = m_image_hdus.begin();
  std::ifstream file;
 
  // open the file and check
  file.open(head_filename.native());
  if (!file.good() || !file.is_open()) {
    throw Elements::Exception() << "Cannot load ascii header file: " << head_filename;
  }
 
  while (file.good() && hdu_iter != m_image_hdus.end()) {
    int current_hdu = *hdu_iter;
 
    std::string line;
    std::getline(file, line);
 
    static boost::regex regex_blank_line("\\s*$");
    line = boost::regex_replace(line, regex_blank_line, std::string(""));
    if (line.size() == 0) {
      continue;
    }
 
    if (boost::to_upper_copy(line) == "END") {
      current_hdu = *(++hdu_iter);
    } else {
      static boost::regex regex("([^=]{1,8})=([^\\/]*)(\\/ (.*))?");
      boost::smatch       sub_matches;
      if (boost::regex_match(line, sub_matches, regex) && sub_matches.size() >= 3) {
        auto keyword = boost::to_upper_copy(sub_matches[1].str());
        auto value   = sub_matches[2].str();
        auto comment = sub_matches[4].str();
        boost::trim(keyword);
        boost::trim(value);
        boost::trim(comment);
        m_headers.at(current_hdu - 1)[keyword] = MetadataEntry{valueAutoCast(value), {{"comment", comment}}};
        ;
      }
    }
  }
}
 
std::vector<int> FitsFile::getDimensions(int hdu) const {
  int status = 0;
 
  // save current HDU
  int original_hdu = 0;
  fits_get_hdu_num(m_fits_ptr.get(), &original_hdu);
 
  // got to requested HDU
  int hdu_type = 0;
  fits_movabs_hdu(m_fits_ptr.get(), hdu, &hdu_type, &status);
 
  // get dimensions
  long naxes[3] = {1, 1, 1};
  int bitpix, naxis;
 
  fits_get_img_param(m_fits_ptr.get(), 3, &bitpix, &naxis, naxes, &status);
  if (status != 0 || (naxis != 2 && naxis != 3)) {
    throw Elements::Exception()
        << "Can't find 2D image or data cube in FITS file: " << m_path << "[" << hdu << "]";
  }
 
  std::vector<int> dims;
  dims.push_back(naxes[0]);
  dims.push_back(naxes[1]);
  if (naxis == 3) {
    dims.push_back(naxes[2]);
  }
 
  // go back to saved HDU
  fits_movabs_hdu(m_fits_ptr.get(), original_hdu, &hdu_type, &status);
 
  return dims;
}
 
 
}  // namespace SourceXtractor