15using Geometry::CSGObject;
19 l1.reserve(numVolumeElements);
20 volume.reserve(numVolumeElements);
26constexpr static V3D X_AXIS{1, 0, 0};
27constexpr static V3D Y_AXIS{0, 1, 0};
28constexpr static V3D Z_AXIS{0, 0, 1};
30struct CylinderParameters {
37struct HollowCylinderParameters {
47V3D createPerpendicular(
const V3D &symmetryAxis) {
48 const std::array<double, 3> scalars = {{
fabs(symmetryAxis.scalar_prod(X_AXIS)),
49 fabs(symmetryAxis.scalar_prod(Y_AXIS)),
50 fabs(symmetryAxis.scalar_prod(Z_AXIS))}};
53 return symmetryAxis.cross_prod(X_AXIS);
54 else if (scalars[1] == 0.)
55 return symmetryAxis.cross_prod(Y_AXIS);
56 else if (scalars[2] == 0.)
57 return symmetryAxis.cross_prod(Z_AXIS);
60 if (scalars[0] < scalars[1] && scalars[0] < scalars[2])
61 return symmetryAxis.cross_prod(X_AXIS);
62 else if (scalars[1] < scalars[0] && scalars[1] < scalars[2])
63 return symmetryAxis.cross_prod(Y_AXIS);
65 return symmetryAxis.cross_prod(Z_AXIS);
68const V3D CalculatePosInCylinder(
const double phi,
const double R,
const double z,
const std::array<V3D, 3> coords) {
69 const auto &x_prime = coords[0];
70 const auto &y_prime = coords[1];
71 const auto &z_prime = coords[2];
73 double rSinPhi = -R * sin(phi);
74 const double rCosPhi = -R * cos(phi);
77 const double xcomp = rCosPhi * x_prime[0] + rSinPhi * y_prime[0] +
z * z_prime[0];
78 const double ycomp = rCosPhi * x_prime[1] + rSinPhi * y_prime[1] +
z * z_prime[1];
79 const double zcomp = rCosPhi * x_prime[2] + rSinPhi * y_prime[2] +
z * z_prime[2];
80 return V3D(xcomp, ycomp, zcomp);
93double calcDistanceInShapeNoCheck(
const V3D &beamDirection,
const IObject &shape,
const V3D &
position) {
95 Track incoming(
position, -beamDirection);
97 if (shape.interceptSurface(incoming) > 0) {
98 return incoming.totalDistInsideObject();
104Raster calculateGeneric(
const V3D &beamDirection,
const IObject &shape,
const double cubeSizeInMetre) {
105 if (cubeSizeInMetre <= 0.)
106 throw std::runtime_error(
"Tried to section shape into zero size elements");
108 const auto bbox = shape.getBoundingBox();
109 assert(bbox.xMax() > bbox.xMin());
110 assert(bbox.yMax() > bbox.yMin());
111 assert(bbox.zMax() > bbox.zMin());
113 const double xLength = bbox.xMax() - bbox.xMin();
114 const double yLength = bbox.yMax() - bbox.yMin();
115 const double zLength = bbox.zMax() - bbox.zMin();
117 const auto numXSlices =
static_cast<size_t>(xLength / cubeSizeInMetre);
118 const auto numYSlices =
static_cast<size_t>(yLength / cubeSizeInMetre);
119 const auto numZSlices =
static_cast<size_t>(zLength / cubeSizeInMetre);
120 const double XSliceThickness = xLength /
static_cast<double>(numXSlices);
121 const double YSliceThickness = yLength /
static_cast<double>(numYSlices);
122 const double ZSliceThickness = zLength /
static_cast<double>(numZSlices);
123 const double elementVolume = XSliceThickness * YSliceThickness * ZSliceThickness;
125 const size_t numVolumeElements = numXSlices * numYSlices * numZSlices;
129 result.reserve(numVolumeElements);
133 throw std::logic_error(
"Too many volume elements requested - try increasing the value "
134 "of the ElementSize property.");
139 for (
size_t i = 0; i < numZSlices; ++i) {
140 const double z = (
static_cast<double>(i) + 0.5) * ZSliceThickness + bbox.xMin();
142 for (
size_t j = 0; j < numYSlices; ++j) {
143 const double y = (
static_cast<double>(j) + 0.5) * YSliceThickness + bbox.yMin();
145 for (
size_t k = 0; k < numXSlices; ++k) {
146 const double x = (
static_cast<double>(k) + 0.5) * XSliceThickness + bbox.zMin();
148 const Kernel::V3D currentPosition = V3D(
x,
y,
z);
150 if (shape.isValid(currentPosition)) {
152 Track incoming(currentPosition, -beamDirection);
158 if (shape.interceptSurface(incoming) > 0) {
159 result.l1.emplace_back(incoming.totalDistInsideObject());
160 result.position.emplace_back(currentPosition);
161 result.volume.emplace_back(elementVolume);
169 result.totalvolume =
static_cast<double>(result.l1.size()) * elementVolume;
182 const auto primitive = shape.
shape();
183 if (hasCustomizedRaster(primitive)) {
186 const auto &shapeInfo = shape.
shapeInfo();
189 const size_t numSlice = std::max<size_t>(1,
static_cast<size_t>(params.height / cubeSizeInMetre));
190 const size_t numAnnuli = std::max<size_t>(1,
static_cast<size_t>(params.radius / cubeSizeInMetre));
193 const auto params = shapeInfo.hollowCylinderGeometry();
194 const size_t numSlice = std::max<size_t>(1,
static_cast<size_t>(params.height / cubeSizeInMetre));
195 const size_t numAnnuli =
196 std::max<size_t>(1,
static_cast<size_t>((params.radius - params.innerRadius) / cubeSizeInMetre));
199 throw std::runtime_error(
"Rasterize::calculate should never get to this point");
202 return calculateGeneric(beamDirection, shape, cubeSizeInMetre);
207 const size_t numAnnuli) {
209 throw std::invalid_argument(
"Given shape is not a cylinder.");
212 throw std::runtime_error(
"Tried to section cylinder into zero slices");
214 throw std::runtime_error(
"Tried to section cylinder into zero annuli");
217 const auto &shapeInfo = shape.
shapeInfo();
221 const V3D center = (params.axis * .5 * params.height) + params.centreOfBottomBase;
223 const double sliceThickness{params.height /
static_cast<double>(numSlices)};
224 const double deltaR{params.radius /
static_cast<double>(numAnnuli)};
231 const size_t numVolumeElements = numSlices * numAnnuli * (numAnnuli + 1) * 3;
234 result.
reserve(numVolumeElements);
235 result.totalvolume = params.height * M_PI * params.radius * params.radius;
240 const V3D x_prime = createPerpendicular(z_prime);
242 const auto coords = std::array<V3D, 3>{{x_prime, y_prime, z_prime}};
246 for (
size_t i = 0; i < numSlices; ++i) {
247 const double z = (
static_cast<double>(i) + 0.5) * sliceThickness - 0.5 * params.height;
252 for (
size_t j = 0; j < numAnnuli; ++j) {
254 const double R = (
static_cast<double>(j) * params.radius /
static_cast<double>(numAnnuli)) + (0.5 * deltaR);
257 const double outerR = R + (deltaR / 2.0);
258 const double innerR = outerR - deltaR;
259 const double elementVolume =
260 M_PI * (outerR * outerR - innerR * innerR) * sliceThickness /
static_cast<double>(Ni);
263 for (
size_t k = 0; k < Ni; ++k) {
264 const double phi = 2. * M_PI *
static_cast<double>(k) /
static_cast<double>(Ni);
265 const auto position = center + CalculatePosInCylinder(phi, R,
z, coords);
269 result.position.emplace_back(
position);
270 result.volume.emplace_back(elementVolume);
272 result.l1.emplace_back(calcDistanceInShapeNoCheck(beamDirection, shape,
position));
281 const size_t numAnnuli) {
283 throw std::invalid_argument(
"Given shape is not a hollow cylinder.");
286 throw std::runtime_error(
"Tried to section cylinder into zero slices");
288 throw std::runtime_error(
"Tried to section cylinder into zero annuli");
291 const auto &shapeInfo = shape.
shapeInfo();
295 const V3D center = (params.axis * .5 * params.height) + params.centreOfBottomBase;
297 const double sliceThickness{params.height /
static_cast<double>(numSlices)};
298 const double deltaR{(params.radius - params.innerRadius) /
static_cast<double>(numAnnuli)};
305 const size_t numVolumeElements = numSlices * numAnnuli * (numAnnuli + 1) * 3;
308 result.
reserve(numVolumeElements);
309 result.totalvolume = params.height * M_PI * (params.radius * params.radius - params.innerRadius * params.innerRadius);
313 V3D z_prime = params.axis;
315 const V3D x_prime = createPerpendicular(z_prime);
317 const auto coords = std::array<V3D, 3>{{x_prime, y_prime, z_prime}};
321 for (
size_t i = 0; i < numSlices; ++i) {
322 const double z = (
static_cast<double>(i) + 0.5) * sliceThickness - 0.5 * params.height;
332 const auto nSteps = params.innerRadius / deltaR;
333 size_t Ni =
static_cast<size_t>(nSteps) * 6;
335 for (
size_t j = 0; j < numAnnuli; ++j) {
339 (
static_cast<double>(j) * (params.radius - params.innerRadius) /
static_cast<double>(numAnnuli)) +
343 const double outerR = R + (deltaR / 2.0);
344 const double innerR = outerR - deltaR;
345 const double elementVolume =
346 M_PI * (outerR * outerR - innerR * innerR) * sliceThickness /
static_cast<double>(Ni);
349 for (
size_t k = 0; k < Ni; ++k) {
350 const double phi = 2. * M_PI *
static_cast<double>(k) /
static_cast<double>(Ni);
352 const auto position = center + CalculatePosInCylinder(phi, R,
z, coords);
356 result.position.emplace_back(
position);
357 result.volume.emplace_back(elementVolume);
359 result.l1.emplace_back(calcDistanceInShapeNoCheck(beamDirection, shape,
position));
IObject : Interface for geometry objects.
virtual const detail::ShapeInfo & shapeInfo() const =0
virtual detail::ShapeInfo::GeometryShape shape() const =0
@ HOLLOWCYLINDER
HOLLOW CYLINDER.
CylinderGeometry cylinderGeometry() const
HollowCylinderGeometry hollowCylinderGeometry() const
double normalize()
Make a normalized vector (return norm value)
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
MANTID_GEOMETRY_DLL Raster calculate(const Kernel::V3D &beamDirection, const IObject &shape, const double cubeSizeInMetre)
MANTID_GEOMETRY_DLL Raster calculateCylinder(const Kernel::V3D &beamDirection, const IObject &shape, const size_t numSlices, const size_t numAnnuli)
MANTID_GEOMETRY_DLL Raster calculateHollowCylinder(const Kernel::V3D &beamDirection, const IObject &shape, const size_t numSlices, const size_t numAnnuli)
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Holds the information used for doing numerical integrations of object in the beam.
std::vector< double > volume
Cached element volumes.
std::vector< double > l1
Cached L1 distances.
void reserve(size_t numVolumeElements)
std::vector< Kernel::V3D > position
Cached element positions.