Blame view

stim/cuda/templates/conv2sep.cuh 8.47 KB
5cc0976c   David Mayerich   added separable c...
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
  #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>
  
  #define pi	3.14159
  
  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
ca99f951   David Mayerich   faster implementa...
33
34
  			//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);
5cc0976c   David Mayerich   added separable c...
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
  
  			//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
bf731970   Laila Saadatifard   fix some bugs in ...
53
  			if(xi >= x || yi >= y)
5cc0976c   David Mayerich   added separable c...
54
55
56
57
58
59
60
61
62
63
64
  				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++){
bf731970   Laila Saadatifard   fix some bugs in ...
65
  				sum += s[si + ki] * g[ki + kr];
5cc0976c   David Mayerich   added separable c...
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
  			}
  			
  			//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
ca99f951   David Mayerich   faster implementa...
92
93
  			//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);
5cc0976c   David Mayerich   added separable c...
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
  
  			//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);
  
59781ee3   Pavel Govyadinov   fixed a stask bug...
218
219
  			cudaDestroyTextureObject(texObj);
  
5cc0976c   David Mayerich   added separable c...
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
263
  		}
  
  		/// 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);
  		}
  		
  	};
  };
  
84eff8b1   Pavel Govyadinov   Merged only the n...
264
  #endif