IDetector.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  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/DetectorContext.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() {}
 
IDetector::IDetector(const IDetector& other)
    : INode()
    , m_explicitROI(other.m_explicitROI)
    , m_axes(other.m_axes)
    , m_polAnalyzer(other.m_polAnalyzer)
    , m_detector_mask(other.m_detector_mask)
{
    if (other.m_detector_resolution)
        setDetectorResolution(*other.m_detector_resolution);
}
 
IDetector::~IDetector() = default;
 
void IDetector::addDetAxis(const IAxis& axis)
{
    m_axes.emplace_back(axis.clone());
    if (rank() == 2)
        m_detector_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());
}
 
OwningVector<IAxis> IDetector::axesClippedToRegionOfInterest() const
{
    OwningVector<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& p_detector_resolution)
{
    m_detector_resolution.reset(p_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* p_intensity_map) const
{
    ASSERT(p_intensity_map);
 
    if (m_detector_resolution) {
        m_detector_resolution->applyDetectorResolution(p_intensity_map);
        if (detectorMask() && detectorMask()->hasMasks()) {
            // sets amplitude in masked areas to zero
            std::unique_ptr<Datafield> buff(new Datafield(p_intensity_map->frame().cloned_axes()));
            iterateOverNonMaskedPoints([&](const_iterator it) {
                (*buff)[it.roiIndex()] = (*p_intensity_map)[it.roiIndex()];
            });
            p_intensity_map->setVector(buff->flatVector());
        }
    }
}
 
const IDetectorResolution* IDetector::detectorResolution() const
{
    return m_detector_resolution.get();
}
 
Datafield* IDetector::createDetectorIntensity(
    const std::vector<std::unique_ptr<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_detector_resolution)
        applyDetectorResolution(detectorMap.get());
 
    return detectorMap.release();
}
 
std::unique_ptr<Datafield> IDetector::createDetectorMap() const
{
    const size_t dim = rank();
    ASSERT(dim != 0);
 
    std::vector<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);
}
 
size_t IDetector::numberOfElements() const
{
    size_t result(0);
    iterateOverNonMaskedPoints([&result](const_iterator) { ++result; });
    return result;
}
 
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_detector_resolution;
}
 
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);
 
    detectorSize = axis.size();
    roiSize = upperIndex - lowerIndex + 1;
}
 
/* -- 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(FixedBinAxis(axisName(0), n_x, x_min, x_max));
    addDetAxis(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;
}
 
std::unique_ptr<DetectorContext> IDetector::createContext() const
{
    return std::make_unique<DetectorContext>(this);
}
 
void IDetector::addMask(const IShape2D& shape, bool mask_value)
{
    m_detector_mask->addMask(shape, mask_value);
}
 
void IDetector::maskAll()
{
    if (rank() != 2)
        return;
    addMask(InfinitePlane(), true);
}
 
const DetectorMask* IDetector::detectorMask() const
{
    return m_detector_mask.get();
}
 
size_t IDetector::getGlobalIndex(size_t x, size_t y) const
{
    if (rank() != 2)
        return totalSize();
    return x * axis(1).size() + y;
}