codebase / stimlib

Commit 0174d8239377c3eada0bf7a09d5182112e7e8db3

Authored by dmayerich
1 parent 6e257ab3

bug fixes and exits

Showing 9 changed files with 120 additions and 405 deletions   Show diff stats
Inline Side-by-side
rts/cuda/error.h
  Show/Hide comments   View file @0174d82
1   -#include <stdio.h>
  1 +#include <stdio.h>
  2 +#include <iostream>
  3 +using namespace std;
2 4 #include "cuda_runtime.h"
3 5 #include "device_launch_parameters.h"
4 6 #include "cufft.h"
... ... @@ -24,15 +26,19 @@ static void CufftError( cufftResult err )
24 26 if (err != CUFFT_SUCCESS)
25 27 {
26 28 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");
  29 + cout<<"The plan parameter is not a valid handle."<<endl;
  30 + else if(err == CUFFT_ALLOC_FAILED)
  31 + cout<<"Allocation failed."<<endl;
  32 + else if(err == CUFFT_INVALID_VALUE)
  33 + cout<<"At least one of the parameters idata, odata, and direction is not valid."<<endl;
  34 + else if(err == CUFFT_INTERNAL_ERROR)
  35 + cout<<"An internal driver error was detected."<<endl;
  36 + else if(err == CUFFT_EXEC_FAILED)
  37 + cout<<"CUFFT failed to execute the transform on the GPU."<<endl;
  38 + else if(err == CUFFT_SETUP_FAILED)
  39 + cout<<"The CUFFT library failed to initialize."<<endl;
  40 + else
  41 + cout<<"Unknown error: "<<err<<endl;
36 42  
37 43 }
38 44 }
... ...
rts/envi/envi.h
  Show/Hide comments   View file @0174d82
... ... @@ -115,8 +115,9 @@ class EnviFile
115 115 f = fopen(filename.c_str(), "rb");
116 116 if(f == NULL)
117 117 {
118   - std::cout<<"ENVI Error: error testing file existance."<<std::endl;
119   - exit(1);
  118 + //std::cout<<"ENVI Error: error testing file existance."<<std::endl;
  119 + //exit(1);
  120 + return false;
120 121 }
121 122  
122 123 fseek(f, 0, SEEK_END);
... ... @@ -135,6 +136,7 @@ class EnviFile
135 136 int i;
136 137 int nBands = header.wavelength.size();
137 138 low = high = 0;
  139 + //std::cout<<"Wavelength to search: "<<wavelength<<" Number of Bands: "<<nBands<<std::endl;
138 140 for(i=0; i<nBands; i++)
139 141 {
140 142 if(header.wavelength[i] < wavelength)
... ... @@ -152,7 +154,9 @@ class EnviFile
152 154 if(header.wavelength[i] < header.wavelength[high])
153 155 high = i;
154 156 }
  157 + //std::cout<<"Low: "<<low<<" High: "<<high<<std::endl;
155 158 }
  159 + //exit(1);
156 160 }
157 161  
158 162 void interpolate_band(void* result, void* A, void* B, double a)
... ... @@ -197,9 +201,14 @@ public:
197 201  
198 202 void open(std::string file_name, std::string header_name, std::string file_mode = "r")
199 203 {
  204 +
200 205 //load the header file
201 206 // if the file doesn't exist, that is okay unless it is read-only
202   - if(!header.load(header_name))
  207 +
  208 + //if the file is being written, create a new header
  209 + if(file_mode == "w")
  210 + header.name = header_name;
  211 + else if(!header.load(header_name))
203 212 {
204 213 if(file_mode == "r")
205 214 {
... ... @@ -224,7 +233,8 @@ public:
224 233  
225 234  
226 235 //open the data file
227   - std::string tmpmode = mode+"b";
  236 + std::string tmpmode = mode+"b";
  237 + //std::cout<<"Mode: "<<tmpmode<<std::endl;
228 238 data = fopen(file_name.c_str(), tmpmode.c_str());
229 239 if(data == NULL)
230 240 {
... ... @@ -363,8 +373,8 @@ public:
363 373 size_t n_read = fread(ptr, getElementSize(), header.samples * header.lines, data);
364 374 if(n_read != header.samples * header.lines)
365 375 {
366   - std::cout<<"ENVI Error: Error reading band data."<<std::endl;
367   - return;
  376 + std::cout<<"ENVI Error -- Error reading band data: "<<n_read<<"/"<<samples() * lines()<<" pixels read"<<std::endl;
  377 + exit(1);
368 378 }
369 379  
370 380 return;
... ... @@ -383,12 +393,13 @@ public:
383 393 //get the low and high band numbers
384 394 unsigned int low, high;
385 395 get_surrounding_bands(wavelength, low, high);
  396 + //std::cout<<"High: "<<high<<" Low: "<<low<<std::endl;
386 397  
387 398 //calculate the position between the two bands
388 399 double a;
389 400 if(low == high)
390 401 a = 1.0;
391   - else
  402 + else
392 403 a = (wavelength - header.wavelength[low]) / (header.wavelength[high] - header.wavelength[low]);
393 404  
394 405 //read both bands
... ...
rts/envi/envi_header.h
  Show/Hide comments   View file @0174d82
... ... @@ -187,16 +187,20 @@ struct EnviHeader
187 187 std::vector<double> get_double_seq(std::string token, std::string sequence)
188 188 {
189 189 //this function returns a sequence of comma-delimited strings
190   - std::vector<double> result;
  190 + std::vector<double> result;
  191 +
  192 + double fentry;
191 193  
192 194 std::string entry;
193 195 size_t i;
194 196 do
195 197 {
196 198 i = sequence.find_first_of(',');
197   - entry = sequence.substr(0, i);
  199 + entry = sequence.substr(0, i);
  200 + fentry = atof(entry.c_str());
198 201 sequence = sequence.substr(i+1);
199   - result.push_back(atof(entry.c_str()));
  202 + result.push_back(atof(entry.c_str()));
  203 + //std::cout<<entry<<" ";
200 204 }while(i != std::string::npos);
201 205  
202 206 return result;
... ... @@ -222,12 +226,14 @@ struct EnviHeader
222 226 }
223 227  
224 228 //for each line in the file, get the token
225   - std::string token;
  229 + std::string token;
  230 +
  231 + //get a line
  232 + getline(file, line);
226 233 while(file)
227 234 {
228 235  
229   - //get a line
230   - getline(file, line);
  236 +
231 237  
232 238 //get the token
233 239 token = get_token(line);
... ... @@ -290,8 +296,19 @@ struct EnviHeader
290 296 else if(token == "y start")
291 297 y_start = atof(get_data_str(line).c_str());
292 298 else if(token == "wavelength units")
293   - wavelength_units = get_data_str(line);
294   - }
  299 + wavelength_units = get_data_str(line);
  300 +
  301 + //get the next line
  302 + getline(file, line);
  303 + }
  304 +
  305 + //make sure the number of bands matches the number of wavelengths
  306 + unsigned int wavelengths = wavelength.size();
  307 + if(bands != wavelengths)
  308 + {
  309 + std::cout<<"ENVI Header Error -- Number of wavelengths doesn't match the number of bands. Bands = "<<bands<<", Wavelengths = "<<wavelength.size()<<std::endl;
  310 + exit(1);
  311 + }
295 312  
296 313 //close the file
297 314 file.close();
... ...
rts/math/point.h
  Show/Hide comments   View file @0174d82
1 1 #ifndef RTS_rtsPoint_H
2 2 #define RTS_rtsPoint_H
3 3  
4   -//#include "rts/vector.h"
  4 +#include "rts/math/vector.h"
5 5 #include <string.h>
  6 +#include "rts/cuda/callable.h"
6 7  
7 8 namespace rts
8 9 {
... ...
rts/math/quad.h
  Show/Hide comments   View file @0174d82
... ... @@ -5,6 +5,7 @@
5 5 #include "rts/cuda/callable.h"
6 6 #include "rts/math/vector.h"
7 7 #include "rts/math/point.h"
  8 +#include "rts/math/triangle.h"
8 9 #include <iostream>
9 10  
10 11 namespace rts{
... ... @@ -78,28 +79,6 @@ struct rtsQuad
78 79  
79 80 }
80 81  
81   - /*CUDA_CALLABLE rtsQuad(rts::rtsPoint<T, N> p, rts::vector<T, N> x, rts::vector<T, N> y, T sx, T sy)
82   - {
83   - //This constructor creates a rtsQuad given a position, orientation, and size
84   - // p = center position of the rtsQuad
85   - // x = x-axis for the rtsQuadangle
86   - // y = y-axis for the rtsQuadangle
87   - // sx = size of the rtsQuad along the A-B axis
88   - // sy = size of the rtsQuad along the A-C axis
89   -
90   - //normalize x and y
91   - rts::vector<T, N> nx = x.norm();
92   - rts::vector<T, N> ny = y.norm();
93   -
94   - //compute X and Y
95   - X = sx * x;
96   - Y = sy * y;
97   -
98   - //compute A
99   - A = p - 0.5 * X - 0.5 * Y;
100   -
101   - }*/
102   -
103 82 CUDA_CALLABLE rts::rtsPoint<T, N> p(T a, T b)
104 83 {
105 84 rts::rtsPoint<T, N> result;
... ... @@ -143,6 +122,24 @@ struct rtsQuad
143 122 return result;
144 123  
145 124 }
  125 +
  126 + CUDA_CALLABLE T dist(rtsPoint<T, N> p)
  127 + {
  128 + //compute the distance between a point and this quad
  129 +
  130 + //first break the quad up into two triangles
  131 + rtsTriangle<T, N> T0(A, A+X, A+Y);
  132 + rtsTriangle<T, N> T1(A+X+Y, A+X, A+Y);
  133 +
  134 +
  135 + ptype d0 = T0.dist(p);
  136 + ptype d1 = T1.dist(p);
  137 +
  138 + if(d0 < d1)
  139 + return d0;
  140 + else
  141 + return d1;
  142 + }
146 143 };
147 144  
148 145 } //end namespace rts
... ...
rts/math/spherical_bessel.h
  Show/Hide comments   View file @0174d82
... ... @@ -5,7 +5,8 @@
5 5  
6 6 namespace rts{
7 7  
8   -#define RTS_BESSEL_CONVERGENCE_FLOAT 0.0145
  8 +#define RTS_BESSEL_CONVERGENCE_MIN 0.0145
  9 +#define RTS_BESSEL_CONVERGENCE_MAX 0.4
9 10 #define RTS_BESSEL_MAXIMUM_FLOAT -1e33
10 11  
11 12 template <typename T>
... ... @@ -109,7 +110,7 @@ CUDA_CALLABLE void shift_sbesselj(int n, T x, T* b)//, T stability = 1.4)
109 110  
110 111 //if(n > stability*x)
111 112 if(n > real(x))
112   - if(real(bnew) < RTS_BESSEL_CONVERGENCE_FLOAT)
  113 + if(real(bnew) < RTS_BESSEL_CONVERGENCE_MIN || real(bnew) > RTS_BESSEL_CONVERGENCE_MAX)
113 114 bnew = 0.0;
114 115  
115 116 //shift and add the new value to the array
... ...
rts/math/vector.h
  Show/Hide comments   View file @0174d82
1 1 #ifndef RTS_VECTOR_H
2 2 #define RTS_VECTOR_H
3 3  
4   -#include <iostream>
  4 +#include <iostream>
  5 +#include <cmath>
  6 +#include <sstream>
5 7 //#include "rts/point.h"
  8 +#include "rts/cuda/callable.h"
6 9  
7 10 namespace rts
8 11 {
... ...
rts/optics/material.h
  Show/Hide comments   View file @0174d82
... ... @@ -113,8 +113,17 @@ namespace rts{
113 113 //read the entry from an input string
114 114 for(int i=0; i<valueList.size(); i++)
115 115 {
  116 +
  117 + char c;
  118 + while(ss.peek() < '0' || ss.peek() > '9')
  119 + {
  120 + //cout<<"bad char"<<endl;
  121 + ss.ignore();
  122 + }
  123 +
116 124 //retrieve the value
117 125 ss>>val;
  126 + //cout<<val<<endl;
118 127  
119 128 //store the value in the appropriate location
120 129 switch(valueList[i])
... ... @@ -184,6 +193,7 @@ namespace rts{
184 193 template <class T>
185 194 class material
186 195 {
  196 + std::string name;
187 197 //dispersion (refractive index as a function of wavelength)
188 198 std::vector< refIndex<T> > dispersion;
189 199  
... ... @@ -218,7 +228,22 @@ namespace rts{
218 228 }
219 229  
220 230 public:
221   -
  231 + std::string getName()
  232 + {
  233 + return name;
  234 + }
  235 + void setName(std::string n)
  236 + {
  237 + name = n;
  238 + }
  239 + T getN0()
  240 + {
  241 + return n0;
  242 + }
  243 + void setN0(T n)
  244 + {
  245 + n0 = n;
  246 + }
222 247 unsigned int nSamples()
223 248 {
224 249 return dispersion.size();
... ... @@ -367,8 +392,10 @@ namespace rts{
367 392 fftw_free(Chi1);
368 393  
369 394 return newM;
370   -#endif
  395 +#else
  396 + std::cout<<"MATERIAL Error: Must have FFTW in order to compute Kramers-Kronig."<<std::endl;
371 397 return material<T>();
  398 +#endif
372 399 }
373 400  
374 401 material(T lambda = 1.0, T n = 1.4, T k = 0.0)
... ... @@ -387,10 +414,12 @@ namespace rts{
387 414 material(std::string filename, std::string format = "", T scaleA = 1.0)
388 415 {
389 416 fromFile(filename, format);
  417 + n0 = 0;
390 418 }
391 419  
392 420 void fromFile(std::string filename, std::string format = "", T scaleA = 1.0)
393 421 {
  422 + name = filename;
394 423 //clear any previous values
395 424 dispersion.clear();
396 425  
... ... @@ -401,7 +430,7 @@ namespace rts{
401 430  
402 431 if(!ifs.is_open())
403 432 {
404   - std::cout<<"Error: material file not found"<<std::endl;
  433 + std::cout<<"Material Error -- file not found: "<<filename<<std::endl;
405 434 exit(1);
406 435 }
407 436  
... ... @@ -449,7 +478,7 @@ namespace rts{
449 478  
450 479 entryType<T> entry(format);
451 480  
452   - std::cout<<"Loading material with format: "<<format<<std::endl;
  481 + //std::cout<<"Loading material with format: "<<format<<std::endl;
453 482  
454 483 T lambda, n, k;
455 484 while(!ss.eof())
... ...
rts/tools/colormap.h deleted
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   -