scalarslice.cu
3.2 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
#include "scalarslice.h"
#include "rts/cuda/error.h"
#include "cublas_v2.h"
#include "rts/envi/envi.h"
scalarslice::scalarslice(int x, int y)
{
//set the resolution
R[0] = x;
R[1] = y;
//allocate memory on the GPU
HANDLE_ERROR(cudaMalloc( (void**)&S, sizeof(ptype) * x * y ));
//std::cout<<"Scalerslice created."<<std::endl;
}
scalarslice::scalarslice()
{
R[0] = R[1] = 0;
S = NULL;
//std::cout<<"Scalerslice created (default)."<<std::endl;
}
scalarslice::~scalarslice()
{
if(S != NULL)
HANDLE_ERROR(cudaFree(S));
S = NULL;
R[0] = R[1] = 0;
//std::cout<<"Scalerslice destroyed."<<std::endl;
}
//assignment operator
scalarslice & scalarslice::operator= (const scalarslice & rhs)
{
//de-allocate any existing GPU memory
if(S != NULL)
HANDLE_ERROR(cudaFree(S));
//copy the slice resolution
R[0] = rhs.R[0];
R[1] = rhs.R[1];
//allocate the necessary memory
HANDLE_ERROR(cudaMalloc(&S, sizeof(ptype) * R[0] * R[1]));
//copy the slice
HANDLE_ERROR(cudaMemcpy(S, rhs.S, sizeof(ptype) * R[0] * R[1], cudaMemcpyDeviceToDevice));
std::cout<<"Assignment operator."<<std::endl;
return *this;
}
void scalarslice::toImage(std::string filename, ptype vmin, ptype vmax, rts::colormapType cmap)
{
rts::gpu2image<ptype>(S, filename, R[0], R[1], vmin, vmax, cmap);
}
void scalarslice::toImage(std::string filename, bool positive, rts::colormapType cmap)
{
cublasStatus_t stat;
cublasHandle_t handle;
//create a CUBLAS handle
stat = cublasCreate(&handle);
if(stat != CUBLAS_STATUS_SUCCESS)
{
std::cout<<"CUBLAS Error: initialization failed"<<std::endl;
exit(1);
}
//find the index of the value with maximum magnitude
int N = R[0] * R[1];
int result;
#ifdef PRECISION_SINGLE
stat = cublasIsamax(handle, N, S, 1, &result);
#elif defined PRECISION_DOUBLE
stat = cublasIdamax(handle, N, S, 1, &result);
#endif
//adjust for 1-based indexing
result -= 1;
if(stat != CUBLAS_STATUS_SUCCESS)
{
std::cout<<"CUBLAS Error: failure finding maximum value."<<std::endl;
exit(1);
}
//retrieve the maximum value
ptype maxVal;
HANDLE_ERROR(cudaMemcpy(&maxVal, S + result, sizeof(ptype), cudaMemcpyDeviceToHost));
//destroy the CUBLAS handle
cublasDestroy(handle);
//output the image
if(positive)
toImage(filename, 0, maxVal, cmap);
else
toImage(filename, -abs(maxVal), abs(maxVal), cmap);
}
void scalarslice::toEnvi(std::string filename, ptype wavelength, bool append)
{
std::string mode;
if(append) mode = "a";
else mode = "w";
//open the ENVI file
EnviFile outfile(filename, mode);
//get the scalar slice from the GPU to the CPU
int memsize = sizeof(ptype) * R[0] * R[1];
ptype* cpuData = (ptype*) malloc( memsize );
HANDLE_ERROR(cudaMemcpy( cpuData, S, memsize, cudaMemcpyDeviceToHost));
//add a band to the ENVI file
outfile.addBand(cpuData, R[0], R[1], wavelength);
outfile.close();
}
void scalarslice::clear()
{
//this function sets the slice to zero
if(S != NULL)
//HANDLE_ERROR(cudaMalloc( (void**)&S, sizeof(ptype) * R[0] * R[1] ));
HANDLE_ERROR(cudaMemset(S, 0, sizeof(ptype) * R[0] * R[1]));
}