IDetector.cpp

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Device/Detector/IDetector.cpp
//! @brief     Implements common detector interface.
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************

#include "Device/Detector/IDetector.h"
#include "Base/Axis/FixedBinAxis.h"
#include "Base/Axis/Frame.h"
#include "Base/Axis/IAxis.h"
#include "Base/Const/Units.h"
#include "Base/Util/Assert.h"
#include "Device/Beam/Beam.h"
#include "Device/Detector/SimulationAreaIterator.h"
#include "Device/Mask/DetectorMask.h"
#include "Device/Mask/InfinitePlane.h"
#include "Device/Resolution/ConvolutionDetectorResolution.h"
#include "Resample/Element/DiffuseElement.h"

namespace {

inline size_t xcoord(size_t index, size_t sizeX, size_t sizeY)
{
    return index / sizeY % sizeX;
}

inline size_t ycoord(size_t index, size_t sizeY)
{
    return index % sizeY;
}

} // namespace


IDetector::IDetector() = default;

IDetector::IDetector(const IDetector& other)
    : INode()
    , m_explicitROI(other.m_explicitROI)
    , m_axes(other.m_axes)
    , m_polAnalyzer(other.m_polAnalyzer)
    , m_mask(other.m_mask)
{
    if (other.m_resolution)
        setDetectorResolution(*other.m_resolution);
}

IDetector::~IDetector() = default;

void IDetector::addDetAxis(const std::shared_ptr<IAxis>& axis)
{
    m_axes.push_back(axis);
    if (rank() == 2)
        m_mask.reset(new DetectorMask(*m_axes[0], *m_axes[1]));
}

size_t IDetector::rank() const
{
    return m_axes.size();
}

void IDetector::clear()
{
    m_axes.clear();
}

const IAxis& IDetector::axis(size_t index) const
{
    ASSERT(index < rank());
    return *m_axes[index];
}

size_t IDetector::axisBinIndex(size_t index, size_t selected_axis) const
{
    const size_t dim = rank();
    size_t remainder(index);
    size_t i_axis = dim;
    for (size_t i = 0; i < dim; ++i) {
        --i_axis;
        if (selected_axis == i_axis)
            return remainder % m_axes[i_axis]->size();
        remainder /= m_axes[i_axis]->size();
    }
    ASSERT(0);
}

size_t IDetector::sizeOfExplicitRegionOfInterest() const
{
    if (m_explicitROI.size() != m_axes.size())
        return 0;

    size_t s = 1;
    for (const auto& roiOfAxis : m_explicitROI)
        s *= roiOfAxis.roiSize;

    return s;
}

std::pair<double, double> IDetector::boundsOfExplicitRegionOfInterest(size_t iAxis) const
{
    ASSERT(iAxis < rank());
    if (iAxis < m_explicitROI.size())
        return {m_explicitROI[iAxis].lower, m_explicitROI[iAxis].upper};
    return {0., 0.};
}

size_t IDetector::totalSize() const
{
    const size_t dim = rank();
    if (dim == 0)
        return 0;
    size_t result = 1;
    for (size_t i_axis = 0; i_axis < dim; ++i_axis)
        result *= m_axes[i_axis]->size();
    return result;
}

size_t IDetector::sizeOfRegionOfInterest() const
{
    const auto explicitSize = sizeOfExplicitRegionOfInterest();
    return (explicitSize != 0) ? explicitSize : totalSize();
}

bool IDetector::hasExplicitRegionOfInterest() const
{
    return !m_axes.empty() && (m_explicitROI.size() == m_axes.size());
}

std::vector<const IAxis*> IDetector::axesClippedToRegionOfInterest() const
{
    std::vector<const IAxis*> result;
    for (size_t iAxis = 0; iAxis < m_axes.size(); ++iAxis) {
        auto* axis = m_axes[iAxis]->clone();
        axis->clip(regionOfInterestBounds(iAxis));
        result.emplace_back(axis);
    }
    return result;
}

void IDetector::setAnalyzer(const R3 direction, double efficiency, double total_transmission)
{
    m_polAnalyzer = PolFilter(direction, efficiency, total_transmission);
}

void IDetector::setDetectorResolution(const IDetectorResolution& detector_resolution)
{
    m_resolution.reset(detector_resolution.clone());
}

// TODO: pass dimension-independent argument to this function
void IDetector::setResolutionFunction(const IResolutionFunction2D& resFunc)
{
    ConvolutionDetectorResolution convFunc(resFunc);
    setDetectorResolution(convFunc);
}

void IDetector::resetRegionOfInterest()
{
    m_explicitROI.clear();
}

void IDetector::applyDetectorResolution(Datafield* intensity_map) const
{
    ASSERT(intensity_map);

    if (m_resolution) {
        m_resolution->applyDetectorResolution(intensity_map);
        if (detectorMask() && detectorMask()->hasMasks()) {
            // sets amplitude in masked areas to zero
            std::unique_ptr<Datafield> buff(new Datafield(intensity_map->frame().cloned_axes()));
            iterateOverNonMaskedPoints([&](const_iterator it) {
                (*buff)[it.roiIndex()] = (*intensity_map)[it.roiIndex()];
            });
            intensity_map->setVector(buff->flatVector());
        }
    }
}

const Datafield*
IDetector::createDetectorIntensity(const std::vector<DiffuseElement>& elements) const
{
    std::unique_ptr<Datafield> detectorMap(createDetectorMap());
    ASSERT(detectorMap);
    size_t elementIndex = 0;
    for (auto it = beginNonMaskedPoints(); it != endNonMaskedPoints(); ++it)
        (*detectorMap)[it.roiIndex()] = elements[elementIndex++].intensity();
    if (m_resolution)
        applyDetectorResolution(detectorMap.get());

    return detectorMap.release();
}

std::unique_ptr<Datafield> IDetector::createDetectorMap() const
{
    const size_t dim = rank();
    ASSERT(dim != 0);

    std::vector<const IAxis*> axes;
    for (size_t iAxis = 0; iAxis < dim; ++iAxis) {
        IAxis* tmp = axis(iAxis).clone();
        tmp->clip(regionOfInterestBounds(iAxis));
        axes.emplace_back(tmp);
    }
    return std::make_unique<Datafield>(axes);
}

std::pair<double, double> IDetector::regionOfInterestBounds(size_t iAxis) const
{
    ASSERT(iAxis < m_axes.size());

    const auto explicitBounds = boundsOfExplicitRegionOfInterest(iAxis);
    if (explicitBounds.first != 0 || explicitBounds.second != 0)
        return explicitBounds;

    return m_axes[iAxis]->bounds();
}

std::vector<const INode*> IDetector::nodeChildren() const
{
    return std::vector<const INode*>() << &m_polAnalyzer << m_resolution.get();
}

void IDetector::iterateOverRegionOfInterest(std::function<void(const_iterator)> func) const
{
    if (this->rank() == 0)
        return;

    for (auto it = beginRegionOfInterestPoints(); it != endRegionOfInterestPoints(); ++it)
        func(it);
}

void IDetector::iterateOverNonMaskedPoints(std::function<void(const_iterator)> func) const
{
    if (this->rank() == 0)
        return;

    for (auto it = beginNonMaskedPoints(); it != endNonMaskedPoints(); ++it)
        func(it);
}

SimulationAreaIterator IDetector::beginNonMaskedPoints() const
{
    return SimulationAreaIterator::createBegin(this, SimulationAreaIterator::notMasked);
}

SimulationAreaIterator IDetector::endNonMaskedPoints() const
{
    return SimulationAreaIterator::createEnd(this, SimulationAreaIterator::notMasked);
}

SimulationAreaIterator IDetector::beginRegionOfInterestPoints() const
{
    return SimulationAreaIterator::createBegin(this, SimulationAreaIterator::regionOfInterest);
}

SimulationAreaIterator IDetector::endRegionOfInterestPoints() const
{
    return SimulationAreaIterator::createEnd(this, SimulationAreaIterator::regionOfInterest);
}

size_t IDetector::regionOfInterestIndexToDetectorIndex(const size_t regionOfInterestIndex) const
{
    if (m_explicitROI.size() != 2)
        return regionOfInterestIndex;

    const auto& x = m_explicitROI[0];
    const auto& y = m_explicitROI[1];

    const size_t globalIndex0 = y.lowerIndex + x.lowerIndex * y.detectorSize;
    return globalIndex0 + ycoord(regionOfInterestIndex, y.roiSize)
           + xcoord(regionOfInterestIndex, x.roiSize, y.roiSize) * y.detectorSize;
}

size_t IDetector::detectorIndexToRegionOfInterestIndex(const size_t detectorIndex) const
{
    if (m_explicitROI.size() != 2)
        return detectorIndex;

    const auto& x = m_explicitROI[0];
    const auto& y = m_explicitROI[1];

    const size_t ny = ycoord(detectorIndex, y.detectorSize);
    if (ny < y.lowerIndex || ny > y.upperIndex)
        throw std::runtime_error("IDetector::detectorIndexToRegionOfInterestIndex() -> Error.");

    const size_t nx = xcoord(detectorIndex, x.detectorSize, y.detectorSize);
    if (nx < x.lowerIndex || nx > x.upperIndex)
        throw std::runtime_error("IDetector::detectorIndexToRegionOfInterestIndex() -> Error.");

    return ny - y.lowerIndex + (nx - x.lowerIndex) * y.roiSize;
}

/* -- IDetector::RoiOfAxis ----------------------------------------------------------- */

IDetector::RoiOfAxis::RoiOfAxis(const IAxis& axis, double _lower, double _upper)
    : lower(_lower)
    , upper(_upper)
    , lowerIndex(axis.findClosestIndex(lower))
    , upperIndex(axis.findClosestIndex(upper))
    , roiSize(upperIndex - lowerIndex + 1)
    , detectorSize(axis.size())
{
}

/* -- from IDetector2D -- */

void IDetector::setDetectorParameters(size_t n_x, double x_min, double x_max, size_t n_y,
                                      double y_min, double y_max)
{
    clear();
    addDetAxis(std::make_shared<FixedBinAxis>(axisName(0), n_x, x_min, x_max));
    addDetAxis(std::make_shared<FixedBinAxis>(axisName(1), n_y, y_min, y_max));
}

void IDetector::setRegionOfInterest(double xlow, double ylow, double xup, double yup)
{
    ASSERT(rank() == 2);

    m_explicitROI.clear();
    m_explicitROI.emplace_back(axis(0), xlow, xup);
    m_explicitROI.emplace_back(axis(1), ylow, yup);
}

std::vector<size_t> IDetector::active_indices() const
{
    std::vector<size_t> result;

    iterateOverNonMaskedPoints([&](const_iterator it) { result.push_back(it.detectorIndex()); });

    return result;
}

void IDetector::addMask(const IShape2D& shape, bool mask_value)
{
    m_mask->addMask(shape, mask_value);
}

void IDetector::maskAll()
{
    if (rank() != 2)
        return;
    addMask(InfinitePlane(), true);
}

const DetectorMask* IDetector::detectorMask() const
{
    return m_mask.get();
}

size_t IDetector::getGlobalIndex(size_t x, size_t y) const
{
    ASSERT(rank() == 2);
    return x * axis(1).size() + y;
}