Blame view

stim/envi/hsi.h 8.55 KB
9d3ba0b1   David Mayerich   added stim::hsi a...
1
2
3
4
5
6
7
8
  #ifndef STIM_HSI_H
  #define STIM_HSI_H
  
  #include "../envi/envi_header.h"
  #include "../envi/binary.h"
  #include <cstring>
  #include <utility>
  
ba51ae6a   David Mayerich   fixed metric calc...
9
10
11
12
13
14
  #ifdef _WIN32
  	#include <float.h>
  #else
  	#include<cmath>
  #endif
  
9d3ba0b1   David Mayerich   added stim::hsi a...
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
  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
  
1a224b6a   David Mayerich   fixed linux compa...
35
36
37
  	using binary<T>::R;
  	using binary<T>::progress;
  
9d3ba0b1   David Mayerich   added stim::hsi a...
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
  	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
  	}
  
  	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
  	}
  
  	/// Gets the two band indices surrounding a given wavelength
  	void band_bounds(double wavelength, unsigned long long& low, unsigned long long& high){
  		unsigned long long 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
1a224b6a   David Mayerich   fixed linux compa...
76
77
  	void band_bounds(std::vector<double> wavelengths,
  					 std::vector<unsigned long long>& low_bands,
9d3ba0b1   David Mayerich   added stim::hsi a...
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
  					 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()
  	/// @param wavelength is the wavelength value to interpolate out
  	T interp_spectrum(T* s, double wavelength){
  		unsigned long long 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
1a224b6a   David Mayerich   fixed linux compa...
98
  
9d3ba0b1   David Mayerich   added stim::hsi a...
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
  		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;
  	}
a2bf1d08   David Mayerich   general bug fixes...
114
115
116
117
118
119
120
121
122
123
124
  
  	/// 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)
  	}
ba51ae6a   David Mayerich   fixed metric calc...
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
  
  	/// 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:
  			/// Get a mask that has all pixels with inf or NaN values masked out (false)
e37ce7bf   David Mayerich   added the ability...
142
  	void mask_finite(unsigned char* out_mask, unsigned char* mask, bool PROGRESS = false){
ba51ae6a   David Mayerich   fixed metric calc...
143
  		size_t XY = X() * Y();
e37ce7bf   David Mayerich   added the ability...
144
  		if(mask == NULL)												//if no mask is provided
41e2f5ab   David Mayerich   fixed finite mask...
145
  			memset(out_mask, 255, XY * sizeof(unsigned char));				//initialize the mask to 255
e37ce7bf   David Mayerich   added the ability...
146
147
  		else															//if a mask is provided
  			memcpy(out_mask, mask, XY * sizeof(unsigned char));			//initialize the current mask to that one
ba51ae6a   David Mayerich   fixed metric calc...
148
  		T* page = (T*)malloc(R[0] * R[1] * sizeof(T));		//allocate space for a page of data
1a224b6a   David Mayerich   fixed linux compa...
149
  
ba51ae6a   David Mayerich   fixed metric calc...
150
151
152
  		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
1a224b6a   David Mayerich   fixed linux compa...
153
  
ba51ae6a   David Mayerich   fixed metric calc...
154
  #ifdef _WIN32
3c714475   David Mayerich   fixed incorrect c...
155
  				if(!_finite(page[i])){						//if the value at index i is not finite
ba51ae6a   David Mayerich   fixed metric calc...
156
  #else
3c714475   David Mayerich   fixed incorrect c...
157
  				if(!std::isfinite(page[i])){				//C++11 implementation
ba51ae6a   David Mayerich   fixed metric calc...
158
159
  #endif
  					size_t x, y, b;
3c714475   David Mayerich   fixed incorrect c...
160
  					xyb(p * R[0] * R[1] + i, x, y, b);		//find the 3D coordinates of the value
e37ce7bf   David Mayerich   added the ability...
161
  					out_mask[ y * X() + x ] = 0;				//remove the pixel (it's bad)
ba51ae6a   David Mayerich   fixed metric calc...
162
163
164
165
166
167
  				}
  			}
  			if(PROGRESS) progress = (double)(p + 1) / (double)R[2] * 100;
  		}
  	}
  
2ce6954b   David Mayerich   added the ability...
168
169
170
171
172
173
174
175
176
177
178
179
  	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
1a224b6a   David Mayerich   fixed linux compa...
180
  	void calc_combined_size(long long xp, long long yp, long long Sx, long long Sy,
2ce6954b   David Mayerich   added the ability...
181
182
183
184
185
186
187
  							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
1a224b6a   David Mayerich   fixed linux compa...
188
  
2ce6954b   David Mayerich   added the ability...
189
190
191
192
193
194
195
196
197
198
199
  		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
1a224b6a   David Mayerich   fixed linux compa...
200
201
  		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
2ce6954b   David Mayerich   added the ability...
202
203
204
205
206
207
  
  		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
  		}
  	}
  
6db250f5   David Mayerich   adaptive streamin...
208
  	size_t read(T* dest, size_t x, size_t y, size_t z, size_t sx, size_t sy, size_t sz){
8b899c24   David Mayerich   fixed asynchronou...
209
210
211
212
213
214
215
216
217
218
219
  		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;
  
6db250f5   David Mayerich   adaptive streamin...
220
  		return binary<T>::read(dest, d[0], d[1], d[2], sd[0], sd[1], sd[2]);
8b899c24   David Mayerich   fixed asynchronou...
221
222
  	}
  
9d3ba0b1   David Mayerich   added stim::hsi a...
223
224
225
226
227
  };
  
  }		//end namespace STIM
  
  #endif