codebase / stimlib

  #ifndef STIM_HSI_H
  #define STIM_HSI_H
  

  #include <stim/envi/envi_header.h>
  #include <stim/envi/binary.h>

  #include <cstring>
  #include <utility>
  

  #ifdef _WIN32
  	#include <float.h>
  #else
  	#include<cmath>
  #endif
  

  namespace stim{
  
  /**
  	The BIL class represents a 3-dimensional binary file stored using band interleaved by line (BIL) image encoding. The binary file is stored
  	such that X-Z "frames" are stored sequentially to form an image stack along the y-axis. When accessing the data sequentially on disk,
  	the dimensions read, from fastest to slowest, are X, Z, Y.
  
  	This class is optimized for data streaming, and therefore supports extremely large (terabyte-scale) files. Data is loaded from disk
  	on request. Functions used to access data are written to support efficient reading.
  */
  
  template <typename T>
  class hsi: public binary<T> {
  
  protected:
  	unsigned char O[3];					//order of dimensions (orientation on disk)
  										//[X Y B]: [0 1 2] = bsq, [0 2 1] = bil, [1 2 0] = bip
  
  	std::vector<double> w;		//band wavelengths
  

  	using binary<T>::R;
  	using binary<T>::progress;
  

  	unsigned long long X(){ return R[O[0]]; }
  	unsigned long long Y(){ return R[O[1]]; }
  	unsigned long long Z(){ return R[O[2]]; }
  
  	/// Initialize axis orders based on common interleave formats
  	void init_bsq(){ O[0] = 0; O[1] = 1; O[2] = 2; }
  	void init_bil(){ O[0] = 0; O[1] = 2; O[2] = 1; }
  	void init_bip(){ O[0] = 1; O[1] = 2; O[2] = 0; }
  
  	/// Calculate the number of masked pixels in a given mask
  	unsigned long long nnz(unsigned char* mask){
  		unsigned long long XY = X() * Y();							//calculate the total number of pixels in the HSI
  		if(mask == NULL) return XY;									//if the mask is null, assume that all of the pixels are masked
  
  		unsigned long long n = 0;									//initialize the number of masked pixels to zero (0)
  		for(unsigned long long xy = 0; xy < XY; xy++)				//for each pixel in the HSI
  			if(mask[xy]) n++;										//if the mask value is nonzero, increment the number of masked pixels
  		return n;													//return the number of masked pixels
  	}
  

  	//perform linear interpolation between two bands

  	T lerp(double w, T low_v, double low_w, T high_v, double high_w){
  		if(low_w == high_w) return low_v;										//if the interval is of zero length, just return one of the bounds
  		double alpha = (w - low_w) / (high_w - low_w);							//calculate the interpolation factor
  		return (T)((1.0 - alpha) * low_v + alpha * high_v);							//interpolate
  	}
  

  	/// Returns the interpolated in the given spectrum based on the given wavelength
  
  	/// @param s is the spectrum in main memory of length Z()
  	/// @param wavelength is the wavelength value to interpolate out
  	T interp_spectrum(T* s, double wavelength){

  		size_t low, high;								//indices for the bands surrounding wavelength

  		band_bounds(wavelength, low, high);							//get the surrounding band indices
  
  		if(high == w.size()) return s[w.size()-1];					//if the high band is above the wavelength range, return the highest wavelength value

  		return lerp(wavelength, s[low], w[low], s[high], w[high]);
  	}
  
  	/// Returns the interpolated value corresponding to the given list of wavelengths
  	std::vector<T> interp_spectrum(T* s, std::vector<double> wavelengths){
  
  		unsigned long long N = wavelengths.size();						//get the number of wavelength measurements
  
  		std::vector<T> v;												//allocate space for the resulting values
  		v.resize(wavelengths.size());
  		for(unsigned long long n = 0; n < N; n++){						//for each measurement
  			v[n] = interp_spectrum(s, wavelengths[n]);					//interpolate the measurement
  		}
  		return v;
  	}

  
  	/// Returns the 1D on-disk index of a specified pixel location and band
  	size_t idx(size_t x, size_t y, size_t b){
  		size_t c[3];										//generate a coefficient list
  
  		c[O[0]] = x;										//assign the coordinates based on the coefficient order
  		c[O[1]] = y;
  		c[O[2]] = b;
  
  		return c[2] * R[0] * R[1] + c[1] * R[0] + c[0];		//calculate and return the index (trust me this works)
  	}

  
  	/// Returns the 3D coordinates of a specified index
  	void xyb(size_t idx, size_t& x, size_t& y, size_t& b){
  
  		size_t c[3];
  
  		c[2] = idx / (R[0] * R[1]);
  		c[1] = (idx - c[2] * R[0] * R[1]) / R[0];
  		c[0] = idx - c[2] * R[0] * R[1] - c[1] * R[0];
  
  		x = c[O[0]];
  		y = c[O[1]];
  		b = c[O[2]];
  	}
  
  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<size_t>& low_bands,
  		std::vector<size_t>& 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 (size_t 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();

  		if(mask == NULL)												//if no mask is provided

  			memset(out_mask, 255, XY * sizeof(unsigned char));				//initialize the mask to 255

  		else															//if a mask is provided
  			memcpy(out_mask, mask, XY * sizeof(unsigned char));			//initialize the current mask to that one

  		T* page = (T*)malloc(R[0] * R[1] * sizeof(T));		//allocate space for a page of data

  		for(size_t p = 0; p < R[2]; p++){					//for each page
  			binary<T>::read_page(page, p);					//read a page
  			for(size_t i = 0; i < R[0] * R[1]; i++){		//for each pixel in that page

  #ifdef _WIN32

  				if(!_finite(page[i])){						//if the value at index i is not finite

  #else

  				if(!std::isfinite(page[i])){				//C++11 implementation

  #endif
  					size_t x, y, b;

  					xyb(p * R[0] * R[1] + i, x, y, b);		//find the 3D coordinates of the value

  					out_mask[ y * X() + x ] = 0;				//remove the pixel (it's bad)

  				}
  			}
  			if(PROGRESS) progress = (double)(p + 1) / (double)R[2] * 100;
  		}
  	}
  

  	void calc_combined_size(long long xp, long long yp, long long Sx, long long Sy,
  							size_t& Sfx, size_t& Sfy){
  		long long left = std::min<long long>(0, xp);										//calculate the left edge of the final image
  		long long right = std::max<long long>((long long)X(), Sx + xp);						//calculate the right edge of the final image
  		long long top = std::min<long long>(0, yp);											//calculate the top edge of the final image
  		long long bottom = std::max<long long>((long long)Y(), Sy + yp);					//calculate the bottom edge of the final image
  
  		Sfx = right - left;
  		Sfy = bottom - top;													//calculate the size of the final image
  	}
  
  	/// Calculates the necessary image size required to combine two images given the specified offset (position) of the second image

  	void calc_combined_size(long long xp, long long yp, long long Sx, long long Sy,

  							size_t& Sfx, size_t& Sfy,
  							size_t& p0_x, size_t& p0_y,
  							size_t& p1_x, size_t& p1_y){
  
  		calc_combined_size(xp, yp, Sx, Sy, Sfx, Sfy);
  
  		p0_x = p0_y = p1_x = p1_y = 0;										//initialize all boundary positions to zero

  		if(xp < 0) p0_x = -xp;												//set the left positions of the current and source image
  		else p1_x = xp;
  		if(yp < 0) p0_y = -yp;
  		else p1_y = yp;
  	}
  
  	/// Inserts an image of a band into a larger image
  	void pad_band(T* padded, T* src, size_t x0, size_t x1, size_t y0, size_t y1){
  
  		size_t w = X();										//calculate the number of pixels in a line
  		size_t wb = w * sizeof(T);							//calculate the number of bytes in a line

  		size_t pw = (X() + x0 + x1);								//calculate the width of the padded image
  		size_t pwb = pw * sizeof(T);								//calculate the number of bytes in a line of the padded image

  
  		for(size_t y = 0; y < Y(); y++){								//for each line in the real image
  			memcpy( &padded[ (y + y0) * pw + x0 ], &src[y * w], wb );	//use memcpy to copy the line to the appropriate place in the padded image
  		}
  	}
  

  	size_t read(T* dest, size_t x, size_t y, size_t z, size_t sx, size_t sy, size_t sz){

  		size_t d[3];					//position in the binary coordinate system
  		size_t sd[3];					//size in the binary coordinate system
  
  		d[O[0]] = x;					//set the point in the binary coordinate system
  		d[O[1]] = y;
  		d[O[2]] = z;
  
  		sd[O[0]] = sx;					//set the size in the binary coordinate system
  		sd[O[1]] = sy;
  		sd[O[2]] = sz;
  

  		return binary<T>::read(dest, d[0], d[1], d[2], sd[0], sd[1], sd[2]);

  	}
  

  };
  
  }		//end namespace STIM
  

  #endif