Mantid
Loading...
Searching...
No Matches
V3D.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#include "MantidKernel/V3D.h"
9#include "MantidKernel/Quat.h"
10#include <algorithm>
11#include <boost/version.hpp>
12#include <sstream>
13
14#if BOOST_VERSION < 106700
15#include <boost/math/common_factor.hpp>
16using boost::math::gcd;
17#else
18#include <boost/integer/common_factor.hpp>
19using boost::integer::gcd;
20#endif
21
22#include <nexus/NeXusFile.hpp>
23
24namespace {
32double nearInt(double val, double eps, double mult) noexcept {
33 if (val > 0.0) {
34 if (val < 1.0) {
35 mult /= val;
36 } else {
37 if (std::abs(val - std::round(val)) > eps) {
38 mult *= std::ceil(val / eps) * eps / val;
39 }
40 }
41 }
42 return mult;
43}
44} // namespace
45
46namespace Mantid::Kernel {
47
57void V3D::spherical(const double R, const double theta, const double phi) noexcept {
58 constexpr double deg2rad = M_PI / 180.0;
59 spherical_rad(R, theta * deg2rad, phi * deg2rad);
60}
61
70void V3D::spherical_rad(const double R, const double polar, const double azimuth) noexcept {
71 m_pt[2] = R * cos(polar);
72 const double ct = R * sin(polar);
73 m_pt[0] = ct * cos(azimuth);
74 m_pt[1] = ct * sin(azimuth);
75
76 // Setting this way can lead to very small values of x & y that should really
77 // be zero.
78 // This can cause confusion for the atan2 function used in getSpherical.
79 if (std::abs(m_pt[0]) < Tolerance)
80 m_pt[0] = 0.0;
81 if (std::abs(m_pt[1]) < Tolerance)
82 m_pt[1] = 0.0;
83}
84
95void V3D::azimuth_polar_SNS(const double R, const double azimuth, const double polar) noexcept {
96 m_pt[1] = R * cos(polar);
97 const double ct = R * sin(polar);
98 m_pt[0] = ct * cos(azimuth);
99 m_pt[2] = ct * sin(azimuth);
100
101 // Setting this way can lead to very small values of x & y that should really
102 // be zero.
103 // This can cause confusion for the atan2 function used in getSpherical.
104 if (std::abs(m_pt[0]) < Tolerance)
105 m_pt[0] = 0.0;
106 if (std::abs(m_pt[1]) < Tolerance)
107 m_pt[1] = 0.0;
108 if (std::abs(m_pt[2]) < Tolerance)
109 m_pt[2] = 0.0;
110}
111
117void V3D::getSpherical(double &R, double &theta, double &phi) const noexcept {
118 constexpr double rad2deg = 180.0 / M_PI;
119 R = norm();
120 theta = 0.0;
121 if (R != 0.0)
122 theta = acos(m_pt[2] / R) * rad2deg;
123 phi = atan2(m_pt[1], m_pt[0]) * rad2deg;
124}
125
131 const double ND(norm());
132 if (ND == 0) {
133 throw std::runtime_error("Unable to normalize a zero length vector.");
134 }
135 this->operator/=(ND);
136 return ND;
137}
138
140void V3D::round() noexcept {
141 for (auto &p : m_pt) {
142 p = std::round(p);
143 }
144}
145
150double V3D::zenith(const V3D &v) const noexcept {
151 const double R = distance(v);
152 if (R != 0.0) {
153 const double zOffset = m_pt[2] - v.m_pt[2];
154 return acos(zOffset / R);
155 } else {
156 return 0.0;
157 }
158}
159
165double V3D::angle(const V3D &v) const {
166 const double ratio = this->cosAngle(v);
167 if (ratio >= 1.0) // NOTE: Due to rounding errors, if v is
168 return 0.0; // is nearly the same as "this" or
169 else if (ratio <= -1.0) // as "-this", ratio can be slightly
170 return M_PI; // more than 1 in absolute value.
171 // That causes acos() to return NaN.
172 return acos(ratio);
173}
174
180double V3D::cosAngle(const V3D &v) const {
181 const double n1 = norm();
182 const double n2 = v.norm();
183 if (n1 == 0. || n2 == 0.) {
184 throw std::runtime_error("Cannot calculate an angle when one of the vectors has zero length.");
185 }
186 return this->scalar_prod(v) / (n1 * n2);
187}
188
199int V3D::reBase(const V3D &A, const V3D &B, const V3D &C) noexcept {
200 Matrix<double> T(3, 3);
201 for (int i = 0; i < 3; i++) {
202 T[i][0] = A[i];
203 T[i][1] = B[i];
204 T[i][2] = C[i];
205 }
206 const double det = T.Invert();
207 if (std::abs(det) < 1e-13) // failed
208 return -1;
209 rotate(T);
210 return 0;
211}
212
217void V3D::rotate(const Kernel::Matrix<double> &A) noexcept {
218 const double xold(m_pt[0]), yold(m_pt[1]), zold(m_pt[2]);
219 for (size_t i = 0; i < m_pt.size(); ++i) {
220 m_pt[i] = A[i][0] * xold + A[i][1] * yold + A[i][2] * zold;
221 }
222}
223
230bool V3D::coLinear(const V3D &Bv, const V3D &Cv) const noexcept {
231 const V3D &Av = *this;
232 const V3D Tmp((Bv - Av).cross_prod(Cv - Av));
233 return Tmp.norm() <= Tolerance;
234}
235
241bool V3D::nullVector(const double tolerance) const noexcept {
242 for (const double p : m_pt) {
243 if (std::abs(p) > tolerance) {
244 return false;
245 }
246 }
247 // Getting to this point means a null vector
248 return true;
249}
250
251bool V3D::unitVector(const double tolerance) const noexcept {
252 const auto l = norm();
253 return std::abs(l - 1.) < tolerance;
254}
255
264int V3D::masterDir(const double tolerance) const noexcept {
265 // Calc max dist
266 double max = m_pt[0] * m_pt[0];
267 double other = max;
268 double u2 = m_pt[1] * m_pt[1];
269 int idx = (m_pt[0] > 0) ? 1 : -1;
270 if (u2 > max) {
271 max = u2;
272 idx = (m_pt[1] > 0) ? 2 : -2;
273 }
274 other += u2;
275 u2 = m_pt[2] * m_pt[2];
276 if (u2 > max) {
277 max = u2;
278 idx = (m_pt[2] > 0) ? 3 : -3;
279 }
280 other += u2;
281 other -= max;
282 if ((other / max) > tolerance) // doesn't have master direction
283 {
284 return 0;
285 }
286 return idx;
287}
288
298std::vector<V3D> V3D::makeVectorsOrthogonal(const std::vector<V3D> &vectors) {
299 if (vectors.size() != 2)
300 throw std::invalid_argument("makeVectorsOrthogonal() only works with 2 vectors");
301
302 const V3D v0 = Kernel::normalize(vectors[0]);
303 V3D v1 = Kernel::normalize(vectors[1]);
304
305 std::vector<V3D> out;
306 out.reserve(3);
307 out.emplace_back(v0);
308
309 // Make a rotation 90 degrees from 0 to 1
310 Quat q(v0, v1);
311 q.setRotation(90);
312 // Rotate v1 so it is 90 deg
313 v1 = v0;
314 q.rotate(v1);
315 out.emplace_back(v1);
316
317 // Finally, the 3rd vector = cross product of 0 and 1
318 V3D v2 = v0.cross_prod(v1);
319 out.emplace_back(v2);
320 return out;
321}
322
328void V3D::read(std::istream &IX) { IX >> m_pt[0] >> m_pt[1] >> m_pt[2]; }
329
330void V3D::write(std::ostream &OX) const
335{
336 OX << m_pt[0] << " " << m_pt[1] << " " << m_pt[2];
337}
338
340std::string V3D::toString() const {
341 std::ostringstream mess;
342 this->write(mess);
343 return mess.str();
344}
345
348void V3D::fromString(const std::string &str) {
349 std::istringstream mess(str);
350 this->read(mess);
351}
352
357void V3D::printSelf(std::ostream &os) const { os << "[" << m_pt[0] << "," << m_pt[1] << "," << m_pt[2] << "]"; }
358
364void V3D::readPrinted(std::istream &IX) {
365 std::string in;
366 std::getline(IX, in);
367 size_t i = in.find_first_of('[');
368 if (i == std::string::npos)
369 throw std::runtime_error("Wrong format for V3D input: " + in);
370 size_t j = in.find_last_of(']');
371 if (j == std::string::npos || j < i + 6)
372 throw std::runtime_error("Wrong format for V3D input: " + in);
373
374 size_t c1 = in.find_first_of(',');
375 size_t c2 = in.find_first_of(',', c1 + 1);
376 if (c1 == std::string::npos || c2 == std::string::npos)
377 throw std::runtime_error("Wrong format for V3D input: [" + in + "]");
378
379 m_pt[0] = std::stod(in.substr(i + 1, c1 - i - 1));
380 m_pt[1] = std::stod(in.substr(c1 + 1, c2 - c1 - 1));
381 m_pt[2] = std::stod(in.substr(c2 + 1, j - c2 - 1));
382}
383
390std::ostream &operator<<(std::ostream &os, const V3D &v) {
391 v.printSelf(os);
392 return os;
393}
394
395std::istream &operator>>(std::istream &IX, V3D &A)
402{
403 A.readPrinted(IX);
404 return IX;
405}
406
407//--------------------------------------------------------------------------------------------
412void V3D::saveNexus(::NeXus::File *file, const std::string &name) const {
413 file->makeData(name, ::NeXus::FLOAT64, 3, true);
414 file->putData(m_pt.data());
415 file->closeData();
416}
417
418//--------------------------------------------------------------------------------------------
423void V3D::loadNexus(::NeXus::File *file, const std::string &name) {
424 std::vector<double> data;
425 file->readData(name, data);
426 if (data.size() != 3)
427 throw std::runtime_error("Unexpected data size when reading a V3D NXS field '" + name + "'. Expected 3.");
428 std::copy(data.cbegin(), data.cend(), m_pt.begin());
429}
430
431double V3D::toMillerIndexes(double eps) {
432 if (eps < 0)
433 eps = -eps;
434 if (eps < FLT_EPSILON)
435 eps = FLT_EPSILON;
436
437 // assuming eps is in 1.e-x form
438
439 double ax = std::abs(m_pt[0]);
440 double ay = std::abs(m_pt[1]);
441 double az = std::abs(m_pt[2]);
442
443 double amax = (ax > ay) ? ax : ay;
444 amax = (az > amax) ? az : amax;
445 if (amax < FLT_EPSILON)
446 throw(std::invalid_argument("vector length is less then accuracy requested"));
447
448 if (ax < eps) {
449 m_pt[0] = 0;
450 ax = 0;
451 }
452 if (ay < eps) {
453 m_pt[1] = 0;
454 ay = 0;
455 }
456 if (az < eps) {
457 m_pt[2] = 0;
458 az = 0;
459 }
460
461 double mult(1);
462 mult = nearInt(ax, eps, mult);
463 mult = nearInt(ay, eps, mult);
464 mult = nearInt(az, eps, mult);
465
466 const size_t iax = std::lround(ax * mult / eps);
467 const size_t iay = std::lround(ay * mult / eps);
468 const size_t iaz = std::lround(az * mult / eps);
469
470 const size_t div = gcd(iax, gcd(iay, iaz));
471 mult /= (static_cast<double>(div) * eps);
472 this->operator*=(mult);
473
474 return mult;
475}
476
483bool V3D::compareMagnitude(const V3D &v1, const V3D &v2) {
484 const double mag_sq_1 = v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2];
485 const double mag_sq_2 = v2[0] * v2[0] + v2[1] * v2[1] + v2[2] * v2[2];
486 return (mag_sq_1 < mag_sq_2);
487}
488
495V3D V3D::directionAngles(bool inDegrees) const {
496 const double conversionFactor = inDegrees ? 180. / M_PI : 1.0;
497 const double divisor = this->norm();
498 if (divisor == 0.) {
499 throw std::runtime_error("Cannot calculate direction angles for zero length vector");
500 }
501 return V3D(conversionFactor * acos(m_pt[0] / divisor), conversionFactor * acos(m_pt[1] / divisor),
502 conversionFactor * acos(m_pt[2] / divisor));
503}
504
510 int MaxOrder = 0;
511 if (abs(static_cast<int>(m_pt[0])) > MaxOrder)
512 MaxOrder = abs(static_cast<int>(m_pt[0]));
513 if (abs(static_cast<int>(m_pt[1])) > MaxOrder)
514 MaxOrder = abs(static_cast<int>(m_pt[1]));
515 if (abs(static_cast<int>(m_pt[2])) > MaxOrder)
516 MaxOrder = abs(static_cast<int>(m_pt[2]));
517 return MaxOrder;
518}
519
524V3D V3D::absoluteValue() const { return V3D(fabs(m_pt[0]), fabs(m_pt[1]), fabs(m_pt[2])); }
525
530double V3D::hklError() const {
531 return fabs(m_pt[0] - std::round(m_pt[0])) + fabs(m_pt[1] - std::round(m_pt[1])) +
532 fabs(m_pt[2] - std::round(m_pt[2]));
533}
534
535} // namespace Mantid::Kernel
#define fabs(x)
Definition: Matrix.cpp:22
double tolerance
Numerical Matrix class.
Definition: Matrix.h:42
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
Class for quaternions.
Definition: Quat.h:39
void rotate(V3D &) const
Rotate a vector.
Definition: Quat.cpp:397
void setRotation(const double deg)
Set the rotation (both don't change rotation axis)
Definition: Quat.cpp:158
Class for 3D vectors.
Definition: V3D.h:34
void loadNexus(::NeXus::File *file, const std::string &name)
Load the object from an open NeXus file.
Definition: V3D.cpp:423
void spherical_rad(const double R, const double polar, const double azimuth) noexcept
Sets the vector position based on spherical coordinates, in radians.
Definition: V3D.cpp:70
std::array< double, 3 > m_pt
Definition: V3D.h:329
static bool compareMagnitude(const Kernel::V3D &v1, const Kernel::V3D &v2)
Convenience method for sorting list of V3D objects based on magnitude.
Definition: V3D.cpp:483
bool coLinear(const V3D &, const V3D &) const noexcept
Determines if this,B,C are colinear.
Definition: V3D.cpp:230
static std::vector< V3D > makeVectorsOrthogonal(const std::vector< V3D > &vectors)
Take a list of 2 vectors and makes a 3D orthogonal system out of them The first vector i0 is taken as...
Definition: V3D.cpp:298
void azimuth_polar_SNS(const double R, const double azimuth, const double polar) noexcept
Sets the vector position based on azimuth and polar angle, in RADIANS, in the SNS instrument coordina...
Definition: V3D.cpp:95
V3D absoluteValue() const
Absolute value.
Definition: V3D.cpp:524
double zenith(const V3D &) const noexcept
Zenith (theta) angle between this and another vector.
Definition: V3D.cpp:150
constexpr double scalar_prod(const V3D &v) const noexcept
Calculates the cross product.
Definition: V3D.h:274
void printSelf(std::ostream &) const
Prints a text representation of itself in format "[x,y,z]".
Definition: V3D.cpp:357
V3D & operator*=(const V3D &v) noexcept
Self-Inner product.
Definition: V3D.h:110
V3D & operator/=(const V3D &v) noexcept
Self-Inner division.
Definition: V3D.h:122
double normalize()
Make a normalized vector (return norm value)
Definition: V3D.cpp:130
constexpr V3D cross_prod(const V3D &v) const noexcept
Cross product (this * argument)
Definition: V3D.h:278
std::string toString() const
Definition: V3D.cpp:340
int maxCoeff()
Maximum absolute integer value.
Definition: V3D.cpp:509
V3D directionAngles(bool inDegrees=true) const
Direction angles.
Definition: V3D.cpp:495
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 read(std::istream &)
Read data from a stream.
Definition: V3D.cpp:328
double angle(const V3D &) const
Angle between this and another vector.
Definition: V3D.cpp:165
void rotate(const Matrix< double > &) noexcept
Rotate a point by a matrix.
Definition: V3D.cpp:217
double norm() const noexcept
Definition: V3D.h:263
int reBase(const V3D &, const V3D &, const V3D &) noexcept
rebase to new basis vector
Definition: V3D.cpp:199
void fromString(const std::string &str)
Sets the vector using a string.
Definition: V3D.cpp:348
bool unitVector(const double tolerance=Kernel::Tolerance) const noexcept
Definition: V3D.cpp:251
double hklError() const
Calculates the error in hkl.
Definition: V3D.cpp:530
void readPrinted(std::istream &)
Read data from a stream in the format returned by printSelf ("[x,y,z]").
Definition: V3D.cpp:364
void write(std::ostream &) const
Write out the point values.
Definition: V3D.cpp:330
void getSpherical(double &R, double &theta, double &phi) const noexcept
Return the vector's position in spherical coordinates.
Definition: V3D.cpp:117
double toMillerIndexes(double eps=1.e-3)
transform vector into form, used to describe directions in crystallogaphical coodinate system
Definition: V3D.cpp:431
double cosAngle(const V3D &) const
cos(Angle) between this and another vector
Definition: V3D.cpp:180
void round() noexcept
Round each component to the nearest integer.
Definition: V3D.cpp:140
int masterDir(const double Tol=1e-3) const noexcept
Determine if there is a master direction.
Definition: V3D.cpp:264
void saveNexus(::NeXus::File *file, const std::string &name) const
Save the object to an open NeXus file.
Definition: V3D.cpp:412
constexpr V3D() noexcept
Definition: V3D.h:36
bool nullVector(const double tolerance=1e-3) const noexcept
Determine if the point is null.
Definition: V3D.cpp:241
constexpr double rad2deg
constexpr double deg2rad
Defines units/enum for Crystal work.
Definition: AngleUnits.h:20
MANTID_KERNEL_DLL std::ostream & operator<<(std::ostream &, CPUTimer &)
Convenience function to provide for easier debug printing.
Definition: CPUTimer.cpp:86
MANTID_KERNEL_DLL std::istream & operator>>(std::istream &, Interpolation &)
Reads in parameter value.
constexpr double Tolerance
Standard tolerance value.
Definition: Tolerance.h:12
MANTID_KERNEL_DLL V3D normalize(V3D v)
Normalizes a V3D.
Definition: V3D.h:341