Workspace2D.cpp

#include "MantidDataObjects/Workspace2D.h"
#include "MantidKernel/Exception.h"
#include "MantidAPI/RefAxis.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/ISpectrum.h"
#include "MantidKernel/VectorHelper.h"

using Mantid::API::ISpectrum;
using Mantid::API::MantidImage;

namespace Mantid {
namespace DataObjects {
using std::size_t;

DECLARE_WORKSPACE(Workspace2D)

/// Constructor
Workspace2D::Workspace2D(): m_noVectors(0) {}

Workspace2D::Workspace2D(const Workspace2D &other)
    : MatrixWorkspace(other), m_noVectors(other.m_noVectors),
      m_monitorList(other.m_monitorList) {
  data.resize(m_noVectors);

  for (size_t i = 0; i < m_noVectors; ++i) {
    // Careful: data holds pointers to ISpectrum, but they point to Histogram1D.
    // There are copy constructors, but we would need a clone() function that is
    // aware of the polymorphism. Use cast + copy constructor for now.
    data[i] = new Histogram1D(*static_cast<Histogram1D *>(other.data[i]));
  }
}

/// Destructor
Workspace2D::~Workspace2D() {
// Clear out the memory
  for (size_t i = 0; i < m_noVectors; i++) {
    delete data[i];
  }
}

/**
 * Sets the size of the workspace and initializes arrays to zero
 *
 * @param NVectors :: The number of vectors/histograms/detectors in the
 * workspace
 *
 * @param XLength :: The number of X data points/bin boundaries in each vector
 * (must all be the same)
 *
 * @param YLength :: The number of data/error points in each vector (must all be
 * the same)
*/
void Workspace2D::init(const std::size_t &NVectors, const std::size_t &XLength,
                       const std::size_t &YLength) {
  m_noVectors = NVectors;
  data.resize(m_noVectors);

  MantidVecPtr t1, t2;
  t1.access().resize(XLength); // this call initializes array to zero
  t2.access().resize(YLength);
  for (size_t i = 0; i < m_noVectors; i++) {
    // Create the spectrum upon init
    Histogram1D *spec = new Histogram1D();
    data[i] = spec;
    // Set the data and X
    spec->setX(t1);
    spec->setDx(t1);
    spec->resetHasDx();
    // Y,E arrays populated
    spec->setData(t2, t2);
    // Default spectrum number = starts at 1, for workspace index 0.
    spec->setSpectrumNo(specid_t(i + 1));
    spec->setDetectorID(detid_t(i + 1));
  }

  // Add axes that reference the data
  m_axes.resize(2);
  m_axes[0] = new API::RefAxis(XLength, this);
  m_axes[1] = new API::SpectraAxis(this);
}

/** Gets the number of histograms
@return Integer
*/
size_t Workspace2D::getNumberHistograms() const {
  return getHistogramNumberHelper();
}

/// get pseudo size
size_t Workspace2D::size() const { return data.size() * blocksize(); }

/// get the size of each vector
size_t Workspace2D::blocksize() const {
  return (data.size() > 0)
             ? static_cast<ISpectrum const *>(data[0])->dataY().size()
             : 0;
}

/**
  * Copy the data (Y's) from an image to this workspace.
  * @param image :: An image to copy the data from.
  * @param start :: Startinf workspace indx to copy data to.
  * @param parallelExecution :: Should inner loop run as parallel operation
  */
void Workspace2D::setImageY(const MantidImage &image, size_t start,
                            bool parallelExecution) {
  MantidImage m;
  setImageYAndE(image, m, start, parallelExecution);
}

/**
  * Copy the data from an image to this workspace's errors.
  * @param image :: An image to copy the data from.
  * @param start :: Startinf workspace indx to copy data to.
  * @param parallelExecution :: Should inner loop run as parallel operation
  */
void Workspace2D::setImageE(const MantidImage &image, size_t start,
                            bool parallelExecution) {
  MantidImage m;
  setImageYAndE(m, image, start, parallelExecution);
}

/**
  * Copy the data from an image to the (Y's) and the errors for this
  * workspace.
  *
  * @param imageY :: An image to copy the data from.
  * @param imageE :: An image to copy the errors from.
  * @param start :: Startinf workspace indx to copy data to.
  *
  * @param loadAsRectImg :: load using one histogram per row and one
  * bin per column, instead of the default one histogram per pixel
  *
  * @param scale_1 :: scale factor for the X axis (norammly
  * representing the inverse of the pixel width or similar.
  *
  * @param parallelExecution :: Should inner loop run as parallel operation
  */
void Workspace2D::setImageYAndE(const API::MantidImage &imageY,
                                const API::MantidImage &imageE, size_t start,
                                bool loadAsRectImg, double scale_1,
                                bool parallelExecution) {
  UNUSED_ARG(parallelExecution) // for parallel for

  if (imageY.empty() && imageE.empty())
    return;
  if (imageY.empty() && imageE[0].empty())
    return;
  if (imageE.empty() && imageY[0].empty())
    return;

  if (!loadAsRectImg && blocksize() != 1) {
    throw std::runtime_error(
        "Cannot set image in workspace: a single bin workspace is "
        "required when initializing a workspace from an "
        "image using a histogram per pixel.");
  }

  size_t height;
  size_t width;
  if (!imageY.empty()) {
    height = imageY.size();
    width = imageY.front().size();
  } else {
    height = imageE.size();
    width = imageE.front().size();
  }
  size_t dataSize = width * height;

  if (start + dataSize > getNumberHistograms() * blocksize()) {
    throw std::runtime_error(
        "Cannot set image: image is bigger than workspace.");
  }

  if (!loadAsRectImg) {
    // 1 pixel - one spectrum
    PARALLEL_FOR_IF(parallelExecution)
    for (int i = 0; i < static_cast<int>(height); ++i) {

      const auto &rowY = imageY[i];
      const auto &rowE = imageE[i];
      size_t spec = start + static_cast<size_t>(i) * width;
      auto pE = rowE.begin();
      for (auto pY = rowY.begin(); pY != rowY.end() && pE != rowE.end();
           ++pY, ++pE, ++spec) {
        data[spec]->dataY()[0] = *pY;
        data[spec]->dataE()[0] = *pE;
      }
    }
  } else {

    if (height != (getNumberHistograms()))
      throw std::runtime_error(
          std::string("To load an image into a workspace with one spectrum per "
                      "row, then number of spectra (") +
          boost::lexical_cast<std::string>(getNumberHistograms()) +
          ") needs to be equal to the height (rows) of the image (" +
          boost::lexical_cast<std::string>(height) + ")");

    if (width != blocksize())
      throw std::runtime_error(
          std::string("To load an image into a workspace with one spectrum per "
                      "row, then number of bins (") +
          boost::lexical_cast<std::string>(blocksize()) +
          ") needs to be equal to the width (columns) of the image (" +
          boost::lexical_cast<std::string>(width) + ")");

    // one spectrum - one row
    PARALLEL_FOR_IF(parallelExecution)
    for (int i = 0; i < static_cast<int>(height); ++i) {

      const auto &rowY = imageY[i];
      const auto &rowE = imageE[i];
      data[i]->dataY() = rowY;
      data[i]->dataE() = rowE;
    }
    // X values. Set first spectrum and copy/propagate that one to all the other
    // spectra
    PARALLEL_FOR_IF(parallelExecution)
    for (int i = 0; i < static_cast<int>(width) + 1; ++i) {
      data[0]->dataX()[i] = i * scale_1;
    }
    PARALLEL_FOR_IF(parallelExecution)
    for (int i = 1; i < static_cast<int>(height); ++i) {
      data[i]->setX(data[0]->dataX());
    }
  }
}

//--------------------------------------------------------------------------------------------
/// Return the underlying ISpectrum ptr at the given workspace index.
ISpectrum *Workspace2D::getSpectrum(const size_t index) {
  if (index >= m_noVectors) {
    std::stringstream ss;
    ss << "Workspace2D::getSpectrum, histogram number " << index
       << " out of range " << m_noVectors;
    throw std::range_error(ss.str());
  }
  invalidateCommonBinsFlag();
  return data[index];
}

const ISpectrum *Workspace2D::getSpectrum(const size_t index) const {
  if (index >= m_noVectors) {
    std::stringstream ss;
    ss << "Workspace2D::getSpectrum, histogram number " << index
       << " out of range " << m_noVectors;
    throw std::range_error(ss.str());
  }
  return data[index];
}

//--------------------------------------------------------------------------------------------
/** Returns the number of histograms.
 *  For some reason Visual Studio couldn't deal with the main
 * getHistogramNumber() method
 *  being virtual so it now just calls this private (and virtual) method which
 * does the work.
 *  @return the number of histograms associated with the workspace
 */
size_t Workspace2D::getHistogramNumberHelper() const { return data.size(); }

//---------------------------------------------------------------------------
/** Rebin a particular spectrum to a new histogram bin boundaries.
 *
 * @param index :: workspace index to generate
 * @param X :: input X vector of the bin boundaries.
 * @param Y :: output vector to be filled with the Y data.
 * @param E :: output vector to be filled with the Error data (optionally)
 * @param skipError :: if true, the error vector is NOT calculated.
 *        CURRENTLY IGNORED, the Error is always calculated.
 */
void Workspace2D::generateHistogram(const std::size_t index, const MantidVec &X,
                                    MantidVec &Y, MantidVec &E,
                                    bool skipError) const {
  UNUSED_ARG(skipError);
  if (index >= this->m_noVectors)
    throw std::range_error(
        "Workspace2D::generateHistogram, histogram number out of range");
  // output data arrays are implicitly filled by function
  const ISpectrum *spec = this->getSpectrum(index);
  const MantidVec &currentX = spec->readX();
  const MantidVec &currentY = spec->readY();
  const MantidVec &currentE = spec->readE();
  if (X.size() <= 1)
    throw std::runtime_error(
        "Workspace2D::generateHistogram(): X vector must be at least length 2");
  Y.resize(X.size() - 1, 0);
  E.resize(X.size() - 1, 0);

  // Perform the rebin from the current bins to the new ones
  if (currentX.size() ==
      currentY.size()) // First, convert to bin boundaries if needed.  The
  {                    // VectorHelper::rebin, assumes bin boundaries, even if
    std::vector<double> histX; // it is a distribution!
    histX.resize(currentX.size() + 1);
    Mantid::Kernel::VectorHelper::convertToBinBoundary(currentX, histX);
    Mantid::Kernel::VectorHelper::rebin(histX, currentY, currentE, X, Y, E,
                                        this->isDistribution());
  } else // assume x_size = y_size + 1
  {
    Mantid::Kernel::VectorHelper::rebin(currentX, currentY, currentE, X, Y, E,
                                        this->isDistribution());
  }
}

} // namespace DataObjects
} // NamespaceMantid

///\cond TEMPLATE
template DLLExport class Mantid::API::WorkspaceProperty<
    Mantid::DataObjects::Workspace2D>;

namespace Mantid {
namespace Kernel {
template <>
DLLExport Mantid::DataObjects::Workspace2D_sptr
IPropertyManager::getValue<Mantid::DataObjects::Workspace2D_sptr>(
    const std::string &name) const {
  PropertyWithValue<Mantid::DataObjects::Workspace2D_sptr> *prop =
      dynamic_cast<PropertyWithValue<Mantid::DataObjects::Workspace2D_sptr> *>(
          getPointerToProperty(name));
  if (prop) {
    return *prop;
  } else {
    std::string message = "Attempt to assign property " + name +
                          " to incorrect type. Expected Workspace2D.";
    throw std::runtime_error(message);
  }
}

template <>
DLLExport Mantid::DataObjects::Workspace2D_const_sptr
IPropertyManager::getValue<Mantid::DataObjects::Workspace2D_const_sptr>(
    const std::string &name) const {
  PropertyWithValue<Mantid::DataObjects::Workspace2D_sptr> *prop =
      dynamic_cast<PropertyWithValue<Mantid::DataObjects::Workspace2D_sptr> *>(
          getPointerToProperty(name));
  if (prop) {
    return prop->operator()();
  } else {
    std::string message = "Attempt to assign property " + name +
                          " to incorrect type. Expected const Workspace2D.";
    throw std::runtime_error(message);
  }
}
} // namespace Kernel
} // namespace Mantid

///\endcond TEMPLATE