#ifndef STIM_CUDA_TPB_CUH
#define STIM_CUDA_TPB_CUH

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


//void array_multiply(float* lhs, float rhs, unsigned int N);
//void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);
//void chi_grad(float* img, float* cpu_copy, unsigned int w, unsigned int h, int r, unsigned int bin_n, unsigned int bin_size, float theta);

/// This function evaluates the tPb given a grayscale image

/// @param img is the multi-channel image
/// @param theta_n is the number of angles used for computing oriented chi-gradient
/// @param r is an array of radii for different scaled discs(filters)
/// @param alpha is is an array of weights for different scaled discs(filters)
/// @param s is the number of scales
/// @param K is the number of clusters

stim::image<float> tPb(stim::image<float> img, int* r, float* alpha, unsigned int theta_n, unsigned int bin_n, int s, unsigned K){

	unsigned int w = img.width();              // get the width of picture
	unsigned int h = img.height();             // get the height of picture
	unsigned int N = w * h;						// get the number of pixels

	stim::image<float> img_textons(w, h, 1, theta_n*2+1);               // allocate space for img_textons
	stim::image<float> img_texture(w, h, 1, 1);               // allocate space for img_texture
	stim::image<float> tPb_theta(w, h, 1, 1);               // allocate space for tPb_theta
	stim::image<float> tPb(w, h, 1, 1);						// allocate space for tPb
	unsigned size = tPb_theta.size();                    // get the size of tPb_theta
	memset (tPb.data(), 0, size * sizeof(float)); // initialize all the pixels of tPb to 0
	stim::image<float> temp(w, h, 1, 1);                    // set the temporary image to store the addtion result

	std::ostringstream ss;                      // (optional) set the stream to designate the test result file

	
	img_textons = textons(img, theta_n);

	img_texture = kmeans(img_textons, K);                // changing kmeans result into float type is required

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

	unsigned int max1 = img_texture.maxv();						// get the maximum of Pb used for normalization
	unsigned int bin_size = (max1 + 1)/bin_n;					// (whether"+1" or not depends on kmeans result)

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

		float theta = 180 * ((float)i/theta_n);              // calculate the even-splited angle for each tPb_theta

		memset (tPb_theta.data(), 0, size * sizeof(float)); // initialize all the pixels of tPb_theta to 0

		//ss << "data_output/0922tPb_theta"<< theta << ".bmp";  // set the name for test result file (optional)
		//std::string sss = ss.str();

		for (int j = 0; j < s; j++){
		    
			// get the chi-gradient by convolving each image slice with the mask
			chi_grad(img_texture.data(), temp.data(), w, h, r[j], bin_n, bin_size, theta);

			float max2 = temp.maxv();						// get the maximum of tPb_theta used for normalization
			array_multiply(temp.data(), 1/max2, N);		    // normalize the tPb_theta

			//output the test result of each slice (optional) 
			//stim::cpu2image(temp.data(), "data_output/tPb_slice0924_2.bmp", w, h, stim::cmBrewer);
			
			// multiply each chi-gradient with its weight
			array_multiply(temp.data(), alpha[j], N);
			
			// add up all the weighted chi-gradients
			array_add(tPb_theta.data(), temp.data(), tPb_theta.data(), N);


		}	
		
		//ss.str("");   //(optional) clear the space for stream

		for(unsigned long ti = 0; ti < N; ti++){

			if(tPb_theta.data()[ti] > tPb.data()[ti]){          //get the maximum value among all tPb_theta for ith pixel
				tPb.data()[ti] = tPb_theta.data()[ti];
			}
			
			else{
			}
		}
	}
	
	float max3 = tPb.maxv();						// get the maximum of tPb used for normalization
	array_multiply(tPb.data(), 1/max3, N);		    // normalize the tPb

	//output the test result of tPb (optional) 
	//stim::cpu2image(tPb.data(),  "data_output/tPb_0922.bmp", w, h, stim::cmBrewer);

	return tPb;
}

#endif