segmentation / ivote3

#ifndef STIM_CUDA_GAUSSIAN_BLUR3_H
#define STIM_CUDA_GAUSSIAN_BLUR3_H

#include <iostream>
#include <cuda.h>
#include <stim/cuda/cudatools.h>
#include <stim/cuda/sharedmem.cuh>
#include <stim/cuda/cudatools/error.h>

#define pi	3.14159

				
		template<typename T>
		__global__ void blur_x(T* out, T* in, T sigma, unsigned int x, unsigned int y, unsigned int z){

			
			//calculate x,y,z coordinates for this thread
			int xi = blockIdx.x * blockDim.x + threadIdx.x;
			int yi = blockIdx.y * blockDim.y + threadIdx.y;
			int zi = blockIdx.z * blockDim.z + threadIdx.z;
			int i = zi * x * y + yi * x + xi;
			
			// calculate the kernel size
			int k_size = sigma * 4;
			
			//if the current pixel is outside of the image
			if(xi >= x || yi >= y || zi>=z)
				return;
			

			int gx, gy, gz, gi;
			T G;
			T sum = 0;		//running weighted sum across the kernel
			out[i] = 0;
				

			//for each element of the kernel
			for(int ki = -k_size; ki <= k_size; ki++){

				//calculate the gaussian value
				G = (1.0 / (sigma * sqrt(2 * pi))) * exp(-(ki*ki) / (2*sigma*sigma));
		

					
				//calculate the global coordinates for this point in the kernel
				gx = xi + ki;
				gy = yi;
				gz = zi;

				//make sure we are inside the bounds of the image
				if(gx >= 0 && gx < x ){
					gi = gz * x * y + gy * x + gx;
					//perform the integration
					sum += G * in[gi];
				}
		
			}
			
			//output the result to global memory
			out[i] = sum;
			
		}

	
		template<typename T>
		__global__ void blur_y(T* out, T* in, T sigma, unsigned int x, unsigned int y, unsigned int z){
			
			//calculate x,y,z coordinates for this thread
			int xi = blockIdx.x * blockDim.x + threadIdx.x;
			int yi = blockIdx.y * blockDim.y + threadIdx.y;
			int zi = blockIdx.z * blockDim.z + threadIdx.z;
			int i = zi * x * y + yi * x + xi;
			
			// calculate the kernel size
			int k_size = sigma * 4;
			
			//if the current pixel is outside of the image
			if(xi >= x || yi >= y || zi>=z)
				return;
			

			int gx, gy, gz, gi;
			T G;
			T sum = 0;		//running weighted sum across the kernel
			out[i] = 0;
				

			//for each element of the kernel
			for(int ki = -k_size; ki <= k_size; ki++){

				//calculate the gaussian value
				G = 1.0 / (sigma * sqrt(2 * pi)) * exp(-(ki*ki) / (2*sigma*sigma));
		

					
				//calculate the global coordinates for this point in the kernel
				gx = xi;
				gy = yi + ki;
				gz = zi;

				//make sure we are inside the bounds of the image
				if(gx >= 0 && gy >= 0 && gz >= 0 && gx < x && gy < y && gz < z){
					gi = gz * x * y + gy * x + gx;
					//perform the integration
					sum += G * in[gi];
				}
		
			}
			
			//output the result to global memory
			out[i] = sum;
		}
		
		template<typename T>
		__global__ void blur_z(T* out, T* in, T sigma, unsigned int x, unsigned int y, unsigned int z){
			
			//calculate x,y,z coordinates for this thread
			int xi = blockIdx.x * blockDim.x + threadIdx.x;
			int yi = blockIdx.y * blockDim.y + threadIdx.y;
			int zi = blockIdx.z * blockDim.z + threadIdx.z;
			int i = zi * x * y + yi * x + xi;
			
			// calculate the kernel size
			int k_size = sigma * 4;
			
			//if the current pixel is outside of the image
			if(xi >= x || yi >= y || zi>=z)
				return;
			

			int gx, gy, gz, gi;
			T G;
			T sum = 0;		//running weighted sum across the kernel
			out[i] = 0;
				

			//for each element of the kernel
			for(int ki = -k_size; ki <= k_size; ki++){

				//calculate the gaussian value
				G = 1.0 / (sigma * sqrt(2 * pi)) * exp(-(ki*ki) / (2*sigma*sigma));
		

					
				//calculate the global coordinates for this point in the kernel
				gx = xi;
				gy = yi;
				gz = zi + ki;

				//make sure we are inside the bounds of the image
				if(gx >= 0 && gy >= 0 && gz >= 0 && gx < x && gy < y && gz < z){
					gi = gz * x * y + gy * x + gx;
					//perform the integration
					sum += G * in[gi];
				}
		
			}
			
			//output the result to global memory
			out[i] = sum;
		}
		
		template<typename T>
		void gpu_gaussian_blur3(T* image, T sigma, unsigned int x, unsigned int y, unsigned int z){

			//get the number of pixels in the image
			unsigned int pixels = x * y * z;
			unsigned int bytes = sizeof(T) * pixels;

			int max_threads = stim::maxThreadsPerBlock();
			dim3 threads(512, 1, 1);
			dim3 blocks(x / threads.x + 1, y, z);

			//allocate temporary space on the GPU
			T* gpuIb_x;
			cudaMalloc(&gpuIb_x, bytes);

			//allocate temporary space on the GPU
			T* gpuIb_y;
			cudaMalloc(&gpuIb_y, bytes);

			// blur the original image in the x direction
			blur_x<T> <<< blocks, threads >>>(gpuIb_x, image, sigma, x, y, z);
			
			// blur the original image in the y direction
			blur_y<float> <<< blocks, threads >>>(gpuIb_y, gpuIb_x, sigma, x, y, z);
			
			// blur the original image in the z direction
			blur_z<float> <<< blocks, threads >>>(image, gpuIb_y, sigma, x, y, z);

			//cudaMemcpy(image, gpuIb_y, bytes, cudaMemcpyDeviceToDevice);

			cudaFree(gpuIb_x);
			cudaFree(gpuIb_y);
			
		}


		/// Applies a Gaussian blur to a 2D image stored on the CPU
		template<typename T>
		void cpu_gaussian_blur3(T* blur, T* image, T sigma, unsigned int x, unsigned int y, unsigned int z){

			
			//get the number of pixels in the image
			unsigned int pixels = x * y *z;
			unsigned int bytes = sizeof(T) * pixels;

			//---------Allocate Image---------
			//allocate space on the GPU for the image
			T* gpuI0;
			cudaMalloc(&gpuI0, bytes);			
			
			//copy the image data to the GPU
			cudaMemcpy(gpuI0, image, bytes, cudaMemcpyHostToDevice);

			//run the GPU-based version of the algorithm
			gpu_gaussian_blur3<T>(gpuI0, sigma, x, y, z);

			//copy the image data from the device
			cudaMemcpy(blur, gpuI0, bytes, cudaMemcpyDeviceToHost);

			//free allocated memory
			cudaFree(gpuI0);

			
		}


#endif