SpecularData container in SimulationElement

Redmine: #1910
Temporary solution to transfer specular data

Core/Binning/SimulationElement.cpp

@@ -25,25 +25,28 @@ SimulationElement::SimulationElement(double wavelength, double alpha_i, double p
     , m_phi_i(phi_i)
     , m_intensity(0.0)
     , mP_pixel(std::move(pixel))
-    , m_contains_specular(false)
 {
     initPolarization();
 }
 
 SimulationElement::SimulationElement(const SimulationElement& other)
     : m_wavelength(other.m_wavelength), m_alpha_i(other.m_alpha_i), m_phi_i(other.m_phi_i)
-    , m_intensity(other.m_intensity), m_contains_specular(other.m_contains_specular)
+    , m_intensity(other.m_intensity)
 {
     mP_pixel.reset(other.mP_pixel->clone());
+    if (other.m_specular_data)
+        m_specular_data.reset(new SpecularData(*other.m_specular_data));
     m_polarization = other.m_polarization;
     m_analyzer_operator = other.m_analyzer_operator;
 }
 
 SimulationElement::SimulationElement(const SimulationElement& other, double x, double y)
     : m_wavelength(other.m_wavelength), m_alpha_i(other.m_alpha_i), m_phi_i(other.m_phi_i)
-    , m_intensity(other.m_intensity), m_contains_specular(other.m_contains_specular)
+    , m_intensity(other.m_intensity)
 {
     mP_pixel.reset(other.mP_pixel->createZeroSizePixel(x, y));
+    if (other.m_specular_data)
+        m_specular_data.reset(new SpecularData(*other.m_specular_data));
     m_polarization = other.m_polarization;
     m_analyzer_operator = other.m_analyzer_operator;
 }
@@ -56,7 +59,7 @@ SimulationElement::SimulationElement(SimulationElement&& other) noexcept
     , m_polarization(std::move(other.m_polarization))
     , m_analyzer_operator(std::move(other.m_analyzer_operator))
     , mP_pixel(std::move(other.mP_pixel))
-    , m_contains_specular(other.m_contains_specular)
+    , m_specular_data(std::move(other.m_specular_data))
 {
 }
 
@@ -108,7 +111,7 @@ void SimulationElement::swapContent(SimulationElement &other)
     std::swap(m_polarization, other.m_polarization);
     std::swap(m_analyzer_operator, other.m_analyzer_operator);
     std::swap(mP_pixel, other.mP_pixel);
-    std::swap(m_contains_specular, other.m_contains_specular);
+    std::swap(m_specular_data, other.m_specular_data);
 }
 
 void SimulationElement::initPolarization()
@@ -127,14 +130,14 @@ double SimulationElement::getPhi(double x, double y) const
     return getKf(x,y).phi();
 }
 
-bool SimulationElement::containsSpecularWavevector() const
+void SimulationElement::setSpecular()
 {
-    return m_contains_specular;
+    m_specular_data.reset(new SpecularData);
 }
 
-void SimulationElement::setSpecular(bool contains_specular)
+void SimulationElement::setSpecular(std::unique_ptr<SpecularData> specular_data)
 {
-    m_contains_specular = contains_specular;
+    m_specular_data = std::move(specular_data);
 }
 
 double SimulationElement::getIntegrationFactor(double x, double y) const {
@@ -152,3 +155,26 @@ void addElementsWithWeight(std::vector<SimulationElement>::const_iterator first,
     for (std::vector<SimulationElement>::const_iterator it = first; it != last; ++it, ++result)
         result->addIntensity(it->getIntensity() * weight);
 }
+
+SpecularData::SpecularData() : data_type_used(DATA_TYPE::Invalid) {}
+
+SpecularData::SpecularData(MatrixVector coefficients)
+    : data(std::move(coefficients))
+    , data_type_used(DATA_TYPE::Matrix)
+{}
+
+SpecularData::SpecularData(ScalarVector coefficients)
+    : data(std::move(coefficients))
+    , data_type_used(DATA_TYPE::Scalar)
+{}
+
+const ILayerRTCoefficients& SpecularData::operator[](size_t index) const
+{
+    if (data_type_used == DATA_TYPE::Invalid)
+        throw std::runtime_error(
+            "Error in SpecularData::operator[]: attempt to access uninitialized data");
+    if (data_type_used == DATA_TYPE::Scalar)
+        return (*boost::get<ScalarVector>(&data))[index];
+    else
+        return (*boost::get<MatrixVector>(&data))[index];
+}

Core/Binning/SimulationElement.h

@@ -20,10 +20,14 @@
 #include "EigenCore.h"
 #include "Vectors3D.h"
 #include "IPixel.h"
+#include "MatrixRTCoefficients.h"
+#include "ScalarRTCoefficients.h"
+#include <boost/variant.hpp>
 #include <memory>
 #include <vector>
 
 class IPixel;
+class SpecularData;
 
 //! Data stucture containing both input and output of a single detector cell.
 //! @ingroup simulation
@@ -80,11 +84,12 @@ public:
     double getAlpha(double x, double y) const;
     double getPhi(double x, double y) const;
 
-    //! check if element contains given wavevector
-    bool containsSpecularWavevector() const;
+    //! check if element corresponds to specular peak
+    SpecularData* specularData() const {return m_specular_data.get();}
 
-    //! indicate that this element contains the specular wavevector
-    void setSpecular(bool contains_specular);
+    //! Turn on specular data
+    void setSpecular();
+    void setSpecular(std::unique_ptr<SpecularData> specular_data);
 
 private:
     void swapContent(SimulationElement &other);
@@ -99,7 +104,7 @@ private:
     Eigen::Matrix2cd m_analyzer_operator; //!< polarization analyzer operator
 #endif
     std::unique_ptr<IPixel> mP_pixel;
-    bool m_contains_specular;
+    std::unique_ptr<SpecularData> m_specular_data;
 };
 
 
@@ -109,4 +114,29 @@ void addElementsWithWeight(std::vector<SimulationElement>::const_iterator first,
                            std::vector<SimulationElement>::iterator result,
                            double weight);
 
+//! Helper class for SimulationElement to carry specular information
+//! @ingroup simulation
+
+class SpecularData
+{
+    // FIXME: find a better way to carry the specular data in SimulationElement
+    using ScalarVector = std::vector<ScalarRTCoefficients>;
+    using MatrixVector = std::vector<MatrixRTCoefficients>;
+public:
+    SpecularData();
+
+    SpecularData(MatrixVector coefficients);
+
+    SpecularData(ScalarVector coefficients);
+
+    const ILayerRTCoefficients& operator[](size_t index) const;
+
+    bool isInited() const {return data_type_used != DATA_TYPE::Invalid;}
+
+private:
+    enum class DATA_TYPE { Invalid = -1, Scalar, Matrix };
+    boost::variant<ScalarVector, MatrixVector> data;
+    DATA_TYPE data_type_used;
+};
+
 #endif // SIMULATIONELEMENT_H

Core/Computation/SpecularComputationTerm.cpp

@@ -17,12 +17,13 @@ void SpecularComputationTerm::eval(ProgressHandler*,
         return;
 
     for (auto it = begin_it; it != end_it; ++it)
-        if (it->containsSpecularWavevector())
+        if (it->specularData())
             evalSingle(it);
 }
 
 void SpecularComputationTerm::evalSingle(const std::vector<SimulationElement>::iterator& iter) const
 {
-    complex_t R = mp_fresnel_map->getInCoefficients(*iter, 0)->getScalarR();
-    iter->setIntensity(std::norm(R));
+    mp_fresnel_map->fillSpecularData(*iter);
+    const ILayerRTCoefficients& layer_data = (*iter->specularData())[0];
+    iter->setIntensity(std::norm(layer_data.getScalarR()));
 }

Core/Instrument/IDetector2D.cpp

@@ -111,7 +111,7 @@ std::vector<SimulationElement> IDetector2D::createSimulationElements(const Beam
         sim_element.setPolarization(beam_polarization);
         sim_element.setAnalyzerOperator(analyzer_operator);
         if (it.index()==spec_index) {
-            sim_element.setSpecular(true);
+            sim_element.setSpecular();
         }
         result.push_back(std::move(sim_element));
     }

Core/Instrument/SpecularDetector1D.cpp

@@ -80,7 +80,7 @@ std::vector<SimulationElement> SpecularDetector1D::createSimulationElements(cons
                                       std::make_unique<SpecularPixel>(alpha_i));
         sim_element.setPolarization(beam_polarization);
         sim_element.setAnalyzerOperator(analyzer_operator);
-        sim_element.setSpecular(true);
+        sim_element.setSpecular();
         result.push_back(std::move(sim_element));
     }
 

Core/Multilayer/IFresnelMap.h

@@ -23,6 +23,7 @@
 class ILayerRTCoefficients;
 class MultiLayer;
 class SimulationElement;
+struct SpecularData;
 
 //! Holds the necessary information to calculate the radiation wavefunction in every layer
 //! for different incoming (outgoing) angles of the beam in the top layer
@@ -43,6 +44,9 @@ public:
     virtual const ILayerRTCoefficients* getInCoefficients(
             const SimulationElement& sim_element, size_t layer_index) const =0;
 
+    //! Fills simulation element specular data
+    virtual void fillSpecularData(SimulationElement& sim_element) const = 0;
+
     //! Sets the multilayer to be used for the Fresnel calculations.
     virtual void setMultilayer(const MultiLayer& multilayer);
 

Core/Multilayer/MatrixFresnelMap.cpp

@@ -46,6 +46,17 @@ MatrixFresnelMap::getInCoefficients(const SimulationElement& sim_element, size_t
     return getCoefficients(sim_element.getKi(), layer_index, *mP_multilayer, m_hash_table_in);
 }
 
+void MatrixFresnelMap::fillSpecularData(SimulationElement& sim_element) const
+{
+    const auto& kvec = sim_element.getKi();
+    std::vector<MatrixRTCoefficients> coef_vector;
+    if (m_use_cache)
+        coef_vector = getCoefficientsFromCache(kvec, *mP_multilayer, m_hash_table_in);
+    else
+        coef_vector = calculateCoefficients(*mP_multilayer, kvec);
+    sim_element.setSpecular(std::make_unique<SpecularData>(std::move(coef_vector)));
+}
+
 const ILayerRTCoefficients* MatrixFresnelMap::getCoefficients(kvector_t kvec, size_t layer_index,
                                                 const MultiLayer& multilayer,
                                                 CoefficientHash& hash_table) const

Core/Multilayer/MatrixFresnelMap.h

@@ -44,6 +44,9 @@ public:
 
     void setMultilayer(const MultiLayer& multilayer) final override;
 
+    //! Fills simulation element specular data
+    virtual void fillSpecularData(SimulationElement& sim_element) const override;
+
     typedef std::unordered_map<kvector_t, std::vector<MatrixRTCoefficients>, HashKVector>
         CoefficientHash;
 

Core/Multilayer/ScalarFresnelMap.cpp

@@ -42,6 +42,17 @@ const ILayerRTCoefficients* ScalarFresnelMap::getInCoefficients(
     return getCoefficients(sim_element.getKi(), layer_index);
 }
 
+void ScalarFresnelMap::fillSpecularData(SimulationElement& sim_element) const
+{
+    const auto& kvec = sim_element.getKi();
+    std::vector<ScalarRTCoefficients> coef_vector;
+    if (m_use_cache)
+        coef_vector = getCoefficientsFromCache(kvec);
+    else
+        coef_vector = calculateCoefficients(*mP_multilayer, kvec);
+    sim_element.setSpecular(std::make_unique<SpecularData>(std::move(coef_vector)));
+}
+
 const ScalarRTCoefficients* ScalarFresnelMap::getCoefficients(
         kvector_t kvec, size_t layer_index) const
 {

Core/Multilayer/ScalarFresnelMap.h

@@ -44,6 +44,9 @@ public:
     virtual const ILayerRTCoefficients* getInCoefficients(
         const SimulationElement& sim_element, size_t layer_index) const final override;
 
+    //! Fills simulation element specular data
+    virtual void fillSpecularData(SimulationElement& sim_element) const override;
+
 private:
     const ScalarRTCoefficients* getCoefficients(kvector_t kvec, size_t layer_index) const;
     const std::vector<ScalarRTCoefficients>& getCoefficientsFromCache(kvector_t kvec) const;