#ifndef STIM_CUDA_DG1_CONV2_CUH
#define STIM_CUDA_DG1_CONV2_CUH

#include <stim/image/image.h>
//#include <cmath>
#include <stim/visualization/colormap.h>
//#include <stim/image/image_contour_detection.h>

#include <stim/cuda/templates/conv2sep.cuh>

#define SIGMA_N 3

/// This function generates the first-order gaussian derivative filter gx gy, 
/// convolves the image with gx gy, 
/// and returns an image class which channel(0) is Ix and channel(1) is Iy 

/// @param img is the one-channel image
/// @param sigma is the parameter for gaussian function

//void conv2_sep(float* img, unsigned int x, unsigned int y, float* kernel0, unsigned int k0, float* kernel1, unsigned int k1);
//void array_abs(float* img, unsigned int N);

stim::image<float> dG1_conv2(stim::image<float> image, int sigma){

	unsigned int w = image.width();    // get the width of picture
	unsigned int h = image.height();   // get the height of picture

	int r = SIGMA_N * sigma;  
	int winsize = 2 * SIGMA_N * sigma + 1;           // set the winsdow size of filter

	stim::image<float> I(w, h, 1, 2);      // allocate space for return image class
	stim::image<float> Ix(w, h);      // allocate space for Ix
	stim::image<float> Iy(w, h);      // allocate space for Iy
	Ix = image;  // initialize Ix
	Iy = image;  // initialize Iy

	float* array_x1;   
	array_x1 = new float[winsize];  //allocate space for the 1D x-oriented gaussian derivative filter array_x1 for gx
	float* array_y1;   
	array_y1 = new float[winsize];  //allocate space for the 1D y-oriented gaussian derivative filter array_y1 for gx
	float* array_x2;   
	array_x2 = new float[winsize];  //allocate space for the 1D x-oriented gaussian derivative filter array_x2 for gy
	float* array_y2;   
	array_y2 = new float[winsize];  //allocate space for the 1D y-oriented gaussian derivative filter array_y2 for gy


	for (int i = 0; i < winsize; i++){	

		int x = i - r;          //range of x
		int y = i - r;          //range of y

		// create the 1D x-oriented gaussian derivative filter array_x1 for gx
		array_x1[i] = (-1) * x * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
		// create the 1D y-oriented gaussian derivative filter array_y1  for gx
		array_y1[i] = exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
		// create the 1D x-oriented gaussian derivative filter array_x2 for gy
		array_x2[i] = exp((-1)*(pow(x, 2))/(2*pow(sigma, 2)));
		// create the 1D y-oriented gaussian derivative filter array_y2  for gy
		array_y2[i] = (-1) * y * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
	}

	//stim::cpu2image(array_x1, "data_output/array_x1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
	//stim::cpu2image(array_y1, "data_output/array_y1_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
	//stim::cpu2image(array_x2, "data_output/array_x2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result
	//stim::cpu2image(array_y2, "data_output/array_y2_0915.bmp", winsize, 1, stim::cmBrewer); // (optional) show the mask result

	// get Ix by convolving the image with gx
	//conv2_sep(Ix.data(), w, h, array_x1, winsize, array_y1, winsize);
	stim::cuda::cpu_conv2sep(Ix.data(), w, h, array_x1, winsize, array_y1, winsize);
	
	//stim::cpu2image(Ix.data(), "data_output/Ix_0915.bmp", w, h, stim::cmBrewer); 
	// get Iy by convolving the image with gy
	//conv2_sep(Iy.data(), w, h, array_x2, winsize, array_y2, winsize);
	stim::cuda::cpu_conv2sep(Iy.data(), w, h, array_x2, winsize, array_y2, winsize);
	
	//stim::cpu2image(Iy.data(), "data_output/Iy_0915.bmp", w, h, stim::cmBrewer); 

	delete [] array_x1;            //free the memory of array_x1
	delete [] array_y1;            //free the memory of array_y1
	delete [] array_x2;            //free the memory of array_x2
	delete [] array_y2;            //free the memory of array_y2

	I.set_channel(0, Ix.data());
	I.set_channel(1, Iy.data());

	return I;

}

#endif