LoadFITS.cpp

#include "MantidAPI/FileProperty.h"
#include "MantidAPI/MultipleFileProperty.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidDataHandling/LoadFITS.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/UnitFactory.h"

#include <boost/algorithm/string.hpp>

#include <Poco/BinaryReader.h>
#include <Poco/FileStream.h>
#include <Poco/Path.h>

using namespace Mantid::DataHandling;
using namespace Mantid::DataObjects;
using namespace Mantid::API;
using namespace Mantid::Kernel;

namespace Mantid {
namespace DataHandling {
// Register the algorithm into the AlgorithmFactory
DECLARE_FILELOADER_ALGORITHM(LoadFITS)

const std::string LoadFITS::g_BIT_DEPTH_NAME = "BitDepthName";
const std::string LoadFITS::g_ROTATION_NAME = "RotationName";
const std::string LoadFITS::g_AXIS_NAMES_NAME = "AxisNames";
const std::string LoadFITS::g_IMAGE_KEY_NAME = "ImageKeyName";
const std::string LoadFITS::g_HEADER_MAP_NAME = "HeaderMapFile";

const std::string LoadFITS::g_defaultImgType = "SAMPLE";

/**
 * Constructor. Just initialize everything to prevent issues.
 */
LoadFITS::LoadFITS()
    : m_headerScaleKey(), m_headerOffsetKey(), m_headerBitDepthKey(),
      m_headerRotationKey(), m_headerImageKeyKey(), m_headerAxisNameKeys(),
      m_mapFile(), m_pixelCount(0) {
  setupDefaultKeywordNames();
}

/**
* Return the confidence with with this algorithm can load the file
* @param descriptor A descriptor for the file
* @returns An integer specifying the confidence level. 0 indicates it will not
* be used
*/
int LoadFITS::confidence(Kernel::FileDescriptor &descriptor) const {
  // Should really improve this to check the file header (of first file at
  // least) to make sure it contains the fields wanted
  return (descriptor.extension() == ".fits" || descriptor.extension() == ".fit")
             ? 80
             : 0;
}

/**
* Initialise the algorithm. Declare properties which can be set before execution
* (input) or
* read from after the execution (output).
*/
void LoadFITS::init() {
  // Specify file extensions which can be associated with a FITS file.
  std::vector<std::string> exts, exts2;

  // Declare the Filename algorithm property. Mandatory. Sets the path to the
  // file to load.
  exts.clear();
  exts.push_back(".fits");
  exts.push_back(".fit");

  exts2.push_back(".*");

  declareProperty(new MultipleFileProperty("Filename", exts),
                  "The name of the input file (note that you can give "
                  "multiple file names separated by commas).");

  declareProperty(new API::WorkspaceProperty<API::Workspace>(
      "OutputWorkspace", "", Kernel::Direction::Output));

  declareProperty(
      new Kernel::PropertyWithValue<bool>("LoadAsRectImg", false,
                                          Kernel::Direction::Input),
      "If enabled (not by default), the output Workspace2D will have "
      "one histogram per row and one bin per pixel, such that a 2D "
      "color plot (color fill plot) will display an image.");

  auto zeroOrPosDbl = boost::make_shared<BoundedValidator<double>>();
  zeroOrPosDbl->setLower(0.0);
  declareProperty("FilterNoiseLevel", 0.0, zeroOrPosDbl,
                  "Threshold to remove noisy pixels. Try 50 for example.");

  auto posInt = boost::make_shared<BoundedValidator<int>>();
  posInt->setLower(1);
  declareProperty("BinSize", 1, posInt,
                  "Rebunch n*n on both axes, generating pixels with sums of "
                  "blocks of n by n original pixels.");

  auto posDbl = boost::make_shared<BoundedValidator<double>>();
  posDbl->setLower(std::numeric_limits<double>::epsilon());
  declareProperty("Scale", 80.0, posDbl, "Pixels per cm.",
                  Kernel::Direction::Input);

  declareProperty(
      new FileProperty(g_HEADER_MAP_NAME, "", FileProperty::OptionalDirectory,
                       "", Kernel::Direction::Input),
      "A file mapping header key names to non-standard names [line separated "
      "values in the format KEY=VALUE, e.g. BitDepthName=BITPIX] - do not use "
      "this if you want to keep compatibility with standard FITS files.");
}

/**
 * Execute the algorithm.
 */
void LoadFITS::exec() {
  // for non-standard headers, by default won't do anything
  mapHeaderKeys();

  std::string fName = getPropertyValue("Filename");
  std::vector<std::string> paths;
  boost::split(paths, fName, boost::is_any_of(","));

  int binSize = getProperty("BinSize");
  double noiseThresh = getProperty("FilterNoiseLevel");
  bool loadAsRectImg = getProperty("LoadAsRectImg");
  const std::string outWSName = getPropertyValue("OutputWorkspace");

  doLoadFiles(paths, outWSName, loadAsRectImg, binSize, noiseThresh);
}

/**
 * Load header(s) from FITS file(s) into FITSInfo header
 * struct(s). This is usually the first step when loading FITS files
 * into workspaces or anything else. In the simplest case, paths has
 * only one string and only one header struct is added in headers.
 *
 * @param paths File name(s)
 * @param headers Vector where to store the header struct(s)
 *
 * @throws std::runtime_error if issues are found in the headers
 */
void LoadFITS::doLoadHeaders(const std::vector<std::string> &paths,
                             std::vector<FITSInfo> &headers) {
  headers.resize(paths.size());

  for (size_t i = 0; i < paths.size(); ++i) {
    headers[i].extension = "";
    headers[i].filePath = paths[i];

    // Get various pieces of information from the file header which are used
    // to
    // create the workspace
    try {
      parseHeader(headers[i]);
    } catch (std::exception &e) {
      // Unable to parse the header, throw.
      throw std::runtime_error(
          "Severe problem found while parsing the header of "
          "this FITS file (" +
          paths[i] +
          "). This file may not be standard FITS. Error description: " +
          e.what());
    }

    // Get and convert specific standard header values which will are
    // needed to know how to load the data: BITPIX, NAXIS, NAXISi (where i =
    // 1..NAXIS, e.g. NAXIS2 for two axis).
    try {
      std::string tmpBitPix = headers[i].headerKeys[m_headerBitDepthKey];
      if (boost::contains(tmpBitPix, "-")) {
        boost::erase_all(tmpBitPix, "-");
        headers[i].isFloat = true;
      } else {
        headers[i].isFloat = false;
      }

      try {
        headers[i].bitsPerPixel = boost::lexical_cast<int>(tmpBitPix);
      } catch (std::exception &e) {
        throw std::runtime_error(
            "Coult not interpret the entry number of bits per pixel (" +
            m_headerBitDepthKey + ") as an integer. Error: " + e.what());
      }
      // Check that the files use bit depths of either 8, 16 or 32
      if (headers[i].bitsPerPixel != 8 && headers[i].bitsPerPixel != 16 &&
          headers[i].bitsPerPixel != 32 && headers[i].bitsPerPixel != 64)
        throw std::runtime_error(
            "This algorithm only supports 8, 16, 32 or 64 "
            "bits per pixel. The header of file '" +
            paths[i] + "' says that its bit depth is: " +
            boost::lexical_cast<std::string>(headers[i].bitsPerPixel));
    } catch (std::exception &e) {
      throw std::runtime_error(
          "Failed to process the '" + m_headerBitDepthKey +
          "' entry (bits per pixel) in the header of this file: " + paths[i] +
          ". Error description: " + e.what());
    }

    try {