RandomPoint.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2019 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidGeometry/RandomPoint.h"
#include "MantidGeometry/Objects/IObject.h"
#include "MantidKernel/PseudoRandomNumberGenerator.h"
#include <boost/math/special_functions/pow.hpp>
 
namespace Mantid::Geometry::RandomPoint {
 
Kernel::V3D localPointInCylinder(const Kernel::V3D &basis, const Kernel::V3D &alongAxis, double polarAngle,
                                 double radialLength) {
  // Use basis to get a second perpendicular vector to define basis2
  Kernel::V3D basis2;
  if (basis.X() == 0) {
    basis2.setX(1.);
  } else if (basis.Y() == 0) {
    basis2.setY(1.);
  } else if (basis.Z() == 0) {
    basis2.setZ(1.);
  } else {
    basis2.setX(-basis.Y());
    basis2.setY(basis.X());
    basis2.normalize();
  }
  const Kernel::V3D basis3{basis.cross_prod(basis2)};
  const Kernel::V3D localPoint{((basis2 * std::cos(polarAngle) + basis3 * std::sin(polarAngle)) * radialLength) +
                               alongAxis};
  return localPoint;
}
 
Kernel::V3D inCuboid(const detail::ShapeInfo &shapeInfo, Kernel::PseudoRandomNumberGenerator &rng) {
  const auto geometry = shapeInfo.cuboidGeometry();
  const double r1{rng.nextValue()};
  const double r2{rng.nextValue()};
  const double r3{rng.nextValue()};
  const Kernel::V3D basis1{geometry.leftFrontTop - geometry.leftFrontBottom};
  const Kernel::V3D basis2{geometry.leftBackBottom - geometry.leftFrontBottom};
  const Kernel::V3D basis3{geometry.rightFrontBottom - geometry.leftFrontBottom};
  return geometry.leftFrontBottom + (basis1 * r1 + basis2 * r2 + basis3 * r3);
}
 
Kernel::V3D inCylinder(const detail::ShapeInfo &shapeInfo, Kernel::PseudoRandomNumberGenerator &rng) {
  const auto geometry = shapeInfo.cylinderGeometry();
  const double r1{rng.nextValue()};
  const double r2{rng.nextValue()};
  const double r3{rng.nextValue()};
  const double polar{2. * M_PI * r1};
  // The sqrt is needed for a uniform distribution of points.
  const double r{geometry.radius * std::sqrt(r2)};
  const double z{geometry.height * r3};
  const Kernel::V3D alongAxis{geometry.axis * z};
  auto localPoint = localPointInCylinder(geometry.axis, alongAxis, polar, r);
  return localPoint + geometry.centreOfBottomBase;
}
 
Kernel::V3D inHollowCylinder(const detail::ShapeInfo &shapeInfo, Kernel::PseudoRandomNumberGenerator &rng) {
  const auto geometry = shapeInfo.hollowCylinderGeometry();
  const double r1{rng.nextValue()};
  const double r2{rng.nextValue()};
  const double r3{rng.nextValue()};
  const double polar{2. * M_PI * r1};
  // We need a random number between the inner radius and outer radius, but also
  // need the square root for a uniform distribution of points
  const double c1 = geometry.innerRadius * geometry.innerRadius;
  const double c2 = geometry.radius * geometry.radius;
  const double r{std::sqrt(c1 + (c2 - c1) * r2)};
  const double z{geometry.height * r3};
  const Kernel::V3D alongAxis{geometry.axis * z};
  auto localPoint = localPointInCylinder(geometry.axis, alongAxis, polar, r);
  return localPoint + geometry.centreOfBottomBase;
}
 
Kernel::V3D inSphere(const detail::ShapeInfo &shapeInfo, Kernel::PseudoRandomNumberGenerator &rng) {
  const auto geometry = shapeInfo.sphereGeometry();
  const double r1{rng.nextValue()};
  const double r2{rng.nextValue()};
  const double r3{rng.nextValue()};
  const double azimuthal{2. * M_PI * r1};
  // The acos is needed for a uniform distribution of points.
  const double polar{std::acos(2. * r2 - 1.)};
  const double r{r3 * geometry.radius};
  const double x{r * std::cos(azimuthal) * std::sin(polar)};
  const double y{r * std::sin(azimuthal) * std::sin(polar)};
  const double z{r * std::cos(polar)};
  return geometry.centre + Kernel::V3D{x, y, z};
}
 
boost::optional<Kernel::V3D> inGenericShape(const IObject &object, Kernel::PseudoRandomNumberGenerator &rng,
                                            size_t maxAttempts) {
  return bounded(object, rng, object.getBoundingBox(), maxAttempts);
}
 
boost::optional<Kernel::V3D> bounded(const IObject &object, Kernel::PseudoRandomNumberGenerator &rng,
                                     const BoundingBox &box, size_t maxAttempts) {
  boost::optional<Kernel::V3D> point{boost::none};
  if (box.isNull()) {
    throw std::invalid_argument("Invalid bounding box. Cannot generate random point.");
  }
  for (size_t attempts{0}; attempts < maxAttempts; ++attempts) {
    const double r1 = rng.nextValue();
    const double r2 = rng.nextValue();
    const double r3 = rng.nextValue();
    auto pt = box.generatePointInside(r1, r2, r3);
    if (object.isValid(pt)) {
      point = pt;
      break;
    }
  };
  return point;
}
} // namespace Mantid::Geometry::RandomPoint