From 49b8b0758cacffc34f71cf5fb3d4ae4cf4641dfc Mon Sep 17 00:00:00 2001
From: Walter Van Herck <w.van.herck@fz-juelich.de>
Date: Wed, 11 Jul 2012 10:55:27 +0200
Subject: [PATCH] Detector resolution function (step 2): finished
---
App/src/TestDetectorResolution.cpp | 6 +--
.../inc/ConvolutionDetectorResolution.h | 14 +++---
Core/Algorithms/inc/GISASExperiment.h | 2 +-
.../src/ConvolutionDetectorResolution.cpp | 46 +++++++++++++++----
Core/Algorithms/src/GISASExperiment.cpp | 2 +-
Core/Tools/inc/MathFunctions.h | 2 +-
Core/Tools/src/MathFunctions.cpp | 7 ++-
7 files changed, 52 insertions(+), 27 deletions(-)
diff --git a/App/src/TestDetectorResolution.cpp b/App/src/TestDetectorResolution.cpp
index 212c1262521..9b4b03169cb 100644
--- a/App/src/TestDetectorResolution.cpp
+++ b/App/src/TestDetectorResolution.cpp
@@ -29,7 +29,7 @@ void TestDetectorResolution::execute()
GISASExperiment experiment;
experiment.setSample(mp_sample);
experiment.setDetectorParameters(0.0*Units::degree, 2.0*Units::degree, 100
- , 0.0*Units::degree, 2.0*Units::degree, 100, true);
+ , 0.0*Units::degree, 2.0*Units::degree, 100);
experiment.setDetectorResolutionFunction(&testResolutionFunction);
experiment.setBeamParameters(1.0*Units::angstrom, -0.2*Units::degree, 0.0*Units::degree);
experiment.runSimulation();
@@ -42,9 +42,7 @@ double testResolutionFunction(double u, double v)
{
double sigma_u = 0.001;
double sigma_v = 0.001;
- double step_u = 0.02*Units::degree;
- double step_v = 0.02*Units::degree;
- return MathFunctions::IntegratedGaussian(u, 0.0, sigma_u, step_u)*MathFunctions::IntegratedGaussian(v, 0.0, sigma_v, step_v);
+ return MathFunctions::IntegratedGaussian(u, 0.0, sigma_u)*MathFunctions::IntegratedGaussian(v, 0.0, sigma_v);
}
void TestDetectorResolution::initializeSample()
diff --git a/Core/Algorithms/inc/ConvolutionDetectorResolution.h b/Core/Algorithms/inc/ConvolutionDetectorResolution.h
index c434df3d01d..4f8ae180727 100644
--- a/Core/Algorithms/inc/ConvolutionDetectorResolution.h
+++ b/Core/Algorithms/inc/ConvolutionDetectorResolution.h
@@ -26,12 +26,12 @@
class ConvolutionDetectorResolution : public IDetectorResolution
{
public:
- typedef double (*resolution_function_1d)(double);
- typedef double (*resolution_function_2d)(double, double);
+ typedef double (*cumulative_DF_1d)(double);
+ typedef double (*cumulative_DF_2d)(double, double);
//! Constructor taking a 1 dimensional resolution function as argument
- ConvolutionDetectorResolution(resolution_function_1d res_function_1d);
+ ConvolutionDetectorResolution(cumulative_DF_1d res_function_1d);
//! Constructor taking a 2 dimensional resolution function as argument
- ConvolutionDetectorResolution(resolution_function_2d res_function_2d);
+ ConvolutionDetectorResolution(cumulative_DF_2d res_function_2d);
//! Destructor
virtual ~ConvolutionDetectorResolution();
@@ -39,10 +39,12 @@ public:
virtual void applyDetectorResolution(OutputData<double> *p_intensity_map) const;
private:
size_t m_dimension;
- resolution_function_1d m_res_function_1d;
- resolution_function_2d m_res_function_2d;
+ cumulative_DF_1d m_res_function_1d;
+ cumulative_DF_2d m_res_function_2d;
void apply1dConvolution(const std::vector<NamedVectorBase *> &axes, OutputData<double> *p_intensity_map) const;
void apply2dConvolution(const std::vector<NamedVectorBase *> &axes, OutputData<double> *p_intensity_map) const;
+ double getIntegratedPDF1d(double x, double step) const;
+ double getIntegratedPDF2d(double x, double step_x, double y, double step_y) const;
};
diff --git a/Core/Algorithms/inc/GISASExperiment.h b/Core/Algorithms/inc/GISASExperiment.h
index 5a1e477f62d..92c92a4961a 100644
--- a/Core/Algorithms/inc/GISASExperiment.h
+++ b/Core/Algorithms/inc/GISASExperiment.h
@@ -28,7 +28,7 @@ public:
void setDetectorParameters(double phi_f_min, double phi_f_max, size_t n_phi,
double alpha_f_min, double alpha_f_max, size_t n_alpha, bool isgisaxs_style=false);
- void setDetectorResolutionFunction(ConvolutionDetectorResolution::resolution_function_2d resolution_function);
+ void setDetectorResolutionFunction(ConvolutionDetectorResolution::cumulative_DF_2d resolution_function);
private:
void initializeAnglesIsgisaxs(NamedVector<double> *p_axis, double start, double end, size_t size);
};
diff --git a/Core/Algorithms/src/ConvolutionDetectorResolution.cpp b/Core/Algorithms/src/ConvolutionDetectorResolution.cpp
index ef2c56b9bce..66accbf7001 100644
--- a/Core/Algorithms/src/ConvolutionDetectorResolution.cpp
+++ b/Core/Algorithms/src/ConvolutionDetectorResolution.cpp
@@ -6,7 +6,7 @@
#include <iostream>
ConvolutionDetectorResolution::ConvolutionDetectorResolution(
- resolution_function_1d res_function_1d)
+ cumulative_DF_1d res_function_1d)
: m_dimension(1)
, m_res_function_1d(res_function_1d)
, m_res_function_2d(0)
@@ -14,7 +14,7 @@ ConvolutionDetectorResolution::ConvolutionDetectorResolution(
}
ConvolutionDetectorResolution::ConvolutionDetectorResolution(
- resolution_function_2d res_function_2d)
+ cumulative_DF_2d res_function_2d)
: m_dimension(2)
, m_res_function_1d(0)
, m_res_function_2d(res_function_2d)
@@ -58,15 +58,18 @@ void ConvolutionDetectorResolution::apply1dConvolution(
// Construct source vector from original intensity map
std::vector<double> source_vector = p_intensity_map->getRawDataVector();
size_t data_size = source_vector.size();
+ if (data_size < 2) {
+ return; // No convolution for sets of zero or one element
+ }
// Construct kernel vector from resolution function
if (p_axis->getSize() != data_size) {
throw LogicErrorException("Size of axis for intensity map does not correspond to the size of the data in the map.");
}
- double step_size = std::abs((*p_axis)[0]-(*p_axis)[p_axis->getSize()-1])/p_axis->getSize();
+ double step_size = std::abs((*p_axis)[0]-(*p_axis)[p_axis->getSize()-1])/(data_size-1);
double mid_value = (*p_axis)[p_axis->getSize()/2]; // Needed because Convolve class expects zero at midpoint
std::vector<double> kernel;
for (size_t index=0; index<data_size; ++index) {
- kernel.push_back(m_res_function_1d((*p_axis)[index] - mid_value)*step_size);
+ kernel.push_back(getIntegratedPDF1d((*p_axis)[index] - mid_value, step_size));
}
// Calculate convolution
std::vector<double> result;
@@ -90,6 +93,9 @@ void ConvolutionDetectorResolution::apply2dConvolution(
NamedVector<double> *p_axis_2 = dynamic_cast<NamedVector<double> *>(axes[1]);
size_t axis_size_1 = p_axis_1->getSize();
size_t axis_size_2 = p_axis_2->getSize();
+ if (axis_size_1 < 2 || axis_size_2 < 2) {
+ return; // No 2d convolution for 1d data
+ }
// Construct source vector array from original intensity map
std::vector<double> raw_source_vector = p_intensity_map->getRawDataVector();
std::vector<std::vector<double> > source;
@@ -102,26 +108,23 @@ void ConvolutionDetectorResolution::apply2dConvolution(
source.push_back(row_vector);
}
// Construct kernel vector from resolution function
- double kernel_total_sum = 0.0; // Temporary sum of the kernel values (should be 1 in total)
std::vector<std::vector<double> > kernel;
kernel.resize(axis_size_1);
double mid_value_1 = (*p_axis_1)[axis_size_1/2]; // Needed because Convolve class expects zero at midpoint
double mid_value_2 = (*p_axis_2)[axis_size_2/2]; // Needed because Convolve class expects zero at midpoint
- double step_size_1 = std::abs((*p_axis_1)[0]-(*p_axis_1)[axis_size_1-1])/axis_size_1;
- double step_size_2 = std::abs((*p_axis_2)[0]-(*p_axis_2)[axis_size_2-1])/axis_size_2;
+ double step_size_1 = std::abs((*p_axis_1)[0]-(*p_axis_1)[axis_size_1-1])/(axis_size_1-1);
+ double step_size_2 = std::abs((*p_axis_2)[0]-(*p_axis_2)[axis_size_2-1])/(axis_size_2-1);
for (size_t index_1=0; index_1<axis_size_1; ++index_1) {
double value_1 = (*p_axis_1)[index_1]-mid_value_1;
std::vector<double> row_vector;
row_vector.resize(axis_size_2);
for (size_t index_2=0; index_2<axis_size_2;++index_2) {
double value_2 = (*p_axis_2)[index_2]-mid_value_2;
- double z_value = m_res_function_2d(value_1, value_2); // *step_size_1*step_size_2;
+ double z_value = getIntegratedPDF2d(value_1, step_size_1, value_2, step_size_2);
row_vector[index_2] = z_value;
- kernel_total_sum += z_value;
}
kernel[index_1] = row_vector;
}
- std::cout << "Total kernel value = " << kernel_total_sum << std::endl;
// Calculate convolution
std::vector<std::vector<double> > result;
MathFunctions::Convolve cv;
@@ -136,3 +139,26 @@ void ConvolutionDetectorResolution::apply2dConvolution(
}
p_intensity_map->setRawDataVector(result_vector);
}
+
+double ConvolutionDetectorResolution::getIntegratedPDF1d(double x,
+ double step) const
+{
+ double halfstep = step/2.0;
+ double xmin = x - halfstep;
+ double xmax = x + halfstep;
+ return m_res_function_1d(xmax) - m_res_function_1d(xmin);
+}
+
+double ConvolutionDetectorResolution::getIntegratedPDF2d(double x,
+ double step_x, double y, double step_y) const
+{
+ double halfstepx = step_x/2.0;
+ double halfstepy = step_y/2.0;
+ double xmin = x - halfstepx;
+ double xmax = x + halfstepx;
+ double ymin = y - halfstepy;
+ double ymax = y + halfstepy;
+ double result = m_res_function_2d(xmax, ymax) - m_res_function_2d(xmax, ymin)
+ - m_res_function_2d(xmin, ymax) + m_res_function_2d(xmin, ymin);
+ return result;
+}
diff --git a/Core/Algorithms/src/GISASExperiment.cpp b/Core/Algorithms/src/GISASExperiment.cpp
index 017626e462d..9f58b7d8b71 100644
--- a/Core/Algorithms/src/GISASExperiment.cpp
+++ b/Core/Algorithms/src/GISASExperiment.cpp
@@ -44,7 +44,7 @@ void GISASExperiment::setDetectorParameters(double phi_f_min, double phi_f_max,
}
void GISASExperiment::setDetectorResolutionFunction(
- ConvolutionDetectorResolution::resolution_function_2d resolution_function)
+ ConvolutionDetectorResolution::cumulative_DF_2d resolution_function)
{
m_detector.setDetectorResolution(new ConvolutionDetectorResolution(resolution_function));
}
diff --git a/Core/Tools/inc/MathFunctions.h b/Core/Tools/inc/MathFunctions.h
index 1ac50026a62..377e04a2668 100644
--- a/Core/Tools/inc/MathFunctions.h
+++ b/Core/Tools/inc/MathFunctions.h
@@ -26,7 +26,7 @@ namespace MathFunctions
{
double Gaussian(double value, double average, double std_dev);
-double IntegratedGaussian(double value, double average, double std_dev, double step);
+double IntegratedGaussian(double value, double average, double std_dev);
double GenerateNormalRandom(double average, double std_dev);
diff --git a/Core/Tools/src/MathFunctions.cpp b/Core/Tools/src/MathFunctions.cpp
index 4c9f8644688..c9d0af01ee9 100644
--- a/Core/Tools/src/MathFunctions.cpp
+++ b/Core/Tools/src/MathFunctions.cpp
@@ -12,12 +12,11 @@ double MathFunctions::Gaussian(double value, double average, double std_dev)
return StandardNormal((value-average)/std_dev)/std_dev;
}
-double MathFunctions::IntegratedGaussian(double value, double average, double std_dev, double step)
+double MathFunctions::IntegratedGaussian(double value, double average, double std_dev)
{
- double left_margin = (value - average - step/2)/std_dev;
- double right_margin = (value - average + step/2)/std_dev;
+ double normalized_value = (value - average)/std_dev;
double root2 = std::sqrt(2.0);
- return (gsl_sf_erf(right_margin/root2) - gsl_sf_erf(left_margin/root2))/2.0;
+ return (gsl_sf_erf(normalized_value/root2) + 1.0)/2.0;
}
double MathFunctions::StandardNormal(double value)
--
GitLab