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

stim::image<float> func_mPb(stim::image<float> lab, unsigned int theta_n, unsigned int w, unsigned int h){

	std::clock_t start;
	start = std::clock();

	//---------------pavel's suggesiton------------------------------------
	std::ostringstream ss;
	unsigned int N = w * h;
	stim::image<float> mPb_theta(w,h), mPb(w,h);
	unsigned size = mPb_theta.size();
	memset ( mPb.data(), 0, size * sizeof(float));

	float* ptr;
	ptr = (float*) malloc(size * sizeof(float) * theta_n);

	for (unsigned int n = 0; n < theta_n; n++){
	
		ss << "data_output/mPb_theta"<< n << "_conv2.bmp";
		float theta = 180 * ((float)n/theta_n); 

		mPb_theta = func_mPb_theta(lab, theta, w, h);
		//mPb_theta.load("101087.bmp");
		float* ptr_n = &ptr[ n * w * h * 1 ];
		mPb_theta.channel(0).data_noninterleaved(ptr_n);

		double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
		std::cout<<"mPb_theta_"<< theta <<" complished   time:"<< duration1 <<"s"<<'\n';


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

		std::string sss = ss.str();
		//stim::cpu2image(&ptr[idx], sss, w, h, stim::cmBrewer);
	 

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

			if(pixel > mPb.data()[i]){
				mPb.data()[i] = pixel;
			}
			
			else{
			}
		}

		

		ss.str("");
	}              

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

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

	return mPb; 

	//---------------my first method------------------------------------
	/*
	std::clock_t start;
	start = std::clock();

	stim::image<float> mPb_stack(w,h,theta_n), mPb(w,h), mPb_theta(w,h), A, B, temp;
	float* ptr[8];

	for (unsigned int n = 0; n < theta_n; n++){

		//int* x = new int(5);
        //int* y = x;
		//*y = 1;
		
		float theta = 180 * ((float)n/theta_n); 
		mPb_theta = func_mPb_theta(lab, theta, w, h);
		mPb_stack.getslice(n) = mPb_theta;
		float* ptr[n] = mPb_stack.getslice(n).data();

		double duration1 = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
		std::cout<<"mPb_theta, theta = "<< theta <<" time:"<< duration1 <<"s"<<'\n';
		

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

			*(ptr[n]+i) = mPb_theta.data()[i];
		
	       
			//float a = mPb_theta.data()[i];
			//float* B = ptr[n]+i;
			//A.data()[i] = mPb_theta.data()[i];
			//float* C = ptr[0]+1;
			//*C = 1;

	//		
		}
		stim::cpu2image(ptr[0], "data_output/mPb_theta.bmp", w, h, stim::cmBrewer);
	}

	for (unsigned long i = 0; i < N; i++){
		
		mPb.data()[i] = 0;
		for (unsigned int n = 0; n < theta_n; n++){

			float* ptr2 = ptr[i]+n;
			float temp = *ptr2;

			if(temp > mPb.data()[i]){
				mPb.data()[i] = temp;
			}	
			else{
			}
	    }
	}                  

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

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

	return mPb; */



}