fixed the crop function to work with wavelength

David Mayerich
1 parent 95f1e985
Showing 3 changed files with 158 additions and 53 deletions Show diff stats
stim/envi/envi.h
stim/envi/hsi.h
stim/math/matrix.h
@@ -1368,6 +1368,39 @@ public:
 		return false;
 	}
+	void band_bounds(double wavelength, size_t& low, size_t& high) {
+		if (header.interleave == envi_header::BSQ) {		//if the infile is bsq file
+			if (header.data_type == envi_header::float32)
+				((bsq<float>*)file)->band_bounds(wavelength, low, high);
+			else if (header.data_type == envi_header::float64)
+				((bsq<double>*)file)->band_bounds(wavelength, low, high);
+			else {
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		else if (header.interleave == envi_header::BIL) {
+			if (header.data_type == envi_header::float32)
+				((bil<float>*)file)->band_bounds(wavelength, low, high);
+			else if (header.data_type == envi_header::float64)
+				((bil<double>*)file)->band_bounds(wavelength, low, high);
+			else {
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+		else if (header.interleave == envi_header::BIP) {
+			if (header.data_type == envi_header::float32)
+				((bip<float>*)file)->band_bounds(wavelength, low, high);
+			else if (header.data_type == envi_header::float64)
+				((bip<double>*)file)->band_bounds(wavelength, low, high);
+			else {
+				std::cout << "ERROR: unidentified data type" << std::endl;
+				exit(1);
+			}
+		}
+	}
+
 	// Retrieve a spectrum at the specified 1D location
 	/// @param ptr is a pointer to pre-allocated memory of size B*sizeof(T)
@@ -62,31 +62,6 @@ protected:
 		return (T)((1.0 - alpha) * low_v + alpha * high_v);							//interpolate
 	}
-	/// Gets the two band indices surrounding a given wavelength
-	void band_bounds(double wavelength, size_t& low, size_t& high){
-		size_t B = Z();
-		for(high = 0; high < B; high++){
-			if(w[high] > wavelength) break;
-		}
-		low = 0;
-		if(high > 0)
-			low = high-1;
-	}
-
-	/// Get the list of band numbers that bound a list of wavelengths
-	void band_bounds(std::vector<double> wavelengths,
-					 std::vector<unsigned long long>& low_bands,
-					 std::vector<unsigned long long>& high_bands){
-
-		unsigned long long W = w.size();									//get the number of wavelengths in the list
-		low_bands.resize(W);												//pre-allocate space for the band lists
-		high_bands.resize(W);
-
-		for(unsigned long long wl = 0; wl < W; wl++){						//for each wavelength
-			band_bounds(wavelengths[wl], low_bands[wl], high_bands[wl]);	//find the low and high bands
-		}
-	}
-
 	/// Returns the interpolated in the given spectrum based on the given wavelength
 	/// @param s is the spectrum in main memory of length Z()
@@ -139,6 +114,31 @@ protected:
 	}
 public:
+
+	/// Gets the two band indices surrounding a given wavelength
+	void band_bounds(double wavelength, size_t& low, size_t& high) {
+		size_t B = Z();
+		for (high = 0; high < B; high++) {
+			if (w[high] > wavelength) break;
+		}
+		low = 0;
+		if (high > 0)
+			low = high - 1;
+	}
+
+	/// Get the list of band numbers that bound a list of wavelengths
+	void band_bounds(std::vector<double> wavelengths,
+		std::vector<unsigned long long>& low_bands,
+		std::vector<unsigned long long>& high_bands) {
+
+		unsigned long long W = w.size();									//get the number of wavelengths in the list
+		low_bands.resize(W);												//pre-allocate space for the band lists
+		high_bands.resize(W);
+
+		for (unsigned long long wl = 0; wl < W; wl++) {						//for each wavelength
+			band_bounds(wavelengths[wl], low_bands[wl], high_bands[wl]);	//find the low and high bands
+		}
+	}
 			/// Get a mask that has all pixels with inf or NaN values masked out (false)
 	void mask_finite(unsigned char* out_mask, unsigned char* mask, bool PROGRESS = false){
 		size_t XY = X() * Y();
@@ -33,32 +33,58 @@ namespace stim{
 		}
 	}
+	//class encapsulates a mat4 file, and can be used to write multiple matrices to a single mat4 file
+	class mat4file {
+		std::ofstream matfile;
+
+	public:
+		/// Constructor opens a mat4 file for writing
+		mat4file(std::string filename) {
+			matfile.open(filename, std::ios::binary);
+		}
+
+		bool is_open() {
+			return matfile.is_open();
+		}
+
+		void close() {
+			matfile.close();
+		}
+
+		bool writemat(char* data, std::string varname, size_t sx, size_t sy, mat4Format format) {
+			//save the matrix file here (use the mat4 function above)
+			//data format: https://maxwell.ict.griffith.edu.au/spl/matlab-page/matfile_format.pdf (page 32)
+
+			int MOPT = 0;									//initialize the MOPT type value to zero
+			int m = 0;										//little endian
+			int o = 0;										//reserved, always 0
+			int p = format;
+			int t = 0;
+			MOPT = m * 1000 + o * 100 + p * 10 + t;			//calculate the type value
+			int mrows = (int)sx;
+			int ncols = (int)sy;
+			int imagf = 0;									//assume real (for now)
+			varname.push_back('\0');									//add a null to the string
+			int namlen = (int)varname.size();						//calculate the name size
+
+			size_t bytes = sx * sy * mat4Format_size(format);
+			matfile.write((char*)&MOPT, 4);
+			matfile.write((char*)&mrows, 4);
+			matfile.write((char*)&ncols, 4);
+			matfile.write((char*)&imagf, 4);
+			matfile.write((char*)&namlen, 4);
+			matfile.write((char*)&varname[0], namlen);
+			matfile.write((char*)data, bytes);				//write the matrix data
+			return is_open();
+		}
+	};
+
 	static void save_mat4(char* data, std::string filename, std::string varname, size_t sx, size_t sy, mat4Format format){
-		//save the matrix file here (use the mat4 function above)
-		//data format: https://maxwell.ict.griffith.edu.au/spl/matlab-page/matfile_format.pdf (page 32)
-		
-		int MOPT = 0;									//initialize the MOPT type value to zero
-		int m = 0;										//little endian
-		int o = 0;										//reserved, always 0
-		int p = format;
-		int t = 0;
-		MOPT = m * 1000 + o * 100 + p * 10 + t;			//calculate the type value
-		int mrows = (int)sx;
-		int ncols = (int)sy;
-		int imagf = 0;									//assume real (for now)
-		varname.push_back('\0');									//add a null to the string
-		int namlen = (int)varname.size();						//calculate the name size
-
-		size_t bytes = sx * sy * mat4Format_size(format);
-		std::ofstream outfile(filename, std::ios::binary);
-		outfile.write((char*)&MOPT, 4);
-		outfile.write((char*)&mrows, 4);
-		outfile.write((char*)&ncols, 4);
-		outfile.write((char*)&imagf, 4);
-		outfile.write((char*)&namlen, 4);
-		outfile.write((char*)&varname[0], namlen);
-		outfile.write((char*)data, bytes);				//write the matrix data
-		outfile.close();
+		mat4file outfile(filename);									//create a mat4 file object
+		if (outfile.is_open()) {									//if the file is open
+			outfile.writemat(data, varname, sx, sy, format);		//write the matrix
+			outfile.close();										//close the file
+		}		
 	}
 template <class T>
@@ -409,8 +435,21 @@ public:
 		}
 	}
-	// saves the matrix as a Level-4 MATLAB file
-	void mat4(std::string filename, std::string name = std::string("unknown"), mat4Format format = mat4_float) {
+	void mat4(stim::mat4file& file, std::string name = std::string("unknown"), mat4Format format = mat4_float) {
+		//make sure the matrix name is valid (only numbers and letters, with a letter at the beginning
+		for (size_t c = 0; c < name.size(); c++) {
+			if (name[c] < 48 ||												//if the character isn't a number or letter, replace it with '_'
+				(name[c] > 57 && name[c] < 65) ||
+				(name[c] > 90 && name[c] < 97) ||
+				(name[c] > 122)) {
+				name[c] = '_';
+			}
+		}
+		if (name[0] < 65 ||
+			(name[0] > 91 && name[0] < 97) ||
+			name[0] > 122) {
+			name = std::string("m") + name;
+		}
 		if (format == mat4_float) {
 			if (sizeof(T) == 4) format = mat4_float32;
 			else if (sizeof(T) == 8) format = mat4_float64;
@@ -419,7 +458,40 @@ public:
 				exit(1);
 			}
 		}
-		stim::save_mat4((char*)M, filename, name, rows(), cols(), format);
+		//the name is now valid
+
+		//if the size of the array is more than 100,000,000 elements, the matrix isn't supported
+		if (rows() * cols() > 100000000) {											//break the matrix up into multiple parts
+			//mat4file out(filename);													//create a mat4 object to write the matrix
+			if (file.is_open()) {
+				if (rows() < 100000000) {												//if the size of the row is less than 100,000,000, split the matrix up by columns
+					size_t ncols = 100000000 / rows();									//calculate the number of columns that can fit in one matrix
+					size_t nmat = (size_t)std::ceil((double)cols() / (double)ncols);			//calculate the number of matrices required
+					for (size_t m = 0; m < nmat; m++) {									//for each matrix
+						std::stringstream ss;
+						ss << name << "_part_" << m + 1;
+						if (m == nmat - 1)
+							file.writemat((char*)(data() + m * ncols * rows()), ss.str(), rows(), cols() - m * ncols, format);
+						else
+							file.writemat((char*)(data() + m * ncols * rows()), ss.str(), rows(), ncols, format);
+					}
+				}
+			}
+		}
+		//call the mat4 subroutine
+		else
+			//stim::save_mat4((char*)M, filename, name, rows(), cols(), format);
+			file.writemat((char*)data(), name, rows(), cols(), format);
+	}
+
+	// saves the matrix as a Level-4 MATLAB file
+	void mat4(std::string filename, std::string name = std::string("unknown"), mat4Format format = mat4_float) {
+		stim::mat4file matfile(filename);
+
+		if (matfile.is_open()) {
+			mat4(matfile, name, format);
+			matfile.close();
+		}
 	}
 };