codebase / stimlib

Commit 6e257ab383f6bc6490c95efd435b8f5156db6667

Authored by dmayerich
1 parent cc8715ef

ENVI and colormap additions

Showing 5 changed files with 966 additions and 418 deletions   Show diff stats
Inline Side-by-side
rts/cuda/error.h
  Show/Hide comments   View file @6e257ab
1 1 #include <stdio.h>
2 2 #include "cuda_runtime.h"
3   -#include "device_launch_parameters.h"
  3 +#include "device_launch_parameters.h"
  4 +#include "cufft.h"
4 5  
5 6 #ifndef CUDA_HANDLE_ERROR_H
6 7 #define CUDA_HANDLE_ERROR_H
... ... @@ -13,10 +14,28 @@ static void HandleError( cudaError_t err, const char *file, int line ) {
13 14 //fclose(outfile);
14 15 printf("%s in %s at line %d\n", cudaGetErrorString( err ), file, line );
15 16 //exit( EXIT_FAILURE );
16   -
  17 +
17 18 }
18 19 }
19   -#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))
  20 +#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))
  21 +
  22 +static void CufftError( cufftResult err )
  23 +{
  24 + if (err != CUFFT_SUCCESS)
  25 + {
  26 + if(err == CUFFT_INVALID_PLAN)
  27 + printf("The plan parameter is not a valid handle.\n");
  28 + if(err == CUFFT_INVALID_VALUE)
  29 + printf("At least one of the parameters idata, odata, and direction is not valid.\n");
  30 + if(err == CUFFT_INTERNAL_ERROR)
  31 + printf("An internal driver error was detected.\n");
  32 + if(err == CUFFT_EXEC_FAILED)
  33 + printf("CUFFT failed to execute the transform on the GPU.\n");
  34 + if(err == CUFFT_SETUP_FAILED)
  35 + printf("The CUFFT library failed to initialize.\n");
  36 +
  37 + }
  38 +}
20 39  
21 40  
22 41  
... ...
rts/envi/envi.h
  Show/Hide comments   View file @6e257ab
... ... @@ -4,401 +4,17 @@
4 4 #include <string>
5 5 #include <iostream>
6 6 #include <fstream>
7   -#include <sstream>
  7 +#include <sstream>
  8 +
  9 +#include "rts/envi/envi_header.h"
8 10  
9 11 /* This is a class for reading and writing ENVI binary files. This class will be updated as needed.
10 12  
11 13 What this class CAN currently do:
12 14 *) Write band images to a BSQ file.
13   - *) Append band images to a BSQ file.
14   -
  15 + *) Append band images to a BSQ file.
  16 + *) Interpolate and query bands in a float32 BSQ file
15 17 */
16   -
17   -//information from an ENVI header file
18   -//A good resource can be found here: http://www.exelisvis.com/docs/enviheaderfiles.html
19   -struct EnviHeader
20   -{
21   - std::string name;
22   -
23   - std::string description;
24   -
25   - unsigned int samples; //x-axis
26   - unsigned int lines; //y-axis
27   - unsigned int bands; //spectral axis
28   - unsigned int header_offset; //header offset for binary file (in bytes)
29   - std::string file_type; //should be "ENVI Standard"
30   - unsigned int data_type; //data representation; common value is 4 (32-bit float)
31   -
32   - enum interleaveType {BIP, BIL, BSQ}; //bip = Z,X,Y; bil = X,Z,Y; bsq = X,Y,Z
33   - interleaveType interleave;
34   -
35   - std::string sensor_type; //not really used
36   -
37   - enum endianType {endianLittle, endianBig};
38   - endianType byte_order; //little = least significant bit first (most systems)
39   -
40   - double x_start, y_start; //coordinates of the upper-left corner of the image
41   - std::string wavelength_units; //stored wavelength units
42   - std::string z_plot_titles[2];
43   -
44   - double pixel_size[2]; //pixel size along X and Y
45   -
46   - std::vector<std::string> band_names; //name for each band in the image
47   - std::vector<double> wavelength; //wavelength for each band
48   -
49   - EnviHeader()
50   - {
51   - name = "";
52   -
53   - //specify default values for a new or empty ENVI file
54   - samples = 0;
55   - lines = 0;
56   - bands = 0;
57   - header_offset = 0;
58   - data_type = 4;
59   - interleave = BSQ;
60   - byte_order = endianLittle;
61   - x_start = y_start = 0;
62   - pixel_size[0] = pixel_size[1] = 1;
63   -
64   - //strings
65   - file_type = "ENVI Standard";
66   - sensor_type = "Unknown";
67   - wavelength_units = "Unknown";
68   - z_plot_titles[0] = z_plot_titles[1] = "Unknown";
69   - }
70   -
71   - std::string trim(std::string line)
72   - {
73   - //trims whitespace from the beginning and end of line
74   - int start_i, end_i;
75   - for(start_i=0; start_i < line.length(); start_i++)
76   - if(line[start_i] != 32)
77   - {
78   - break;
79   - }
80   -
81   - for(end_i = line.length()-1; end_i >= start_i; end_i--)
82   - if(line[end_i] != ' ')
83   - {
84   - break;
85   - }
86   -
87   - return line.substr(start_i, end_i - start_i+1);
88   - }
89   -
90   -
91   - std::string get_token(std::string line)
92   - {
93   - //returns a variable name; in this case we look for the '=' sign
94   - size_t i = line.find_first_of('=');
95   -
96   - std::string result;
97   - if(i != std::string::npos)
98   - result = trim(line.substr(0, i-1));
99   -
100   - return result;
101   - }
102   -
103   - std::string get_data_str(std::string line)
104   - {
105   - size_t i = line.find_first_of('=');
106   -
107   - std::string result;
108   - if(i != std::string::npos)
109   - result = trim(line.substr(i+1));
110   - else
111   - {
112   - std::cout<<"ENVI Header error - data not found for token: "<<get_token(line)<<std::endl;
113   - exit(1);
114   - }
115   - return result;
116   - }
117   -
118   - std::string get_brace_str(std::string token, std::string line, std::ifstream &file)
119   - {
120   - //this function assembles all of the characters between curly braces
121   - //this is how strings are defined within an ENVI file
122   -
123   - std::string result;
124   -
125   - //first, find the starting brace
126   - size_t i;
127   - do
128   - {
129   - i = line.find_first_of('{');
130   - if(i != std::string::npos)
131   - break;
132   - }while(file);
133   -
134   - //if i is still npos, we have reached the end of the file without a brace...something is wrong
135   - if(i == std::string::npos)
136   - {
137   - std::cout<<"ENVI Header error - string token declared without being defined: "<<token<<std::endl;
138   - exit(1);
139   - }
140   - line = line.substr(i+1);
141   -
142   - //copy character data into the result string until we find a closing brace
143   - while(file)
144   - {
145   - i = line.find_first_of('}');
146   -
147   -
148   - if(i != std::string::npos)
149   - {
150   - result += line.substr(0, i);
151   - break;
152   - }
153   - else
154   - result += line;
155   -
156   - getline(file, line);
157   - }
158   -
159   - if(i == std::string::npos)
160   - {
161   - std::cout<<"ENVI Header error - string token declared without a terminating '}': "<<token<<std::endl;
162   - exit(1);
163   - }
164   -
165   - return trim(result);
166   - }
167   -
168   - std::vector<std::string> get_string_seq(std::string token, std::string sequence)
169   - {
170   - //this function returns a sequence of comma-delimited strings
171   - std::vector<std::string> result;
172   -
173   - std::string entry;
174   - size_t i;
175   - do
176   - {
177   - i = sequence.find_first_of(',');
178   - entry = sequence.substr(0, i);
179   - sequence = sequence.substr(i+1);
180   - result.push_back(trim(entry));
181   - }while(i != std::string::npos);
182   -
183   - return result;
184   - }
185   -
186   - std::vector<double> get_double_seq(std::string token, std::string sequence)
187   - {
188   - //this function returns a sequence of comma-delimited strings
189   - std::vector<double> result;
190   -
191   - std::string entry;
192   - size_t i;
193   - do
194   - {
195   - i = sequence.find_first_of(',');
196   - entry = sequence.substr(0, i);
197   - sequence = sequence.substr(i+1);
198   - result.push_back(atof(entry.c_str()));
199   - }while(i != std::string::npos);
200   -
201   - return result;
202   - }
203   -
204   - void load(std::string filename, std::string mode = "r")
205   - {
206   - name = filename;
207   -
208   -
209   - //open the header file
210   - std::ifstream file(name.c_str());
211   -
212   - if(!file)
213   - {
214   - //if we are in read mode, throw an exception
215   - if(mode == "r")
216   - {
217   - std::cout<<"ENVI header file not found: "<<name<<std::endl;
218   - exit(1);
219   - }
220   - else
221   - {
222   - //there is no file, so just return the default structure
223   - return;
224   - }
225   - }
226   -
227   - //the first line should just be "ENVI"
228   - std::string line;
229   - getline(file, line);
230   - if(line != "ENVI")
231   - {
232   - std::cout<<"The header doesn't appear to be an ENVI file. The first line should be 'ENVI'."<<std::endl;
233   - exit(1);
234   - }
235   -
236   - //for each line in the file, get the token
237   - std::string token;
238   - while(file)
239   - {
240   -
241   - //get a line
242   - getline(file, line);
243   -
244   - //get the token
245   - token = get_token(line);
246   -
247   - if(token == "description")
248   - description = get_brace_str(token, line, file);
249   - else if(token == "band names")
250   - {
251   - std::string string_sequence = get_brace_str(token, line, file);
252   - band_names = get_string_seq(token, string_sequence);
253   - }
254   - else if(token == "wavelength")
255   - {
256   - std::string string_sequence = get_brace_str(token, line, file);
257   - wavelength = get_double_seq(token, string_sequence);
258   - }
259   - else if(token == "pixel size")
260   - {
261   - std::string string_sequence = get_brace_str(token, line, file);
262   - std::vector<double> pxsize = get_double_seq(token, string_sequence);
263   - pixel_size[0] = pxsize[0];
264   - pixel_size[1] = pxsize[1];
265   - }
266   - else if(token == "z plot titles")
267   - {
268   - std::string string_sequence = get_brace_str(token, line, file);
269   - std::vector<std::string> titles = get_string_seq(token, string_sequence);
270   - z_plot_titles[0] = titles[0];
271   - z_plot_titles[1] = titles[1];
272   - }
273   -
274   - else if(token == "samples")
275   - samples = atoi(get_data_str(line).c_str());
276   - else if(token == "lines")
277   - {
278   - lines = atoi(get_data_str(line).c_str());
279   - std::cout<<"Number of lines loaded: "<<lines<<std::endl;
280   - }
281   - else if(token == "bands")
282   - bands = atoi(get_data_str(line).c_str());
283   - else if(token == "header offset")
284   - header_offset = atoi(get_data_str(line).c_str());
285   - else if(token == "file type")
286   - file_type = get_data_str(line);
287   - else if(token == "data type")
288   - data_type = atoi(get_data_str(line).c_str());
289   - else if(token == "interleave")
290   - {
291   - std::string interleave_str = get_data_str(line);
292   - if(interleave_str == "bip")
293   - interleave = BIP;
294   - else if(interleave_str == "bil")
295   - interleave = BIL;
296   - else if(interleave_str == "bsq")
297   - interleave = BSQ;
298   - }
299   - else if(token == "sensor type")
300   - sensor_type = get_data_str(line);
301   - else if(token == "byte order")
302   - byte_order = (endianType)atoi(get_data_str(line).c_str());
303   - else if(token == "x start")
304   - x_start = atof(get_data_str(line).c_str());
305   - else if(token == "y start")
306   - y_start = atof(get_data_str(line).c_str());
307   - else if(token == "wavelength units")
308   - wavelength_units = get_data_str(line);
309   - }
310   -
311   - //close the file
312   - file.close();
313   - }
314   -
315   - void save(std::string filename)
316   - {
317   - //open a file
318   - std::ofstream outfile(filename.c_str());
319   -
320   - //write the ENVI type identifier
321   - outfile<<"ENVI"<<std::endl;
322   -
323   - //output all of the data
324   - outfile<<"description = {"<<std::endl;
325   - outfile<<" "<<description<<"}"<<std::endl;
326   -
327   - outfile<<"samples = "<<samples<<std::endl;
328   - outfile<<"lines = "<<lines<<std::endl;
329   - outfile<<"bands = "<<bands<<std::endl;
330   - outfile<<"header offset = "<<header_offset<<std::endl;
331   - outfile<<"file type = "<<file_type<<std::endl;
332   - outfile<<"data type = "<<data_type<<std::endl;
333   - outfile<<"interleave = ";
334   - if(interleave == BIP)
335   - outfile<<"bip";
336   - if(interleave == BIL)
337   - outfile<<"bil";
338   - if(interleave == BSQ)
339   - outfile<<"bsq";
340   - outfile<<std::endl;
341   - outfile<<"sensor type = "<<sensor_type<<std::endl;
342   - outfile<<"byte order = "<<byte_order<<std::endl;
343   - outfile<<"x start = "<<x_start<<std::endl;
344   - outfile<<"y start = "<<y_start<<std::endl;
345   - outfile<<"wavelength units = "<<wavelength_units<<std::endl;
346   - outfile<<"z plot titles = {";
347   - outfile<<z_plot_titles[0]<<", "<<z_plot_titles[1]<<"}"<<std::endl;
348   - outfile<<"pixel size = {"<<pixel_size[0]<<", "<<pixel_size[1]<<", units=Meters}"<<std::endl;
349   - if(band_names.size() > 0)
350   - {
351   - outfile<<"band names = {"<<std::endl;
352   - for(int i=0; i<band_names.size(); i++)
353   - {
354   - outfile<<band_names[i];
355   - if(i < band_names.size() - 1)
356   - outfile<<", ";
357   - }
358   - outfile<<"}"<<std::endl;
359   - }
360   - outfile<<"wavelength = {"<<std::endl;
361   - for(int i=0; i<wavelength.size()-1; i++)
362   - outfile<<wavelength[i]<<", ";
363   - outfile<<wavelength.back()<<"}"<<std::endl;
364   -
365   - outfile.close();
366   - }
367   -
368   - void save()
369   - {
370   - //std::cout<<"ENVI Header Name: "<<name<<std::endl;
371   - save(name);
372   - }
373   -
374   - //returns the size of a single value (in bytes)
375   - unsigned int valsize()
376   - {
377   - switch(data_type)
378   - {
379   - case 1: //1 = 8-bit byte
380   - return 1;
381   - case 2: //16-bit signed integer
382   - case 12: //16-bit unsigned integer
383   - return 2;
384   - case 3: //32-bit signed long integer
385   - case 4: //32-bit floating point
386   - case 13: //32-bit unsigned long integer
387   - return 4;
388   - case 5: //64-bit double-precision floating point
389   - case 6: //32-bit complex value
390   - case 14: //64-bit signed long integer
391   - case 15: //64-bit unsigned long integer
392   - return 8;
393   - case 9: //64-bit complex value
394   - return 16;
395   - default:
396   - return 0;
397   - }
398   -
399   - }
400   -}; //end EnviHeader
401   -
402 18 class EnviFile
403 19 {
404 20 EnviHeader header;
... ... @@ -430,7 +46,7 @@ class EnviFile
430 46 readonly(); //if the file is read-only, throw an exception
431 47  
432 48 //we currently only support BSQ
433   - if(new_header.interleave == EnviHeader::BIP || new_header.interleave == EnviHeader::BIL)
  49 + if(new_header.interleave == Envi::BIP || new_header.interleave == Envi::BIL)
434 50 {
435 51 std::cout<<"This software only supports BSQ.";
436 52 return false;
... ... @@ -440,7 +56,7 @@ class EnviFile
440 56 if(header.bands != new_header.bands)
441 57 {
442 58 //the new header must be a BSQ file
443   - if(new_header.interleave != EnviHeader::BSQ)
  59 + if(new_header.interleave != Envi::BSQ)
444 60 {
445 61 std::cout<<"ENVI Error: can only add bands to a BSQ file."<<std::endl;
446 62 return false;
... ... @@ -476,10 +92,10 @@ class EnviFile
476 92  
477 93 }
478 94  
479   - EnviHeader load_header(std::string header_file, std::string file_mode)
  95 + EnviHeader load_header(std::string header_file)
480 96 {
481 97 EnviHeader new_header;
482   - new_header.load(header_file, file_mode);
  98 + new_header.load(header_file);
483 99 return new_header;
484 100 }
485 101  
... ... @@ -512,12 +128,64 @@ class EnviFile
512 128 else
513 129 return false;
514 130  
  131 + }
  132 +
  133 + void get_surrounding_bands(double wavelength, unsigned int &low, unsigned int &high)
  134 + {
  135 + int i;
  136 + int nBands = header.wavelength.size();
  137 + low = high = 0;
  138 + for(i=0; i<nBands; i++)
  139 + {
  140 + if(header.wavelength[i] < wavelength)
  141 + {
  142 + if(header.wavelength[low] > wavelength)
  143 + low = i;
  144 + if(header.wavelength[i] > header.wavelength[low])
  145 + low = i;
  146 + }
  147 +
  148 + if(header.wavelength[i] > wavelength)
  149 + {
  150 + if(header.wavelength[high] < wavelength)
  151 + high = i;
  152 + if(header.wavelength[i] < header.wavelength[high])
  153 + high = i;
  154 + }
  155 + }
  156 + }
  157 +
  158 + void interpolate_band(void* result, void* A, void* B, double a)
  159 + {
  160 + //interpolate between two bands independent of data type
  161 + // CURRENTLY ONLY HANDLES FLOAT
  162 + int X = header.samples;
  163 + int Y = header.lines;
  164 + int N = X * Y;
  165 +
  166 + if(header.data_type != Envi::float32)
  167 + {
  168 + std::cout<<"ENVI Error: this software only handles interpolation of 32-bit floats."<<std::endl;
  169 + exit(1);
  170 + }
  171 +
  172 + float r, v0, v1;
  173 + for(int n=0; n<N; n++)
  174 + {
  175 + v0 = ((float*)A)[n];
  176 + v1 = ((float*)B)[n];
  177 + float r = v0 * a + v1 * (1.0 - a);
  178 + //((float*)result)[n] = ((float*)A)[n] * a + ((float*)B)[n] * (1.0 - a);
  179 + ((float*)result)[n] = r;
  180 + }
  181 +
515 182 }
516 183  
517 184 public:
518 185  
519 186 EnviFile()
520   - {
  187 + {
  188 + mode = "w";
521 189 init();
522 190 }
523 191  
... ... @@ -528,30 +196,41 @@ public:
528 196 }
529 197  
530 198 void open(std::string file_name, std::string header_name, std::string file_mode = "r")
531   - {
  199 + {
  200 + //load the header file
  201 + // if the file doesn't exist, that is okay unless it is read-only
  202 + if(!header.load(header_name))
  203 + {
  204 + if(file_mode == "r")
  205 + {
  206 + std::cout<<"ENVI Error: header file not found: "<<header_name<<std::endl;
  207 + exit(1);
  208 + }
  209 + //otherwise use the default header and save the name
  210 + header.name = header_name;
  211 + }
  212 +
532 213 mode = file_mode;
533 214  
534   - header.name = header_name;
535   -
536   - //lock the file if we are loading or appending
537   - if(mode == "r" || mode == "a" || mode == "r+" || mode == "a+")
  215 + //lock the file if it is read-only
  216 + if(mode == "r")
  217 + locked = true;
  218 + else if(mode == "a" || mode == "r+" || mode == "a+")
538 219 {
539   - //load the ENVI header - failure will cause an exit
540   - set_header(load_header(header_name, file_mode));
541   -
542   - if(mode == "r" || test_exist(file_name))
543   - {
544   - locked = true;
545   - }
546   -
  220 + //if the file can be edited, lock it if data has already been written
  221 + if(test_exist(file_name))
  222 + locked = true;
547 223 }
548   - else
549   - {
550   - header.name = header_name;
551   - }
  224 +
552 225  
553 226 //open the data file
554   - data = fopen(file_name.c_str(), (mode + "b").c_str());
  227 + std::string tmpmode = mode+"b";
  228 + data = fopen(file_name.c_str(), tmpmode.c_str());
  229 + if(data == NULL)
  230 + {
  231 + std::cout<<"ENVI Error: unable to open binary file: "<<file_name<<std::endl;
  232 + exit(1);
  233 + }
555 234  
556 235  
557 236 }
... ... @@ -602,14 +281,35 @@ public:
602 281 EnviHeader newHeader = header;
603 282 newHeader.lines = lines;
604 283 set_header(newHeader);
605   - }
  284 + }
  285 +
  286 + unsigned int lines() {return header.lines;}
  287 + unsigned int samples() {return header.samples;}
  288 + unsigned int bands() {return header.bands;}
  289 + double wavelength(unsigned int b) {return header.wavelength[b];}
606 290  
607 291 void setBands(unsigned int bands)
608 292 {
609 293 EnviHeader newHeader = header;
610 294 newHeader.bands = bands;
611 295 set_header(newHeader);
612   - }
  296 + }
  297 +
  298 + unsigned int getElementSize()
  299 + {
  300 + //return the size of each element in bytes
  301 + return header.valsize();
  302 + }
  303 + unsigned int getBandSize()
  304 + {
  305 + //returns the size of a band in bytes
  306 + return header.valsize() * header.lines * header.samples;
  307 + }
  308 +
  309 + Envi::dataType getDataType()
  310 + {
  311 + return header.data_type;
  312 + }
613 313  
614 314 //DATA MANIPULATION
615 315 void addBand(void* band, unsigned int samples, unsigned int lines, double wavelength, std::string band_name ="")
... ... @@ -630,7 +330,83 @@ public:
630 330  
631 331 //since data was written, the file must be locked
632 332 locked = true;
633   - }
  333 + }
  334 +
  335 + //DATA RETRIEVAL
  336 + void getBand(void* ptr, unsigned int n)
  337 + {
  338 + //reads the n'th band and writes it to the given pointer
  339 + // memory has been assumed to be allocated
  340 +
  341 + if(n > header.wavelength.size())
  342 + {
  343 + std::cout<<"ENVI Error: Invalid band number."<<std::endl;
  344 + return;
  345 + }
  346 +
  347 +
  348 + if(header.interleave == Envi::BSQ)
  349 + {
  350 +
  351 + //compute the position of the desired band (relative to the beginning of the file)
  352 + int bandpos = header.header_offset + getBandSize() * n;
  353 +
  354 + //seek to that position
  355 + int seek_result = fseek(data, bandpos, SEEK_SET);
  356 + if(seek_result)
  357 + {
  358 + std::cout<<"ENVI Error: Error seeking through data file."<<std::endl;
  359 + return;
  360 + }
  361 +
  362 + //read the band data from the file
  363 + size_t n_read = fread(ptr, getElementSize(), header.samples * header.lines, data);
  364 + if(n_read != header.samples * header.lines)
  365 + {
  366 + std::cout<<"ENVI Error: Error reading band data."<<std::endl;
  367 + return;
  368 + }
  369 +
  370 + return;
  371 + }
  372 + else
  373 + {
  374 + std::cout<<"ENVI Error: Only BSQ operations are currently supported."<<std::endl;
  375 + return;
  376 + }
  377 + }
  378 +
  379 + void getWavelength(void* band, double wavelength)
  380 + {
  381 + //this function retrieves a band via interpolation of wavelength values
  382 +
  383 + //get the low and high band numbers
  384 + unsigned int low, high;
  385 + get_surrounding_bands(wavelength, low, high);
  386 +
  387 + //calculate the position between the two bands
  388 + double a;
  389 + if(low == high)
  390 + a = 1.0;
  391 + else
  392 + a = (wavelength - header.wavelength[low]) / (header.wavelength[high] - header.wavelength[low]);
  393 +
  394 + //read both bands
  395 + void* A;
  396 + A = malloc( getBandSize() );
  397 + getBand(A, low);
  398 +
  399 + void* B;
  400 + B = malloc( getBandSize() );
  401 + getBand(B, high);
  402 +
  403 + //interpolate between the bands
  404 + interpolate_band(band, A, B, a);
  405 +
  406 + //free the surrounding bands
  407 + free(A);
  408 + free(B);
  409 + }
634 410  
635 411 //close file
636 412 void close()
... ...
rts/envi/envi_header.h 0 → 100644
  Show/Hide comments   View file @6e257ab
  1 +#ifndef ENVI_HEADER_H
  2 +#define ENVI_HEADER_H
  3 +
  4 +#include <string>
  5 +#include <iostream>
  6 +#include <fstream>
  7 +#include <sstream>
  8 +#include <vector>
  9 +#include <stdlib.h>
  10 +
  11 +//information from an ENVI header file
  12 +//A good resource can be found here: http://www.exelisvis.com/docs/enviheaderfiles.html
  13 +
  14 +namespace Envi
  15 +{
  16 + enum dataType {dummy0, int8, int16, int32, float32, float64, complex32, dummy7, dummy8, complex64, dummy10, dummy11, uint16, uint32, int64, uint64};
  17 + enum interleaveType {BIP, BIL, BSQ}; //bip = Z,X,Y; bil = X,Z,Y; bsq = X,Y,Z
  18 + enum endianType {endianLittle, endianBig};
  19 +}
  20 +struct EnviHeader
  21 +{
  22 + std::string name;
  23 +
  24 + std::string description;
  25 +
  26 + unsigned int samples; //x-axis
  27 + unsigned int lines; //y-axis
  28 + unsigned int bands; //spectral axis
  29 + unsigned int header_offset; //header offset for binary file (in bytes)
  30 + std::string file_type; //should be "ENVI Standard"
  31 +
  32 + Envi::dataType data_type; //data representation; common value is 4 (32-bit float)
  33 +
  34 +
  35 + Envi::interleaveType interleave;
  36 +
  37 + std::string sensor_type; //not really used
  38 +
  39 + Envi::endianType byte_order; //little = least significant bit first (most systems)
  40 +
  41 + double x_start, y_start; //coordinates of the upper-left corner of the image
  42 + std::string wavelength_units; //stored wavelength units
  43 + std::string z_plot_titles[2];
  44 +
  45 + double pixel_size[2]; //pixel size along X and Y
  46 +
  47 + std::vector<std::string> band_names; //name for each band in the image
  48 + std::vector<double> wavelength; //wavelength for each band
  49 +
  50 + EnviHeader()
  51 + {
  52 + name = "";
  53 +
  54 + //specify default values for a new or empty ENVI file
  55 + samples = 0;
  56 + lines = 0;
  57 + bands = 0;
  58 + header_offset = 0;
  59 + data_type = Envi::float32;
  60 + interleave = Envi::BSQ;
  61 + byte_order = Envi::endianLittle;
  62 + x_start = y_start = 0;
  63 + pixel_size[0] = pixel_size[1] = 1;
  64 +
  65 + //strings
  66 + file_type = "ENVI Standard";
  67 + sensor_type = "Unknown";
  68 + wavelength_units = "Unknown";
  69 + z_plot_titles[0] = z_plot_titles[1] = "Unknown";
  70 + }
  71 +
  72 + std::string trim(std::string line)
  73 + {
  74 + //trims whitespace from the beginning and end of line
  75 + int start_i, end_i;
  76 + for(start_i=0; start_i < line.length(); start_i++)
  77 + if(line[start_i] != 32)
  78 + {
  79 + break;
  80 + }
  81 +
  82 + for(end_i = line.length()-1; end_i >= start_i; end_i--)
  83 + if(line[end_i] != ' ')
  84 + {
  85 + break;
  86 + }
  87 +
  88 + return line.substr(start_i, end_i - start_i+1);
  89 + }
  90 +
  91 +
  92 + std::string get_token(std::string line)
  93 + {
  94 + //returns a variable name; in this case we look for the '=' sign
  95 + size_t i = line.find_first_of('=');
  96 +
  97 + std::string result;
  98 + if(i != std::string::npos)
  99 + result = trim(line.substr(0, i-1));
  100 +
  101 + return result;
  102 + }
  103 +
  104 + std::string get_data_str(std::string line)
  105 + {
  106 + size_t i = line.find_first_of('=');
  107 +
  108 + std::string result;
  109 + if(i != std::string::npos)
  110 + result = trim(line.substr(i+1));
  111 + else
  112 + {
  113 + std::cout<<"ENVI Header error - data not found for token: "<<get_token(line)<<std::endl;
  114 + exit(1);
  115 + }
  116 + return result;
  117 + }
  118 +
  119 + std::string get_brace_str(std::string token, std::string line, std::ifstream &file)
  120 + {
  121 + //this function assembles all of the characters between curly braces
  122 + //this is how strings are defined within an ENVI file
  123 +
  124 + std::string result;
  125 +
  126 + //first, find the starting brace
  127 + size_t i;
  128 + do
  129 + {
  130 + i = line.find_first_of('{');
  131 + if(i != std::string::npos)
  132 + break;
  133 + }while(file);
  134 +
  135 + //if i is still npos, we have reached the end of the file without a brace...something is wrong
  136 + if(i == std::string::npos)
  137 + {
  138 + std::cout<<"ENVI Header error - string token declared without being defined: "<<token<<std::endl;
  139 + exit(1);
  140 + }
  141 + line = line.substr(i+1);
  142 +
  143 + //copy character data into the result string until we find a closing brace
  144 + while(file)
  145 + {
  146 + i = line.find_first_of('}');
  147 +
  148 +
  149 + if(i != std::string::npos)
  150 + {
  151 + result += line.substr(0, i);
  152 + break;
  153 + }
  154 + else
  155 + result += line;
  156 +
  157 + getline(file, line);
  158 + }
  159 +
  160 + if(i == std::string::npos)
  161 + {
  162 + std::cout<<"ENVI Header error - string token declared without a terminating '}': "<<token<<std::endl;
  163 + exit(1);
  164 + }
  165 +
  166 + return trim(result);
  167 + }
  168 +
  169 + std::vector<std::string> get_string_seq(std::string token, std::string sequence)
  170 + {
  171 + //this function returns a sequence of comma-delimited strings
  172 + std::vector<std::string> result;
  173 +
  174 + std::string entry;
  175 + size_t i;
  176 + do
  177 + {
  178 + i = sequence.find_first_of(',');
  179 + entry = sequence.substr(0, i);
  180 + sequence = sequence.substr(i+1);
  181 + result.push_back(trim(entry));
  182 + }while(i != std::string::npos);
  183 +
  184 + return result;
  185 + }
  186 +
  187 + std::vector<double> get_double_seq(std::string token, std::string sequence)
  188 + {
  189 + //this function returns a sequence of comma-delimited strings
  190 + std::vector<double> result;
  191 +
  192 + std::string entry;
  193 + size_t i;
  194 + do
  195 + {
  196 + i = sequence.find_first_of(',');
  197 + entry = sequence.substr(0, i);
  198 + sequence = sequence.substr(i+1);
  199 + result.push_back(atof(entry.c_str()));
  200 + }while(i != std::string::npos);
  201 +
  202 + return result;
  203 + }
  204 +
  205 + bool load(std::string filename)
  206 + {
  207 + //open the header file
  208 + std::ifstream file(filename.c_str());
  209 +
  210 + if(!file)
  211 + {
  212 + return false;
  213 + }
  214 +
  215 + //the first line should just be "ENVI"
  216 + std::string line;
  217 + file>>line;
  218 + if(line != "ENVI")
  219 + {
  220 + std::cout<<"ENVI Header Error: The header doesn't appear to be an ENVI file. The first line should be 'ENVI'."<<std::endl;
  221 + exit(1);
  222 + }
  223 +
  224 + //for each line in the file, get the token
  225 + std::string token;
  226 + while(file)
  227 + {
  228 +
  229 + //get a line
  230 + getline(file, line);
  231 +
  232 + //get the token
  233 + token = get_token(line);
  234 +
  235 + if(token == "description")
  236 + description = get_brace_str(token, line, file);
  237 + else if(token == "band names")
  238 + {
  239 + std::string string_sequence = get_brace_str(token, line, file);
  240 + band_names = get_string_seq(token, string_sequence);
  241 + }
  242 + else if(token == "wavelength")
  243 + {
  244 + std::string string_sequence = get_brace_str(token, line, file);
  245 + wavelength = get_double_seq(token, string_sequence);
  246 + }
  247 + else if(token == "pixel size")
  248 + {
  249 + std::string string_sequence = get_brace_str(token, line, file);
  250 + std::vector<double> pxsize = get_double_seq(token, string_sequence);
  251 + pixel_size[0] = pxsize[0];
  252 + pixel_size[1] = pxsize[1];
  253 + }
  254 + else if(token == "z plot titles")
  255 + {
  256 + std::string string_sequence = get_brace_str(token, line, file);
  257 + std::vector<std::string> titles = get_string_seq(token, string_sequence);
  258 + z_plot_titles[0] = titles[0];
  259 + z_plot_titles[1] = titles[1];
  260 + }
  261 +
  262 + else if(token == "samples")
  263 + samples = atoi(get_data_str(line).c_str());
  264 + else if(token == "lines")
  265 + lines = atoi(get_data_str(line).c_str());
  266 + else if(token == "bands")
  267 + bands = atoi(get_data_str(line).c_str());
  268 + else if(token == "header offset")
  269 + header_offset = atoi(get_data_str(line).c_str());
  270 + else if(token == "file type")
  271 + file_type = get_data_str(line);
  272 + else if(token == "data type")
  273 + data_type = (Envi::dataType)atoi(get_data_str(line).c_str());
  274 + else if(token == "interleave")
  275 + {
  276 + std::string interleave_str = get_data_str(line);
  277 + if(interleave_str == "bip")
  278 + interleave = Envi::BIP;
  279 + else if(interleave_str == "bil")
  280 + interleave = Envi::BIL;
  281 + else if(interleave_str == "bsq")
  282 + interleave = Envi::BSQ;
  283 + }
  284 + else if(token == "sensor type")
  285 + sensor_type = get_data_str(line);
  286 + else if(token == "byte order")
  287 + byte_order = (Envi::endianType)atoi(get_data_str(line).c_str());
  288 + else if(token == "x start")
  289 + x_start = atof(get_data_str(line).c_str());
  290 + else if(token == "y start")
  291 + y_start = atof(get_data_str(line).c_str());
  292 + else if(token == "wavelength units")
  293 + wavelength_units = get_data_str(line);
  294 + }
  295 +
  296 + //close the file
  297 + file.close();
  298 +
  299 + //set the file name
  300 + name = filename;
  301 +
  302 + return true;
  303 + }
  304 +
  305 + void save(std::string filename)
  306 + {
  307 + //open a file
  308 + std::ofstream outfile(filename.c_str());
  309 +
  310 + //write the ENVI type identifier
  311 + outfile<<"ENVI"<<std::endl;
  312 +
  313 + //output all of the data
  314 + outfile<<"description = {"<<std::endl;
  315 + outfile<<" "<<description<<"}"<<std::endl;
  316 +
  317 + outfile<<"samples = "<<samples<<std::endl;
  318 + outfile<<"lines = "<<lines<<std::endl;
  319 + outfile<<"bands = "<<bands<<std::endl;
  320 + outfile<<"header offset = "<<header_offset<<std::endl;
  321 + outfile<<"file type = "<<file_type<<std::endl;
  322 + outfile<<"data type = "<<data_type<<std::endl;
  323 + outfile<<"interleave = ";
  324 + if(interleave == Envi::BIP)
  325 + outfile<<"bip";
  326 + if(interleave == Envi::BIL)
  327 + outfile<<"bil";
  328 + if(interleave == Envi::BSQ)
  329 + outfile<<"bsq";
  330 + outfile<<std::endl;
  331 + outfile<<"sensor type = "<<sensor_type<<std::endl;
  332 + outfile<<"byte order = "<<byte_order<<std::endl;
  333 + outfile<<"x start = "<<x_start<<std::endl;
  334 + outfile<<"y start = "<<y_start<<std::endl;
  335 + outfile<<"wavelength units = "<<wavelength_units<<std::endl;
  336 + outfile<<"z plot titles = {";
  337 + outfile<<z_plot_titles[0]<<", "<<z_plot_titles[1]<<"}"<<std::endl;
  338 + outfile<<"pixel size = {"<<pixel_size[0]<<", "<<pixel_size[1]<<", units=Meters}"<<std::endl;
  339 + if(band_names.size() > 0)
  340 + {
  341 + outfile<<"band names = {"<<std::endl;
  342 + for(int i=0; i<band_names.size(); i++)
  343 + {
  344 + outfile<<band_names[i];
  345 + if(i < band_names.size() - 1)
  346 + outfile<<", ";
  347 + }
  348 + outfile<<"}"<<std::endl;
  349 + }
  350 + outfile<<"wavelength = {"<<std::endl;
  351 + for(int i=0; i<wavelength.size()-1; i++)
  352 + outfile<<wavelength[i]<<", ";
  353 + outfile<<wavelength.back()<<"}"<<std::endl;
  354 +
  355 + outfile.close();
  356 + }
  357 +
  358 + void save()
  359 + {
  360 + //std::cout<<"ENVI Header Name: "<<name<<std::endl;
  361 + save(name);
  362 + }
  363 +
  364 + //returns the size of a single value (in bytes)
  365 + unsigned int valsize()
  366 + {
  367 + switch(data_type)
  368 + {
  369 + case Envi::int8: //1 = 8-bit byte
  370 + return 1;
  371 + case Envi::int16: //16-bit signed integer
  372 + case Envi::uint16: //16-bit unsigned integer
  373 + return 2;
  374 + case Envi::int32: //32-bit signed long integer
  375 + case Envi::float32: //32-bit floating point
  376 + case Envi::uint32: //32-bit unsigned long integer
  377 + return 4;
  378 + case Envi::float64: //64-bit double-precision floating point
  379 + case Envi::complex32: //32-bit complex value
  380 + case Envi::int64: //64-bit signed long integer
  381 + case Envi::uint64: //64-bit unsigned long integer
  382 + return 8;
  383 + case Envi::complex64: //64-bit complex value
  384 + return 16;
  385 + default:
  386 + return 0;
  387 + }
  388 +
  389 + }
  390 +}; //end EnviHeader
  391 +
  392 +#endif
... ...
rts/optics/material.h
  Show/Hide comments   View file @6e257ab
... ... @@ -558,6 +558,17 @@ namespace rts{
558 558 }
559 559  
560 560 }
  561 +
  562 + //returns the scattering efficiency at wavelength l
  563 + rtsComplex<T> eta(T l)
  564 + {
  565 + //get the refractive index
  566 + rtsComplex<T> ri = getN(l);
  567 +
  568 + //convert the refractive index to scattering efficiency
  569 + return ri*ri - 1.0;
  570 +
  571 + }
561 572 //interpolate the given lambda value and return the index of refraction
562 573 rtsComplex<T> operator()(T l)
563 574 {
... ...
rts/tools/colormap.h 0 → 100644
  Show/Hide comments   View file @6e257ab
  1 +#ifndef RTS_COLORMAP_H
  2 +#define RTS_COLORMAP_H
  3 +
  4 +#include <string>
  5 +#include <qimage.h>
  6 +#include <qcolor.h>
  7 +#include "rts/cuda/error.h"
  8 +
  9 +
  10 +#define BREWER_CTRL_PTS 11
  11 +
  12 +static float BREWERCP[BREWER_CTRL_PTS*4] = {0.192157f, 0.211765f, 0.584314f, 1.0f,
  13 + 0.270588f, 0.458824f, 0.705882f, 1.0f,
  14 + 0.454902f, 0.678431f, 0.819608f, 1.0f,
  15 + 0.670588f, 0.85098f, 0.913725f, 1.0f,
  16 + 0.878431f, 0.952941f, 0.972549f, 1.0f,
  17 + 1.0f, 1.0f, 0.74902f, 1.0f,
  18 + 0.996078f, 0.878431f, 0.564706f, 1.0f,
  19 + 0.992157f, 0.682353f, 0.380392f, 1.0f,
  20 + 0.956863f, 0.427451f, 0.262745f, 1.0f,
  21 + 0.843137f, 0.188235f, 0.152941f, 1.0f,
  22 + 0.647059f, 0.0f, 0.14902f, 1.0f};
  23 +
  24 +
  25 +#ifdef __CUDACC__
  26 +texture<float4, cudaTextureType1D> cudaTexBrewer;
  27 +static cudaArray* gpuBrewer;
  28 +#endif
  29 +
  30 +
  31 +
  32 +namespace rts{
  33 +
  34 +enum colormapType {cmBrewer, cmGrayscale};
  35 +
  36 +static void buffer2image(unsigned char* buffer, std::string filename, unsigned int x_size, unsigned int y_size)
  37 +{
  38 + //x_size = 256;
  39 + //y_size = 256;
  40 + //create an image object
  41 + QImage image(x_size, y_size, QImage::Format_RGB32);
  42 + if(image.isNull())
  43 + {
  44 + std::cout<<"Error creating QImage."<<std::endl;
  45 + return;
  46 + }
  47 +
  48 + int i;
  49 + unsigned char r, g, b;
  50 + unsigned int x, y;
  51 + for(y=0; y<y_size; y++)
  52 + for(x=0; x<x_size; x++)
  53 + {
  54 + //calculate the 1D index
  55 + i = y * x_size + x;
  56 +
  57 + r = buffer[i * 3 + 0];
  58 + g = buffer[i * 3 + 1];
  59 + b = buffer[i * 3 + 2];
  60 +
  61 + //set the image pixel
  62 + QColor color(r, g, b);
  63 + image.setPixel(x, y, color.rgb());
  64 + }
  65 +
  66 + if(!image.save(filename.c_str()))
  67 + std::cout<<"Error saving QImage."<<std::endl;
  68 +}
  69 +
  70 +#ifdef __CUDACC__
  71 +static void initBrewer()
  72 +{
  73 + //initialize the Brewer colormap
  74 +
  75 + //allocate CPU space
  76 + float4 cpuColorMap[BREWER_CTRL_PTS];
  77 +
  78 + //define control rtsPoints
  79 + cpuColorMap[0] = make_float4(0.192157f, 0.211765f, 0.584314f, 1.0f);
  80 + cpuColorMap[1] = make_float4(0.270588f, 0.458824f, 0.705882f, 1.0f);
  81 + cpuColorMap[2] = make_float4(0.454902f, 0.678431f, 0.819608f, 1.0f);
  82 + cpuColorMap[3] = make_float4(0.670588f, 0.85098f, 0.913725f, 1.0f);
  83 + cpuColorMap[4] = make_float4(0.878431f, 0.952941f, 0.972549f, 1.0f);
  84 + cpuColorMap[5] = make_float4(1.0f, 1.0f, 0.74902f, 1.0f);
  85 + cpuColorMap[6] = make_float4(0.996078f, 0.878431f, 0.564706f, 1.0f);
  86 + cpuColorMap[7] = make_float4(0.992157f, 0.682353f, 0.380392f, 1.0f);
  87 + cpuColorMap[8] = make_float4(0.956863f, 0.427451f, 0.262745f, 1.0f);
  88 + cpuColorMap[9] = make_float4(0.843137f, 0.188235f, 0.152941f, 1.0f);
  89 + cpuColorMap[10] = make_float4(0.647059f, 0.0f, 0.14902f, 1.0f);
  90 +
  91 +
  92 + int width = BREWER_CTRL_PTS;
  93 + int height = 0;
  94 +
  95 +
  96 + // allocate array and copy colormap data
  97 + cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32, 32, 32, 32, cudaChannelFormatKindFloat);
  98 +
  99 + HANDLE_ERROR(cudaMallocArray(&gpuBrewer, &channelDesc, width, height));
  100 +
  101 + HANDLE_ERROR(cudaMemcpyToArray(gpuBrewer, 0, 0, cpuColorMap, sizeof(float4)*width, cudaMemcpyHostToDevice));
  102 +
  103 + // set texture parameters
  104 + cudaTexBrewer.addressMode[0] = cudaAddressModeClamp;
  105 + //texBrewer.addressMode[1] = cudaAddressModeClamp;
  106 + cudaTexBrewer.filterMode = cudaFilterModeLinear;
  107 + cudaTexBrewer.normalized = true; // access with normalized texture coordinates
  108 +
  109 + // Bind the array to the texture
  110 + HANDLE_ERROR(cudaBindTextureToArray( cudaTexBrewer, gpuBrewer, channelDesc));
  111 +
  112 +}
  113 +
  114 +static void destroyBrewer()
  115 +{
  116 + HANDLE_ERROR(cudaFreeArray(gpuBrewer));
  117 +
  118 +}
  119 +
  120 +template<class T>
  121 +__global__ static void applyBrewer(T* gpuSource, unsigned char* gpuDest, unsigned int N, T minVal = 0, T maxVal = 1)
  122 +{
  123 +
  124 + int i = blockIdx.y * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x;
  125 + if(i >= N) return;
  126 +
  127 + //compute the normalized value on [minVal maxVal]
  128 + float a = (gpuSource[i] - minVal) / (maxVal - minVal);
  129 +
  130 + //lookup the color
  131 + float shift = 1.0/BREWER_CTRL_PTS;
  132 + float4 color = tex1D(cudaTexBrewer, a+shift);
  133 +
  134 + gpuDest[i * 3 + 0] = 255 * color.x;
  135 + gpuDest[i * 3 + 1] = 255 * color.y;
  136 + gpuDest[i * 3 + 2] = 255 * color.z;
  137 +}
  138 +
  139 +template<class T>
  140 +__global__ static void applyGrayscale(T* gpuSource, unsigned char* gpuDest, unsigned int N, T minVal = 0, T maxVal = 1)
  141 +{
  142 + int i = blockIdx.y * gridDim.x * blockDim.x + blockIdx.x * blockDim.x + threadIdx.x;
  143 + if(i >= N) return;
  144 +
  145 + //compute the normalized value on [minVal maxVal]
  146 + float a = (gpuSource[i] - minVal) / (maxVal - minVal);
  147 +
  148 + //threshold
  149 + if(a > 1.0)
  150 + a = 1.0;
  151 + if(a < 0.0)
  152 + a = 0.0;
  153 +
  154 + gpuDest[i * 3 + 0] = 255 * a;
  155 + gpuDest[i * 3 + 1] = 255 * a;
  156 + gpuDest[i * 3 + 2] = 255 * a;
  157 +}
  158 +
  159 +template<class T>
  160 +static void gpu2gpu(T* gpuSource, unsigned char* gpuDest, unsigned int nVals, T minVal = 0, T maxVal = 1, colormapType cm = cmGrayscale, int blockDim = 128)
  161 +{
  162 + //This function converts a scalar field on the GPU to a color image on the GPU
  163 + int gridX = (nVals + blockDim - 1)/blockDim;
  164 + int gridY = 1;
  165 + if(gridX > 65535)
  166 + {
  167 + gridY = (gridX + 65535 - 1) / 65535;
  168 + gridX = 65535;
  169 + }
  170 + dim3 dimGrid(gridX, gridY);
  171 + //int gridDim = (nVals + blockDim - 1)/blockDim;
  172 + if(cm == cmGrayscale)
  173 + applyGrayscale<<<dimGrid, blockDim>>>(gpuSource, gpuDest, nVals, minVal, maxVal);
  174 + else if(cm == cmBrewer)
  175 + {
  176 + initBrewer();
  177 + applyBrewer<<<dimGrid, blockDim>>>(gpuSource, gpuDest, nVals, minVal, maxVal);
  178 + //HANDLE_ERROR(cudaMemset(gpuDest, 0, sizeof(unsigned char) * nVals * 3));
  179 + destroyBrewer();
  180 + }
  181 +
  182 +}
  183 +
  184 +template<class T>
  185 +static void gpu2cpu(T* gpuSource, unsigned char* cpuDest, unsigned int nVals, T minVal, T maxVal, colormapType cm = cmGrayscale)
  186 +{
  187 + //this function converts a scalar field on the GPU to a color image on the CPU
  188 +
  189 + //first create the color image on the GPU
  190 +
  191 + //allocate GPU memory for the color image
  192 + unsigned char* gpuDest;
  193 + HANDLE_ERROR(cudaMalloc( (void**)&gpuDest, sizeof(unsigned char) * nVals * 3 ));
  194 +
  195 + //HANDLE_ERROR(cudaMemset(gpuSource, 0, sizeof(T) * nVals));
  196 +
  197 + //create the image on the gpu
  198 + gpu2gpu(gpuSource, gpuDest, nVals, minVal, maxVal, cm);
  199 +
  200 + //HANDLE_ERROR(cudaMemset(gpuDest, 0, sizeof(unsigned char) * nVals * 3));
  201 +
  202 + //copy the image from the GPU to the CPU
  203 + HANDLE_ERROR(cudaMemcpy(cpuDest, gpuDest, sizeof(unsigned char) * nVals * 3, cudaMemcpyDeviceToHost));
  204 +
  205 + HANDLE_ERROR(cudaFree( gpuDest ));
  206 +
  207 +}
  208 +
  209 +template<typename T>
  210 +static void gpu2image(T* gpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, T valMin, T valMax, colormapType cm = cmGrayscale)
  211 +{
  212 + //allocate a color buffer
  213 + unsigned char* cpuBuffer = NULL;
  214 + cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size);
  215 +
  216 + //memset(cpuBuffer, 255, sizeof(unsigned char) * 3 * x_size * y_size);
  217 +
  218 + //do the mapping
  219 + gpu2cpu<T>(gpuSource, cpuBuffer, x_size * y_size, valMin, valMax, cm);
  220 + std::cout<<"cpu scalar buffer: "<<(unsigned int)cpuBuffer[0]<<std::endl;
  221 +
  222 + //copy the buffer to an image
  223 + buffer2image(cpuBuffer, fileDest, x_size, y_size);
  224 +
  225 + free(cpuBuffer);
  226 +}
  227 +
  228 +#endif
  229 +
  230 +template<class T>
  231 +static void cpuApplyBrewer(T* cpuSource, unsigned char* cpuDest, unsigned int N, T minVal = 0, T maxVal = 1)
  232 +{
  233 + for(int i=0; i<N; i++)
  234 + {
  235 + //compute the normalized value on [minVal maxVal]
  236 + T v = cpuSource[i];
  237 + float a = (cpuSource[i] - minVal) / (maxVal - minVal);
  238 + if(a < 0) a = 0;
  239 + if(a > 1) a = 1;
  240 +
  241 + float c = a * (float)(BREWER_CTRL_PTS-1);
  242 + int ptLow = (int)c;
  243 + float m = c - (float)ptLow;
  244 +
  245 + float r, g, b;
  246 + if(ptLow == BREWER_CTRL_PTS - 1)
  247 + {
  248 + r = BREWERCP[ptLow * 4 + 0];
  249 + g = BREWERCP[ptLow * 4 + 1];
  250 + b = BREWERCP[ptLow * 4 + 2];
  251 + }
  252 + else
  253 + {
  254 + r = BREWERCP[ptLow * 4 + 0] * m + BREWERCP[ (ptLow+1) * 4 + 0] * (1.0 - m);
  255 + g = BREWERCP[ptLow * 4 + 1] * m + BREWERCP[ (ptLow+1) * 4 + 1] * (1.0 - m);
  256 + b = BREWERCP[ptLow * 4 + 2] * m + BREWERCP[ (ptLow+1) * 4 + 2] * (1.0 - m);
  257 + }
  258 +
  259 +
  260 + cpuDest[i * 3 + 0] = 255 * r;
  261 + cpuDest[i * 3 + 1] = 255 * g;
  262 + cpuDest[i * 3 + 2] = 255 * b;
  263 +
  264 + }
  265 +}
  266 +
  267 +template<class T>
  268 +static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T valMin, T valMax, colormapType cm = cmGrayscale)
  269 +{
  270 +
  271 + if(cm == cmBrewer)
  272 + cpuApplyBrewer(cpuSource, cpuDest, nVals, valMin, valMax);
  273 + else if(cm == cmGrayscale)
  274 + {
  275 + int i;
  276 + float a;
  277 + float range = valMax - valMin;
  278 + for(i = 0; i<nVals; i++)
  279 + {
  280 + //normalize to the range [valMin valMax]
  281 + a = (cpuSource[i] - valMin) / range;
  282 +
  283 + if(a < 0) a = 0.0;
  284 + if(a > 1) a = 1.0;
  285 +
  286 + cpuDest[i * 3 + 0] = 255 * a;
  287 + cpuDest[i * 3 + 1] = 255 * a;
  288 + cpuDest[i * 3 + 2] = 255 * a;
  289 + }
  290 + }
  291 +}
  292 +
  293 +template<class T>
  294 +static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale, bool positive = false)
  295 +{
  296 + //computes the max and min range automatically
  297 +
  298 + //find the largest magnitude value
  299 + T maxVal = abs(cpuSource[0]);
  300 + for(int i=0; i<nVals; i++)
  301 + {
  302 + if(abs(cpuSource[i]) > maxVal)
  303 + maxVal = abs(cpuSource[i]);
  304 + }
  305 +
  306 + if(positive)
  307 + cpu2cpu(cpuSource, cpuDest, nVals, (T)0.0, maxVal, cm);
  308 + else
  309 + cpu2cpu(cpuSource, cpuDest, nVals, -maxVal, maxVal, cm);
  310 +
  311 +}
  312 +
  313 +
  314 +
  315 +template<typename T>
  316 +static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, T valMin, T valMax, colormapType cm = cmGrayscale)
  317 +{
  318 + //allocate a color buffer
  319 + unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size);
  320 +
  321 + //do the mapping
  322 + cpu2cpu<T>(cpuSource, cpuBuffer, x_size * y_size, valMin, valMax, cm);
  323 +
  324 + //copy the buffer to an image
  325 + buffer2image(cpuBuffer, fileDest, x_size, y_size);
  326 +
  327 + free(cpuBuffer);
  328 +
  329 +}
  330 +
  331 +template<typename T>
  332 +static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, colormapType cm = cmGrayscale, bool positive = false)
  333 +{
  334 + //allocate a color buffer
  335 + unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size);
  336 +
  337 + //do the mapping
  338 + cpu2cpu<T>(cpuSource, cpuBuffer, x_size * y_size, cm, positive);
  339 +
  340 + //copy the buffer to an image
  341 + buffer2image(cpuBuffer, fileDest, x_size, y_size);
  342 +
  343 + free(cpuBuffer);
  344 +
  345 +}
  346 +
  347 +} //end namespace colormap and rts
  348 +
  349 +#endif
  350 +
... ...