From f368a6dd500d7d3e0a5d5c4be1a89382f342af4e Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Wed, 2 Nov 2022 14:23:26 +0100
Subject: [PATCH] implemented uniformly random rotation quaternion

---
 quaternion.cc | 27 +++++++++++++++++++++++++++
 quaternion.h  |  1 +
 2 files changed, 28 insertions(+)

diff --git a/quaternion.cc b/quaternion.cc
index 2e6f5db..26fee64 100644
--- a/quaternion.cc
+++ b/quaternion.cc
@@ -408,6 +408,33 @@ r *= ::pow(xnorm,y);
 return r;
 }
 
+
+#ifndef QUAT_NO_RANDOM
+//method to efficiently generate a uniformly random point on a unit 4-sphere by G. Marsaglia, Ann. Math. Stat. 43, 645 (1972)
+template<typename T>
+void Quaternion<T>::random_rotation()
+{
+T  s1,s2,s;
+do
+	{
+	q[0]=2.*random()/(1. + RAND_MAX) - 1.;
+	q[1]=2.*random()/(1. + RAND_MAX) - 1.;
+	s1 = q[0]*q[0] + q[1]*q[1];
+	}
+	while(s1>1);
+do
+        {
+        q[2]=2.*random()/(1. + RAND_MAX) - 1.;
+        q[3]=2.*random()/(1. + RAND_MAX) - 1.;
+        s2 = q[2]*q[2] + q[3]*q[3];
+        }
+        while(s2>1);
+s = sqrt((1-s1)/s2);
+q[2] *= s;
+q[3] *= s;
+}
+#endif
+
 //force instantization
 #define INSTANTIZE(T) \
 template class Quaternion<T>; \
diff --git a/quaternion.h b/quaternion.h
index 08f80bc..bdc674e 100644
--- a/quaternion.h
+++ b/quaternion.h
@@ -124,6 +124,7 @@ public:
 	void axis2normquat(const T *axis, const T &angle);
 	void normquat2axis(T *axis,  T &angle) const;
 	
+	void random_rotation(); //generate uniformly random unit quaternion
 
 	//C-style IO
 	int fprintf(FILE *f, const char *format) const {return ::fprintf(f,format,q[0],q[1],q[2],q[3]);};