Mantid
Loading...
Searching...
No Matches
MeshObject.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 +
16
17#include <algorithm>
18#include <memory>
19#include <utility>
20
21namespace Mantid::Geometry {
22
23MeshObject::MeshObject(std::vector<uint32_t> faces, std::vector<Kernel::V3D> vertices, const Kernel::Material &material)
24 : m_boundingBox(), m_id("MeshObject"), m_triangles(std::move(faces)), m_vertices(std::move(vertices)),
25 m_material(material) {
26
27 initialize();
28}
29
30MeshObject::MeshObject(std::vector<uint32_t> &&faces, std::vector<Kernel::V3D> &&vertices,
31 const Kernel::Material &&material)
32 : m_boundingBox(), m_id("MeshObject"), m_triangles(std::move(faces)), m_vertices(std::move(vertices)),
33 m_material(material) {
34
35 initialize();
36}
37
38// Do things that need to be done in constructor
40
42 m_handler = std::make_shared<GeometryHandler>(*this);
43}
44
49
54
61 // May enclose volume if there are at
62 // at least 4 triangles and 4 vertices (Tetrahedron)
63 return (numberOfTriangles() >= 4 && numberOfVertices() >= 4);
64}
65
71bool MeshObject::isValid(const Kernel::V3D &point) const {
72
74 if (!bb.isPointInside(point)) {
75 return false;
76 }
77
78 Kernel::V3D direction(0.0, 0.0, 1.0); // direction to look for intersections
79 std::vector<Kernel::V3D> intersectionPoints;
80 std::vector<TrackDirection> entryExitFlags;
81
82 getIntersections(point, direction, intersectionPoints, entryExitFlags);
83
84 if (intersectionPoints.empty()) {
85 return false;
86 }
87
88 // True if point is on surface
89 const auto it = std::find_if(
90 intersectionPoints.cbegin(), intersectionPoints.cend(),
91 [this, &point](const auto &intersectionPoint) { return point.distance(intersectionPoint) < M_TOLERANCE; });
92 if (it != intersectionPoints.cend())
93 return true;
94
95 // Look for nearest point then check its entry-exit flag
96 double nearestPointDistance = point.distance(intersectionPoints[0]);
97 size_t nearestPointIndex = 0;
98 for (size_t i = 1; i < intersectionPoints.size(); ++i) {
99 if (point.distance(intersectionPoints[i]) < nearestPointDistance) {
100 nearestPointDistance = point.distance(intersectionPoints[i]);
101 nearestPointIndex = i;
102 }
103 }
104 return (entryExitFlags[nearestPointIndex] == TrackDirection::LEAVING);
105}
106
112bool MeshObject::isOnSide(const Kernel::V3D &point) const {
113
115 if (!bb.isPointInside(point)) {
116 return false;
117 }
118
119 const std::vector<Kernel::V3D> directions = {Kernel::V3D{0, 0, 1}, Kernel::V3D{0, 1, 0},
120 Kernel::V3D{1, 0, 0}}; // directions to look for intersections
121 // We have to look in several directions in case a point is on a face
122 // or edge parallel to the first direction or also the second direction.
123 for (const auto &direction : directions) {
124 std::vector<Kernel::V3D> intersectionPoints;
125 std::vector<TrackDirection> entryExitFlags;
126
127 getIntersections(point, direction, intersectionPoints, entryExitFlags);
128
129 if (intersectionPoints.empty()) {
130 return false;
131 }
132
133 const auto it = std::find_if(
134 intersectionPoints.cbegin(), intersectionPoints.cend(),
135 [this, &point](const auto &intersectionPoint) { return point.distance(intersectionPoint) < M_TOLERANCE; });
136 if (it != intersectionPoints.cend())
137 return true;
138 }
139 return false;
140}
141
148 int originalCount = UT.count(); // Number of intersections original track
150 if (!bb.doesLineIntersect(UT)) {
151 return 0;
152 }
153
154 std::vector<Kernel::V3D> intersectionPoints;
155 std::vector<TrackDirection> entryExit;
156
157 getIntersections(UT.startPoint(), UT.direction(), intersectionPoints, entryExit);
158 if (intersectionPoints.empty())
159 return 0; // Quit if no intersections found
160
161 // For a 3D mesh, a ray may intersect several segments
162 for (size_t i = 0; i < intersectionPoints.size(); ++i) {
163 UT.addPoint(entryExit[i], intersectionPoints[i], *this);
164 }
165 UT.buildLink();
166
167 return UT.count() - originalCount;
168}
169
176double MeshObject::distance(const Track &track) const {
177 Kernel::V3D vertex1, vertex2, vertex3, intersection;
178 TrackDirection unused;
179 for (size_t i = 0; getTriangle(i, vertex1, vertex2, vertex3); ++i) {
180 if (MeshObjectCommon::rayIntersectsTriangle(track.startPoint(), track.direction(), vertex1, vertex2, vertex3,
181 intersection, unused)) {
182 return track.startPoint().distance(intersection);
183 }
184 }
185 std::ostringstream os;
186 os << "Unable to find intersection with object with track starting at " << track.startPoint() << " in direction "
187 << track.direction() << "\n";
188 throw std::runtime_error(os.str());
189}
190
198void MeshObject::getIntersections(const Kernel::V3D &start, const Kernel::V3D &direction,
199 std::vector<Kernel::V3D> &intersectionPoints,
200 std::vector<TrackDirection> &entryExitFlags) const {
201
202 Kernel::V3D vertex1, vertex2, vertex3, intersection;
203 TrackDirection entryExit;
204 for (size_t i = 0; getTriangle(i, vertex1, vertex2, vertex3); ++i) {
205 if (MeshObjectCommon::rayIntersectsTriangle(start, direction, vertex1, vertex2, vertex3, intersection, entryExit)) {
206 intersectionPoints.emplace_back(intersection);
207 entryExitFlags.emplace_back(entryExit);
208 }
209 }
210 // still need to deal with edge cases
211}
212
213/*
214 * Get a triangle - useful for iterating over triangles
215 * @param index :: Index of triangle in MeshObject
216 * @param v1 :: First vertex of triangle
217 * @param v2 :: Second vertex of triangle
218 * @param v3 :: Third vertex of triangle
219 * @returns true if the specified triangle exists
220 */
221bool MeshObject::getTriangle(const size_t index, Kernel::V3D &vertex1, Kernel::V3D &vertex2,
222 Kernel::V3D &vertex3) const {
223 bool triangleExists = index < m_triangles.size() / 3;
224 if (triangleExists) {
225 vertex1 = m_vertices[m_triangles[3 * index]];
226 vertex2 = m_vertices[m_triangles[3 * index + 1]];
227 vertex3 = m_vertices[m_triangles[3 * index + 2]];
228 }
229 return triangleExists;
230}
231
240int MeshObject::calcValidType(const Kernel::V3D &point, const Kernel::V3D &uVec) const {
241 const Kernel::V3D shift(uVec * Kernel::Tolerance * 25.0);
242 const Kernel::V3D testA(point - shift);
243 const Kernel::V3D testB(point + shift);
244 const int flagA = isValid(testA);
245 const int flagB = isValid(testB);
246 if (!(flagA ^ flagB))
247 return 0;
248 return (flagA) ? -1 : 1;
249}
250
263void MeshObject::getBoundingBox(double &xmax, double &ymax, double &zmax, double &xmin, double &ymin,
264 double &zmin) const {
265 return MeshObjectCommon::getBoundingBox(m_vertices, m_boundingBox, xmax, ymax, zmax, xmin, ymin, zmin);
266}
267
273double MeshObject::solidAngle(const Kernel::V3D &observer) const {
274 double solidAngleSum(0), solidAngleNegativeSum(0);
275 Kernel::V3D vertex1, vertex2, vertex3;
276 for (size_t i = 0; this->getTriangle(i, vertex1, vertex2, vertex3); ++i) {
277 double sa = MeshObjectCommon::getTriangleSolidAngle(vertex1, vertex2, vertex3, observer);
278 if (sa > 0.0) {
279 solidAngleSum += sa;
280 } else {
281 solidAngleNegativeSum += sa;
282 }
283 }
284 /*
285 Same implementation as CSGObject. Assumes a convex closed mesh with
286 solidAngleSum == -solidAngleNegativeSum
287
288 Average is used to bypass issues with winding order. Surface normal
289 affects magnitude of solid angle. See CSGObject.
290 */
291 return 0.5 * (solidAngleSum - solidAngleNegativeSum);
292}
293
301double MeshObject::solidAngle(const Kernel::V3D &observer, const Kernel::V3D &scaleFactor) const {
302 std::vector<Kernel::V3D> scaledVertices;
303 scaledVertices.reserve(m_vertices.size());
304 std::transform(m_vertices.cbegin(), m_vertices.cend(), std::back_inserter(scaledVertices),
305 [&scaleFactor](const auto &vertex) { return scaleFactor * vertex; });
306 MeshObject meshScaled(m_triangles, scaledVertices, m_material);
307 return meshScaled.solidAngle(observer);
308}
309
314double MeshObject::volume() const {
315 // Select centre of bounding box as centre point.
316 // For each triangle calculate the signed volume of
317 // the tetrahedron formed by the triangle and the
318 // centre point. Then add to total.
319
321 double cX = 0.5 * (bb.xMax() + bb.xMin());
322 double cY = 0.5 * (bb.yMax() + bb.yMin());
323 double cZ = 0.5 * (bb.zMax() + bb.zMin());
324 Kernel::V3D centre(cX, cY, cZ);
325
326 double volumeTimesSix(0.0);
327
328 Kernel::V3D vertex1, vertex2, vertex3;
329 for (size_t i = 0; getTriangle(i, vertex1, vertex2, vertex3); ++i) {
330 Kernel::V3D a = vertex1 - centre;
331 Kernel::V3D b = vertex2 - centre;
332 Kernel::V3D c = vertex3 - centre;
333 volumeTimesSix += a.scalar_prod(b.cross_prod(c));
334 }
335
336 return volumeTimesSix / 6.0;
337}
338
345}
346
353
354 Kernel::V3D testPt(0, 0, 0);
355 // Try centre of bounding box as initial guess, if we have one.
356 const BoundingBox &boundingBox = getBoundingBox();
357 if (boundingBox.isNonNull()) {
358 testPt = boundingBox.centrePoint();
359 if (searchForObject(testPt)) {
360 point = testPt;
361 return 1;
362 }
363 }
364
365 return 0;
366}
367
378 const size_t maxAttempts) const {
379 const auto &bbox = getBoundingBox();
380 if (bbox.isNull()) {
381 throw std::runtime_error("Object::generatePointInObject() - Invalid "
382 "bounding box. Cannot generate new point.");
383 }
384 return generatePointInObject(rng, bbox, maxAttempts);
385}
386
398 const BoundingBox &activeRegion,
399 const size_t maxAttempts) const {
400
401 const auto point = RandomPoint::bounded(*this, rng, activeRegion, maxAttempts);
402
403 return point;
404}
405
412 //
413 // Method - check if point in object, if not search directions along
414 // principle axes using interceptSurface
415 //
416 if (isValid(point))
417 return true;
418 for (const auto &dir : {Kernel::V3D(1., 0., 0.), Kernel::V3D(-1., 0., 0.), Kernel::V3D(0., 1., 0.),
419 Kernel::V3D(0., -1., 0.), Kernel::V3D(0., 0., 1.), Kernel::V3D(0., 0., -1.)}) {
420 Geometry::Track tr(point, dir);
421 if (this->interceptSurface(tr) > 0) {
422 point = tr.cbegin()->entryPoint;
423 return true;
424 }
425 }
426 return false;
427}
428
434void MeshObject::setGeometryHandler(const std::shared_ptr<GeometryHandler> &h) {
435 if (h == nullptr)
436 return;
437 m_handler = h;
438}
439
444void MeshObject::draw() const {
445 if (m_handler == nullptr)
446 return;
447 // Render the Object
448 m_handler->render();
449}
450
457 if (m_handler == nullptr)
458 return;
459 // Render the Object
460 m_handler->initialize();
461}
462
466std::shared_ptr<GeometryHandler> MeshObject::getGeometryHandler() const {
467 // Check if the geometry handler is upto dated with the cache, if not then
468 // cache it now.
469 return m_handler;
470}
471
476void MeshObject::rotate(const Kernel::Matrix<double> &rotationMatrix) {
477 std::for_each(m_vertices.begin(), m_vertices.end(),
478 [&rotationMatrix](auto &vertex) { vertex.rotate(rotationMatrix); });
479}
480
485void MeshObject::translate(const Kernel::V3D &translationVector) {
486 std::transform(m_vertices.cbegin(), m_vertices.cend(), m_vertices.begin(),
487 [&translationVector](const auto &vertex) { return vertex + translationVector; });
488}
489
494void MeshObject::scale(const double scaleFactor) {
495 std::transform(m_vertices.cbegin(), m_vertices.cend(), m_vertices.begin(),
496 [&scaleFactor](const auto &vertex) { return vertex * scaleFactor; });
497}
498
505 if ((matrix.numCols() != 4) || (matrix.numRows() != 4)) {
506 throw "Transformation matrix must be 4 x 4";
507 }
508
509 // create homogenous coordinates for the input vector with 4th element
510 // equal to 1 (position)
511 for (Kernel::V3D &vertex : m_vertices) {
512 std::vector<double> vertexin(4);
513 vertexin[0] = vertex.X();
514 vertexin[1] = vertex.Y();
515 vertexin[2] = vertex.Z();
516 vertexin[3] = 1;
517 std::vector<double> vertexout(4);
518 matrix.multiplyPoint(vertexin, vertexout);
519 Kernel::V3D newvertex(vertexout[0], vertexout[1], vertexout[2]);
520 vertex = newvertex;
521 }
522}
523
528 return; // Hopefully nothing necessary here
529}
530
534size_t MeshObject::numberOfTriangles() const { return m_triangles.size() / 3; }
535
539std::vector<uint32_t> MeshObject::getTriangles() const { return m_triangles; }
540
544size_t MeshObject::numberOfVertices() const { return m_vertices.size(); }
545
550
554const std::vector<Kernel::V3D> &MeshObject::getV3Ds() const { return m_vertices; }
555
557
559 throw std::runtime_error("MeshObject::shapeInfo() is not implemented");
560}
561
565void MeshObject::GetObjectGeom(detail::ShapeInfo::GeometryShape &type, std::vector<Kernel::V3D> &vectors,
566 double &innerRadius, double &radius, double &height) const {
567 // In practice, this outputs type = -1,
568 // to indicate not a "standard" object (cuboid/cone/cyl/sphere).
569 // Retained for possible future use.
571 if (m_handler == nullptr)
572 return;
573 m_handler->GetObjectGeom(type, vectors, innerRadius, radius, height);
574}
575
576} // namespace Mantid::Geometry
double height
Definition: GetAllEi.cpp:155
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
double radius
Definition: Rasterize.cpp:31
double innerRadius
Definition: Rasterize.cpp:39
A simple structure that defines an axis-aligned cuboid shaped bounding box for a geometrical object.
Definition: BoundingBox.h:34
bool isPointInside(const Kernel::V3D &point) const
Is the given point within the bounding box?
Definition: BoundingBox.cpp:28
double xMax() const
Return the maximum value of X.
Definition: BoundingBox.h:80
bool isNonNull() const
Is the box considered valid. Convenience for !isNull()
Definition: BoundingBox.h:106
double zMin() const
Return the minimum value of Z.
Definition: BoundingBox.h:86
double zMax() const
Return the maximum value of Z.
Definition: BoundingBox.h:88
double yMax() const
Return the maximum value of Y.
Definition: BoundingBox.h:84
double xMin() const
Return the minimum value of X.
Definition: BoundingBox.h:78
bool doesLineIntersect(const Track &track) const
Does a specified track intersect the bounding box.
Definition: BoundingBox.cpp:44
Kernel::V3D centrePoint() const
Returns the centre of the bounding box.
Definition: BoundingBox.h:94
double yMin() const
Return the minimum value of Y.
Definition: BoundingBox.h:82
Triangular Mesh Object.
Definition: MeshObject.h:50
Kernel::Material m_material
material composition
Definition: MeshObject.h:178
boost::optional< Kernel::V3D > generatePointInObject(Kernel::PseudoRandomNumberGenerator &rng, const size_t) const override
Select a random point within the object.
Definition: MeshObject.cpp:377
void translate(const Kernel::V3D &)
Translate the mesh according to the supplied x, y, z vector.
Definition: MeshObject.cpp:485
double solidAngle(const Kernel::V3D &observer) const override
Find solid angle of object wrt the observer.
Definition: MeshObject.cpp:273
int getPointInObject(Kernel::V3D &point) const override
Try to find a point that lies within (or on) the object.
Definition: MeshObject.cpp:352
size_t numberOfTriangles() const
Output functions for rendering, may also be used internally.
Definition: MeshObject.cpp:534
int interceptSurface(Geometry::Track &) const override
Given a track, fill the track with valid section.
Definition: MeshObject.cpp:147
const detail::ShapeInfo & shapeInfo() const override
Definition: MeshObject.cpp:558
std::vector< double > getVertices() const
get vertices
Definition: MeshObject.cpp:549
detail::ShapeInfo::GeometryShape shape() const override
Definition: MeshObject.cpp:556
std::vector< Kernel::V3D > m_vertices
Definition: MeshObject.h:176
bool hasValidShape() const override
Return whether this object has a valid shape.
Definition: MeshObject.cpp:60
bool isValid(const Kernel::V3D &) const override
Check if a point is inside.
Definition: MeshObject.cpp:71
void setMaterial(const Kernel::Material &material) override
Definition: MeshObject.cpp:53
void draw() const override
Draws the Object using geometry handler, If the handler is not set then this function does nothing.
Definition: MeshObject.cpp:444
MeshObject(std::vector< uint32_t > faces, std::vector< Kernel::V3D > vertices, const Kernel::Material &material)
Constructor.
Definition: MeshObject.cpp:23
const BoundingBox & getBoundingBox() const override
Return cached value of axis-aligned bounding box.
Definition: MeshObject.cpp:343
void getIntersections(const Kernel::V3D &start, const Kernel::V3D &direction, std::vector< Kernel::V3D > &intersectionPoints, std::vector< Mantid::Geometry::TrackDirection > &entryExitFlags) const
Get intersections.
Definition: MeshObject.cpp:198
std::vector< uint32_t > getTriangles() const
get faces
Definition: MeshObject.cpp:539
const std::vector< Kernel::V3D > & getV3Ds() const
get vertices in V3D form
Definition: MeshObject.cpp:554
double volume() const override
Calculates the volume of this object.
Definition: MeshObject.cpp:314
std::vector< uint32_t > m_triangles
Contents Triangles are specified by indices into a list of vertices.
Definition: MeshObject.h:175
void scale(const double scaleFactor)
Scale the mesh according to the supplied scale factor.
Definition: MeshObject.cpp:494
bool searchForObject(Kernel::V3D &point) const
Search object for valid point.
Definition: MeshObject.cpp:411
void initDraw() const override
Initializes/prepares the object to be rendered, this will generate geometry for object,...
Definition: MeshObject.cpp:456
size_t numberOfVertices() const
Read access to mesh object for rendering.
Definition: MeshObject.cpp:544
void updateGeometryHandler()
Updates the geometry handler if needed.
Definition: MeshObject.cpp:527
const Kernel::Material & material() const override
Definition: MeshObject.cpp:48
std::shared_ptr< GeometryHandler > m_handler
Geometry Handle for rendering.
Definition: MeshObject.h:160
std::shared_ptr< GeometryHandler > getGeometryHandler() const override
Returns the geometry handler.
Definition: MeshObject.cpp:466
void multiply(const Kernel::Matrix< double > &)
Transform the mesh (scale, translate, rotate) according to the supplied transformation matrix.
Definition: MeshObject.cpp:504
bool getTriangle(const size_t index, Kernel::V3D &v1, Kernel::V3D &v2, Kernel::V3D &v3) const
Get triangle.
Definition: MeshObject.cpp:221
void rotate(const Kernel::Matrix< double > &)
Rotate the mesh according to the supplied rotation matrix.
Definition: MeshObject.cpp:476
int calcValidType(const Kernel::V3D &Pt, const Kernel::V3D &uVec) const
Calculate if a point PT is a valid point on the track.
Definition: MeshObject.cpp:240
BoundingBox m_boundingBox
Cache for object's bounding box.
Definition: MeshObject.h:154
bool isOnSide(const Kernel::V3D &) const override
Determines wither point is on the surface.
Definition: MeshObject.cpp:112
void GetObjectGeom(detail::ShapeInfo::GeometryShape &type, std::vector< Kernel::V3D > &vectors, double &innerRadius, double &radius, double &height) const override
get info on standard shapes (none for Mesh Object)
Definition: MeshObject.cpp:565
void setGeometryHandler(const std::shared_ptr< GeometryHandler > &h)
Set Geometry Handler.
Definition: MeshObject.cpp:434
double distance(const Track &track) const override
Compute the distance to the first point of intersection with the surface.
Definition: MeshObject.cpp:176
Defines a track as a start point and a direction.
Definition: Track.h:165
const Kernel::V3D & startPoint() const
Returns the starting point.
Definition: Track.h:191
void addPoint(const TrackDirection direction, const Kernel::V3D &endPoint, const IObject &obj, const ComponentID compID=nullptr)
Adds a point of intersection to the track.
Definition: Track.cpp:124
void buildLink()
Construct links between added points.
Definition: Track.cpp:184
int count() const
Returns the number of links.
Definition: Track.h:219
const Kernel::V3D & direction() const
Returns the direction as a unit vector.
Definition: Track.h:193
LType::const_iterator cbegin() const
Returns an interator to the start of the set of links (const version)
Definition: Track.h:206
A material is defined as being composed of a given element, defined as a PhysicalConstants::NeutronAt...
Definition: Material.h:50
Numerical Matrix class.
Definition: Matrix.h:42
void multiplyPoint(const std::vector< T > &in, std::vector< T > &out) const
Multiply M*Vec.
Definition: Matrix.cpp:375
size_t numRows() const
Return the number of rows in the matrix.
Definition: Matrix.h:144
size_t numCols() const
Return the number of columns in the matrix.
Definition: Matrix.h:147
Defines a 1D pseudo-random number generator, i.e.
Class for 3D vectors.
Definition: V3D.h:34
double distance(const V3D &v) const noexcept
Calculates the distance between two vectors.
Definition: V3D.h:287
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
Definition: V3D.h:274
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
Definition: V3D.h:278
MANTID_GEOMETRY_DLL bool rayIntersectsTriangle(const Kernel::V3D &start, const Kernel::V3D &direction, const Kernel::V3D &v1, const Kernel::V3D &v2, const Kernel::V3D &v3, Kernel::V3D &intersection, TrackDirection &entryExit)
Get intersection points and their in out directions on the given ray.
MANTID_GEOMETRY_DLL void checkVertexLimit(size_t nVertices)
MANTID_GEOMETRY_DLL std::vector< double > getVertices(const std::vector< Kernel::V3D > &vertices)
getVertices converts vector Kernel::V3D to vector doubles.
MANTID_GEOMETRY_DLL double getTriangleSolidAngle(const Kernel::V3D &a, const Kernel::V3D &b, const Kernel::V3D &c, const Kernel::V3D &observer)
Find the solid angle of a triangle defined by vectors a,b,c from point "observer".
MANTID_GEOMETRY_DLL const BoundingBox & getBoundingBox(const std::vector< Kernel::V3D > &vertices, BoundingBox &cacheBB)
Takes input vertices and calculates bounding box.
boost::optional< Kernel::V3D > bounded(const detail::ShapeInfo &shapeInfo, Kernel::PseudoRandomNumberGenerator &rng, const BoundingBox &box, size_t maxAttempts)
Return a random point in a known shape restricted by a bounding box.
Definition: RandomPoint.h:62
bool MANTID_GEOMETRY_DLL intersection(const ConvexPolygon &P, const ConvexPolygon &Q, ConvexPolygon &out)
Compute the instersection of two convex polygons.
constexpr double Tolerance
Standard tolerance value.
Definition: Tolerance.h:12
STL namespace.