#ifndef STIM_CUDA_DOWN_SAMPLE_H
#define STIM_CUDA_DOWN_SAMPLE_H

#include <iostream>
#include <cuda.h>
#include <stim/cuda/cudatools.h>
#include <stim/cuda/templates/gaussian_blur.cuh>

namespace stim{
	namespace cuda{

		template<typename T>
		__global__ void down_sample(T* gpuI, T* gpuI0, T resize, unsigned int x, unsigned int y){

			unsigned int sigma_ds = 1/resize;
			unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));
			unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));
			
			
			// calculate the 2D coordinates for this current thread.
			int xi = blockIdx.x * blockDim.x + threadIdx.x;
			int yi = blockIdx.y;
			// convert 2D coordinates to 1D
			int i = yi * x_ds + xi;
			
			if(xi< x_ds && yi< y_ds){

				int x_org = xi * sigma_ds ;
				int y_org = yi * sigma_ds ;
				int i_org = y_org * x + x_org;
				gpuI[i] = gpuI0[i_org];
			}

		}


		/// Applies a Gaussian blur to a 2D image stored on the GPU
		template<typename T>
		void gpu_down_sample(T* gpuI, T* gpuI0, T resize, unsigned int x, unsigned int y){

			
			unsigned int sigma_ds = 1/resize;
			unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));
			unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));
			
			//get the number of pixels in the image
//			unsigned int pixels_ds = x_ds * y_ds;
			
			unsigned int max_threads = stim::maxThreadsPerBlock();
			dim3 threads(max_threads, 1);
			dim3 blocks(x_ds/threads.x + (x_ds %threads.x == 0 ? 0:1) , y_ds);
			
			stim::cuda::gpu_gaussian_blur2<float>(gpuI0, sigma_ds,x ,y);
			
			//resample the image
			down_sample<float> <<< blocks, threads >>>(gpuI, gpuI0, resize, x, y);

		}

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

			//get the number of pixels in the image
			unsigned int pixels = x * y;
			unsigned int bytes = sizeof(T) * pixels;
			
			unsigned int sigma_ds = 1/resize;
			unsigned int x_ds = (x/sigma_ds + (x %sigma_ds == 0 ? 0:1));
			unsigned int y_ds = (y/sigma_ds + (y %sigma_ds == 0 ? 0:1));
			unsigned int bytes_ds = sizeof(T) * x_ds * y_ds;
			


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

			//allocate space on the GPU for the down sampled image
			T* gpuI;
			cudaMalloc(&gpuI, bytes_ds);

			//run the GPU-based version of the algorithm
			gpu_down_sample<T>(gpuI, gpuI0, resize, x, y);

			//copy the image data to the GPU
			cudaMemcpy(re_img, gpuI, bytes_ds, cudaMemcpyHostToDevice);

			cudaFree(gpuI0);
			cudeFree(gpuI);
		}
	
	}
}

#endif