96f9b10f
Laila Saadatifard
change the header...
|
1
2
|
#ifndef STIM_CUDA_CONV2SEP_H
#define STIM_CUDA_CONV2SEP_H
|
13fe3c84
Laila Saadatifard
update the stimli...
|
3
4
5
|
#include <iostream>
#include <cuda.h>
|
96f9b10f
Laila Saadatifard
change the header...
|
6
7
|
#include <stim/cuda/cudatools/devices.h>
#include <stim/cuda/cudatools/timer.h>
|
13fe3c84
Laila Saadatifard
update the stimli...
|
8
|
#include <stim/cuda/sharedmem.cuh>
|
96f9b10f
Laila Saadatifard
change the header...
|
9
|
#include <stim/cuda/cudatools/error.h>
|
13fe3c84
Laila Saadatifard
update the stimli...
|
10
11
12
13
14
15
16
|
#define pi 3.14159
namespace stim{
namespace cuda{
template<typename T>
|
96f9b10f
Laila Saadatifard
change the header...
|
17
18
|
__global__ void conv2sep_0(T* out, cudaTextureObject_t in, unsigned int x, unsigned int y,
T* kernel0, unsigned int k0){
|
13fe3c84
Laila Saadatifard
update the stimli...
|
19
|
|
96f9b10f
Laila Saadatifard
change the header...
|
20
21
|
//generate a pointer to shared memory (size will be specified as a kernel parameter)
extern __shared__ T s[];
|
13fe3c84
Laila Saadatifard
update the stimli...
|
22
|
|
96f9b10f
Laila Saadatifard
change the header...
|
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
|
int kr = k0/2; //calculate the kernel radius
//get a pointer to the gaussian in memory
T* g = (T*)&s[blockDim.x + 2 * kr];
//calculate the start point for this block
int bxi = blockIdx.x * blockDim.x;
int byi = blockIdx.y;
//copy the portion of the image necessary for this block to shared memory
stim::cuda::sharedMemcpy_tex2D(s, in, bxi - kr, byi, 2 * kr + blockDim.x, 1, threadIdx, blockDim);
//calculate the thread index
int ti = threadIdx.x;
//calculate the spatial coordinate for this thread
int xi = bxi + ti;
int yi = byi;
//use the first 2kr+1 threads to transfer the kernel to shared memory
if(ti < k0){
g[ti] = kernel0[ti];
}
//make sure that all writing to shared memory is done before continuing
__syncthreads();
//if the current pixel is outside of the image
if(xi > x || yi > y)
return;
//calculate the coordinates of the current thread in shared memory
int si = ti + kr;
T sum = 0; //running weighted sum across the kernel
//for each element of the kernel
for(int ki = -kr; ki <= kr; ki++){
sum += g[ki + kr] * s[si + ki];
|
13fe3c84
Laila Saadatifard
update the stimli...
|
65
|
}
|
96f9b10f
Laila Saadatifard
change the header...
|
66
67
68
|
//calculate the 1D image index for this thread
unsigned int i = byi * x + xi;
|
13fe3c84
Laila Saadatifard
update the stimli...
|
69
|
|
96f9b10f
Laila Saadatifard
change the header...
|
70
71
|
//output the result to global memory
out[i] = sum;
|
13fe3c84
Laila Saadatifard
update the stimli...
|
72
73
74
|
}
template<typename T>
|
96f9b10f
Laila Saadatifard
change the header...
|
75
76
77
78
79
80
81
82
83
84
85
86
87
88
|
__global__ void conv2sep_1(T* out, cudaTextureObject_t in, unsigned int x, unsigned int y,
T* kernel0, unsigned int k0){
//generate a pointer to shared memory (size will be specified as a kernel parameter)
extern __shared__ T s[];
int kr = k0/2; //calculate the kernel radius
//get a pointer to the gaussian in memory
T* g = (T*)&s[blockDim.y + 2 * kr];
//calculate the start point for this block
int bxi = blockIdx.x;
int byi = blockIdx.y * blockDim.y;
|
13fe3c84
Laila Saadatifard
update the stimli...
|
89
|
|
96f9b10f
Laila Saadatifard
change the header...
|
90
91
|
//copy the portion of the image necessary for this block to shared memory
stim::cuda::sharedMemcpy_tex2D(s, in, bxi, byi - kr, 1, 2 * kr + blockDim.y, threadIdx, blockDim);
|
13fe3c84
Laila Saadatifard
update the stimli...
|
92
|
|
96f9b10f
Laila Saadatifard
change the header...
|
93
94
|
//calculate the thread index
int ti = threadIdx.y;
|
13fe3c84
Laila Saadatifard
update the stimli...
|
95
|
|
96f9b10f
Laila Saadatifard
change the header...
|
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
//calculate the spatial coordinate for this thread
int xi = bxi;
int yi = byi + ti;
//use the first 2kr+1 threads to transfer the kernel to shared memory
if(ti < k0){
g[ti] = kernel0[ti];
}
//make sure that all writing to shared memory is done before continuing
__syncthreads();
//if the current pixel is outside of the image
if(xi > x || yi > y)
return;
|
13fe3c84
Laila Saadatifard
update the stimli...
|
113
|
|
96f9b10f
Laila Saadatifard
change the header...
|
114
115
|
//calculate the coordinates of the current thread in shared memory
int si = ti + kr;
|
13fe3c84
Laila Saadatifard
update the stimli...
|
116
|
|
96f9b10f
Laila Saadatifard
change the header...
|
117
|
T sum = 0; //running weighted sum across the kernel
|
13fe3c84
Laila Saadatifard
update the stimli...
|
118
|
|
96f9b10f
Laila Saadatifard
change the header...
|
119
120
121
122
123
124
125
126
127
128
129
|
//for each element of the kernel
for(int ki = -kr; ki <= kr; ki++){
sum += g[ki + kr] * s[si + ki];
}
//calculate the 1D image index for this thread
unsigned int i = yi * x + xi;
//output the result to global memory
out[i] = sum;
|
13fe3c84
Laila Saadatifard
update the stimli...
|
130
131
|
}
|
13fe3c84
Laila Saadatifard
update the stimli...
|
132
|
template<typename T>
|
96f9b10f
Laila Saadatifard
change the header...
|
133
134
135
136
|
void tex_conv2sep(T* out, unsigned int x, unsigned int y,
cudaTextureObject_t texObj, cudaArray* cuArray,
T* kernel0, unsigned int k0,
T* kernel1, unsigned int k1){
|
13fe3c84
Laila Saadatifard
update the stimli...
|
137
|
|
96f9b10f
Laila Saadatifard
change the header...
|
138
139
140
|
//get the maximum number of threads per block for the CUDA device
int max_threads = stim::maxThreadsPerBlock();
dim3 threads(max_threads, 1);
|
13fe3c84
Laila Saadatifard
update the stimli...
|
141
|
|
96f9b10f
Laila Saadatifard
change the header...
|
142
143
|
//calculate the number of blocks
dim3 blocks(x / threads.x + 1, y);
|
13fe3c84
Laila Saadatifard
update the stimli...
|
144
|
|
96f9b10f
Laila Saadatifard
change the header...
|
145
146
147
148
149
|
//calculate the shared memory used in the kernel
unsigned int pixel_bytes = max_threads * sizeof(T); //bytes devoted to pixel data being processed
unsigned int apron_bytes = k0/2 * sizeof(T); //bytes devoted to the apron on each side of the window
unsigned int gaussian_bytes = k0 * sizeof(T); //bytes devoted to memory used to store the pre-computed Gaussian window
unsigned int shared_bytes = pixel_bytes + 2 * apron_bytes + gaussian_bytes; //total number of bytes shared memory used
|
5cc0976c
David Mayerich
added separable c...
|
150
|
|
96f9b10f
Laila Saadatifard
change the header...
|
151
152
|
//blur the image along the x-axis
conv2sep_0<T> <<< blocks, threads, shared_bytes >>>(out, texObj, x, y, kernel0, k0);
|
13fe3c84
Laila Saadatifard
update the stimli...
|
153
|
|
96f9b10f
Laila Saadatifard
change the header...
|
154
155
156
157
158
159
160
161
162
163
164
165
|
// Copy the x-blurred data from global memory to the texture
cudaMemcpyToArray(cuArray, 0, 0, out, x * y * sizeof(T),
cudaMemcpyDeviceToDevice);
//transpose the block and thread dimensions
threads.x = 1;
threads.y = max_threads;
blocks.x = x;
blocks.y = y / threads.y + 1;
//blur the image along the y-axis
conv2sep_1<T> <<< blocks, threads, shared_bytes >>>(out, texObj, x, y, kernel1, k1);
|
13fe3c84
Laila Saadatifard
update the stimli...
|
166
167
168
|
}
|
13fe3c84
Laila Saadatifard
update the stimli...
|
169
|
template<typename T>
|
96f9b10f
Laila Saadatifard
change the header...
|
170
171
172
|
void gpu_conv2sep(T* image, unsigned int x, unsigned int y,
T* kernel0, unsigned int k0,
T* kernel1, unsigned int k1){
|
13fe3c84
Laila Saadatifard
update the stimli...
|
173
|
|
96f9b10f
Laila Saadatifard
change the header...
|
174
175
176
|
//get the number of pixels in the image
unsigned int pixels = x * y;
unsigned int bytes = sizeof(T) * pixels;
|
13fe3c84
Laila Saadatifard
update the stimli...
|
177
|
|
96f9b10f
Laila Saadatifard
change the header...
|
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
|
// Allocate CUDA array in device memory
//define a channel descriptor for a single 32-bit channel
cudaChannelFormatDesc channelDesc =
cudaCreateChannelDesc(32, 0, 0, 0,
cudaChannelFormatKindFloat);
cudaArray* cuArray; //declare the cuda array
cudaMallocArray(&cuArray, &channelDesc, x, y); //allocate the cuda array
// Copy the image data from global memory to the array
cudaMemcpyToArray(cuArray, 0, 0, image, bytes,
cudaMemcpyDeviceToDevice);
// Specify texture
struct cudaResourceDesc resDesc; //create a resource descriptor
memset(&resDesc, 0, sizeof(resDesc)); //set all values to zero
resDesc.resType = cudaResourceTypeArray; //specify the resource descriptor type
resDesc.res.array.array = cuArray; //add a pointer to the cuda array
|
13fe3c84
Laila Saadatifard
update the stimli...
|
196
|
|
96f9b10f
Laila Saadatifard
change the header...
|
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
|
// Specify texture object parameters
struct cudaTextureDesc texDesc; //create a texture descriptor
memset(&texDesc, 0, sizeof(texDesc)); //set all values in the texture descriptor to zero
texDesc.addressMode[0] = cudaAddressModeWrap; //use wrapping (around the edges)
texDesc.addressMode[1] = cudaAddressModeWrap;
texDesc.filterMode = cudaFilterModePoint; //use linear filtering
texDesc.readMode = cudaReadModeElementType; //reads data based on the element type (32-bit floats)
texDesc.normalizedCoords = 0; //not using normalized coordinates
// Create texture object
cudaTextureObject_t texObj = 0;
cudaCreateTextureObject(&texObj, &resDesc, &texDesc, NULL);
//call the texture version of the separable convolution function
tex_conv2sep(image, x, y, texObj, cuArray, kernel0, k0, kernel1, k1);
|
5cc0976c
David Mayerich
added separable c...
|
212
|
|
96f9b10f
Laila Saadatifard
change the header...
|
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
|
//free allocated memory
cudaFree(cuArray);
}
/// Applies a Gaussian blur to a 2D image stored on the CPU
template<typename T>
void cpu_conv2sep(T* image, unsigned int x, unsigned int y,
T* kernel0, unsigned int k0,
T* kernel1, unsigned int k1){
//get the number of pixels in the image
unsigned int pixels = x * y;
unsigned int bytes = sizeof(T) * pixels;
//---------Allocate Image---------
//allocate space on the GPU for the image
T* gpuI0;
HANDLE_ERROR(cudaMalloc(&gpuI0, bytes));
//copy the image data to the GPU
HANDLE_ERROR(cudaMemcpy(gpuI0, image, bytes, cudaMemcpyHostToDevice));
//---------Allocate Kernel--------
//allocate and copy the 0 (x) kernel
T* gpuK0;
HANDLE_ERROR(cudaMalloc(&gpuK0, k0 * sizeof(T)));
HANDLE_ERROR(cudaMemcpy(gpuK0, kernel0, k0 * sizeof(T), cudaMemcpyHostToDevice));
//allocate and copy the 1 (y) kernel
T* gpuK1;
HANDLE_ERROR(cudaMalloc(&gpuK1, k1 * sizeof(T)));
HANDLE_ERROR(cudaMemcpy(gpuK1, kernel1, k1 * sizeof(T), cudaMemcpyHostToDevice));
//run the GPU-based version of the algorithm
gpu_conv2sep<T>(gpuI0, x, y, gpuK0, k0, gpuK1, k1);
//copy the image data from the device
cudaMemcpy(image, gpuI0, bytes, cudaMemcpyDeviceToHost);
//free allocated memory
cudaFree(gpuI0);
|
13fe3c84
Laila Saadatifard
update the stimli...
|
255
256
|
}
|
6d30a707
Tianshu Cheng
add cuda/array_ad...
|
257
258
|
};
};
|
13fe3c84
Laila Saadatifard
update the stimli...
|
259
260
|
#endif
|