Fixed the pdf implementation in spharmonics classes, added the random.h class in…

… stim/math for access to randomized variables and to sample a solid sphere (or angle), the functions are static but an instance can be initialized if a specific seed is needed

Fixed the pdf implementation in spharmonics classes, added the random.h class in…
… stim/math for access to randomized variables and to sample a solid sphere (or angle), the functions are static but an instance can be initialized if a specific seed is needed
Pavel Govyadinov
1 parent 9bfae4d0
Showing 2 changed files with 152 additions and 2 deletions Show diff stats
stim/math/random.h
stim/math/spharmonics.h
+#ifndef STIM_RANDOM
+#define STIM_RANDOM
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include <stim/math/vec3.h>
+#include <stim/math/constants.h>
+
+namespace stim{
+
+template<class T>
+class Random{
+protected:
+	void init(int seed = 0)
+	{
+		if(seed <= 0)
+			srand(time(NULL));
+		else if(seed > 0)
+			srand(time(seed));
+		else
+			std::cout << "Error: Unknown value: in STIM::RANDOM"  << std::endl;
+	}
+
+public:
+	/// Default Constructor
+	Random(){
+		init(-1);
+	}
+
+	/// Constructor from a seed.
+	/// A positive seed sets, 0 or negative yeilds the 
+	Random(int seed){
+		init(seed);
+	}
+
+	///Returns a random number uniformly sampled between 0 and 1
+	static
+	T uniformRandom()
+	{
+	     return (  (T)(rand()))/(  (T)(RAND_MAX));  ///generates a random number between 0 and 1 using the uniform distribution.
+	}
+
+	///Returns a random number from a normal distribution between 0 to 1.
+	static
+	T normalRandom()
+	{
+		T u1 = uniformRandom();
+		T u2 = uniformRandom();
+		return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1));                                                             ///generate a random number using the normal distribution between      0 and 1.
+	}
+	///Return a random vec3 each value between 0 and 1 from a uniform distribution.
+	static
+	stim::vec3<T> uniformRandVector()
+	{
+		stim::vec3<T> r(uniformRandom(), uniformRandom(), 1.0);                                                  ///generate a random vector using the uniform distribution between 0 and 1.
+		return r;
+	}
+	///Return a random vec3, each value between 0 and 1 from a normal distribution.
+	static
+	stim::vec3<T> normalRandVector()
+	{
+		stim::vec3<float> r(normalRandom(), normalRandom(), 1.0);                                                    ///generate a random vector using the normal distribution between      0 and 1.
+		return r;
+	}
+
+	///place num_samples of samples on the surface of a sphere of radius r.
+	///returns an std::vector of vec3's in cartisian coordinates.
+	static std::vector<stim::vec3 <T> >
+	sample_sphere(unsigned int num_samples, T radius = 1, T solidAngle = stim::TAU)
+	{
+		std::cout << "did this" << std::endl;
+		T PHI[2], Z[2], range;      ///Range of angles in cylinderical coordinates
+		PHI[0] = solidAngle/2;          ///project the solid angle into spherical coords
+		PHI[1] = asin(0);               ///
+		Z[0] = cos(PHI[0]);             ///project the z into spherical coordinates
+		Z[1] = cos(PHI[1]);             ///
+		range = Z[0] - Z[1];            ///the range of all possible z values.
+
+		T z, theta, phi;            /// temporary individual
+
+		std::vector<stim::vec3<T> > samples;
+
+		//srand(100);                     ///set random seed
+
+		for(int i = 0; i < num_samples; i++)		///for each sample
+		{
+			z = uniformRandom()*range + Z[1];	///find a random z based on the solid angle
+			theta = uniformRandom() * stim::TAU;	///find theta
+			phi = acos(z);				///project into spherical coord phi
+			stim::vec3<T> sph(radius, theta, phi);	///assume spherical
+			stim::vec3<T> cart = sph.sph2cart();	///conver to cartesisn
+			samples.push_back(cart);		///push into list
+		}
+
+///THIS IS DEBUGGING CODE, UNCOMMENT TO CHECK WHETHER THE SURFACE IS WELL SAMPLED!
+/*
+		std::stringstream name;
+	      	for(int i = 0; i < num_samples; i++)
+			name << samples[i].str() << std::endl;
+		
+	    
+	      	std::ofstream outFile;
+	      	outFile.open("Sampled Surface.txt");
+	      	outFile << name.str().c_str();                                                              
+*/
+
+		return samples;					///return full list.
+	}
+};
+
+}
+
+#endif
@@ -5,6 +5,7 @@
 #include <stim/math/vector.h>
 #include <boost/math/special_functions/spherical_harmonic.hpp>
 #include <stim/math/constants.h>
+#include <stim/math/random.h>
 #include <vector>
  
 #define WIRE_SCALE 1.001
@@ -26,7 +27,6 @@ protected:
 	double SH(int l, int m, double theta, double phi){
 		std::complex<T> result = boost::math::spherical_harmonic(l, m, phi, theta);
 		return result.imag() + result.real();
-//		return boost::math::spherical_harmonic_i(l, m, phi, theta);
 	}
  
 	unsigned int coeff_1d(unsigned int l, int m){
@@ -117,10 +117,10 @@ public:
 	}
  
 	/// Generates a PDF describing the probability distribution of points on a spherical surface
-
 	/// @param sph_pts is a list of points in spherical coordinates (theta, phi) where theta = [0, 2pi] and phi = [0, pi]
 	/// @param l is the maximum degree of the spherical harmonic function
 	/// @param m is the maximum order
+
 	void pdf(std::vector<stim::vec<double> > sph_pts, unsigned int l, int m){
  
 		mcBegin( l, m );		//begin spherical harmonic sampling
@@ -134,6 +134,42 @@ public:
 		mcEnd();
 	}
  
+	/// Generates a PDF describing the probability distribution of points on a spherical surface
+	/// @param sph_pts is a list of points in cartesian coordinates 
+	/// @param l is the maximum degree of the spherical harmonic function
+	/// @param m is the maximum order
+	/// @param c is the centroid of the points in sph_pts. DEFAULT 0,0,0
+	/// @param n is the number of points of the surface of the sphere used to create the PDF. DEFAULT 1000
+	void pdf(std::vector<stim::vec3<double> > sph_pts, unsigned int l, int m, stim::vec3<double> c = stim::vec3<double>(0,0,0), unsigned int n = 1000)
+	{
+		std::vector<double> weights;		///the weight at each point on the surface of the sphere.
+//		weights.resize(n);
+		unsigned int nP = sph_pts.size();
+		std::vector<stim::vec3<double> > sphere = stim::Random<double>::sample_sphere(n, 1.0, stim::TAU);
+		for(int i = 0; i < n; i++)
+		{
+			double val = 0;
+			for(int j = 0; j < nP; j++)
+			{
+				stim::vec3<double> temp = sph_pts[j] - c;
+				if(temp.dot(sphere[i]) > 0)
+					val += pow(temp.dot(sphere[i]),4);
+			}
+			weights.push_back(val);
+		}
+		
+		
+		mcBegin(l, m);		//begin spherical harmonic sampling
+		
+		for(unsigned int i = 0; i < n; i++)
+		{
+			stim::vec3<double> sph = sphere[i].cart2sph();
+			mcSample(sph[1], sph[2], weights[i]);
+		}
+
+		mcEnd();
+	}
+
 	std::string str(){
  
 		std::stringstream ss;