From 55c20c406637f93ba08ddd20bde31ea2b0c13302 Mon Sep 17 00:00:00 2001
From: Gennady Pospelov <g.pospelov@fz-juelich.de>
Date: Tue, 18 Oct 2016 10:30:24 +0200
Subject: [PATCH] IDetector::getIntensityData now returns ROI-clipped map in
desired units.
---
Core/Instrument/IDetector2D.cpp | 69 ++++++++++++++++---
Core/Instrument/IDetector2D.h | 5 ++
.../UnitTests/Core/3/SphericalDetectorTest.h | 69 +++++++++++++++++++
3 files changed, 135 insertions(+), 8 deletions(-)
diff --git a/Core/Instrument/IDetector2D.cpp b/Core/Instrument/IDetector2D.cpp
index 6136d0836da..ebc3aeda49f 100644
--- a/Core/Instrument/IDetector2D.cpp
+++ b/Core/Instrument/IDetector2D.cpp
@@ -117,21 +117,20 @@ OutputData<double> *IDetector2D::getDetectorIntensity(const OutputData<double> &
const Beam& beam,
IDetector2D::EAxesUnits units_type) const
{
- std::unique_ptr<OutputData<double> > result (data.clone());
+ std::unique_ptr<OutputData<double>> result(data.clone());
applyDetectorResolution(result.get());
- if(units_type == IDetector2D::DEFAULT) {
+ if(units_type == IDetector2D::DEFAULT && regionOfInterest() == nullptr)
return result.release();
- }
- OutputData<double>* detectorMap = createDetectorMap(beam, units_type);
+ std::unique_ptr<OutputData<double>> detectorMap(createDetectorMap(beam, units_type));
if(!detectorMap)
throw Exceptions::RuntimeErrorException("Instrument::getDetectorIntensity() -> Error."
"Can't create detector map.");
- detectorMap->setRawDataVector(result->getRawDataVector());
- return detectorMap;
+ setDataToDetectorMap(*detectorMap.get(), *result.get());
+ return detectorMap.release();
}
OutputData<double>* IDetector2D::createDetectorMap(const Beam& beam, EAxesUnits units) const
@@ -296,8 +295,35 @@ std::unique_ptr<IAxis> IDetector2D::constructAxis(size_t axis_index, const Beam
{
double amin(0), amax(0);
calculateAxisRange(axis_index, beam, units, amin, amax);
- return std::unique_ptr<IAxis>(
- new FixedBinAxis(getAxisName(axis_index), getAxis(axis_index).getSize(), amin, amax));
+
+ std::unique_ptr<IAxis> result(new FixedBinAxis(getAxisName(axis_index),
+ getAxis(axis_index).getSize(), amin, amax));
+
+ if(m_region_of_interest) {
+ return clipAxisToRoi(axis_index, *result.get());
+ } else {
+ return result;
+ }
+}
+
+std::unique_ptr<IAxis> IDetector2D::clipAxisToRoi(size_t axis_index, const IAxis &axis) const
+{
+ if(!m_region_of_interest)
+ throw Exceptions::RuntimeErrorException("IDetector2D::clipAxisToRoi() -> Error. ROI is "
+ "absent");
+
+ size_t nbin1(0), nbin2(axis.getSize()-1);
+ if(axis_index == BornAgain::X_AXIS_INDEX) {
+ nbin1 = getAxis(axis_index).findClosestIndex(m_region_of_interest->getXlow());
+ nbin2 = getAxis(axis_index).findClosestIndex(m_region_of_interest->getXup());
+
+ }else if(axis_index == BornAgain::Y_AXIS_INDEX) {
+ nbin1 = getAxis(axis_index).findClosestIndex(m_region_of_interest->getYlow());
+ nbin2 = getAxis(axis_index).findClosestIndex(m_region_of_interest->getYup());
+ }
+
+ return std::unique_ptr<IAxis>(new FixedBinAxis(axis.getName(), nbin2-nbin1+1,
+ axis.getBin(nbin1).m_lower, axis.getBin(nbin2).m_upper));
}
void IDetector2D::calculateAxisRange(size_t axis_index, const Beam &beam,
@@ -395,6 +421,33 @@ Eigen::Matrix2cd IDetector2D::calculateAnalyzerOperator(
return result;
}
+void IDetector2D::setDataToDetectorMap(OutputData<double> &detectorMap,
+ const OutputData<double> &data) const
+{
+ if(detectorMap.getAllocatedSize() == data.getAllocatedSize()) {
+ detectorMap.setRawDataVector(data.getRawDataVector());
+ return;
+ }
+
+ if(const Geometry::Rectangle *roi = regionOfInterest()) {
+ size_t roi_x = getAxis(BornAgain::X_AXIS_INDEX).findClosestIndex(roi->getXlow());
+ size_t roi_y = getAxis(BornAgain::Y_AXIS_INDEX).findClosestIndex(roi->getYlow());
+ size_t index0 = getGlobalIndex(roi_x, roi_y);
+ const IAxis &yAxisOfDetector = getAxis(BornAgain::Y_AXIS_INDEX);
+
+ const IAxis &xAxisOfMap = *detectorMap.getAxis(BornAgain::X_AXIS_INDEX);
+ const IAxis &yAxisOfMap = *detectorMap.getAxis(BornAgain::Y_AXIS_INDEX);
+ for(size_t ix=0; ix<xAxisOfMap.getSize(); ++ix) {
+ for(size_t iy=0; iy<yAxisOfMap.getSize(); ++iy) {
+ size_t mapIndex = iy + ix*yAxisOfMap.getSize();
+ size_t globalIndex = index0 + iy + ix*yAxisOfDetector.getSize();
+ detectorMap[mapIndex] = data[globalIndex];
+ }
+ }
+ }
+
+}
+
void IDetector2D::addAxis(const IAxis& axis)
{
m_axes.push_back(axis.clone());
diff --git a/Core/Instrument/IDetector2D.h b/Core/Instrument/IDetector2D.h
index 3481c64a3ff..7a99d61f6e7 100644
--- a/Core/Instrument/IDetector2D.h
+++ b/Core/Instrument/IDetector2D.h
@@ -174,6 +174,8 @@ protected:
std::unique_ptr<IAxis> constructAxis(size_t axis_index, const Beam& beam,
EAxesUnits units) const;
+ std::unique_ptr<IAxis> clipAxisToRoi(size_t axis_index, const IAxis &axis) const;
+
//! Calculates axis range from original detector axes in given units (mm, rad, etc)
virtual void calculateAxisRange(size_t axis_index, const Beam& beam, EAxesUnits units,
double &amin, double &amax) const;
@@ -216,6 +218,9 @@ private:
const kvector_t direction, double efficiency, double total_transmission = 1.0) const;
#endif
+ void setDataToDetectorMap(OutputData<double> &detectorMap,
+ const OutputData<double> &data) const;
+
std::unique_ptr<Geometry::Rectangle> m_region_of_interest;
};
diff --git a/Tests/UnitTests/Core/3/SphericalDetectorTest.h b/Tests/UnitTests/Core/3/SphericalDetectorTest.h
index 2a827a03685..52dab56ca67 100644
--- a/Tests/UnitTests/Core/3/SphericalDetectorTest.h
+++ b/Tests/UnitTests/Core/3/SphericalDetectorTest.h
@@ -224,6 +224,75 @@ TEST_F(SphericalDetectorTest, regionOfInterestAndData)
EXPECT_EQ(data.getAxis(1)->getMax(), 4.0);
}
+//! Create detector map in the presence of region of interest.
+
+TEST_F(SphericalDetectorTest, regionOfInterestAndDetectorMap)
+{
+ SphericalDetector detector(6, -1.0*Units::deg, 5.0*Units::deg,
+ 4, 0.0*Units::deg, 4.0*Units::deg);
+
+ detector.setRegionOfInterest(0.1*Units::deg, 1.1*Units::deg, 3.0*Units::deg, 2.9*Units::deg);
+ Beam beam;
+ beam.setCentralK(1.0*Units::angstrom, 0.4*Units::deg, 0.0);
+
+ // Creating map in default units, which are radians and checking that axes are clipped
+ // to region of interest.
+ std::unique_ptr<OutputData<double>> data(
+ detector.createDetectorMap(beam, IDetector2D::DEFAULT));
+ EXPECT_EQ(data->getAxis(0)->getSize(), 4);
+ EXPECT_EQ(data->getAxis(0)->getMin(), 0.0*Units::deg);
+ EXPECT_EQ(data->getAxis(0)->getMax(), 4.0*Units::deg);
+ EXPECT_EQ(data->getAxis(1)->getSize(), 2);
+ EXPECT_EQ(data->getAxis(1)->getMin(), 1.0*Units::deg);
+ EXPECT_EQ(data->getAxis(1)->getMax(), 3.0*Units::deg);
+
+ // Creating map with axes in degrees, and checking that it is clipped to the region of interest
+ data.reset(detector.createDetectorMap(beam, IDetector2D::DEGREES));
+ EXPECT_EQ(data->getAxis(0)->getSize(), 4);
+ EXPECT_EQ(data->getAxis(0)->getMin(), 0.0);
+ EXPECT_EQ(data->getAxis(0)->getMax(), 4.0);
+ EXPECT_EQ(data->getAxis(1)->getSize(), 2);
+ EXPECT_EQ(data->getAxis(1)->getMin(), 1.0);
+ EXPECT_EQ(data->getAxis(1)->getMax(), 3.0);
+}
+
+//! Checking IDetector2D::getIntensityData in the presence of region of interest.
+
+TEST_F(SphericalDetectorTest, getIntensityData)
+{
+ SphericalDetector detector(6, -1.0*Units::deg, 5.0*Units::deg,
+ 4, 0.0*Units::deg, 4.0*Units::deg);
+ detector.setRegionOfInterest(0.1*Units::deg, 1.1*Units::deg, 3.0*Units::deg, 2.9*Units::deg);
+ Beam beam;
+ beam.setCentralK(1.0*Units::angstrom, 0.4*Units::deg, 0.0);
+
+ // Initializing data (no region of interest involved yet) and filling with amplitudes
+ OutputData<double> intensityData;
+ detector.initOutputData(intensityData);
+ EXPECT_EQ(intensityData.getAllocatedSize(), 6*4);
+ for(size_t i=0; i<intensityData.getAllocatedSize(); ++i) {
+ intensityData[i] = static_cast<double>(i);
+ }
+ EXPECT_EQ(intensityData[intensityData.getAllocatedSize()-1], 23.0);
+
+ // Getting detectorIntensity and checking that amplitudes are correct and it is clipped to
+ // region of interest.
+
+ std::unique_ptr<OutputData<double>> detectorIntensity(
+ detector.getDetectorIntensity(intensityData, beam, IDetector2D::DEGREES));
+
+ EXPECT_EQ(detectorIntensity->getAllocatedSize(), 8);
+ EXPECT_EQ((*detectorIntensity)[0], 5.0);
+ EXPECT_EQ((*detectorIntensity)[7], 18.0);
+ EXPECT_EQ(detectorIntensity->getAxis(0)->getSize(), 4);
+ EXPECT_EQ(detectorIntensity->getAxis(0)->getMin(), 0.0);
+ EXPECT_EQ(detectorIntensity->getAxis(0)->getMax(), 4.0);
+ EXPECT_EQ(detectorIntensity->getAxis(1)->getSize(), 2);
+ EXPECT_EQ(detectorIntensity->getAxis(1)->getMin(), 1.0);
+ EXPECT_EQ(detectorIntensity->getAxis(1)->getMax(), 3.0);
+}
+
+
TEST_F(SphericalDetectorTest, DetectorCopying)
{
--
GitLab