MaxentTransformMultiFourier.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidAlgorithms/MaxEnt/MaxentTransformMultiFourier.h"
#include <gsl/gsl_fft_complex.h>
 
#include <iterator>
#include <stdexcept>
#include <utility>
 
namespace Mantid::Algorithms {
 
MaxentTransformMultiFourier::MaxentTransformMultiFourier(const MaxentSpaceComplex_sptr &dataSpace,
                                                         MaxentSpace_sptr imageSpace, size_t numSpec)
    : MaxentTransformFourier(dataSpace, std::move(imageSpace)), m_numSpec(numSpec), m_linearAdjustments(),
      m_constAdjustments() {}
 
std::vector<double> MaxentTransformMultiFourier::imageToData(const std::vector<double> &image) {
 
  std::vector<double> dataOneSpec = MaxentTransformFourier::imageToData(image);
 
  // Create concatenated copies of transformed data (one for each spectrum)
  std::vector<double> data;
  data.reserve(m_numSpec * dataOneSpec.size());
  for (size_t s = 0; s < m_numSpec; s++) {
    std::copy(dataOneSpec.begin(), dataOneSpec.end(), std::back_inserter(data));
  }
 
  // Apply adjustments (we assume there are sufficient adjustments supplied)
  double dataR = 0.123456789;
  double dataI = 0.987654321;
  if (!m_linearAdjustments.empty() && !m_constAdjustments.empty()) {
    for (size_t i = 0; i < data.size(); i++) {
      if (i % 2 == 0) { // Real part
        dataR = data[i];
        dataI = data[i + 1];
        data[i] = m_linearAdjustments[i] * dataR - m_linearAdjustments[i + 1] * dataI + m_constAdjustments[i];
      } else { // Imaginary part
        data[i] = m_linearAdjustments[i] * dataR + m_linearAdjustments[i - 1] * dataI + m_constAdjustments[i];
      }
    }
  } else if (!m_linearAdjustments.empty() && m_constAdjustments.empty()) {
    for (size_t i = 0; i < data.size(); i++) {
      if (i % 2 == 0) { // Real part
        dataR = data[i];
        dataI = data[i + 1];
        data[i] = m_linearAdjustments[i] * dataR - m_linearAdjustments[i + 1] * dataI;
      } else { // Imaginary part
        data[i] = m_linearAdjustments[i] * dataR + m_linearAdjustments[i - 1] * dataI;
      }
    }
  } else if (m_linearAdjustments.empty() && !m_constAdjustments.empty()) {
    for (size_t i = 0; i < data.size(); i++) {
      data[i] = data[i] + m_constAdjustments[i];
    }
  }
 
  return data;
}
 
std::vector<double> MaxentTransformMultiFourier::dataToImage(const std::vector<double> &data) {
 
  if (data.size() % m_numSpec)
    throw std::invalid_argument("Size of data vector must be a multiple of number of spectra.");
 
  // Sum the concatenated spectra in data
  size_t nData = data.size() / (m_numSpec);
  std::vector<double> dataSum(nData, 0.0);
  for (size_t s = 0; s < m_numSpec; s++) {
    for (size_t i = 0; i < nData; i++) {
      dataSum[i] += data[s * nData + i];
    }
  }
  // Then apply forward fourier transform to sum
  std::vector<double> image = MaxentTransformFourier::dataToImage(dataSum);
 
  return image;
}
 
void MaxentTransformMultiFourier::setAdjustments(const std::vector<double> &linAdj,
                                                 const std::vector<double> &constAdj) {
  m_linearAdjustments = linAdj;
  m_constAdjustments = constAdj;
}
 
} // namespace Mantid::Algorithms