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

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

/// @param img is the multi-channel image
/// @param theta_n is the number of angles used for computing the gradient
/// @param r is an array of radii for different scaled discs(filters)
/// @param alpha is an array of weights for different scaled discs(filters)
/// @param s is the number of scales
stim::image<float> func_mPb(stim::image<float> img, unsigned int theta_n, int* r, float* alpha, int s){

	std::clock_t start;      // (optional) set the timer to calculate the total time 
	start = std::clock();    // (optional) set timer start point


	std::ostringstream ss;                              // (optional) set the stream to designate the test result file
	
	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> mPb_theta(w,h);                  // allocate space for theta-dependent multicale Pb (mPb_theta)
	stim::image<float> mPb(w,h);                        // allocate space for multicale Pb (mPb)

	unsigned size = mPb.size();							// get the size of mPb
	memset ( mPb.data(), 0, size * sizeof(float));      // initialize all the pixels of mPb to 0

	float* ptr;											   // set a pointer 
	ptr = (float*) malloc(size * sizeof(float) * theta_n); // this pointer points to a continuous space allocated to store all the mPb_theta

	for (unsigned int n = 0; n < theta_n; n++){
	
		ss << "data_output/mPb_theta"<< n << "_0911.bmp";    // (optional) set the name for test result file
		std::string sss = ss.str();							// (optional) 
		float theta = 180 * ((float)n/theta_n);              // calculate the even-splited angle for each mPb_theta

		mPb_theta = func_mPb_theta(img, theta, r, alpha, s); // calculate the mPb_theta

		float* ptr_n = &ptr[ n * w * h * 1 ];       // set a pointer which points to the space for each mPb_theta
		mPb_theta.data_noninterleaved(ptr_n);       // set this pointer to point to the each mPb_theta

		double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;             // (optional) calculate the time for generating each mPb_theta
		std::cout<<"mPb_theta_"<< theta <<" complished   time:"<< duration1 <<"s"<<'\n';   // (optional) show this time 


		unsigned long idx = n * w * h * 1;  //index for the nth mPb_theta

		
		stim::cpu2image(mPb_theta.data(), sss, w, h, stim::cmBrewer);	 // (optional) output the nth mPb_theta

		for(unsigned long i = 0; i < N; i++){
		
			float pixel = ptr[i+idx];           //get the ith pixel in nth mPb_theta

			if(pixel > mPb.data()[i]){          //get the maximum value among all mPb_theta for ith pixel
				mPb.data()[i] = pixel;
			}
			
			else{
			}
		}

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

	}              

	stim::cpu2image(mPb.data(), "data_output/mPb_500_0911_neat.bmp", w, h, stim::cmBrewer);  // output the mPb

	double duration2 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time 
	std::cout<<"total time:"<< duration2 <<"s"<<'\n';                      // (optional) show the total time 

	return mPb; 

}