hsi.h
4.64 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
#ifndef STIM_HSI_H
#define STIM_HSI_H
#include "../envi/envi_header.h"
#include "../envi/binary.h"
#include <cstring>
#include <utility>
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
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
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()
/// @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
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)
}
};
} //end namespace STIM
#endif