Mantid
Loading...
Searching...
No Matches
Rasterize.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
11#include <array>
12
13namespace Mantid::Geometry {
14
15using Geometry::CSGObject;
16using Kernel::V3D;
17
18void Raster::reserve(size_t numVolumeElements) {
19 l1.reserve(numVolumeElements);
20 volume.reserve(numVolumeElements);
21 position.reserve(numVolumeElements);
22}
23
24namespace { // anonymous
25
26constexpr static V3D X_AXIS{1, 0, 0};
27constexpr static V3D Y_AXIS{0, 1, 0};
28constexpr static V3D Z_AXIS{0, 0, 1};
29
30struct CylinderParameters {
31 double radius;
32 double height;
35};
36
37struct HollowCylinderParameters {
38 double radius;
40 double height;
42 V3D symmetryaxis;
43};
44
45// since cylinders are symmetric around the main axis, choose a random
46// perpendicular to have as the second axis
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))}};
51 // check against the cardinal axes
52 if (scalars[0] == 0.)
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);
58
59 // see which one is smallest
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);
64 else
65 return symmetryAxis.cross_prod(Z_AXIS);
66}
67
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];
72
73 double rSinPhi = -R * sin(phi);
74 const double rCosPhi = -R * cos(phi);
75
76 // Calculate the current position in the shape in Cartesian coordinates
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);
81}
82
83// This must be updated when new primitives become available
84bool hasCustomizedRaster(detail::ShapeInfo::GeometryShape shape) {
86 return true;
88 return true;
89 // nothing else has been customized
90 return false;
91}
92
93double calcDistanceInShapeNoCheck(const V3D &beamDirection, const IObject &shape, const V3D &position) {
94 // Create track for distance in cylinder before scattering point
95 Track incoming(position, -beamDirection);
96
97 if (shape.interceptSurface(incoming) > 0) {
98 return incoming.totalDistInsideObject();
99 } else {
100 return 0;
101 }
102}
103
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");
107
108 const auto bbox = shape.getBoundingBox();
109 assert(bbox.xMax() > bbox.xMin());
110 assert(bbox.yMax() > bbox.yMin());
111 assert(bbox.zMax() > bbox.zMin());
112
113 const double xLength = bbox.xMax() - bbox.xMin();
114 const double yLength = bbox.yMax() - bbox.yMin();
115 const double zLength = bbox.zMax() - bbox.zMin();
116
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;
124
125 const size_t numVolumeElements = numXSlices * numYSlices * numZSlices;
126
127 Raster result;
128 try {
129 result.reserve(numVolumeElements);
130 } catch (...) {
131 // Typically get here if the number of volume elements is too large
132 // Provide a bit more information
133 throw std::logic_error("Too many volume elements requested - try increasing the value "
134 "of the ElementSize property.");
135 }
136
137 // go through the bounding box generating cubes and seeing if they are
138 // inside the shape
139 for (size_t i = 0; i < numZSlices; ++i) {
140 const double z = (static_cast<double>(i) + 0.5) * ZSliceThickness + bbox.xMin();
141
142 for (size_t j = 0; j < numYSlices; ++j) {
143 const double y = (static_cast<double>(j) + 0.5) * YSliceThickness + bbox.yMin();
144
145 for (size_t k = 0; k < numXSlices; ++k) {
146 const double x = (static_cast<double>(k) + 0.5) * XSliceThickness + bbox.zMin();
147 // Set the current position in the sample in Cartesian coordinates.
148 const Kernel::V3D currentPosition = V3D(x, y, z);
149 // Check if the current point is within the object. If not, skip.
150 if (shape.isValid(currentPosition)) {
151 // Create track for distance in sample before scattering point
152 Track incoming(currentPosition, -beamDirection);
153 // We have an issue where occasionally, even though a point is
154 // within the object a track segment to the surface isn't correctly
155 // created. In the context of this algorithm I think it's safe to
156 // just chuck away the element in this case. This will also throw
157 // away points that are inside a gauge volume but outside the sample
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);
162 }
163 }
164 }
165 }
166 }
167
168 // Record the number of elements we ended up with
169 result.totalvolume = static_cast<double>(result.l1.size()) * elementVolume;
170
171 return result;
172}
173
174} // namespace
175
176namespace Rasterize {
177
178// -------------------
179// collection of calculations that convert to CSGObjects and pass the work on
180
181Raster calculate(const V3D &beamDirection, const IObject &shape, const double cubeSizeInMetre) {
182 const auto primitive = shape.shape();
183 if (hasCustomizedRaster(primitive)) {
184 // convert to the underlying primitive type - this assumes that there are
185 // only customizations for CSGObjects
186 const auto &shapeInfo = shape.shapeInfo();
188 const auto params = shapeInfo.cylinderGeometry();
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));
191 return calculateCylinder(beamDirection, shape, numSlice, numAnnuli);
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));
197 return calculateHollowCylinder(beamDirection, shape, numSlice, numAnnuli);
198 } else {
199 throw std::runtime_error("Rasterize::calculate should never get to this point");
200 }
201 } else {
202 return calculateGeneric(beamDirection, shape, cubeSizeInMetre);
203 }
204}
205
206Raster calculateCylinder(const V3D &beamDirection, const IObject &shape, const size_t numSlices,
207 const size_t numAnnuli) {
209 throw std::invalid_argument("Given shape is not a cylinder.");
210
211 if (numSlices == 0)
212 throw std::runtime_error("Tried to section cylinder into zero slices");
213 if (numAnnuli == 0)
214 throw std::runtime_error("Tried to section cylinder into zero annuli");
215
216 // convert to the underlying primitive type
217 const auto &shapeInfo = shape.shapeInfo();
218
219 // get the geometry for the volume elements
220 const auto params = shapeInfo.cylinderGeometry();
221 const V3D center = (params.axis * .5 * params.height) + params.centreOfBottomBase;
222
223 const double sliceThickness{params.height / static_cast<double>(numSlices)};
224 const double deltaR{params.radius / static_cast<double>(numAnnuli)};
225
226 /* The number of volume elements is
227 * numslices*(1+2+3+.....+numAnnuli)*6
228 * Since the first annulus is separated in 6 segments, the next one in 12 and
229 * so on.....
230 */
231 const size_t numVolumeElements = numSlices * numAnnuli * (numAnnuli + 1) * 3;
232
233 Raster result;
234 result.reserve(numVolumeElements);
235 result.totalvolume = params.height * M_PI * params.radius * params.radius;
236
237 // Assume that z' = axis. Then select whatever has the smallest dot product
238 // with axis to be the x' direction
239 const V3D z_prime = normalize(params.axis);
240 const V3D x_prime = createPerpendicular(z_prime);
241 const V3D y_prime = z_prime.cross_prod(x_prime);
242 const auto coords = std::array<V3D, 3>{{x_prime, y_prime, z_prime}};
243
244 // loop over the elements of the shape and create everything
245 // loop over slices
246 for (size_t i = 0; i < numSlices; ++i) {
247 const double z = (static_cast<double>(i) + 0.5) * sliceThickness - 0.5 * params.height;
248
249 // Number of elements in 1st annulus
250 size_t Ni = 0;
251 // loop over annuli
252 for (size_t j = 0; j < numAnnuli; ++j) {
253 Ni += 6;
254 const double R = (static_cast<double>(j) * params.radius / static_cast<double>(numAnnuli)) + (0.5 * deltaR);
255
256 // all the volume elements in the ring/slice are the same
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);
261
262 // loop over elements in current annulus
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);
266
267 assert(shape.isValid(position));
268
269 result.position.emplace_back(position);
270 result.volume.emplace_back(elementVolume);
271 // TODO should be customized for cylinder
272 result.l1.emplace_back(calcDistanceInShapeNoCheck(beamDirection, shape, position));
273 } // loop over k
274 } // loop over j
275 } // loop over i
276
277 return result;
278}
279
280Raster calculateHollowCylinder(const V3D &beamDirection, const IObject &shape, const size_t numSlices,
281 const size_t numAnnuli) {
283 throw std::invalid_argument("Given shape is not a hollow cylinder.");
284
285 if (numSlices == 0)
286 throw std::runtime_error("Tried to section cylinder into zero slices");
287 if (numAnnuli == 0)
288 throw std::runtime_error("Tried to section cylinder into zero annuli");
289
290 // convert to the underlying primitive type
291 const auto &shapeInfo = shape.shapeInfo();
292
293 // get the geometry for the volume elements
294 const auto params = shapeInfo.hollowCylinderGeometry();
295 const V3D center = (params.axis * .5 * params.height) + params.centreOfBottomBase;
296
297 const double sliceThickness{params.height / static_cast<double>(numSlices)};
298 const double deltaR{(params.radius - params.innerRadius) / static_cast<double>(numAnnuli)};
299
300 /* The number of volume elements is
301 * numslices*(1+2+3+.....+numAnnuli)*6
302 * Since the first annulus is separated in 6 segments, the next one in 12 and
303 * so on.....
304 */
305 const size_t numVolumeElements = numSlices * numAnnuli * (numAnnuli + 1) * 3;
306
307 Raster result;
308 result.reserve(numVolumeElements);
309 result.totalvolume = params.height * M_PI * (params.radius * params.radius - params.innerRadius * params.innerRadius);
310
311 // Assume that z' = axis. Then select whatever has the smallest dot product
312 // with axis to be the x' direction
313 V3D z_prime = params.axis;
314 z_prime.normalize();
315 const V3D x_prime = createPerpendicular(z_prime);
316 const V3D y_prime = z_prime.cross_prod(x_prime);
317 const auto coords = std::array<V3D, 3>{{x_prime, y_prime, z_prime}};
318
319 // loop over the elements of the shape and create everything
320 // loop over slices
321 for (size_t i = 0; i < numSlices; ++i) {
322 const double z = (static_cast<double>(i) + 0.5) * sliceThickness - 0.5 * params.height;
323
324 // Number of elements in 1st annulus
325 // NOTE:
326 // For example, if the hollow cylinder consist of an inner cylinder surface with
327 // two annulus, and two annulus for the hollow ring (i.e. total four annulus for
328 // the outter cylinder surface). We have
329 // Ni = [6, 12, 18, 24]
330 // ^
331 // inner outter
332 const auto nSteps = params.innerRadius / deltaR;
333 size_t Ni = static_cast<size_t>(nSteps) * 6;
334 // loop over annuli
335 for (size_t j = 0; j < numAnnuli; ++j) {
336 Ni += 6;
337 const double R =
338 params.innerRadius +
339 (static_cast<double>(j) * (params.radius - params.innerRadius) / static_cast<double>(numAnnuli)) +
340 (0.5 * deltaR);
341
342 // all the volume elements in the ring/slice are the same
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);
347
348 // loop over elements in current annulus
349 for (size_t k = 0; k < Ni; ++k) {
350 const double phi = 2. * M_PI * static_cast<double>(k) / static_cast<double>(Ni);
351
352 const auto position = center + CalculatePosInCylinder(phi, R, z, coords);
353
354 assert(shape.isValid(position));
355
356 result.position.emplace_back(position);
357 result.volume.emplace_back(elementVolume);
358 // TODO should be customized for hollow cylinder
359 result.l1.emplace_back(calcDistanceInShapeNoCheck(beamDirection, shape, position));
360 } // loop over k
361 } // loop over j
362 } // loop over i
363
364 return result;
365}
366
367} // namespace Rasterize
368} // namespace Mantid::Geometry
double height
Definition: GetAllEi.cpp:155
double position
Definition: GetAllEi.cpp:154
#define fabs(x)
Definition: Matrix.cpp:22
double radius
Definition: Rasterize.cpp:31
V3D symmetryaxis
Definition: Rasterize.cpp:34
double innerRadius
Definition: Rasterize.cpp:39
V3D centerBottomBase
Definition: Rasterize.cpp:33
IObject : Interface for geometry objects.
Definition: IObject.h:41
virtual const detail::ShapeInfo & shapeInfo() const =0
virtual detail::ShapeInfo::GeometryShape shape() const =0
CylinderGeometry cylinderGeometry() const
Definition: ShapeInfo.cpp:53
HollowCylinderGeometry hollowCylinderGeometry() const
Definition: ShapeInfo.cpp:58
Class for 3D vectors.
Definition: V3D.h:34
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
Definition: V3D.h:278
MANTID_GEOMETRY_DLL Raster calculate(const Kernel::V3D &beamDirection, const IObject &shape, const double cubeSizeInMetre)
Definition: Rasterize.cpp:181
MANTID_GEOMETRY_DLL Raster calculateCylinder(const Kernel::V3D &beamDirection, const IObject &shape, const size_t numSlices, const size_t numAnnuli)
Definition: Rasterize.cpp:206
MANTID_GEOMETRY_DLL Raster calculateHollowCylinder(const Kernel::V3D &beamDirection, const IObject &shape, const size_t numSlices, const size_t numAnnuli)
Definition: Rasterize.cpp:280
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Definition: V3D.h:341
Holds the information used for doing numerical integrations of object in the beam.
Definition: Rasterize.h:21
std::vector< double > volume
Cached element volumes.
Definition: Rasterize.h:26
std::vector< double > l1
Cached L1 distances.
Definition: Rasterize.h:25
void reserve(size_t numVolumeElements)
Definition: Rasterize.cpp:18
std::vector< Kernel::V3D > position
Cached element positions.
Definition: Rasterize.h:27