#ifndef STIM_CUDA_DG2_D2X_THETA_CONV2_CUH
#define STIM_CUDA_DG2_D2X_THETA_CONV2_CUH

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

#define SIGMA_N 3

/// This function evaluates the theta-dependent even-symmetric gaussian derivative gradient of an one-channel image

/// @param img is the one-channel image
/// @param r is an array of radii for different scaled discs(filters)
/// @param sigma_n is the number of standard deviations used to define the sigma
/// @param theta is angle used for computing the gradient

//void conv2(float* img, float* mask, float* cpu_copy, unsigned int w, unsigned int h, unsigned int M);
//void array_abs(float* img, unsigned int N);

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

	unsigned int w = image.width();    // get the width of picture
	unsigned int h = image.height();   // get the height of picture
	unsigned N = w * h;				   // get the number of pixels of picture
	
	int r = SIGMA_N * sigma;                         // set the radius of filter
	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> dG2_d2x_theta(w, h);      // allocate space for dG2_d2x_theta
	stim::image<float> mask_x(winsize, winsize);    // allocate space for x-axis-oriented filter
	stim::image<float> mask_r(winsize, winsize);    // allocate space for theta-oriented filter

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

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

			// create the x-oriented gaussian derivative filter mask_x
			mask_x.data()[j*winsize + i] = (-1) * (1 - pow(x, 2)) * exp((-1)*(pow(x, 2))/(2*pow(sigma, 2))) * exp((-1)*(pow(y, 2))/(2*pow(sigma, 2)));
		
		}
	}

	mask_r = mask_x.rotate(theta, r, r);
	//mask_r = mask_x.rotate(45, r, r);
	//stim::cpu2image(mask_r.data(), "data_output/mask_r_0919.bmp", winsize, winsize, stim::cmBrewer);  

	// do the 2D convolution with image and mask
	conv2(image.data(), mask_r.data(), dG2_d2x_theta.data(), w, h, winsize);
	array_abs(dG2_d2x_theta.data(), N);

	//stim::cpu2image(dG2_d2x_theta.data(), "data_output/dG2_d2x_theta_0919.bmp", w, h, stim::cmGrayscale); 
	
	return dG2_d2x_theta;
}

#endif