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
142
143
144
145
|
/// 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)
void mask_finite(unsigned char* mask, bool PROGRESS = false){
size_t XY = X() * Y();
memset(mask, 255, XY * sizeof(unsigned char)); //initialize the mask to zero (0)
T* page = (T*)malloc(R[0] * R[1] * sizeof(T)); //allocate space for a page of data
|
1a224b6a
David Mayerich
fixed linux compa...
|
146
|
|
ba51ae6a
David Mayerich
fixed metric calc...
|
147
148
149
|
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...
|
150
|
|
ba51ae6a
David Mayerich
fixed metric calc...
|
151
|
#ifdef _WIN32
|
3c714475
David Mayerich
fixed incorrect c...
|
152
|
if(!_finite(page[i])){ //if the value at index i is not finite
|
ba51ae6a
David Mayerich
fixed metric calc...
|
153
|
#else
|
3c714475
David Mayerich
fixed incorrect c...
|
154
|
if(!std::isfinite(page[i])){ //C++11 implementation
|
ba51ae6a
David Mayerich
fixed metric calc...
|
155
156
|
#endif
size_t x, y, b;
|
3c714475
David Mayerich
fixed incorrect c...
|
157
158
|
xyb(p * R[0] * R[1] + i, x, y, b); //find the 3D coordinates of the value
mask[ y * X() + x ] = 0; //remove the pixel (it's bad)
|
ba51ae6a
David Mayerich
fixed metric calc...
|
159
160
161
162
163
164
|
}
}
if(PROGRESS) progress = (double)(p + 1) / (double)R[2] * 100;
}
}
|
2ce6954b
David Mayerich
added the ability...
|
165
166
167
168
169
170
171
172
173
174
175
176
|
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...
|
177
|
void calc_combined_size(long long xp, long long yp, long long Sx, long long Sy,
|
2ce6954b
David Mayerich
added the ability...
|
178
179
180
181
182
183
184
|
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...
|
185
|
|
2ce6954b
David Mayerich
added the ability...
|
186
187
188
189
190
191
192
193
194
195
196
|
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...
|
197
198
|
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...
|
199
200
201
202
203
204
|
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
}
}
|
846b3f18
David Mayerich
fixed asynchronou...
|
205
|
void 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...
|
206
207
208
209
210
211
212
213
214
215
216
|
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;
|
846b3f18
David Mayerich
fixed asynchronou...
|
217
|
if(!binary<T>::read(dest, d[0], d[1], d[2], sd[0], sd[1], sd[2])){
|
8b899c24
David Mayerich
fixed asynchronou...
|
218
219
220
221
222
|
std::cout<<"error reading block in stim::hsi: ("<<d[0]<<", "<<d[1]<<", "<<d[2]<<") - ["<<sd[0]<<", "<<sd[1]<<", "<<sd[2]<<"]"<<std::endl;
exit(1);
}
}
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
223
224
225
226
227
|
};
} //end namespace STIM
#endif
|