#include "bbUtilitiesBuildBinaryImage.h"
#include "bbUtilitiesPackage.h"
namespace bbUtilities
{

BBTK_ADD_BLACK_BOX_TO_PACKAGE(Utilities,BuildBinaryImage)
//BBTK_USER_BLACK_BOX_IMPLEMENTATION(BuildBinaryImage,bbtk::AtomicBlackBox);
BBTK_BLACK_BOX_IMPLEMENTATION(BuildBinaryImage,bbtk::AtomicBlackBox);
void BuildBinaryImage::Process()
{
	//Getting the vectors
	/*
	std::vector< double > xv = bbGetInputInX();
	std::vector< double > yv = bbGetInputInY();
	std::vector< double > zv = bbGetInputInZ();
    */
	//thresholds
	//upper = bbGetInputTUpper();
	//lower = bbGetInputTLower();

	//-MANUAL AXE POINTS-----------------------------------------------
	std::vector< double > xv;
	std::vector< double > yv;
	std::vector< double > zv;
	
	//-MANUAL AXE POINTS-----------------------------------------------
	int RADIO=8;
	std::vector< int > firstpoint;

	
	//-PATH DEFINITION---------------------------------------------------
	std::string path= "C:/CAT08/Desarrollo/Utilidades/Utilities/data/";
    std::string TR= "TR";
	//std::string LC= "Eje_LC_MIDDLE";
	std::string dataset= "D00";
    std::string vessel= "V1";

    //-SPLINE POINTS FILE---------------------------------------------
	//std::string cpath=path + RC + "_"+dataset + "_"+serie +".txt";
	std::string cpath=path +TR +"_"+dataset+"_"+ vessel +".txt";
    //-------------------------------------------------------------

	//-READING THE FILE---------------------------------------------------------
     
	std::string line;
	std::ifstream file;

	file.open( cpath.c_str() );	
	if(file.is_open())
	{
		
		std::getline(file,line);
		//NUMBER OF SPLINE POINTS-------------------
		int nPoints=atoi(line.c_str());

		int i=0;
		while(!file.eof() && i<nPoints)
		{
			std::getline(file,line);
			int pos=line.find(" ");
			int pos2=line.find_last_of(" ");
			int size=line.size();
			std::string xs=line.substr(0,pos);
			std::string ys=line.substr(pos+1,pos2);
			std::string zs=line.substr(pos2+1,size);
			double x0=atof(xs.c_str());
			double y0=atof(ys.c_str());
			double z0=atof(zs.c_str());
			//FIRST POINT FOR MARACAS ALGORITHM
			if(i >= nPoints/2){
				firstpoint.push_back(x0);
				firstpoint.push_back(y0);
				firstpoint.push_back(z0);
			}
			xv.push_back(x0);
			yv.push_back(y0);
			zv.push_back(z0);

			i++;

		}
	}
	file.close();
	
	//--GETTING THE  IMAGE----------------
	vtkImageData* img = bbGetInputImageIn();
    
	//---BINARY IMAGE INITIALIZATION
	if(imgbin==NULL)
	   imgbin = vtkImageData::New();

	//TYPE
    imgbin->SetScalarTypeToUnsignedShort();
    imgbin->SetNumberOfScalarComponents(1);
	//SPACING
	double* spc = img->GetSpacing();
    imgbin->SetSpacing(spc[0], spc[1], spc[2]);
    imgbin->SetOrigin(0, 0, 0);
    //SIZE
	int ext[6];
	img->GetExtent(ext);
    imgbin->SetExtent(ext);
    //GETTING MEMORY
    imgbin->AllocateScalars();

	//IMAGE SIZE
	int extx = ext[1]-ext[0]+1, exty = ext[3]-ext[2]+1, extz = ext[5]-ext[4]+1;

	//GETTING THE IMAGE IN BLACK
	unsigned short* imgbinpoint = (unsigned short*)imgbin->GetScalarPointer();
	//The whole volume is going black
	for(int i = 0; i < extx*exty*extz; i++){

		imgbinpoint[i]=0;

	}
	
	
	unsigned short* imgpoint = (unsigned short*)img->GetScalarPointer();
	//Drawing the cells according to the vector points x, y , z
	
	//CREATING SPHERES IN EACH POINT
	for(int i = 0; i < xv.size();i++)
	
	{
		//---SPLINE POINT---------------------------
		double xp = xv[i], yp = yv[i], zp = zv[i];
		
		//---CUBE COORDINATES--------------------
		int extint[6];		
		extint[0]=(xp-RADIO); 
		extint[1]=(xp+RADIO); 
		extint[2]=(yp-RADIO); 
		extint[3]=(yp+RADIO); 
		extint[4]=(zp-RADIO); 
		extint[5]=(zp+RADIO);
		
		//---EXTRACTING THE VOI----------
		vtkExtractVOI* evoi = vtkExtractVOI::New();
		evoi->SetInput(img);
		evoi->SetVOI(extint);	
		evoi->UpdateWholeExtent();
		evoi->Update();
		evoi->GetOutput()->UpdateInformation();
		evoi->GetOutput()->SetUpdateExtent(evoi->GetOutput()->GetWholeExtent());
		evoi->GetOutput()->Update();
		
		//---THE VOI----------
		vtkImageData* imgVOI = evoi->GetOutput();

		
		//------CALCUALTING THE ADAPTATIVE THRESHOLD AND STANDAR DESVIATION---------------
		//Calculate adaptative Threshold, average and standard deviation
		unsigned short* ptrimgvoi;
		double greyavg=0, greystd = 0;
		int greypoints=0;
        

		//-------SPHERE SIZE FOR THE ADAPTATIVE THRESHOLD---------------
		int rad = 2;

		//------ADAPTATIVE THRESHOLD---------------
		
		for(int i = xp-rad; i < xp+rad; i++){
			for(int j = yp-rad; j < yp+rad; j++){
				for(int k = zp-rad; k < zp+rad; k++){
					if(pow(i-xp,2)+pow(j-yp,2)+pow(k-zp,2)<pow(rad, 2)){
						ptrimgvoi = (unsigned short*)imgVOI->GetScalarPointer(i, j, k);
						//borrame
						//std::cout << "NG: " << (double)*ptrimgvoi << std::endl;
						greyavg+= (double)*ptrimgvoi;
						greypoints++;
					}
				}
			}
		}
		/*
		for(int i=0 ; i <= 8*rad*rad*rad ; i++)
		{
           if(pow(i-xp,2)+pow(j-yp,2)+pow(k-zp,2)<pow(rad, 2))
		   {
						ptrimgvoi = (unsigned short*)imgVOI->GetScalarPointer(i, j, k);
						greyavg+= (double)*ptrimgvoi;
						greypoints++;
		   }
		}
		*/
       //------ADAPTATIVE STANDAR DESVIATION---------------
		if(greypoints!=0){
			greyavg=greyavg/greypoints;
			for(int i = xp-rad; i < xp+rad; i++){
				for(int j = yp-rad; j < yp+rad; j++){
					for(int k = zp-rad; k < zp+rad; k++){
						if(pow(i-xp,2)+pow(j-yp,2)+pow(k-zp,2)<pow(rad, 2)){
							ptrimgvoi = (unsigned short*)imgVOI->GetScalarPointer(i, j, k);
							greystd+= pow((double)*ptrimgvoi-greyavg, 2);						
						}
					}
				}
			}            
			greystd = sqrt(greystd/greypoints);
		}	
		//------APPLYING THE  THRESHOLD---------------

		vtkImageThreshold *thresh2 = vtkImageThreshold::New();
		thresh2->SetInValue(255.0);
		thresh2->SetOutputScalarTypeToUnsignedShort();
  		thresh2->SetOutValue(0.0);
		thresh2->ReplaceInOn();
		thresh2->ReplaceOutOn();
		thresh2->SetInput(evoi->GetOutput());	
		//thresh2->ThresholdBetween(lower, upper);

		//ptrimgvoi = (unsigned short*)imgVOI->GetScalarPointer(xp, yp, zp);
		//greyavg = (double)*ptrimgvoi;
		thresh2->ThresholdBetween(greyavg-(greystd), greyavg+ (greystd));
		
		//thresh2->ThresholdBetween(1174, 2000);
	
		thresh2->UpdateWholeExtent();	
		thresh2->Update();


		//----CASTING-----------------
		vtkImageCast* cast2 = vtkImageCast::New();
 		cast2->SetInput(thresh2->GetOutput());
		cast2->SetOutputScalarTypeToUnsignedChar();	
		cast2->Update();

		//---SEGMENTATION CONNECTIVITY-------------

		vtkImageSeedConnectivity* connect2 = vtkImageSeedConnectivity::New();
		connect2->SetInput(cast2->GetOutput());
		connect2->SetInputConnectValue(255);
		connect2->SetOutputConnectedValue(255);
		connect2->SetOutputUnconnectedValue(0);
		connect2->RemoveAllSeeds();
		connect2->AddSeed(xp, yp, zp );
		connect2->Update();
        
		//----CASTING-----------------
		vtkImageCast* cast4 = vtkImageCast::New();
		cast4->SetInput(connect2->GetOutput());
		cast4->SetOutputScalarTypeToUnsignedShort();
		cast4->Update();
	
		//----COPYING THE VOLUME TO THE FINAL IMAGE-----------------
		for (int i=extint[0];i <= extint[1];i++){
			for (int j=extint[2];j<=extint[3];j++){
				for (int k=extint[4];k<=extint[5];k++){
					if(i >= ext[0] && i <= ext[1] && j >= ext[2] && j <= ext[3] && k >= ext[4] && k <= ext[5]){
						unsigned short* ptr=(unsigned short *) cast4->GetOutput()->GetScalarPointer(i,j,k); //en *ptr esta el nivel de gris del punto
						unsigned short* NEWptr=(unsigned short *)imgbin->GetScalarPointer(i,j,k);		
						if(*NEWptr==0){
							*NEWptr = *ptr;					
						}						 
					}
				}
			}
		}		
	}
	bbSetOutputImageOut( imgbin );  
	bbSetOutputFirstPoint(firstpoint);
    //----SETTING THE OUTPUT-----------------
	
  
}
void BuildBinaryImage::bbUserSetDefaultValues()
{
	imgbin =NULL;

}
void BuildBinaryImage::bbUserInitializeProcessing()
{

}
void BuildBinaryImage::bbUserFinalizeProcessing()
{

}

}
// EO namespace bbUtilities