added median spectra to ENVI files, updated mean spectra in interleave formats

David Mayerich
1 parent 1b2858a2
Showing 9 changed files with 140 additions and 79 deletions Show diff stats
matlab/enviLoadRaw.m
matlab/enviSaveRaw.m
stim/cuda/cudatools/error.h
stim/envi/bil.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
stim/gl/gl_spider.h
stim/math/vector.h
 %loads an ENVI file without any manipulation (changing orientation)
+% enviLoadRaw(filename, headername)
 function M = enviLoadRaw(filename, headername)
     %if a header isn't provided, assume it's just the filename
 %saves an ENVI file without any manipulation, assumes (X, Y, S)
+% enviSaveRaw(M, filename, headername)
 function enviSaveRaw(M, filename, headername)
     %if a header isn't provided, assume it's just the filename
@@ -11,12 +11,7 @@ using namespace std;
 //handle error macro
 static void HandleError( cudaError_t err, const char *file,  int line ) {
    	if (err != cudaSuccess) {
-			//FILE* outfile = fopen("cudaErrorLog.txt", "w");
-      		//fprintf(outfile,  "%s in %s at line %d\n", cudaGetErrorString( err ),  file, line );
-			//fclose(outfile);
             printf("%s in %s at line %d\n", cudaGetErrorString( err ),  file, line );
-       		//exit( EXIT_FAILURE );
-
    	}
 }
 #define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))
@@ -992,32 +992,35 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
+	bool mean_spectrum(double* m, double* std, unsigned char* mask = NULL, bool PROGRESS = false){
 		unsigned long long XZ = X() * Z();
 		unsigned long long XY = X() * Y();
 		T* temp = (T*)malloc(sizeof(T) * XZ);
-		for (unsigned long long j = 0; j < Z(); j++){
-			p[j] = 0;
-		}
+		memset(m, 0, Z() * sizeof(double));							//initialize the mean to zero
+		double* e_x2 = (double*)malloc(Z() * sizeof(double));		//allocate space for E[x^2]
+		memset(e_x2, 0, Z() * sizeof(double));						//initialize E[x^2] to zero
 		//calculate vaild number in a band
-		unsigned long long count = 0;
-		for (unsigned long long j = 0; j < XY; j++){
-			if (mask == NULL || mask[j] != 0){
-				count++;
-			}
-		}
+		size_t count = nnz(mask);							//count the number of pixels in the mask
+
+		double x;											//create a register to store the pixel value
 		for (unsigned long long k = 0; k < Y(); k++){
 			read_plane_y(temp, k);
 			unsigned long long kx = k * X();
 			for (unsigned long long i = 0; i < X(); i++){
 				if (mask == NULL || mask[kx + i] != 0){
 					for (unsigned long long j = 0; j < Z(); j++){
-						p[j] += temp[j * X() + i] / (double)count;
+						x = temp[j * X() + i];
+						m[j] += x / (double)count;
+						e_x2[j] += x*x / (double)count;
 					}
 				}
 			}
 			if(PROGRESS) progress = (double)(k+1) / Y() * 100;
 		}
+
+		for(size_t i = 0; i < Z(); i++)							//calculate the standard deviation
+			std[i] = sqrt(e_x2[i] - m[i] * m[i]);
+
 		free(temp);
 		return true;
 	}
@@ -954,26 +954,31 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
+	bool mean_spectrum(double* m, double* std, unsigned char* mask = NULL, bool PROGRESS = false){
 		unsigned long long XY = X() * Y();							//calculate the total number of pixels in the HSI
 		T* temp = (T*)malloc(sizeof(T) * Z());						//allocate space for the current spectrum to be read
-		memset(p, 0, sizeof(double) * Z());							//initialize the average spectrum to zero (0)
-		//for (unsigned j = 0; j < Z(); j++){
-		//	p[j] = 0;
-		//}
+		memset(m, 0, Z() * sizeof(double));							//set the mean spectrum to zero
+		double* e_x2 = (double*)malloc(Z() * sizeof(double));		//allocate space for E[x^2]
+		memset(e_x2, 0, Z() * sizeof(double));						//set all values for E[x^2] to zero
 		unsigned long long count = nnz(mask);									//calculate the number of masked pixels
-
+		double x;
 		for (unsigned long long i = 0; i < XY; i++){							//for each pixel in the HSI
 			if (mask == NULL || mask[i] != 0){						//if the pixel is masked
 				pixel(temp, i);										//get the spectrum
 				for (unsigned long long j = 0; j < Z(); j++){					//for each spectral component
-					p[j] += (double)temp[j] / (double)count;		//add the weighted value to the average
+					x = temp[j];
+					m[j] += x / (double)count;		//add the weighted value to the average
+					e_x2[j] += x*x / (double)count;
 				}
 			}
 			if(PROGRESS) progress = (double)(i+1) / XY * 100;		//increment the progress
 		}
+		//calculate the standard deviation
+		for(size_t i = 0; i < Z(); i++)
+			std[i] = sqrt(e_x2[i] - m[i] * m[i]);
+
 		free(temp);
 		return true;
 	}
@@ -561,12 +561,12 @@ public:
 		free(src[1]);
 		free(dst[0]);
 		free(dst[1]);
-		//if(VERBOSE){
+		if(VERBOSE){
 			std::cout<<"total time to execute bsq::bip(): "<<t_total<<" ms"<<std::endl;
 			std::cout<<"     total time spent processing: "<<pt_total<<" ms"<<std::endl;
 			std::cout<<"        total time spent reading: "<<rt_total<<" ms"<<std::endl;
 			std::cout<<"        total time spent writing: "<<wt_total<<" ms"<<std::endl;
-		//}
+		}
 		return true;															//return true
 	}
@@ -1120,27 +1120,61 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool avg_band(double* p, unsigned char* mask = NULL, bool PROGRESS = false){
+	bool mean_spectrum(double* m, double* std, unsigned char* mask = NULL, bool PROGRESS = false){
 		unsigned long long XY = X() * Y();
-		unsigned long long count = 0;						//count will store the number of masked pixels
+		unsigned long long count = nnz(mask);						//count will store the number of masked pixels
 		T* temp = (T*)malloc(sizeof(T) * XY);
-		//calculate this loop counts the number of true pixels in the mask
-		for (unsigned j = 0; j < XY; j++){
-			if (mask == NULL || mask[j] != 0){
-				count++;
-			}
-		}
+		
 		//this loops goes through each band in B (Z())
 		//	masked (or valid) pixels from that band are averaged and the average is stored in p
+		double e_x;										//stores E[x]^2
+		double e_x2;										//stores E[x^2]
+		double x;
 		for (unsigned long long i = 0; i < Z(); i++){
-			p[i] = 0;
+			e_x = 0;
+			e_x2 = 0;
 			band_index(temp, i);				//get the band image and store it in temp
 			for (unsigned long long j = 0; j < XY; j++){	//loop through temp, averaging valid pixels
 				if (mask == NULL || mask[j] != 0){
-					p[i] += (double)temp[j] / (double)count;
+					x = (double)temp[j];
+					e_x += x / (double)count;				//sum the expected value of x
+					e_x2 += (x * x) / (double)count;		//sum the expected value of x^2
+				}
+			}
+			m[i] = e_x;												//store the mean
+			std[i] = sqrt(e_x2 - e_x * e_x);						//calculate the standard deviation
+			if(PROGRESS) progress = (double)(i+1) / Z() * 100;		//update the progress counter
+		}
+		free(temp);
+		return true;
+	}
+
+	/// Calculate the median value for all masked (or valid) pixels in a band and returns the median spectrum
+
+	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
+	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
+	bool median_spectrum(double* m, unsigned char* mask = NULL, bool PROGRESS = false){
+		size_t XY = X() * Y();
+		size_t count = nnz(mask);						//count will store the number of masked pixels
+		T* temp = (T*)malloc(sizeof(T) * XY);
+		
+		std::vector<T> band_values(count);				//create an STD vector of band values
+
+		//this loops goes through each band in B (Z())
+		//	masked (or valid) pixels from that band are averaged and the average is stored in p
+		size_t k;
+		for (size_t i = 0; i < Z(); i++){							//for each band
+			band_index(temp, i);									//get the band image and store it in temp
+			k = 0;													//initialize the band_value index to zero
+			for (size_t j = 0; j < XY; j++){						//loop through temp, averaging valid pixels
+				if (mask == NULL || mask[j] != 0){
+					band_values[k] = temp[j];				//store the value in the band_values array
+					k++;											//increment the band_values index
 				}
 			}
-			if(PROGRESS) progress = (double)(i+1) / Z() * 100;
+			std::sort(band_values.begin(), band_values.end());		//sort all of the values in the band
+			m[i] = band_values[ count/2 ];							//store the center value in the array
+			if(PROGRESS) progress = (double)(i+1) / Z() * 100;		//update the progress counter
 		}
 		free(temp);
 		return true;
@@ -1380,12 +1380,12 @@ public:
 	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
 	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
-	bool avg_band(double * p, unsigned char* mask, bool PROGRESS = false){
+	bool mean_spectrum(double * p, double* std, unsigned char* mask, bool PROGRESS = false){
 		if (header.interleave == envi_header::BSQ){
 			if (header.data_type == envi_header::float32)
-				return ((bsq<float>*)file)->avg_band(p, mask, PROGRESS);
+				return ((bsq<float>*)file)->mean_spectrum(p, std, mask, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->avg_band(p,  mask, PROGRESS);
+				return ((bsq<double>*)file)->mean_spectrum(p, std, mask, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1393,9 +1393,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIL){
 			if (header.data_type == envi_header::float32)
-				return ((bil<float>*)file)->avg_band(p,  mask, PROGRESS);
+				return ((bil<float>*)file)->mean_spectrum(p, std, mask, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->avg_band(p,  mask, PROGRESS);
+				return ((bil<double>*)file)->mean_spectrum(p, std, mask, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1403,9 +1403,9 @@ public:
 		}
 		else if (header.interleave == envi_header::BIP){
 			if (header.data_type == envi_header::float32)
-				return ((bip<float>*)file)->avg_band(p, mask, PROGRESS);
+				return ((bip<float>*)file)->mean_spectrum(p, std, mask, PROGRESS);
 			else if (header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->avg_band(p, mask, PROGRESS);
+				return ((bip<double>*)file)->mean_spectrum(p, std, mask, PROGRESS);
 			else{
 				std::cout << "ERROR: unidentified data type" << std::endl;
 				exit(1);
@@ -1414,6 +1414,28 @@ public:
 		return false;
 	}
+	/// Calculate the mean value for all masked (or valid) pixels in a band and returns the average spectrum
+
+	/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
+	/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
+	bool median_spectrum(double* m, unsigned char* mask, bool PROGRESS = false){
+		if (header.interleave == envi_header::BSQ){
+			if (header.data_type == envi_header::float32)
+				return ((bsq<float>*)file)->median_spectrum(m, mask, PROGRESS);
+			else if (header.data_type == envi_header::float64)
+				return ((bsq<double>*)file)->median_spectrum(m, mask, PROGRESS);
+			else{
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		else{
+			std::cout<<"ERROR: median calculation is only supported for BSQ interleave types. Convert to process."<<std::endl;
+			exit(1);
+		}
+		return false;
+	}
+
 	/// Calculate the covariance matrix for all masked pixels in the image.
 	/// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
@@ -473,7 +473,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 			glEndList();					///finilize the display list.
 			#ifdef DEBUG
 				for(int i = 0; i < numSamplesPos; i++)
-					std::cout << pV[i] << std::endl;
+					std::cout << pV[i].str() << std::endl;
 			#endif
 		}
@@ -1145,8 +1145,8 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 				out[3] = temp[2];
 			}
 			#ifdef DEBUG
-				std::cout << "out is " << out << std::endl;
-				std::cout << "when rotating from " << from << " to " << dir << std::endl;
+				std::cout << "out is " << out.str() << std::endl;
+				std::cout << "when rotating from " << from.str() << " to " << dir.str() << std::endl;
 			#endif
 			return out;
 		}
@@ -1531,7 +1531,7 @@ class gl_spider // : public virtual gl_texture&lt;T&gt;
 				setMagnitude(curSeedMag);
 				#ifdef DEBUG
-					std::cout << "The new seed " << curSeed << curSeedVec << curSeedMag << std::endl;
+					std::cout << "The new seed " << curSeed.str() << curSeedVec.str() << curSeedMag << std::endl;
 				#endif
 //			Bind(direction_texID, direction_buffID, numSamples, n_pixels);
@@ -74,11 +74,11 @@ struct vec : public std::vector&lt;T&gt;
 			at(i) = other[i];
 		}
 	}
-
+
 //	vec( vec3<T>& other){
 //		resize(3);							//resize the current vector to match the copy
 //		for(size_t i=0; i<3; i++){	//copy each element
-//			at(i) = other[i];
+//			at(i) = other[i];
 //		}
 //	}
@@ -139,16 +139,16 @@ struct vec : public std::vector&lt;T&gt;
 	}
-	
-	vec<T> cyl2cart() const
-	{
-		vec<T> cyl;
-		cyl.push_back(at(0)*std::sin(at(1)));
-		cyl.push_back(at(0)*std::cos(at(1)));
-		cyl.push_back(at(2));
-		return(cyl);
-		
-	}
+	
+	vec<T> cyl2cart() const
+	{
+		vec<T> cyl;
+		cyl.push_back(at(0)*std::sin(at(1)));
+		cyl.push_back(at(0)*std::cos(at(1)));
+		cyl.push_back(at(2));
+		return(cyl);
+		
+	}
 	/// Convert the vector from cartesian to spherical coordinates (x, y, z -> r, theta, phi where theta = [0, 2*pi])
 	vec<T> cart2sph() const
 	{
@@ -335,16 +335,16 @@ struct vec : public std::vector&lt;T&gt;
 		return *this;
 	}
-	/// Cast to a vec3
-	operator stim::vec3<T>(){
-		stim::vec3<T> r;
-		size_t N = std::min<size_t>(size(), 3);
-		for(size_t i = 0; i < N; i++)
-			r[i] = at(i);
-		return r;
-	}
-
-
+	/// Cast to a vec3
+	operator stim::vec3<T>(){
+		stim::vec3<T> r;
+		size_t N = min(size(), (size_t)3);
+		for(size_t i = 0; i < N; i++)
+			r[i] = at(i);
+		return r;
+	}
+
+
 	/// Casting and assignment
 	template<typename Y>
 	vec<T> & operator=(vec<Y> rhs){
@@ -355,16 +355,16 @@ struct vec : public std::vector&lt;T&gt;
 			at(i) = rhs[i];
 		return *this;
 	}
-
-	/// Assign a vec = vec3
-	template<typename Y>
-	vec<T> & operator=(vec3<Y> rhs)
-	{
-		resize(3);
-		for(size_t i=0; i<3; i++)
-			at(i) = rhs[i];
-		return *this;
-	}
+
+	/// Assign a vec = vec3
+	template<typename Y>
+	vec<T> & operator=(vec3<Y> rhs)
+	{
+		resize(3);
+		for(size_t i=0; i<3; i++)
+			at(i) = rhs[i];
+		return *this;
+	}
 	/// Unary minus (returns the negative of the vector)
 	vec<T> operator-() const{