#include "ContourExtractData.h"


//----------------------------------------------------------------------
ContourExtractData::ContourExtractData( bool okImagesResults)
{
	this->imagedata			= NULL;
	imagedataValueResult	= NULL;
	imagedataMaskResult		= NULL;
	this->okImagesResults	= okImagesResults;
	_typeOperation			= 0;
}

// ------------------------------------------------------------------------

ContourExtractData::~ContourExtractData()
{
}


//----------------------------------------------------------------------
void ContourExtractData::SetImage( vtkImageData* imagedata)
{
	this->imagedata			= imagedata;
	this->imagedata->GetScalarRange(scalarRange);

	// RaC 20-11-09 Changes in InitLstContoursLinesYPoints
	int ext[6];
	this->imagedata->GetWholeExtent(ext);
	_sizeImageY = ext[3]-ext[2]+1;

	// init vtk image result : valuesImage maskImage  
	if (this->okImagesResults==true){  InitVtkImagesResult(); }
}
//----------------------------------------------------------------------
void ContourExtractData::SetZtoBeAnalys( int z )
{
	this->zImage			= z;
}

//------------------------------------------------------------------------
void ContourExtractData::SetLstManualContourModel( std::vector<manualBaseModel*> lstManConMod)
{
	this->lstManConMod = lstManConMod;
}


//------------------------------------------------------------------------
void ContourExtractData::GetMinMaxPoint(int *minPoint, 
										int *maxPoint, 
										manualBaseModel *manualcontourmodel
										)
{
	int i;
	//int	np		= manualcontourmodel->GetSizeLstPoints( );  // number of control points // JPRx

	// JSTG 26-02-08 ---------------------------------------------------------------------------------------
	//int nps = manualviewbaseecontour->GetNumberOfPointsSpline(); // number of points in the spline
	int nps = manualcontourmodel->GetNumberOfPointsSpline(); // number of points in the spline
	//------------------------------------------------------------------------------------------------------

	// JSTG 26-02-08 ---------------------------------------------------------------------------------------
	//double x,y,z,t;
	double x,y,z;
	//double delta=( double ) ( np  ) / ( double ) ( nps  );
	manualcontourmodel->UpdateSpline();
	for (i=0; i<nps; i++)
	{
		//t= delta * (double)i;
		//manualcontourmodel->GetSplinePoint(t,x,y,z);
		manualcontourmodel->GetSpline_i_Point(i,&x,&y,&z);
		if (x<minPoint[0]){ minPoint[0]=(int)x; }
		if (y<minPoint[1]){ minPoint[1]=(int)y; }
		if (x>maxPoint[0]){ maxPoint[0]=(int)x; }
		if (y>maxPoint[1]){ maxPoint[1]=(int)y; }
	}
	minPoint[0]--;
	minPoint[1]--;
	maxPoint[0]++;
	maxPoint[1]++;
	//------------------------------------------------------------------------------------------------------
}

//------------------------------------------------------------------------
void ContourExtractData::GetMinMaxPoint_Of_LstManConMod(	int *minPoint, 
														int *maxPoint
														)

{
	int i,size = lstManConMod.size();

	for(i=0 ; i<size ; i++)
	{
		GetMinMaxPoint(minPoint,maxPoint,lstManConMod[i]);
	}
}


//------------------------------------------------------------------------
int ContourExtractData::AnalisisContourInsideV2(int x, int y, int iContour )
{
	bool inBorder=false;
	int result	= 0;
	int i;

	//cout<<"RaC ContourExtractData::AnalisisContourInsideV2-iContour:"<<iContour<<" y:"<<y<<endl; 

	int nps=_lstlstlstVecX1[iContour][y].size();

	double x1,y1,x2,y2;
	double borderX, borderY;
	double xx1, yy1,xx2, yy2;
	double xx=x, yy=y;
	double d;

	for (i=0; i<nps; i++)
	{
		x1=_lstlstlstVecX1[iContour][y][i];
		y1=_lstlstlstVecY1[iContour][y][i];
		x2=_lstlstlstVecX2[iContour][y][i];
		y2=_lstlstlstVecY2[iContour][y][i];

		borderX=x1;
		borderY=y1;

		if (y1<y2) 
		{ 
			xx1=x1; yy1=y1; xx2=x2; yy2=y2;
		} else {
			xx1=x2; yy1=y2; xx2=x1; yy2=y1; 
		} 

		double difxx2xx1=fabs(xx2-xx1);
		if (difxx2xx1==0) difxx2xx1=0.0000000001;

		// Finding border looking in vertical direction  AND  verifing if pixel is at right of the line 
		if  ( (yy>=yy1)&&(yy<=yy2) ) 
		{
			//by triangle similarity
			d = ( fabs(xx2-xx1)*(yy-yy1) ) / (yy2-yy1) ;
			if ( (xx1<=xx2)&&(x<(xx1+d)) )	{	result++;	}
			if ( (xx1>xx2)&&(x<(xx1-d)) )	{	result++;	} 

			if ( (yy2-yy1)/difxx2xx1 >= 1.0)
			{
				if (xx1<=xx2)   
				{ 
					borderX = xx1+d;
					borderY = y;
				} else { 
					borderX = xx1-d;
					borderY = y;
				} 
			} 
		} // if point inside y 


		// Finding border looking in vertical direction
		if ( ((xx1<=xx2)&&(xx>=xx1)&&(xx<xx2)) || ((xx1>xx2)&&(xx>=xx2)&&(xx<xx1)) )
		{
			if ( (yy2-yy1)/difxx2xx1 <= 1.0)
			{
				//by triangle similarity
				d = ( fabs(xx1-xx)*(yy2-yy1) ) / difxx2xx1;
				if (yy1+d<=yy2) {
					borderX=x;
					borderY=yy1+d;
				}
			} 			
		} // if point inside x			


		//Border verification
		if (   (x==(int)borderX) && (y==(int)borderY)   )  { inBorder=true; }// if point in border

		// Verification : border in horizontal line 
		if ( ((int)y1==(int)y2)  &&  ((int)y1==y) && (x1<x2) && (x>=x1) && (x<=x2))		{ inBorder=true;	}
		if ( ((int)y1==(int)y2)  &&  ((int)y1==y) && (x2<x1) && (x>=x2) && (x<=x1))		{ inBorder=true;	}

		if (inBorder==true){ i=nps; }		
	} // for i

	if (inBorder==true)	{ result=1;	}

	return result;
}



//------------------------------------------------------------------------
// typeOperation=0 AND
// typeOperation=1 OR
// typeOperation=2 XOR
bool ContourExtractData::isInside(int x, int y, int typeOperation)
{
	bool result				= false;
	int numberLeft			= 0;
	int i,size				= this->lstManConMod.size();
	int numberInside        = 0;

	/* RaC 20-11-09 (C1) Changes to use the method without the image.
	int ext[6];
	imagedata->GetExtent(ext);

	if ((x>=0) && (x<=ext[1]) && (y>=0) && (y<=ext[3]))
	{
	*/

	if (typeOperation==0)  // AND  Intersection
	{
		for (i=0;i<size;i++)
		{
			// To process the statistics of the Points contour the procedure is different
			// RaC 19-09-09
			manualBaseModel *mbm = lstManConMod[i];
			if(mbm->GetTypeModel()==7)
			{
				if(mbm->IsPoint(x,y)==true)
				{
					numberLeft=1;
				}
			}//if
			else
			{
				numberLeft =  AnalisisContourInsideV2(x,y, i );
			}//else

			if ( (numberLeft % 2) ==1){         numberInside++;  }
		}
		if ( numberInside == (size) ){ result=true; }
	} // AND  Intersection

	numberLeft=0;

	if (typeOperation==1)  // OR  All
	{
		for (i=0;i<size;i++)
		{
			// To process the statistics of the Points contour the procedure is different
			// RaC 19-09-09
			manualBaseModel *mbm = lstManConMod[i];
			if(mbm->GetTypeModel()==7)
			{
				if(mbm->IsPoint(x,y)==true)
				{
					numberLeft=1;
				}
			}//if
			else
			{
				numberLeft =  AnalisisContourInsideV2(x,y, i );
			}//else
			if ( (numberLeft % 2) ==1){ result=true;  }
		}
	} // OR  All

	numberLeft=0;

	if (typeOperation==2)  // XOR  crown
	{
		for (i=0;i<size;i++)
		{
			// To process the statistics of the Points contour the procedure is different
			// RaC 19-09-09
			manualBaseModel *mbm = lstManConMod[i];
			if(mbm->GetTypeModel()==7)
			{
				if(mbm->IsPoint(x,y)==true)
				{
					numberLeft=1;
				}
			}//if
			else
			{
				numberLeft =  numberLeft + AnalisisContourInsideV2(x,y, i );
			}//else

		}
		if ( numberLeft % 2 ==1){ 	result = true;  } 
	}// XOR  crown



	//RaC }

	return result;
}

//------------------------------------------------------------------------

double ContourExtractData::GetDataValue(int x, int y, int z)
{
	//	wxVtk2DBaseView *wxvtk2dbaseview = (wxVtk2DBaseView*)wxvtkbaseview;
	//	int z = (int)wxvtk2dbaseview->GetVtkBaseData()->GetZ();
	//JSTG 13-03-08-----
	//EED OJO avec JS	_zz = z;
	//------------------
	double result;
	void *p;
	p = imagedata->GetScalarPointer(x,y,z);

	if (imagedata->GetScalarType()==VTK_CHAR)
	{
		char *pp = (char*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_SIGNED_CHAR)
	{
		signed char *pp = (signed char*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_UNSIGNED_CHAR)
	{
		unsigned char *pp = (unsigned char*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_SHORT)
	{
		short *pp = (short*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_UNSIGNED_SHORT)
	{
		unsigned short *pp = (unsigned short*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_INT)
	{
		int *pp = (int*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_UNSIGNED_INT)
	{
		unsigned int *pp = (unsigned int*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_LONG)
	{
		long *pp = (long*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_UNSIGNED_LONG)
	{
		unsigned long *pp = (unsigned long*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_FLOAT)
	{	
		float *pp = (float*)p;
		result = (double)(*pp);
	}
	else if (imagedata->GetScalarType()==VTK_DOUBLE)
	{
		double *pp = (double*)p;
		result = (double)(*pp);
	}

	return result;
}

//------------------------------------------------------------------------

void ContourExtractData::PutVtkImageDataResultValue( int x, int y, int z, double value )
{
	unsigned short *pValue;
	unsigned short *pMask;
	pValue	= (unsigned short *)imagedataValueResult->GetScalarPointer(x,y,z);
	pMask	= (unsigned short *)imagedataMaskResult->GetScalarPointer(x,y,z);
	*pMask	= 255;
	*pValue = (unsigned short)value;
}

//------------------------------------------------------------------------
void ContourExtractData::ResetImageResult(int z)
{
	if (okImagesResults==true)
	{
		unsigned short *pValue;
		unsigned short *pMask;
		pValue	= (unsigned short *)imagedataValueResult->GetScalarPointer(0,0,z);
		pMask	= (unsigned short *)imagedataMaskResult->GetScalarPointer(0,0,z);

		int ext[6];
		imagedataValueResult->GetExtent(ext);

		int size = (ext[1]-ext[0]+1) * (ext[3]-ext[2]+1); 
		memset(pValue,0,size*2);
		memset(pMask,0,size*2);
	} // if
}


//------------------------------------------------------------------------
void ContourExtractData::CalculateImageResult()
{
	if (okImagesResults==true)
	{
		ResetImageResult(zImage);

		int minPoint[2];
		int maxPoint[2];
		int i,j;
		double value;

		minPoint[0] = 999999;
		minPoint[1] = 999999;
		maxPoint[0] = -999999;
		maxPoint[1] = -999999;

		GetMinMaxPoint_Of_LstManConMod(minPoint,maxPoint);
		InitLstContoursLinesYPoints();

		for (j=minPoint[1]; j<=maxPoint[1]; j++)
		{
			for (i=minPoint[0]; i<=maxPoint[0]; i++)
			{

				//RaC 20-11-09 Changes specified in isInside (C1)
				int ext[6];
				imagedata->GetExtent(ext);

				if ((i>=0) && (i<=ext[1]) && (j>=0) && (j<=ext[3]))
				{
                  if (isInside(i,j,_typeOperation)==true)
					{
  					  value = GetDataValue(i,j,zImage);
					  if ( (value>=scalarRange[0]) && (value<=scalarRange[1]) )
					  {
						PutVtkImageDataResultValue(i,j,zImage,  value );
					  } // scalarRange
					} // if isInside

				} // if ext
			} // for i
		} // for j


		imagedataValueResult->Modified();
		imagedataMaskResult->Modified();
		imagedataValueResult->Update();
		imagedataMaskResult->Update();
	} // if

}

//------------------------------------------------------------------------
void ContourExtractData::GetValuesInsideCrown(	int *numberOfPixels, 
												std::vector<double> *pLstValue,
												std::vector<double> *pLstValuePosX,
												std::vector<double> *pLstValuePosY,
												std::vector<double> *pLstValuePosZ)
{
	pLstValue->clear();
	pLstValuePosX->clear();
	pLstValuePosY->clear();
	pLstValuePosZ->clear();

	//	if (okImagesResults==true)
	//	{
	//		ResetImageResult(zImage);
	//	}

	int minPoint[2];
	int maxPoint[2];
	int i,j;
	double value;
	int acum=0;

	minPoint[0] = 999999;
	minPoint[1] = 999999;
	maxPoint[0] = -999999;
	maxPoint[1] = -999999;

	GetMinMaxPoint_Of_LstManConMod(minPoint,maxPoint);
	InitLstContoursLinesYPoints();

	for (j=minPoint[1]; j<=maxPoint[1]; j++)
	{
		for (i=minPoint[0]; i<=maxPoint[0]; i++)
		{
			//RaC 20-11-09 Changes specified in isInside (C1)
			int ext[6];
			imagedata->GetExtent(ext);

			if ((i>=0) && (i<=ext[1]) && (j>=0) && (j<=ext[3]))
			{
	            if (isInside(i,j,_typeOperation)==true)
				{

			  	   	acum++;
				   	value = GetDataValue(i,j,zImage);
					if ( (value>=scalarRange[0]) && (value<=scalarRange[1]) )
					{
						pLstValue		-> push_back( value );
						pLstValuePosX	-> push_back( i );
						pLstValuePosY	-> push_back( j );
						pLstValuePosZ	-> push_back( -1 );
					} // scalarRange
				} // if isInside
			} // if ext
		} // for
	} // for

	*numberOfPixels = acum;
}

//------------------------------------------------------------------------

vtkImageData *ContourExtractData::GetVtkImageValueResult()
{
	return imagedataValueResult;
}
//------------------------------------------------------------------------
vtkImageData *ContourExtractData::GetVtkImageMaskResult()
{
	return imagedataMaskResult;
}
// ------------------------------------------------------------------------
void ContourExtractData::InitVtkImagesResult()
{
	int ext[6];
	int newDim[3];
	double spc[3];
	int scalartype;

	imagedata->GetSpacing(spc);
	imagedata->GetExtent(ext);
	newDim[0]=ext[1]-ext[0]+1;
	newDim[1]=ext[3]-ext[2]+1;
	newDim[2]=ext[5]-ext[4]+1;
	scalartype = imagedata->GetScalarType();

	if (imagedataValueResult!=NULL)
	{
		imagedataValueResult->Delete();
	}
	imagedataValueResult = vtkImageData::New();
	//	  imagedataValueResult->SetScalarType(scalartype);
	imagedataValueResult->SetScalarTypeToUnsignedShort();
	imagedataValueResult->SetSpacing(spc);
	imagedataValueResult->SetDimensions( newDim );
	imagedataValueResult->AllocateScalars();

	if (imagedataMaskResult!=NULL)
	{
		imagedataMaskResult->Delete();
	}
	imagedataMaskResult  = vtkImageData::New();
	//	  imagedataMaskResult->SetScalarType(scalartype);
	imagedataMaskResult->SetScalarTypeToUnsignedShort();
	imagedataMaskResult->SetSpacing(spc);
	imagedataMaskResult->SetDimensions( newDim );
	imagedataMaskResult->AllocateScalars();
}


//------------------------------------------------------------------------
void ContourExtractData::InitVolumeStatistics()
{
	vol_rCountRange				= 0;
	vol_rsize					= 0;
	vol_minValue				= 9999999;
	vol_maxValue				=-9999999;
	vol_acum_average			= 0;
	vol_acum_standardeviation	= 0;
}

//------------------------------------------------------------------------
void ContourExtractData::SetVolumeStatistics(int rCountRange, 
											 int rsize,
											 double minValue,
											 double maxValue,
											 double acum_average,
											 double acum_standardeviation)
{
	vol_rCountRange				= vol_rCountRange + rCountRange; 
	vol_rsize					= vol_rsize + rsize; 

	if (minValue<vol_minValue){ vol_minValue = minValue;  }
	if (maxValue>vol_maxValue){ vol_maxValue = maxValue;  }

	vol_acum_average			= vol_acum_average + acum_average; 
	vol_acum_standardeviation	= vol_acum_standardeviation + acum_standardeviation; 
}

//------------------------------------------------------------------------
void ContourExtractData::GetVolumeStatistics(int *vol_rCountRange, 
											 int *vol_rsize,
											 double *vol_minValue,
											 double *vol_maxValue,
											 double *vol_average,
											 double *vol_standardeviation)
{
	*vol_rCountRange		= this->vol_rCountRange; 
	*vol_rsize				= this->vol_rsize; 
	*vol_minValue			= this->vol_minValue; 
	*vol_maxValue			= this->vol_maxValue; 
	*vol_average			= this->vol_acum_average / this->vol_rsize;
	*vol_standardeviation	= sqrt(this->vol_acum_standardeviation / this->vol_rsize);
}


//------------------------------------------------------------------------
void ContourExtractData::Statistics( std::vector<double> *inputLstValue, 
									int     grayRangeMin,
									int     grayRangeMax,
									int		*rCountRange, 
									int		*rsize, 
									double	*rmin, 
									double	*rmax,
									double	*raverage,
									double	*rstandardeviation
									)
{
	double min						= 0;
	double max						= 0;
	double average					= 0;
	double standardeviation			= 0;
	double acum_average				= 0;
	double acum_standardeviation		= 0;
	int	 size						= 0;
	int	 countRange					= 0;
	double ng;

	if (inputLstValue!=NULL)
	{
		size=inputLstValue->size();
		if (size>0){
			max=(*inputLstValue)[0];
			min=(*inputLstValue)[0];
			// Average , countRange
			int i;
			for ( i=0; i<size; i++ )
			{
				ng=(*inputLstValue)[i];
				acum_average = acum_average + ng;
				if (max<ng) max=ng;		// Max
				if (min>ng) min=ng;     // Min
				if ((ng>=grayRangeMin) && (ng<=grayRangeMax)) countRange++;  // countRange
			} // for average
			average = acum_average / size;

			// Standar Deviation
			acum_standardeviation=0;
			double tmp;
			for ( i=0; i<size; i++ )
			{
				tmp = (*inputLstValue)[i] - average;
				acum_standardeviation = acum_standardeviation + tmp*tmp;
			} // for standar deviation
			standardeviation = sqrt(acum_standardeviation/size);
			SetVolumeStatistics(countRange, (*rsize),
								min,max,
							    acum_average,acum_standardeviation);
		} // if size
	} // if NULL

	// OUTPUT
	*rsize				= size; 
	*rCountRange		= countRange;
	*rmin				= min; 
	*rmax				= max;
	*raverage			= average;
	*rstandardeviation	= standardeviation;
}

//------------------------------------------------------------------------
void ContourExtractData::SetTypeOperation(int type)
{
	_typeOperation=type;
}

//------------------------------------------------------------------------
void ContourExtractData::Fill_lstlstlstVecXY(int iContour, int sizeY)
{
	int i,y;
	double x1,y1,z1,x2,y2,z2;
	manualBaseModel *manualcontourmodel= lstManConMod[iContour];	
	int nps = manualcontourmodel->GetNumberOfPointsSpline(); // number of points in the spline
	manualcontourmodel->UpdateSpline();
	//------------------------------------------------------------------------------------------------------

	for (y=0;y<sizeY;y++)
	{
		manualcontourmodel->GetSpline_i_Point(0,&x1,&y1,&z1);
		x1=x1+0.5; y1=y1+0.5;
		for (i=1; i<nps; i++)
		{
			manualcontourmodel->GetSpline_i_Point(i,&x2,&y2,&z2);
			x2=x2+0.5; y2=y2+0.5;
			if ( ((y1<y2)&&(y>=y1)&&(y<=y2)) || ((y1>y2)&&(y>=y2)&&(y<=y1)) || ((int)y1==y) || ((int)y2==y)  )
			{
				_lstlstlstVecX1[iContour][y].push_back(x1);
				_lstlstlstVecY1[iContour][y].push_back(y1);
				_lstlstlstVecX2[iContour][y].push_back(x2);
				_lstlstlstVecY2[iContour][y].push_back(y2);
			} 
			x1=x2; y1=y2; z1=z2;
		} // for i  Points in spline
	} // y

}

void ContourExtractData::InitLstContoursLinesYPoints()
{
	// init InInside Optimisation  
	int i;

	_lstlstlstVecX1.clear();
	_lstlstlstVecY1.clear();
	_lstlstlstVecX2.clear();
	_lstlstlstVecY2.clear();

	/* RaC 20-11-09
	int ext[6];
	this->imagedata->GetWholeExtent(ext);
	int sizeY = ext[3]-ext[2]+1;
	*/
	std::vector<double> vecDouble;
	std::vector< std::vector<double> > vecVecDouble;
	for ( i=0 ; i<_sizeImageY ; i++ )
	{
		vecVecDouble.push_back( vecDouble );
	}

	//Fill structure with points
	int sizeContours = lstManConMod.size();
	for( i=0 ; i<sizeContours ; i++ )
	{
		_lstlstlstVecX1.push_back( vecVecDouble );
		_lstlstlstVecY1.push_back( vecVecDouble );
		_lstlstlstVecX2.push_back( vecVecDouble );
		_lstlstlstVecY2.push_back( vecVecDouble );
		Fill_lstlstlstVecXY(i,_sizeImageY);
	}

}

void ContourExtractData::SetScalarRange(double min, double max)
{
	scalarRange[0]=min;
	scalarRange[1]=max;
}

void ContourExtractData::SetSizeImageY(int pSizeImageY)
{
	_sizeImageY=pSizeImageY;
}