Hi,
I have implemented a function to rotate a 3D vector a given angle around a
given axis (basically wrapping the functionality provided by Eigen::AngleAxis)
as an Rcpp function.
Below is an extract from the source file:
#include <Rcpp.h>
#include <RcppEigen.h>
#include <Eigen/Eigen>
#include <Eigen/Geometry>
using namespace Rcpp;
using Rcpp::as;
using Eigen::Map;
using Eigen::VectorXd;
using Eigen::Vector3d;
using Eigen::Matrix;
using Eigen::MatrixXd;
using Eigen::Matrix3d;
using Eigen::Dynamic;
using Eigen::AngleAxisd;
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
NumericVector rotate3DVectorAngleAxis(NumericVector x, NumericVector axis,
double angle) {
// Note: x should be 1 single vector to rotate
Map<VectorXd> xEigen(as<Map<VectorXd> >(x));
Map<VectorXd> axisEigen(as<Map<VectorXd> >(axis));
Matrix3d rotation;
rotation = AngleAxisd(angle, axisEigen.normalized());
Vector3d rotatedVector;
rotatedVector = rotation.matrix() * xEigen;
return wrap(rotatedVector);
}
The function executes with no problems and gives the expected output.
However, I have noticed that there is times where a value of 0 would be
expected for some vector components, but instead a very small value is
returned. For example, executing the following on my machine:
rotate3DVectorAngleAxis(c(0,0,1), c(1,0,0), pi/2)
Which represents a rotation of 90 degrees around the X axis, and therefore the
expected output value would be c(0,1,0)
However, it instead results in an output vector of c(0.000000e+00,
-1.000000e+00, 6.123234e-17)
I.e., the value for the Z component is 6.123234e-17.
I guess this is somehow related to machine precision, but is there an exact
cause that could be possibly fixed?
Additionally, why does the deviation from 0 only seem to happen for the Z
component, and not for the X component?
Varying the amount of rotation also seems to lead to a cumulative error on the
same component. E.g., rotating by 4001 times pi/2 (i.e., 1000 complete 360
degrees rotations plus a 90 degree rotation):
rotate3DVectorAngleAxis(c(0,0,1), c(1,0,0), 14401*pi/2)
results in c(0.000000e+00, -1.000000e+00, 4.776933e-13)
Which seems again to be accumulating in the same Z component.
Is there anything I could improve in my implementation to avoid these numerical
errors?
For reference, .Machine$double.eps on my machine is 2.220446e-16
Thanks a lot in advance
Best wishes,
Rafa
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
