Mantid
Loading...
Searching...
No Matches
ShapeFactory.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 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
10
23
24#include "MantidKernel/Logger.h"
25#include "MantidKernel/Quat.h"
26
27#include <Poco/AutoPtr.h>
28#include <Poco/DOM/DOMParser.h>
29#include <Poco/DOM/DOMWriter.h>
30#include <Poco/DOM/Document.h>
31#include <Poco/DOM/Element.h>
32#include <Poco/DOM/NodeList.h>
33
34#include "boost/make_shared.hpp"
35
36#include <algorithm>
37
38using Poco::XML::Document;
39using Poco::XML::DOMParser;
40using Poco::XML::DOMWriter;
41using Poco::XML::Element;
42using Poco::XML::Node;
43using Poco::XML::NodeList;
44
45namespace Mantid::Geometry {
46
47using namespace Kernel;
48
49namespace {
50const V3D DEFAULT_CENTRE(0, 0, 0);
51const V3D DEFAULT_AXIS(0, 0, 1);
52
54Logger g_log("ShapeFactory");
55} // namespace
56
57namespace {
58std::vector<double> DegreesToRadians(const std::vector<double> &anglesDegrees) {
59 std::vector<double> anglesRadians;
60 std::transform(anglesDegrees.cbegin(), anglesDegrees.cend(), std::back_inserter(anglesRadians),
61 [](const auto angle) { return angle * M_PI / 180.0; });
62 return anglesRadians;
63}
64} // namespace
65
74std::shared_ptr<CSGObject> ShapeFactory::createShape(std::string shapeXML, bool addTypeTag) {
75 // wrap in a type tag
76 if (addTypeTag)
77 shapeXML = "<type name=\"userShape\"> " + shapeXML + " </type>";
78
79 // Set up the DOM parser and parse xml string
80 DOMParser pParser;
81 Poco::AutoPtr<Document> pDoc;
82 try {
83 pDoc = pParser.parseString(shapeXML);
84 } catch (...) {
85 g_log.warning("Unable to parse XML string " + shapeXML + " . Empty geometry Object is returned.");
86
87 return std::make_shared<CSGObject>();
88 }
89 // Get pointer to root element
90 Element *pRootElem = pDoc->documentElement();
91
92 // convert into a Geometry object
93 return createShape(pRootElem);
94}
95
104std::shared_ptr<CSGObject> ShapeFactory::createShape(Poco::XML::Element *pElem) {
105 // Write the definition to a string to store in the final object
106 std::stringstream xmlstream;
107 DOMWriter writer;
108 writer.writeNode(xmlstream, pElem);
109 std::string shapeXML = xmlstream.str();
110 auto retVal = std::make_shared<CSGObject>(shapeXML);
111
112 // if no <algebra> element then use default algebra
113 bool defaultAlgebra(false);
114 // get algebra string
115 Poco::AutoPtr<NodeList> pNL_algebra = pElem->getElementsByTagName("algebra");
116 std::string algebraFromUser;
117 if (pNL_algebra->length() == 0) {
118 defaultAlgebra = true;
119 } else if (pNL_algebra->length() == 1) {
120 auto *pElemAlgebra = static_cast<Element *>(pNL_algebra->item(0));
121 algebraFromUser = pElemAlgebra->getAttribute("val");
122 } else {
123 g_log.warning() << "More than one algebra string defined for this shape. "
124 << "Maximum one allowed. Therefore empty shape is returned.";
125 return retVal;
126 }
127
128 Poco::AutoPtr<NodeList> pNL_gonio = pElem->getElementsByTagName("goniometer");
129 auto *pElemGonio = static_cast<Element *>(pNL_gonio->item(0));
131 if (pElemGonio) {
132 // Parse the rotate matrix, defined in units of radians
133 for (size_t i = 0; i < 3; ++i) {
134 for (size_t j = 0; j < 3; ++j) {
135 m_gonioRotateMatrix[i][j] = getDoubleAttribute(pElemGonio, "a" + std::to_string(i + 1) + std::to_string(j + 1));
136 }
137 }
138 }
139
140 Poco::AutoPtr<NodeList> pNL_rotate_all = pElem->getElementsByTagName("rotate-all");
141 auto *pElemRotateAll = static_cast<Element *>(pNL_rotate_all->item(0));
143 if (pElemRotateAll) {
144 std::vector<double> rotateAngles = DegreesToRadians(parsePosition(pElemRotateAll));
145 m_rotateAllMatrix = generateMatrix(rotateAngles[0], rotateAngles[1], rotateAngles[2]);
146 }
147
148 // match id given to a shape by the user to
149 // id understandable by Mantid code
150 std::map<std::string, std::string> idMatching;
151
152 // loop over all the sub-elements of pElem
153 Poco::AutoPtr<NodeList> pNL = pElem->childNodes(); // get all child nodes
154 unsigned long pNL_length = pNL->length();
155 int numPrimitives = 0;
156 // stores the primitives that will be
157 // used to build final shape
158 std::map<int, std::shared_ptr<Surface>> primitives;
159 // used to build up unique id's for each shape added. Must start
160 // from int > zero.
161 int l_id = 1;
162 // Element of fixed complete object
163 Element *lastElement = nullptr;
164 for (unsigned int i = 0; i < pNL_length; i++) {
165 if ((pNL->item(i))->nodeType() == Node::ELEMENT_NODE) {
166 auto *pE = static_cast<Element *>(pNL->item(i));
167
168 // assume for now that if sub-element has attribute id then it is a shape
169 // element
170 if (pE->hasAttribute("id")) {
171 std::string idFromUser = pE->getAttribute("id"); // get id
172
173 std::string primitiveName = pE->tagName(); // get name of primitive
174
175 // if there are any error thrown while parsing the XML string for a
176 // given shape
177 // write out a warning to the user that this shape is ignored. If all
178 // shapes are ignored
179 // this way an empty object is returned to the user.
180 try {
181 if (primitiveName == "sphere") {
182 lastElement = pE;
183 idMatching[idFromUser] = parseSphere(pE, primitives, l_id);
184 numPrimitives++;
185 } else if (primitiveName == "infinite-plane") {
186 idMatching[idFromUser] = parseInfinitePlane(pE, primitives, l_id);
187 retVal->setFiniteGeometryFlag(false);
188 numPrimitives++;
189 } else if (primitiveName == "infinite-cylinder") {
190 idMatching[idFromUser] = parseInfiniteCylinder(pE, primitives, l_id);
191 retVal->setFiniteGeometryFlag(false);
192 numPrimitives++;
193 } else if (primitiveName == "cylinder") {
194 lastElement = pE;
195 idMatching[idFromUser] = parseCylinder(pE, primitives, l_id);
196 numPrimitives++;
197 } else if (primitiveName == "hollow-cylinder") {
198 lastElement = pE;
199 idMatching[idFromUser] = parseHollowCylinder(pE, primitives, l_id);
200 numPrimitives++;
201 } else if (primitiveName == "cuboid") {
202 lastElement = pE;
203 idMatching[idFromUser] = parseCuboid(pE, primitives, l_id);
204 numPrimitives++;
205 } else if (primitiveName == "infinite-cone") {
206 idMatching[idFromUser] = parseInfiniteCone(pE, primitives, l_id);
207 retVal->setFiniteGeometryFlag(false);
208 numPrimitives++;
209 } else if (primitiveName == "cone") {
210 lastElement = pE;
211 idMatching[idFromUser] = parseCone(pE, primitives, l_id);
212 numPrimitives++;
213 } else if (primitiveName == "hexahedron") {
214 lastElement = pE;
215 idMatching[idFromUser] = parseHexahedron(pE, primitives, l_id);
216 numPrimitives++;
217 } else if (primitiveName == "tapered-guide") {
218 idMatching[idFromUser] = parseTaperedGuide(pE, primitives, l_id);
219 numPrimitives++;
220 } else if (primitiveName == "torus") {
221 idMatching[idFromUser] = parseTorus(pE, primitives, l_id);
222 numPrimitives++;
223 } else if (primitiveName == "slice-of-cylinder-ring") {
224 idMatching[idFromUser] = parseSliceOfCylinderRing(pE, primitives, l_id);
225 numPrimitives++;
226 } else {
227 g_log.warning(primitiveName + " not a recognised geometric shape. This shape is ignored.");
228 }
229 } catch (std::invalid_argument &e) {
230 g_log.warning() << e.what() << " <" << primitiveName << "> shape is ignored.";
231 } catch (std::runtime_error &e) {
232 g_log.warning() << e.what() << " <" << primitiveName << "> shape is ignored.";
233 } catch (...) {
234 g_log.warning() << " Problem with parsing XML string for <" << primitiveName << ">. This shape is ignored.";
235 }
236 }
237 }
238 }
239
240 if (!defaultAlgebra) {
241 // Translate algebra string defined by the user into something Mantid can
242 // understand
243
244 // std::string algebra; // to hold algebra in a way Mantid can understand
245 std::map<std::string, std::string>::iterator iter;
246 std::map<size_t, std::string, std::greater<size_t>> allFound;
247 for (iter = idMatching.begin(); iter != idMatching.end(); ++iter) {
248 size_t found = algebraFromUser.find(iter->first);
249
250 if (found == std::string::npos) {
251 defaultAlgebra = true;
252 g_log.warning() << "Algebra shape Warning: " + iter->first +
253 " not found in algebra string: " + algebraFromUser + "\n" +
254 ". Default to equal shape to intersection of those defined.";
255 break;
256 } else {
257 allFound[found] = iter->first;
258 }
259 }
260
261 // Here do the actually swapping of strings
262 // but only if the algebra containes all the shapes
263 if (allFound.size() == idMatching.size()) {
264 std::map<size_t, std::string, std::greater<size_t>>::iterator iter2;
265 for (iter2 = allFound.begin(); iter2 != allFound.end(); ++iter2) {
266 // std::string kuse = iter2->second;
267 algebraFromUser.replace(iter2->first, (iter2->second).size(), idMatching[iter2->second]);
268 }
269 }
270 }
271 if (defaultAlgebra) {
272 algebraFromUser = ""; // reset in case we are overwriten invalid string
273 std::map<std::string, std::string>::iterator iter;
274 for (iter = idMatching.begin(); iter != idMatching.end(); ++iter) {
275 algebraFromUser.append(iter->second + " "); // default is intersection
276 }
277 }
278
279 // check to see if there actually were any primitives in 'type' xml element
280 // and if yes then return empty Object. Otherwise populate Object with the
281 // primitives
282
283 if (numPrimitives == 0)
284 return retVal;
285 else {
286 retVal->setObject(21, algebraFromUser);
287 retVal->populate(primitives);
288 // check whether there is only one surface/closed surface
289 if (numPrimitives == 1 && lastElement != nullptr) // special case
290 {
291 // parse the primitive and create a Geometry handler for the object
292 createGeometryHandler(lastElement, retVal);
293 }
294
295 // get bounding box string
296 Poco::AutoPtr<NodeList> pNL_boundingBox = pElem->getElementsByTagName("bounding-box");
297 if (pNL_boundingBox->length() != 1) // should probably throw an error if
298 // more than 1 bounding box is
299 // defined...
300 return retVal;
301
302 double xmin = std::stod((getShapeElement(pElem, "x-min"))->getAttribute("val"));
303 double ymin = std::stod((getShapeElement(pElem, "y-min"))->getAttribute("val"));
304 double zmin = std::stod((getShapeElement(pElem, "z-min"))->getAttribute("val"));
305 double xmax = std::stod((getShapeElement(pElem, "x-max"))->getAttribute("val"));
306 double ymax = std::stod((getShapeElement(pElem, "y-max"))->getAttribute("val"));
307 double zmax = std::stod((getShapeElement(pElem, "z-max"))->getAttribute("val"));
308
309 retVal->defineBoundingBox(xmax, ymax, zmax, xmin, ymin, zmin);
310
311 return retVal;
312 }
313}
314
326std::string ShapeFactory::parseSphere(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
327 int &l_id) {
328 Element *pElemCentre = getOptionalShapeElement(pElem, "centre");
329 Element *pElemRadius = getShapeElement(pElem, "radius");
330 const double radius = getDoubleAttribute(pElemRadius, "val");
331 V3D centre = pElemCentre ? parsePosition(pElemCentre) : DEFAULT_CENTRE;
332
333 // Only rotate the normal vector by the rotate-all and goniometer tags
334 // Special case: do not obey rotate tag as rotation around the centre of sphere is unneccessary
335 if (m_rotateAllMatrix != Kernel::Matrix<double>(3, 3, true)) {
336 // Apply automatic rotation due to the rotate-all tag
338 }
339 if (m_gonioRotateMatrix != Kernel::Matrix<double>(3, 3, true)) {
340 // Apply automatic rotation due to the goniometer that should NOT be manually defined by the user
342 }
343
344 prim[l_id] = std::make_shared<Sphere>(centre, radius);
345 const auto algebra = sphereAlgebra(l_id);
346 l_id++;
347 return algebra;
348}
349
355std::string ShapeFactory::sphereAlgebra(const int surfaceID) { return "(-" + std::to_string(surfaceID) + ")"; }
356
368std::string ShapeFactory::parseInfinitePlane(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
369 int &l_id) {
370 Element *pElemPip = getShapeElement(pElem, "point-in-plane");
371 Element *pElemNormal = getShapeElement(pElem, "normal-to-plane");
372 Element *pElem_rot = getOptionalShapeElement(pElem, "rotate");
373
374 V3D normVec = normalize(parsePosition(pElemNormal));
375 V3D centre = parsePosition(pElemPip);
376
377 // Rotate the normal vector by the rotate, rotate-all and goniometer tags
378 // Rotate the centre by the rotate-all and goniometer tags
379 if (pElem_rot) {
380 // Apply manual rotation supplied to rotate tag
381 std::vector<double> rotateAngles = DegreesToRadians(parsePosition(pElem_rot));
382 const std::vector<double> rotateMatrix = generateMatrix(rotateAngles[0], rotateAngles[1], rotateAngles[2]);
383 normVec.rotate(rotateMatrix);
384 }
385 if (m_rotateAllMatrix != Kernel::Matrix<double>(3, 3, true)) {
386 // Apply automatic rotation due to the rotate-all tag
388 normVec.rotate(m_rotateAllMatrix);
389 }
390 if (m_gonioRotateMatrix != Kernel::Matrix<double>(3, 3, true)) {
391 // Apply automatic rotation due to the goniometer that should NOT be manually defined by the user
394 }
395
396 // create infinite-plane
397 auto pPlane = std::make_shared<Plane>();
398 pPlane->setPlane(centre, normVec);
399 prim[l_id] = pPlane;
400
401 std::stringstream retAlgebraMatch;
402 retAlgebraMatch << "(" << l_id << ")";
403 l_id++;
404 return retAlgebraMatch.str();
405}
406
418std::string ShapeFactory::parseInfiniteCylinder(Poco::XML::Element *pElem,
419 std::map<int, std::shared_ptr<Surface>> &prim, int &l_id) {
420 Element *pElemCentre = getShapeElement(pElem, "centre");
421 Element *pElemAxis = getShapeElement(pElem, "axis");
422 Element *pElemRadius = getShapeElement(pElem, "radius");
423 Element *pElem_rot = getOptionalShapeElement(pElem, "rotate");
424
425 // getDoubleAttribute can throw - put the calls above any new
426 const double radius = getDoubleAttribute(pElemRadius, "val");
427 V3D normVec = normalize(parsePosition(pElemAxis));
428 V3D centre = parsePosition(pElemCentre);
429
430 // Rotate the normal vector by the rotate, rotate-all and goniometer tags
431 // Rotate the centre by the rotate-all and goniometer tags
432 if (pElem_rot) {
433 // Apply manual rotation supplied to rotate tag
434 std::vector<double> rotateAngles = DegreesToRadians(parsePosition(pElem_rot));
435 const std::vector<double> rotateMatrix = generateMatrix(rotateAngles[0], rotateAngles[1], rotateAngles[2]);
436 normVec.rotate(rotateMatrix);
437 }
438 if (m_rotateAllMatrix != Kernel::Matrix<double>(3, 3, true)) {
439 // Apply automatic rotation due to the rotate-all tag
441 normVec.rotate(m_rotateAllMatrix);
442 }
443 if (m_gonioRotateMatrix != Kernel::Matrix<double>(3, 3, true)) {
444 // Apply automatic rotation due to the goniometer that should NOT be manually defined by the user
447 }
448
449 // create infinite-cylinder
450 auto pCylinder = std::make_shared<Cylinder>();
451 pCylinder->setNorm(normVec);
452 pCylinder->setCentre(centre);
453
454 pCylinder->setRadius(radius);
455 prim[l_id] = pCylinder;
456
457 std::stringstream retAlgebraMatch;
458 retAlgebraMatch << "(-" << l_id << ")";
459 l_id++;
460 return retAlgebraMatch.str();
461}
462
474std::string ShapeFactory::parseCylinder(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
475 int &l_id) {
476 Element *pElemBase = getShapeElement(pElem, "centre-of-bottom-base");
477 Element *pElemAxis = getShapeElement(pElem, "axis");
478 Element *pElemRadius = getShapeElement(pElem, "radius");
479 Element *pElemHeight = getShapeElement(pElem, "height");
480 Element *pElem_rot = getOptionalShapeElement(pElem, "rotate");
481
482 V3D normVec = normalize(parsePosition(pElemAxis));
483
484 // getDoubleAttribute can throw - put the calls above any new
485 const double radius = getDoubleAttribute(pElemRadius, "val");
486 const double height = getDoubleAttribute(pElemHeight, "val");
487
488 // add infinite cylinder
489 auto pCylinder = std::make_shared<Cylinder>();
490 V3D centreOfBottomBase = parsePosition(pElemBase);
491 V3D centre = centreOfBottomBase + normVec * (0.5 * height);
492 pCylinder->setRadius(radius);
493
494 // Rotate the normal vector by the rotate, rotate-all and goniometer tags
495 // Rotate the centre by the rotate-all and goniometer tags
496 if (pElem_rot) {
497 // Apply manual rotation supplied to rotate tag
498 std::vector<double> rotateAngles = DegreesToRadians(parsePosition(pElem_rot));
499 const std::vector<double> rotateMatrix = generateMatrix(rotateAngles[0], rotateAngles[1], rotateAngles[2]);
500 normVec.rotate(rotateMatrix);
501 }
502 if (m_rotateAllMatrix != Kernel::Matrix<double>(3, 3, true)) {
503 // Apply automatic rotation due to the rotate-all tag
505 normVec.rotate(m_rotateAllMatrix);
506 }
507 if (m_gonioRotateMatrix != Kernel::Matrix<double>(3, 3, true)) {
508 // Apply automatic rotation due to the goniometer that should NOT be manually defined by the user
511 }
512
513 pCylinder->setNorm(normVec);
514 pCylinder->setCentre(centre);
515
516 prim[l_id] = pCylinder;
517 std::stringstream retAlgebraMatch;
518 retAlgebraMatch << "(-" << l_id << " ";
519 l_id++;
520
521 // add top plane
522 auto pPlaneTop = std::make_shared<Plane>();
523 // to get point in top plane
524 V3D pointInPlaneTop = centre + (normVec * height * 0.5);
525 pPlaneTop->setPlane(pointInPlaneTop, normVec);
526 prim[l_id] = pPlaneTop;
527 retAlgebraMatch << "-" << l_id << " ";
528 l_id++;
529
530 // add bottom plane
531 auto pPlaneBottom = std::make_shared<Plane>();
532 V3D pointInPlaneBottom = centre - (normVec * height * 0.5);
533 pPlaneBottom->setPlane(pointInPlaneBottom, normVec);
534 prim[l_id] = pPlaneBottom;
535 retAlgebraMatch << "" << l_id << ")";
536 l_id++;
537
538 return retAlgebraMatch.str();
539}
540
553std::string ShapeFactory::parseHollowCylinder(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
554 int &l_id) {
555 Element *pElemBase = getShapeElement(pElem, "centre-of-bottom-base");
556 Element *pElemAxis = getShapeElement(pElem, "axis");
557 Element *pElemInnerRadius = getShapeElement(pElem, "inner-radius");
558 Element *pElemOuterRadius = getShapeElement(pElem, "outer-radius");
559 Element *pElemHeight = getShapeElement(pElem, "height");
560 Element *pElem_rot = getOptionalShapeElement(pElem, "rotate");
561
562 const double innerRadius = getDoubleAttribute(pElemInnerRadius, "val");
563 if (innerRadius <= 0.0) {
564 throw std::runtime_error("ShapeFactory::parseHollowCylinder(): inner-radius < 0.0");
565 }
566 const double outerRadius = getDoubleAttribute(pElemOuterRadius, "val");
567 if (outerRadius <= 0.0) {
568 throw std::runtime_error("ShapeFactory::parseHollowCylinder(): outer-radius < 0.0");
569 }
570 if (innerRadius > outerRadius) {
571 throw std::runtime_error("ShapeFactory::parseHollowCylinder(): inner-radius > outer-radius.");
572 }
573 const double height = getDoubleAttribute(pElemHeight, "val");
574 if (height <= 0.0) {
575 throw std::runtime_error("ShapeFactory::parseHollowCylinder(): height < 0.0");
576 }
577 V3D centreOfBottomBase = parsePosition(pElemBase);
578 V3D normVec = normalize(parsePosition(pElemAxis));
579 V3D centre = centreOfBottomBase + normVec * (0.5 * height);
580
581 // Rotate the normal vector by the rotate, rotate-all and goniometer tags
582 // Rotate the centre by the rotate-all and goniometer tags
583 if (pElem_rot) {
584 // Apply manual rotation supplied to rotate tag
585 std::vector<double> rotateAngles = DegreesToRadians(parsePosition(pElem_rot));
586 const std::vector<double> rotateMatrix = generateMatrix(rotateAngles[0], rotateAngles[1], rotateAngles[2]);
587 normVec.rotate(rotateMatrix);
588 }
589 if (m_rotateAllMatrix != Kernel::Matrix<double>(3, 3, true)) {
590 // Apply automatic rotation due to the rotate-all tag
592 normVec.rotate(m_rotateAllMatrix);
593 }
594 if (m_gonioRotateMatrix != Kernel::Matrix<double>(3, 3, true)) {
595 // Apply automatic rotation due to the goniometer that should NOT be manually defined by the user
598 }
599
600 // add outer infinite cylinder surface
601 auto outerCylinder = std::make_shared<Cylinder>();
602 outerCylinder->setCentre(centre);
603 outerCylinder->setNorm(normVec);
604 outerCylinder->setRadius(outerRadius);
605 prim[l_id] = outerCylinder;
606
607 std::stringstream retAlgebraMatch;
608 retAlgebraMatch << "(-" << l_id << " ";
609 l_id++;
610
611 // add inner infinite cylinder surface
612 auto innerCylinder = std::make_shared<Cylinder>();
613 innerCylinder->setCentre(centre);
614 innerCylinder->setNorm(normVec);
615 innerCylinder->setRadius(innerRadius);
616 prim[l_id] = innerCylinder;
617 retAlgebraMatch << l_id << " ";
618 l_id++;
619
620 // add top plane
621 auto pPlaneTop = std::make_shared<Plane>();
622 // to get point in top plane
623 V3D pointInPlaneTop = centre + (normVec * height * 0.5);
624 pPlaneTop->setPlane(pointInPlaneTop, normVec);
625 prim[l_id] = pPlaneTop;
626 retAlgebraMatch << "-" << l_id << " ";
627 l_id++;
628
629 // add bottom plane
630 auto pPlaneBottom = std::make_shared<Plane>();
631 V3D pointInPlaneBottom = centre - (normVec * height * 0.5);
632 pPlaneBottom->setPlane(pointInPlaneBottom, normVec);
633 prim[l_id] = pPlaneBottom;
634 retAlgebraMatch << "" << l_id << ")";
635 l_id++;
636
637 return retAlgebraMatch.str();
638}
639
649CuboidCorners ShapeFactory::parseCuboid(Poco::XML::Element *pElem) {
650 // Users have two syntax options when defining cuboids:
651
652 // A - "Point" syntax.
653 Element *pElem_lfb = getOptionalShapeElement(pElem, "left-front-bottom-point");
654 Element *pElem_lft = getOptionalShapeElement(pElem, "left-front-top-point");
655 Element *pElem_lbb = getOptionalShapeElement(pElem, "left-back-bottom-point");
656 Element *pElem_rfb = getOptionalShapeElement(pElem, "right-front-bottom-point");
657
658 // B - "Alternate" syntax.
659 Element *pElem_height = getOptionalShapeElement(pElem, "height");
660 Element *pElem_width = getOptionalShapeElement(pElem, "width");
661 Element *pElem_depth = getOptionalShapeElement(pElem, "depth");
662 Element *pElem_centre = getOptionalShapeElement(pElem, "centre");
663 Element *pElem_axis = getOptionalShapeElement(pElem, "axis");
664 Element *pElem_rot = getOptionalShapeElement(pElem, "rotate");
665
666 const bool usingPointSyntax = pElem_lfb && pElem_lft && pElem_lbb && pElem_rfb;
667 const bool usingAlternateSyntax = pElem_height && pElem_width && pElem_depth;
668
669 const bool usedPointSyntaxField = pElem_lfb || pElem_lft || pElem_lbb || pElem_rfb;
670 const bool usedAlternateSyntaxField =
671 pElem_height || pElem_width || pElem_depth || pElem_centre || pElem_axis || pElem_rot;
672
673 const std::string SYNTAX_ERROR_MSG = "XML element: <" + pElem->tagName() +
674 "> may contain EITHER corner points (LFB, LFT, LBB and RFB) OR " +
675 "height, width, depth, centre and axis values.";
676
677 CuboidCorners result;
678
679 if (usingPointSyntax && !usingAlternateSyntax) {
680 if (usedAlternateSyntaxField)
681 throw std::invalid_argument(SYNTAX_ERROR_MSG);
682
683 result.lfb = parsePosition(pElem_lfb);
684 result.lft = parsePosition(pElem_lft);
685 result.lbb = parsePosition(pElem_lbb);
686 result.rfb = parsePosition(pElem_rfb);
687 } else if (usingAlternateSyntax && !usingPointSyntax) {
688 if (usedPointSyntaxField)
689 throw std::invalid_argument(SYNTAX_ERROR_MSG);
690
691 const double deltaH = getDoubleAttribute(pElem_height, "val") / 2;
692 const double deltaW = getDoubleAttribute(pElem_width, "val") / 2;
693 const double deltaD = getDoubleAttribute(pElem_depth, "val") / 2;
694
695 V3D centre = pElem_centre ? parsePosition(pElem_centre) : DEFAULT_CENTRE;
696
697 result.lfb = V3D(-deltaW, -deltaH, -deltaD);
698 result.lft = V3D(-deltaW, deltaH, -deltaD);
699 result.lbb = V3D(-deltaW, -deltaH, deltaD);
700 result.rfb = V3D(deltaW, -deltaH, -deltaD);
701
702 if (pElem_axis) {
703 // Use a quarternion to do a rotate for us, with respect to the default
704 // axis. Our "Quat" implementation requires that the vectors passed to
705 // it be normalised.
706 const V3D axis = normalize(parsePosition(pElem_axis));
707 const Quat rotate(DEFAULT_AXIS, axis);
708
709 rotate.rotate(result.lfb);
710 rotate.rotate(result.lft);
711 rotate.rotate(result.lbb);
712 rotate.rotate(result.rfb);
713 }
714
715 // Rotate the points by the rotate, rotate-all and goniometer tags
716 // Rotate the centre by the rotate-all and goniometer tags
717 if (pElem_rot) {
718 // Apply manual rotation supplied to rotate tag
719 std::vector<double> rotateAngles = DegreesToRadians(parsePosition(pElem_rot));
720 const std::vector<double> rotateMatrix = generateMatrix(rotateAngles[0], rotateAngles[1], rotateAngles[2]);
721 result.rotatePoints(rotateMatrix);
722 }
723 if (m_rotateAllMatrix != Kernel::Matrix<double>(3, 3, true)) {
724 // Apply automatic rotation due to the rotate-all tag
727 }
728 if (m_gonioRotateMatrix != Kernel::Matrix<double>(3, 3, true)) {
729 // Apply automatic rotation due to the goniometer that should NOT be manually defined by the user
732 }
733
734 result.lfb += centre;
735 result.lft += centre;
736 result.lbb += centre;
737 result.rfb += centre;
738 } else {
739 throw std::invalid_argument(SYNTAX_ERROR_MSG);
740 }
741
742 return result;
743}
744
756std::string ShapeFactory::parseCuboid(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
757 int &l_id) {
758 auto corners = parseCuboid(pElem);
759
760 const V3D pointTowardBack = normalize(corners.lbb - corners.lfb);
761
762 // add front plane cutoff
763 auto pPlaneFrontCutoff = std::make_shared<Plane>();
764 try {
765 pPlaneFrontCutoff->setPlane(corners.lfb, pointTowardBack);
766 } catch (std::invalid_argument &) {
767 throw;
768 }
769 prim[l_id] = pPlaneFrontCutoff;
770
771 std::stringstream retAlgebraMatch;
772 retAlgebraMatch << "(" << l_id << " ";
773 l_id++;
774
775 // add back plane cutoff
776 auto pPlaneBackCutoff = std::make_shared<Plane>();
777 try {
778 pPlaneBackCutoff->setPlane(corners.lbb, pointTowardBack);
779 } catch (std::invalid_argument &) {
780 throw;
781 }
782 prim[l_id] = pPlaneBackCutoff;
783 retAlgebraMatch << "-" << l_id << " ";
784 l_id++;
785
786 const V3D pointTowardRight = normalize(corners.rfb - corners.lfb);
787
788 // add left plane cutoff
789 auto pPlaneLeftCutoff = std::make_shared<Plane>();
790 try {
791 pPlaneLeftCutoff->setPlane(corners.lfb, pointTowardRight);
792 } catch (std::invalid_argument &) {
793 throw;
794 }
795 prim[l_id] = pPlaneLeftCutoff;
796 retAlgebraMatch << "" << l_id << " ";
797 l_id++;
798
799 // add right plane cutoff
800 auto pPlaneRightCutoff = std::make_shared<Plane>();
801 try {
802 pPlaneRightCutoff->setPlane(corners.rfb, pointTowardRight);
803 } catch (std::invalid_argument &) {
804 throw;
805 }
806 prim[l_id] = pPlaneRightCutoff;
807 retAlgebraMatch << "-" << l_id << " ";
808 l_id++;
809
810 const V3D pointTowardTop = normalize(corners.lft - corners.lfb);
811
812 // add bottom plane cutoff
813 auto pPlaneBottomCutoff = std::make_shared<Plane>();
814 try {
815 pPlaneBottomCutoff->setPlane(corners.lfb, pointTowardTop);
816 } catch (std::invalid_argument &) {
817 throw;
818 }
819 prim[l_id] = pPlaneBottomCutoff;
820 retAlgebraMatch << "" << l_id << " ";
821 l_id++;
822
823 // add top plane cutoff
824 auto pPlaneTopCutoff = std::make_shared<Plane>();
825 try {
826 pPlaneTopCutoff->setPlane(corners.lft, pointTowardTop);
827 } catch (std::invalid_argument &) {
828 throw;
829 }
830 prim[l_id] = pPlaneTopCutoff;
831 retAlgebraMatch << "-" << l_id << ")";
832 l_id++;
833
834 return retAlgebraMatch.str();
835}
836
848std::string ShapeFactory::parseInfiniteCone(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
849 int &l_id) {
850 Element *pElemTipPoint = getShapeElement(pElem, "tip-point");
851 Element *pElemAxis = getShapeElement(pElem, "axis");
852 Element *pElemAngle = getShapeElement(pElem, "angle");
853
854 const V3D normVec = normalize(parsePosition(pElemAxis));
855
856 // getDoubleAttribute can throw - put the calls above any new
857 const double angle = getDoubleAttribute(pElemAngle, "val");
858
859 // add infinite double cone
860 auto pCone = std::make_shared<Cone>();
861 pCone->setCentre(parsePosition(pElemTipPoint));
862 pCone->setNorm(normVec);
863 pCone->setAngle(angle);
864 prim[l_id] = pCone;
865
866 std::stringstream retAlgebraMatch;
867 retAlgebraMatch << "(" << l_id << " ";
868 l_id++;
869
870 // plane top cut of top part of double cone
871 auto pPlaneBottom = std::make_shared<Plane>();
872 pPlaneBottom->setPlane(parsePosition(pElemTipPoint), normVec);
873 prim[l_id] = pPlaneBottom;
874 retAlgebraMatch << "-" << l_id << ")";
875 l_id++;
876
877 return retAlgebraMatch.str();
878}
879
891std::string ShapeFactory::parseCone(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
892 int &l_id) {
893 Element *pElemTipPoint = getShapeElement(pElem, "tip-point");
894 Element *pElemAxis = getShapeElement(pElem, "axis");
895 Element *pElemAngle = getShapeElement(pElem, "angle");
896 Element *pElemHeight = getShapeElement(pElem, "height");
897
898 const V3D normVec = normalize(parsePosition(pElemAxis));
899
900 // getDoubleAttribute can throw - put the calls above any new
901 const double angle = getDoubleAttribute(pElemAngle, "val");
902 const double height = getDoubleAttribute(pElemHeight, "val");
903
904 // add infinite double cone
905 auto pCone = std::make_shared<Cone>();
906 pCone->setCentre(parsePosition(pElemTipPoint));
907 pCone->setNorm(normVec);
908 pCone->setAngle(angle);
909 prim[l_id] = pCone;
910
911 std::stringstream retAlgebraMatch;
912 retAlgebraMatch << "(" << l_id << " ";
913 l_id++;
914
915 // Plane to cut off cone from below
916 auto pPlaneTop = std::make_shared<Plane>();
917 V3D pointInPlane = parsePosition(pElemTipPoint);
918 pointInPlane -= (normVec * height);
919 pPlaneTop->setPlane(pointInPlane, normVec);
920 prim[l_id] = pPlaneTop;
921 retAlgebraMatch << "" << l_id << " ";
922 l_id++;
923
924 // plane top cut of top part of double cone
925 auto pPlaneBottom = std::make_shared<Plane>();
926 pPlaneBottom->setPlane(parsePosition(pElemTipPoint), normVec);
927 prim[l_id] = pPlaneBottom;
928 retAlgebraMatch << "-" << l_id << ")";
929 l_id++;
930
931 return retAlgebraMatch.str();
932}
933
941 std::map<int, std::shared_ptr<Surface>> &prim, int &l_id) {
942 // add front face
943 auto pPlaneFrontCutoff = std::make_shared<Plane>();
944 auto normal = (hex.rfb - hex.lfb).cross_prod(hex.lft - hex.lfb);
945
946 // V3D jjj = (normal*(rfb-rbb));
947 if (normal.scalar_prod(hex.rfb - hex.rbb) < 0)
948 normal *= -1.0;
949 pPlaneFrontCutoff->setPlane(hex.lfb, normal);
950 prim[l_id] = pPlaneFrontCutoff;
951 std::stringstream retAlgebraMatch;
952 retAlgebraMatch << "(-" << l_id << " ";
953 l_id++;
954
955 // add back face
956 auto pPlaneBackCutoff = std::make_shared<Plane>();
957 normal = (hex.rbb - hex.lbb).cross_prod(hex.lbt - hex.lbb);
958 if (normal.scalar_prod(hex.rfb - hex.rbb) < 0)
959 normal *= -1.0;
960 pPlaneBackCutoff->setPlane(hex.lbb, normal);
961 prim[l_id] = pPlaneBackCutoff;
962 retAlgebraMatch << "" << l_id << " ";
963 l_id++;
964
965 // add left face
966 auto pPlaneLeftCutoff = std::make_shared<Plane>();
967 normal = (hex.lbb - hex.lfb).cross_prod(hex.lft - hex.lfb);
968 if (normal.scalar_prod(hex.rfb - hex.lfb) < 0)
969 normal *= -1.0;
970 pPlaneLeftCutoff->setPlane(hex.lfb, normal);
971 prim[l_id] = pPlaneLeftCutoff;
972 retAlgebraMatch << "" << l_id << " ";
973 l_id++;
974
975 // add right face
976 auto pPlaneRightCutoff = std::make_shared<Plane>();
977 normal = (hex.rbb - hex.rfb).cross_prod(hex.rft - hex.rfb);
978 if (normal.scalar_prod(hex.rfb - hex.lfb) < 0)
979 normal *= -1.0;
980 pPlaneRightCutoff->setPlane(hex.rfb, normal);
981 prim[l_id] = pPlaneRightCutoff;
982 retAlgebraMatch << "-" << l_id << " ";
983 l_id++;
984
985 // add top face
986 auto pPlaneTopCutoff = std::make_shared<Plane>();
987 normal = (hex.rft - hex.lft).cross_prod(hex.lbt - hex.lft);
988 if (normal.scalar_prod(hex.rft - hex.rfb) < 0)
989 normal *= -1.0;
990 pPlaneTopCutoff->setPlane(hex.lft, normal);
991 prim[l_id] = pPlaneTopCutoff;
992 retAlgebraMatch << "-" << l_id << " ";
993 l_id++;
994
995 // add bottom face
996 auto pPlaneBottomCutoff = std::make_shared<Plane>();
997 normal = (hex.rfb - hex.lfb).cross_prod(hex.lbb - hex.lfb);
998 if (normal.scalar_prod(hex.rft - hex.rfb) < 0)
999 normal *= -1.0;
1000 pPlaneBottomCutoff->setPlane(hex.lfb, normal);
1001 prim[l_id] = pPlaneBottomCutoff;
1002 retAlgebraMatch << "" << l_id << ")";
1003 l_id++;
1004
1005 return retAlgebraMatch.str();
1006}
1007
1016Hexahedron ShapeFactory::parseHexahedron(Poco::XML::Element *pElem) {
1017 Element *pElem_lfb = getShapeElement(pElem, "left-front-bottom-point");
1018 Element *pElem_lft = getShapeElement(pElem, "left-front-top-point");
1019 Element *pElem_lbb = getShapeElement(pElem, "left-back-bottom-point");
1020 Element *pElem_lbt = getShapeElement(pElem, "left-back-top-point");
1021 Element *pElem_rfb = getShapeElement(pElem, "right-front-bottom-point");
1022 Element *pElem_rft = getShapeElement(pElem, "right-front-top-point");
1023 Element *pElem_rbb = getShapeElement(pElem, "right-back-bottom-point");
1024 Element *pElem_rbt = getShapeElement(pElem, "right-back-top-point");
1025 Element *pElem_rot = getOptionalShapeElement(pElem, "rotate");
1026
1027 const bool isValid =
1028 pElem_lfb && pElem_lft && pElem_lbb && pElem_lbt && pElem_rfb && pElem_rft && pElem_rbb && pElem_rbt;
1029
1030 std::ostringstream ERROR_MSG;
1031 ERROR_MSG << "XML element: <" + pElem->tagName() + ""
1032 << "> contains invalid syntax for defining hexahedron. The "
1033 "following points have not been defined:\n\n";
1034
1035 if (!pElem_lfb)
1036 ERROR_MSG << "left-front-bottom-point\n";
1037 if (!pElem_lft)
1038 ERROR_MSG << "left-front-top-point\n";
1039 if (!pElem_lbb)
1040 ERROR_MSG << "left-back-bottom-point\n";
1041 if (!pElem_lbt)
1042 ERROR_MSG << "left-back-top-point\n";
1043 if (!pElem_rfb)
1044 ERROR_MSG << "right-front-bottom-point\n";
1045 if (!pElem_rft)
1046 ERROR_MSG << "right-front-top-point\n";
1047 if (!pElem_rbb)
1048 ERROR_MSG << "right-back-bottom-point\n";
1049 if (!pElem_rbt)
1050 ERROR_MSG << "right-back-top-point\n";
1051
1052 if (!isValid)
1053 throw std::invalid_argument(ERROR_MSG.str());
1054
1055 Hexahedron hex;
1056 hex.lfb = parsePosition(pElem_lfb);
1057 hex.lft = parsePosition(pElem_lft);
1058 hex.lbb = parsePosition(pElem_lbb);
1059 hex.lbt = parsePosition(pElem_lbt);
1060 hex.rfb = parsePosition(pElem_rfb);
1061 hex.rft = parsePosition(pElem_rft);
1062 hex.rbb = parsePosition(pElem_rbb);
1063 hex.rbt = parsePosition(pElem_rbt);
1064
1065 // Rotate the points by the rotate, rotate-all and goniometer tags
1066 if (pElem_rot) {
1067 // Apply manual rotation supplied to rotate tag
1068 std::vector<double> rotateAngles = DegreesToRadians(parsePosition(pElem_rot));
1069 const std::vector<double> rotateMatrix = generateMatrix(rotateAngles[0], rotateAngles[1], rotateAngles[2]);
1070 hex.rotatePoints(rotateMatrix);
1071 }
1072 if (m_rotateAllMatrix != Kernel::Matrix<double>(3, 3, true)) {
1073 // Apply automatic rotation due to the rotate-all tag
1075 }
1076 if (m_gonioRotateMatrix != Kernel::Matrix<double>(3, 3, true)) {
1077 // Apply automatic rotation due to the goniometer that should NOT be manually defined by the user
1079 }
1080
1081 return hex;
1082}
1083
1095std::string ShapeFactory::parseHexahedron(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
1096 int &l_id) {
1097 Hexahedron hex = parseHexahedron(pElem);
1098
1099 return parseHexahedronFromStruct(hex, prim, l_id);
1100}
1101
1114std::string ShapeFactory::parseTaperedGuide(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
1115 int &l_id) {
1116 Element *pElemApertureStart = getShapeElement(pElem, "aperture-start");
1117 Element *pElemLength = getShapeElement(pElem, "length");
1118 Element *pElemApertureEnd = getShapeElement(pElem, "aperture-end");
1119 Element *pElemCentre = getOptionalShapeElement(pElem, "centre");
1120 Element *pElemAxis = getOptionalShapeElement(pElem, "axis");
1121 Element *pElem_rot = getOptionalShapeElement(pElem, "rotate");
1122
1123 // For centre and axis we allow defaults.
1124 V3D centre = pElemCentre ? parsePosition(pElemCentre) : DEFAULT_CENTRE;
1125 // Quat requires normalised axes.
1126 const V3D axis = normalize(pElemAxis ? parsePosition(pElemAxis) : DEFAULT_AXIS);
1127
1128 const double apertureStartWidth = getDoubleAttribute(pElemApertureStart, "width");
1129 const double apertureStartHeight = getDoubleAttribute(pElemApertureStart, "height");
1130 const double length = getDoubleAttribute(pElemLength, "val");
1131 const double apertureEndWidth = getDoubleAttribute(pElemApertureEnd, "width");
1132 const double apertureEndHeight = getDoubleAttribute(pElemApertureEnd, "height");
1133
1134 const double halfSW = apertureStartWidth / 2;
1135 const double halfSH = apertureStartHeight / 2;
1136 const double halfEW = apertureEndWidth / 2;
1137 const double halfEH = apertureEndHeight / 2;
1138
1139 // Build the basic shape.
1140 Hexahedron hex;
1141 hex.lfb = V3D(-halfSW, -halfSH, 0);
1142 hex.lft = V3D(-halfSW, halfSH, 0);
1143 hex.lbb = V3D(-halfEW, -halfEH, length);
1144 hex.lbt = V3D(-halfEW, halfEH, length);
1145 hex.rfb = V3D(halfSW, -halfSH, 0);
1146 hex.rft = V3D(halfSW, halfSH, 0);
1147 hex.rbb = V3D(halfEW, -halfEH, length);
1148 hex.rbt = V3D(halfEW, halfEH, length);
1149
1150 // Point it along the defined axis.
1151 if (axis != DEFAULT_AXIS) {
1152 const Quat q(DEFAULT_AXIS, axis);
1153
1154 q.rotate(hex.lfb);
1155 q.rotate(hex.lft);
1156 q.rotate(hex.lbb);
1157 q.rotate(hex.lbt);
1158 q.rotate(hex.rfb);
1159 q.rotate(hex.rft);
1160 q.rotate(hex.rbb);
1161 q.rotate(hex.rbt);
1162 }
1163
1164 // Rotate the points by the rotate, rotate-all and goniometer tags
1165 if (pElem_rot) {
1166 // Apply manual rotation supplied to rotate tag
1167 std::vector<double> rotateAngles = DegreesToRadians(parsePosition(pElem_rot));
1168 const std::vector<double> rotateMatrix = generateMatrix(rotateAngles[0], rotateAngles[1], rotateAngles[2]);
1169 hex.rotatePoints(rotateMatrix);
1170 }
1171 if (m_rotateAllMatrix != Kernel::Matrix<double>(3, 3, true)) {
1172 // Apply automatic rotation due to the rotate-all tag
1174 centre.rotate(m_rotateAllMatrix);
1175 }
1176 if (m_gonioRotateMatrix != Kernel::Matrix<double>(3, 3, true)) {
1177 // Apply automatic rotation due to the goniometer that should NOT be manually defined by the user
1180 }
1181
1182 // Move it to the defined centre.
1183 hex.lfb += centre;
1184 hex.lft += centre;
1185 hex.lbb += centre;
1186 hex.lbt += centre;
1187 hex.rfb += centre;
1188 hex.rft += centre;
1189 hex.rbb += centre;
1190 hex.rbt += centre;
1191
1192 return parseHexahedronFromStruct(hex, prim, l_id);
1193}
1194
1206std::string ShapeFactory::parseTorus(Poco::XML::Element *pElem, std::map<int, std::shared_ptr<Surface>> &prim,
1207 int &l_id) {
1208 Element *pElemCentre = getShapeElement(pElem, "centre");
1209 Element *pElemAxis = getShapeElement(pElem, "axis");
1210 Element *pElemRadiusFromCentre = getShapeElement(pElem, "radius-from-centre-to-tube");
1211 Element *pElemRadiusTube = getShapeElement(pElem, "radius-tube");
1212
1213 const V3D normVec = normalize(parsePosition(pElemAxis));
1214
1215 // getDoubleAttribute can throw - put the calls above any new
1216 const double radiusCentre = getDoubleAttribute(pElemRadiusFromCentre, "val");
1217 const double radiusTube = getDoubleAttribute(pElemRadiusTube, "val");
1218
1219 // add torus
1220 auto pTorus = std::make_shared<Torus>();
1221 pTorus->setCentre(parsePosition(pElemCentre));
1222 pTorus->setNorm(normVec);
1223 pTorus->setDistanceFromCentreToTube(radiusCentre);
1224 pTorus->setTubeRadius(radiusTube);
1225 prim[l_id] = pTorus;
1226
1227 std::stringstream retAlgebraMatch;
1228 retAlgebraMatch << "(-" << l_id << ")";
1229 l_id++;
1230
1231 return retAlgebraMatch.str();
1232}
1233
1246std::string ShapeFactory::parseSliceOfCylinderRing(Poco::XML::Element *pElem,
1247 std::map<int, std::shared_ptr<Surface>> &prim, int &l_id) {
1248 Element *pElemArc = getShapeElement(pElem, "arc");
1249 Element *pElemInnerRadius = getShapeElement(pElem, "inner-radius");
1250 Element *pElemOuterRadius = getShapeElement(pElem, "outer-radius");
1251 Element *pElemDepth = getShapeElement(pElem, "depth");
1252 Element *pElem_rot = getOptionalShapeElement(pElem, "rotate");
1253
1254 const double outerRadius = getDoubleAttribute(pElemOuterRadius, "val");
1255 double innerRadius = getDoubleAttribute(pElemInnerRadius, "val");
1256 if (innerRadius <= 0.0) {
1257 innerRadius = outerRadius / 1000;
1258 g_log.warning() << "ShapeFactory::parseSliceOfCylinderRing(): inner-radius cannot be < 0.0 or = 0.0, so has been "
1259 "automatically set.";
1260 }
1261 const double middleRadius = (outerRadius + innerRadius) / 2.0;
1262 const double depth = getDoubleAttribute(pElemDepth, "val");
1263 const double arc = (M_PI / 180.0) * getDoubleAttribute(pElemArc, "val");
1264
1265 V3D normVec(0, 0, 1);
1266 V3D centrePoint(-middleRadius, 0, 0);
1267 V3D planeSlice1 = V3D(cos(arc / 2.0 + M_PI / 2.0), sin(arc / 2.0 + M_PI / 2.0), 0);
1268 V3D planeSlice2 = V3D(cos(-arc / 2.0 + M_PI / 2.0), sin(-arc / 2.0 + M_PI / 2.0), 0);
1269
1270 // Rotate the normal vector and both place slices by the rotate, rotate-all and goniometer tags
1271 // Rotate the centre by the rotate-all and goniometer tags
1272 if (pElem_rot) {
1273 // Apply manual rotation supplied to rotate tag
1274 std::vector<double> rotateAngles = DegreesToRadians(parsePosition(pElem_rot));
1275 const std::vector<double> rotateMatrix = generateMatrix(rotateAngles[0], rotateAngles[1], rotateAngles[2]);
1276 normVec.rotate(rotateMatrix);
1277 planeSlice1.rotate(rotateMatrix);
1278 planeSlice2.rotate(rotateMatrix);
1279 }
1280 if (m_rotateAllMatrix != Kernel::Matrix<double>(3, 3, true)) {
1281 // Apply automatic rotation due to the rotate-all tag
1282 normVec.rotate(m_rotateAllMatrix);
1283 planeSlice1.rotate(m_rotateAllMatrix);
1284 planeSlice2.rotate(m_rotateAllMatrix);
1285 centrePoint.rotate(m_rotateAllMatrix);
1286 }
1287 if (m_gonioRotateMatrix != Kernel::Matrix<double>(3, 3, true)) {
1288 // Apply automatic rotation due to the goniometer that should NOT be manually defined by the user
1289 normVec.rotate(m_gonioRotateMatrix);
1290 planeSlice1.rotate(m_gonioRotateMatrix);
1291 planeSlice2.rotate(m_gonioRotateMatrix);
1292 centrePoint.rotate(m_gonioRotateMatrix);
1293 }
1294
1295 // add inner infinite cylinder
1296 auto pCylinder1 = std::make_shared<Cylinder>();
1297 pCylinder1->setCentre(centrePoint);
1298 pCylinder1->setNorm(normVec);
1299 pCylinder1->setRadius(innerRadius);
1300 prim[l_id] = pCylinder1;
1301 std::stringstream retAlgebraMatch;
1302 retAlgebraMatch << "(" << l_id << " ";
1303 l_id++;
1304
1305 // add outer infinite cylinder
1306 auto pCylinder2 = std::make_shared<Cylinder>();
1307 pCylinder2->setCentre(centrePoint);
1308 pCylinder2->setNorm(normVec);
1309 pCylinder2->setRadius(outerRadius);
1310 prim[l_id] = pCylinder2;
1311 retAlgebraMatch << "-" << l_id << " ";
1312 l_id++;
1313
1314 // add top cutoff plane of infinite cylinder ring
1315 auto pPlaneTop = std::make_shared<Plane>();
1316 V3D pointInPlaneTop = centrePoint + (normVec * depth * 0.5);
1317 pPlaneTop->setPlane(pointInPlaneTop, normVec);
1318 prim[l_id] = pPlaneTop;
1319 retAlgebraMatch << "-" << l_id << " ";
1320 l_id++;
1321
1322 // add bottom cutoff plane (which is assumed to face the sample)
1323 // which at this point will result in a cylinder ring
1324 auto pPlaneBottom = std::make_shared<Plane>();
1325 V3D pointInPlaneBottom = centrePoint - (normVec * depth * 0.5);
1326 pPlaneBottom->setPlane(pointInPlaneBottom, normVec);
1327 prim[l_id] = pPlaneBottom;
1328 retAlgebraMatch << "" << l_id << " ";
1329 l_id++;
1330
1331 // the two planes that are going to cut a slice of the cylinder ring
1332
1333 auto pPlaneSlice1 = std::make_shared<Plane>();
1334 pPlaneSlice1->setPlane(centrePoint, planeSlice1);
1335 prim[l_id] = pPlaneSlice1;
1336 retAlgebraMatch << "-" << l_id << " ";
1337 l_id++;
1338
1339 auto pPlaneSlice2 = std::make_shared<Plane>();
1340 pPlaneSlice2->setPlane(centrePoint, planeSlice2);
1341 prim[l_id] = pPlaneSlice2;
1342 retAlgebraMatch << "" << l_id << ")";
1343 l_id++;
1344
1345 return retAlgebraMatch.str();
1346}
1347
1358Poco::XML::Element *ShapeFactory::getShapeElement(Poco::XML::Element *pElem, const std::string &name) {
1359 // check if this shape element contain an element with name specified by the
1360 // 2nd function argument
1361 Poco::AutoPtr<NodeList> pNL = pElem->getElementsByTagName(name);
1362 if (pNL->length() != 1) {
1363 throw std::invalid_argument("XML element: <" + pElem->tagName() +
1364 "> must contain exactly one sub-element with name: <" + name + ">.");
1365 }
1366 auto *retVal = static_cast<Element *>(pNL->item(0));
1367 return retVal;
1368}
1369
1381Poco::XML::Element *ShapeFactory::getOptionalShapeElement(Poco::XML::Element *pElem, const std::string &name) {
1382 // Allow zero or one occurances of subelements with the given name.
1383 Poco::AutoPtr<NodeList> pNL = pElem->getElementsByTagName(name);
1384 if (pNL->length() == 0)
1385 return nullptr;
1386 else if (pNL->length() > 1)
1387 throw std::invalid_argument("XML element: <" + pElem->tagName() +
1388 "> may contain at most one sub-element with name: <" + name + ">.");
1389
1390 auto *retVal = static_cast<Element *>(pNL->item(0));
1391 return retVal;
1392}
1393
1404double ShapeFactory::getDoubleAttribute(Poco::XML::Element *pElem, const std::string &name) {
1405 if (pElem->hasAttribute(name)) {
1406 return std::stod(pElem->getAttribute(name));
1407 } else {
1408 throw std::invalid_argument("XML element: <" + pElem->tagName() + "> does not have the attribute: " + name + ".");
1409 }
1410}
1411
1417V3D ShapeFactory::parsePosition(Poco::XML::Element *pElem) {
1418 V3D retVal;
1419
1420 if (pElem->hasAttribute("R") || pElem->hasAttribute("theta") || pElem->hasAttribute("phi")) {
1421 double R = 0.0, theta = 0.0, phi = 0.0;
1422
1423 if (pElem->hasAttribute("R"))
1424 R = std::stod(pElem->getAttribute("R"));
1425 if (pElem->hasAttribute("theta"))
1426 theta = std::stod(pElem->getAttribute("theta"));
1427 if (pElem->hasAttribute("phi"))
1428 phi = std::stod(pElem->getAttribute("phi"));
1429
1430 retVal.spherical(R, theta, phi);
1431 } else if (pElem->hasAttribute("r") || pElem->hasAttribute("t") || pElem->hasAttribute("p"))
1432 // This is alternative way a user may specify sphecical coordinates
1433 // which may be preferred in the long run to the more verbose of
1434 // using R, theta and phi.
1435 {
1436 double R = 0.0, theta = 0.0, phi = 0.0;
1437
1438 if (pElem->hasAttribute("r"))
1439 R = std::stod(pElem->getAttribute("r"));
1440 if (pElem->hasAttribute("t"))
1441 theta = std::stod(pElem->getAttribute("t"));
1442 if (pElem->hasAttribute("p"))
1443 phi = std::stod(pElem->getAttribute("p"));
1444
1445 retVal.spherical(R, theta, phi);
1446 } else {
1447 double x = 0.0, y = 0.0, z = 0.0;
1448
1449 if (pElem->hasAttribute("x"))
1450 x = std::stod(pElem->getAttribute("x"));
1451 if (pElem->hasAttribute("y"))
1452 y = std::stod(pElem->getAttribute("y"));
1453 if (pElem->hasAttribute("z"))
1454 z = std::stod(pElem->getAttribute("z"));
1455
1456 retVal(x, y, z);
1457 }
1458
1459 return retVal;
1460}
1461
1468std::shared_ptr<CSGObject> ShapeFactory::createSphere(const Kernel::V3D &centre, double radius) {
1469 const int surfaceID = 1;
1470 const std::map<int, std::shared_ptr<Surface>> primitives{{surfaceID, std::make_shared<Sphere>(centre, radius)}};
1471
1472 auto shape = std::make_shared<CSGObject>();
1473 shape->setObject(21, sphereAlgebra(surfaceID));
1474 shape->populate(primitives);
1475
1476 auto handler = std::make_shared<GeometryHandler>(shape);
1477 shape->setGeometryHandler(handler);
1478 detail::ShapeInfo shapeInfo;
1479 shapeInfo.setSphere(centre, radius);
1480 handler->setShapeInfo(std::move(shapeInfo));
1481
1482 shape->defineBoundingBox(radius, radius, radius, -radius, -radius, -radius);
1483 return shape;
1484}
1485
1498std::shared_ptr<CSGObject> ShapeFactory::createHexahedralShape(double xlb, double xlf, double xrf, double xrb,
1499 double ylb, double ylf, double yrf, double yrb) {
1500 Hexahedron hex;
1501 static const double ZDEPTH = 0.001;
1502 hex.lbb = V3D(xlb, ylb, 0);
1503 hex.lbt = V3D(xlb, ylb, ZDEPTH);
1504 hex.lfb = V3D(xlf, ylf, 0);
1505 hex.lft = V3D(xlf, ylf, ZDEPTH);
1506 hex.rbb = V3D(xrb, yrb, 0);
1507 hex.rbt = V3D(xrb, yrb, ZDEPTH);
1508 hex.rfb = V3D(xrf, yrf, 0);
1509 hex.rft = V3D(xrf, yrf, ZDEPTH);
1510
1511 std::map<int, std::shared_ptr<Surface>> prim;
1512 int l_id = 1;
1513 auto algebra = parseHexahedronFromStruct(hex, prim, l_id);
1514
1515 auto shape = std::make_shared<CSGObject>();
1516 shape->setObject(21, algebra);
1517 shape->populate(prim);
1518
1519 auto handler = std::make_shared<GeometryHandler>(shape);
1520 detail::ShapeInfo shapeInfo;
1521 shape->setGeometryHandler(handler);
1522
1523 shapeInfo.setHexahedron(hex.lbb, hex.lfb, hex.rfb, hex.rbb, hex.lbt, hex.lft, hex.rft, hex.rbt);
1524
1525 handler->setShapeInfo(std::move(shapeInfo));
1526
1527 shape->defineBoundingBox(std::max(xrb, xrf), yrf, ZDEPTH, std::min(xlf, xlb), ylb, 0);
1528
1529 return shape;
1530}
1531
1533void ShapeFactory::createGeometryHandler(Poco::XML::Element *pElem, const std::shared_ptr<CSGObject> &Obj) {
1534
1535 auto geomHandler = std::make_shared<GeometryHandler>(Obj);
1536 detail::ShapeInfo shapeInfo;
1537 Obj->setGeometryHandler(geomHandler);
1538
1539 if (pElem->tagName() == "cuboid") {
1540 auto corners = parseCuboid(pElem);
1541 shapeInfo.setCuboid(corners.lfb, corners.lft, corners.lbb, corners.rfb);
1542 } else if (pElem->tagName() == "hexahedron") {
1543 auto corners = parseHexahedron(pElem);
1544 shapeInfo.setHexahedron(corners.lbb, corners.lfb, corners.rfb, corners.rbb, corners.lbt, corners.lft, corners.rft,
1545 corners.rbt);
1546 } else if (pElem->tagName() == "sphere") {
1547 Element *pElemCentre = getOptionalShapeElement(pElem, "centre");
1548 Element *pElemRadius = getShapeElement(pElem, "radius");
1549 V3D centre;
1550 if (pElemCentre)
1551 centre = parsePosition(pElemCentre);
1552 shapeInfo.setSphere(centre, std::stod(pElemRadius->getAttribute("val")));
1553 } else if (pElem->tagName() == "cylinder") {
1554 Element *pElemCentre = getShapeElement(pElem, "centre-of-bottom-base");
1555 Element *pElemAxis = getShapeElement(pElem, "axis");
1556 Element *pElemRadius = getShapeElement(pElem, "radius");
1557 Element *pElemHeight = getShapeElement(pElem, "height");
1558 const V3D normVec = normalize(parsePosition(pElemAxis));
1559 shapeInfo.setCylinder(parsePosition(pElemCentre), normVec, std::stod(pElemRadius->getAttribute("val")),
1560 std::stod(pElemHeight->getAttribute("val")));
1561 } else if (pElem->tagName() == "hollow-cylinder") {
1562 Element *pElemCentre = getShapeElement(pElem, "centre-of-bottom-base");
1563 Element *pElemAxis = getShapeElement(pElem, "axis");
1564 Element *pElemInnerRadius = getShapeElement(pElem, "inner-radius");
1565 Element *pElemOuterRadius = getShapeElement(pElem, "outer-radius");
1566 Element *pElemHeight = getShapeElement(pElem, "height");
1567 V3D normVec = parsePosition(pElemAxis);
1568 normVec.normalize();
1569 shapeInfo.setHollowCylinder(parsePosition(pElemCentre), normVec, std::stod(pElemInnerRadius->getAttribute("val")),
1570 std::stod(pElemOuterRadius->getAttribute("val")),
1571 std::stod(pElemHeight->getAttribute("val")));
1572 } else if (pElem->tagName() == "cone") {
1573 Element *pElemTipPoint = getShapeElement(pElem, "tip-point");
1574 Element *pElemAxis = getShapeElement(pElem, "axis");
1575 Element *pElemAngle = getShapeElement(pElem, "angle");
1576 Element *pElemHeight = getShapeElement(pElem, "height");
1577
1578 const V3D normVec = normalize(parsePosition(pElemAxis));
1579 const double height = std::stod(pElemHeight->getAttribute("val"));
1580 const double radius = height * tan(M_PI * std::stod(pElemAngle->getAttribute("val")) / 180.0);
1581 shapeInfo.setCone(parsePosition(pElemTipPoint), normVec, radius, height);
1582 }
1583
1584 geomHandler->setShapeInfo(std::move(shapeInfo));
1585}
1586
1595Kernel::Matrix<double> ShapeFactory::generateMatrix(double xrotate, double yrotate, double zrotate) {
1596 Kernel::Matrix<double> xMatrix = generateXRotation(xrotate);
1597 Kernel::Matrix<double> yMatrix = generateYRotation(yrotate);
1598 Kernel::Matrix<double> zMatrix = generateZRotation(zrotate);
1599
1600 return zMatrix * yMatrix * xMatrix;
1601}
1602
1609 const double sinX = sin(xrotate);
1610 const double cosX = cos(xrotate);
1611 std::vector<double> matrixList = {1, 0, 0, 0, cosX, -sinX, 0, sinX, cosX};
1612 return Kernel::Matrix<double>(matrixList);
1613}
1614
1621 const double sinY = sin(yrotate);
1622 const double cosY = cos(yrotate);
1623 std::vector<double> matrixList = {cosY, 0, sinY, 0, 1, 0, -sinY, 0, cosY};
1624 return Kernel::Matrix<double>(matrixList);
1625}
1626
1633 const double sinZ = sin(zrotate);
1634 const double cosZ = cos(zrotate);
1635 std::vector<double> matrixList = {cosZ, -sinZ, 0, sinZ, cosZ, 0, 0, 0, 1};
1636 return Kernel::Matrix<double>(matrixList);
1637}
1638
1639std::string ShapeFactory::addGoniometerTag(const Kernel::Matrix<double> &rotateMatrix, std::string xml) {
1640
1641 // Delete previous goniometer from xml
1642 std::size_t foundGonioTag = xml.find("<goniometer");
1643 if (foundGonioTag != std::string::npos) {
1644 std::size_t gonioTagLength = xml.find(">", foundGonioTag + 1) - foundGonioTag;
1645 xml.erase(foundGonioTag, gonioTagLength);
1646 }
1647
1648 // Put goniometer tag in correct place in xml
1649 std::size_t gonioPlace;
1650 std::size_t foundType = xml.find("</type>");
1651 std::size_t foundSampleGeometry = xml.find("</samplegeometry");
1652
1653 if (foundType != std::string::npos) {
1654 // Add goniometer BEFORE Type end tag
1655 gonioPlace = foundType;
1656 } else if (foundSampleGeometry != std::string::npos) {
1657 // If no type tag, add goniometer BEFORE SampleGeometry end tag
1658 gonioPlace = foundSampleGeometry;
1659 } else {
1660 // If no Type or SampleGeometry tag, add goniometer to the end
1661 gonioPlace = xml.size();
1662 }
1663
1664 const std::vector<std::string> matrixElementNames = {"a11", "a12", "a13", "a21", "a22", "a23", "a31", "a32", "a33"};
1665 std::string goniometerRotate = " <goniometer ";
1666 for (size_t i = 0; i < rotateMatrix.numRows(); ++i) {
1667 for (size_t j = 0; j < rotateMatrix.numCols(); ++j) {
1668 goniometerRotate += matrixElementNames[3 * i + j] + " = '" + std::to_string(rotateMatrix[i][j]) + "' ";
1669 }
1670 }
1671 goniometerRotate += "/>";
1672 xml.insert(gonioPlace, goniometerRotate);
1673
1674 return xml;
1675}
1676} // namespace Mantid::Geometry
double height
Definition: GetAllEi.cpp:155
double radius
Definition: Rasterize.cpp:31
double innerRadius
Definition: Rasterize.cpp:39
std::string parseCone(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'cone' element.
static std::shared_ptr< CSGObject > createHexahedralShape(double xlb, double xlf, double xrf, double xrb, double ylb, double ylf, double yrf, double yrb)
Create a hexahedral shape object.
static Kernel::Matrix< double > generateYRotation(double yRotation)
Generates the y component of the rotate matrix.
std::string parseHollowCylinder(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'hollow-cylinder' element.
std::shared_ptr< CSGObject > createShape(Poco::XML::Element *pElem)
Creates a geometric object from a DOM-element-node pointing to an element whose child nodes contain t...
static std::string parseHexahedronFromStruct(const Hexahedron &hex, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
The "tapered-guide" shape is actually a special case of hexahedron; once we have the 8 points that ma...
void createGeometryHandler(Poco::XML::Element *, const std::shared_ptr< CSGObject > &)
create a special geometry handler for the known finite primitives
static Kernel::Matrix< double > generateMatrix(double xRotation, double yRotation, double zRotation)
Generates a rotate Matrix applying the x rotate then y rotate, then z rotate.
CuboidCorners parseCuboid(Poco::XML::Element *pElem)
Get the four corners of a cuboid from an XML element.
std::string addGoniometerTag(const Kernel::Matrix< double > &rotateMatrix, std::string xml)
Kernel::V3D parsePosition(Poco::XML::Element *pElem)
Get position coordinates from XML element.
std::string parseSliceOfCylinderRing(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'slice-of-cylinder-ring' element.
static std::shared_ptr< CSGObject > createSphere(const Kernel::V3D &centre, double radius)
Create a Sphere.
static Kernel::Matrix< double > generateZRotation(double zRotation)
Generates the z component of the rotate matrix.
double getDoubleAttribute(Poco::XML::Element *pElem, const std::string &name)
Return value of attribute to XML element.
Kernel::Matrix< double > m_rotateAllMatrix
Definition: ShapeFactory.h:133
Poco::XML::Element * getShapeElement(Poco::XML::Element *pElem, const std::string &name)
Return a subelement of an XML element, but also checks that there exist exactly one entry of this sub...
static std::string sphereAlgebra(const int surfaceID)
Create the algebra string for a Sphere.
std::string parseInfiniteCone(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'infinite-cone' element.
std::string parseCylinder(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'cylinder' element.
Hexahedron parseHexahedron(Poco::XML::Element *pElem)
Get all corners of a hexahedron from an XML element.
Kernel::Matrix< double > m_gonioRotateMatrix
Definition: ShapeFactory.h:132
std::string parseTorus(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'torus' element.
std::string parseSphere(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'sphere' element.
std::string parseInfinitePlane(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'infinite-plane' element.
Poco::XML::Element * getOptionalShapeElement(Poco::XML::Element *pElem, const std::string &name)
Return a subelement of an XML element.
std::string parseTaperedGuide(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'tapered-guide' element, which is a special case of the XML 'hexahedron' element.
std::string parseInfiniteCylinder(Poco::XML::Element *pElem, std::map< int, std::shared_ptr< Surface > > &prim, int &l_id)
Parse XML 'infinite-cylinder' element.
static Kernel::Matrix< double > generateXRotation(double xRotation)
Generates the x component of the rotate matrix.
void setSphere(const Kernel::V3D &center, double radius)
sets the geometry handler for a sphere
Definition: ShapeInfo.cpp:85
void setCylinder(const Kernel::V3D &centerBottomBase, const Kernel::V3D &symmetryAxis, double radius, double height)
sets the geometry handler for a cylinder
Definition: ShapeInfo.cpp:93
void setCone(const Kernel::V3D &center, const Kernel::V3D &symmetryAxis, double radius, double height)
sets the geometry handler for a cone
Definition: ShapeInfo.cpp:101
void setHollowCylinder(const Kernel::V3D &centreBottomBase, const Kernel::V3D &symmetryAxis, double innerRadius, double outerRadius, double height)
sets the geometry handler for a hollow cylinder
Definition: ShapeInfo.cpp:109
void setCuboid(const Kernel::V3D &, const Kernel::V3D &, const Kernel::V3D &, const Kernel::V3D &)
sets the geometry handler for a cuboid
Definition: ShapeInfo.cpp:68
void setHexahedron(const Kernel::V3D &, const Kernel::V3D &, const Kernel::V3D &, const Kernel::V3D &, const Kernel::V3D &, const Kernel::V3D &, const Kernel::V3D &, const Kernel::V3D &)
sets the geometry handler for a hexahedron
Definition: ShapeInfo.cpp:76
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
Numerical Matrix class.
Definition: Matrix.h:42
void identityMatrix()
Makes the matrix an idenity matrix.
Definition: Matrix.cpp:661
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
Class for quaternions.
Definition: Quat.h:39
void rotate(V3D &) const
Rotate a vector.
Definition: Quat.cpp:397
Class for 3D vectors.
Definition: V3D.h:34
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
void spherical(const double R, const double theta, const double phi) noexcept
Sets the vector position based on spherical coordinates.
Definition: V3D.cpp:57
void rotate(const Matrix< double > &) noexcept
Rotate a point by a matrix.
Definition: V3D.cpp:217
Mantid::Kernel::Logger g_log("Goniometer")
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Definition: V3D.h:341
std::string to_string(const wide_integer< Bits, Signed > &n)
void rotatePoints(const std::vector< double > &rotationMatrix)
Definition: ShapeFactory.h:44
void rotatePoints(const std::vector< double > &rotationMatrix)
Definition: ShapeFactory.h:62