#include <math.h>
#include <iostream>
#include <fstream>
#include <boost/math/special_functions/round.hpp>
#include <boost/make_shared.hpp>
#include "MantidDataObjects/NoShape.h"
#include "MantidDataObjects/PeakShapeEllipsoid.h"
#include "MantidMDAlgorithms/Integrate3DEvents.h"
extern "C" {
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
}
using namespace Mantid::DataObjects;
namespace Mantid {
namespace MDAlgorithms {
using namespace std;
using Mantid::Kernel::DblMatrix;
using Mantid::Kernel::V3D;
/**
* Construct an object to store events that correspond to a peak an are
* within the specified radius of the specified peak centers, and to
* integrate the peaks.
*
* @param peak_q_list List of Q-vectors for peak centers.
* @param UBinv The matrix that maps Q-vectors to h,k,l
* @param radius The maximum distance from a peak's Q-vector, for
* an event to be stored in the list associated with
* that peak.
*/
Integrate3DEvents::Integrate3DEvents(
std::vector<std::pair<double, V3D>> const &peak_q_list,
DblMatrix const &UBinv, double radius) {
this->UBinv = UBinv;
this->radius = radius;
int64_t hkl_key;
for (size_t it = 0; it != peak_q_list.size(); ++it) {
hkl_key = getHklKey(peak_q_list[it].second);
if (hkl_key != 0) // only save if hkl != (0,0,0)
peak_qs[hkl_key] = peak_q_list[it].second;
}
}
/**
* Destructor.
*/
Integrate3DEvents::~Integrate3DEvents() {}
/**
* Add the specified event Q's to lists of events near peaks. An event is
* added to at most one list. First the nearest h,k,l for that event Q vector
* is calculated. If a peak with that h,k,l was specified when this object
* was constructed and if the distance from the specified event Q to that
* peak is less than the radius that was specified at construction time,
* then the event Q vector is added to the list of event Q vectors for that
* peak.
* NOTE: The Q-vectors passed in to this method will be shifted by the center
* Q for it's associated peak, so that the list of Q-vectors for a peak
* are centered around 0,0,0 and represent offsets in Q from the peak
* center.
*
* @param event_qs List of event Q vectors to add to lists of Q's associated
* with peaks.
* @param hkl_integ
*/
void Integrate3DEvents::addEvents(
std::vector<std::pair<double, V3D>> const &event_qs, bool hkl_integ) {
for (size_t i = 0; i < event_qs.size(); i++) {
addEvent(event_qs[i], hkl_integ);
}
}
/**
* Integrate the events around the specified peak Q-vector. The principal
* axes of the events near this Q-vector and the standard deviations in the
* directions of these principal axes determine ellipsoidal regions
* for integrating the peak and estimating the background. Alternatively,
* if peak and background radii are specified, then those will be used for
* half the major axis length of the ellipsoids, and the other axes of the
* ellipsoids will be set proportionally, based on the standard deviations.
*
* @param E1Vec Vector of values for calculating edge of detectors
* @param peak_q The Q-vector for the peak center.
* @param specify_size If true the integration will be done using the
* ellipsoids with major axes determined by the
* peak, back_inner and back_outer radii
* parameters. If false, the integration will be
* done using a peak region with major axis chosen
* so that it covers +- three standard deviations
* of the data in each direction. In this case,
* the background ellipsoidal shell is chosen to
* have the same VOLUME as the peak ellipsoid, and
* to use the peak ellipsoid for the inner
* radius.
* @param peak_radius Size of half the major axis of the ellipsoidal
* peak region.
* @param back_inner_radius Size of half the major axis of the INNER
* ellipsoidal boundary of the background region
* @param back_outer_radius Size of half the major axis of the OUTER
* ellipsoidal boundary of the background region
*
* @param axes_radii The radii used for integration in the
* directions of the three principal axes.
* @param inti Returns the net integrated intensity
* @param sigi Returns an estimate of the standard deviation
* of the net integrated intensity
*
*/
Mantid::Geometry::PeakShape_const_sptr
Integrate3DEvents::ellipseIntegrateEvents(
std::vector<Kernel::V3D> E1Vec, V3D const &peak_q, bool specify_size,
double peak_radius, double back_inner_radius, double back_outer_radius,
std::vector<double> &axes_radii, double &inti, double &sigi) {
inti = 0.0; // default values, in case something
sigi = 0.0; // is wrong with the peak.
int64_t hkl_key = getHklKey(peak_q);
if (hkl_key == 0) {
return boost::make_shared<NoShape>();
}
std::vector<std::pair<double, V3D>> &some_events = event_lists[hkl_key];
if (some_events.size() < 3) // if there are not enough events to
{ // find covariance matrix, return
return boost::make_shared<NoShape>();
}
DblMatrix cov_matrix(3, 3);
makeCovarianceMatrix(some_events, cov_matrix, radius);
std::vector<V3D> eigen_vectors;
getEigenVectors(cov_matrix, eigen_vectors);
std::vector<double> sigmas;
for (int i = 0; i < 3; i++) {
sigmas.push_back(stdDev(some_events, eigen_vectors[i], radius));
}
bool invalid_peak = false;
for (int i = 0; i < 3; i++) {
if ((boost::math::isnan)(sigmas[i])) {
invalid_peak = true;
} else if (sigmas[i] <= 0) {
invalid_peak = true;
}
}
if (invalid_peak) // if data collapses to a line or
{ // to a plane, the volume of the
return boost::make_shared<NoShape>(); // ellipsoids will be zero.
}
return ellipseIntegrateEvents(E1Vec, peak_q, some_events, eigen_vectors,
sigmas, specify_size, peak_radius,
back_inner_radius, back_outer_radius,
axes_radii, inti, sigi);
}
/**
* Calculate the number of events in an ellipsoid centered at 0,0,0 with
* the three specified axes and the three specified sizes in the direction
* of those axes. NOTE: The three axes must be mutually orthogonal unit
* vectors.
*
* @param events List of 3D events centered at 0,0,0
* @param directions List of 3 orthonormal directions for the axes of
* the ellipsoid.
* @param sizes List of three values a,b,c giving half the length
* of the three axes of the ellisoid.
* @return Then number of events that are in or on the specified ellipsoid.
*/
double Integrate3DEvents::numInEllipsoid(
std::vector<std::pair<double, V3D>> const &events,
std::vector<V3D> const &directions, std::vector<double> const &sizes) {
double count = 0;
for (size_t i = 0; i < events.size(); i++) {
double sum = 0;
for (size_t k = 0; k < 3; k++) {
double comp = events[i].second.scalar_prod(directions[k]) / sizes[k];
sum += comp * comp;
}
if (sum <= 1)
count += events[i].first;
}
return count;
}
/**
* Calculate the number of events in an ellipsoid centered at 0,0,0 with
* the three specified axes and the three specified sizes in the direction
* of those axes. NOTE: The three axes must be mutually orthogonal unit
* vectors.
*
* @param events List of 3D events centered at 0,0,0
* @param directions List of 3 orthonormal directions for the axes of
* the ellipsoid.
* @param sizes List of three values a,b,c giving half the length
* of the three axes of the ellisoid.
* @param sizesIn List of three values a,b,c giving half the length
* of the three inner axes of the ellisoid.
* @return Then number of events that are in or on the specified ellipsoid.
*/
double Integrate3DEvents::numInEllipsoidBkg(
std::vector<std::pair<double, V3D>> const &events,
std::vector<V3D> const &directions, std::vector<double> const &sizes,
std::vector<double> const &sizesIn) {
double count = 0;
std::vector<double> eventVec;
for (size_t i = 0; i < events.size(); i++) {
double sum = 0;
double sumIn = 0;
for (size_t k = 0; k < 3; k++) {
double comp = events[i].second.scalar_prod(directions[k]) / sizes[k];
sum += comp * comp;
comp = events[i].second.scalar_prod(directions[k]) / sizesIn[k];
sumIn += comp * comp;
}
if (sum <= 1 && sumIn >= 1)
eventVec.push_back(events[i].first);
}
std::sort(eventVec.begin(), eventVec.end());
// Remove top 1% of background
for (size_t k = 0;
k < static_cast<size_t>(0.99 * static_cast<double>(eventVec.size()));
k++) {
count += eventVec[k];
}
return count;
}
/**
* Given a list of events, associated with a particular peak
* and already SHIFTED to be centered at (0,0,0), calculate the 3x3
* covariance matrix for finding the principal axes of that
* local event data. Only events within the specified radius
* of (0,0,0) will be used.
*
* The covariance matrix can be easily constructed. X, Y, Z of each peak
*position are the variables we wish to determine
* the covariance. The mean position in each dimension has already been
*calculated on subtracted, since this corresponds to the centre position of
*each
* peak, which we knew aprori. The expected values of each correlation test X,X
*X,Y X,Z e.t.c form the elements of this 3 by 3 matrix, but since the
* probabilities are equal, we can remove them from the sums of the expected
*values, and simply divide by the number of events for each matrix element.
* Note that the diagonal elements form the variance X,X, Y,Y, Z,Z
*
* @param events Vector of V3D objects containing the
* Q vectors for a peak, with mean at (0,0,0).
* @param matrix A 3x3 matrix that will be filled out with
* the covariance matrix for the list of
* events.
* @param radius Only events within this radius of the
* peak center (0,0,0) will be used for
* calculating the covariance matrix.
*/
void Integrate3DEvents::makeCovarianceMatrix(
std::vector<std::pair<double, V3D>> const &events, DblMatrix &matrix,
double radius) {
for (int row = 0; row < 3; row++) {
for (int col = 0; col < 3; col++) {
double sum = 0;
for (size_t i = 0; i < events.size(); i++) {
auto event = events[i].second;
if (event.norm() <= radius) {
sum += event[row] * event[col];
}
}
if (events.size() > 1)
matrix[row][col] = sum / (double)(events.size() - 1);
else
matrix[row][col] = sum;
}
}
}
/**
* Calculate the eigen vectors of a 3x3 real symmetric matrix using the GSL.
*
* @param cov_matrix 3x3 real symmetric matrix.
* @param eigen_vectors The eigen vectors for the matrix are returned
* in this list.
*/
void Integrate3DEvents::getEigenVectors(DblMatrix const &cov_matrix,
std::vector<V3D> &eigen_vectors) {
unsigned int size = 3;
gsl_matrix *matrix = gsl_matrix_alloc(size, size);
gsl_vector *eigen_val = gsl_vector_alloc(size);
gsl_matrix *eigen_vec = gsl_matrix_alloc(size, size);
gsl_eigen_symmv_workspace *wkspace = gsl_eigen_symmv_alloc(size);
// copy the matrix data into the gsl matrix
for (size_t row = 0; row < size; row++)
for (size_t col = 0; col < size; col++) {
gsl_matrix_set(matrix, row, col, cov_matrix[row][col]);
}
gsl_eigen_symmv(matrix, eigen_val, eigen_vec, wkspace);
// copy the resulting eigen vectors to output vector
for (size_t col = 0; col < size; col++) {
eigen_vectors.push_back(V3D(gsl_matrix_get(eigen_vec, 0, col),
gsl_matrix_get(eigen_vec, 1, col),
gsl_matrix_get(eigen_vec, 2, col)));
}
gsl_matrix_free(matrix);
gsl_vector_free(eigen_val);
gsl_matrix_free(eigen_vec);
gsl_eigen_symmv_free(wkspace);
}
/**
* Calculate the standard deviation of the given list of 3D events in
* the direction of the specified vector. Only events that are within
* the specified radius of 0,0,0 will be considered.
*
* @param events List of 3D events centered at 0,0,0
* @param direction Unit vector giving the direction vector on which
* the 3D events will be projected.
* @param radius Maximun size of event vectors that will be used
* in calculating the standard deviation.
*/
double
Integrate3DEvents::stdDev(std::vector<std::pair<double, V3D>> const &events,
V3D const &direction, double radius) {
double sum = 0;
double sum_sq = 0;
double stdev = 0;
int count = 0;
for (size_t i = 0; i < events.size(); i++) {
auto event = events[i].second;
if (event.norm() <= radius) {
double dot_prod = event.scalar_prod(direction);
sum += dot_prod;
sum_sq += dot_prod * dot_prod;
count++;
}
}
if (count > 1) {
double ave = sum / count;
stdev = sqrt((sum_sq / count - ave * ave) * (double)count / (count - 1.0));
}
return stdev;
}
/**
* Form a map key as 10^12*h + 10^6*k + l from the integers h,k,l.
*
* @param h The first Miller index
* @param k The second Miller index
* @param l The third Miller index
*/
int64_t Integrate3DEvents::getHklKey(int h, int k, int l) {
int64_t key(0);
if (h != 0 || k != 0 || l != 0)
key = 1000000000000 * h + 1000000 * k + l;
return key;
}
/**
* Form a map key for the specified q_vector. The q_vector is mapped to
* h,k,l by UBinv and the map key is then formed from those rounded h,k,l
* values.
*
* @param q_vector The q_vector to be mapped to h,k,l
*/
int64_t Integrate3DEvents::getHklKey2(V3D const &hkl) {
int h = boost::math::iround<double>(hkl[0]);
int k = boost::math::iround<double>(hkl[1]);
int l = boost::math::iround<double>(hkl[2]);
return getHklKey(h, k, l);
}
/**
* Form a map key for the specified q_vector. The q_vector is mapped to
* h,k,l by UBinv and the map key is then formed from those rounded h,k,l
* values.
*
* @param q_vector The q_vector to be mapped to h,k,l
*/
int64_t Integrate3DEvents::getHklKey(V3D const &q_vector) {
V3D hkl = UBinv * q_vector;
int h = boost::math::iround<double>(hkl[0]);
int k = boost::math::iround<double>(hkl[1]);
int l = boost::math::iround<double>(hkl[2]);
return getHklKey(h, k, l);
}
/**
* Add an event to the appropriate vector of events for the closest h,k,l,
* if it is within the required radius of the corresponding peak in the
* PeakQMap.
*
* NOTE: The event passed in may be modified by this method. In particular,
* if it corresponds to one of the specified peak_qs, the corresponding peak q
* will be subtracted from the event and the event will be added to that
* peak's vector in the event_lists map.
*
* @param event_Q The Q-vector for the event that may be added to the
* event_lists map, if it is close enough to some peak
* @param hkl_integ
*/
void Integrate3DEvents::addEvent(std::pair<double, V3D> event_Q,
bool hkl_integ) {
int64_t hkl_key;
if (hkl_integ)
hkl_key = getHklKey2(event_Q.second);
else
hkl_key = getHklKey(event_Q.second);
if (hkl_key == 0) // don't keep events associated with 0,0,0
return;
auto peak_it = peak_qs.find(hkl_key);
if (peak_it != peak_qs.end()) {
if (!peak_it->second.nullVector()) {
if (hkl_integ)
event_Q.second = event_Q.second - UBinv * peak_it->second;
else
event_Q.second = event_Q.second - peak_it->second;
if (event_Q.second.norm() < radius) {
event_lists[hkl_key].push_back(event_Q);
}
}
}
}
/**
* Integrate a list of events, centered about (0,0,0) given the principal
* axes for the events and the standard deviations in the the directions
* of the principal axes.
*
* @param E1Vec Vector of values for calculating edge of detectors
* @param peak_q The Q-vector for the peak center.
* @param ev_list List of events centered at (0,0,0) for a
* particular peak.
* @param directions The three principal axes of the list of events
* @param sigmas The standard deviations of the events in the
* directions of the three principal axes.
* @param specify_size If true the integration will be done using the
* ellipsoids with major axes determined by the
* peak, back_inner and back_outer radii
* parameters. If false, the integration will be
* done using a peak region with major axis chosen
* so that it covers +- three standard deviations
* of the data in each direction. In this case,
* the background ellipsoidal shell is chosen to
* have the same VOLUME as the peak ellipsoid, and
* to use the peak ellipsoid for the inner
* radius.
* @param peak_radius Size of half the major axis of the ellipsoidal
* peak region.
* @param back_inner_radius Size of half the major axis of the INNER
* ellipsoidal boundary of the background region
* @param back_outer_radius Size of half the major axis of the OUTER
* ellipsoidal boundary of the background region
*
* @param axes_radii The radii used for integration in the
* directions of the three principal axes.
* @param inti Returns the net integrated intensity
* @param sigi Returns an estimate of the standard deviation
* of the net integrated intensity
*
*/
PeakShapeEllipsoid_const_sptr Integrate3DEvents::ellipseIntegrateEvents(
std::vector<Kernel::V3D> E1Vec, V3D const &peak_q,
std::vector<std::pair<double, Mantid::Kernel::V3D>> const &ev_list,
std::vector<Mantid::Kernel::V3D> const &directions,
std::vector<double> const &sigmas, bool specify_size, double peak_radius,
double back_inner_radius, double back_outer_radius,
std::vector<double> &axes_radii, double &inti, double &sigi) {
// r1, r2 and r3 will give the sizes of the major axis of
// the peak ellipsoid, and of the inner and outer surface
// of the background ellipsoidal shell, respectively.
// They will specify the size as the number of standard
// deviations in the direction of each of the pricipal
// axes that the ellipsoid will extend from the center.
double r1, r2, r3;
double max_sigma = sigmas[0];
for (int i = 1; i < 3; i++) {
if (sigmas[i] > max_sigma) {
max_sigma = sigmas[i];
}
}
if (specify_size) {
r1 = peak_radius / max_sigma; // scale specified sizes by 1/max_sigma
r2 = back_inner_radius / max_sigma; // so when multiplied by the individual
r3 = back_outer_radius / max_sigma; // sigmas in different directions, the
} // major axis has the specified size
else {
r1 = 3;
r2 = 3;
r3 = r2 * 1.25992105; // A factor of 2 ^ (1/3) will make the background
// shell volume equal to the peak region volume.
// if necessary restrict the background ellipsoid
// to lie within the specified sphere, and adjust
// the other sizes, proportionally
if (r3 * max_sigma > radius) {
r3 = radius / max_sigma;
r1 = r3 * 0.79370053f; // This value for r1 and r2 makes the background
r2 = r1; // shell volume equal to the peak region volume.
}
}
axes_radii.clear();
std::vector<double> abcBackgroundOuterRadii;
std::vector<double> abcBackgroundInnerRadii;
std::vector<double> abcRadii;
for (int i = 0; i < 3; i++) {
abcBackgroundOuterRadii.push_back(r3 * sigmas[i]);
abcBackgroundInnerRadii.push_back(r2 * sigmas[i]);
abcRadii.push_back(r1 * sigmas[i]);
axes_radii.push_back(r1 * sigmas[i]);
}
if (E1Vec.size() > 0) {
double h3 = 1.0 - detectorQ(E1Vec, peak_q, abcBackgroundOuterRadii);
// scaled from area of circle minus segment when r normalized to 1
double m3 = std::sqrt(
1.0 -
(std::acos(1.0 - h3) - (1.0 - h3) * std::sqrt(2.0 * h3 - h3 * h3)) /
M_PI);
double h1 = 1.0 - detectorQ(E1Vec, peak_q, axes_radii);
// Do not use peak if edge of detector is inside integration radius
if (h1 > 0.0)
return boost::make_shared<const PeakShapeEllipsoid>(
directions, abcRadii, abcBackgroundInnerRadii,
abcBackgroundOuterRadii, Mantid::Kernel::QLab, "IntegrateEllipsoids");
r3 *= m3;
if (r2 != r1) {
double h2 = 1.0 - detectorQ(E1Vec, peak_q, abcBackgroundInnerRadii);
// scaled from area of circle minus segment when r normalized to 1
double m2 = std::sqrt(
1.0 -
(std::acos(1.0 - h2) - (1.0 - h2) * std::sqrt(2.0 * h2 - h2 * h2)) /
M_PI);
r2 *= m2;
}
}
double backgrd = numInEllipsoidBkg(
ev_list, directions, abcBackgroundOuterRadii, abcBackgroundInnerRadii);
double peak_w_back = numInEllipsoid(ev_list, directions, axes_radii);
double ratio = pow(r1, 3) / (pow(r3, 3) - pow(r2, 3));
inti = peak_w_back - ratio * backgrd;
sigi = sqrt(peak_w_back + ratio * ratio * backgrd);
// Make the shape and return it.
return boost::make_shared<const PeakShapeEllipsoid>(
directions, abcRadii, abcBackgroundInnerRadii, abcBackgroundOuterRadii,
Mantid::Kernel::QLab, "IntegrateEllipsoids");
}
/** Calculate if this Q is on a detector
* The distance from C to OE is given by dv=C-E*(C.scalar_prod(E))
* If dv.norm<integration_radius, one of the detector trajectories on the edge
*is too close to the peak
* This method is applied to all masked pixels. If there are masked pixels
*trajectories inside an integration volume, the peak must be rejected.
*
* @param E1Vec Vector of values for calculating edge of detectors
* @param QLabFrame: The Peak center.
* @param r: Peak radius.
*/
double Integrate3DEvents::detectorQ(std::vector<Kernel::V3D> E1Vec,
const Mantid::Kernel::V3D QLabFrame,
std::vector<double> &r) {
double quot = 1.0;
for (auto E1 = E1Vec.begin(); E1 != E1Vec.end(); ++E1) {
V3D distv = QLabFrame -
*E1 * (QLabFrame.scalar_prod(
*E1)); // distance to the trajectory as a vector
double quot0 = distv.norm() / *(std::min_element(r.begin(), r.end()));
if (quot0 < quot) {
quot = quot0;
}
}
return quot;
}
} // namespace MDAlgorithms
} // namespace Mantid