Commit 228d4819 authored by Antti Soininen's avatar Antti Soininen
Browse files

Fix issues regarding zero length V3D

Most issues stemmed from trying to normalize a zero length V3D in the
unit tests. We now throw in a few places where an illegal vector
operation is tried out.

Re #24898
parent 3da2a5c6
...@@ -75,7 +75,12 @@ void DetectorSearcher::createDetectorCache() { ...@@ -75,7 +75,12 @@ void DetectorSearcher::createDetectorCache() {
auto pos = m_detInfo.position(pointNo); auto pos = m_detInfo.position(pointNo);
pos.normalize(); pos.normalize();
auto E1 = (pos - beam) * -m_crystallography_convention; auto E1 = (pos - beam) * -m_crystallography_convention;
E1.normalize(); const auto norm = E1.norm();
if (norm == 0.) {
E1 = V3D(up) * -m_crystallography_convention;
} else {
E1 /= norm;
}
Eigen::Vector3d point(E1[0], E1[1], E1[2]); Eigen::Vector3d point(E1[0], E1[1], E1[2]);
......
...@@ -441,28 +441,33 @@ void LoadIsawDetCal::doRotation(V3D rX, V3D rY, ComponentInfo &componentInfo, ...@@ -441,28 +441,33 @@ void LoadIsawDetCal::doRotation(V3D rX, V3D rY, ComponentInfo &componentInfo,
rY.normalize(); rY.normalize();
// These are the original axes // These are the original axes
const V3D oX(1., 0., 0.); constexpr V3D oX(1., 0., 0.);
const V3D oY(0., 1., 0.); constexpr V3D oY(0., 1., 0.);
// Axis that rotates X // Axis that rotates X
V3D ax1 = oX.cross_prod(rX); const V3D ax1 = oX.cross_prod(rX);
// Rotation angle from oX to rX Quat Q1;
double angle1 = oX.angle(rX) * DegreesPerRadian; if (!ax1.nullVector(1e-12)) {
if (doWishCorrection) // Rotation angle from oX to rX
angle1 += 180.0; double angle1 = oX.angle(rX) * DegreesPerRadian;
// Create the first quaternion if (doWishCorrection)
Quat Q1(angle1, ax1); angle1 += 180.0;
// Create the first quaternion
Q1.setAngleAxis(angle1, ax1);
}
// Now we rotate the original Y using Q1 // Now we rotate the original Y using Q1
V3D roY = oY; V3D roY = oY;
Q1.rotate(roY); Q1.rotate(roY);
// Find the axis that rotates oYr onto rY // Find the axis that rotates oYr onto rY
V3D ax2 = roY.cross_prod(rY); const V3D ax2 = roY.cross_prod(rY);
const double angle2 = roY.angle(rY) * DegreesPerRadian; Quat Q2;
Quat Q2(angle2, ax2); if (!ax2.nullVector(1e-12)) {
const double angle2 = roY.angle(rY) * DegreesPerRadian;
Q2.setAngleAxis(angle2, ax2);
}
// Final = those two rotations in succession; Q1 is done first. // Final = those two rotations in succession; Q1 is done first.
Quat Rot = Q2 * Q1; const Quat Rot = Q2 * Q1;
// Then find the corresponding relative position // Then find the corresponding relative position
const auto componentIndex = componentInfo.indexOf(comp->getComponentID()); const auto componentIndex = componentInfo.indexOf(comp->getComponentID());
......
...@@ -55,7 +55,6 @@ public: ...@@ -55,7 +55,6 @@ public:
if (det != nullptr) { if (det != nullptr) {
const auto detPos = det->getPos(); const auto detPos = det->getPos();
const V3D testPos(x, y, z); const V3D testPos(x, y, z);
std::cout << detPos << " " << testPos << std::endl;
TS_ASSERT_EQUALS(detPos, testPos); TS_ASSERT_EQUALS(detPos, testPos);
} else { } else {
throw std::runtime_error("In checkPosition IComponent is nullptr"); throw std::runtime_error("In checkPosition IComponent is nullptr");
...@@ -67,7 +66,6 @@ public: ...@@ -67,7 +66,6 @@ public:
if (det != nullptr) { if (det != nullptr) {
const auto detRot = det->getRotation(); const auto detRot = det->getRotation();
const Quat testRot(w, a, b, c); const Quat testRot(w, a, b, c);
std::cout << detRot << " " << testRot << std::endl;
TS_ASSERT_EQUALS(detRot, testRot); TS_ASSERT_EQUALS(detRot, testRot);
} else { } else {
throw std::runtime_error("In checkRotation IComponent is nullptr"); throw std::runtime_error("In checkRotation IComponent is nullptr");
......
...@@ -191,11 +191,14 @@ const DblMatrix &OrientedLattice::setUFromVectors(const V3D &u, const V3D &v) { ...@@ -191,11 +191,14 @@ const DblMatrix &OrientedLattice::setUFromVectors(const V3D &u, const V3D &v) {
throw std::invalid_argument("|B.v|~0"); throw std::invalid_argument("|B.v|~0");
buVec.normalize(); // 1st unit vector, along Bu buVec.normalize(); // 1st unit vector, along Bu
V3D bwVec = buVec.cross_prod(bvVec); V3D bwVec = buVec.cross_prod(bvVec);
if (bwVec.normalize() < 1e-5) const auto norm = bwVec.norm();
throw std::invalid_argument( if (norm < 1e-5) {
"u and v are parallel"); // 3rd unit vector, perpendicular to Bu,Bv // 3rd unit vector, perpendicular to Bu,Bv
bvVec = bwVec.cross_prod( throw std::invalid_argument("u and v are parallel");
buVec); // 2nd unit vector, perpendicular to Bu, in the Bu,Bv plane }
bwVec /= norm;
// 2nd unit vector, perpendicular to Bu, in the Bu,Bv plane
bvVec = bwVec.cross_prod(buVec);
DblMatrix tau(3, 3), lab(3, 3), U(3, 3); DblMatrix tau(3, 3), lab(3, 3), U(3, 3);
/*lab = U tau /*lab = U tau
/ 0 1 0 \ /bu[0] bv[0] bw[0]\ / 0 1 0 \ /bu[0] bv[0] bw[0]\
......
...@@ -2021,11 +2021,9 @@ void InstrumentDefinitionParser::makeXYplaneFaceComponent( ...@@ -2021,11 +2021,9 @@ void InstrumentDefinitionParser::makeXYplaneFaceComponent(
//----------------------------------------------------------------------------------------------------------------------- //-----------------------------------------------------------------------------------------------------------------------
/** Make the shape defined in 1st argument face the position in the second /** Make the shape defined in 1st argument face the position in the second
*argument, * argument, by rotating the z-axis of the component passed in 1st argument so
* by rotating the z-axis of the component passed in 1st argument so that it * that it points in the direction: from the position (as specified 2nd
*points in the * argument) to the component (1st argument).
* direction: from the position (as specified 2nd argument) to the component
*(1st argument).
* *
* @param in :: Component to be rotated * @param in :: Component to be rotated
* @param facingPoint :: position to face * @param facingPoint :: position to face
...@@ -2036,22 +2034,22 @@ void InstrumentDefinitionParser::makeXYplaneFaceComponent( ...@@ -2036,22 +2034,22 @@ void InstrumentDefinitionParser::makeXYplaneFaceComponent(
// vector from facing object to component we want to rotate // vector from facing object to component we want to rotate
Kernel::V3D facingDirection = pos - facingPoint; Kernel::V3D facingDirection = pos - facingPoint;
facingDirection.normalize();
if (facingDirection.norm() == 0.0) if (facingDirection.norm() == 0.0)
return; return;
facingDirection.normalize();
// now aim to rotate shape such that the z-axis of of the object we want to // now aim to rotate shape such that the z-axis of of the object we want to
// rotate // rotate points in the direction of facingDirection. That way the XY plane
// points in the direction of facingDirection. That way the XY plane faces // faces the 'facing object'.
// the constexpr Kernel::V3D z(0, 0, 1);
// 'facing object'.
Kernel::V3D z = Kernel::V3D(0, 0, 1);
Kernel::Quat R = in->getRotation(); Kernel::Quat R = in->getRotation();
R.inverse(); R.inverse();
R.rotate(facingDirection); R.rotate(facingDirection);
Kernel::V3D normal = facingDirection.cross_prod(z); Kernel::V3D normal = facingDirection.cross_prod(z);
if (normal.norm() == 0.) {
normal = -facingDirection;
}
normal.normalize(); normal.normalize();
double theta = (180.0 / M_PI) * facingDirection.angle(z); double theta = (180.0 / M_PI) * facingDirection.angle(z);
...@@ -2059,8 +2057,7 @@ void InstrumentDefinitionParser::makeXYplaneFaceComponent( ...@@ -2059,8 +2057,7 @@ void InstrumentDefinitionParser::makeXYplaneFaceComponent(
in->rotate(Kernel::Quat(-theta, normal)); in->rotate(Kernel::Quat(-theta, normal));
else { else {
// To take into account the case where the facing direction is in the // To take into account the case where the facing direction is in the
// (0,0,1) // (0,0,1) or (0,0,-1) direction.
// or (0,0,-1) direction.
in->rotate(Kernel::Quat(-theta, Kernel::V3D(0, 1, 0))); in->rotate(Kernel::Quat(-theta, Kernel::V3D(0, 1, 0)));
} }
} }
......
...@@ -115,11 +115,11 @@ int ObjComponent::interceptSurface(Track &track) const { ...@@ -115,11 +115,11 @@ int ObjComponent::interceptSurface(Track &track) const {
// TODO: If scaling parameters are ever enabled, would they need need to be // TODO: If scaling parameters are ever enabled, would they need need to be
// used here? // used here?
V3D trkStart = factorOutComponentPosition(track.startPoint()); const V3D trkStart = factorOutComponentPosition(track.startPoint());
V3D trkDirection = takeOutRotation(track.direction()); const V3D trkDirection = takeOutRotation(track.direction());
Track probeTrack(trkStart, trkDirection); Track probeTrack(trkStart, trkDirection);
int intercepts = shape()->interceptSurface(probeTrack); const int intercepts = shape()->interceptSurface(probeTrack);
Track::LType::const_iterator it; Track::LType::const_iterator it;
for (it = probeTrack.cbegin(); it != probeTrack.cend(); ++it) { for (it = probeTrack.cbegin(); it != probeTrack.cend(); ++it) {
......
...@@ -22,7 +22,7 @@ using Kernel::V3D; ...@@ -22,7 +22,7 @@ using Kernel::V3D;
/** /**
* Default constructor * Default constructor
*/ */
Track::Track() : m_startPoint(), m_unitVector() {} Track::Track() : m_startPoint(), m_unitVector(0., 0., 1.) {}
/** /**
* Constructor * Constructor
...@@ -30,7 +30,11 @@ Track::Track() : m_startPoint(), m_unitVector() {} ...@@ -30,7 +30,11 @@ Track::Track() : m_startPoint(), m_unitVector() {}
* @param unitVector :: Directional vector. It must be unit vector. * @param unitVector :: Directional vector. It must be unit vector.
*/ */
Track::Track(const V3D &startPt, const V3D &unitVector) Track::Track(const V3D &startPt, const V3D &unitVector)
: m_startPoint(startPt), m_unitVector(unitVector) {} : m_startPoint(startPt), m_unitVector(unitVector) {
if (!unitVector.unitVector()) {
throw std::invalid_argument("Failed to construct track: direction is not a unit vector.");
}
}
/** /**
* Resets the track starting point and direction. * Resets the track starting point and direction.
...@@ -38,6 +42,9 @@ Track::Track(const V3D &startPt, const V3D &unitVector) ...@@ -38,6 +42,9 @@ Track::Track(const V3D &startPt, const V3D &unitVector)
* @param direction :: The new direction. Must be a unit vector! * @param direction :: The new direction. Must be a unit vector!
*/ */
void Track::reset(const V3D &startPoint, const V3D &direction) { void Track::reset(const V3D &startPoint, const V3D &direction) {
if (!direction.unitVector()) {
throw std::invalid_argument("Failed to reset track: direction is not a unit vector.");
}
m_startPoint = startPoint; m_startPoint = startPoint;
m_unitVector = direction; m_unitVector = direction;
} }
......
...@@ -154,17 +154,25 @@ public: ...@@ -154,17 +154,25 @@ public:
true); true);
// Mix // Mix
V3D dir(1., 1., 0);
dir.normalize();
TS_ASSERT_EQUALS( TS_ASSERT_EQUALS(
bbox.doesLineIntersect(Track(V3D(-5.0, -1.0, 0.0), V3D(1.0, 1.0, 0.0))), bbox.doesLineIntersect(Track(V3D(-5.0, -1.0, 0.0), dir)),
true); true);
dir = {1., 1., 1.,};
dir.normalize();
TS_ASSERT_EQUALS(bbox.doesLineIntersect( TS_ASSERT_EQUALS(bbox.doesLineIntersect(
Track(V3D(-5.0, -1.0, -0.5), V3D(1.0, 1.0, 1.0))), Track(V3D(-5.0, -1.0, -0.5), dir)),
true); true);
dir = {-1., -0.4, 0.};
dir.normalize();
TS_ASSERT_EQUALS(bbox.doesLineIntersect( TS_ASSERT_EQUALS(bbox.doesLineIntersect(
Track(V3D(10.0, 10.0, 0.0), V3D(-1.0, -0.4, 0.0))), Track(V3D(10.0, 10.0, 0.0), dir)),
false); false);
dir = {1., 1., 0.};
dir.normalize();
TS_ASSERT_EQUALS(bbox.doesLineIntersect( TS_ASSERT_EQUALS(bbox.doesLineIntersect(
Track(V3D(-10.0, -10.0, 0.0), V3D(1.0, 1.0, 0.0))), Track(V3D(-10.0, -10.0, 0.0), dir)),
true); // Hits box at edge true); // Hits box at edge
} }
......
...@@ -390,7 +390,9 @@ public: ...@@ -390,7 +390,9 @@ public:
std::vector<Link> std::vector<Link>
expectedResults; // left empty as there are no expected results expectedResults; // left empty as there are no expected results
auto geom_obj = createCappedCylinder(); auto geom_obj = createCappedCylinder();
Track track(V3D(-10, 0, 0), V3D(1, 1, 0)); V3D dir(1., 1., 0.);
dir.normalize();
Track track(V3D(-10, 0, 0), dir);
checkTrackIntercept(geom_obj, track, expectedResults); checkTrackIntercept(geom_obj, track, expectedResults);
} }
......
...@@ -189,10 +189,10 @@ public: ...@@ -189,10 +189,10 @@ public:
Cylinder A; Cylinder A;
A.setSurface("cx 5"); A.setSurface("cx 5");
TS_ASSERT_EQUALS(A.surfaceNormal(V3D(10, 0, 0)), V3D(0, 0, 0)); TS_ASSERT_THROWS_ANYTHING(A.surfaceNormal(V3D(10, 0, 0)));
TS_ASSERT_EQUALS(A.surfaceNormal(V3D(0, 10, 0)), V3D(0, 1, 0)); TS_ASSERT_EQUALS(A.surfaceNormal(V3D(0, 10, 0)), V3D(0, 1, 0));
TS_ASSERT_EQUALS(A.surfaceNormal(V3D(0, 0, 10)), V3D(0, 0, 1)); TS_ASSERT_EQUALS(A.surfaceNormal(V3D(0, 0, 10)), V3D(0, 0, 1));
TS_ASSERT_EQUALS(A.surfaceNormal(V3D(-10, 0, 0)), V3D(0, 0, 0)); TS_ASSERT_THROWS_ANYTHING(A.surfaceNormal(V3D(-10, 0, 0)));
TS_ASSERT_EQUALS(A.surfaceNormal(V3D(0, -10, 0)), V3D(0, -1, 0)); TS_ASSERT_EQUALS(A.surfaceNormal(V3D(0, -10, 0)), V3D(0, -1, 0));
TS_ASSERT_EQUALS(A.surfaceNormal(V3D(0, 0, -10)), V3D(0, 0, -1)); TS_ASSERT_EQUALS(A.surfaceNormal(V3D(0, 0, -10)), V3D(0, 0, -1));
......
...@@ -258,8 +258,8 @@ public: ...@@ -258,8 +258,8 @@ public:
TS_ASSERT(!ptrDet4->isValid(V3D(-0.02, 0.0, 0.0) + ptrDet4->getPos())); TS_ASSERT(!ptrDet4->isValid(V3D(-0.02, 0.0, 0.0) + ptrDet4->getPos()));
TS_ASSERT(!ptrDet4->isValid(V3D(0.0, 0.02, 0.0) + ptrDet4->getPos())); TS_ASSERT(!ptrDet4->isValid(V3D(0.0, 0.02, 0.0) + ptrDet4->getPos()));
TS_ASSERT(!ptrDet4->isValid(V3D(0.0, -0.02, 0.0) + ptrDet4->getPos())); TS_ASSERT(!ptrDet4->isValid(V3D(0.0, -0.02, 0.0) + ptrDet4->getPos()));
TS_ASSERT(!ptrDet4->isValid(V3D(0.0, 0.0, 0.02) + ptrDet4->getPos())); TS_ASSERT(ptrDet4->isValid(V3D(0.0, 0.0, 0.02) + ptrDet4->getPos()));
TS_ASSERT(ptrDet4->isValid(V3D(0.0, 0.0, -0.02) + ptrDet4->getPos())); TS_ASSERT(!ptrDet4->isValid(V3D(0.0, 0.0, -0.02) + ptrDet4->getPos()));
// test of facing as a sub-element of location // test of facing as a sub-element of location
boost::shared_ptr<const IDetector> ptrDet5 = i->getDetector(5); boost::shared_ptr<const IDetector> ptrDet5 = i->getDetector(5);
......
...@@ -120,6 +120,7 @@ public: ...@@ -120,6 +120,7 @@ public:
Instrument_sptr testInst = setupInstrument(); Instrument_sptr testInst = setupInstrument();
InstrumentRayTracer tracker(testInst); InstrumentRayTracer tracker(testInst);
V3D testDir(0.010, 0.0, 15.004); V3D testDir(0.010, 0.0, 15.004);
testDir.normalize();
tracker.trace(testDir); tracker.trace(testDir);
Links results = tracker.getResults(); Links results = tracker.getResults();
TS_ASSERT_EQUALS(results.size(), 1); TS_ASSERT_EQUALS(results.size(), 1);
...@@ -156,9 +157,13 @@ public: ...@@ -156,9 +157,13 @@ public:
*/ */
void doTestRectangularDetector(std::string message, Instrument_sptr inst, void doTestRectangularDetector(std::string message, Instrument_sptr inst,
V3D testDir, int expectX, int expectY) { V3D testDir, int expectX, int expectY) {
// std::cout << message << '\n';
InstrumentRayTracer tracker(inst); InstrumentRayTracer tracker(inst);
testDir.normalize(); // Force to be unit vector const auto norm = testDir.norm();
if (norm == 0.) {
TS_ASSERT_THROWS_ANYTHING(tracker.traceFromSample(testDir));
return;
}
testDir /= norm;
tracker.traceFromSample(testDir); tracker.traceFromSample(testDir);
Links results = tracker.getResults(); Links results = tracker.getResults();
......
...@@ -320,7 +320,9 @@ public: ...@@ -320,7 +320,9 @@ public:
std::vector<Link> std::vector<Link>
expectedResults; // left empty as there are no expected results expectedResults; // left empty as there are no expected results
auto geom_obj = createCube(4.0); auto geom_obj = createCube(4.0);
Track track(V3D(-10, 0, 0), V3D(1, 1, 0)); V3D dir(1., 1., 0.);
dir.normalize();
Track track(V3D(-10, 0, 0), dir);
checkTrackIntercept(std::move(geom_obj), track, expectedResults); checkTrackIntercept(std::move(geom_obj), track, expectedResults);
} }
...@@ -361,7 +363,9 @@ public: ...@@ -361,7 +363,9 @@ public:
void testInterceptLShapeTwoPass() { void testInterceptLShapeTwoPass() {
std::vector<Link> expectedResults; std::vector<Link> expectedResults;
auto geom_obj = createLShape(); auto geom_obj = createLShape();
Track track(V3D(0, 2.5, 0.5), V3D(0.707, -0.707, 0)); V3D dir(1., -1., 0.);
dir.normalize();
Track track(V3D(0, 2.5, 0.5), dir);
// format = startPoint, endPoint, total distance so far // format = startPoint, endPoint, total distance so far
expectedResults.emplace_back( expectedResults.emplace_back(
......
...@@ -164,8 +164,10 @@ public: ...@@ -164,8 +164,10 @@ public:
ocyl.setParent(&parent); ocyl.setParent(&parent);
// Check original track misses // Check original track misses
TS_ASSERT_EQUALS(ocyl.interceptSurface(track), 0); TS_ASSERT_EQUALS(ocyl.interceptSurface(track), 0);
// Create a new test track going from the origin down the line y = -x // Create a new test track going from the origin down the line z = -y
Track track2(V3D(0, 0, 0), V3D(0, 1, -1)); V3D dir(0., 1., -1.);
dir.normalize();
Track track2(V3D(0, 0, 0), dir);
TS_ASSERT_EQUALS(ocyl.interceptSurface(track2), 1); TS_ASSERT_EQUALS(ocyl.interceptSurface(track2), 1);
Track::LType::const_iterator it2 = track2.cbegin(); Track::LType::const_iterator it2 = track2.cbegin();
if (it2 == track2.cend()) if (it2 == track2.cend())
...@@ -347,14 +349,14 @@ public: ...@@ -347,14 +349,14 @@ public:
TS_ASSERT_EQUALS(itscale->entryPoint, V3D(-6.4, 0.0, 0.0)); TS_ASSERT_EQUALS(itscale->entryPoint, V3D(-6.4, 0.0, 0.0));
TS_ASSERT_EQUALS(itscale->exitPoint, V3D(2.4, 0.0, 0.0)); TS_ASSERT_EQUALS(itscale->exitPoint, V3D(2.4, 0.0, 0.0));
Track trackScaleY(V3D(0.0, -2, 0), V3D(0, 2.0, 0)); Track trackScaleY(V3D(0.0, -2, 0), V3D(0, 1.0, 0));
TS_ASSERT_EQUALS(ocyl->interceptSurface(trackScaleY), 1); TS_ASSERT_EQUALS(ocyl->interceptSurface(trackScaleY), 1);
Track::LType::const_iterator itscaleY = trackScaleY.cbegin(); Track::LType::const_iterator itscaleY = trackScaleY.cbegin();
TS_ASSERT_EQUALS(itscaleY->distFromStart, 2.5); TS_ASSERT_EQUALS(itscaleY->distFromStart, 2.5);
TS_ASSERT_EQUALS(itscaleY->entryPoint, V3D(0.0, -0.5, 0.0)); TS_ASSERT_EQUALS(itscaleY->entryPoint, V3D(0.0, -0.5, 0.0));
TS_ASSERT_EQUALS(itscaleY->exitPoint, V3D(0.0, 0.5, 0.0)); TS_ASSERT_EQUALS(itscaleY->exitPoint, V3D(0.0, 0.5, 0.0));
Track trackScaleW(V3D(0, 0, -5), V3D(0, 0, 5)); Track trackScaleW(V3D(0, 0, -5), V3D(0, 0, 1));
TS_ASSERT_EQUALS(ocyl->interceptSurface(trackScaleW), 1); TS_ASSERT_EQUALS(ocyl->interceptSurface(trackScaleW), 1);
Track::LType::const_iterator itscaleW = trackScaleW.cbegin(); Track::LType::const_iterator itscaleW = trackScaleW.cbegin();
TS_ASSERT_DELTA(itscaleW->distFromStart, 6.5, 1e-6); TS_ASSERT_DELTA(itscaleW->distFromStart, 6.5, 1e-6);
......
...@@ -176,8 +176,10 @@ public: ...@@ -176,8 +176,10 @@ public:
ocyl.setParent(&parent); ocyl.setParent(&parent);
// Check original track misses // Check original track misses
TS_ASSERT_EQUALS(pocyl.interceptSurface(track), 0); TS_ASSERT_EQUALS(pocyl.interceptSurface(track), 0);
// Create a new test track going from the origin down the line y = -x // Create a new test track going from the origin down the line z = -y
Track track2(V3D(0, 0, 0), V3D(0, 1, -1)); V3D dir(0, 1., -1.);
dir.normalize();
Track track2(V3D(0, 0, 0), dir);
TS_ASSERT_EQUALS(pocyl.interceptSurface(track2), 1); TS_ASSERT_EQUALS(pocyl.interceptSurface(track2), 1);
Track::LType::const_iterator it2 = track2.cbegin(); Track::LType::const_iterator it2 = track2.cbegin();
if (it2 == track2.cend()) if (it2 == track2.cend())
......
...@@ -113,7 +113,7 @@ public: ...@@ -113,7 +113,7 @@ public:
Quat &operator*=(const Quat &); Quat &operator*=(const Quat &);
bool operator==(const Quat &) const; bool operator==(const Quat &) const;
bool operator!=(const Quat &) const; bool operator!=(const Quat &) const;
const double &operator[](int) const; double operator[](int) const;
double &operator[](int); double &operator[](int);
/** @name Element access. */ /** @name Element access. */
......
...@@ -351,6 +351,7 @@ public: ...@@ -351,6 +351,7 @@ public:
1e-3) const; ///< Determine if there is a master direction 1e-3) const; ///< Determine if there is a master direction
bool bool
nullVector(const double Tol = 1e-3) const; ///< Determine if the point is null nullVector(const double Tol = 1e-3) const; ///< Determine if the point is null
bool unitVector(const double tolerance = Kernel::Tolerance) const noexcept;
bool coLinear(const V3D &, const V3D &) const; bool coLinear(const V3D &, const V3D &) const;
void saveNexus(::NeXus::File *file, const std::string &name) const; void saveNexus(::NeXus::File *file, const std::string &name) const;
......
...@@ -213,27 +213,31 @@ void Quat::operator()(const V3D &rX, const V3D &rY, const V3D &rZ) { ...@@ -213,27 +213,31 @@ void Quat::operator()(const V3D &rX, const V3D &rY, const V3D &rZ) {
(void)rZ; // Avoid compiler warning (void)rZ; // Avoid compiler warning
// These are the original axes // These are the original axes
V3D oX = V3D(1., 0., 0.); constexpr V3D oX = V3D(1., 0., 0.);
V3D oY = V3D(0., 1., 0.); constexpr V3D oY = V3D(0., 1., 0.);
// Axis that rotates X // Axis that rotates X
V3D ax1 = oX.cross_prod(rX); const V3D ax1 = oX.cross_prod(rX);
// Rotation angle from oX to rX
double angle1 = oX.angle(rX);
// Create the first quaternion // Create the first quaternion
Quat Q1(angle1 * 180.0 / M_PI, ax1); Quat Q1;
if (!ax1.nullVector()) {
// Rotation angle from oX to rX
const double angle1 = oX.angle(rX);
Q1.setAngleAxis(angle1 * 180.0 / M_PI, ax1);
}
// Now we rotate the original Y using Q1 // Now we rotate the original Y using Q1
V3D roY = oY; V3D roY = oY;
Q1.rotate(roY); Q1.rotate(roY);
// Find the axis that rotates oYr onto rY // Find the axis that rotates oYr onto rY
V3D ax2 = roY.cross_prod(rY); const V3D ax2 = roY.cross_prod(rY);
double angle2 = roY.angle(rY); Quat Q2;
Quat Q2(angle2 * 180.0 / M_PI, ax2); if (!ax2.nullVector()) {
const double angle2 = roY.angle(rY);
Q2.setAngleAxis(angle2 * 180.0 / M_PI, ax2);
}
// Final = those two rotations in succession; Q1 is done first. // Final = those two rotations in succession; Q1 is done first.
Quat final = Q2 * Q1; const Quat final = Q2 * Q1;