6aa04ba2
David Mayerich
interleave types
|
1
2
3
|
#ifndef STIM_BIL_H
#define STIM_BIL_H
|
e8eb202f
David Mayerich
added a new ENVI ...
|
4
|
#include "../envi/envi_header.h"
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
5
|
#include "../envi/hsi.h"
|
1a224b6a
David Mayerich
fixed linux compa...
|
6
|
#include "../math/fd_coefficients.h"
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
7
|
#include <stim/cuda/cudatools/error.h>
|
c25e7d0d
heziqi
speed of bip.base...
|
8
9
|
#include <cstring>
#include <utility>
|
1a224b6a
David Mayerich
fixed linux compa...
|
10
|
#include <deque>
|
c25e7d0d
heziqi
speed of bip.base...
|
11
|
|
8a86bd56
David Mayerich
changed rts names...
|
12
|
namespace stim{
|
c25e7d0d
heziqi
speed of bip.base...
|
13
|
|
a23c4132
David Mayerich
Doxygen comments ...
|
14
15
16
17
18
19
20
21
22
|
/**
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.
*/
|
c25e7d0d
heziqi
speed of bip.base...
|
23
24
|
template <typename T>
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
25
|
class bil: public hsi<T> {
|
c25e7d0d
heziqi
speed of bip.base...
|
26
27
|
protected:
|
c25e7d0d
heziqi
speed of bip.base...
|
28
|
|
1a224b6a
David Mayerich
fixed linux compa...
|
29
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
30
31
|
using hsi<T>::w; //use the wavelength array in stim::hsi
using hsi<T>::nnz;
|
63fc1df8
David Mayerich
added an inverse ...
|
32
|
using binary<T>::progress;
|
1a224b6a
David Mayerich
fixed linux compa...
|
33
34
35
|
using hsi<T>::X;
using hsi<T>::Y;
using hsi<T>::Z;
|
cbce396e
David Mayerich
added support for...
|
36
|
|
c25e7d0d
heziqi
speed of bip.base...
|
37
38
39
40
|
public:
using binary<T>::open;
using binary<T>::file;
|
6708cc25
heziqi
Ziqi added envi c...
|
41
|
using binary<T>::R;
|
c25e7d0d
heziqi
speed of bip.base...
|
42
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
43
44
|
bil(){ hsi<T>::init_bil(); }
|
a23c4132
David Mayerich
Doxygen comments ...
|
45
46
47
48
49
50
51
52
|
/// Open a data file for reading using the class interface.
/// @param filename is the name of the binary file on disk
/// @param X is the number of samples along dimension 1
/// @param Y is the number of samples (lines) along dimension 2
/// @param B is the number of samples (bands) along dimension 3
/// @param header_offset is the number of bytes (if any) in the binary header
/// @param wavelengths is an optional STL vector of size B specifying a numerical label for each band
|
1a224b6a
David Mayerich
fixed linux compa...
|
53
54
55
56
57
|
bool open(std::string filename,
unsigned long long X,
unsigned long long Y,
unsigned long long B,
unsigned long long header_offset,
|
bb457896
David Mayerich
improved error me...
|
58
59
|
std::vector<double> wavelengths,
stim::iotype io = stim::io_in){
|
c25e7d0d
heziqi
speed of bip.base...
|
60
|
|
6708cc25
heziqi
Ziqi added envi c...
|
61
|
w = wavelengths;
|
c25e7d0d
heziqi
speed of bip.base...
|
62
|
|
bb457896
David Mayerich
improved error me...
|
63
|
return open(filename, vec<unsigned long long>(X, B, Y), header_offset, io);
|
6708cc25
heziqi
Ziqi added envi c...
|
64
|
|
c25e7d0d
heziqi
speed of bip.base...
|
65
66
|
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
67
68
69
70
|
/// Retrieve a single band (based on index) and stores it in pre-allocated memory.
/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
/// @param page <= B is the integer number of the band to be copied.
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
71
|
bool band_index( T * p, unsigned long long page, bool PROGRESS = false){
|
724ec347
David Mayerich
simplified the EN...
|
72
|
//return binary<T>::read_plane_1(p, page);
|
c25e7d0d
heziqi
speed of bip.base...
|
73
|
|
63fc1df8
David Mayerich
added an inverse ...
|
74
|
if(PROGRESS) progress = 0;
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
75
76
|
unsigned long long L = X() * sizeof(T); //caculate the number of bytes in a sample line
unsigned long long jump = X() * (Z() - 1) * sizeof(T);
|
c25e7d0d
heziqi
speed of bip.base...
|
77
|
|
724ec347
David Mayerich
simplified the EN...
|
78
|
if (page >= Z()){ //make sure the bank number is right
|
c25e7d0d
heziqi
speed of bip.base...
|
79
80
81
82
|
std::cout<<"ERROR: page out of range"<<std::endl;
return false;
}
|
724ec347
David Mayerich
simplified the EN...
|
83
|
file.seekg(X() * page * sizeof(T), std::ios::beg);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
84
|
for (unsigned long long i = 0; i < Y(); i++)
|
c25e7d0d
heziqi
speed of bip.base...
|
85
|
{
|
724ec347
David Mayerich
simplified the EN...
|
86
|
file.read((char *)(p + i * X()), L);
|
c25e7d0d
heziqi
speed of bip.base...
|
87
|
file.seekg( jump, std::ios::cur);
|
63fc1df8
David Mayerich
added an inverse ...
|
88
|
if(PROGRESS) progress = (double)(i+1) / Y() * 100;
|
c25e7d0d
heziqi
speed of bip.base...
|
89
90
91
92
93
|
}
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
94
95
96
97
|
/// Retrieve a single band (by numerical label) and stores it in pre-allocated memory.
/// @param p is a pointer to an allocated region of memory at least X * Y * sizeof(T) in size.
/// @param wavelength is a floating point value (usually a wavelength in spectral data) used as a label for the band to be copied.
|
63fc1df8
David Mayerich
added an inverse ...
|
98
|
bool band( T * p, double wavelength, bool PROGRESS = false){
|
c25e7d0d
heziqi
speed of bip.base...
|
99
|
|
ee4dea28
David Mayerich
fixed errors in c...
|
100
101
|
//if there are no wavelengths in the BSQ file
if(w.size() == 0)
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
102
|
return band_index(p, (unsigned long long)wavelength);
|
ee4dea28
David Mayerich
fixed errors in c...
|
103
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
104
105
|
unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
unsigned long long S = XY * sizeof(T); //calculate the number of bytes of a band
|
c25e7d0d
heziqi
speed of bip.base...
|
106
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
107
|
unsigned long long page=0; //bands around the wavelength
|
1a224b6a
David Mayerich
fixed linux compa...
|
108
|
|
c25e7d0d
heziqi
speed of bip.base...
|
109
110
111
112
|
//get the bands numbers around the wavelength
//if wavelength is smaller than the first one in header file
|
6708cc25
heziqi
Ziqi added envi c...
|
113
|
if ( w[page] > wavelength ){
|
63fc1df8
David Mayerich
added an inverse ...
|
114
|
band_index(p, page, PROGRESS);
|
c25e7d0d
heziqi
speed of bip.base...
|
115
116
117
|
return true;
}
|
6708cc25
heziqi
Ziqi added envi c...
|
118
|
while( w[page] < wavelength )
|
c25e7d0d
heziqi
speed of bip.base...
|
119
120
121
|
{
page++;
//if wavelength is larger than the last wavelength in header file
|
724ec347
David Mayerich
simplified the EN...
|
122
|
if (page == Z()) {
|
6e9cf1ac
David Mayerich
fixed image bug i...
|
123
|
band_index(p, Z()-1, PROGRESS);
|
c25e7d0d
heziqi
speed of bip.base...
|
124
125
126
|
return true;
}
}
|
6708cc25
heziqi
Ziqi added envi c...
|
127
|
if ( wavelength < w[page] )
|
c25e7d0d
heziqi
speed of bip.base...
|
128
|
{
|
517876d6
heziqi
metrics finished ...
|
129
130
|
T * p1;
T * p2;
|
c25e7d0d
heziqi
speed of bip.base...
|
131
132
133
|
p1=(T*)malloc(S); //memory allocation
p2=(T*)malloc(S);
band_index(p1, page - 1);
|
63fc1df8
David Mayerich
added an inverse ...
|
134
|
band_index(p2, page, PROGRESS);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
135
136
137
|
for(unsigned long long i=0; i < XY; i++){
double r = (wavelength - w[page-1]) / (w[page] - w[page-1]);
p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);
|
c25e7d0d
heziqi
speed of bip.base...
|
138
|
}
|
517876d6
heziqi
metrics finished ...
|
139
140
|
free(p1);
free(p2);
|
c25e7d0d
heziqi
speed of bip.base...
|
141
142
143
|
}
else //if the wavelength is equal to a wavelength in header file
{
|
63fc1df8
David Mayerich
added an inverse ...
|
144
|
band_index(p, page, PROGRESS);
|
c25e7d0d
heziqi
speed of bip.base...
|
145
146
147
148
149
|
}
return true;
}
|
724ec347
David Mayerich
simplified the EN...
|
150
151
152
153
154
155
|
/// Retrieves a band of x values from a given xz plane.
/// @param p is a pointer to pre-allocated memory at least Z * sizeof(T) in size
/// @param c is a pointer to an existing XZ plane (size X*Z*sizeof(T))
/// @param wavelength is the wavelength to retrieve
bool read_x_from_xz(T* p, T* c, double wavelength)
|
c25e7d0d
heziqi
speed of bip.base...
|
156
|
{
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
157
158
|
unsigned long long B = Z();
unsigned long long L = X() * sizeof(T);
|
c25e7d0d
heziqi
speed of bip.base...
|
159
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
160
|
unsigned long long page=0; //samples around the wavelength
|
c25e7d0d
heziqi
speed of bip.base...
|
161
162
163
164
165
166
|
T * p1;
T * p2;
//get the bands numbers around the wavelength
//if wavelength is smaller than the first one in header file
|
6708cc25
heziqi
Ziqi added envi c...
|
167
|
if ( w[page] > wavelength ){
|
1a224b6a
David Mayerich
fixed linux compa...
|
168
|
memcpy(p, c, L);
|
c25e7d0d
heziqi
speed of bip.base...
|
169
170
171
|
return true;
}
|
6708cc25
heziqi
Ziqi added envi c...
|
172
|
while( w[page] < wavelength )
|
c25e7d0d
heziqi
speed of bip.base...
|
173
174
175
176
|
{
page++;
//if wavelength is larger than the last wavelength in header file
if (page == B) {
|
724ec347
David Mayerich
simplified the EN...
|
177
|
memcpy(p, c + (B - 1) * X(), L);
|
c25e7d0d
heziqi
speed of bip.base...
|
178
179
180
|
return true;
}
}
|
6708cc25
heziqi
Ziqi added envi c...
|
181
|
if ( wavelength < w[page] )
|
c25e7d0d
heziqi
speed of bip.base...
|
182
183
184
185
|
{
p1=(T*)malloc( L ); //memory allocation
p2=(T*)malloc( L );
|
724ec347
David Mayerich
simplified the EN...
|
186
187
|
memcpy(p1, c + (page - 1) * X(), L);
memcpy(p2, c + page * X(), L);
|
1a224b6a
David Mayerich
fixed linux compa...
|
188
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
189
|
for(unsigned long long i=0; i < X(); i++){
|
6708cc25
heziqi
Ziqi added envi c...
|
190
|
double r = (double) (wavelength - w[page-1]) / (double) (w[page] - w[page-1]);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
191
|
p[i] = (T)(((double)p2[i] - (double)p1[i]) * r + (double)p1[i]);
|
c25e7d0d
heziqi
speed of bip.base...
|
192
193
|
}
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
194
|
else //if the wavelength is equal to a wavelength in header file
|
724ec347
David Mayerich
simplified the EN...
|
195
|
memcpy(p, c + page * X(), L);
|
1a224b6a
David Mayerich
fixed linux compa...
|
196
197
|
return true;
|
c25e7d0d
heziqi
speed of bip.base...
|
198
199
|
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
200
201
202
203
204
|
/// Retrieve a single spectrum (B-axis line) at a given (x, y) location and stores it in pre-allocated memory.
/// @param p is a pointer to pre-allocated memory at least B * sizeof(T) in size.
/// @param x is the x-coordinate (dimension 1) of the spectrum.
/// @param y is the y-coordinate (dimension 2) of the spectrum.
|
735a2a24
David Mayerich
started testing o...
|
205
206
207
208
209
|
bool spectrum(T* p, size_t n, bool PROGRESS = false){
size_t y = n / X();
size_t x = n - y * X();
return binary<T>::read_line_1(p, x, y, PROGRESS);
}
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
210
|
bool spectrum(T * p, unsigned long long x, unsigned long long y, bool PROGRESS = false){
|
63fc1df8
David Mayerich
added an inverse ...
|
211
|
return binary<T>::read_line_1(p, x, y, PROGRESS);
|
c25e7d0d
heziqi
speed of bip.base...
|
212
|
}
|
bfe9c6db
heziqi
added feature_mat...
|
213
|
|
a23c4132
David Mayerich
Doxygen comments ...
|
214
215
216
217
|
/// Retrieve a single pixel and stores it in pre-allocated memory.
/// @param p is a pointer to pre-allocated memory at least sizeof(T) in size.
/// @param n is an integer index to the pixel using linear array indexing.
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
218
|
bool pixel(T * p, unsigned long long n){
|
bfe9c6db
heziqi
added feature_mat...
|
219
220
|
//calculate the corresponding x, y
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
221
222
|
unsigned long long x = n % X();
unsigned long long y = n / X();
|
bfe9c6db
heziqi
added feature_mat...
|
223
224
225
226
|
//get the pixel
return spectrum(p, x, y);
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
227
|
|
c25e7d0d
heziqi
speed of bip.base...
|
228
|
//given a Y ,return a XZ slice
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
229
|
bool read_plane_xz(T * p, size_t y){
|
724ec347
David Mayerich
simplified the EN...
|
230
|
return binary<T>::read_plane_2(p, y);
|
c25e7d0d
heziqi
speed of bip.base...
|
231
|
}
|
e843658b
Brad Deutsch
Previous push did...
|
232
|
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
233
234
235
236
237
238
239
240
241
242
243
244
245
246
|
//given a Y, return ZX slice (transposed such that the spectrum is the leading dimension)
int read_plane_zx(T* p, size_t y){
T* temp = (T*) malloc(X() * Z() * sizeof(T)); //allocate space to store the temporary xz plane
binary<T>::read_plane_2(temp, y); //load the plane from disk
size_t z, x;
for(z = 0; z < Z(); z++){
for(x = 0; x <= z; x++){
p[x * Z() + z] = temp[z * X() + x]; //copy to the destination frame
}
}
}
//load a frame y into a pre-allocated double-precision array
int read_plane_xzd(double* f, size_t y){
|
ebd02452
David Mayerich
added cuBLAS erro...
|
247
|
size_t XB = X() * Z();
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
|
T* temp = (T*) malloc(XB * sizeof(T)); //create a temporary location to store the plane at current precision
if(!read_plane_y(temp, y)) return 1; //read the plane in its native format, if it fails return a 1
for(size_t i = 0; i < XB; i++) f[i] = temp[i]; //convert the plane to a double
return 0;
}
//given a Y, return ZX slice (transposed such that the spectrum is the leading dimension)
int read_plane_zxd(double* p, size_t y){
T* temp = (T*) malloc(X() * Z() * sizeof(T)); //allocate space to store the temporary xz plane
binary<T>::read_plane_2(temp, y); //load the plane from disk
size_t z, x;
for(z = 0; z < Z(); z++){
for(x = 0; x < X(); x++){
p[x * Z() + z] = (double)temp[z * X() + x]; //copy to the destination frame
}
}
return 0;
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
267
|
|
a23c4132
David Mayerich
Doxygen comments ...
|
268
269
270
271
|
/// Perform baseline correction given a list of baseline points and stores the result in a new BSQ file.
/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
/// @param wls is the list of baseline points based on band labels.
|
798cea97
David Mayerich
added progress ba...
|
272
|
bool baseline(std::string outname, std::vector<double> wls, unsigned char* mask = NULL, bool PROGRESS = false){
|
1a224b6a
David Mayerich
fixed linux compa...
|
273
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
274
|
unsigned long long N = wls.size(); //get the number of baseline points
|
1a224b6a
David Mayerich
fixed linux compa...
|
275
|
|
c25e7d0d
heziqi
speed of bip.base...
|
276
|
std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
|
1a224b6a
David Mayerich
fixed linux compa...
|
277
278
|
std::string headername = outname + ".hdr"; //the header file name
|
c25e7d0d
heziqi
speed of bip.base...
|
279
|
//simplify image resolution
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
280
281
282
|
unsigned long long ZX = Z() * X(); //calculate the number of points in a Y slice
unsigned long long L = ZX * sizeof(T); //calculate the number of bytes of a Y slice
unsigned long long B = Z();
|
724ec347
David Mayerich
simplified the EN...
|
283
|
|
c25e7d0d
heziqi
speed of bip.base...
|
284
285
|
T* c; //pointer to the current Y slice
c = (T*)malloc(L); //memory allocation
|
1a224b6a
David Mayerich
fixed linux compa...
|
286
|
|
c25e7d0d
heziqi
speed of bip.base...
|
287
288
|
T* a; //pointer to the two YZ lines surrounding the current YZ line
T* b;
|
1a224b6a
David Mayerich
fixed linux compa...
|
289
|
|
724ec347
David Mayerich
simplified the EN...
|
290
291
|
a = (T*)malloc(X() * sizeof(T));
b = (T*)malloc(X() * sizeof(T));
|
c25e7d0d
heziqi
speed of bip.base...
|
292
293
294
295
|
double ai, bi; //stores the two baseline points wavelength surrounding the current band
double ci; //stores the current band's wavelength
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
296
|
unsigned long long control;
|
c25e7d0d
heziqi
speed of bip.base...
|
297
298
299
300
301
302
|
if (a == NULL || b == NULL || c == NULL){
std::cout<<"ERROR: error allocating memory";
exit(1);
}
// loop start correct every y slice
|
1a224b6a
David Mayerich
fixed linux compa...
|
303
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
304
|
for (unsigned long long k =0; k < Y(); k++)
|
c25e7d0d
heziqi
speed of bip.base...
|
305
306
|
{
//get the current y slice
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
307
|
read_plane_xz(c, k);
|
1a224b6a
David Mayerich
fixed linux compa...
|
308
309
|
//initialize lownum, highnum, low, high
|
6708cc25
heziqi
Ziqi added envi c...
|
310
|
ai = w[0];
|
f6169dea
heziqi
Ziqi completed bi...
|
311
|
control=0;
|
c25e7d0d
heziqi
speed of bip.base...
|
312
313
|
//if no baseline point is specified at band 0,
//set the baseline point at band 0 to 0
|
6708cc25
heziqi
Ziqi added envi c...
|
314
|
if(wls[0] != w[0]){
|
1a224b6a
David Mayerich
fixed linux compa...
|
315
|
bi = wls[control];
|
724ec347
David Mayerich
simplified the EN...
|
316
|
memset(a, (char)0, X() * sizeof(T) );
|
c25e7d0d
heziqi
speed of bip.base...
|
317
318
319
320
|
}
//else get the low band
else{
control++;
|
724ec347
David Mayerich
simplified the EN...
|
321
|
read_x_from_xz(a, c, ai);
|
c25e7d0d
heziqi
speed of bip.base...
|
322
323
324
|
bi = wls[control];
}
//get the high band
|
724ec347
David Mayerich
simplified the EN...
|
325
|
read_x_from_xz(b, c, bi);
|
1a224b6a
David Mayerich
fixed linux compa...
|
326
|
|
c25e7d0d
heziqi
speed of bip.base...
|
327
|
//correct every YZ line
|
1a224b6a
David Mayerich
fixed linux compa...
|
328
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
329
|
for(unsigned long long cii = 0; cii < B; cii++){
|
c25e7d0d
heziqi
speed of bip.base...
|
330
331
|
//update baseline points, if necessary
|
6708cc25
heziqi
Ziqi added envi c...
|
332
|
if( w[cii] >= bi && cii != B - 1) {
|
c25e7d0d
heziqi
speed of bip.base...
|
333
334
|
//if the high band is now on the last BL point
if (control != N-1) {
|
1a224b6a
David Mayerich
fixed linux compa...
|
335
|
|
c25e7d0d
heziqi
speed of bip.base...
|
336
|
control++; //increment the index
|
1a224b6a
David Mayerich
fixed linux compa...
|
337
|
|
c25e7d0d
heziqi
speed of bip.base...
|
338
|
std::swap(a, b); //swap the baseline band pointers
|
1a224b6a
David Mayerich
fixed linux compa...
|
339
|
|
c25e7d0d
heziqi
speed of bip.base...
|
340
341
|
ai = bi;
bi = wls[control];
|
724ec347
David Mayerich
simplified the EN...
|
342
|
read_x_from_xz(b, c, bi);
|
1a224b6a
David Mayerich
fixed linux compa...
|
343
|
|
c25e7d0d
heziqi
speed of bip.base...
|
344
345
|
}
//if the last BL point on the last band of the file?
|
6708cc25
heziqi
Ziqi added envi c...
|
346
|
else if ( wls[control] < w[B - 1]) {
|
1a224b6a
David Mayerich
fixed linux compa...
|
347
|
|
c25e7d0d
heziqi
speed of bip.base...
|
348
|
std::swap(a, b); //swap the baseline band pointers
|
1a224b6a
David Mayerich
fixed linux compa...
|
349
|
|
724ec347
David Mayerich
simplified the EN...
|
350
|
memset(b, (char)0, X() * sizeof(T) ); //clear the high band
|
1a224b6a
David Mayerich
fixed linux compa...
|
351
|
|
c25e7d0d
heziqi
speed of bip.base...
|
352
|
ai = bi;
|
6708cc25
heziqi
Ziqi added envi c...
|
353
|
bi = w[B - 1];
|
c25e7d0d
heziqi
speed of bip.base...
|
354
355
|
}
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
356
|
|
6708cc25
heziqi
Ziqi added envi c...
|
357
|
ci = w[cii];
|
1a224b6a
David Mayerich
fixed linux compa...
|
358
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
359
|
unsigned long long jump = cii * X();
|
c25e7d0d
heziqi
speed of bip.base...
|
360
|
//perform the baseline correction
|
724ec347
David Mayerich
simplified the EN...
|
361
|
for(unsigned i=0; i < X(); i++)
|
c25e7d0d
heziqi
speed of bip.base...
|
362
363
|
{
double r = (double) (ci - ai) / (double) (bi - ai);
|
92e4cf05
heziqi
include file fixed
|
364
|
c[i + jump] =(T) ( c[i + jump] - (b[i] - a[i]) * r - a[i] );
|
c25e7d0d
heziqi
speed of bip.base...
|
365
|
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
366
367
368
|
}//loop for YZ line end
target.write(reinterpret_cast<const char*>(c), L); //write the corrected data into destination
|
cbce396e
David Mayerich
added support for...
|
369
|
|
63fc1df8
David Mayerich
added an inverse ...
|
370
|
if(PROGRESS) progress = (double)(k+1) / Y() * 100;
|
1a224b6a
David Mayerich
fixed linux compa...
|
371
372
|
}//loop for Y slice end
|
c25e7d0d
heziqi
speed of bip.base...
|
373
374
375
376
|
free(a);
free(b);
free(c);
target.close();
|
cbce396e
David Mayerich
added support for...
|
377
|
|
1a224b6a
David Mayerich
fixed linux compa...
|
378
379
|
return true;
|
c25e7d0d
heziqi
speed of bip.base...
|
380
|
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
381
|
|
a23c4132
David Mayerich
Doxygen comments ...
|
382
383
384
385
386
|
/// Normalize all spectra based on the value of a single band, storing the result in a new BSQ file.
/// @param outname is the name of the output file used to store the resulting baseline-corrected data.
/// @param w is the label specifying the band that the hyperspectral image will be normalized to.
/// @param t is a threshold specified such that a spectrum with a value at w less than t is set to zero. Setting this threshold allows the user to limit division by extremely small numbers.
|
c30dd3e6
David Mayerich
added vector norm...
|
387
|
bool ratio(std::string outname, double w, unsigned char* mask = NULL, bool PROGRESS = false)
|
c25e7d0d
heziqi
speed of bip.base...
|
388
|
{
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
389
390
391
|
unsigned long long B = Z(); //calculate the number of bands
unsigned long long ZX = Z() * X();
unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
|
1a224b6a
David Mayerich
fixed linux compa...
|
392
|
unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
393
|
unsigned long long L = ZX * sizeof(T);
|
c25e7d0d
heziqi
speed of bip.base...
|
394
395
396
397
398
399
400
401
|
std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
std::string headername = outname + ".hdr"; //the header file name
T * c; //pointer to the current ZX slice
T * b; //pointer to the standard band
b = (T*)malloc( S ); //memory allocation
|
1a224b6a
David Mayerich
fixed linux compa...
|
402
|
c = (T*)malloc( L );
|
c25e7d0d
heziqi
speed of bip.base...
|
403
|
|
e933c14e
David Mayerich
cleaned up the code
|
404
|
band(b, w); //get the certain band into memory
|
c25e7d0d
heziqi
speed of bip.base...
|
405
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
406
|
for(unsigned long long j = 0; j < Y(); j++)
|
c25e7d0d
heziqi
speed of bip.base...
|
407
|
{
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
408
|
read_plane_xz(c, j);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
409
|
for(unsigned long long i = 0; i < B; i++)
|
c25e7d0d
heziqi
speed of bip.base...
|
410
|
{
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
411
|
for(unsigned long long m = 0; m < X(); m++)
|
c25e7d0d
heziqi
speed of bip.base...
|
412
|
{
|
798cea97
David Mayerich
added progress ba...
|
413
|
if( mask != NULL && !mask[m + j * X()] )
|
724ec347
David Mayerich
simplified the EN...
|
414
|
c[m + i * X()] = (T)0.0;
|
93503ec6
David Mayerich
allowed masking d...
|
415
|
else
|
724ec347
David Mayerich
simplified the EN...
|
416
|
c[m + i * X()] = c[m + i * X()] / b[m + j * X()];
|
1a224b6a
David Mayerich
fixed linux compa...
|
417
|
}
|
c25e7d0d
heziqi
speed of bip.base...
|
418
419
|
}
target.write(reinterpret_cast<const char*>(c), L); //write normalized data into destination
|
cbce396e
David Mayerich
added support for...
|
420
|
|
63fc1df8
David Mayerich
added an inverse ...
|
421
|
if(PROGRESS) progress = (double)(j+1) / Y() * 100;
|
c25e7d0d
heziqi
speed of bip.base...
|
422
|
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
423
|
|
c25e7d0d
heziqi
speed of bip.base...
|
424
425
426
|
free(b);
free(c);
target.close();
|
c25e7d0d
heziqi
speed of bip.base...
|
427
428
|
return true;
}
|
c30dd3e6
David Mayerich
added vector norm...
|
429
|
|
904614c3
David Mayerich
added normalizati...
|
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
|
void normalize(std::string outfile, unsigned char* mask = NULL, bool PROGRESS = false){
std::cout<<"ERROR: algorithm not implemented"<<std::endl;
exit(1);
//This code is almost done but has to be debugged
/*std::ofstream target(outfile.c_str(), std::ios::binary); //open the target binary file
size_t B = Z(); //calculate the number of bands
size_t L = X(); //get the number of items in a line
size_t Lb = L * sizeof(T); //calculate the number of bytes in a line
size_t XY = X() * Y(); //calculate the number of pixels in an image
T* line = (T*) malloc(Lb); //allocate space for a line
T* len = (T*) malloc(Lb); //allocate space for the lengths in a line
for(size_t y = 0; y < Y(); y++){ //for each line in the image
file.seekg(y * Lb, std::ios::beg); //move the pointer to the current file to the beginning
for(size_t b = 0; b < B; b++){ //for each band in this line
file.read((char*)line, Lb); //read a line of pixels
for(size_t l = 0; l < L; l++) //for each element in the line for this band
len[l] += line[l] * line[l]; //add the square of the spectral component
}
for(size_t l = 0; l < L; l++) //for each length element in the line
len[l] = sqrt(len[l]); //calculate the square root of the sum of squares
file.seekg(y * Lb, std::ios::beg); //move the pointer to the current file to the beginning
for(size_t b = 0; b < B; b++){ //for each band in this line
file.read((char*)line, Lb); //read a line of pixels
for(size_t l = 0; l < L; l++) //for each element in the line for this band
line[l] /= len[l]; //divide each element by the length
}
target.write((char*)line, Lb); //write the normalized line to the target file
if(PROGRESS) progress = (double)(y+1) / (double)Y();
}*/
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
464
|
|
b31c02eb
David Mayerich
added a --select ...
|
465
466
467
468
469
|
bool select(std::string outfile, std::vector<double> bandlist, unsigned char* mask = NULL, bool PROGRESS = NULL) {
std::cout << "ERROR: select() not implemented for BIL" << std::endl;
exit(1);
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
470
471
472
|
/// Convert the current BIL file to a BSQ file with the specified file name.
/// @param outname is the name of the output BSQ file to be saved to disk.
|
63fc1df8
David Mayerich
added an inverse ...
|
473
|
bool bsq(std::string outname, bool PROGRESS = false)
|
c25e7d0d
heziqi
speed of bip.base...
|
474
|
{
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
475
|
unsigned long long S = X() * Y() * sizeof(T); //calculate the number of bytes in a band
|
1a224b6a
David Mayerich
fixed linux compa...
|
476
|
|
c25e7d0d
heziqi
speed of bip.base...
|
477
478
|
std::ofstream target(outname.c_str(), std::ios::binary);
std::string headername = outname + ".hdr";
|
1a224b6a
David Mayerich
fixed linux compa...
|
479
|
|
c25e7d0d
heziqi
speed of bip.base...
|
480
481
|
T * p; //pointer to the current band
p = (T*)malloc(S);
|
1a224b6a
David Mayerich
fixed linux compa...
|
482
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
483
|
for ( unsigned long long i = 0; i < Z(); i++)
|
1a224b6a
David Mayerich
fixed linux compa...
|
484
|
{
|
c25e7d0d
heziqi
speed of bip.base...
|
485
|
band_index(p, i);
|
cbce396e
David Mayerich
added support for...
|
486
487
|
target.write(reinterpret_cast<const char*>(p), S); //write a band data into target file
|
1a224b6a
David Mayerich
fixed linux compa...
|
488
|
if(PROGRESS) progress = (double)(i+1) / Z() * 100; //store the progress for the current operation
|
c25e7d0d
heziqi
speed of bip.base...
|
489
|
}
|
6708cc25
heziqi
Ziqi added envi c...
|
490
|
|
c25e7d0d
heziqi
speed of bip.base...
|
491
492
493
494
495
|
free(p);
target.close();
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
496
497
498
|
/// Convert the current BIL file to a BIP file with the specified file name.
/// @param outname is the name of the output BIP file to be saved to disk.
|
63fc1df8
David Mayerich
added an inverse ...
|
499
|
bool bip(std::string outname, bool PROGRESS = false)
|
f6169dea
heziqi
Ziqi completed bi...
|
500
|
{
|
c1078e19
David Mayerich
HSIproc bug fixes...
|
501
|
unsigned long long S = X() * Z() * sizeof(T); //calculate the number of bytes in a ZX slice
|
1a224b6a
David Mayerich
fixed linux compa...
|
502
|
|
f6169dea
heziqi
Ziqi completed bi...
|
503
504
|
std::ofstream target(outname.c_str(), std::ios::binary);
std::string headername = outname + ".hdr";
|
1a224b6a
David Mayerich
fixed linux compa...
|
505
|
|
f6169dea
heziqi
Ziqi completed bi...
|
506
507
508
509
|
T * p; //pointer to the current XZ slice for bil file
p = (T*)malloc(S);
T * q; //pointer to the current ZX slice for bip file
q = (T*)malloc(S);
|
1a224b6a
David Mayerich
fixed linux compa...
|
510
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
511
|
for ( unsigned long long i = 0; i < Y(); i++)
|
1a224b6a
David Mayerich
fixed linux compa...
|
512
|
{
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
513
|
read_plane_xz(p, i);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
514
|
for ( unsigned long long k = 0; k < Z(); k++)
|
f6169dea
heziqi
Ziqi completed bi...
|
515
|
{
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
516
517
|
unsigned long long ks = k * X();
for ( unsigned long long j = 0; j < X(); j++)
|
724ec347
David Mayerich
simplified the EN...
|
518
|
q[k + j * Z()] = p[ks + j];
|
cbce396e
David Mayerich
added support for...
|
519
|
|
1a224b6a
David Mayerich
fixed linux compa...
|
520
|
if(PROGRESS) progress = (double)((i+1) * Z() + k+1) / (Z() * Y()) * 100; //store the progress for the current operation
|
f6169dea
heziqi
Ziqi completed bi...
|
521
|
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
522
523
|
target.write(reinterpret_cast<const char*>(q), S); //write a band data into target file
|
f6169dea
heziqi
Ziqi completed bi...
|
524
525
|
}
|
f6169dea
heziqi
Ziqi completed bi...
|
526
527
528
529
530
531
|
free(p);
free(q);
target.close();
return true;
}
|
20c212c0
heziqi
Ziqi added functi...
|
532
|
|
a23c4132
David Mayerich
Doxygen comments ...
|
533
534
535
536
537
538
539
540
|
/// Return a baseline corrected band given two adjacent baseline points and their bands. The result is stored in a pre-allocated array.
/// @param lb is the label value for the left baseline point
/// @param rb is the label value for the right baseline point
/// @param lp is a pointer to an array holding the band image for the left baseline point
/// @param rp is a pointer to an array holding the band image for the right baseline point
/// @param wavelength is the label value for the requested baseline-corrected band
/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
|
20c212c0
heziqi
Ziqi added functi...
|
541
542
|
bool baseline_band(double lb, double rb, T* lp, T* rp, double wavelength, T* result){
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
543
|
unsigned long long XY = X() * Y();
|
20c212c0
heziqi
Ziqi added functi...
|
544
545
546
547
|
band(result, wavelength); //get band
//perform the baseline correction
double r = (double) (wavelength - lb) / (double) (rb - lb);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
548
|
for(unsigned long long i=0; i < XY; i++){
|
20c212c0
heziqi
Ziqi added functi...
|
549
550
551
552
|
result[i] =(T) (result[i] - (rp[i] - lp[i]) * r - lp[i] );
}
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
553
554
555
556
557
558
559
|
/// Return a baseline corrected band given two adjacent baseline points. The result is stored in a pre-allocated array.
/// @param lb is the label value for the left baseline point
/// @param rb is the label value for the right baseline point
/// @param bandwavelength is the label value for the desired baseline-corrected band
/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size.
|
70407ea9
heziqi
Ziqi completed he...
|
560
|
bool height(double lb, double rb, double bandwavelength, T* result){
|
20c212c0
heziqi
Ziqi added functi...
|
561
562
563
|
T* lp;
T* rp;
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
564
565
|
unsigned long long XY = X() * Y();
unsigned long long S = XY * sizeof(T);
|
70407ea9
heziqi
Ziqi completed he...
|
566
567
|
lp = (T*) malloc(S); //memory allocation
rp = (T*) malloc(S);
|
20c212c0
heziqi
Ziqi added functi...
|
568
|
|
70407ea9
heziqi
Ziqi completed he...
|
569
|
band(lp, lb);
|
1a224b6a
David Mayerich
fixed linux compa...
|
570
|
band(rp, rb);
|
20c212c0
heziqi
Ziqi added functi...
|
571
572
573
|
baseline_band(lb, rb, lp, rp, bandwavelength, result);
|
517876d6
heziqi
metrics finished ...
|
574
575
|
free(lp);
free(rp);
|
20c212c0
heziqi
Ziqi added functi...
|
576
577
578
|
return true;
}
|
70407ea9
heziqi
Ziqi completed he...
|
579
|
|
a23c4132
David Mayerich
Doxygen comments ...
|
580
581
582
583
584
585
586
587
|
/// Calculate the area under the spectrum between two specified points and stores the result in a pre-allocated array.
/// @param lb is the label value for the left baseline point
/// @param rb is the label value for the right baseline point
/// @param lab is the label value for the left bound (start of the integration)
/// @param rab is the label value for the right bound (end of the integration)
/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
|
70407ea9
heziqi
Ziqi completed he...
|
588
|
bool area(double lb, double rb, double lab, double rab, T* result){
|
20c212c0
heziqi
Ziqi added functi...
|
589
590
591
592
593
|
T* lp; //left band pointer
T* rp; //right band pointer
T* cur; //current band 1
T* cur2; //current band 2
|
20c212c0
heziqi
Ziqi added functi...
|
594
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
595
596
|
unsigned long long XY = X() * Y();
unsigned long long S = XY * sizeof(T);
|
20c212c0
heziqi
Ziqi added functi...
|
597
598
599
600
601
|
lp = (T*) malloc(S); //memory allocation
rp = (T*) malloc(S);
cur = (T*) malloc(S);
cur2 = (T*) malloc(S);
|
20c212c0
heziqi
Ziqi added functi...
|
602
603
604
605
|
memset(result, (char)0, S);
//find the wavelenght position in the whole band
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
606
607
608
|
unsigned long long n = w.size();
unsigned long long ai = 0; //left bound position
unsigned long long bi = n - 1; //right bound position
|
20c212c0
heziqi
Ziqi added functi...
|
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
|
//to make sure the left and the right bound are in the bandwidth
if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
exit(1);
}
//to make sure rigth bound is bigger than left bound
else if(lb > rb){
std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
exit(1);
}
//get the position of lb and rb
|
70407ea9
heziqi
Ziqi completed he...
|
624
|
while (lab >= w[ai]){
|
20c212c0
heziqi
Ziqi added functi...
|
625
626
|
ai++;
}
|
70407ea9
heziqi
Ziqi completed he...
|
627
|
while (rab <= w[bi]){
|
20c212c0
heziqi
Ziqi added functi...
|
628
629
630
631
632
633
634
|
bi--;
}
band(lp, lb);
band(rp, rb);
//calculate the beginning and the ending part
|
70407ea9
heziqi
Ziqi completed he...
|
635
636
|
baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
baseline_band(lb, rb, lp, rp, w[bi], cur);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
637
638
|
for(unsigned long long j = 0; j < XY; j++){
result[j] += (T)((rab - w[bi]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
|
20c212c0
heziqi
Ziqi added functi...
|
639
|
}
|
70407ea9
heziqi
Ziqi completed he...
|
640
641
|
baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
baseline_band(lb, rb, lp, rp, w[ai], cur);
|
1a224b6a
David Mayerich
fixed linux compa...
|
642
|
for(unsigned long long j = 0; j < XY; j++){
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
643
|
result[j] += (T)((w[ai] - lab) * ((double)cur[j] + (double)cur2[j]) / 2.0);
|
20c212c0
heziqi
Ziqi added functi...
|
644
|
}
|
517876d6
heziqi
metrics finished ...
|
645
|
|
20c212c0
heziqi
Ziqi added functi...
|
646
647
|
//calculate the area
ai++;
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
648
|
for(unsigned long long i = ai; i <= bi ;i++)
|
20c212c0
heziqi
Ziqi added functi...
|
649
650
|
{
baseline_band(lb, rb, lp, rp, w[ai], cur2);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
651
|
for(unsigned long long j = 0; j < XY; j++)
|
20c212c0
heziqi
Ziqi added functi...
|
652
|
{
|
1a224b6a
David Mayerich
fixed linux compa...
|
653
|
result[j] += (T)((w[ai] - w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 2.0);
|
20c212c0
heziqi
Ziqi added functi...
|
654
655
656
|
}
std::swap(cur,cur2); //swap the band pointers
}
|
70407ea9
heziqi
Ziqi completed he...
|
657
|
|
517876d6
heziqi
metrics finished ...
|
658
659
660
661
|
free(lp);
free(rp);
free(cur);
free(cur2);
|
20c212c0
heziqi
Ziqi added functi...
|
662
663
664
|
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
665
666
667
668
669
670
671
672
673
|
/// Compute the ratio of two baseline-corrected peaks. The result is stored in a pre-allocated array.
/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
/// @param rb1 is the label value for the right baseline point for the first peak (numerator)
/// @param pos1 is the label value for the first peak (numerator) position
/// @param lb2 is the label value for the left baseline point for the second peak (denominator)
/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
/// @param pos2 is the label value for the second peak (denominator) position
/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
|
ba51ae6a
David Mayerich
fixed metric calc...
|
674
|
bool ph_to_ph(T* result, double lb1, double rb1, double pos1, double lb2, double rb2, double pos2, unsigned char* mask = NULL){
|
70407ea9
heziqi
Ziqi completed he...
|
675
|
|
724ec347
David Mayerich
simplified the EN...
|
676
677
|
T* p1 = (T*)malloc(X() * Y() * sizeof(T));
T* p2 = (T*)malloc(X() * Y() * sizeof(T));
|
1a224b6a
David Mayerich
fixed linux compa...
|
678
|
|
70407ea9
heziqi
Ziqi completed he...
|
679
680
681
682
|
//get the two peak band
height(lb1, rb1, pos1, p1);
height(lb2, rb2, pos2, p2);
//calculate the ratio in result
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
683
|
for(unsigned long long i = 0; i < X() * Y(); i++){
|
70407ea9
heziqi
Ziqi completed he...
|
684
685
686
687
688
|
if(p1[i] == 0 && p2[i] ==0)
result[i] = 1;
else
result[i] = p1[i] / p2[i];
}
|
517876d6
heziqi
metrics finished ...
|
689
690
691
|
free(p1);
free(p2);
|
1a224b6a
David Mayerich
fixed linux compa...
|
692
|
return true;
|
70407ea9
heziqi
Ziqi completed he...
|
693
|
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
694
|
|
a23c4132
David Mayerich
Doxygen comments ...
|
695
696
697
698
699
700
701
702
703
|
/// Compute the ratio between a peak area and peak height.
/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
/// @param rb1 is the label value for the right baseline point for the first peak (numerator)
/// @param pos1 is the label value for the first peak (numerator) position
/// @param lb2 is the label value for the left baseline point for the second peak (denominator)
/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
/// @param pos2 is the label value for the second peak (denominator) position
/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
|
ba51ae6a
David Mayerich
fixed metric calc...
|
704
705
|
bool pa_to_ph(T* result, double lb1, double rb1, double lab1, double rab1,
double lb2, double rb2, double pos, unsigned char* mask = NULL){
|
1a224b6a
David Mayerich
fixed linux compa...
|
706
|
|
724ec347
David Mayerich
simplified the EN...
|
707
708
|
T* p1 = (T*)malloc(X() * Y() * sizeof(T));
T* p2 = (T*)malloc(X() * Y() * sizeof(T));
|
1a224b6a
David Mayerich
fixed linux compa...
|
709
|
|
70407ea9
heziqi
Ziqi completed he...
|
710
711
712
713
|
//get the area and the peak band
area(lb1, rb1, lab1, rab1, p1);
height(lb2, rb2, pos, p2);
//calculate the ratio in result
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
714
|
for(unsigned long long i = 0; i < X() * Y(); i++){
|
70407ea9
heziqi
Ziqi completed he...
|
715
716
717
718
719
|
if(p1[i] == 0 && p2[i] ==0)
result[i] = 1;
else
result[i] = p1[i] / p2[i];
}
|
517876d6
heziqi
metrics finished ...
|
720
721
722
|
free(p1);
free(p2);
|
1a224b6a
David Mayerich
fixed linux compa...
|
723
724
725
|
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
726
727
728
729
730
731
732
733
734
|
/// Compute the ratio between two peak areas.
/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
/// @param rb1 is the label value for the right baseline point for the first peak (numerator)
/// @param lab1 is the label value for the left bound (start of the integration) of the first peak (numerator)
/// @param rab1 is the label value for the right bound (end of the integration) of the first peak (numerator)
/// @param lb2 is the label value for the left baseline point for the second peak (denominator)
/// @param rb2 is the label value for the right baseline point for the second peak (denominator)
/// @param lab2 is the label value for the left bound (start of the integration) of the second peak (denominator)
|
1a224b6a
David Mayerich
fixed linux compa...
|
735
|
/// @param rab2 is the label value for the right bound (end of the integration) of the second peak (denominator)
|
a23c4132
David Mayerich
Doxygen comments ...
|
736
|
/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
|
ba51ae6a
David Mayerich
fixed metric calc...
|
737
738
|
bool pa_to_pa(T* result, double lb1, double rb1, double lab1, double rab1,
double lb2, double rb2, double lab2, double rab2, unsigned char* mask = NULL){
|
1a224b6a
David Mayerich
fixed linux compa...
|
739
|
|
724ec347
David Mayerich
simplified the EN...
|
740
741
|
T* p1 = (T*)malloc(X() * Y() * sizeof(T));
T* p2 = (T*)malloc(X() * Y() * sizeof(T));
|
1a224b6a
David Mayerich
fixed linux compa...
|
742
|
|
70407ea9
heziqi
Ziqi completed he...
|
743
744
745
746
|
//get the area and the peak band
area(lb1, rb1, lab1, rab1, p1);
area(lb2, rb2, lab2, rab2, p2);
//calculate the ratio in result
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
747
|
for(unsigned long long i = 0; i < X() * Y(); i++){
|
70407ea9
heziqi
Ziqi completed he...
|
748
749
750
751
752
|
if(p1[i] == 0 && p2[i] ==0)
result[i] = 1;
else
result[i] = p1[i] / p2[i];
}
|
517876d6
heziqi
metrics finished ...
|
753
754
755
|
free(p1);
free(p2);
|
1a224b6a
David Mayerich
fixed linux compa...
|
756
757
|
return true;
}
|
70407ea9
heziqi
Ziqi completed he...
|
758
|
|
a23c4132
David Mayerich
Doxygen comments ...
|
759
760
761
762
763
764
765
|
/// Compute the definite integral of a baseline corrected peak.
/// @param lb is the label value for the left baseline point
/// @param rb is the label value for the right baseline point
/// @param lab is the label for the start of the definite integral
/// @param rab is the label for the end of the definite integral
/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
|
517876d6
heziqi
metrics finished ...
|
766
767
768
769
770
771
|
bool x_area(double lb, double rb, double lab, double rab, T* result){
T* lp; //left band pointer
T* rp; //right band pointer
T* cur; //current band 1
T* cur2; //current band 2
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
772
773
|
unsigned long long XY = X() * Y();
unsigned long long S = XY * sizeof(T);
|
517876d6
heziqi
metrics finished ...
|
774
775
776
777
778
779
780
781
782
|
lp = (T*) malloc(S); //memory allocation
rp = (T*) malloc(S);
cur = (T*) malloc(S);
cur2 = (T*) malloc(S);
memset(result, (char)0, S);
//find the wavelenght position in the whole band
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
783
784
785
|
unsigned long long n = w.size();
unsigned long long ai = 0; //left bound position
unsigned long long bi = n - 1; //right bound position
|
517876d6
heziqi
metrics finished ...
|
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
|
//to make sure the left and the right bound are in the bandwidth
if (lb < w[0] || rb < w[0] || lb > w[n-1] || rb >w[n-1]){
std::cout<<"ERROR: left bound or right bound out of bandwidth"<<std::endl;
exit(1);
}
//to make sure rigth bound is bigger than left bound
else if(lb > rb){
std::cout<<"ERROR: right bound should be bigger than left bound"<<std::endl;
exit(1);
}
//get the position of lb and rb
while (lab >= w[ai]){
ai++;
}
while (rab <= w[bi]){
bi--;
}
band(lp, lb);
band(rp, rb);
//calculate the beginning and the ending part
baseline_band(lb, rb, lp, rp, rab, cur2); //ending part
baseline_band(lb, rb, lp, rp, w[bi], cur);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
812
813
|
for(unsigned long long j = 0; j < XY; j++){
result[j] += (T)((rab - w[bi]) * (rab + w[bi]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
|
517876d6
heziqi
metrics finished ...
|
814
815
816
|
}
baseline_band(lb, rb, lp, rp, lab, cur2); //beginnning part
baseline_band(lb, rb, lp, rp, w[ai], cur);
|
1a224b6a
David Mayerich
fixed linux compa...
|
817
|
for(unsigned long long j = 0; j < XY; j++){
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
818
|
result[j] += (T)((w[ai] - lab) * (w[ai] + lab) * ((double)cur[j] + (double)cur2[j]) / 4.0);
|
517876d6
heziqi
metrics finished ...
|
819
820
821
822
|
}
//calculate f(x) times x
ai++;
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
823
|
for(unsigned long long i = ai; i <= bi ;i++)
|
517876d6
heziqi
metrics finished ...
|
824
825
|
{
baseline_band(lb, rb, lp, rp, w[ai], cur2);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
826
|
for(unsigned long long j = 0; j < XY; j++)
|
517876d6
heziqi
metrics finished ...
|
827
|
{
|
1a224b6a
David Mayerich
fixed linux compa...
|
828
|
result[j] += (T)((w[ai] - w[ai-1]) * (w[ai] + w[ai-1]) * ((double)cur[j] + (double)cur2[j]) / 4.0);
|
517876d6
heziqi
metrics finished ...
|
829
830
831
832
833
834
835
836
837
838
839
|
}
std::swap(cur,cur2); //swap the band pointers
}
free(lp);
free(rp);
free(cur);
free(cur2);
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
840
841
842
843
844
845
846
|
/// Compute the centroid of a baseline corrected peak.
/// @param lb is the label value for the left baseline point
/// @param rb is the label value for the right baseline point
/// @param lab is the label for the start of the peak
/// @param rab is the label for the end of the peak
/// @param result is a pointer to a pre-allocated array at least X * Y * sizeof(T) in size
|
ba51ae6a
David Mayerich
fixed metric calc...
|
847
|
bool centroid(T* result, double lb, double rb, double lab, double rab, unsigned char* mask = NULL){
|
724ec347
David Mayerich
simplified the EN...
|
848
849
|
T* p1 = (T*)malloc(X() * Y() * sizeof(T));
T* p2 = (T*)malloc(X() * Y() * sizeof(T));
|
1a224b6a
David Mayerich
fixed linux compa...
|
850
|
|
517876d6
heziqi
metrics finished ...
|
851
852
853
854
|
//get the area and the peak band
x_area(lb, rb, lab, rab, p1);
area(lb, rb, lab, rab, p2);
//calculate the ratio in result
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
855
|
for(unsigned long long i = 0; i < X() * Y(); i++){
|
ba51ae6a
David Mayerich
fixed metric calc...
|
856
|
if(mask == NULL || mask[i])
|
517876d6
heziqi
metrics finished ...
|
857
858
859
860
861
|
result[i] = p1[i] / p2[i];
}
free(p1);
free(p2);
|
1a224b6a
David Mayerich
fixed linux compa...
|
862
|
return true;
|
517876d6
heziqi
metrics finished ...
|
863
864
|
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
865
866
867
868
869
870
871
|
/// Create a mask based on a given band and threshold value.
/// All pixels in the
/// specified band greater than the threshold are true and all pixels less than the threshold are false.
/// @param mask_band is the band used to specify the mask
/// @param threshold is the threshold used to determine if the mask value is true or false
/// @param p is a pointer to a pre-allocated array at least X * Y in size
|
e37ce7bf
David Mayerich
added the ability...
|
872
|
bool build_mask(unsigned char* out_mask, double mask_band, double lower, double upper, unsigned char* mask = NULL, bool PROGRESS = false){
|
4a6f666c
heziqi
Added mask method
|
873
|
|
724ec347
David Mayerich
simplified the EN...
|
874
|
T* temp = (T*)malloc(X() * Y() * sizeof(T)); //allocate memory for the certain band
|
798cea97
David Mayerich
added progress ba...
|
875
|
band(temp, mask_band, PROGRESS);
|
4a6f666c
heziqi
Added mask method
|
876
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
877
|
for (unsigned long long i = 0; i < X() * Y(); i++) {
|
e37ce7bf
David Mayerich
added the ability...
|
878
879
880
881
882
883
884
|
if(mask == NULL || mask[i] != 0){
if(temp[i] > lower && temp[i] < upper){
out_mask[i] = 255;
}
else
out_mask[i] = 0;
}
|
4a6f666c
heziqi
Added mask method
|
885
886
|
}
|
517876d6
heziqi
metrics finished ...
|
887
|
free(temp);
|
4a6f666c
heziqi
Added mask method
|
888
889
890
|
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
891
892
893
894
|
/// Apply a mask file to the BSQ image, setting all values outside the mask to zero.
/// @param outfile is the name of the masked output file
/// @param p is a pointer to memory of size X * Y, where p(i) = 0 for pixels that will be set to zero.
|
6a46c8ff
David Mayerich
fixed bugs in app...
|
895
|
bool apply_mask(std::string outfile, unsigned char* p, bool PROGRESS = false){
|
740f8cd2
heziqi
added apply_mask
|
896
897
898
|
std::ofstream target(outfile.c_str(), std::ios::binary);
|
1a55a328
David Mayerich
Added comments fo...
|
899
|
//I THINK THIS IS WRONG
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
900
901
|
unsigned long long XZ = X() * Z(); //calculate the number of values in a page on disk
unsigned long long L = XZ * sizeof(T); //calculate the size of the page (in bytes)
|
740f8cd2
heziqi
added apply_mask
|
902
|
|
1a55a328
David Mayerich
Added comments fo...
|
903
|
T * temp = (T*)malloc(L); //allocate memory for a temporary page
|
740f8cd2
heziqi
added apply_mask
|
904
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
905
|
for (unsigned long long i = 0; i < Y(); i++) //for each value in Y() (BIP should be X)
|
740f8cd2
heziqi
added apply_mask
|
906
|
{
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
907
|
read_plane_xz(temp, i); //retrieve an ZX slice, stored in temp
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
908
|
for ( unsigned long long j = 0; j < Z(); j++) //for each Z() (Y)
|
740f8cd2
heziqi
added apply_mask
|
909
|
{
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
910
|
for (unsigned long long k = 0; k < X(); k++) //for each band
|
740f8cd2
heziqi
added apply_mask
|
911
|
{
|
724ec347
David Mayerich
simplified the EN...
|
912
913
|
if(p[i * X() + k] == 0)
temp[j * X() + k] = 0;
|
740f8cd2
heziqi
added apply_mask
|
914
915
916
917
|
else
continue;
}
}
|
1a224b6a
David Mayerich
fixed linux compa...
|
918
|
target.write(reinterpret_cast<const char*>(temp), L); //write a band data into target file
|
6a46c8ff
David Mayerich
fixed bugs in app...
|
919
|
if(PROGRESS) progress = (double)(i+1) / (double)Y() * 100;
|
740f8cd2
heziqi
added apply_mask
|
920
921
922
923
924
925
|
}
target.close();
free(temp);
return true;
}
|
ba51ae6a
David Mayerich
fixed metric calc...
|
926
927
928
929
930
|
/// Copies all spectra corresponding to nonzero values of a mask into a pre-allocated matrix of size (B x P)
/// where P is the number of masked pixels and B is the number of bands. The allocated memory can be accessed
/// using the following indexing: i = p*B + b
/// @param matrix is the destination for the pixel data
/// @param mask is the mask
|
814eb271
David Mayerich
fixed problems wi...
|
931
|
bool sift(T* matrix, unsigned char* mask = NULL, bool PROGRESS = false){
|
ba51ae6a
David Mayerich
fixed metric calc...
|
932
933
934
935
|
size_t Lbytes = sizeof(T) * X();
T* line = (T*) malloc( Lbytes ); //allocate space for a line
file.seekg(0, std::ios::beg); //seek to the beginning of the file
|
1a224b6a
David Mayerich
fixed linux compa...
|
936
|
|
ba51ae6a
David Mayerich
fixed metric calc...
|
937
938
939
940
941
942
943
|
size_t pl;
size_t p = 0; //create counter variables
for(size_t y = 0; y < Y(); y++){ //for each line in the data set
for(size_t b = 0; b < Z(); b++){ //for each band in the data set
pl = 0; //initialize the pixel offset for the current line to zero (0)
file.read( (char*)line, Lbytes ); //read the current line
for(size_t x = 0; x < X(); x++){
|
d4aa4cf7
David Mayerich
fixed the in-memo...
|
944
|
if(mask == NULL || mask[y * X() + x]){ //if the current pixel is masked
|
ba51ae6a
David Mayerich
fixed metric calc...
|
945
946
947
948
949
|
size_t i = (p + pl) * Z() + b; //calculate the index in the sifted matrix
matrix[i] = line[x]; //store the current value in the line at the correct matrix location
pl++; //increment the pixel pointer
}
}
|
814eb271
David Mayerich
fixed problems wi...
|
950
|
if(PROGRESS) progress = (double)( (y+1)*Z() + 1) / (double)(Y() * Z()) * 100;
|
ba51ae6a
David Mayerich
fixed metric calc...
|
951
952
953
954
955
956
|
}
p += pl; //add the line increment to the running pixel index
}
return true;
}
|
f7e99517
David Mayerich
added sifting for...
|
957
958
|
/// Saves to disk only those spectra corresponding to mask values != 0
/// Unlike the BIP and BSQ versions of this function, this version saves a different format: the BIL file is saved as a BIP
|
c1078e19
David Mayerich
HSIproc bug fixes...
|
959
|
bool sift(std::string outfile, unsigned char* p, bool PROGRESS = false){
|
724ec347
David Mayerich
simplified the EN...
|
960
|
// Assume X() = X, Y() = Y, Z() = Z.
|
e843658b
Brad Deutsch
Previous push did...
|
961
962
963
|
std::ofstream target(outfile.c_str(), std::ios::binary);
//for loading pages:
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
964
965
966
|
unsigned long long XZ = X() * Z(); //calculate the number of values in an XZ page on disk
unsigned long long B = Z(); //calculate the number of bands
unsigned long long L = XZ * sizeof(T); //calculate the size of the page (in bytes)
|
e843658b
Brad Deutsch
Previous push did...
|
967
|
|
f7e99517
David Mayerich
added sifting for...
|
968
969
970
971
972
|
//allocate temporary memory for a XZ slice
T* slice = (T*) malloc(L);
//allocate temporary memory for a spectrum
T* spec = (T*) malloc(B * sizeof(T));
|
e843658b
Brad Deutsch
Previous push did...
|
973
|
|
f7e99517
David Mayerich
added sifting for...
|
974
|
//for each slice along the y axis
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
975
|
for (unsigned long long y = 0; y < Y(); y++) //Select a page by choosing Y coordinate, Y()
|
e843658b
Brad Deutsch
Previous push did...
|
976
|
{
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
977
|
read_plane_xz(slice, y); //retrieve an ZX page, store in "slice"
|
f7e99517
David Mayerich
added sifting for...
|
978
|
|
1a224b6a
David Mayerich
fixed linux compa...
|
979
|
//for each sample along X
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
980
|
for (unsigned long long x = 0; x < X(); x++) //Select a pixel by choosing X coordinate in the page, X()
|
e843658b
Brad Deutsch
Previous push did...
|
981
|
{
|
f7e99517
David Mayerich
added sifting for...
|
982
983
|
//if the mask != 0 at that xy pixel
if (p[y * X() + x] != 0) //if the mask != 0 at that XY pixel
|
e843658b
Brad Deutsch
Previous push did...
|
984
|
{
|
f7e99517
David Mayerich
added sifting for...
|
985
|
//for each band at that pixel
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
986
|
for (unsigned long long b = 0; b < B; b++) //Select a voxel by choosing Z coordinate at the pixel
|
e843658b
Brad Deutsch
Previous push did...
|
987
|
{
|
f7e99517
David Mayerich
added sifting for...
|
988
|
spec[b] = slice[b*X() + x]; //Pass the correct spectral value from XZ page into the spectrum to be saved.
|
e843658b
Brad Deutsch
Previous push did...
|
989
|
}
|
f7e99517
David Mayerich
added sifting for...
|
990
|
target.write((char*)spec, B * sizeof(T)); //write that spectrum to disk. Size is L2.
|
e843658b
Brad Deutsch
Previous push did...
|
991
|
}
|
3150c537
David Mayerich
added further pro...
|
992
|
|
c1078e19
David Mayerich
HSIproc bug fixes...
|
993
|
if(PROGRESS) progress = (double) ((y+1) * X() + x+1) / (Y() * X()) * 100;
|
e843658b
Brad Deutsch
Previous push did...
|
994
995
996
|
}
}
target.close();
|
f7e99517
David Mayerich
added sifting for...
|
997
998
|
free(slice);
free(spec);
|
3150c537
David Mayerich
added further pro...
|
999
|
|
e843658b
Brad Deutsch
Previous push did...
|
1000
1001
1002
|
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
1003
1004
1005
|
/// Calculate the mean band value (average along B) at each pixel location.
/// @param p is a pointer to memory of size X * Y * sizeof(T) that will store the band averages.
|
a43c4fe1
heziqi
Added crop in env...
|
1006
|
bool band_avg(T* p){
|
724ec347
David Mayerich
simplified the EN...
|
1007
|
unsigned long long XZ = X() * Z();
|
a43c4fe1
heziqi
Added crop in env...
|
1008
|
T* temp = (T*)malloc(sizeof(T) * XZ);
|
724ec347
David Mayerich
simplified the EN...
|
1009
|
T* line = (T*)malloc(sizeof(T) * X());
|
a43c4fe1
heziqi
Added crop in env...
|
1010
|
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
1011
|
for (unsigned long long i = 0; i < Y(); i++){
|
a43c4fe1
heziqi
Added crop in env...
|
1012
1013
|
getY(temp, i);
//initialize x-line
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
1014
|
for (unsigned long long j = 0; j < X(); j++){
|
a43c4fe1
heziqi
Added crop in env...
|
1015
1016
|
line[j] = 0;
}
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
1017
1018
1019
|
unsigned long long c = 0;
for (unsigned long long j = 0; j < Z(); j++){
for (unsigned long long k = 0; k < X(); k++){
|
724ec347
David Mayerich
simplified the EN...
|
1020
|
line[k] += temp[c] / (T)Z();
|
a43c4fe1
heziqi
Added crop in env...
|
1021
1022
1023
|
c++;
}
}
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
1024
|
for (unsigned long long j = 0; j < X(); j++){
|
724ec347
David Mayerich
simplified the EN...
|
1025
|
p[j + i * X()] = line[j];
|
a43c4fe1
heziqi
Added crop in env...
|
1026
1027
1028
1029
1030
1031
|
}
}
free(temp);
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
1032
1033
1034
1035
|
/// Calculate the mean value for all masked (or valid) pixels in a band and returns the average spectrum
/// @param p is a pointer to pre-allocated memory of size [B * sizeof(T)] that stores the mean spectrum
/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
|
35a7195f
David Mayerich
added median spec...
|
1036
|
bool mean_spectrum(double* m, double* std, unsigned char* mask = NULL, bool PROGRESS = false){
|
724ec347
David Mayerich
simplified the EN...
|
1037
1038
|
unsigned long long XZ = X() * Z();
unsigned long long XY = X() * Y();
|
a43c4fe1
heziqi
Added crop in env...
|
1039
|
T* temp = (T*)malloc(sizeof(T) * XZ);
|
35a7195f
David Mayerich
added median spec...
|
1040
1041
1042
|
memset(m, 0, Z() * sizeof(double)); //initialize the mean to zero
double* e_x2 = (double*)malloc(Z() * sizeof(double)); //allocate space for E[x^2]
memset(e_x2, 0, Z() * sizeof(double)); //initialize E[x^2] to zero
|
a43c4fe1
heziqi
Added crop in env...
|
1043
|
//calculate vaild number in a band
|
35a7195f
David Mayerich
added median spec...
|
1044
1045
1046
|
size_t count = nnz(mask); //count the number of pixels in the mask
double x; //create a register to store the pixel value
|
2aa68315
David Mayerich
major bug fixes f...
|
1047
|
for (unsigned long long k = 0; k < Y(); k++){
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
1048
|
read_plane_xz(temp, k);
|
2aa68315
David Mayerich
major bug fixes f...
|
1049
1050
1051
1052
|
unsigned long long kx = k * X();
for (unsigned long long i = 0; i < X(); i++){
if (mask == NULL || mask[kx + i] != 0){
for (unsigned long long j = 0; j < Z(); j++){
|
35a7195f
David Mayerich
added median spec...
|
1053
1054
1055
|
x = temp[j * X() + i];
m[j] += x / (double)count;
e_x2[j] += x*x / (double)count;
|
a43c4fe1
heziqi
Added crop in env...
|
1056
1057
1058
|
}
}
}
|
63fc1df8
David Mayerich
added an inverse ...
|
1059
|
if(PROGRESS) progress = (double)(k+1) / Y() * 100;
|
a43c4fe1
heziqi
Added crop in env...
|
1060
|
}
|
35a7195f
David Mayerich
added median spec...
|
1061
1062
1063
1064
|
for(size_t i = 0; i < Z(); i++) //calculate the standard deviation
std[i] = sqrt(e_x2[i] - m[i] * m[i]);
|
a43c4fe1
heziqi
Added crop in env...
|
1065
1066
1067
1068
|
free(temp);
return true;
}
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
1069
|
int co_matrix_cublas(double* co, double* avg, unsigned char *mask, bool PROGRESS = false){
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
|
cublasStatus_t stat;
cublasHandle_t handle;
progress = 0; //initialize the progress to zero (0)
size_t XY = X() * Y(); //calculate the number of elements in a band image
size_t XB = X() * Z();
size_t B = Z(); //calculate the number of spectral elements
double* F = (double*)malloc(sizeof(double) * B * X()); //allocate space for the frame that will be pulled from the file
double* F_dev;
HANDLE_ERROR(cudaMalloc(&F_dev, X() * B * sizeof(double))); //allocate space for the frame on the GPU
double* s_dev; //declare a device pointer that will store the spectrum on the GPU
double* A_dev; //declare a device pointer that will store the covariance matrix on the GPU
double* avg_dev; //declare a device pointer that will store the average spectrum
HANDLE_ERROR(cudaMalloc(&s_dev, B * sizeof(double))); //allocate space on the CUDA device for a spectrum
HANDLE_ERROR(cudaMalloc(&A_dev, B * B * sizeof(double))); //allocate space on the CUDA device for the covariance matrix
HANDLE_ERROR(cudaMemset(A_dev, 0, B * B * sizeof(double))); //initialize the covariance matrix to zero (0)
HANDLE_ERROR(cudaMalloc(&avg_dev, XB * sizeof(double))); //allocate space on the CUDA device for the average spectrum
for(size_t x = 0; x < X(); x++) //make multiple copies of the average spectrum in order to build a matrix
HANDLE_ERROR(cudaMemcpy(&avg_dev[x * B], avg, B * sizeof(double), cudaMemcpyHostToDevice));
//stat = cublasSetVector((int)B, sizeof(double), avg, 1, avg_dev, 1); //copy the average spectrum to the CUDA device
double ger_alpha = 1.0/(double)XY; //scale the outer product by the inverse of the number of samples (mean outer product)
double axpy_alpha = -1; //multiplication factor for the average spectrum (in order to perform a subtraction)
CUBLAS_HANDLE_ERROR(stat = cublasCreate(&handle)); //create a cuBLAS instance
if (stat != CUBLAS_STATUS_SUCCESS) return 1; //test the cuBLAS instance to make sure it is valid
else std::cout<<"Using cuBLAS to calculate the mean covariance matrix..."<<std::endl;
double beta = 1.0;
size_t x, y;
for(y = 0; y < Y(); y++){ //for each line
read_plane_zxd(F, y); //read a frame from the file
HANDLE_ERROR(cudaMemcpy(F_dev, F, XB * sizeof(double), cudaMemcpyHostToDevice)); //copy the frame to the GPU
CUBLAS_HANDLE_ERROR(cublasDgeam(handle, CUBLAS_OP_N, CUBLAS_OP_N, (int)B, (int)X(), &axpy_alpha, avg_dev, (int)B, &beta, F_dev, (int)B, F_dev, (int)B));//subtract the mean spectrum
for(x = 0; x < X(); x++)
CUBLAS_HANDLE_ERROR(cublasDsyr(handle, CUBLAS_FILL_MODE_UPPER, (int)B, &ger_alpha, &F_dev[x*B], 1, A_dev, (int)B)); //perform an outer product
if(PROGRESS) progress = (double)(y + 1) / Y() * 100;
}
cublasGetMatrix((int)B, (int)B, sizeof(double), A_dev, (int)B, co, (int)B); //copy the result from the GPU to the CPU
cudaFree(A_dev); //clean up allocated device memory
cudaFree(s_dev);
cudaFree(avg_dev);
for(unsigned long long i = 0; i < B; i++){ //copy the upper triangular portion to the lower triangular portion
for(unsigned long long j = i+1; j < B; j++){
co[B * i + j] = co[B * j + i];
}
}
return 0;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
1130
1131
1132
1133
1134
|
/// Calculate the covariance matrix for all masked pixels in the image.
/// @param co is a pointer to pre-allocated memory of size [B * B] that stores the resulting covariance matrix
/// @param avg is a pointer to memory of size B that stores the average spectrum
/// @param mask is a pointer to memory of size [X * Y] that stores the mask value at each pixel location
|
dd5aab2f
David Mayerich
fixed errors in s...
|
1135
|
bool co_matrix(double* co, double* avg, unsigned char *mask, int cuda_device = 0, bool PROGRESS = false){
|
63fc1df8
David Mayerich
added an inverse ...
|
1136
|
progress = 0;
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
1137
|
|
dd5aab2f
David Mayerich
fixed errors in s...
|
1138
|
if (cuda_device >= 0) { //if a CUDA device is specified
|
dc99699e
David Mayerich
implemented a no-...
|
1139
1140
|
int dev_count;
HANDLE_ERROR(cudaGetDeviceCount(&dev_count)); //get the number of CUDA devices
|
dd5aab2f
David Mayerich
fixed errors in s...
|
1141
|
std::cout << "Number of CUDA devices: " << dev_count << std::endl; //output the number of CUDA devices
|
dc99699e
David Mayerich
implemented a no-...
|
1142
|
cudaDeviceProp prop;
|
dd5aab2f
David Mayerich
fixed errors in s...
|
1143
1144
1145
1146
|
//int best_device_id = 0; //stores the best CUDA device
//float best_device_cc = 0.0f; //stores the compute capability of the best device
std::cout << "CUDA devices----" << std::endl;
for (int d = 0; d < dev_count; d++) { //for each CUDA device
|
dc99699e
David Mayerich
implemented a no-...
|
1147
|
cudaGetDeviceProperties(&prop, d); //get the property of the first device
|
dd5aab2f
David Mayerich
fixed errors in s...
|
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
|
//float cc = prop.major + prop.minor / 10.0f; //calculate the compute capability
std::cout << d << ": [" << prop.major << "." << prop.minor << "] " << prop.name << std::endl; //display the device information
//if(cc > best_device_cc){
// best_device_cc = cc; //if this is better than the previous device, use it
// best_device_id = d;
//}
}
if (dev_count > 0 && dev_count > cuda_device) { //if the first device is not an emulator
cudaGetDeviceProperties(&prop, cuda_device); //get the property of the requested CUDA device
if (prop.major != 9999) {
std::cout << "Using device " << cuda_device << std::endl;
HANDLE_ERROR(cudaSetDevice(cuda_device));
int status = co_matrix_cublas(co, avg, mask, PROGRESS); //use cuBLAS to calculate the covariance matrix
if (status == 0) return true; //if the cuBLAS function returned correctly, we're done
|
dc99699e
David Mayerich
implemented a no-...
|
1162
|
}
|
dd5aab2f
David Mayerich
fixed errors in s...
|
1163
|
}
|
dc99699e
David Mayerich
implemented a no-...
|
1164
|
}
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
1165
|
|
a43c4fe1
heziqi
Added crop in env...
|
1166
|
//memory allocation
|
724ec347
David Mayerich
simplified the EN...
|
1167
|
unsigned long long xy = X() * Y();
|
2aa68315
David Mayerich
major bug fixes f...
|
1168
|
unsigned long long B = Z();
|
a43c4fe1
heziqi
Added crop in env...
|
1169
1170
|
T* temp = (T*)malloc(sizeof(T) * B);
//count vaild pixels in a band
|
2aa68315
David Mayerich
major bug fixes f...
|
1171
1172
1173
|
unsigned long long count = 0;
for (unsigned long long j = 0; j < xy; j++){
if (mask == NULL || mask[j] != 0){
|
a43c4fe1
heziqi
Added crop in env...
|
1174
1175
1176
1177
|
count++;
}
}
//initialize correlation matrix
|
2aa68315
David Mayerich
major bug fixes f...
|
1178
1179
|
for (unsigned long long i = 0; i < B; i++){
for (unsigned long long k = 0; k < B; k++){
|
a43c4fe1
heziqi
Added crop in env...
|
1180
1181
1182
1183
|
co[i * B + k] = 0;
}
}
//calculate correlation coefficient matrix
|
2aa68315
David Mayerich
major bug fixes f...
|
1184
1185
|
for (unsigned long long j = 0; j < xy; j++){
if (mask == NULL || mask[j] != 0){
|
a43c4fe1
heziqi
Added crop in env...
|
1186
|
pixel(temp, j);
|
2aa68315
David Mayerich
major bug fixes f...
|
1187
1188
1189
|
for (unsigned long long i = 0; i < B; i++){
for (unsigned long long k = i; k < B; k++){
co[i * B + k] += ((double)temp[i] - (double)avg[i]) * ((double)temp[k] - (double)avg[k]) / (double)count;
|
a43c4fe1
heziqi
Added crop in env...
|
1190
1191
1192
|
}
}
}
|
63fc1df8
David Mayerich
added an inverse ...
|
1193
|
if(PROGRESS) progress = (double)(j+1) / xy * 100;
|
a43c4fe1
heziqi
Added crop in env...
|
1194
1195
|
}
//because correlation matrix is symmetric
|
2aa68315
David Mayerich
major bug fixes f...
|
1196
1197
|
for (unsigned long long i = 0; i < B; i++){
for (unsigned long long k = i + 1; k < B; k++){
|
a43c4fe1
heziqi
Added crop in env...
|
1198
1199
1200
1201
1202
|
co[k * B + i] = co[i * B + k];
}
}
free(temp);
|
a43c4fe1
heziqi
Added crop in env...
|
1203
1204
1205
|
return true;
}
|
0df38ff3
heziqi
fixed head detached
|
1206
|
|
a23c4132
David Mayerich
Doxygen comments ...
|
1207
1208
1209
1210
1211
1212
1213
|
/// Crop a region of the image and save it to a new file.
/// @param outfile is the file name for the new cropped image
/// @param x0 is the lower-left x pixel coordinate to be included in the cropped image
/// @param y0 is the lower-left y pixel coordinate to be included in the cropped image
/// @param x1 is the upper-right x pixel coordinate to be included in the cropped image
/// @param y1 is the upper-right y pixel coordinate to be included in the cropped image
|
e9bddc57
David Mayerich
added band croppi...
|
1214
1215
1216
1217
1218
|
bool crop(std::string outfile, unsigned long long x0,
unsigned long long y0,
unsigned long long x1,
unsigned long long y1,
unsigned long long b0,
|
63fc1df8
David Mayerich
added an inverse ...
|
1219
1220
|
unsigned long long b1,
bool PROGRESS = false){
|
a43c4fe1
heziqi
Added crop in env...
|
1221
|
|
3033b886
David Mayerich
implemented spect...
|
1222
|
//calculate the new image parameters
|
be1064b1
David Mayerich
fixed cropping er...
|
1223
1224
1225
|
unsigned long long samples = x1 - x0 + 1;
unsigned long long lines = y1 - y0 + 1;
unsigned long long bands = b1 - b0 + 1;
|
3033b886
David Mayerich
implemented spect...
|
1226
1227
1228
1229
1230
1231
1232
1233
1234
|
//calculate the size of a line
unsigned long long L = samples * sizeof(T);
//allocate space for a line
T* temp = (T*)malloc(bands * L);
//create an output stream to store the cropped file
|
a43c4fe1
heziqi
Added crop in env...
|
1235
|
std::ofstream out(outfile.c_str(), std::ios::binary);
|
3033b886
David Mayerich
implemented spect...
|
1236
1237
1238
1239
1240
|
//calculate the distance between bands
unsigned long long jumpb = (X() - samples) * sizeof(T);
//distance needed to jump from the previous line of the last band to the next line of the first band
|
be1064b1
David Mayerich
fixed cropping er...
|
1241
|
unsigned long long longjump = ((Z() - bands) * X()) * sizeof(T);
|
3033b886
David Mayerich
implemented spect...
|
1242
1243
1244
1245
|
//set the start position for the cropped region
file.seekg((y0 * X() * Z() + b0 * X() + x0) * sizeof(T), std::ios::beg);
|
9d3ba0b1
David Mayerich
added stim::hsi a...
|
1246
|
for (unsigned long long x = 0; x < lines; x++)
|
a43c4fe1
heziqi
Added crop in env...
|
1247
|
{
|
be1064b1
David Mayerich
fixed cropping er...
|
1248
|
for (unsigned long long z = b0; z <= b1; z++)
|
a43c4fe1
heziqi
Added crop in env...
|
1249
|
{
|
3033b886
David Mayerich
implemented spect...
|
1250
|
file.read((char *)(temp + z * samples), sizeof(T) * samples);
|
1a224b6a
David Mayerich
fixed linux compa...
|
1251
|
file.seekg(jumpb, std::ios::cur); //go to the next band
|
3150c537
David Mayerich
added further pro...
|
1252
|
|
be1064b1
David Mayerich
fixed cropping er...
|
1253
|
if(PROGRESS) progress = (double)(x * (b1 - b0 + 1) + z + 1) / (lines * (b1 - b0 + 1)) * 100;
|
a43c4fe1
heziqi
Added crop in env...
|
1254
|
}
|
3033b886
David Mayerich
implemented spect...
|
1255
1256
1257
1258
1259
1260
|
//write slice data into target file
out.write(reinterpret_cast<const char*>(temp), bands * L);
//seek to the beginning of the next X-Z slice
file.seekg(longjump, std::ios::cur);
|
a43c4fe1
heziqi
Added crop in env...
|
1261
|
}
|
3033b886
David Mayerich
implemented spect...
|
1262
1263
|
//free the temporary frame
|
a43c4fe1
heziqi
Added crop in env...
|
1264
|
free(temp);
|
3150c537
David Mayerich
added further pro...
|
1265
|
|
a43c4fe1
heziqi
Added crop in env...
|
1266
1267
1268
|
return true;
}
|
dc8eb8aa
David Mayerich
added band trimmi...
|
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
|
/// Remove a list of bands from the ENVI file
/// @param outfile is the file name for the output hyperspectral image (with trimmed bands)
/// @param b is an array of bands to be eliminated
void trim(std::string outfile, std::vector<size_t> band_array, bool PROGRESS = false){
std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
file.seekg(0, std::ios::beg); //move to the beginning of the input file
size_t Xb = X() * sizeof(T); //calculate the number of bytes in a line
T* line = (T*)malloc(Xb); //allocate space for a line
size_t i; //create an index into the band array
for(size_t y = 0; y < Y(); y++){ //for each Y plane
i = 0;
for(size_t b = 0; b < Z(); b++){ //for each band
if(b != band_array[i]){ //if this band isn't trimmed
file.read((char*)line, Xb); //read a line
out.write((char*)line, Xb); //write the line
}
else{
file.seekg(Xb, std::ios::cur); //if this band is trimmed, skip it
i++;
}
}
if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100;
}
free(line);
}
|
2ce6954b
David Mayerich
added the ability...
|
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
|
/// Combine two BSQ images along the Y axis
/// @param outfile is the combined file to be output
/// @param infile is the input file stream for the image to combine with this one
/// @param Sx is the size of the second image along X
/// @param Sy is the size of the second image along Y
/// @param offset is a shift (negative or positive) in the combined image to the left or right
void combine(std::string outfile, bil<T>* C, long long xp, long long yp, bool PROGRESS = false){
std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
file.seekg(0, std::ios::beg); //move to the beginning of both files
C->file.seekg(0, std::ios::beg);
size_t S[2]; //size of the output band image
size_t p0[2]; //position of the current image in the output
size_t p1[2]; //position of the source image in the output
hsi<T>::calc_combined_size(xp, yp, C->X(), C->Y(), S[0], S[1], p0[0], p0[1], p1[0], p1[1]); //calculate the image placement parameters
size_t line_bytes = X() * sizeof(T);
size_t band_bytes = X() * Y() * sizeof(T);
T* cur = (T*)malloc(line_bytes); //allocate space for a band of the current image
size_t line_src_bytes = C->X() * sizeof(T);
size_t band_src_bytes = C->X() * C->Y() * sizeof(T);
T* src = (T*)malloc(line_src_bytes); //allocate space for a band of the source image
size_t line_dst_bytes = S[0] * sizeof(T);
size_t band_dst_bytes = S[0] * S[1] * sizeof(T);
T* dst = (T*)malloc(line_dst_bytes); //allocate space for a band of the destination image
for(size_t y = 0; y < S[1]; y++){ //for each line in the destination file
memset(dst, 0, line_dst_bytes); //set all values to zero (0) in the destination image
for(size_t b = 0; b < Z(); b++){ //for each band in both images
if(y >= p0[1] && y < p0[1] + Y()){ //if the current image crosses this line
file.read((char*)cur, line_bytes); //read a line from the current image
memcpy(&dst[p0[0]], cur, line_bytes); //copy the current line to the correct spot in the destination line
}
if(y >= p1[1] && y < p1[1] + C->Y()){ //if the source image crosses this line
C->file.read((char*)src, line_src_bytes); //read a line from the source image
memcpy(&dst[p1[0]], src, line_src_bytes); //copy the source line into the correct spot in the destination line
}
out.write((char*)dst, line_dst_bytes);
if(PROGRESS) progress = (double)((y + 1)*Z() + b + 1) / (double) (Z() * S[1]) * 100;
}
}
}
|
82a7c08c
David Mayerich
added an append c...
|
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
|
///Append two files together along the band dimension
void append(std::string outfile, bil<T>* C, bool PROGRESS = false) {
std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
file.seekg(0, std::ios::beg); //move to the beginning of both files
C->file.seekg(0, std::ios::beg);
size_t a_bytes = X() * Z() * sizeof(T); //calculate the number of bytes in a single plane of this file
size_t b_bytes = C->X() * C->Z() * sizeof(T); //calculate the number of bytes in a single plane of the appending file
T* a = (T*)malloc(a_bytes); //allocate space for a plane of the current file
T* b = (T*)malloc(b_bytes); //allocate space for a plane of the appended file
if (PROGRESS) progress = 0;
for (size_t y = 0; y < Y(); y++) {
read_plane_xz(a, y); //read a plane from the current file
out.write((char*)a, a_bytes); //write the plane to disk
C->read_plane_xz(b, y); //read a plane from the appending file
out.write((char*)b, b_bytes);
if (PROGRESS) progress = (double)(y + 1) / (double)(Y()) * 100;
}
out.close();
}
|
9b2a5d71
David Mayerich
implemented finit...
|
1366
|
/// Convolve the given band range with a kernel specified by a vector of coefficients.
|
2ce6954b
David Mayerich
added the ability...
|
1367
|
|
9b2a5d71
David Mayerich
implemented finit...
|
1368
1369
1370
1371
1372
1373
|
/// @param outfile is an already open stream to the output file
/// @param C is an array of coefficients
/// @param start is the band to start processing (the first coefficient starts here)
/// @param nbands is the number of bands to process
/// @param center is the index for the center coefficient for the kernel (used to set the wavelengths in the output file)
|
9d77bbd9
David Mayerich
updates for HSIvi...
|
1374
|
void convolve(std::string outfile, std::vector<double> C, size_t start, size_t end, unsigned char* mask = NULL, bool PROGRESS = false){
|
9b2a5d71
David Mayerich
implemented finit...
|
1375
1376
1377
1378
1379
1380
1381
1382
|
std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
size_t B = end - start + 1;
size_t Xb = X() * sizeof(T); //calculate the size of a band (frame) in bytes
size_t XBb = Xb * Z();
size_t N = C.size(); //get the number of bands that the kernel spans
std::deque<T*> frame(N, NULL); //create an array to store pointers to each frame
for(size_t f = 0; f < N; f++)
|
1a224b6a
David Mayerich
fixed linux compa...
|
1383
|
frame[f] = (T*)malloc(Xb); //allocate space for the frame
|
9b2a5d71
David Mayerich
implemented finit...
|
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
|
T* outline = (T*)malloc(Xb); //allocate space for the output band
//Implementation: In order to minimize reads from secondary storage, each band is only loaded once into the 'frame' deque.
// When a new band is loaded, the last band is popped, a new frame is copied to the pointer, and it is
// re-inserted into the deque.
for(size_t y = 0; y < Y(); y++){
file.seekg(y * XBb + start * Xb, std::ios::beg); //move to the beginning of the 'start' band
for(size_t f = 0; f < N; f++) //for each frame
file.read((char*)frame[f], Xb); //load the frame
for(size_t b = 0; b < B; b++){ //for each band
memset(outline, 0, Xb); //set the output band to zero (0)
for(size_t c = 0; c < N; c++){ //for each frame (corresponding to each coefficient)
for(size_t x = 0; x < X(); x++){ //for each pixel
|
9d77bbd9
David Mayerich
updates for HSIvi...
|
1398
1399
1400
|
if(mask == NULL || mask[y * X() + x]){
outline[x] += (T)(C[c] * frame[c][x]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
}
|
9b2a5d71
David Mayerich
implemented finit...
|
1401
1402
1403
1404
1405
1406
1407
1408
|
}
}
out.write((char*)outline, Xb); //output the band
file.read((char*)frame[0], Xb); //read the next band
frame.push_back(frame.front()); //put the first element in the back
frame.pop_front(); //pop the first element
}
if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100;
|
1a224b6a
David Mayerich
fixed linux compa...
|
1409
|
}
|
9b2a5d71
David Mayerich
implemented finit...
|
1410
1411
1412
|
}
/// Approximate the spectral derivative of the image
|
9d77bbd9
David Mayerich
updates for HSIvi...
|
1413
|
void deriv(std::string outfile, size_t d, size_t order, const std::vector<double> w = std::vector<double>(), unsigned char* mask = NULL, bool PROGRESS = false){
|
9b2a5d71
David Mayerich
implemented finit...
|
1414
|
std::ofstream out(outfile.c_str(), std::ios::binary); //open the output file for writing
|
1a224b6a
David Mayerich
fixed linux compa...
|
1415
|
|
9b2a5d71
David Mayerich
implemented finit...
|
1416
1417
1418
1419
|
size_t Xb = X() * sizeof(T); //calculate the size of a line (frame) in bytes
size_t XBb = Xb * Z();
size_t B = Z();
|
1a224b6a
David Mayerich
fixed linux compa...
|
1420
|
file.seekg(0, std::ios::beg); //seek to the beginning of the file
|
9b2a5d71
David Mayerich
implemented finit...
|
1421
1422
1423
1424
|
size_t N = order + d; //approximating a derivative requires order + d samples
std::deque<T*> frame(N, NULL); //create an array to store pointers to each frame
for(size_t f = 0; f < N; f++) //for each frame
|
1a224b6a
David Mayerich
fixed linux compa...
|
1425
|
frame[f] = (T*)malloc(Xb); //allocate space for the frame
|
9b2a5d71
David Mayerich
implemented finit...
|
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
|
T* outline = (T*)malloc(Xb); //allocate space for the output band
//Implementation: In order to minimize reads from secondary storage, each band is only loaded once into the 'frame' deque.
// When a new band is loaded, the last band is popped, a new frame is copied to the pointer, and it is
// re-inserted into the deque.
size_t mid = (size_t)(N / 2); //calculate the mid point of the kernel
size_t iw; //index to the first wavelength used to evaluate the derivative at this band
for(size_t y = 0; y < Y(); y++){
//file.seekg(y * XBb, std::ios::beg); //seek to the beginning of the current Y plane
for(size_t f = 0; f < N; f++) //for each frame
file.read((char*)frame[f], Xb); //load a line
for(size_t b = 0; b < B; b++){ //for each band
if(b < mid) //calculate the first wavelength used to evaluate the derivative at this band
iw = 0;
else if(b > B - (N - mid + 1))
iw = B - N;
else{
iw = b - mid;
file.read((char*)frame[0], Xb); //read the next band
frame.push_back(frame.front()); //put the first element in the back
frame.pop_front(); //pop the first element
}
std::vector<double> w_pts(w.begin() + iw, w.begin() + iw + N); //get the wavelengths corresponding to each sample
std::vector<double> C = diff_coefficients(w[b], w_pts, d); //get the optimal sample weights
memset(outline, 0, Xb); //set the output band to zero (0)
for(size_t c = 0; c < N; c++){ //for each frame (corresponding to each coefficient)
for(size_t x = 0; x < X(); x++){ //for each pixel
|
9d77bbd9
David Mayerich
updates for HSIvi...
|
1457
1458
1459
|
if(mask == NULL || mask[y * X() + x]){
outline[x] += (T)(C[c] * frame[c][x]); //calculate the contribution of the current frame (scaled by the corresponding coefficient)
}
|
9b2a5d71
David Mayerich
implemented finit...
|
1460
1461
1462
|
}
}
out.write((char*)outline, Xb); //output the band
|
1a224b6a
David Mayerich
fixed linux compa...
|
1463
|
|
9b2a5d71
David Mayerich
implemented finit...
|
1464
1465
1466
1467
1468
|
}
if(PROGRESS) progress = (double)(y+1) / (double)Y() * 100;
}
}
|
0df38ff3
heziqi
fixed head detached
|
1469
|
|
278c622b
David Mayerich
finished Windows ...
|
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
|
bool multiply(std::string outname, double v, unsigned char* mask = NULL, bool PROGRESS = false){
unsigned long long B = Z(); //calculate the number of bands
unsigned long long ZX = Z() * X();
unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band
unsigned long long L = ZX * sizeof(T);
std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
T * c; //pointer to the current ZX slice
c = (T*)malloc( L ); //allocate space for the slice
for(unsigned long long j = 0; j < Y(); j++){ //for each line
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
1483
|
read_plane_xz(c, j); //load the line into memory
|
278c622b
David Mayerich
finished Windows ...
|
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
|
for(unsigned long long i = 0; i < B; i++){ //for each band
for(unsigned long long m = 0; m < X(); m++){ //for each sample
if( mask == NULL && mask[m + j * X()] ) //if the pixel is masked
c[m + i * X()] *= (T)v;
}
}
target.write(reinterpret_cast<const char*>(c), L); //write normalized data into destination
if(PROGRESS) progress = (double)(j+1) / Y() * 100; //update the progress
}
free(c); //free the slice memory
target.close(); //close the output file
return true;
}
bool add(std::string outname, double v, unsigned char* mask = NULL, bool PROGRESS = false){
unsigned long long B = Z(); //calculate the number of bands
unsigned long long ZX = Z() * X();
unsigned long long XY = X() * Y(); //calculate the number of pixels in a band
unsigned long long S = XY * sizeof(T); //calculate the number of bytes in a band
unsigned long long L = ZX * sizeof(T);
std::ofstream target(outname.c_str(), std::ios::binary); //open the target binary file
T * c; //pointer to the current ZX slice
c = (T*)malloc( L ); //allocate space for the slice
for(unsigned long long j = 0; j < Y(); j++){ //for each line
|
2781a48c
David Mayerich
fixed cuBLAS func...
|
1513
|
read_plane_xz(c, j); //load the line into memory
|
278c622b
David Mayerich
finished Windows ...
|
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
|
for(unsigned long long i = 0; i < B; i++){ //for each band
for(unsigned long long m = 0; m < X(); m++){ //for each sample
if( mask == NULL && mask[m + j * X()] ) //if the pixel is masked
c[m + i * X()] += (T)v;
}
}
target.write(reinterpret_cast<const char*>(c), L); //write normalized data into destination
if(PROGRESS) progress = (double)(j+1) / Y() * 100; //update the progress
}
free(c); //free the slice memory
target.close(); //close the output file
return true;
}
|
a23c4132
David Mayerich
Doxygen comments ...
|
1530
|
/// Close the file.
|
f6169dea
heziqi
Ziqi completed bi...
|
1531
1532
1533
1534
1535
|
bool close(){
file.close();
return true;
}
|
c25e7d0d
heziqi
speed of bip.base...
|
1536
|
};
|
cac62fd3
David Mayerich
modified to work ...
|
1537
|
}
|
6aa04ba2
David Mayerich
interleave types
|
1538
1539
|
#endif
|