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

#define PI 3.1415926

void array_multiply(float* lhs, float rhs, unsigned int N);
void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
void array_abs(float* img, unsigned int N);

/// This function evaluates the theta-dependent odd 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

stim::image<float> Gd_odd(stim::image<float> image, int r, unsigned int sigma_n, float theta){

	float theta_r = (theta * PI)/180; //change angle unit from degree to rad

	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 winsize = 2 * r + 1;           // set the winsdow size of disc(filter)

	stim::image<float> I(w, h, 1, 2);       // allocate space for return image of Gd1
	stim::image<float> Ix(w, h);       // allocate space for Ix
	stim::image<float> Iy(w, h);       // allocate space for Iy
	stim::image<float> Gd_odd_theta(w, h);       // allocate space for Pb

	I = Gd1(image, r, sigma_n);  // calculate the Ix, Iy
	Ix = I.channel(0);
	Iy = I.channel(1);

	array_multiply(Ix.data(), cos(theta_r), N);								//Ix = Ix*cos(theta_r)
	array_multiply(Iy.data(), sin(theta_r), N);								//Iy = Iy*sin(theta_r)
	array_add(Ix.data(), Iy.data(), Gd_odd_theta.data(), N);				//Gd_odd_theta = Ix + Iy;
	array_abs(Gd_odd_theta.data(), N);		

	//stim::cpu2image(I.channel(0).data(), "data_output/Gd_odd_x_0919.bmp", w, h, stim::cmBrewer); 
	//stim::cpu2image(I.channel(1).data(), "data_output/Gd_odd_y_0919.bmp", w, h, stim::cmBrewer); 
	//stim::cpu2image(Gd_odd_theta.data(), "data_output/Gd_odd_theta_0919.bmp", w, h, stim::cmBrewer); 

	return Gd_odd_theta;

}