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

BBTK_ADD_BLACK_BOX_TO_PACKAGE(Utilities,PolyDataToImageData)
BBTK_BLACK_BOX_IMPLEMENTATION(PolyDataToImageData,bbtk::AtomicBlackBox);
void PolyDataToImageData::Process()
{
 

	vtkPolyData* poly = bbGetInputIn();
	//vtkImageData* img = bbGetInputImage();
	//_img = NULL;
	double* bounds = NULL;
	if(poly!=NULL){//&&img!=NULL){

		bounds = poly->GetBounds();
		//img = createImageData(poly);
		//img = createImageData1(poly);
		_img = createImageData2(poly);
		//img = createImageData3(poly, img);
	}

	bbSetOutputBounds(bounds);
	bbSetOutputOut(_img);

      
}
vtkImageData* PolyDataToImageData::createImageData3(vtkPolyData* poly, vtkImageData* img){
	
	/*vtkImageData* img = vtkImageData::New();
	img->SetScalarTypeToUnsignedChar();
	img->SetOrigin( 0,0,0 );
	img->SetSpacing(0.26483,0.26483,0.700012);
	img->SetDimensions(189,171,73);
	img->AllocateScalars();*/
	/*int iDim[3];
	img->GetDimensions(iDim);
	
	for(int n=0;n<iDim[0]*iDim[1]*iDim[2];n++)
	{
		img->GetPointData()->GetScalars()->SetComponent(n,0,0);
	}*/

	//int iDim[3];
	//img->GetDimensions(iDim);
	//img->SetDimensions(iDim[0]*2,iDim[1]*2,iDim[2]*2 );
	//img->AllocateScalars();

	vtkPolyDataToImageStencil* imgStencil = vtkPolyDataToImageStencil::New();
	imgStencil->SetInput( poly );
	//imgStencil->SetInput( normalsFilter->GetOutput() );
	//imgStencil->SetInformationInput(img);
	//imgStencil->Update();
 
	
    vtkImageStencil * stencil = vtkImageStencil::New();	
	stencil->SetInput( img );			
	stencil->SetStencil( imgStencil->GetOutput() );
	stencil->ReverseStencilOff();
	stencil->SetBackgroundValue( 0 );
	stencil->Update();

	

	vtkImageData* img1 = stencil->GetOutput();
		
	//return getImageBorder(img1);
	return img1;

}
vtkImageData* PolyDataToImageData::getImageBorder(vtkImageData* img){

	vtkImageCast* cast4 = vtkImageCast::New();
	cast4->SetInput(img);
	cast4->SetOutputScalarTypeToUnsignedShort();
	cast4->Update();


    vtkImageData* img1 = cast4->GetOutput();

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

	for(int x = ext[0]; x <= ext[1];x++){
		for(int y = ext[2]; y <= ext[3]  ; y++){
			for(int z = ext[4]; z <= ext[5] ; z++){	
				unsigned short* ptr=(unsigned short*) img1->GetScalarPointer(x,y,z);
				if(*ptr!=0){
					*ptr = 50;
				}
			}
		}
	}
	
	for(int x = ext[0]; x <= ext[1];x++){
		for(int y = ext[2]; y <= ext[3]  ; y++){
			for(int z = ext[4]; z <= ext[5] ; z++){	
                bool color = false;
				unsigned short* ptr=(unsigned short*) img1->GetScalarPointer(x,y,z);

				if(*ptr == 50){
					unsigned short* ptr1;
					for(int i = x - 1; i <= x + 1; i++){
						for(int j = y - 1; j <= y + 1; j++){
							for(int k = z - 1; k <= z + 1; k++){
								if(ext[0]<= i && i <= ext[1]&&ext[2]<= j && j <= ext[3]&&ext[4]<= k && k <= ext[5]){
									ptr1=(unsigned short*) img1->GetScalarPointer(i,j,k);
									if(*ptr1==0){
										color = true;
									}
								}						
							}
						}						
					}
				}
				if(color){
					*ptr= 100;	
				}
			}
		}
	}
	return img1;
}
vtkImageData* PolyDataToImageData::createImageData2(vtkPolyData* poly){

	double dBounds[6];
	int iDim[3];
	double dSpacing[3];
//	double dPointLocation[3];
	int iPointLocation[3];

	double point0[3];
	double point1[3];
	double point2[3];
	double *vect0;
	double *vect1;
	double *vectCross;
	double temppoint[3], lamb, pointproy[3];
	
//	double distNorm, distPoint, d;
	double  distPoint, d;
	
	std::vector<double*> vectpoints;

	int minx,maxx, miny, maxy, minz, maxz;

	poly->GetBounds(dBounds);

	dSpacing[0]=0.1;
	dSpacing[1]=0.1;
	dSpacing[2]=0.1;

	iDim[0]=1+(int)((dBounds[1]-dBounds[0])/dSpacing[0]);
	iDim[1]=1+(int)((dBounds[3]-dBounds[2])/dSpacing[1]);
	iDim[2]=1+(int)((dBounds[5]-dBounds[4])/dSpacing[2]);
	
	vtkImageData* img = vtkImageData::New();
	img->SetScalarTypeToUnsignedShort();
	img->SetNumberOfScalarComponents(1);
	img->SetSpacing(1,1,1);		
	img->SetDimensions(iDim);
	img->AllocateScalars();

	int n, tmpMax=iDim[0]*iDim[1]*iDim[2];
	
	for( n=0 ; n<tmpMax ; n++ )
	{
		img->GetPointData()->GetScalars()->SetComponent(n,0,0);
	}

	n = poly->GetNumberOfCells();
	
	for(int i = 0; i < n; i++)
	{			
		vtkPoints* points = poly->GetCell(i)->GetPoints();			

		points->GetPoint(0,point0);
		points->GetPoint(1,point1);
		points->GetPoint(2,point2);			

		point0[0] = (point0[0]-dBounds[0])/dSpacing[0];
		point0[1] = (point0[1]-dBounds[2])/dSpacing[1];
		point0[2] = (point0[2]-dBounds[4])/dSpacing[2];
		//roundPoint(point0);

		point1[0] = (point1[0]-dBounds[0])/dSpacing[0];
		point1[1] = (point1[1]-dBounds[2])/dSpacing[1];
		point1[2] = (point1[2]-dBounds[4])/dSpacing[2];
		//roundPoint(point1);

		point2[0] = (point2[0]-dBounds[0])/dSpacing[0];
		point2[1] = (point2[1]-dBounds[2])/dSpacing[1];
		point2[2] = (point2[2]-dBounds[4])/dSpacing[2];
		//roundPoint(point2);

		vectpoints.clear();
		vectpoints.push_back(point0);
		vectpoints.push_back(point1);
		vectpoints.push_back(point2);

		int size = vectpoints.size();
		if(size>3){
			size++;
		}

		vect0 = makeVector(point0, point1);
		vect1 = makeVector(point0, point2);
		vectCross = getCrossProduct(vect0, vect1);
		//vectCross[0]*x + vectCross[1]*y + vectCross[2]*z + d = 0
		d = -1*getPointProduct(vectCross, point0);

		/*minx = (int)((minval(minval(point0[0], point1[0]), point2[0])-dBounds[0])/dSpacing[0]);
		miny = (int)((minval(minval(point0[1], point1[1]), point2[1])-dBounds[2])/dSpacing[1]);
		minz = (int)((minval(minval(point0[2], point1[2]), point2[2])-dBounds[4])/dSpacing[2]);
		maxx = (int)((maxval(maxval(point0[0], point1[0]), point2[0])-dBounds[0])/dSpacing[0]);
		maxy = (int)((maxval(maxval(point0[1], point1[1]), point2[1])-dBounds[2])/dSpacing[1]);
		maxz = (int)((maxval(maxval(point0[2], point1[2]), point2[2])-dBounds[4])/dSpacing[2]);*/

		minx = minval(minval(point0[0], point1[0]), point2[0]);
		miny = minval(minval(point0[1], point1[1]), point2[1]);
		minz = minval(minval(point0[2], point1[2]), point2[2]);
		maxx = maxval(maxval(point0[0], point1[0]), point2[0]);
		maxy = maxval(maxval(point0[1], point1[1]), point2[1]);
		maxz = maxval(maxval(point0[2], point1[2]), point2[2]);

		        
		for(int x = minx; x <= maxx;x++){
			for(int y = miny; y <= maxy  ; y++){
				for(int z = minz; z <= maxz ; z++){				             
					
					temppoint[0] = x;
					temppoint[1] = y;
					temppoint[2] = z;
					
					lamb = getLambda(vectCross, temppoint, d);

					distPoint = abs(lamb)*getMagnitud(vectCross);

					pointproy[0] = lamb*vectCross[0]+temppoint[0];
					pointproy[1] = lamb*vectCross[1]+temppoint[1];
					pointproy[2] = lamb*vectCross[2]+temppoint[2];

					
					if(distPoint < 1) {
						//if( pointInsideAngles(pointproy, vectpoints)){
						if(pointInsideDistance(pointproy, point0, point1)||
							pointInsideDistance(pointproy, point1, point2)||
							pointInsideDistance(pointproy, point2, point0)){

							iPointLocation[0] = x;
							iPointLocation[1] = y;
							iPointLocation[2] = z;
							colorPoint(iPointLocation,img,100);

						}else{

							iPointLocation[0] = x;
							iPointLocation[1] = y;
							iPointLocation[2] = z;
							//colorPoint(iPointLocation,img,200);

						}
					}				                    
				}				
			}
		}	

		delete vect0;
		delete vect1;
		delete vectCross;

	}	
	return img;

}
bool PolyDataToImageData::pointInsideDistance(double* pointproy,double point0[3], double point1[3]){

	bool ret = false;
	double* vect0 = makeVector(point0, point1);
	double* vect1 = makeVector(point0, pointproy);
	double* v0xv1 = getCrossProduct(vect0, vect1);
	double* vectcross = getCrossProduct(vect0, v0xv1);

	//vectCross[0]*x + vectCross[1]*y + vectCross[2]*z + d = 0
	double d = -1*getPointProduct(vectcross, point0);	
	    				
	double lamb = getLambda(vectcross, pointproy, d);

	double pointproytemp[3];
	pointproytemp[0] = lamb*vectcross[0]+pointproy[0];
	pointproytemp[1] = lamb*vectcross[1]+pointproy[1];
	pointproytemp[2] = lamb*vectcross[2]+pointproy[2];

    double dist = getDistance(point0, point1);
	double dist1 = getDistance(pointproytemp, point0);
	double dist2 = getDistance(pointproytemp, point1);

	if(dist1<=dist && dist2 <= dist){		
		ret = true;
	}
	delete vect0;
	delete vect1;
	delete v0xv1;
	delete vectcross;
	return ret;
}
/**
**	Returns if the point pointproy is inside the plane described by point0, point1, point2 forming
**	a triangle
**/
bool PolyDataToImageData::pointInsideAngles(double* pointproy,std::vector<double*> points){

	double a, b, c;
	double alfa, beta, gama;

	double a1, b1, c1;
	double alfa1, beta1, gama1;

	int index1, index2;

	std::ofstream file1;
		

	if(points.size()>2){
		
		for(int i = 0; i < (int)(points.size()); i++){

			if(i == (int)(points.size())-1){
				index1 = 0;
			}else{
				index1 = i + 1;
			}

			if(i == 0){
				index2 = points.size()-1;
				
			}else{
				index2 = i - 1;
			}

			a = getDistance(points[index1], points[index2]);
			b = getDistance(points[index2], points[i] );
			c = getDistance(points[i], points[index1]);
			
			alfa = getAngleTheoCos(b, c, a);
			beta = getAngleTheoCos(c, a, b);
			gama = getAngleTheoCos(a, b, c);

			if(alfa+beta+gama > 3.1416 || alfa+beta+gama < 3.1415 ){
				std::cout<<"Error "<<alfa+beta+gama<<std::endl;
			}

						
			a1 = getDistance(points[index1], pointproy);
			b1 = getDistance(pointproy, points[i]  );
			c1 = getDistance(points[i], points[index1]);

			alfa1 = getAngleTheoCos(b1, c1, a1);
			beta1 = getAngleTheoCos(c1, a1, b1);
			gama1 = getAngleTheoCos(a1, b1, c1);	

			if(alfa1+beta1+gama1 > 3.1416 || alfa1+beta1+gama1 < 3.1415 ){
				std::cout<<"Error "<<alfa1+beta1+gama1<<std::endl;
			}


			if(alfa1 > alfa){
				//std::cout<<"Out alfa= "<<alfa<<" alfa1= "<<alfa1<<std::endl;
				/*file1.open( "D:/creaNefsev/Desarrollo/Test/test.txt", ios::out | ios::app| ios::ate );
				if(file1.is_open())
				{
					file1 << "<Group>"<<std::endl;
					file1 << "<Triangle>"<<std::endl;
					file1 << "\t<Points Triangle>"<<std::endl;
					file1 << "\t\t<pointA>"<<"<x>"<<points[i][0]<<"</x><y>"<<points[i][1]<<"</y><z>"<<points[i][2]<<"</z></pointA>"<<std::endl;
					file1 << "\t\t<pointB>"<<"<x>"<<points[index1][0]<<"</x><y>"<<points[index1][1]<<"</y><z>"<<points[index1][2]<<"</z></pointB>"<<std::endl;
					file1 << "\t\t<pointC>"<<"<x>"<<points[index2][0]<<"</x><y>"<<points[index2][1]<<"</y><z>"<<points[index2][2]<<"</z></pointC>"<<std::endl;
					file1 << "\t</Points Triangle>"<<std::endl;

					file1 << "\t<Angles Triangle>"<<std::endl;
					file1 << "\t\t<alfa>"<<alfa<<"</alfa>"<<std::endl;
					file1 << "\t\t<beta>"<<beta<<"</beta>"<<std::endl;
					file1 << "\t\t<gama>"<<gama<<"</gama>"<<std::endl;
					file1 << "\t</Angles Triangle>"<<std::endl;

					file1 << "\t<Sides Triangle>"<<std::endl;
					file1 << "\t\t<lenghta>"<<a<<"</lenghta>"<<std::endl;
					file1 << "\t\t<lenghtb>"<<b<<"</lenghtb>"<<std::endl;
					file1 << "\t\t<lenghtc>"<<c<<"</lenghtc>"<<std::endl;
					file1 << "\t</Sides Triangle>"<<std::endl;
					file1 << "</Triangle>"<<std::endl;

					file1 << "<Triangle>"<<std::endl;
					file1 << "\t<Points Triangle>"<<std::endl;
					file1 << "\t\t<pointA>"<<"<x>"<<points[i][0]<<"</x><y>"<<points[i][1]<<"</y><z>"<<points[i][2]<<"</z></pointA>"<<std::endl;
					file1 << "\t\t<pointB>"<<"<x>"<<points[index1][0]<<"</x><y>"<<points[index1][1]<<"</y><z>"<<points[index1][2]<<"</z></pointB>"<<std::endl;
					file1 << "\t\t<pointC>"<<"<x>"<<pointproy[0]<<"</x><y>"<<pointproy[1]<<"</y><z>"<<pointproy[2]<<"</z></pointC>"<<std::endl;
					file1 << "\t</Points Triangle>"<<std::endl;

					file1 << "\t<Angles Triangle>"<<std::endl;
					file1 << "\t\t<alfa>"<<alfa1<<"</alfa>"<<std::endl;
					file1 << "\t\t<beta>"<<beta1<<"</beta>"<<std::endl;
					file1 << "\t\t<gama>"<<gama1<<"</gama>"<<std::endl;
					file1 << "\t</Angles Triangle>"<<std::endl;

					file1 << "\t<Sides Triangle>"<<std::endl;
					file1 << "\t\t<lenghta>"<<a1<<"</lenghta>"<<std::endl;
					file1 << "\t\t<lenghtb>"<<b1<<"</lenghtb>"<<std::endl;
					file1 << "\t\t<lenghtc>"<<c1<<"</lenghtc>"<<std::endl;
					file1 << "\t</Sides Triangle>"<<std::endl;
					file1 << "</Triangle>"<<std::endl;

					file1 << "</Group>"<<std::endl;
											
					
				}
				file1.close();*/
				return false;
			}
		}
	}	
	return true;

}
/**
**	Get's the angle between the 3 distances using the Cos theorem
**	c^2 = a^2 + b^2 - 2*a*b cos g
**
**/
double PolyDataToImageData::getAngleTheoCos(double a, double b, double c){

	double val;

	val = pow(c, 2) - (pow(a, 2) + pow(b, 2));
	val = -(val/(2*a*b));

    return 	acos(val);
    
}
/**
**	Returns the value of the distance between the 2 given points
**/
double PolyDataToImageData::getDistance(double* point0, double* point1){
	double dist = sqrt(pow(point0[0]-point1[0],2) + pow(point0[1]-point1[1],2) + pow(point0[2]-point1[2],2));
	return dist;
}
/**
**	Returns the magnitud of the given vector
**/
double PolyDataToImageData::getMagnitud(double* vect){

	double mag;

	mag = sqrt(pow(vect[0],2) + pow(vect[1],2) + pow(vect[2],2));

	return mag;

}
/**
**	Return the escalar product between the two vectors [x0, y0, z0] . [x1, y1, z1]
**	[x0, y0, z0] . [x1, y1, z1] = x0*x1 + y0*y1 + z0*z1 
**/
double PolyDataToImageData::getPointProduct(double* vect0,double* vect1){
	double val;

	val = vect0[0]*vect1[0] + vect0[1]*vect1[1] + vect0[2]*vect1[2];

	return val;
}
/**
**	returns the unitary vector of the given vector
**	u = 1/|vect| . vect
**/
double* PolyDataToImageData::getNormal(double* vect){

	double* vectnorm;
	double mag = getMagnitud(vect);

	vectnorm = new double[3];
	

	if(mag!=0){
		vectnorm[0] = vect[0]/mag;
		vectnorm[1] = vect[1]/mag;
		vectnorm[2] = vect[2]/mag;
	}else{
		vectnorm[0] = 0;
		vectnorm[1] = 0;
		vectnorm[2] = 0;
	}

	return vectnorm;


}
/**
**	Returns lambda 
**	lambda = -(a*x1 + b*y1 + c*z1 + d)/(a^2 + b^2 + c^2)
**	where [a, b, c] is a vector normal to the plane and x1, y1, z1 is the given point outside the plane
**	d is the constant of the plane in the ecuation ax + by + cz + d = 0 describing the plane
**/
double PolyDataToImageData::getLambda(double* vectnorm, double* point, double d){
	double lamb=0, xn, yn, zn, px, py,pz, mag = 0;

	xn = vectnorm[0];
	yn = vectnorm[1];
	zn = vectnorm[2];
	px = point[0];
	py = point[1];
	pz = point[2];


	mag = pow(xn, 2) + pow(yn, 2) + pow(zn, 2);

	if(mag != 0){
		lamb = -(xn*px + yn*py + zn*pz + d)/mag;
	}
	return lamb;

}
vtkImageData* PolyDataToImageData::createImageData1(vtkPolyData* poly){

	int iDim[3];
	double dBounds[6];
	double dSpacing[3];

	dSpacing[0]=0.003;
	dSpacing[1]=0.003;
	dSpacing[2]=0.003;

	poly->GetBounds(dBounds);

	iDim[0]=2+(int)((dBounds[1]-dBounds[0])/dSpacing[0]);
	iDim[1]=2+(int)((dBounds[3]-dBounds[2])/dSpacing[1]);
	iDim[2]=2+(int)((dBounds[5]-dBounds[4])/dSpacing[2]);


	vtkImageData* img = vtkImageData::New();
	img->SetScalarTypeToUnsignedShort();
	img->SetNumberOfScalarComponents(1);
	img->SetSpacing(1,1,1);		
	img->SetDimensions(iDim);
	img->AllocateScalars();


	vtkPolyDataToImageStencil* polydataToImageStencil = vtkPolyDataToImageStencil::New();
	polydataToImageStencil->SetInput( poly );
	polydataToImageStencil->Update();

	vtkImageStencil* imageStencil = vtkImageStencil::New();
	imageStencil->SetStencil( polydataToImageStencil->GetOutput());
	imageStencil->SetInput( img);
	imageStencil->Update();

	return imageStencil->GetOutput();

}
vtkImageData* PolyDataToImageData::createImageData(vtkPolyData* poly){
	
	double dBounds[6];
	int iDim[3];
	double dSpacing[3];
//	double dPointLocation[3];
	int iPointLocation[3];

	double point0[3];
	double point1[3];
	double point2[3];
	double *vect0;
	double *vect1;
	double *vectCross;

	int minx,maxx, miny, maxy, minz, maxz;

	poly->GetBounds(dBounds);

	dSpacing[0]=0.0003;
	dSpacing[1]=0.0003;
	dSpacing[2]=0.0003;

	iDim[0]=2+(int)((dBounds[1]-dBounds[0])/dSpacing[0]);
	iDim[1]=2+(int)((dBounds[3]-dBounds[2])/dSpacing[1]);
	iDim[2]=2+(int)((dBounds[5]-dBounds[4])/dSpacing[2]);
	
	vtkImageData* img = vtkImageData::New();
	img->SetScalarTypeToUnsignedShort();
	img->SetNumberOfScalarComponents(1);
	img->SetSpacing(dSpacing[0],dSpacing[1],dSpacing[2]);		
	img->SetDimensions(iDim);
	img->AllocateScalars();

	int n, tmpMax=iDim[0]*iDim[1]*iDim[2];
	for(n=0 ; n<tmpMax ; n++ )
	{
		img->GetPointData()->GetScalars()->SetComponent(n,0,0);
	}

	n = poly->GetNumberOfCells();
	
	for(int i = 0; i < n; i++)
	{			
		vtkPoints* points = poly->GetCell(i)->GetPoints();			

		points->GetPoint(0,point0);
		points->GetPoint(1,point1);
		points->GetPoint(2,point2);			

		//int iIndex=(iPointLocation[2] *iDim[0]*iDim[1]) + ( iPointLocation[1] *iDim[0]) + iPointLocation[0];
		point0[0] = (point0[0]-dBounds[0])/dSpacing[0];
		point0[1] = (point0[1]-dBounds[2])/dSpacing[1];
		point0[2] = (point0[2]-dBounds[4])/dSpacing[2];
		roundPoint(point0);

		point1[0] = (point1[0]-dBounds[0])/dSpacing[0];
		point1[1] = (point1[1]-dBounds[2])/dSpacing[1];
		point1[2] = (point1[2]-dBounds[4])/dSpacing[2];
		roundPoint(point1);

		point2[0] = (point2[0]-dBounds[0])/dSpacing[0];
		point2[1] = (point2[1]-dBounds[2])/dSpacing[1];
		point2[2] = (point2[2]-dBounds[4])/dSpacing[2];
		roundPoint(point2);

		vect0 = makeVector(point0, point1);
		vect1 = makeVector(point0, point2);
		vectCross = getCrossProduct(vect0, vect1);

		/*minx = (int)((minval(minval(point0[0], point1[0]), point2[0])-dBounds[0])/dSpacing[0]);
		miny = (int)((minval(minval(point0[1], point1[1]), point2[1])-dBounds[2])/dSpacing[1]);
		minz = (int)((minval(minval(point0[2], point1[2]), point2[2])-dBounds[4])/dSpacing[2]);
		maxx = (int)((maxval(maxval(point0[0], point1[0]), point2[0])-dBounds[0])/dSpacing[0]);
		maxy = (int)((maxval(maxval(point0[1], point1[1]), point2[1])-dBounds[2])/dSpacing[1]);
		maxz = (int)((maxval(maxval(point0[2], point1[2]), point2[2])-dBounds[4])/dSpacing[2]);*/

		minx = minval(minval(point0[0], point1[0]), point2[0]);
		miny = minval(minval(point0[1], point1[1]), point2[1]);
		minz = minval(minval(point0[2], point1[2]), point2[2]);
		maxx = maxval(maxval(point0[0], point1[0]), point2[0]);
		maxy = maxval(maxval(point0[1], point1[1]), point2[1]);
		maxz = maxval(maxval(point0[2], point1[2]), point2[2]);

		        
		for(int x = minx; x <= maxx;x++){
			for(int y = miny; y <= maxy  ; y++){
				for(int z = minz; z <= maxz ; z++){
                    double temppoint[3];
					double *tempvect;
					/*temppoint[0] = (x*dSpacing[0])+dBounds[0];
					temppoint[1] = (y*dSpacing[1])+dBounds[2];
					temppoint[2] = (z*dSpacing[2])+dBounds[4];*/

					temppoint[0] = x;
					temppoint[1] = y;
					temppoint[2] = z;

					tempvect = makeVector(point0, temppoint);			
					
					if(isOnPlane(vectCross, tempvect)){
						iPointLocation[0] = x;
						iPointLocation[1] = y;
						iPointLocation[2] = z;
						colorPoint(iPointLocation,img,100);
					}
				
					delete tempvect;
				}				
			}
		}	

		delete vect0;
		delete vect1;
		delete vectCross;

	}	
	return img;
}

void PolyDataToImageData::roundPoint(double* point){

	double fractpart, intpart;  
  
  	for(int i = 0; i < 3; i++){

		fractpart = modf (point[i] , &intpart);

		if(fractpart >= 0.5){
			point[i] = ceil(point[i]);
		}else{
			point[i] = floor(point[i]);
		}
		
	}
}

bool PolyDataToImageData::isOnPlane(double* vect0, double* vect1){
	int x;
	x =  vect0[0]*vect1[0] + (vect0[1]*vect1[1]) + vect0[2]*vect1[2];

	if(x == 0){
		return true;
	}
	return false;
}

double PolyDataToImageData::minval(double x0, double x1){
	if(x0 < x1)
		return x0;
	else
		return x1;
}
double PolyDataToImageData::maxval(double x0, double x1){
	if(x0>x1)
		return x0;
	else
		return x1;
}

double* PolyDataToImageData::makeVector(double point0[3], double point1[3]){
	double *vect;
	vect = new double[3];

	vect[0]= point0[0]-point1[0];
	vect[1]= point0[1]-point1[1];
	vect[2]= point0[2]-point1[2];

	return vect;

}

void PolyDataToImageData::colorPoint(int* point, vtkImageData* img, int grey){

	int x, y, z, ext[6];
	x = point[0];
	y = point[1];
	z = point[2];

	img->GetExtent(ext);	

	if(ext[0]<= x && x <= ext[1]&&ext[2]<= y && y <= ext[3]&&ext[4]<= z && z <= ext[5]){
		unsigned short* ptr=(unsigned short*) img->GetScalarPointer(x,y,z); 		
		*ptr= grey;	
	}	
	
}

double* PolyDataToImageData::getCrossProduct(double* vect0,double* vect1){
	double* vectCross;
	vectCross = new double[3];
	vectCross[0] = vect0[1]*vect1[2]-(vect0[2]*vect1[1]);
	vectCross[1] = -(vect0[0]*vect1[2]-(vect0[2]*vect1[0]));
	vectCross[2] = vect0[0]*vect1[1]-(vect0[1]*vect1[0]);

	return vectCross;
}

void PolyDataToImageData::bbUserSetDefaultValues()
{
	_img = NULL;
    bbSetInputIn(NULL); 
	bbSetInputImage(NULL); 

}
void PolyDataToImageData::bbUserInitializeProcessing()
{

}
void PolyDataToImageData::bbUserFinalizeProcessing()
{
	if(_img!=NULL){
		_img->Delete();
	}
}
}
// EO namespace bbUtilities