codebase / stimlib

  #ifndef STIM_CUDA_CONV2_H
  #define STIM_CUDA_CONV2_H
  
  #include <iostream>
  #include <cuda.h>

  #include <stim/cuda/cudatools.h>

  #include <cmath>
  #include <algorithm>
  
  namespace stim{
  	namespace cuda{
  
  		template<typename T>

  		__global__ void cuda_conv2(T* mask, T* copy, cudaTextureObject_t texObj, unsigned int w, unsigned int h, unsigned int M){

  
  
  			//the radius of mask

  			int r = (M - 1)/2;

  
  
  			//calculate the 1D index for this thread
  			//int idx = blockIdx.x * blockDim.x + threadIdx.x;
  
  			//change 1D index to 2D cordinates
  			int i = blockIdx.x * blockDim.x + threadIdx.x;
  			int j = blockIdx.y;
  
  			int idx = j * w + i;
  			//unsigned long N = w * h;
  
  			if(i < w && j < h){
  

  				//copy[idx] = tex2D<float>(texObj, i+100, j+100);

  				//return;

  				tex2D<float>(texObj, (float)i/w, (float)j/h);

  
  				//allocate memory for result
  				T sum = 0;
  
  				//for (unsigned int y = max(j - r, 0); y <= min(j + r, h - 1); y++){
  
  					//for (unsigned int x = max(i - r, 0); x <= min(i + r, w - 1); x++){
  
  				for (int y = j - r; y <= j + r; y++){
  
  					for (int x = i - r; x <= i + r; x++){
  
  						//idx to mask cordinates(xx, yy)

  						int xx = x - (i - r);
  						int yy = y - (j - r);

  						sum += tex2D<T>(texObj, (float)x/w, (float)y/h) * mask[yy * M + xx];

  					}		
  				}
  				copy[idx] = sum;
  			 }
  		}
  		
  
  		template<typename T>
  		void gpu_conv2(T* img, T* mask, T* copy, unsigned int w, unsigned int h, unsigned M){
  
  			unsigned long N = w * h;
  
  			// 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, w, h);			//allocate the cuda array
  
  			// Copy the image data from global memory to the array
  			cudaMemcpyToArray(cuArray, 0, 0, img, N * sizeof(T),
  							  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]   = cudaAddressModeClamp;			//use wrapping (around the edges)
  			texDesc.addressMode[1]   = cudaAddressModeClamp;

  			texDesc.filterMode       = cudaFilterModePoint;		//use linear filtering
  			texDesc.readMode         = cudaReadModeElementType;		//reads data based on the element type (32-bit floats)

  			texDesc.normalizedCoords = 1;							//using normalized coordinates

  
  			// Create texture object
  			cudaTextureObject_t texObj = 0;
  			cudaCreateTextureObject(&texObj, &resDesc, &texDesc, NULL);
  
  			//get the maximum number of threads per block for the CUDA device
  			int threads = stim::maxThreadsPerBlock();
  
  			//calculate the number of blocks
  			dim3 blocks(w / threads + 1, h);
  
  			//call the kernel to do the multiplication

  			cuda_conv2 <<< blocks, threads >>>(mask, copy, texObj, w, h, M);

  
  		}
  
  		template<typename T>
  		void cpu_conv2(T* img, T* mask, T* cpu_copy, unsigned int w, unsigned int h, unsigned M){
  
  			unsigned long N = w * h;
  			//allocate memory on the GPU for the array
  			T* gpu_img; 
  			T* gpu_mask; 
  			T* gpu_copy;
  			HANDLE_ERROR( cudaMalloc( &gpu_img, N * sizeof(T) ) );
  			HANDLE_ERROR( cudaMalloc( &gpu_mask, pow(M, 2) * sizeof(T) ) );
  			HANDLE_ERROR( cudaMalloc( &gpu_copy, N * sizeof(T) ) );
  
  			//copy the array to the GPU
  			HANDLE_ERROR( cudaMemcpy( gpu_img, img, N * sizeof(T), cudaMemcpyHostToDevice) );
  			HANDLE_ERROR( cudaMemcpy( gpu_mask, mask, pow(M, 2) * sizeof(T), cudaMemcpyHostToDevice) );
  
  			//call the GPU version of this function
  			gpu_conv2<T>(gpu_img, gpu_mask ,gpu_copy, w, h, M);
  
  			//copy the array back to the CPU
  			HANDLE_ERROR( cudaMemcpy( cpu_copy, gpu_copy, N * sizeof(T), cudaMemcpyDeviceToHost) );
  
  			//free allocated memory
  			cudaFree(gpu_img);
  			cudaFree(gpu_mask);
  			cudaFree(gpu_copy);
  
  		}
  		
  	}
  }
  
  
  #endif