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>

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

namespace {}

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_baseName(), m_pixelCount(0), m_rebin(1),
      m_noiseThresh(0.0), m_progress(NULL) {
  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;
}

/**
 * Sets several keyword names with default (and standard) values. You
 * don't want to change these unless you want to break compatibility
 * with the FITS standard.
 */
void LoadFITS::setupDefaultKeywordNames() {
  // Inits all the absolutely necessary keywords
  // standard headers (If SIMPLE=T)
  m_headerScaleKey = "BSCALE";
  m_headerOffsetKey = "BZERO";
  m_headerBitDepthKey = "BITPIX";
  m_headerImageKeyKey = "IMAGE_TYPE"; // This is a "HIERARCH Image/Type= "
  m_headerRotationKey = "ROTATION";

  m_headerNAxisNameKey = "NAXIS";
  m_headerAxisNameKeys.push_back("NAXIS1");
  m_headerAxisNameKeys.push_back("NAXIS2");

  m_mapFile = "";

  // extensions
  m_sampleRotation = "HIERARCH Sample/Tomo_Angle";
  m_imageType = "HIERARCH Image/Type";
}

/**
* 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(","));

  m_rebin = getProperty("BinSize");
  m_noiseThresh = getProperty("FilterNoiseLevel");
  bool loadAsRectImg = getProperty("LoadAsRectImg");
  const std::string outWSName = getPropertyValue("OutputWorkspace");

  doLoadFiles(paths, outWSName, loadAsRectImg);
}

/**
 * Create FITS file information for each file selected. Loads headers
 * and data from the files and creates and fills the output
 * workspace(s).
 *
 * @param paths File names as given in the algorithm input property
 *
 * @param outWSName name of the output (group) workspace to create
 *
 * @param loadAsRectImg Load files with 1 spectrum per row and 1 bin
 * per column, so a color fill plot displays the image
 *
 * @throw std::runtime_error when load fails (for example a memory
 * allocation issue, wrong rebin requested, etc.)
 */
void LoadFITS::doLoadFiles(const std::vector<std::string> &paths,
                           const std::string &outWSName, bool loadAsRectImg) {
  std::vector<FITSInfo> headers;

  doLoadHeaders(paths, headers);

  // No extension is set -> it's the standard format which we can parse.
  if (headers[0].numberOfAxis > 0)
    m_pixelCount += headers[0].axisPixelLengths[0];

  // Presumably 2 axis, but futureproofing.
  for (int i = 1; i < headers[0].numberOfAxis; ++i) {
    m_pixelCount *= headers[0].axisPixelLengths[i];
  }

  // Check consistency of m_rebin asap
  for (int i = 0; i < headers[0].numberOfAxis; ++i) {
    if (0 != (headers[0].axisPixelLengths[i] % m_rebin)) {
      throw std::runtime_error(
          "Cannot rebin this image in blocks of " +
          boost::lexical_cast<std::string>(m_rebin) + " x " +
          boost::lexical_cast<std::string>(m_rebin) +
          " pixels as requested because the size of dimension " +
          boost::lexical_cast<std::string>(i + 1) + " (" +
          boost::lexical_cast<std::string>(headers[0].axisPixelLengths[i]) +