#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);

/// This function evaluates the theta-dependent multicale Pb given an multi-channel image

/// @param img is the multi-channel image
/// @param theta is the angle used for computing the 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

stim::image<float> func_mPb_theta(stim::image<float> img, float theta, int* r, float* alpha, int s){

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


	stim::image<float> mPb_theta(w, h, 1);               // allocate space for theta-dependent multicale Pb
	unsigned size = mPb_theta.size();                    // get the size of theta-dependent multicale Pb
	memset ( mPb_theta.data(), 0, size * sizeof(float)); // initialize all the pixels of theta-dependent multicale Pb to 0


	unsigned int N = w * h;						// get the number of pixels
	int sigma_n = 3; 							// set the number of standard deviations used to define the sigma

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

	stim::image<float> temp;                    // set the temporary image to store the addtion result

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

			//ss << "data_output/mPb_theta_slice"<< i*s + j << ".bmp";  // set the name for test result file
			//std::string sss = ss.str();
		    
			// get the gaussian gradient by convolving each image slice with the mask
			temp = gaussian_derivative_filter_odd(img.channel(i), r[i*s + j], sigma_n, theta);

			// (optional) output the test result of each slice
			//stim::cpu2image(temp.data(), sss, w, h, stim::cmBrewer);

			// multiply each gaussian gradient with its weight
			array_multiply(temp.data(), alpha[i*s + j], N);

			// add up all the weighted gaussian gradients
			array_add(mPb_theta.data(), temp.data(), mPb_theta.data(), N);

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

		}	
	}


	return mPb_theta;
}