/*# ---------------------------------------------------------------------
#
# Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
#                        pour la Sant�)
# Authors : Eduardo Davila, Frederic Cervenansky, Claire Mouton
# Previous Authors : Laurent Guigues, Jean-Pierre Roux
# CreaTools website : www.creatis.insa-lyon.fr/site/fr/creatools_accueil
#
#  This software is governed by the CeCILL-B license under French law and
#  abiding by the rules of distribution of free software. You can  use,
#  modify and/ or redistribute the software under the terms of the CeCILL-B
#  license as circulated by CEA, CNRS and INRIA at the following URL
#  http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
#  or in the file LICENSE.txt.
#
#  As a counterpart to the access to the source code and  rights to copy,
#  modify and redistribute granted by the license, users are provided only
#  with a limited warranty  and the software's author,  the holder of the
#  economic rights,  and the successive licensors  have only  limited
#  liability.
#
#  The fact that you are presently reading this means that you have had
#  knowledge of the CeCILL-B license and that you accept its terms.
# ------------------------------------------------------------------------ */

#include "Propagation.h"

//------------------------------------------------------------------------------------------------------------------------------------------
//CLASS: Vector ------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
//Constructor
Vector::Vector()
{
	_plane = -1;
	_var = -1;
}
//Destructor
Vector::~Vector()
{
	_vec.clear();
	_vecX.clear();
	_vecY.clear();
	_vecZ.clear();
}

//------------------------------------------------------------------------------------------------------------------------------------------
void Vector::set_vec(double val)
{
	_vec.push_back(val);
}
//------------------------------------------------------------------------------------------------------------------------------------------
void Vector::set_var(double val)
{
	_var = val;
}
//------------------------------------------------------------------------------------------------------------------------------------------
double Vector::get_vec(int id)
{
	if(_vec.size() != 0)
	{
		return _vec[id];
	}
	return -1;
}
//------------------------------------------------------------------------------------------------------------------------------------------
double Vector::get_var()
{
	return _var;
}
//------------------------------------------------------------------------------------------------------------------------------------------
int Vector::getsize_vec()
{
	return _vec.size();
}
//------------------------------------------------------------------------------------------------------------------------------------------
void Vector::copyVector( std::vector<Vector>*vec1, std::vector<Vector>*vec2 )
{
	int i,j;
	if( vec1->size() != 0 )
	{
		vec2->clear();
		for(i=0; i<(int)(vec1->size()); i++)
		{
			Vector *temp = new Vector();
			temp->set_var( (*vec1)[i].get_var() );
			for(j=0; j<(*vec1)[i].getsize_vec(); j++)
			{
				temp->set_vec( (*vec1)[i].get_vec(j) );
			}
			vec2->push_back(*temp);
			delete temp;
		}
	}
}
//------------------------------------------------------------------------------------------------------------------------------------------
void Vector::printVector(std::vector<Vector>*vec1)
{
	int i,j;
	for(i=0; i<(int)(vec1->size()); i++)
	{
		printf("\n Pos (%d) => var = %f",i,(*vec1)[i].get_var());
		for(j=0; j<(*vec1)[i].getsize_vec(); j++)
		{
			printf("\n vec(%d) = %f",j,(*vec1)[i].get_vec(j));
		}
	}
}
//------------------------------------------------------------------------------------------------------------------------------------------
void Vector::set_plane(int val)
{
	_plane = val;
}
//------------------------------------------------------------------------------------------------------------------------------------------
void Vector::set_x(double val)
{
	_vecX.push_back(val);
}
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
void Vector::set_y(double val)
{
	_vecY.push_back(val);
}
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
void Vector::set_z(double val)
{
	_vecZ.push_back(val);
}
//------------------------------------------------------------------------------------------------------------------------------------------
int Vector::get_plane()
{
	return _plane;
}
//------------------------------------------------------------------------------------------------------------------------------------------
double Vector::get_x(int id)
{
	if( (-1<id) && (id<(int)(_vecX.size())  ) )
	{
		return _vecX[id];
	}
	return -1;
}
//------------------------------------------------------------------------------------------------------------------------------------------
int Vector::getsize_x()
{
	if(_vecX.size() != 0)
	{
		return _vecX.size();
	}
	return -1;
}
//------------------------------------------------------------------------------------------------------------------------------------------
double Vector::get_y(int id)
{
	if( (-1<id) && (id<(int)(_vecY.size())) )
	{
		return _vecY[id];
	}
	return -1;
}
//------------------------------------------------------------------------------------------------------------------------------------------
int Vector::getsize_y()
{
	if(_vecY.size() != 0)
	{
		return _vecY.size();
	}
	return -1;
}
//------------------------------------------------------------------------------------------------------------------------------------------
int Vector::getsize_z()
{
	if(_vecZ.size() != 0)
	{
		return _vecZ.size();
	}
	return -1;
}
//------------------------------------------------------------------------------------------------------------------------------------------
double Vector::get_z(int id)
{
	if( (-1<id) && (id<(int)(_vecZ.size())) )
	{
		return _vecZ[id];
	}
	return -1;
}
//------------------------------------------------------------------------------------------------------------------------------------------
std::vector<double> Vector::getVec()
{
	return _vec;
}
//------------------------------------------------------------------------------------------------------------------------------------------
void Vector::resetVec()
{
	_vec.clear();
}
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------


//Constructor
PropContour::PropContour()
{
	_interpnumber = 100;
}

//Destructor
PropContour::~PropContour()
{
	ResetPlaneVector();
	ResetKeyContours();
}

/* EED 03/07/2008
//------------------------------------------------------------------------------------
double PropContour::RBF_WendLand(double norm, double m_rad)
{
	double y;
	  norm = norm / m_rad;
	  y = pow( 1-norm,4 ) * ( (4*norm) + 1 );
	  if(norm >= 1)
	  {
		y = 0;
	  }
	return y; 
}
//------------------------------------------------------------------------------------
double PropContour::RBF_ThinPlate(double norm)
{
	double y;
	if(norm == 0)
	{
		y = 0;
	}
	else
	{
		y = pow(norm,2)*log(norm);
	}
	return y;
}
//------------------------------------------------------------------------------------
vtkImageData* PropContour::method_RBF ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY, 
										std::vector<double>*CoordZ )
{
	_dimImage[0] = 100;			// X Axis
	_dimImage[1] = 100;			// Y Axis
	_dimImage[2] = 1;			// Z Axis
	int pointsSize = CoordX->size();
	double spc[3]={1,1,1};
	double norm = 0, val = 0;

	vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
	vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
	vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
	vnl_vector<double> H( pointsSize, 1 );
	vnl_vector<double> D( pointsSize, 0.0 );
	vnl_vector<double> Q( 2, 0.0 );
	vnl_matrix<double> Impl( _dimImage[0],_dimImage[1], 0.0 );

	unsigned short *pValue;
	imagedataValue = vtkImageData::New();
	imagedataValue->SetScalarTypeToUnsignedShort();
	imagedataValue->SetSpacing(spc);
	imagedataValue->SetDimensions(_dimImage);
	imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
	imagedataValue->AllocateScalars();
	imagedataValue->Update();

    int i,j,h;
	for(i=0; i<pointsSize; i++)
	{
		for(j=0; j<pointsSize; j++)
		{
			norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) );
			A(i,j) = RBF_WendLand(norm,rad);
			//A(i,j) = RBF_ThinPlate(norm);
		}
	}

	Ainv = vnl_matrix_inverse<double>(A);
	D = Ainv* H;

	for(i=0; i < _dimImage[1]; i++)
	{
		for(j=0; j < _dimImage[0]; j++)
		{
			Q(0) = i;
			Q(1) = j;
			val = 0;
			for(h=0; h<pointsSize; h++)
			{
				norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) );
				//val = val + ( D(h) * RBF_ThinPlate(norm) );
				val = val + ( D(h) * RBF_WendLand(norm,rad) );
			}
			Impl(i,j) = val - 1.0;

//			if ( (Impl(i,j)>=-0.001) && (Impl(i,j) <=0.001)  )
//			{Impl(i,j)=128;
//			}
//			else
//			{Impl(i,j)=0;}


//			pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0);
//			vvalue = (Impl(i,j)*256+256);
//			if (vvalue) < 0) 
//			{ 
//				*pValue = 0;
//			} else {
//				*pValue = vvalue;
//			}

			pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0);
			if ( Impl(i,j) >=0 ) 
			{ 
				*pValue = 128;
			} 
			else 
			{
				*pValue = 0;
			}
		}
	}

	for(i=0; i<pointsSize; i=i+5)
	{
         pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,0);
		 *pValue=255;
	}

	return imagedataValue;
}
//----------------------------------------------------------------------------------------------------
double PropContour::RBF_ThinPlate_3D(double norm)
{
	return norm = pow( norm,3 );
}
//----------------------------------------------------------------------------------------------------
vtkImageData* PropContour::method_RBF_3D (  double rad, std::vector<double>*CoordX, std::vector<double>*CoordY, 
											std::vector<double>*CoordZ )
{
	long interval = wxGetElapsedTime(TRUE);

	int i,j,k,h;
	int pointsSize = CoordX->size();
	double max = -1, min = 1000, minz = 100000, maxz = -10000;

	for( i=0; i<pointsSize; i++ )
	{
		if( (*CoordX)[i] > max )
		{
			max = (*CoordX)[i]; 
		}
		if( (*CoordY)[i] > max )
		{
			max = (*CoordY)[i]; 
		}
		if( (*CoordX)[i] < min )
		{
			min = (*CoordX)[i]; 
		}
		if( (*CoordY)[i] < min )
		{
			min = (*CoordY)[i]; 
		}
		if( (*CoordZ)[i] < minz )
		{
			minz = (*CoordZ)[i];
		}
		if( (*CoordZ)[i] > maxz )
		{
			maxz = (*CoordZ)[i];
		}
	}

	_dimImage[0] = 200;			// X axis
	_dimImage[1] = 200;			// Y axis
	_dimImage[2] = 200;			// Z axis

	double spc[3]={1,1,1};
	double norm = 0, val, vvalue;

	vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
	vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
	vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
	vnl_vector<double> H( pointsSize, 1 );
	vnl_vector<double> D( pointsSize, 0.0 );
	vnl_vector<double> Q( 3, 0.0 );

	unsigned short *pValue;
	imagedataValue = vtkImageData::New();
	imagedataValue->SetScalarTypeToUnsignedShort();
	imagedataValue->SetSpacing(spc);
	imagedataValue->SetDimensions(_dimImage);
	imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
	imagedataValue->AllocateScalars();
	imagedataValue->Update();

	for(i=0; i<pointsSize; i++)
	{
		for(j=0; j<pointsSize; j++)
		{
			norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) );
			A(i,j) = RBF_WendLand(norm,rad);
		}
	}
	D = vnl_matrix_inverse<double>(A)* H;

//Inicialization
//	for(k=0; k<_dimImage[2]; k++)
//	{
//		for(j=0; j<_dimImage[1]; j++)
//		{
//			for(i=0; i<_dimImage[0]; i++)
//			{
//				pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
//				*pValue = 0;
//			}
//		}
//	}

//Filling
//	pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0);
//	i = 0;
//	j = 0;
//	k = 0;
//	for(h=0; h<pointsSize; h++)
//	{
//		norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
//		//val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
//		val = val + ( D(h) * RBF_WendLand(norm,rad) );
//		if( h == pointsSize-1 )
//		{
//			val = val - 1;
//			pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
//			if ( val >=0 ) 
//			{ 
//				*pValue = 128;
//			} 
//			else 
//			{
//				*pValue = 0;
//			}
//			//pValue++;
//			j++;
//			h = -1;
//			val = 0;
//			if( j == _dimImage[0] )
//			{
//				i++;
//				j = 0;
//				h = -1;
//				if( i == _dimImage[1])
//				{
//					k++;
//					i = 0;
//					j = 0;
//					h = -1;
//					if(k == _dimImage[2])
//					{
//						h = pointsSize;
//					}
//				}
//			}
//		}
//	}

	i = (int)min-10;
	j = (int)min-10;
	k = (int)minz-10;
	val = 0;

	for(h=0; h<pointsSize; h++)
	{
		norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
		val = val + ( D(h) * RBF_WendLand(norm,rad) );
		if( h == pointsSize-1 )
		{
			val = val - 1;
			pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
			if ( val >=0 ) 
			{ 
				*pValue = 128;
			} 
			else 
			{
				*pValue = 0;
			}
			//pValue++;
			j++;
			h = -1;
			val = 0;
			if( j == (int)max+10)
			{
				i++;
				j = min-10;
				h = -1;
				if( i == (int)max+10)
				{
					k++;
					i = min-10;
					j = min-10;
					h = -1;
					if(k == maxz+10)
					{
						h = pointsSize;
					}
				}
			}
		}
	}

	interval = wxGetElapsedTime(FALSE);
	long interPlane = interval/_dimImage[2];

	for(i=0; i<pointsSize; i++)
	{
		 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]);
		*pValue=255;
	}

	long intervalPC = wxGetElapsedTime();

	printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------");
	printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval);
	printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
	printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane);
	printf("\n NUMBER OF PLANES................................. %d",k);
	printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
	printf("\n ------------------------------------------------------------");

	return imagedataValue;
}
//----------------------------------------------------------------------------------------------------

vtkImageData* PropContour::method_RBF_3D_ThinPlate (  double rad, std::vector<double>*CoordX, std::vector<double>*CoordY, 
													  std::vector<double>*CoordZ )
{
	long interval = wxGetElapsedTime(TRUE);

	int i,j,k,h;
	int pointsSize = CoordX->size();
	double max = -1, min = 1000, minz = 100000, maxz = -10000;

	for( i=0; i<pointsSize; i++ )
	{
		if( (*CoordX)[i] > max )
		{
			max = (*CoordX)[i]; 
		}
		if( (*CoordY)[i] > max )
		{
			max = (*CoordY)[i]; 
		}
		if( (*CoordX)[i] < min )
		{
			min = (*CoordX)[i]; 
		}
		if( (*CoordY)[i] < min )
		{
			min = (*CoordY)[i]; 
		}
		if( (*CoordZ)[i] < minz )
		{
			minz = (*CoordZ)[i];
		}
		if( (*CoordZ)[i] > maxz )
		{
			maxz = (*CoordZ)[i];
		}
	}

	_dimImage[0] = 190;			// X axis
	_dimImage[1] = 190;			// Y axis
	_dimImage[2] = 190;			// Z axis

	double spc[3]={1,1,1};
	double norm = 0, val, vvalue;

	vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
	vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
	vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
	vnl_vector<double> H( pointsSize, 1 );
	vnl_vector<double> D( pointsSize, 0.0 );
	vnl_vector<double> Q( 3, 0.0 );

	unsigned short *pValue;
	imagedataValue = vtkImageData::New();
	imagedataValue->SetScalarTypeToUnsignedShort();
	imagedataValue->SetSpacing(spc);
	imagedataValue->SetDimensions(_dimImage);
	imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,-10,_dimImage[2]-1);
	imagedataValue->AllocateScalars();
	imagedataValue->Update();

	for(i=0; i<pointsSize; i++)
	{
		for(j=0; j<pointsSize; j++)
		{
			norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) );
			A(i,j) = RBF_ThinPlate_3D(norm);
		}
	}
	D = vnl_matrix_inverse<double>(A)* H;

//Inicialization
//	for(k=0; k<_dimImage[2]; k++)
//	{
//		for(j=0; j<_dimImage[1]; j++)
//		{
//			for(i=0; i<_dimImage[0]; i++)
//			{
//				pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
//				*pValue = 0;
//			}
//		}
//	}

//Filling

//	pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0);
//	i = 0;
//	j = 0;
//	k = 0;
//	for(h=0; h<pointsSize; h++)
//	{
//		norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
//		//val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
//		val = val + ( D(h) * RBF_WendLand(norm,rad) );
//		if( h == pointsSize-1 )
//		{
//			val = val - 1;
//			pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
//			if ( val >=0 ) 
//			{ 
//				*pValue = 128;
//			} 
//			else 
//			{
//				*pValue = 0;
//			}
//			//pValue++;
//			j++;
//			h = -1;
//			val = 0;
//			if( j == _dimImage[0] )
//			{
//				i++;
//				j = 0;
//				h = -1;
//				if( i == _dimImage[1])
//				{
//					k++;
//					i = 0;
//					j = 0;
//					h = -1;
//					if(k == _dimImage[2])
//					{
//						h = pointsSize;
//					}
//				}
//			}
//		}
//	}

	i = (int)min-10;
	j = (int)min-10;
	k = (int)minz-10;
	val = 0;

	for(h=0; h<pointsSize; h++)
	{
		norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
		val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
		if( h == pointsSize-1 )
		{
			val = val - 1;
			pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
			if ( val >=0 ) 
			{ 
				*pValue = 128;
			} 
			else 
			{
				*pValue = 0;
			}
			//pValue++;
			j++;
			h = -1;
			val = 0;
			if( j == (int)max+10)
			{
				i++;
				j = min-10;
				h = -1;
				if( i == (int)max+10)
				{
					k++;
					i = min-10;
					j = min-10;
					h = -1;
					if(k == maxz+10)
					{
						h = pointsSize;
					}
				}
			}
		}
	}

	interval = wxGetElapsedTime(FALSE);
	long interPlane = interval/_dimImage[2];

	for(i=0; i<pointsSize; i++)
	{
		 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]);
		*pValue=255;
	}

	long intervalPC = wxGetElapsedTime();

	printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------");
	printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval);
	printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
	printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane);
	printf("\n NUMBER OF PLANES................................. %d",k);
	printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
	printf("\n ------------------------------------------------------------");

	return imagedataValue;
}

*/

//---------------------------------------------------------------------------------------------------------
void PropContour::ReadKeyContour(FILE* fd)
{
	char firstline[30];
	int	size;
	double x,y;		
	int z;
	std::vector<double> tempX;
	std::vector<double> tempY;
	std::vector<double> tempZ;
	tempX.clear();
	tempY.clear();
	tempZ.clear();
	while(!feof(fd))
	{
		//fscanf(fd," %s %d",&firstline,&size); // JPRx
		fscanf(fd," %s %d",firstline,&size);
		for(int i=0; i<size; i++)
		{
			fscanf(fd,"%lf %lf %d",&x,&y,&z);
			tempX.push_back(x);
			tempY.push_back(y);
			tempZ.push_back(z);
		}
		SetKeyContours(&tempX,&tempY,&tempZ);
		tempX.clear();
		tempY.clear();
		tempZ.clear();
	}
}
//---------------------------------------------------------------------------------------------------------
int PropContour::VectorDirection(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
{
	int dir,i;
	double SumX = 0,SumY = 0;
	double ax,ay,bx,by,axb;
	int size = InX->size();
	for(i=0; i<size; i++)
	{
		SumX = SumX + (*InX)[i];
		SumY = SumY + (*InY)[i];
	}
	SumX = SumX/size;		//Mass Center: X coord
	SumY = SumY/size;		//Mass Center: Y coord
	
	int positive = 0;
	int negative = 0;
	for(i=0; i<size; i++)
	{
		ax = (*InX)[i]-SumX;
		ay = (*InY)[i]-SumY;
		bx = (*InX)[(i+1)%size]-SumX;
		by = (*InY)[(i+1)%size]-SumY;
		axb = (ax*by) - (bx*ay);
		if(axb > 0)
		{
			positive++;
		}
		if(axb < 0)
		{
			negative++;
		}
	}
	if(positive >= negative)
	{
		dir = 1;
	}
	else
	{
		dir = -1;
	}

	return dir;
}
//----------------------------------------------------------------------------------------------------
void PropContour::VectorOrder(int dir, int posini, std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
{
	int size = InX->size();
	int i;
	std::vector<double> tempX;
	std::vector<double> tempY;
	std::vector<double> tempZ;
	
	tempX.clear();
	tempY.clear();
	tempZ.clear();

	for(i=0; i<size; i++)
	{
		if(dir == 1)
		{
			tempX.push_back((*InX)[posini]);
			tempY.push_back((*InY)[posini]);
			tempZ.push_back((*InZ)[posini]);
			posini++;
			if(posini == size)
			{
				posini = 0;
			}
		}
		if(dir == -1)
		{
			tempX.push_back((*InX)[posini]);
			tempY.push_back((*InY)[posini]);
			tempZ.push_back((*InZ)[posini]);
			posini--;
			if(posini < 0)
			{
				posini = size-1;
			}
		}
	}
	InX->clear();
	InY->clear();
	InZ->clear();
	for(i=0; i<size; i++)
	{
		InX->push_back(tempX[i]);
		InY->push_back(tempY[i]);
		InZ->push_back(tempZ[i]);
	}

}
//----------------------------------------------------------------------------------------------------
void PropContour::PreparePointsForSpline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ, 
										 std::vector<int>*Sizes)
{
	int sizeS = Sizes->size();
	//int sizeV = InX->size(); // JPRx
	int i,j,mem,posinic,dir,cont;
	double leX;

	std::vector<double> tempX;
	std::vector<double> tempY;
	std::vector<double> tempZ;
	std::vector<double> lstX;
	std::vector<double> lstY;
	std::vector<double> lstZ;
	
	lstX.clear();
	lstY.clear();
	lstZ.clear();

	mem = 0;
	cont = 0;
	for(i=0; i<sizeS; i++)
	{
		leX=1000;
		tempX.clear();
		tempY.clear();
		tempZ.clear();
		for(j=0; j<(*Sizes)[i]; j++)
		{
			tempX.push_back((*InX)[j+mem]);
			tempY.push_back((*InY)[j+mem]);
			tempZ.push_back((*InZ)[j+mem]);
			if( (*InX)[j] < leX )
			{
				posinic = j;
				leX = (*InX)[j];
			}
		}
		mem = mem + (*Sizes)[i];
		dir = VectorDirection(&tempX,&tempY,&tempZ);
		VectorOrder(dir,posinic,&tempX,&tempY,&tempZ);

		for(j=0; j<(*Sizes)[i]; j++)
		{
			lstX.push_back(tempX[j]);
			lstY.push_back(tempY[j]);
			lstZ.push_back((*InZ)[j+cont]);
		}
		cont = cont + (*Sizes)[i];
	}

//Fill the Finally lst in X,Y,Z ---------------
	int sizetemp = lstX.size();
	//printf("\nJSTG-PropContour::PreparePointsForSpline");
	InX->clear();
	InY->clear();
	InZ->clear();
	for(i=0; i<sizetemp; i++)
	{
		//printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
		InX->push_back(lstX[i]);
		InY->push_back(lstY[i]);
		InZ->push_back(lstZ[i]);
	}

	
}
//----------------------------------------------------------------------------------------------------
vtkImageData* PropContour::method_Spline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ, std::vector<int>*Sizes)
{
	//long interval = wxGetElapsedTime(TRUE); // JPRx

	int i,j,k,sizeX,sizeS,sizeInS;
	int numspline;
	double x,y,z;
//EED	double spc[3]={1,1,1};
	std::vector<double> lstX;
	std::vector<double> lstY;
	std::vector<double> lstZ;
	std::vector<double> tempX;
	std::vector<double> tempY;
	std::vector<double> tempZ;
	std::vector<double> tempS;
	_mContourModel = new manualContourModel();
	_mContourModel->SetNumberOfPointsSpline(_interpnumber);

	imagedataValue=NULL;
//EED
//	_dimImage[0] = 200;			// X axis
//	_dimImage[1] = 200;			// Y axis
//	_dimImage[2] = 200;			// Z axis
//
//	unsigned short *pValue;
//	imagedataValue = vtkImageData::New();
//	imagedataValue->SetScalarTypeToUnsignedShort();
//	imagedataValue->SetSpacing(spc);
//	imagedataValue->SetDimensions(_dimImage);
//	imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
//	imagedataValue->AllocateScalars();
//	imagedataValue->Update();
	
//	lstX.clear();
//	lstY.clear();
//	lstZ.clear();

	sizeX = InX->size();
	for(i=0; i<sizeX; i++)
	{
		lstX.push_back((*InX)[i]);
		lstY.push_back((*InY)[i]);
		lstZ.push_back((*InZ)[i]);
	}
	
	PreparePointsForSpline(&lstX,&lstY,&lstZ,Sizes);
	/*int sizetemp = lstX.size();
	printf("\nJSTG-PropContour::method_Spline");
	for(i=0; i<sizetemp; i++)
	{
		printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
	}*/

	sizeS = Sizes->size();
	int cont = 0;
	for(i=0; i<sizeS; i++)
	{
		_mContourModel->DeleteAllPoints();
		_mContourModel->SetCloseContour(true);
		sizeInS = (*Sizes)[i];
		for(j=0; j<sizeInS; j++)
		{
			_mContourModel->AddPoint(lstX[cont],lstY[cont],lstZ[cont]);
			cont ++;
		}
		_mContourModel->UpdateSpline();
		numspline = _mContourModel->GetNumberOfPointsSpline();
		for(k=0; k<numspline; k++)
		{
			_mContourModel->GetSpline_i_Point(k,&x,&y,&z);
//EED			pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
//			*pValue = 128;
			tempX.push_back(x);
			tempY.push_back(y);
			tempZ.push_back(z);
		}
	}
	
	int tam = numspline;
	std::vector<Vector>::iterator it;
	
	ResetPlaneVector();

	for(i=0; i<numspline; i++)
	{
		Vector *vec = new Vector();
		_planevector.push_back(*vec);
	}// for i

	for(j=0; j<tam; j++)
	{
		_mContourModel->DeleteAllPoints();
		_mContourModel->SetCloseContour(false);
		cont = 0;
		for(i=0; i<sizeS; i++)
		{
			//double hh=tempZ[cont+j]; // JPRx
			_mContourModel->AddPoint(tempX[cont+j],tempY[cont+j],tempZ[cont+j]);
			cont = cont + tam;
		} // for i
		_mContourModel->UpdateSpline();
		numspline = _mContourModel->GetNumberOfPointsSpline();
		for(k=0; k<numspline; k++)
		{
//EED002
			_mContourModel->GetSpline_i_Point(k,&x,&y,&z);
			_planevector[k].set_x(x);
			_planevector[k].set_y(y);
			_planevector[k].set_z(z);
			_planevector[k].set_plane(k);
//EED			pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
//EED			*pValue = 128;
		} // for k
	} // for j

/*
	int tempsize = _planevector.size();
	int sizexx;
	for(i=0; i<tempsize; i++)
	{
		sizexx = _planevector[i].getsize_x();
	}
*/
/*
	interval = wxGetElapsedTime(FALSE);
	int zplanes = abs(tempZ[0]-tempZ[tempZ.size()-1]);
	long interPlane = interval/zplanes;
	int pointsSize = 0;
	for(i=0; i<sizeS; i++)
	{
		pointsSize = pointsSize + (*Sizes)[i]; 
	}
	for(i=0; i<pointsSize; i++)
	{
		 pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*InX)[i] , (int)(*InY)[i] ,(int)(*InZ)[i]);
		*pValue=255;
	}

	printf("\n\n JSTG - PropContour::method_Spline ------------------------");
	printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %d (ms)",interval);
	//printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
	printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %d (ms)",interPlane);
	printf("\n NUMBER OF PLANES................................. %d",zplanes);
	printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
	printf("\n ------------------------------------------------------------");
*/
	return imagedataValue;
}

//---------------------------------------------------------------------------------------------------
void PropContour::ResetPlaneVector()
{
//	Vector *vec;
//	int ii,iiSize = _planevector.size();
//	for (ii=0 ; ii<iiSize ; ii++) 
//	{
//		vec = &(_planevector[ii]);
//		delete vec;
//	}
	_planevector.clear();
}


//---------------------------------------------------------------------------------------------------
void PropContour::ResetKeyContours()
{
	_KeyContourSizes.clear();
	_KeyContourX.clear();
	_KeyContourY.clear();
	_KeyContourZ.clear();
}
//---------------------------------------------------------------------------------------------------
void PropContour::SetKeyContours(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
{

	int idKeyContour		= 0;
	int idKeyContourSizes	= 0;
	int tmpIdKeyContSizes	= 0;
	bool okFind				= false;
	int i;
	int sizeKeyContour,Z=(int)(*InZ)[0];
	sizeKeyContour = _KeyContourZ.size();
	for (i=0; i<sizeKeyContour; i++)
	{ 
		if (i>0)
		{
			if ( (_KeyContourZ[i-1]<Z) && (_KeyContourZ[i]>=Z) )
			{
				idKeyContour		= i;
				idKeyContourSizes	= tmpIdKeyContSizes;
				okFind=true;
				i=sizeKeyContour;
			}
			if ( (i<sizeKeyContour) && (_KeyContourZ[i-1] != _KeyContourZ[i]) )
			{
				tmpIdKeyContSizes++;
			}
		} else {
			if  (_KeyContourZ[0]>Z) 
			{
				idKeyContour		= 0;
				idKeyContourSizes	= 0;
				okFind				= true;
				i					= sizeKeyContour;
			}
		} // if >O
	} // for 

	if (okFind==false)
	{
		idKeyContour		= _KeyContourX.size();
		idKeyContourSizes	= _KeyContourSizes.size();
		okFind=true;
	}

	_KeyContourSizes.insert( _KeyContourSizes.begin() + idKeyContourSizes , InX->size() );
	for(i=0; i<(int)(InX->size()); i++)
	{
		_KeyContourX.insert( _KeyContourX.begin() + idKeyContour, (*InX)[i] );
		_KeyContourY.insert( _KeyContourY.begin() + idKeyContour, (*InY)[i] );
		_KeyContourZ.insert( _KeyContourZ.begin() + idKeyContour, (*InZ)[i] );
	}


//EED
//	_KeyContourSizes.push_back( InX->size() );
//	for(i=0; i<InX->size(); i++)
//	{
//		_KeyContourX.push_back( (*InX)[i] );
//		_KeyContourY.push_back( (*InY)[i] );
//		_KeyContourZ.push_back( (*InZ)[i] );
//	}

}
//---------------------------------------------------------------------------------------------------
vtkImageData* PropContour::CalculeSplinePropagation()
{
	if(_KeyContourSizes.size() <= 0)
	{
		printf("\n There would be at last 1 contour");
		return NULL;
	}
	if(_KeyContourSizes.size() == 1)
	{
		return NULL;
	}
	if(_KeyContourSizes.size() >= 2)
	{
		return method_Spline(&_KeyContourX,&_KeyContourY,&_KeyContourZ,&_KeyContourSizes);
	}
	return NULL;
}
//---------------------------------------------------------------------------------------------------
void PropContour::GetKeyContours(std::vector<double>*KeyX, std::vector<double>*KeyY, std::vector<double>*KeyZ, std::vector<int>*KeyS)
{
	int i;
	KeyX->clear();
	KeyY->clear();
	KeyZ->clear();
	KeyS->clear();

	for(i=0; i<(int)(_KeyContourSizes.size()); i++)
	{
		KeyS->push_back( _KeyContourSizes[i] );
	}
	for(i=0; i<(int)(_KeyContourX.size()); i++)
	{
		KeyX->push_back( _KeyContourX[i] );
		KeyY->push_back( _KeyContourY[i] );
		KeyZ->push_back( _KeyContourZ[i] );
	}
}
//---------------------------------------------------------------------------------------------------
void PropContour::GetPropagatedContours( std::vector<Vector>*planevec )
{
	int i,j;
	planevec->clear();
	for(i=0; i<(int)(_planevector.size()); i++)
	{
		Vector *temp = new Vector();
		temp->set_plane( _planevector[i].get_plane() );
		for(j=0; j<_planevector[i].getsize_x(); j++)
		{
			temp->set_x( _planevector[i].get_x(j) );
			temp->set_y( _planevector[i].get_y(j) );
			temp->set_z( _planevector[i].get_z(j) );
		}
		planevec->push_back(*temp);
		delete temp;
	}
}
//---------------------------------------------------------------------------------------------------
void PropContour::SetInterpNumber(int val)
{
	_interpnumber = val;
}

int  PropContour::FindIdWithZ(double z)
{
	int result_ID=0;
	int k,size=_planevector.size();
	double dist,minDist=99999999;
	for ( k=0 ; k<size ; k++)
	{
		dist=fabs( z - _planevector[k].get_z(0) );
		if (dist<minDist)
		{
			minDist=dist;
			result_ID = k;
		}
	}// for i
	return result_ID;
}

//---------------------------------------------------------------------------------------------------
void PropContour::GetIdContour(int id, std::vector<double>*vecX, std::vector<double>*vecY, std::vector<double>*vecZ)
{
	int i;
	vecX->clear();
	vecY->clear();
	vecZ->clear();
	//int sizeplane = _planevector[id].getsize_x();
	double tempx;
	for(i=0; i<_planevector[id].getsize_x(); i++)
	{
		vecX->push_back( _planevector[id].get_x(i) );
		tempx = _planevector[id].get_x(i);
		vecY->push_back( _planevector[id].get_y(i) );
		vecZ->push_back( _planevector[id].get_z(i) );
	}
}
//---------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------
//---------------------------------------------------------------------------------------------------