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

stim::image<float> func_mPb_theta(stim::image<float> lab, float theta, unsigned int w, unsigned int h){

	stim::image<float> mPb_theta(w, h, 1);

	stim::image<float> pic_light, pic_colora, pic_colorb;	
	pic_light = lab.channel(0);
	pic_colora = lab.channel(1);
	pic_colorb = lab.channel(2);

	unsigned int N = w * h;
	float sigma = 2; 
	unsigned int sigma_n = 3;
	unsigned r1 = 3;
	unsigned r2 = 5;
	unsigned r3 = 10;
	unsigned r4 = 20;
	float alpha[9] = {1,1,1,1,1,1,1,1,1};


	stim::image<float> l1,l2,l3,a1,a2,a3,b1,b2,b3;

	l1 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r1 * 2, theta, w, h);
	l2 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r2 * 2, theta, w, h);
	l3 = gaussian_derivative_filter_odd(pic_light, sigma, sigma_n, r3 * 2, theta, w, h);
	a1 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r2 * 2, theta, w, h);
	a2 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r3 * 2, theta, w, h);
	a3 = gaussian_derivative_filter_odd(pic_colora, sigma, sigma_n, r4 * 2, theta, w, h);
	b1 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r2 * 2, theta, w, h);
	b2 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r3 * 2, theta, w, h);
	b3 = gaussian_derivative_filter_odd(pic_colorb, sigma, sigma_n, r4 * 2, theta, w, h);

	for (unsigned i = 0; i<N; i++){

		mPb_theta.data()[i] = l1.data()[i] * alpha[0] +
							  l2.data()[i] * alpha[1] +
							  l3.data()[i] * alpha[2] +
							  a1.data()[i] * alpha[3] +
							  a2.data()[i] * alpha[4] +
							  a3.data()[i] * alpha[5] +
							  b1.data()[i] * alpha[6] +
							  b2.data()[i] * alpha[7] +
							  b3.data()[i] * alpha[8] ;
	
	}

	//stim::cpu2image(mPb_theta.data(), "data_output/cmap_mPb_theta0.bmp", w, h, stim::cmBrewer);


	//getch();

	return mPb_theta;
}