conv2sep.cuh
8.45 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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
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
255
256
257
258
259
260
261
262
#ifndef STIM_CUDA_CONV2SEP_H
#define STIM_CUDA_CONV2SEP_H
#include <iostream>
#include <cuda.h>
#include <stim/cuda/cudatools/devices.h>
#include <stim/cuda/cudatools/timer.h>
#include <stim/cuda/sharedmem.cuh>
#include <stim/cuda/cudatools/error.h>
namespace stim{
namespace cuda{
template<typename T>
__global__ void conv2sep_0(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.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<float, unsigned char>(s, in, bxi - kr, byi, 2 * kr + blockDim.x, 1, threadIdx, blockDim);
stim::cuda::sharedMemcpy_tex2D<float>(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 += s[si + ki] * g[ki + kr];
}
//calculate the 1D image index for this thread
unsigned int i = byi * x + xi;
//output the result to global memory
out[i] = sum;
}
template<typename T>
__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;
//copy the portion of the image necessary for this block to shared memory
//stim::cuda::sharedMemcpy_tex2D<float, unsigned char>(s, in, bxi, byi - kr, 1, 2 * kr + blockDim.y, threadIdx, blockDim);
stim::cuda::sharedMemcpy_tex2D<float>(s, in, bxi, byi - kr, 1, 2 * kr + blockDim.y, threadIdx, blockDim);
//calculate the thread index
int ti = threadIdx.y;
//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;
//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];
}
//calculate the 1D image index for this thread
unsigned int i = yi * x + xi;
//output the result to global memory
out[i] = sum;
}
template<typename T>
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){
//get the maximum number of threads per block for the CUDA device
int max_threads = stim::maxThreadsPerBlock();
dim3 threads(max_threads, 1);
//calculate the number of blocks
dim3 blocks(x / threads.x + 1, y);
//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
//blur the image along the x-axis
conv2sep_0<T> <<< blocks, threads, shared_bytes >>>(out, texObj, x, y, kernel0, k0);
// 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);
}
template<typename T>
void gpu_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 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
// 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);
//free allocated memory
cudaFree(cuArray);
cudaDestroyTextureObject(texObj);
}
/// 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);
}
};
};
#endif