Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

Jiaming Guo
2 parents eceb6571 e702aa4e
Showing 4 changed files with 42 additions and 6 deletions Show diff stats
stim/cuda/templates/conv2sep.cuh
stim/envi/agilent_binary.h
stim/image/image.h
stim/math/spharmonics.h
@@ -8,8 +8,6 @@
 #include <stim/cuda/sharedmem.cuh>
 #include <stim/cuda/cudatools/error.h>
-#define pi	3.14159
-
 namespace stim{
 	namespace cuda{
@@ -125,7 +125,7 @@ public:
 		infile.read((char*)ptr, bytes());				//read the data		
 		infile.close();
 		Z[0] = 1;
-		Z[1] = R[2];
+		Z[1] = (double)R[2];
 	}
 	void save(std::string filename){
@@ -480,7 +480,6 @@ public:
 				max_val = img[n];
 			}
 		}
-
 		return max_val;
 	}
@@ -508,6 +507,33 @@ public:
 		return r;								//return the inverted image
 	}
+	image<T> crop(size_t x0, size_t y0, size_t w, size_t h){
+		image<T> result(w, h, C());								//create the output cropped image
+
+		size_t srci;
+		size_t dsti;
+		size_t line_bytes = w * C();							//calculate the number of bytes in a line
+		for (size_t yi = 0; yi < h; yi++) {						//for each row in the cropped image
+			srci = (y0 + yi) * X() * C() + x0 * C();			//calculate the source index
+			dsti = yi * w * C();								//calculate the destination index
+			memcpy(&result.img[dsti], &img[srci], line_bytes);	//copy the data
+		}
+		return result;
+	}
+
+	//crop regions given by an array of 1D index values
+	std::vector<image<T>> crop_idx(size_t w, size_t h, std::vector<size_t> idx) {
+		std::vector<image<T>> result(idx.size());										//create an array of image files to return
+		for (size_t i = 0; i < idx.size(); i++) {										//for each specified index point
+			size_t y = idx[i] / X();													//calculate the y coordinate from the 1D index (center of ROI)
+			size_t x = idx[i] - y * X();												//calculate the x coordinate (center of ROI)
+			y -= w / 2;																	//update x and y values to reflect the lower corner of the ROI
+			x -= h / 2;
+			result[i] = crop(x, y, w, h);												//get the cropped image and store it in the result array
+		}
+		return result;
+	}
+
 	image<T> srgb2lab(){
 		std::cout<<"ERROR stim::image::srgb2lab - function has been broken, re-implement."<<std::endl;
 		exit(1);
@@ -25,8 +25,20 @@ protected:
 	//calculate the value of the SH basis function (l, m) at (theta, phi)
 		//here, theta = [0, PI], phi = [0, 2*PI]
 	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();
+		//std::complex<T> result = boost::math::spherical_harmonic(l, m, phi, theta);
+		//return result.imag() + result.real();
+
+		//this calculation is based on calculating the real spherical harmonics:
+		//		https://en.wikipedia.org/wiki/Spherical_harmonics#Addition_theorem
+		if (m < 0) {
+			return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, abs(m), phi, theta).imag();
+		}
+		else if (m == 0) {
+			return boost::math::spherical_harmonic(l, m, phi, theta).real();
+		}
+		else {
+			return sqrt(2.0) * pow(-1, m) * boost::math::spherical_harmonic(l, m, phi, theta).real();
+		}
 	}
 	unsigned int coeff_1d(unsigned int l, int m){