From bfaba5dba7595086af2cc61e241da7a8589bc200 Mon Sep 17 00:00:00 2001
From: bkaradzic <branimirkaradzic@gmail.com>
Date: Wed, 30 Oct 2013 20:18:08 -0700
Subject: [PATCH] Added generating random point on sphere and hemisphere.

---
 include/bx/rng.h | 72 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 72 insertions(+)

diff --git a/include/bx/rng.h b/include/bx/rng.h
index db819cb..9889e49 100644
--- a/include/bx/rng.h
+++ b/include/bx/rng.h
@@ -8,6 +8,9 @@
 
 #include "bx.h"
 
+#define _USE_MATH_DEFINES
+#include <math.h>
+
 namespace bx
 {
 	// George Marsaglia's MWC
@@ -107,6 +110,75 @@ namespace bx
 		return 2.0f * bx::frnd(_rng) - 1.0f;
 	}
 
+	/// Generate random point on unit sphere.
+	template <typename Ty>
+	static inline void randUnitSphere(float _result[3], Ty* _rng)
+	{
+		float rand0 = frnd(_rng) * 2.0f - 1.0f;
+		float rand1 = frnd(_rng) * float(M_PI) * 2.0f;
+
+		float sqrtf1 = sqrtf(1.0f - rand0*rand0);
+		_result[0] = sqrtf1 * cosf(rand1);
+		_result[1] = sqrtf1 * sinf(rand1);
+		_result[2] = rand0;
+	}
+
+	/// Generate random point on unit hemisphere.
+	template <typename Ty>
+	static inline void randUnitHemisphere(float _result[3], Ty* _rng, const float _normal[3])
+	{
+		float dir[3];
+		randUnitSphere(dir, _rng);
+
+		float DdotN = dir[0]*_normal[0]
+					+ dir[1]*_normal[1]
+					+ dir[2]*_normal[2]
+					;
+
+		if (0.0f > DdotN)
+		{
+			dir[0] = -dir[0];
+			dir[1] = -dir[1];
+			dir[2] = -dir[2];
+		}
+
+		_result[0] = dir[0];
+		_result[1] = dir[1];
+		_result[2] = dir[2];
+	}
+
+	/// Sampling with Hammersley and Halton Points
+	/// http://www.cse.cuhk.edu.hk/~ttwong/papers/udpoint/udpoints.html
+	///
+	static inline void generateSphereHammersley(void* _data, uint32_t _stride, uint32_t _num, float _scale = 1.0f)
+	{
+		uint8_t* data = (uint8_t*)_data;
+
+		for (uint32_t ii = 0; ii < _num; ii++)
+		{
+			float tt = 0.0f;
+			float pp = 0.5;
+			for (uint32_t jj = ii; jj; jj >>= 1)
+			{
+				tt += (jj & 1) ? pp : 0.0f;
+				pp *= 0.5f;
+			}
+
+			tt = 2.0f * tt - 1.0f;
+
+			const float phi = (ii + 0.5f) / _num;
+			const float phirad =  phi * 2.0f * float(M_PI);
+			const float st = sqrtf(1.0f-tt*tt) * _scale;
+
+			float* xyz = (float*)data;
+			data += _stride;
+
+			xyz[0] = st * cosf(phirad);
+			xyz[1] = st * sinf(phirad);
+			xyz[2] = tt * _scale;
+		}
+	}
+
 } // namespace bx
 
 #endif // __BX_RNG_H__