On 10 January 2024 at 09:52, Rafael Ayala Hernandez wrote:
| 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
>From the looks it seems like an instance of R FAQ 7.31:
https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f
and the 'classic' referenced therein
https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html
It is probably your responsibility to identify such values and explicitly
round or set them to zero. I know nothing about the rotation here but in
other field (esp with iterations) a common rule of thumb is 'less than
sqrt(eps) is zero' meaning 1e-8. Of course that is too general and may be
too broad as "it always depends".
Hope this helps a little. Maybe contact some field experts too.
Cheers, Dirk
| 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
--
dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org
_______________________________________________
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
