Open stla opened 6 years ago
RobinHankin
commented
6 years ago Hello stla, thanks for this! You're absolutely right, this would be a Good Thing to do, and I've never got round to it. But your functions look great and I'll include them in the onion package when I get a minute.
Best wishes, Robin
RobinHankin
commented
6 years ago Hi, thanks for this! I've replied on the github issues page [I am new to github, still feeling my way around]. I'll include your functionality in the package when I get a minute. Who should I credit?'
best wishes
Robin hankin.robin@gmail.com hankin.robin@gmail.com
On Fri, Jul 27, 2018 at 9:04 PM stla notifications@github.com wrote:
Hello Robin,
It would be a nice functionality if the package dealt with the relationship between quaternions and rotation matrices in 3D. Below are some functions of mine.
' @description Quaternion to rotation matrix
' @details http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/index.htm
quaternion2matrix <- function(q){ rbind( c(1 - 2q[3]^2 - 2q[4]^2, 2q[2]q[3] - 2q[4]q[1], 2q[2]q[4] + 2q[3]q[1]), c(2q[2]q[3] + 2q[4]q[1], 1 - 2q[2]^2 - 2q[4]^2, 2q[3]q[4] - 2q[2]q[1]), c(2q[2]q[4] - 2q[3]q[1], 2q[3]q[4] + 2q[2]q[1], 1 - 2q[2]^2 - 2q[3]^2) ) }
' @description Rotation matrix to quaternion
' @details http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm
matrix2quaternion <- function(M){ tr <- M[1,1] + M[2,2] + M[3,3] if (tr > 0) { S <- sqrt(tr+1.0) 2 qw <- 0.25 S qx <- (M[3,2] - M[2,3]) / S qy <- (M[1,3] - M[3,1]) / S qz <- (M[2,1] - M[1,2]) / S } else if ((M[1,1] > M[2,2]) && (M[1,1] > M[3,3])) { S <- sqrt(1.0 + M[1,1] - M[2,2] - M[3,3]) 2 qw <- (M[3,2] - M[2,3]) / S qx <- 0.25 S qy <- (M[1,2] + M[2,1]) / S qz <- (M[1,3] + M[3,1]) / S } else if (M[2,2] > M[3,3]) { S <- sqrt(1.0 + M[2,2] - M[1,1] - M[3,3]) 2 qw <- (M[1,3] - M[3,1]) / S qx <- (M[1,2] + M[2,1]) / S qy <- 0.25 S qz <- (M[2,3] + M[3,2]) / S } else { S <- sqrt(1.0 + M[3,3] - M[1,1] - M[2,2]) 2 qw <- (M[2,1] - M[1,2]) / S qx <- (M[1,3] + M[3,1]) / S qy <- (M[2,3] + M[3,2]) / S qz <- 0.25 S } return(c(qw,qx,qy,qz)) }
' @description Rotation matrix from axis and angle
axisAngle2matrix <- function(u, theta){ u <- u/sqrt(c(crossprod(u))) M1 <- cos(theta) diag(3) M2 <- sin(theta) rbind( c(0, -u[3], u[2]), c(u[3], 0, -u[1]), c(-u[2], u[1], 0) ) M3 <- (1-cos(theta)) * kronecker(u,t(u)) M1 + M2 + M3 }
' @description Quaternion from axis and angle
axisAngle2quaternion <- function(u, theta){ u <- u/sqrt(c(crossprod(u))) sine <- sin(theta/2) c(cos(theta/2), u[1]sine, u[2]sine, u[3]*sine) }
' @description Axis and angle from quaternion
quaternion2axisAngle <- function(q){ den <- sqrt(1-q[1]q[1]) list(axis = q[2:4]/den, angle = 2acos(q[1])) }
' @description Get the rotation transforming U to V
getRotation <- function(U,V){ ma <- sqrt(c(crossprod(U))) mb <- sqrt(c(crossprod(V))) # ?? no sense if ma /= mb ... d <- c(tcrossprod(U, t(V))) c <- sqrt(mamb+d) ma2 <- sqrt(2)ma r <- 1/ma2/c W <- geometry::extprod3d(U,V) c(c/ma2, r*W) }
tests
u <- c(1,1,1); theta <- 1 axisAngle2quaternion(u,theta) ( q <- matrix2quaternion(axisAngle2matrix(u,theta)) ) RSpincalc::DCM2Q(axisAngle2matrix(u,theta)) quaternion2axisAngle(q)
U <- c(1,0,0) V <- c(1,1,1)/sqrt(3) q <- getRotation(U,V) quaternion2matrix(q) %*% U
— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/RobinHankin/onion/issues/1, or mute the thread https://github.com/notifications/unsubscribe-auth/AMpc0vqjcgc37m6qiGILN_spskzzEYhDks5uKteIgaJpZM4VjKPg .
stla
commented
6 years ago Hi Robin, Put my name if you want but don't feel obliged. I mostly copied the formulas from http://www.euclideanspace.com/maths/geometry/rotations/. The formula for getRotation is copied from a Haskell package, but I'm not sure this function covers all the cases. Please check the code. For example I don't remember whether the quaternion must be normalized in quaternion2matrix, and this is not checked by the code. Best,Stéphane
Le vendredi 27 juillet 2018 à 11:26:01 UTC+2, Robin Hankin <notifications@github.com> a écrit : Hi, thanks for this! I've replied on the github issues page [I am new to github, still feeling my way around]. I'll include your functionality in the package when I get a minute. Who should I credit?'
best wishes
Robin hankin.robin@gmail.com hankin.robin@gmail.com
On Fri, Jul 27, 2018 at 9:04 PM stla notifications@github.com wrote:
Hello Robin,
It would be a nice functionality if the package dealt with the relationship between quaternions and rotation matrices in 3D. Below are some functions of mine.
' @description Quaternion to rotation matrix
' @details http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/index.htm
quaternion2matrix <- function(q){ rbind( c(1 - 2q[3]^2 - 2q[4]^2, 2q[2]q[3] - 2q[4]q[1], 2q[2]q[4] + 2q[3]q[1]), c(2q[2]q[3] + 2q[4]q[1], 1 - 2q[2]^2 - 2q[4]^2, 2q[3]q[4] - 2q[2]q[1]), c(2q[2]q[4] - 2q[3]q[1], 2q[3]q[4] + 2q[2]q[1], 1 - 2q[2]^2 - 2q[3]^2) ) }
' @description Rotation matrix to quaternion
' @details http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm
matrix2quaternion <- function(M){ tr <- M[1,1] + M[2,2] + M[3,3] if (tr > 0) { S <- sqrt(tr+1.0) 2 qw <- 0.25 S qx <- (M[3,2] - M[2,3]) / S qy <- (M[1,3] - M[3,1]) / S qz <- (M[2,1] - M[1,2]) / S } else if ((M[1,1] > M[2,2]) && (M[1,1] > M[3,3])) { S <- sqrt(1.0 + M[1,1] - M[2,2] - M[3,3]) 2 qw <- (M[3,2] - M[2,3]) / S qx <- 0.25 S qy <- (M[1,2] + M[2,1]) / S qz <- (M[1,3] + M[3,1]) / S } else if (M[2,2] > M[3,3]) { S <- sqrt(1.0 + M[2,2] - M[1,1] - M[3,3]) 2 qw <- (M[1,3] - M[3,1]) / S qx <- (M[1,2] + M[2,1]) / S qy <- 0.25 S qz <- (M[2,3] + M[3,2]) / S } else { S <- sqrt(1.0 + M[3,3] - M[1,1] - M[2,2]) 2 qw <- (M[2,1] - M[1,2]) / S qx <- (M[1,3] + M[3,1]) / S qy <- (M[2,3] + M[3,2]) / S qz <- 0.25 S } return(c(qw,qx,qy,qz)) }
' @description Rotation matrix from axis and angle
axisAngle2matrix <- function(u, theta){ u <- u/sqrt(c(crossprod(u))) M1 <- cos(theta) diag(3) M2 <- sin(theta) rbind( c(0, -u[3], u[2]), c(u[3], 0, -u[1]), c(-u[2], u[1], 0) ) M3 <- (1-cos(theta)) * kronecker(u,t(u)) M1 + M2 + M3 }
' @description Quaternion from axis and angle
axisAngle2quaternion <- function(u, theta){ u <- u/sqrt(c(crossprod(u))) sine <- sin(theta/2) c(cos(theta/2), u[1]sine, u[2]sine, u[3]*sine) }
' @description Axis and angle from quaternion
quaternion2axisAngle <- function(q){ den <- sqrt(1-q[1]q[1]) list(axis = q[2:4]/den, angle = 2acos(q[1])) }
' @description Get the rotation transforming U to V
getRotation <- function(U,V){ ma <- sqrt(c(crossprod(U))) mb <- sqrt(c(crossprod(V))) # ?? no sense if ma /= mb ... d <- c(tcrossprod(U, t(V))) c <- sqrt(mamb+d) ma2 <- sqrt(2)ma r <- 1/ma2/c W <- geometry::extprod3d(U,V) c(c/ma2, r*W) }
tests
u <- c(1,1,1); theta <- 1 axisAngle2quaternion(u,theta) ( q <- matrix2quaternion(axisAngle2matrix(u,theta)) ) RSpincalc::DCM2Q(axisAngle2matrix(u,theta)) quaternion2axisAngle(q)
U <- c(1,0,0) V <- c(1,1,1)/sqrt(3) q <- getRotation(U,V) quaternion2matrix(q) %*% U
— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/RobinHankin/onion/issues/1, or mute the thread https://github.com/notifications/unsubscribe-auth/AMpc0vqjcgc37m6qiGILN_spskzzEYhDks5uKteIgaJpZM4VjKPg .
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.
RobinHankin
commented
6 years ago thanks for this, I hadn't connected the permutations bug-finder with the quaternion rotation writer!
Best wishes
Robin
hankin.robin@gmail.com hankin.robin@gmail.com
On Sat, Jul 28, 2018 at 1:59 AM stla notifications@github.com wrote:
Hi Robin, Put my name if you want but don't feel obliged. I mostly copied the formulas from http://www.euclideanspace.com/maths/geometry/rotations/. The formula for getRotation is copied from a Haskell package, but I'm not sure this function covers all the cases. Please check the code. For example I don't remember whether the quaternion must be normalized in quaternion2matrix, and this is not checked by the code. Best,Stéphane
Le vendredi 27 juillet 2018 à 11:26:01 UTC+2, Robin Hankin < notifications@github.com> a écrit :
Hi, thanks for this! I've replied on the github issues page [I am new to github, still feeling my way around]. I'll include your functionality in the package when I get a minute. Who should I credit?'
best wishes
Robin hankin.robin@gmail.com hankin.robin@gmail.com
On Fri, Jul 27, 2018 at 9:04 PM stla notifications@github.com wrote:
Hello Robin,
It would be a nice functionality if the package dealt with the relationship between quaternions and rotation matrices in 3D. Below are some functions of mine.
' @description Quaternion to rotation matrix
' @details
http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/index.htm quaternion2matrix <- function(q){ rbind( c(1 - 2q[3]^2 - 2q[4]^2, 2q[2]q[3] - 2q[4]q[1], 2q[2]q[4] + 2q[3]q[1]), c(2q[2]q[3] + 2q[4]q[1], 1 - 2q[2]^2 - 2q[4]^2, 2q[3]q[4] - 2q[2]q[1]), c(2q[2]q[4] - 2q[3]q[1], 2q[3]q[4] + 2q[2]q[1], 1 - 2q[2]^2 - 2q[3]^2) ) }
' @description Rotation matrix to quaternion
' @details
http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/index.htm matrix2quaternion <- function(M){ tr <- M[1,1] + M[2,2] + M[3,3] if (tr > 0) { S <- sqrt(tr+1.0) 2 qw <- 0.25 S qx <- (M[3,2] - M[2,3]) / S qy <- (M[1,3] - M[3,1]) / S qz <- (M[2,1] - M[1,2]) / S } else if ((M[1,1] > M[2,2]) && (M[1,1] > M[3,3])) { S <- sqrt(1.0 + M[1,1] - M[2,2] - M[3,3]) 2 qw <- (M[3,2] - M[2,3]) / S qx <- 0.25 S qy <- (M[1,2] + M[2,1]) / S qz <- (M[1,3] + M[3,1]) / S } else if (M[2,2] > M[3,3]) { S <- sqrt(1.0 + M[2,2] - M[1,1] - M[3,3]) 2 qw <- (M[1,3] - M[3,1]) / S qx <- (M[1,2] + M[2,1]) / S qy <- 0.25 S qz <- (M[2,3] + M[3,2]) / S } else { S <- sqrt(1.0 + M[3,3] - M[1,1] - M[2,2]) 2 qw <- (M[2,1] - M[1,2]) / S qx <- (M[1,3] + M[3,1]) / S qy <- (M[2,3] + M[3,2]) / S qz <- 0.25 S } return(c(qw,qx,qy,qz)) }
' @description Rotation matrix from axis and angle
axisAngle2matrix <- function(u, theta){ u <- u/sqrt(c(crossprod(u))) M1 <- cos(theta) diag(3) M2 <- sin(theta) rbind( c(0, -u[3], u[2]), c(u[3], 0, -u[1]), c(-u[2], u[1], 0) ) M3 <- (1-cos(theta)) * kronecker(u,t(u)) M1 + M2 + M3 }
' @description Quaternion from axis and angle
axisAngle2quaternion <- function(u, theta){ u <- u/sqrt(c(crossprod(u))) sine <- sin(theta/2) c(cos(theta/2), u[1]sine, u[2]sine, u[3]*sine) }
' @description Axis and angle from quaternion
quaternion2axisAngle <- function(q){ den <- sqrt(1-q[1]q[1]) list(axis = q[2:4]/den, angle = 2acos(q[1])) }
' @description Get the rotation transforming U to V
getRotation <- function(U,V){ ma <- sqrt(c(crossprod(U))) mb <- sqrt(c(crossprod(V))) # ?? no sense if ma /= mb ... d <- c(tcrossprod(U, t(V))) c <- sqrt(mamb+d) ma2 <- sqrt(2)ma r <- 1/ma2/c W <- geometry::extprod3d(U,V) c(c/ma2, r*W) }
tests
u <- c(1,1,1); theta <- 1 axisAngle2quaternion(u,theta) ( q <- matrix2quaternion(axisAngle2matrix(u,theta)) ) RSpincalc::DCM2Q(axisAngle2matrix(u,theta)) quaternion2axisAngle(q)
U <- c(1,0,0) V <- c(1,1,1)/sqrt(3) q <- getRotation(U,V) quaternion2matrix(q) %*% U
— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/RobinHankin/onion/issues/1, or mute the thread < https://github.com/notifications/unsubscribe-auth/AMpc0vqjcgc37m6qiGILN_spskzzEYhDks5uKteIgaJpZM4VjKPg
.
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.
— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/RobinHankin/onion/issues/1#issuecomment-408427139, or mute the thread https://github.com/notifications/unsubscribe-auth/AMpc0kvhw4kij3qBItSGbF6XJxkLj-Itks5uKxzWgaJpZM4VjKPg .
stla
commented
6 years ago Hello Robin,
I have also implemented these functions with RcppEigen.
// we only include RcppEigen.h which pulls Rcpp.h in for us
#include <RcppEigen.h>
#include <Eigen/Geometry>
using namespace std;
// [[Rcpp::export]]
Rcpp::NumericVector multQuat(const float & w1, const float & w2,
const float & x1, const float & x2,
const float & y1, const float & y2,
const float & z1, const float & z2)
{
Eigen::Quaternionf q1(w1, x1, y1, z1);
Eigen::Quaternionf q2(w2, x2, y2, z2);
Eigen::Quaternionf q1q2 = q1*q2;
return Rcpp::NumericVector::create(
q1q2.w(), q1q2.x(), q1q2.y(), q1q2.z() );
}
// [[Rcpp::export]]
Rcpp::NumericVector squaredQuat(const float & w,
const float & x,
const float & y,
const float & z)
{
Eigen::Quaternionf q(w, x, y, z);
Eigen::Quaternionf qsquared = q*q;
return Rcpp::NumericVector::create(
qsquared.w(), qsquared.x(), qsquared.y(), qsquared.z() );
}
// [[Rcpp::export]]
Rcpp::NumericVector fromAxisAngle(const float & angle,
Eigen::VectorXf & axis)
{
axis.normalize();
Eigen::AngleAxisf AA(angle, axis);
Eigen::Quaternionf q(AA);
return Rcpp::NumericVector::create(
q.w(), q.x(), q.y(), q.z());
}
// [[Rcpp::export]]
Eigen::MatrixXf fromQuaternion(const float & w,
const float & x,
const float & y,
const float & z)
{
Eigen::Quaternionf q(w, x, y, z);
Eigen::Quaternionf qn = q.normalized();
//Eigen::Quaternionf (q.w(), q.x(), q.y(), q.z());
return qn.toRotationMatrix();
}
// [[Rcpp::export]]
Rcpp::NumericVector toQuaternion(Eigen::MatrixXf & rotmatrix)
{
float angle; Eigen::Vector3f axis;
Eigen::Quaternionf q; q = Eigen::AngleAxisf(angle, axis);
return Rcpp::NumericVector::create(q.w(), q.x(), q.y(), q.z());
}
// [[Rcpp::export]]
Rcpp::NumericVector getRotation(Eigen::VectorXf & a, Eigen::VectorXf & b)
{
Eigen::Quaternionf q;
Eigen::Quaternionf quat = q.setFromTwoVectors(a,b);
return Rcpp::NumericVector::create(quat.w(), quat.x(), quat.y(), quat.z());
}
Take a look at the benchmarks:
library(onion)
library(microbenchmark)
q1 <- quaternion(Re=1, i=2, j=3, k=4)
q2 <- quaternion(Re=10, i=12, j=13, k=14)> microbenchmark(
+ onion = q1*q2,
+ RcppEigen = multQuat(1,10,2,12,3,13,4,14)
+ )
Unit: microseconds
expr min lq mean median uq max neval cld
onion 352.089 357.4445 485.40968 384.442 407.6465 9399.279 100 b
RcppEigen 2.679 3.5710 5.82021 4.910 6.6945 17.405 100 a
hi again, thanks for these. Did you set the use.R flag in the onion evaluation? The package has the option to use C routines or R routines (I can't remember why now). Best wishes, Robin
stla
commented
6 years ago Sorry there were some mistakes. Below is the corrected code.
// we only include RcppEigen.h which pulls Rcpp.h in for us
#include <RcppEigen.h>
#include <Eigen/Geometry>
#include <Eigen/Core>
using namespace std;
// [[Rcpp::export]]
Rcpp::NumericVector multQuat(const float & w1, const float & w2,
const float & x1, const float & x2,
const float & y1, const float & y2,
const float & z1, const float & z2)
{
Eigen::Quaternionf q1(w1, x1, y1, z1);
Eigen::Quaternionf q2(w2, x2, y2, z2);
Eigen::Quaternionf q1q2 = q1*q2;
return Rcpp::NumericVector::create(
q1q2.w(), q1q2.x(), q1q2.y(), q1q2.z() );
}
// [[Rcpp::export]]
Rcpp::NumericVector squaredQuat(const float & w,
const float & x,
const float & y,
const float & z)
{
Eigen::Quaternionf q(w, x, y, z);
Eigen::Quaternionf qsquared = q*q;
return Rcpp::NumericVector::create(
qsquared.w(), qsquared.x(), qsquared.y(), qsquared.z() );
}
// [[Rcpp::export]]
Rcpp::NumericVector fromAxisAngle(const float & angle,
Eigen::VectorXf & axis)
{
axis.normalize();
Eigen::AngleAxisf AA(angle, axis);
Eigen::Quaternionf q(AA);
return Rcpp::NumericVector::create(
q.w(), q.x(), q.y(), q.z());
}
// [[Rcpp::export]]
Eigen::MatrixXf fromQuaternion(const float & w,
const float & x,
const float & y,
const float & z)
{
Eigen::Quaternionf q(w, x, y, z);
Eigen::Quaternionf qn = q.normalized();
//Eigen::Quaternionf (q.w(), q.x(), q.y(), q.z());
return qn.toRotationMatrix();
}
// [[Rcpp::export]]
Rcpp::List toAxisAngle(const Rcpp::NumericVector & v)
{
Eigen::Quaternionf q(v[0],v[1],v[2],v[3]);
Eigen::AngleAxisf AA(q);
Eigen::Vector3f AnAx = AA.axis();
float ang = AA.angle();
return Rcpp::List::create(Rcpp::Named("axis") = AnAx,
Rcpp::Named("angle") = ang);
}
// [[Rcpp::export]]
//Rcpp::NumericVector toQuaternion(Eigen::MatrixXf & rotmatrix)
//{
// Eigen::Quaternionf q;
// Eigen::AngleAxisf AnAX = q(rotmatrix);
// return fromAxisAngle(AnAx.angle(), AnaX.axis());
//}
// [[Rcpp::export]]
Rcpp::NumericVector getRotation(Eigen::VectorXf & a, Eigen::VectorXf & b)
{
Eigen::Quaternionf q;
Eigen::Quaternionf quat = q.setFromTwoVectors(a,b);
return Rcpp::NumericVector::create(quat.w(), quat.x(), quat.y(), quat.z());
}
ur=TRUE only slightly increases the speed:
> microbenchmark(
+ onion = AprodA(q1,q2, ur=TRUE),
+ RcppEigen = multQuat(1,10,2,12,3,13,4,14)
+ )
Unit: microseconds
expr min lq mean median uq max neval cld
onion 232.941 242.9825 523.18000 278.682 323.9755 22543.897 100 b
RcppEigen 2.679 3.5720 6.11484 4.910 7.5875 24.099 100 a
> microbenchmark(
+ onion = AprodA(q1,q2, ur=FALSE),
+ RcppEigen = multQuat(1,10,2,12,3,13,4,14)
+ )
Unit: microseconds
expr min lq mean median uq max neval cld
onion 340.487 363.469 546.98272 377.7485 442.0075 11040.571 100 b
RcppEigen 2.678 4.017 6.28438 5.3565 6.6950 17.405 100 a
no worries. I've made a start on implementing some of the quaternion-to-orthogonal matrix work (a new branch on the onion package repo), but it's gonna take me some time to finish. Also, you might be interested to know that I've been using the gyrogroup package to geneate "random" orthogonal matrices. Best wishes, Robin
stla
commented
6 years ago Sounds good. There's also the rstiefel package to generate random orthogonal matrices.
RobinHankin
commented
6 years ago I suspect the onion functionality is slow because it uses S3 and so operators have to be dispatched via Ops.onion(), and the package has to do housekeeping tasks such as getting and setting object classes and fiddling with extractor methods. But the basic default multiplication routine is written in C which ought to be pretty close to optimal.
But the package is vectorized. Could you try your benchmark on a quaternionic vector which is say 10^5 elements?
> library(onion)
> o1 <- rquat(1e5)
> o2 <- rquat(1e5)
> ignore <- o1*o2Indeed. My code is not vectorized and it is slower if I use vapply.
> library(onion)
> n <- 2000
> o1 <- rquat(n)
> o2 <- rquat(n)
>
> oo1 <- as.matrix(o1)
> oo2 <- as.matrix(o2)
>
> library(microbenchmark)
> microbenchmark(
+ onion = o1*o2,
+ Eigen = vapply(1:n, function(i) quatProd(oo1[,i],oo2[,i]), numeric(4))
+ )
Unit: microseconds
expr min lq mean median uq max neval cld
onion 801.459 958.3145 3817.184 1057.158 1231.193 242040.20 100 a
Eigen 8458.592 9160.5365 10973.875 10476.516 12435.535 17794.95 100 bI've been using the gyrogroup package to geneate "random" orthogonal matrices
FYI, I've found a method here to sample a versor (unit quaternion) such that the corresponding rotation matrix is uniformly sampled.
And the uniform generation of versors is implemented in (Rcpp)Eigen:
Rcpp::NumericVector rversor_(){
Eigen::Quaterniond q = Eigen::Quaterniond::UnitRandom();
return Rcpp::NumericVector::create(q.w(), q.x(), q.y(), q.z());
} For a vectorized sampling:
Rcpp::NumericMatrix rversors_(const unsigned n){
Rcpp::NumericMatrix out(n, 4);
for(unsigned i=0; i<n; i++){
Eigen::Quaterniond q = Eigen::Quaterniond::UnitRandom();
out(i,0) = q.w();
out(i,1) = q.x();
out(i,2) = q.y();
out(i,3) = q.z();
}
return out;
}
Hello Robin,
It would be a nice functionality if the package dealt with the relationship between quaternions and rotation matrices in 3D. Below are some functions of mine.