#include <stim/image/image.h>
#include <cmath>
#include <stim/visualization/colormap.h>
#include <stim/image/image_contour_detection.h>
#include <iostream>
#include <stim/visualization/colormap.h>
#include <sstream>
/// calculate the cPb, tPb and mPb given a multi-channel image

void array_multiply(float* lhs, float rhs, unsigned int N);
void array_add(float* ptr1, float* ptr2, float* sum, unsigned int N);

int main()
{
	stim::image<float> image;                           // generate an image object

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

	//for(unsigned int i = 27; i < 31; i++){

	//ss1 << "3d_sample/b0"<< i << ".bmp";    // (optional) set the name for test result file
	//std::string sss1 = ss1.str();							// (optional) 

	//image.load("bone1001_3.bmp");	
	//image.load(sss1);	
	//image.load("101085.bmp");					// load the input image
	image.load("slice00_500_500.bmp");					// load the input image
	image = image.channel(0);                           // get the first channel of black-and-white image


	unsigned int w = image.width();		  // get the width of picture
	unsigned int h = image.height();      // get the height of picture
	unsigned int N = w * h;				  // get the number of pixels
	int c = image.channels();             // get the number if channels of picture
	int s = 3;						      // set the number of scales
	int sigma[3] = {1,2,3};              // set an array of sigmas for different scaled gaussian filters for cPb, the length of array is c*s
	int r[3] = {5,10,20};                // set an array of radii for different scaled discs(filters) for tPb, the length of array is c*s
	float alpha[3] = {1,1,1};             // set an array of weights for different scaled discs(filters)
	unsigned int theta_n = 8;             // set the number of angles used in filter orientation
	unsigned int bin_n = 16;              // set the number of bins used in chi-distance
	unsigned int K = 16;				  // set the number of cludters, K should be multiple of bin_n
	
	stim::image<float> img_cPb(w, h);				// allocate the space for cPb
	stim::image<float> img_tPb(w, h);				// allocate the space for tPb
	stim::image<float> img_mPb(w, h);				// allocate the space for tPb

	std::cout<<"imagesize: "<< w <<"*"<< h <<'\n';

	//*******************cPb********************
	std::clock_t start1;      // (optional) set the timer to calculate the total time 
	start1 = std::clock();    // (optional) set timer start point
	
	std::cout<<"begin cPb"<<'\n';     

	img_cPb = cPb(image, sigma, alpha, s);      
	
	// show the cPb (optional)	
    stim::cpu2image(img_cPb.data(), "data_output/img_cPb_1008.bmp", w, h, stim::cmBrewer); 
	//ss2 << "3d_sample/cPb0"<< i << ".bmp";    // (optional) set the name for test result file
	//std::string sss2 = ss2.str();

	//stim::cpu2image(img_cPb.data(), "data_output/bone_cPb3.1.bmp", w, h, stim::cmBrewer); 
	
	double duration1 = ( std::clock() - start1 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time 
	std::cout<<"cPb time: "<< duration1 <<"s"<<'\n';                      // (optional) show the total time 
	
	//ss1.str("");                  //(optional) clear the space for stream
	//ss2.str("");                  //(optional) clear the space for stream
	

	//*******************tPb********************
	
	std::cout<<"begin tPb"<<'\n';     
	

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

	img_tPb = tPb(image, r, alpha, theta_n, bin_n, s, K);
	
	// show the tPb (optional)
	// stim::cpu2image(img_tPb.data(), "data_output/bone_tPb_3.1.bmp", w, h, stim::cmBrewer); 
	stim::cpu2image(img_tPb.data(), "data_output/img_tPb_1008.bmp", w, h, stim::cmBrewer); 

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

	
	//*******************************************

	double duration3 = ( std::clock() - start1 ) / (double) CLOCKS_PER_SEC; // (optional) calculate the total time 
	std::cout<<"total running time: "<< duration3 <<"s"<<'\n';                      // (optional) show the total time 
	//}
	//******************mPb**********************
	// set parameters for linear combination
	
	float a = 1;
	float b = 0.5;

	// create mPb by linearly combined cPb and tPb
	array_multiply(img_cPb.data(), a, N);
	array_multiply(img_tPb.data(), b, N);
	array_add(img_cPb.data(), img_tPb.data(), img_mPb.data(), N);

	//ss2 << "3d_sample/mPb0"<< i << ".bmp";    // (optional) set the name for test result file
	//std::string sss2 = ss2.str();
	// show the mPb (optional)
	stim::cpu2image(img_mPb.data(), "data_output/img_mPb_1008.bmp", w, h, stim::cmBrewer); 

	//stim::cpu2image(img_mPb.data(), sss2, w, h, stim::cmBrewer);

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

	//}
	return 0;

}