/*=========================================================================


=========================================================================*/
#include "axisExtractor02.h"
#include "vtkPolyData.h"
#include "vtkObjectFactory.h"
#include "vtkImageThreshold.h"
#include "vtkImageCast.h"
#include "vtkImageSeedConnectivity.h"
#include "vtkImageData.h"
#include "vtkMarchingCubes.h"
#include "vtkDoubleArray.h"
#include "vtkPointData.h"
#include "vtkExtractVOI.h"
#include "vtkImageConstantPad.h"
#include "vtkImageTranslateExtent.h"
#include "vtkMath.h"
#include "vtkTransformPolyDataFilter.h"
#include "vtkTransform.h"
#include <time.h>





vtkStandardNewMacro(axisExtractor02);


//----------------------------------------------------------------------------
void axisExtractor02::SetParam(double value)
{
	this->param=value;
}



//----------------------------------------------------------------------------
double axisExtractor02::GetParam()
{
	return this->param;
}




//----------------------------------------------------------------------------
void axisExtractor02::SetParam2(double value)
{
	this->param2=value;
}



//----------------------------------------------------------------------------
double axisExtractor02::GetParam2()
{
	return this->param2;
}




//----------------------------------------------------------------------------
void axisExtractor02::SetParam3(double value)
{
	this->param3=value;
}



//----------------------------------------------------------------------------
double axisExtractor02::GetParam3()
{
	return this->param3;
}



//----------------------------------------------------------------------------
void axisExtractor02::SetMaxant(int value)
{
	this->maxant=value;
}



//----------------------------------------------------------------------------
int axisExtractor02::GetMaxant()
{
	return this->maxant;
}



//----------------------------------------------------------------------------
void axisExtractor02::SetMinant(int value)
{
	this->minant=value;
}



//----------------------------------------------------------------------------
int axisExtractor02::GetMinant()
{
	return this->minant;
}


//----------------------------------------------------------------------------
void axisExtractor02::SetPoint(double value[3])
{


	double spa[3];

	this->GetInput()->GetSpacing (spa);

	value[0]=value[0]+maxant*spa[0];
	value[1]=value[1]+maxant*spa[1];
	value[2]=value[2]+maxant*spa[2];

	this->m_Stack0.push_front(value[0]);
	this->m_Stack1.push_front(value[1]);
	this->m_Stack2.push_front(value[2]);
	this->m_Stack3.push_front(value[0]);
	this->m_Stack4.push_front(value[1]);
	this->m_Stack5.push_front(value[2]);
	this->m_Stack6.push_front(value[0]);
	this->m_Stack7.push_front(value[1]);
	this->m_Stack8.push_front(value[2]);
	this->m_Stack.push_front(0);
	this->m_Stackr.push_front(0);
	this->m_Stackra.push_front(0);
	this->m_Stackrp.push_front(0);
}



//----------------------------------------------------------------------------
vtkImageData *axisExtractor02::GetVolumen()
{

	vtkImageTranslateExtent *trans;

	trans = vtkImageTranslateExtent::New();

	trans->SetInput(this->data4);
	
	trans->SetTranslation (maxant, maxant, maxant);

	trans->Update();

	return trans->GetOutput();
}





//----------------------------------------------------------------------------
vtkPolyData  *axisExtractor02::GetOutput ()
{

		

	
	vtkTransform *transL1;
	vtkTransformPolyDataFilter *trans;

	transL1 = vtkTransform::New();
	trans = vtkTransformPolyDataFilter::New();

// EED 30 Oct 2006
	double spa[3];
	this->GetInput()->GetSpacing (spa);
	transL1->Translate(-maxant*spa[0], -maxant*spa[1], -maxant*spa[2]);
//	transL1->Translate(-maxant, -maxant, -maxant);

	
	trans->SetInput((vtkPolyData *)(this->Outputs[0]));

	trans->SetTransform (transL1);

	trans->Update();

	return trans->GetOutput();
}







//----------------------------------------------------------------------------
axisExtractor02::axisExtractor02() {


	this->NumberOfRequiredInputs = 1;
	this->param=1;
	this->param2=1;
	this->param3=0.5;
	this->param4=0;
	
	/*this->resample= vtkImageResample::New();
	this->resample->SetDimensionality (3);*/

	this->data4=vtkImageData::New();	

	this->data1=vtkImageData::New();
	this->data2=vtkImageData::New();
	this->data3=vtkImageData::New();
	this->data6=vtkImageData::New();
		
	this->extrac = vtkExtractVOI::New();
//	this->extrac->SetInput(resample->GetOutput());
	this->extrac->SetSampleRate(1, 1, 1);
	
	this->connect = vtkImageSeedConnectivity::New();
	this->connect->SetInput(this->data1);
	this->connect->SetInputConnectValue(255);
	this->connect->SetOutputConnectedValue(255);
	this->connect->SetOutputUnconnectedValue(0);

	this->distance= vtkImageEuclideanDistance::New();
	this->distance->SetInput(this->connect->GetOutput());
	this->distance->InitializeOn();
	this->distance->ConsiderAnisotropyOff();
			

}



//----------------------------------------------------------------------------
void axisExtractor02::SetInput(vtkImageData *input)
{

	vtkImageConstantPad *pad;
	vtkImageTranslateExtent *trans;

	pad = vtkImageConstantPad::New();
	trans = vtkImageTranslateExtent::New();

	pad->SetInput(input);
	trans->SetInput(pad->GetOutput());
	
	
	pad->SetConstant(0);

	pad->SetInput(input);

	int ext[6];

	input->GetExtent(ext);

	ext[0]=ext[0]-maxant;
	ext[1]=ext[1]+maxant;
	ext[2]=ext[2]-maxant;
	ext[3]=ext[3]+maxant;
	ext[4]=ext[4]-maxant;
	ext[5]=ext[5]+maxant;

	pad->SetOutputWholeExtent(ext);

	trans->SetTranslation (maxant, maxant, maxant);

	trans->Update();


	//this->vtkProcessObject::SetNthInput(0, input);
	this->vtkProcessObject::SetNthInput(0, trans->GetOutput());

	//this->GetInput()->SetOrigin(0,0,0);
	//this->resample->SetInput(this->GetInput());
	this->extrac->SetInput(this->GetInput());
  
}



//----------------------------------------------------------------------------
void axisExtractor02::distanciaejes(vtkPolyData *eje1,  vtkPolyData *eje2)
{
	vtkIdType totalpuntos1=eje1->GetNumberOfPoints();
	vtkIdType totalpuntos2=eje2->GetNumberOfPoints();
	vtkIdType totallineas=eje2->GetNumberOfLines();

	vtkCell * lineas;

	int i, j;

	double point[3], point2[3], point3[3], point4[3];

	double v[3], v2[3];

	double x, dist, min, y;
	FILE *ff;

	ff=fopen("comparacion.txt","w");

	

	for(i=0;i<totalpuntos1;i++){
		eje1->GetPoint(i, point);
		min =5000;
		for(j=0;j<totallineas;j++){
			lineas=eje2->GetCell(j);
					
					lineas->GetPoints()->GetPoint(0, point3);
					lineas->GetPoints()->GetPoint(1, point2);
					v[0]=point3[0]-point2[0];
					v[1]=point3[1]-point2[1];
					v[2]=point3[2]-point2[2];
					
					x=((v[0]*point[0])+(v[1]*point[1])+(v[2]*point[2])-(v[0]*point2[0])-(v[1]*point2[1])-(v[2]*point2[2]))/((v[0]*v[0])+(v[1]*v[1])+(v[2]*v[2]));
					
					
					point4[0]=point2[0]+(x*v[0]);
					point4[1]=point2[1]+(x*v[1]);
					point4[2]=point2[2]+(x*v[2]);

					v2[0]=point4[0]-point[0];
					v2[1]=point4[1]-point[1];
					v2[2]=point4[2]-point[2];

					y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);

					if(x>=0 && x<=1){
						dist=sqrt(((point4[0]-point[0])*(point4[0]-point[0]))+((point4[1]-point[1])*(point4[1]-point[1]))+((point4[2]-point[2])*(point4[2]-point[2])));
						if(dist<min){
							min=dist;


							/*printf("punto: %f %f %f\n", point[0], point[1], point[2]);
							printf("punto2: %f %f %f\n", point2[0], point2[1], point2[2]);
							printf("punto3: %f %f %f\n", point3[0], point3[1], point3[2]);
							printf("punto4: %f %f %f\n", point4[0], point4[1], point4[2]);
							printf("v: %f %f %f\n", v[0], v[1], v[2]);
							printf("v2: %f %f %f\n", v2[0], v2[1], v2[2]);
							printf("%f\n", min);
							printf("%f\n", x);
							printf("%f\n", y);*/



						}
					}
				
			/*}*/
		}
		if(min<5){
			fprintf(ff,"%f\n", min);
		}

			

	}

	fclose(ff);

	
	
}




//----------------------------------------------------------------------------
void axisExtractor02::Execute()
{
	flagg=0;
	flagg2=0;
	frama=0;
	fseg=0;
	mejordst=0;
	mejorrad=0;
	mejorcant=0;
	
	//double minspac;

	time_t start,end;
	double dif;
//	FILE *ff;
	
	this->GetInput()->GetExtent(extprin0);
	this->GetInput()->GetExtent(extprin);
	this->GetInput()->GetSpacing(espprin);

// EED
//	stream  = fopen( "resultado.txt", "w" );
//	fprintf(stream, "accion\tindexp\tbuenos\tradiotrabajo\tradioactual\tradiopred\tradioanterior\tmejordst\tmejorrad\tmejorcant\tflag2\tflag3\tflag4\tflag5\tflag6\tflag7\tflag8\tflag9\tflag10\tflag11\tflag12\tflag13\tflag14\tflag15\tflag16\tflag17\tflag18\tflag19\tflag20\tflagg\tflagg2\tcantidad\tcantidadb\tburned\tpuntoactual\tpuntoanterior\tpuntopre\tmejor\n");

	/*minspac=min(espprin[0],min(espprin[1],espprin[2]));

	this->resample->SetAxisOutputSpacing( 0, minspac);
	this->resample->SetAxisOutputSpacing( 1, minspac);
	this->resample->SetAxisOutputSpacing( 2, minspac);
	resample->Update();*/
	

	this->data4->SetScalarTypeToUnsignedChar();
	this->data4->SetExtent(this->GetInput()->GetExtent());
	this->data4->SetSpacing(this->GetInput()->GetSpacing());

	this->blanquear(this->data4);

	this->points = vtkPoints::New();
	this->buenos=0;
	this->lineas = vtkCellArray::New();


	time (&start);
	this->todo();

// EED
//	fclose( stream );
	
	((vtkPolyData *)this->Outputs[0])->SetPoints (this->points);
	((vtkPolyData *)this->Outputs[0])->SetLines(this->lineas);
	time (&end);
	dif = difftime (end,start);

//	ff=fopen("time.txt","w");
//	fprintf(ff,"%d \n",this->buenos);
//	fprintf(ff,"%.2lf \n",dif);
//	fprintf(ff,"%.2lf \n",(double)this->buenos/dif);
//	fclose(ff);
	
}






//----------------------------------------------------------------------------
vtkImageData *axisExtractor02::GetInput()
{
	if (this->NumberOfInputs < 1){
		return NULL;
    }
  
	return (vtkImageData *)(this->Inputs[0]);
}







//----------------------------------------------------------------------------
void axisExtractor02::PrintSelf(ostream& os, vtkIndent indent)
{
	this->Superclass::PrintSelf(os,indent);
}




//----------------------------------------------------------------------------
void axisExtractor02::realtoreal(double a[3], double b[3] )
{
		
	b[0]=a[0]+(espprin[0]*extprin[0]);
	b[1]=a[1]+(espprin[1]*extprin[2]);
	b[2]=a[2];

	

}



//----------------------------------------------------------------------------
void axisExtractor02::realtoreal2(double a[3], double b[3] )
{
		
	b[0]=a[0]-(espprin[0]*extprin[0]);
	b[1]=a[1]-(espprin[1]*extprin[2]);
	b[2]=a[2];

}


//----------------------------------------------------------------------------
void axisExtractor02::realtoindex(double a[3], int b[3] )
{
	double c[3];

	double minspac;

// EED 15 Mai 2007 .NET 
//	minspac=min(espprin[0],min(espprin[1],espprin[2]));

	minspac=espprin[0];
	if (espprin[1]<minspac) { minspac=espprin[1]; }
	if (espprin[2]<minspac) { minspac=espprin[2]; }


	c[0]=(a[0]/minspac);
	c[1]=(a[1]/minspac);
	c[2]=(a[2]/minspac);

	b[0]=(int)(a[0]/minspac);
	b[1]=(int)(a[1]/minspac);
	b[2]=(int)(a[2]/minspac);

	if(c[0]-b[0]>0.5){
		b[0]+=1;
	}
	if(c[1]-b[1]>0.5){
		b[1]+=1;
	}
	if(c[2]-b[2]>0.5){
		b[2]+=1;
	}
	

}	







//----------------------------------------------------------------------------
void axisExtractor02::indextoreal(int a[3], double b[3] )
{	
	double minspac;

// EED 15 Mai 2007 .NET 
//	minspac=min(espprin[0],min(espprin[1],espprin[2]));
	minspac=espprin[0];
	if (espprin[1]<minspac) { minspac=espprin[1]; }
	if (espprin[2]<minspac) { minspac=espprin[2]; }
	
	b[0]=(a[0])*minspac;
	b[1]=(a[1])*minspac;
	b[2]=a[2]*minspac;
}



//----------------------------------------------------------------------------
void axisExtractor02::indextoreal(double a[3], double b[3] )
{	
	double minspac;

// EED 15 Mai 2007 .NET 
//	minspac=min(espprin[0],min(espprin[1],espprin[2]));

	minspac=espprin[0];
	if (espprin[1]<minspac) { minspac=espprin[1]; }
	if (espprin[2]<minspac) { minspac=espprin[2]; }

	b[0]=(a[0])*minspac;
	b[1]=(a[1])*minspac;
	b[2]=a[2]*minspac;
}



//----------------------------------------------------------------------------
double axisExtractor02::distanciaejepunto(double point[3], double point2[3], double point3[3])
{
	
	
	double point4[3];

	double v[3], v2[3], v3[3];

	double x, dist=0,  y;

		
	v[0]=point3[0]-point2[0];
	v[1]=point3[1]-point2[1];
	v[2]=point3[2]-point2[2];

	v3[0]=point[0]-point2[0];
	v3[1]=point[1]-point2[1];
	v3[2]=point[2]-point2[2];
					
	x=((v[0]*v3[0])+(v[1]*v3[1])+(v[2]*v3[2]))/((v[0]*v[0])+(v[1]*v[1])+(v[2]*v[2]));
					
					
	point4[0]=point2[0]+(x*v[0]);
	point4[1]=point2[1]+(x*v[1]);
	point4[2]=point2[2]+(x*v[2]);

	v2[0]=point4[0]-point[0];
	v2[1]=point4[1]-point[1];
	v2[2]=point4[2]-point[2];

	y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);

	dist=sqrt((v2[0]*v2[0])+(v2[1]*v2[1])+(v2[2]*v2[2]));
					
	return dist;
	
}



//----------------------------------------------------------------------------
double axisExtractor02::proporcioejepunto(double point[3], double point2[3], double point3[3])
{
	
	
	double v[3], v3[3];

	double x;

	v[0]=point3[0]-point2[0];
	v[1]=point3[1]-point2[1];
	v[2]=point3[2]-point2[2];

	v3[0]=point[0]-point2[0];
	v3[1]=point[1]-point2[1];
	v3[2]=point[2]-point2[2];
					
	x=((v[0]*v3[0])+(v[1]*v3[1])+(v[2]*v3[2]))/((v[0]*v[0])+(v[1]*v[1])+(v[2]*v[2]));
					
	return x;
				
			
	
	
}






	

//----------------------------------------------------------------------------
void axisExtractor02::searc(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
{

	unsigned char *ptr;
	unsigned char *ptr2;	

	int radio;
	
	int i2, j2, k2;
	int i3, j3, k3;
		
	ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
	ptr2=(unsigned char *)	data2->GetScalarPointer(i,j,k);

	int ext[6];

	data->GetExtent(ext);

	radio=(ext[1]-ext[0])/2;

	i2=i-ext[0]-radio;
	j2=j-ext[2]-radio;
	k2=k-ext[4]-radio;

	*ptr2=label;
	vector[label-1][0]+=1;
	vector[label-1][1]+=i;
	vector[label-1][2]+=j;
	vector[label-1][3]+=k;

	for(i3=-1;i3<=1;i3++){
			for(j3=-1;j3<=1;j3++){
				for(k3=-1;k3<=1;k3++){
					if(i+i3>=ext[0] && i+i3<=ext[1] && j+j3>=ext[2] && j+j3<=ext[3] && k+k3>=ext[4] &&  k+k3<=ext[5] ){
						ptr=(unsigned char *)	data->GetScalarPointer(i+i3,j+j3,k+k3);
						ptr2=(unsigned char *)	data2->GetScalarPointer(i+i3,j+j3,k+k3);
						if(*ptr==255 && *ptr2==0 && (radio*radio)>=(((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3))) && ((radio-1)*(radio-1))<=(((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3)))){
							searc(i+i3, j+j3, k+k3, data, data2, label, vector);
						}
					}
				}
			}
		}
		
}




//----------------------------------------------------------------------------
void axisExtractor02::searcb(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
{

	unsigned char *ptr;
	unsigned char *ptr2;	

	int radio;

	int i2, j2, k2;
	int i3, j3, k3;
		
	ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
	ptr2=(unsigned char *)	data2->GetScalarPointer(i,j,k);

	int ext[6];

	data->GetExtent(ext);

	radio=(ext[1]-ext[0])/2;

	i2=i-ext[0]-radio;
	j2=j-ext[2]-radio;
	k2=k-ext[4]-radio;

	*ptr2=label;
	vector[label-1][0]+=1;
	vector[label-1][1]+=i;
	vector[label-1][2]+=j;
	vector[label-1][3]+=k;

	for(i3=-1;i3<=1;i3++){
			for(j3=-1;j3<=1;j3++){
				for(k3=-1;k3<=1;k3++){
					if(i+i3>=ext[0] && i+i3<=ext[1] && j+j3>=ext[2] && j+j3<=ext[3] && k+k3>=ext[4] &&  k+k3<=ext[5] ){
						ptr=(unsigned char *)	data->GetScalarPointer(i+i3,j+j3,k+k3);
						ptr2=(unsigned char *)	data2->GetScalarPointer(i+i3,j+j3,k+k3);
						if(*ptr==0 && *ptr2==0 && (radio*radio)>=(((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3))) && ((radio-1)*(radio-1))<=(((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3)))){
							searcb(i+i3, j+j3, k+k3, data, data2, label, vector);
						}
					}
				}
			}
		}
		
}




//----------------------------------------------------------------------------
unsigned char axisExtractor02::find_components(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
{
	int ext[6];
	int i, j, k, radio;
	int i2, j2, k2;

	data->GetExtent(ext);
	double *ff=  data->GetOrigin();


	unsigned char *ptr;
	unsigned char *ptr2;

	radio=(ext[1]-ext[0])/2;

	for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
			for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
				for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
					ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
					ptr2=(unsigned char *)	data2->GetScalarPointer(i,j,k);

					if(*ptr==255 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2))  ){
						label=label+1;
						vector[label-1][0]=0;
						vector[label-1][1]=0;
						vector[label-1][2]=0;
						vector[label-1][3]=0;
						searc(i, j, k, data, data2, label, vector);
						
					}
				}
			}
		}

		return label;
}




//----------------------------------------------------------------------------
unsigned char axisExtractor02::find_componentsb(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
{
	int ext[6];
	int i, j, k, radio;
	int i2, j2, k2;

	data->GetExtent(ext);

	unsigned char *ptr;
	unsigned char *ptr2;

	radio=(ext[1]-ext[0])/2;

	for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
			for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
				for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
					ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
					ptr2=(unsigned char *)	data2->GetScalarPointer(i,j,k);

					if(*ptr==0 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2))  ){
						label=label+1;
						vector[label-1][0]=0;
						vector[label-1][1]=0;
						vector[label-1][2]=0;
						vector[label-1][3]=0;
						searcb(i, j, k, data, data2, label, vector);
						
					}
				}
			}
		}

		return label;
}





//----------------------------------------------------------------------------
int axisExtractor02::proporcion(vtkImageData *data )
{
	int ext[6];
	int i, j, k;
	int max=0;
	int total=0;

	data->GetExtent(ext);

	unsigned short *ptr;
	
	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
				if(*ptr!=0){
					max++;
				}
				total++;
			}
		}
	}

	return (max*100)/total;
}








//----------------------------------------------------------------------------
bool axisExtractor02::border(vtkImageData *data, int p1[3] )
{
	int ext[6];
	int i, j, k;
	short val;
	int total1, total2;
	int centro[3];


	int p2[3];
	
	data->GetExtent(ext);

	centro[0]=(ext[1]+ext[0])/2;
	centro[1]=(ext[3]+ext[2])/2;
	centro[2]=(ext[5]+ext[4])/2;

	unsigned char *ptr;

	ptr=(unsigned char *)	data->GetScalarPointer(p1);
	
	val=*ptr;

	bool res=false;
	total1=0;
	total2=0;

	for(i=-1;i<=1;i++){
		for(j=-1;j<=1;j++){
			for(k=-1;k<=1;k++){
				p2[0]=p1[0]+i;
				p2[1]=p1[1]+j;
				p2[2]=p1[2]+k;
				if(ext[0]<=p2[0] && ext[1]>=p2[0] && ext[2]<=p2[1] && ext[3]>=p2[1] && ext[4]<=p2[2] && ext[5]>=p2[2] ){
					ptr=(unsigned char *)	data->GetScalarPointer(p2);
					if(*ptr!=val){
						total1++;
					}
					else{
						total2++;
					}						
				}
			}
		}
	}

	if(total1>0){
		res=true;
	}
	/*if(p1[0]==centro[0] && p1[1]==centro[1] && p1[2]==centro[2]){
		res=false;
	}*/
			
	return res;
}






//----------------------------------------------------------------------------
void axisExtractor02::optim(vtkImageData *data, vtkImageData *data2 )
{
	int ext[6];
	int i, j, k, radio;
	int centro[3];

	double w0, w1;

	double w0p, w1p;

	double inerciaxx, inerciayy, inerciazz, inerciaxy, inerciayz, inerciazx;
	double sumx, sumy, sumz;
	double sumxx, sumyy, sumzz, sumxy, sumyz, sumzx;

	double sumix, sumiy, sumiz;
	double sumixx, sumiyy, sumizz, sumixy, sumiyz, sumizx;

	double inerciaixx, inerciaiyy, inerciaizz, inerciaixy, inerciaiyz, inerciaizx;

	double inerciaxx2, inerciayy2, inerciazz2, inerciaxy2, inerciayz2, inerciazx2;
	double sumx2, sumy2, sumz2;
	double sumxx2, sumyy2, sumzz2, sumxy2, sumyz2, sumzx2;

	double sumix2, sumiy2, sumiz2;
	double sumixx2, sumiyy2, sumizz2, sumixy2, sumiyz2, sumizx2;

	//double sumixp, sumiyp, sumizp;
	//double sumix2p, sumiy2p, sumiz2p;

	double inerciaxxp, inerciayyp, inerciazzp, inerciaxyp, inerciayzp, inerciazxp;
	double sumxp, sumyp, sumzp;
	double sumxxp, sumyyp, sumzzp, sumxyp, sumyzp, sumzxp;

	double inerciaxx2p, inerciayy2p, inerciazz2p, inerciaxy2p, inerciayz2p, inerciazx2p;
	double sumx2p, sumy2p, sumz2p;
	double sumxx2p, sumyy2p, sumzz2p, sumxy2p, sumyz2p, sumzx2p;

	//double inerciaixx2p, inerciaiyy2p, inerciaizz2p, inerciaixy2p, inerciaiyz2p, inerciaizx2p;
	
	unsigned long cantp;
	double sumip;
	double sumiip;
	double inerciaiip;
	double centip;
	
	unsigned long cant2p;
	double sumi2p;
	double sumii2p;
	double inerciaii2p;
	double centi2p;
	
	int point[3];

	
	double sumitotal=0;
	double distcent=0;

	double im, jm, km, lm;
	

	max2=0;
	cant=0;
	sumi=0;
	sumii=0;
	des2=0;
	min2=3200;
	centi=0;

	max3=0;
	cant2=0;
	sumi2=0;
	sumii2=0;
	des3=0;
	min3=3200;
	centi2=0;

	sumx=0;
	sumy=0;
	sumz=0;

	sumxx=0;
	sumyy=0;
	sumzz=0;

	sumxy=0;
	sumyz=0;
	sumzx=0;

	sumix=0;
	sumiy=0;
	sumiz=0;

	sumixx=0;
	sumiyy=0;
	sumizz=0;

	sumixy=0;
	sumiyz=0;
	sumizx=0;

	sumx2=0;
	sumy2=0;
	sumz2=0;

	sumxx2=0;
	sumyy2=0;
	sumzz2=0;

	sumxy2=0;
	sumyz2=0;
	sumzx2=0;

	sumix2=0;
	sumiy2=0;
	sumiz2=0;

	sumixx2=0;
	sumiyy2=0;
	sumizz2=0;

	sumixy2=0;
	sumiyz2=0;
	sumizx2=0;

	max=0;
	min=3200;
		

	data->GetExtent(ext);

	unsigned short *ptr;
	unsigned char *ptr2;

	radio=(ext[1]-ext[0])/2;

	centro[0]=(ext[1]+ext[0])/2;
	centro[1]=(ext[3]+ext[2])/2;
	centro[2]=(ext[5]+ext[4])/2;
	
	double tmpsqrt;

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
				if( radio>=sqrt( tmpsqrt) ){
					ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
					if(*ptr>max){
						max=*ptr;
							
					}
					if(*ptr<min){
						min=*ptr;
							
					}
				}
			}
		}
	}

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
				if( radio>=sqrt(tmpsqrt) ){
					ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
					ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
					point[0]=i;
					point[1]=j;
					point[2]=k;
						
					im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
					jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
					km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
					lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
							
					if(*ptr2!=0){
						sumi+=lm;
						sumii+=lm*lm;
						cant+=1;

						sumx+=im;
						sumy+=jm;
						sumz+=km;

						sumxx+=im*im;
						sumyy+=jm*jm;
						sumzz+=km*km;
						sumxy+=im*jm;
						sumyz+=jm*km;
						sumzx+=km*im;

						sumix+=lm*im;
						sumiy+=lm*jm;
						sumiz+=lm*km;

						sumixx+=lm*im*im;
						sumiyy+=lm*jm*jm;
						sumizz+=lm*km*km;
						sumixy+=lm*im*jm;
						sumiyz+=lm*jm*km;
						sumizx+=lm*km*im;
					}
					else{
						sumi2+=lm;
						sumii2+=lm*lm;
						cant2+=1;
					
						sumx2+=im;
						sumy2+=jm;
						sumz2+=km;

						sumxx2+=im*im;
						sumyy2+=jm*jm;
						sumzz2+=km*km;
						sumxy2+=im*jm;
						sumyz2+=jm*km;
						sumzx2+=km*im;

						sumix2+=lm*im;
						sumiy2+=lm*jm;
						sumiz2+=lm*km;

						sumixx2+=lm*im*im;
						sumiyy2+=lm*jm*jm;
						sumizz2+=lm*km*km;
						sumixy2+=lm*im*jm;
						sumiyz2+=lm*jm*km;
						sumizx2+=lm*km*im;
					}
				}
			}
		}
	}

	centi=(double)sumi/(double)cant;
	inerciaii= ((double)sumii/(double)cant)-(centi*centi);
						
	centi2=(double)sumi2/(double)cant2;
	inerciaii2= ((double)sumii2/(double)cant2)-(centi2*centi2);

	centx=(double)sumx/(double)cant;
	centy=(double)sumy/(double)cant;
	centz=(double)sumz/(double)cant;

	centx2=(double)sumx2/(double)cant2;
	centy2=(double)sumy2/(double)cant2;
	centz2=(double)sumz2/(double)cant2;
	
	inerciaxx= ((double)sumxx/(double)cant)-centx*centx;
	inerciayy= ((double)sumyy/(double)cant)-centy*centy;
	inerciazz= ((double)sumzz/(double)cant)-centz*centz;

	inerciaxx2= ((double)sumxx2/(double)cant2)-centx2*centx2;
	inerciayy2= ((double)sumyy2/(double)cant2)-centy2*centy2;
	inerciazz2= ((double)sumzz2/(double)cant2)-centz2*centz2;

	inerciaxy= ((double)sumxy/(double)cant)-centx*centy;
	inerciayz= ((double)sumyz/(double)cant)-centy*centz;
	inerciazx= ((double)sumzx/(double)cant)-centz*centx;

	inerciaxy2= ((double)sumxy2/(double)cant2)-centx2*centy2;
	inerciayz2= ((double)sumyz2/(double)cant2)-centy2*centz2;
	inerciazx2= ((double)sumzx2/(double)cant2)-centz2*centx2;


	centix=((double)sumix+(double)sumix2)/((double)sumi+(double)sumi2);
	centiy=((double)sumiy+(double)sumiy2)/((double)sumi+(double)sumi2);
	centiz=((double)sumiz+(double)sumiz2)/((double)sumi+(double)sumi2);
		
	inerciaixx= ( ((double)sumixx+(double)sumixx2) / ((double)sumi+(double)sumi2)  )-(centix*centix);
	inerciaiyy= ( ((double)sumiyy+(double)sumiyy2) / ((double)sumi+(double)sumi2)  )-(centiy*centiy);
	inerciaizz= ( ((double)sumizz+(double)sumizz2) / ((double)sumi+(double)sumi2)  )-(centiz*centiz);


	inerciaixy= ( ((double)sumixy+(double)sumixy2) / ((double)sumi+(double)sumi2)  )-(centix*centiy);
	inerciaiyz= ( ((double)sumiyz+(double)sumiyz2) / ((double)sumi+(double)sumi2)  )-(centiy*centiz);
	inerciaizx= ( ((double)sumizx+(double)sumizx2) / ((double)sumi+(double)sumi2)  )-(centiz*centix);


	A[0][0]=inerciaxx;
	A[1][1]=inerciayy;
	A[2][2]=inerciazz;
	A[0][1]=A[1][0]=inerciaxy;
	A[0][2]=A[2][0]=inerciazx;
	A[1][2]=A[2][1]=inerciayz;

	A2[0][0]=inerciaxx2;
	A2[1][1]=inerciayy2;
	A2[2][2]=inerciazz2;
	A2[0][1]=A2[1][0]=inerciaxy2;
	A2[0][2]=A2[2][0]=inerciazx2;
	A2[1][2]=A2[2][1]=inerciayz2;

	Ai[0][0]=inerciaixx;
	Ai[1][1]=inerciaiyy;
	Ai[2][2]=inerciaizz;
	Ai[0][1]=Ai[1][0]=inerciaixy;
	Ai[0][2]=Ai[2][0]=inerciaizx;
	Ai[1][2]=Ai[2][1]=inerciaiyz;

	vtkMath::Diagonalize3x3  (A, w, V);
	vtkMath::Diagonalize3x3  (A2, w2, V2);
	vtkMath::Diagonalize3x3  (Ai, wi, Vi);

	
	if(w[0]>w[1]){
		if(w[0]>w[2]){
			ejemin=0;
			inerciar=w[1]+w[2];
			if(w[1]>w[2]){
				inerciary=w[0]+w[2];
				inerciarz=w[0]+w[1];
			}
			else{
				inerciary=w[0]+w[1];
				inerciarz=w[0]+w[2];
			}
		}
		else{
			ejemin=2;
			inerciar=w[0]+w[1];
			if(w[1]>w[0]){
				inerciary=w[0]+w[2];
				inerciarz=w[1]+w[2];
			}
			else{
				inerciary=w[1]+w[2];
				inerciarz=w[0]+w[2];
			}
		}
	}
	else{
		if(w[1]>w[2]){
			ejemin=1;
			inerciar=w[2]+w[0];
			if(w[2]>w[0]){
				inerciary=w[1]+w[0];
				inerciarz=w[1]+w[2];
			}
			else{
				inerciary=w[1]+w[2];
				inerciarz=w[1]+w[0];
			}
		}
		else{
			ejemin=2;
			inerciar=w[1]+w[0];
			if(w[1]>w[0]){
				inerciary=w[0]+w[2];
				inerciarz=w[1]+w[2];
			}
			else{
				inerciary=w[1]+w[2];
				inerciarz=w[0]+w[2];
			}
		}
	}


	if(w2[0]>w2[1]){
		if(w2[0]>w2[2]){
			ejemin2=0;
			inercia2r=w2[1]+w2[2];
		}
		else{
			ejemin2=2;
			inercia2r=w2[1]+w2[0];
		}
	}
	else{
		if(w2[1]>w2[2]){
			ejemin2=1;
			inercia2r=w2[2]+w2[0];
		}
		else{
			ejemin2=2;
			inercia2r=w2[1]+w2[0];
		}
	}


	if(wi[0]>wi[1]){
		if(wi[0]>wi[2]){
			ejemini=0;
			inerciari=wi[1]+wi[2];
			if(wi[1]>wi[2]){
				inerciariy=wi[0]+wi[2];
				inerciariz=wi[0]+wi[1];
			}
			else{
				inerciariy=wi[0]+wi[1];
				inerciariz=wi[0]+wi[2];
			}
		}
		else{
			ejemini=2;
			inerciari=wi[1]+wi[0];
			if(wi[1]>wi[0]){
				inerciariy=wi[0]+wi[2];
				inerciariz=wi[1]+wi[2];
			}
			else{
				inerciariy=wi[1]+wi[2];
				inerciariz=wi[0]+wi[2];
			}

		}
	}
	else{
		if(wi[1]>wi[2]){
			ejemini=1;
			inerciari=wi[2]+wi[0];
			if(wi[2]>wi[0]){
				inerciariy=wi[1]+wi[0];
				inerciariz=wi[1]+wi[2];
			}
			else{
				inerciariy=wi[1]+wi[2];
				inerciariz=wi[1]+wi[0];
			}
		}
		else{
			ejemini=2;
			inerciari=wi[1]+wi[0];
			if(wi[1]>wi[0]){
				inerciariy=wi[0]+wi[2];
				inerciariz=wi[1]+wi[2];
			}
			else{
				inerciariy=wi[1]+wi[2];
				inerciariz=wi[0]+wi[2];
			}
		}
	}

		
	/*
	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				if( radio>=sqrt(((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]))) ){
					ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
					ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
					
					
					im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
					jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
					km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);

					point1[0]=im;
					point1[1]=jm;
					point1[2]=km;

					point2[0]=centx+V[0][ejemin];
					point2[1]=centy+V[1][ejemin];
					point2[2]=centz+V[2][ejemin];

					point3[0]=centx-V[0][ejemin];
					point3[1]=centy-V[1][ejemin];
					point3[2]=centz-V[2][ejemin];

					//distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
					distcent=1;
						
					if(*ptr2==0){				
						
						sumitotal+=(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
						
					}
				}
			}
		}
	}

	inercia2r=sumitotal/cant2;
	*/




	/*printf("centx: %f\n", centx);
	printf("centy: %f\n", centy);
	printf("centz: %f\n", centz);

	printf("w[0]: %f\n", w[0]);
	printf("w[1]: %f\n", w[1]);
	printf("w[2]: %f\n", w[2]);
	printf("inerciar: %f\n", inerciar);

	printf("w[ejemin]: %f\n", w[ejemin]);

	printf("w[ejemin]/inerciar: %f\n", w[ejemin]/inerciar);

	printf("V[0][0]: %f\n", V[0][0]);
	printf("V[1][0]: %f\n", V[1][0]);
	printf("V[2][0]: %f\n", V[2][0]);
	
	printf("V[0][1]: %f\n", V[0][1]);
	printf("V[1][1]: %f\n", V[1][1]);
	printf("V[2][1]: %f\n", V[2][1]);
	
	printf("V[0][2]: %f\n", V[0][2]);
	printf("V[1][2]: %f\n", V[1][2]);
	printf("V[2][2]: %f\n", V[2][2]);*/

	w0=(double)cant/((double)cant+(double)cant2);
	w1=(double)cant2/((double)cant+(double)cant2);

	
	/*
	ut=w0*centi+w1*centi2;
	intervar=w0*inerciaii+w1*inerciaii2;
	*/

	//w0=1;
	//w1=1;

	

	//costo=w0*inerciaii+w1*inerciaii2+param*inerciar;
	//costo=param*(w0*param3*inerciaii+w1*(1-param3)*inerciaii2)+(1-param)*(w0*param2*inerciar+w1*(1-param2)*inercia2r);
	costo=param*w0*inerciaii+w1*param2*inerciaii2+w0*param3*inerciar+w1*param4*inercia2r;

	int changes=4;

	int total=0;
	
	while(changes>3 && total<500){
		changes=0;
		total++;
		for(i=ext[0];i<=ext[1];i++){
			for(j=ext[2];j<=ext[3];j++){
				for(k=ext[4];k<=ext[5];k++){
					tmpsqrt = ((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
					if( radio>=sqrt(tmpsqrt) ){
						ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
						ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
						point[0]=i;
						point[1]=j;
						point[2]=k;

						im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
						jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
						km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
						lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
						
						if(border(data2, point)){
							if(*ptr2!=0){
								
								sumip=sumi-lm;
								sumiip=sumii-lm*lm;
							
								cantp=cant-1;

								sumi2p=sumi2+lm;
								sumii2p=sumii2+lm*lm;
								
								cant2p=cant2+1;

								sumxp=sumx-im;
								sumyp=sumy-jm;
								sumzp=sumz-km;
								sumxxp=sumxx-im*im;
								sumyyp=sumyy-jm*jm;
								sumzzp=sumzz-km*km;
								sumxyp=sumxy-im*jm;
								sumyzp=sumyz-jm*km;
								sumzxp=sumzx-km*im;

								sumx2p=sumx2+im;
								sumy2p=sumy2+jm;
								sumz2p=sumz2+km;
								sumxx2p=sumxx2+im*im;
								sumyy2p=sumyy2+jm*jm;
								sumzz2p=sumzz2+km*km;
								sumxy2p=sumxy2+im*jm;
								sumyz2p=sumyz2+jm*km;
								sumzx2p=sumzx2+km*im;

								centip=(double)sumip/(double)cantp;
								inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
								
								centi2p=(double)sumi2p/(double)cant2p;
								inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
								
								centxp=(double)sumxp/(double)cantp;
								centyp=(double)sumyp/(double)cantp;
								centzp=(double)sumzp/(double)cantp;

								centx2p=(double)sumx2p/(double)cant2p;
								centy2p=(double)sumy2p/(double)cant2p;
								centz2p=(double)sumz2p/(double)cant2p;
																
								inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
								inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
								inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;

								inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
								inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
								inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;

								inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
								inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
								inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;

								inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
								inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
								inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;

								Ap[0][0]=inerciaxxp;
								Ap[1][1]=inerciayyp;
								Ap[2][2]=inerciazzp;
								Ap[0][1]=Ap[1][0]=inerciaxyp;
								Ap[0][2]=Ap[2][0]=inerciazxp;
								Ap[1][2]=Ap[2][1]=inerciayzp;


								A2p[0][0]=inerciaxx2p;
								A2p[1][1]=inerciayy2p;
								A2p[2][2]=inerciazz2p;
								A2p[0][1]=A2p[1][0]=inerciaxy2p;
								A2p[0][2]=A2p[2][0]=inerciazx2p;
								A2p[1][2]=A2p[2][1]=inerciayz2p;

								

								w0p=(double)cantp/((double)cantp+(double)cant2p);
								w1p=(double)cant2p/((double)cantp+(double)cant2p);

								//w0p=1;
								//w1p=1;
								
								vtkMath::Diagonalize3x3  (Ap, wp, Vp);
								vtkMath::Diagonalize3x3  (A2p, w2p, V2p);

							
								if(wp[0]>wp[1]){
									if(wp[0]>wp[2]){
										ejeminp=0;
										inerciarp=wp[2]+wp[1];
										if(wp[1]>wp[2]){
											inerciarpy=wp[0]+wp[2];
											inerciarpz=wp[0]+wp[1];
										}
										else{
											inerciarpy=wp[0]+wp[1];
											inerciarpz=wp[0]+wp[2];
										}
									}
									else{
										ejeminp=2;
										inerciarp=wp[1]+wp[0];
										if(wp[1]>wp[0]){
											inerciarpy=wp[0]+wp[2];
											inerciarpz=wp[1]+wp[2];
										}
										else{
											inerciarpy=wp[1]+wp[2];
											inerciarpz=wp[0]+wp[2];
										}
									}
								}
								else{
									if(wp[1]>wp[2]){
										ejeminp=1;
										inerciarp=wp[2]+wp[0];
										if(wp[2]>wp[0]){
											inerciarpy=wp[1]+wp[0];
											inerciarpz=wp[1]+wp[2];
										}
										else{
											inerciarpy=wp[1]+wp[2];
											inerciarpz=wp[1]+wp[0];
										}
									}
									else{
										ejeminp=2;
										inerciarp=wp[1]+wp[0];
										if(wp[1]>wp[0]){
											inerciarpy=wp[0]+wp[2];
											inerciarpz=wp[1]+wp[2];
										}
										else{
											inerciarpy=wp[1]+wp[2];
											inerciarpz=wp[0]+wp[2];
										}
									}
								}


								if(w2p[0]>w2p[1]){
									if(w2p[0]>w2p[2]){
										ejemin2p=0;
										inercia2rp=w2p[2]+w2p[1];
									}
									else{
										ejemin2p=2;
										inercia2rp=w2p[1]+w2p[0];
									}
								}
								else{
									if(w2p[1]>w2p[2]){
										ejemin2p=1;
										inercia2rp=w2p[2]+w2p[0];
									}
									else{
										ejemin2p=2;
										inercia2rp=w2p[1]+w2p[0];
									}
								}

							/*	point1[0]=im;
								point1[1]=jm;
								point1[2]=km;

								point2[0]=centxp+Vp[0][ejeminp];
								point2[1]=centyp+Vp[1][ejeminp];
								point2[2]=centzp+Vp[2][ejeminp];

								point3[0]=centxp-Vp[0][ejeminp];
								point3[1]=centyp-Vp[1][ejeminp];
								point3[2]=centzp-Vp[2][ejeminp];*/

								//distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
							/*	distcent=1;

								sumitotalp=sumitotal+(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
							
								inercia2rp=sumitotalp/cant2p;*/

								//costop=w0p*inerciaiip+w1p*inerciaii2p+param*inerciarp;

								//costop=param*(w0p*param3*inerciaiip+w1p*(1-param3)*inerciaii2p)+(1-param)*(w0p*param2*inerciarp+w1p*(1-param2)*inercia2rp);
								
								
								costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;

								if(costop<costo  && centip>centi2p){
																		
									*ptr2=0;
									changes++;

									sumi=sumip;
									sumii=sumiip;
									cant=cantp;

									sumi2=sumi2p;
									sumii2=sumii2p;
									cant2=cant2p;

									sumx=sumxp;
									sumy=sumyp;
									sumz=sumzp;
									sumxx=sumxxp;
									sumyy=sumyyp;
									sumzz=sumzzp;
									sumxy=sumxyp;
									sumyz=sumyzp;
									sumzx=sumzxp;

									sumx2=sumx2p;
									sumy2=sumy2p;
									sumz2=sumz2p;
									sumxx2=sumxx2p;
									sumyy2=sumyy2p;
									sumzz2=sumzz2p;
									sumxy2=sumxy2p;
									sumyz2=sumyz2p;
									sumzx2=sumzx2p;

									centi=centip;
									inerciaii= inerciaiip;
									
									centi2=centi2p;
									inerciaii2p= inerciaii2p;
									
									centx=centxp;
									centy=centyp;
									centz=centzp;

									centx2=centx2p;
									centy2=centy2p;
									centz2=centz2p;
																		
									inerciaxx= inerciaxxp;
									inerciayy= inerciayyp;
									inerciazz= inerciazzp;

									inerciaxx2= inerciaxx2p;
									inerciayy2= inerciayy2p;
									inerciazz2= inerciazz2p;

									inerciaxy= inerciaxyp;
									inerciayz= inerciayzp;
									inerciazx= inerciazxp;

									inerciaxy2= inerciaxy2p;
									inerciayz2= inerciayz2p;
									inerciazx2= inerciazx2p;

								//	sumitotal=sumitotalp;

									costo=costop;

								}
							}
							else{
								sumip=sumi+lm;
								sumiip=sumii+lm*lm;
							
								cantp=cant+1;

								sumi2p=sumi2-lm;
								sumii2p=sumii2-lm*lm;
								
								cant2p=cant2-1;

								sumxp=sumx+im;
								sumyp=sumy+jm;
								sumzp=sumz+km;
								sumxxp=sumxx+im*im;
								sumyyp=sumyy+jm*jm;
								sumzzp=sumzz+km*km;
								sumxyp=sumxy+im*jm;
								sumyzp=sumyz+jm*km;
								sumzxp=sumzx+km*im;

								sumx2p=sumx2-im;
								sumy2p=sumy2-jm;
								sumz2p=sumz2-km;
								sumxx2p=sumxx2-im*im;
								sumyy2p=sumyy2-jm*jm;
								sumzz2p=sumzz2-km*km;
								sumxy2p=sumxy2-im*jm;
								sumyz2p=sumyz2-jm*km;
								sumzx2p=sumzx2-km*im;

								centip=(double)sumip/(double)cantp;
								inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
								
								centi2p=(double)sumi2p/(double)cant2p;
								inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
								
								centxp=(double)sumxp/(double)cantp;
								centyp=(double)sumyp/(double)cantp;
								centzp=(double)sumzp/(double)cantp;

								centx2p=(double)sumx2p/(double)cant2p;
								centy2p=(double)sumy2p/(double)cant2p;
								centz2p=(double)sumz2p/(double)cant2p;
								
								
								inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
								inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
								inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;

								inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
								inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
								inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;

								inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
								inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
								inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;

								inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
								inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
								inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;

								Ap[0][0]=inerciaxxp;
								Ap[1][1]=inerciayyp;
								Ap[2][2]=inerciazzp;
								Ap[0][1]=Ap[1][0]=inerciaxyp;
								Ap[0][2]=Ap[2][0]=inerciazxp;
								Ap[1][2]=Ap[2][1]=inerciayzp;

								A2p[0][0]=inerciaxx2p;
								A2p[1][1]=inerciayy2p;
								A2p[2][2]=inerciazz2p;
								A2p[0][1]=A2p[1][0]=inerciaxy2p;
								A2p[0][2]=A2p[2][0]=inerciazx2p;
								A2p[1][2]=A2p[2][1]=inerciayz2p;

								w0p=(double)cantp/((double)cantp+(double)cant2p);
								w1p=(double)cant2p/((double)cantp+(double)cant2p);

								//w0p=1;
								//w1p=1;
								
								vtkMath::Diagonalize3x3  (Ap, wp, Vp);
								vtkMath::Diagonalize3x3  (A2p, w2p, V2p);
							
								if(wp[0]>wp[1]){
									if(wp[0]>wp[2]){
										ejeminp=0;
										inerciarp=wp[2]+wp[1];
										if(wp[1]>wp[2]){
											inerciarpy=wp[0]+wp[2];
											inerciarpz=wp[0]+wp[1];
										}
										else{
											inerciarpy=wp[0]+wp[1];
											inerciarpz=wp[0]+wp[2];
										}
									}
									else{
										ejeminp=2;
										inerciarp=wp[1]+wp[0];
										if(wp[1]>wp[0]){
											inerciarpy=wp[0]+wp[2];
											inerciarpz=wp[1]+wp[2];
										}
										else{
											inerciarpy=wp[1]+wp[2];
											inerciarpz=wp[0]+wp[2];
										}
									}
								}
								else{
									if(wp[1]>wp[2]){
										ejeminp=1;
										inerciarp=wp[2]+wp[0];
										if(wp[2]>wp[0]){
											inerciarpy=wp[1]+wp[0];
											inerciarpz=wp[1]+wp[2];
										}
										else{
											inerciarpy=wp[1]+wp[2];
											inerciarpz=wp[1]+wp[0];
										}
									}
									else{
										ejeminp=2;
										inerciarp=wp[1]+wp[0];
										if(wp[1]>wp[0]){
											inerciarpy=wp[0]+wp[2];
											inerciarpz=wp[1]+wp[2];
										}
										else{
											inerciarpy=wp[1]+wp[2];
											inerciarpz=wp[0]+wp[2];
										}
									}
								}


								if(w2p[0]>w2p[1]){
									if(w2p[0]>w2p[2]){
										ejemin2p=0;
										inercia2rp=w2p[2]+w2p[1];
									}
									else{
										ejemin2p=2;
										inercia2rp=w2p[1]+w2p[0];
									}
								}
								else{
									if(w2p[1]>w2p[2]){
										ejemin2p=1;
										inercia2rp=w2p[2]+w2p[0];
									}
									else{
										ejemin2p=2;
										inercia2rp=w2p[1]+w2p[0];
									}
								}


							/*	point1[0]=im;
								point1[1]=jm;
								point1[2]=km;

								point2[0]=centxp+Vp[0][ejeminp];
								point2[1]=centyp+Vp[1][ejeminp];
								point2[2]=centzp+Vp[2][ejeminp];

								point3[0]=centxp-Vp[0][ejeminp];
								point3[1]=centyp-Vp[1][ejeminp];
								point3[2]=centzp-Vp[2][ejeminp];*/

								//distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
							/*	distcent=1;
							
								sumitotalp=sumitotal-(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
								
								inercia2rp=sumitotalp/cant2p;*/


								//costop=w0p*inerciaiip+w1p*inerciaii2p+param*inerciarp;
								//costop=param*(w0p*param3*inerciaiip+w1p*(1-param3)*inerciaii2p)+(1-param)*(w0p*param2*inerciarp+w1p*(1-param2)*inercia2rp);
								
								
								costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;

								if(costop<costo && centip>centi2p){
																		
									*ptr2=255;
									changes++;

									sumi=sumip;
									sumii=sumiip;
									cant=cantp;

									sumi2=sumi2p;
									sumii2=sumii2p;
									cant2=cant2p;

									sumx=sumxp;
									sumy=sumyp;
									sumz=sumzp;
									sumxx=sumxxp;
									sumyy=sumyyp;
									sumzz=sumzzp;
									sumxy=sumxyp;
									sumyz=sumyzp;
									sumzx=sumzxp;

									sumx2=sumx2p;
									sumy2=sumy2p;
									sumz2=sumz2p;
									sumxx2=sumxx2p;
									sumyy2=sumyy2p;
									sumzz2=sumzz2p;
									sumxy2=sumxy2p;
									sumyz2=sumyz2p;
									sumzx2=sumzx2p;

									centi=centip;
									inerciaii= inerciaiip;
									
									centi2=centi2p;
									inerciaii2p= inerciaii2p;
									
									centx=centxp;
									centy=centyp;
									centz=centzp;

									centx2=centx2p;
									centy2=centy2p;
									centz2=centz2p;
																		
									inerciaxx= inerciaxxp;
									inerciayy= inerciayyp;
									inerciazz= inerciazzp;

									inerciaxx2= inerciaxx2p;
									inerciayy2= inerciayy2p;
									inerciazz2= inerciazz2p;

									inerciaxy= inerciaxyp;
									inerciayz= inerciayzp;
									inerciazx= inerciazxp;

									inerciaxy2= inerciaxy2p;
									inerciayz2= inerciayz2p;
									inerciazx2= inerciazx2p;

								//	sumitotal=sumitotalp;

									costo=costop;

								}
								
							}
						}
					}
				}
			}
		}

	//	printf("changes: %d\n", changes);
		

	}

	//printf("total: %d\n", total);

		
}





//----------------------------------------------------------------------------
void axisExtractor02::costominimo(vtkImageData *data,  vtkImageData *data2 )
{
	int ext[6];
	int i, j, k, radio;
	int centro[3];
	
	
	data->GetExtent(ext);

	double *ptr;
	unsigned char *ptr2;

	int ind;

	radio=(ext[1]-ext[0])/2;

	centro[0]=(ext[1]+ext[0])/2;
	centro[1]=(ext[3]+ext[2])/2;
	centro[2]=(ext[5]+ext[4])/2;
	

	for(i=0;i<=10;i++){
		minis[i]=0;
	}

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(double *)data->GetScalarPointer(i,j,k);
				ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);

				ind=(int)*ptr2;

				if(*ptr2!=0){	
					if(minis[ind-1]<=*ptr && *ptr!=0 && ind<11){
						minis[ind-1]=*ptr;
						candit[ind-1][0]=i;
						candit[ind-1][1]=j;
						candit[ind-1][2]=k;
					}
				}
			}
		}
	}
}





//----------------------------------------------------------------------------
void axisExtractor02::costominimo2(vtkImageData *data,  vtkImageData *data3, int p1[3], int p2[3], int p3[3])
{
	int ext[6];
	int i, j, k, radio;
	int punto[3];
	int puntob[3];
	int punto2[3];
	int puntob2[3];		
	int puntomin[3];
	int puntobmin[3];
	
	double *ptr;
	double *ptr2; 
	unsigned char *ptr3; 
	unsigned char *ptr4; 
	double *ptr5; 
	unsigned char *ptr6; 
	double *ptr7; 
	double *ptr8;

	vtkImageData *data4;
		
	data4=vtkImageData::New();
	data4->SetScalarTypeToUnsignedChar();
	data4->SetExtent(data->GetExtent());
	data4->SetSpacing(data->GetSpacing());

	vtkImageData *data6;
		
	data6=vtkImageData::New();
	data6->SetScalarTypeToUnsignedChar();
	data6->SetExtent(data->GetExtent());
	data6->SetSpacing(data->GetSpacing());
	

	vtkImageData *data2;

	data2=vtkImageData::New();
	data2->SetScalarTypeToDouble();
	data2->SetExtent(data->GetExtent());
	data2->SetSpacing(data->GetSpacing());

	vtkImageData *data8;

	data8=vtkImageData::New();
	data8->SetScalarTypeToDouble();
	data8->SetExtent(data->GetExtent());
	data8->SetSpacing(data->GetSpacing());

	
	double mincst, mincstb;
	
	bool flag=true;

	double costoactual;

	data->GetExtent(ext);

	radio=(ext[1]-ext[0])/2;

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr2=(double *)data2->GetScalarPointer(i,j,k);
				ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
				ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
				ptr8=(double *)data8->GetScalarPointer(i,j,k);
				
				*ptr2 = 1000000000;
				*ptr4 = 0;
				*ptr6 = 0;
				*ptr8 = 1000000000;
			}
		}
	}



	punto[0]=p1[0];
	punto[1]=p1[1];
	punto[2]=p1[2];
	puntob[0]=p2[0];
	puntob[1]=p2[1];
	puntob[2]=p2[2];
	
	
	ptr2=(double *)data2->GetScalarPointer(p1[0],p1[1],p1[2]);
	*ptr2=0;

	ptr4=(unsigned char *)data4->GetScalarPointer(p1[0],p1[1],p1[2]);
	*ptr4=1;

	ptr8=(double *)data8->GetScalarPointer(p2[0],p2[1],p2[2]);
	*ptr8=0;

	ptr6=(unsigned char *)data6->GetScalarPointer(p2[0],p2[1],p2[2]);
	*ptr6=1;

	


		

	while(flag){
		
		ptr5=(double *)data2->GetScalarPointer(punto[0],punto[1],punto[2]);
		ptr7=(double *)data8->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
			
			for(i=-1;i<=1;i++){
				for(j=-1;j<=1;j++){
					for(k=-1;k<=1;k++){
						if(i!=0 || j!=0 || k!=0){
							punto2[0]=punto[0]+i;
							punto2[1]=punto[1]+j;
							punto2[2]=punto[2]+k;

							puntob2[0]=puntob[0]+i;
							puntob2[1]=puntob[1]+j;
							puntob2[2]=puntob[2]+k;

							if(punto2[0]>=ext[0] && punto2[0]<=ext[1] && punto2[1]>=ext[2] && punto2[1]<=ext[3] && punto2[2]>=ext[4] &&  punto2[2]<=ext[5] ){

								ptr=(double *)data->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
								ptr2=(double *)data2->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
								ptr3=(unsigned char *)data3->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
								ptr4=(unsigned char *)data4->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
								

								if(*ptr3!=0 && *ptr4!=2){
									costoactual=*ptr5+(1/(*ptr));
									if(*ptr2>costoactual){
										*ptr2=costoactual;									
										*ptr4=1;
									}									
								}
							}

							if(puntob2[0]>=ext[0] && puntob2[0]<=ext[1] && puntob2[1]>=ext[2] && puntob2[1]<=ext[3] && puntob2[2]>=ext[4] &&  puntob2[2]<=ext[5] ){

								ptr=(double *)data->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
								ptr8=(double *)data8->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
								ptr3=(unsigned char *)data3->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
								ptr6=(unsigned char *)data6->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);

								if(*ptr3!=0 && *ptr6!=2){
									costoactual=*ptr7+(1/(*ptr));
									if(*ptr8>costoactual){
										*ptr8=costoactual;									
										*ptr6=1;
									}									
								}
							}
						}
					}
				}
			}

			ptr4=(unsigned char *)data4->GetScalarPointer(punto[0],punto[1],punto[2]);
			*ptr4=2;

			ptr6=(unsigned char *)data6->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
			*ptr6=2;

			mincst=1000000000;
			mincstb=1000000000;
			
			for(i=ext[0];i<=ext[1];i++){
				for(j=ext[2];j<=ext[3];j++){
					for(k=ext[4];k<=ext[5];k++){
						ptr2=(double *)data2->GetScalarPointer(i,j,k);
						ptr8=(double *)data8->GetScalarPointer(i,j,k);
						ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
						ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
							
						if(*ptr4==1){
							if(mincst>=*ptr2){
								mincst=*ptr2;
								puntomin[0]=i;
								puntomin[1]=j;
								puntomin[2]=k;
							}
						}
						
						if(*ptr6==1){
							if(mincstb>=*ptr8){
								mincstb=*ptr8;
								puntobmin[0]=i;
								puntobmin[1]=j;
								puntobmin[2]=k;
							}
						}	
					}
				}
			}

			punto[0]=puntomin[0];
			punto[1]=puntomin[1];
			punto[2]=puntomin[2];

			puntob[0]=puntobmin[0];
			puntob[1]=puntobmin[1];
			puntob[2]=puntobmin[2];

			ptr4=(unsigned char *)data4->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
			ptr6=(unsigned char *)data6->GetScalarPointer(punto[0],punto[1],punto[2]);
		

			if(*ptr4==2){
				p3[0]=puntob[0];
				p3[1]=puntob[1];
				p3[2]=puntob[2];
				flag=false;
			}

			if(*ptr6==2){
				p3[0]=punto[0];
				p3[1]=punto[1];
				p3[2]=punto[2];
				flag=false;
			}
		}

		data2->Delete();
		data4->Delete();
		data6->Delete();
		data8->Delete();
}





//----------------------------------------------------------------------------
void axisExtractor02::invertir(vtkImageData *data )
{
	int ext[6];
	int i, j, k, radio;
	int i2, j2, k2;

	data->GetExtent(ext);

	unsigned char *ptr;

	radio=(ext[1]-ext[0])/2;

	double tmpsqrt;

	for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
		for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
			for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
				ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
				tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);		
				if( radio<sqrt(tmpsqrt) ){
					if(*ptr==0){
						*ptr=255;
					}
					else{
						*ptr=0;
					}
				}
			}
		}
	}
}



//----------------------------------------------------------------------------
void axisExtractor02::redondear(vtkImageData *data )
{
	int ext[6];
	int i, j, k, radio;
	int i2, j2, k2;

	data->GetExtent(ext);

	unsigned char *ptr;

	radio=(ext[1]-ext[0])/2;

	double tmpsqrt;

	for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
		for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
			for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
				ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
				tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);			
				if( radio<sqrt(tmpsqrt) ){
					*ptr=0;
				}
				/*if( i2==0 && j2==0 && k2==0 ){
					*ptr=255;
				}*/
			}
		}
	}
}





//----------------------------------------------------------------------------
void axisExtractor02::redondear2(vtkImageData *data )
{
	int ext[6];
	int i, j, k, radio;
	int i2, j2, k2;

	data->GetExtent(ext);

	double *ptr;

	radio=(ext[1]-ext[0])/2;

	double tmpsqrt;

	for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
		for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
			for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
				ptr=(double *)	data->GetScalarPointer(i,j,k);
				tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);	
				if( radio<sqrt(tmpsqrt) ){
					*ptr=0;
					
				}
				
			}
		}
	}
}




//----------------------------------------------------------------------------
void axisExtractor02::redondear3(vtkImageData *data )
{
	int ext[6];
	int i, j, k, radio;
	int i2, j2, k2;

	data->GetExtent(ext);

	unsigned char *ptr;

	radio=(ext[1]-ext[0])/2;

	double tmpsqrt;

	for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
		for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
			for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
				ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
				tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
				if( radio<sqrt(tmpsqrt) ){
					*ptr=255;
					
				}
				
			}
		}
	}
}




//----------------------------------------------------------------------------
double axisExtractor02::distancia(double a[3], double b[3] )
{
	double tmpsqrt=((a[0]-b[0])*(a[0]-b[0]))+((a[1]-b[1])*(a[1]-b[1]))+((a[2]-b[2])*(a[2]-b[2]));
	return sqrt(tmpsqrt);
}





//----------------------------------------------------------------------------
double axisExtractor02::distancia(int a[3], int b[3] )
{
	double tmpsqrt=((a[0]-b[0])*(a[0]-b[0]))+((a[1]-b[1])*(a[1]-b[1]))+((a[2]-b[2])*(a[2]-b[2]));
	return sqrt(tmpsqrt);
}







//----------------------------------------------------------------------------
void axisExtractor02::blanquear(vtkImageData *data )
{
	int ext[6];
	int i, j, k;
	
	data->GetExtent(ext);

	unsigned char *ptr;

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
				*ptr=0;
			}
		}
	}
}



//----------------------------------------------------------------------------
void axisExtractor02::blanquear3(vtkImageData *data )
{
	int ext[6];
	int i, j, k;
	
	data->GetExtent(ext);

	double *ptr;

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(double *)	data->GetScalarPointer(i,j,k);
				*ptr=0;
			}
		}
	}
}




//----------------------------------------------------------------------------
void axisExtractor02::blanquear2(vtkImageData *data )
{
	int ext[6];
	int centro[3];
	int i, j, k;
	
	data->GetExtent(ext);

	centro[0]=(ext[1]+ext[0])/2;
	centro[1]=(ext[3]+ext[2])/2;
	centro[2]=(ext[5]+ext[4])/2;

	unsigned char *ptr;

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
				if(i==centro[0] && j==centro[1] && k==centro[2]){
					*ptr=255;
				}
				else{
					*ptr=0;
				}
			}
		}
	}
}




//----------------------------------------------------------------------------
void axisExtractor02::cilindro(vtkImageData *data, double vector[3] )
{
	int ext[6];
	int i, j, k, radio;
	int i2, j2, k2;
	int centro[3];
	int punto[3];	
	double centrof[3];
	double puntof[3];
	double radiof;
	double puntof2[3];
	

	data->GetExtent(ext);

	unsigned char *ptr;

	radio=(ext[1]-ext[0])/2;
	radiof=espprin[0]*(((double)ext[1]-(double)ext[0])/4);
	centro[0]=(ext[1]+ext[0])/2;
	centro[1]=(ext[3]+ext[2])/2;
	centro[2]=(ext[5]+ext[4])/2;

	indextoreal(centro, centrof );
	
	if(vector[0]==0 && vector[1]==0 && vector[2]==0){
		blanquear2(data);
	}
	else{
		puntof2[0]=centrof[0]+vector[0];
		puntof2[1]=centrof[1]+vector[1];
		puntof2[2]=centrof[2]+vector[2];

		double tmpsqrt;
		for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
			for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
				for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
					ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
					punto[0]=i;
					punto[1]=j;
					punto[2]=k;
					indextoreal(punto, puntof );
					tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
					if( radio<sqrt(tmpsqrt) ){
						*ptr=0;
						
					}
					else if(radiof<distanciaejepunto(puntof, centrof, puntof2)){
						*ptr=0;
					}
					else{
						*ptr=255;
					}					
				}
			}
		}
	}
}





//----------------------------------------------------------------------------
void axisExtractor02::modelo(vtkImageData *data, unsigned char cantidad,  unsigned long vector[50][4], int candit[10][3], double radioactual, double minis[10])
{
	int ext[6];
	int i, j, k, radio, radio2;
	int i2, j2, k2;
	int t;
	int centro[3];
	int punto[3];
	int punto2[3];
	double centrof[3];
	double puntof[3];
	double radiof;
	double puntof2[3];

	double dist, prop, rest;
	

	data->GetExtent(ext);

	unsigned char *ptr;

	radio=(ext[1]-ext[0])/2;
	radio2=(int)radioactual;
	
	centro[0]=(ext[1]+ext[0])/2;
	centro[1]=(ext[3]+ext[2])/2;
	centro[2]=(ext[5]+ext[4])/2;

	indextoreal(centro, centrof );
	
	double tmpsqrt;

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
				punto[0]=i;
				punto[1]=j;
				punto[2]=k;
				
				if( radioactual<distancia(punto, centro) ){
					*ptr=0;						
				}
				else{
					*ptr=255;
					
				}
			}
		}
	}

	for(t=0;t<cantidad;t++){
		radiof=sqrt(minis[t]);

		punto2[0]=candit[t][0];
		punto2[1]=candit[t][1];
		punto2[2]=candit[t][2];

		indextoreal(punto2, puntof2 );
	

		for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
			for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
				for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
					ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
					punto[0]=i;
					punto[1]=j;
					punto[2]=k;
					indextoreal(punto, puntof );

					dist=distanciaejepunto(puntof, centrof, puntof2);
					prop=proporcioejepunto(puntof, centrof, puntof2);

					if(prop>=0 && prop<=1){
						rest=(radiof-radioactual)*prop+radioactual;
						rest=rest+0.25;
						rest=rest*espprin[0];
					}
					else{
						rest=dist-1;
					}
					tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
					if( radio>sqrt(tmpsqrt) && rest>dist){
						*ptr=255;
					}					
				}
			}
		}
	}
}




//----------------------------------------------------------------------------
void axisExtractor02::comparacion(vtkImageData *data, vtkImageData *data2, unsigned long minis[4])
{
	int ext[6];
	int i, j, k;
	double radio;
	
	data->GetExtent(ext);

	unsigned char *ptr;
	unsigned char *ptr2;

	
	minis[0]=0;
	minis[1]=0;
	minis[2]=0;
	minis[3]=0;

	int centro[3];
	int punto[3];
	
	radio=((double)ext[1]-(double)ext[0])/2;
	
	centro[0]=(ext[1]+ext[0])/2;
	centro[1]=(ext[3]+ext[2])/2;
	centro[2]=(ext[5]+ext[4])/2;
	
	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
				ptr2=(unsigned char *)	data2->GetScalarPointer(i,j,k);

				punto[0]=i;
				punto[1]=j;
				punto[2]=k;
				
				if( radio>distancia(punto, centro) ){
									
					if( *ptr==0 && *ptr2==0 ){
						minis[0]++;						
					}
					else if( *ptr!=0 && *ptr2!=0 ){
						minis[1]++;						
					}
					else if(*ptr==0 && *ptr2!=0){
						minis[2]++;						
					}
					else{
						minis[3]++;						
					}
				}
			}
		}
	}
}



//----------------------------------------------------------------------------
void axisExtractor02::copiar(vtkImageData *data, vtkImageData *data2 )
{
	int ext[6];
	int i, j, k;
	
	data->GetExtent(ext);

	unsigned char *ptr, *ptr2;

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(unsigned char *)	data->GetScalarPointer(i,j,k);
				ptr2=(unsigned char *)	data2->GetScalarPointer(i,j,k);
				if(*ptr!=0 && *ptr2==0){
					*ptr2=*ptr;
				}
			}
		}
	}

	data2->Modified();
}




//----------------------------------------------------------------------------
void axisExtractor02::copiar2(vtkImageData *data, vtkImageData *data2 )
{
	int ext[6];
	int i, j, k;
	
	data->GetExtent(ext);

	double *ptr, *ptr2;

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(double *)	data->GetScalarPointer(i,j,k);
				ptr2=(double *)	data2->GetScalarPointer(i,j,k);
				if(*ptr!=0 && *ptr2==0){
					*ptr2=*ptr;
				}
			}
		}
	}

	data2->Modified();
}





//----------------------------------------------------------------------------
double axisExtractor02::angulo(double a[3], double b[3] )
{
	double m1=sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]));
	double m2=sqrt((b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]));
	double res=(acos(((a[0]*b[0])+(a[1]*b[1])+(a[2]*b[2]))/(m1*m2))*180)/3.1415;

	if(res<0){
		res=-res;
	}

	if(res>=0 && res<90){
		res=res;
	}
	else if(res>=90 && res<180){
		res=180-res;
	}
	else if(res>=180 && res<270){
		res=res-180;
	}
	else{
		res=360-res;
	}
	
	return res;
}




//----------------------------------------------------------------------------
double axisExtractor02::angulo(double i1, double j1, double k1, double i2, double j2, double k2 )
{
	double m1=sqrt((i1*i1)+(j1*j1)+(k1*k1));
	double m2=sqrt((i2*i2)+(j2*j2)+(k2*k2));
	double res=(acos(((i1*i2)+(j1*j2)+(k1*k2))/(m1*m2))*180)/3.1415;

	if(res<0){
		res=-res;
	}

	if(res>=0 && res<90){
		res=res;
	}
	else if(res>=90 && res<180){
		res=180-res;
	}
	else if(res>=180 && res<270){
		res=res-180;
	}
	else{
		res=360-res;
	}
	
	return res;
}





//----------------------------------------------------------------------------
int axisExtractor02::envolumen(int a[3], vtkImageData *datae )
{
	
	int res;

	unsigned short *ptr;

	ptr=(unsigned short *)	datae->GetScalarPointer(a);

	if(*ptr==0){
		res=0;
	}
	else{
		res=1;
	}
	
	return res;
}




//----------------------------------------------------------------------------
int axisExtractor02::mincandit(int candit[10][3], int cantidad, double puntoanterior[3])
{
	double distentpu;
	double distmientrp;
	int	inminenpu=0;

	double canditreal[3];
	int i;
	

	distmientrp=64000;
	inminenpu=0;

	for(i=0;i<cantidad && i<10;i++){
		indextoreal(candit[i], canditreal);
		distentpu=distancia(canditreal, puntoanterior);
		if(distentpu<distmientrp){
			distmientrp=distentpu;
			inminenpu=i;
		}
	}

	return inminenpu;

}



//----------------------------------------------------------------------------
int axisExtractor02::maxareacandit(unsigned long vector[50][4], int cantidad)
{
	
	unsigned long max;
	int i;
	int indmax=0;
	

	max=0;

	for(i=0;i<cantidad && i<10;i++){
		if(max<vector[i][0]){
			max=vector[i][0];
			indmax=i;
		}
	}

	return indmax;

}




//----------------------------------------------------------------------------
unsigned long axisExtractor02::totalarea(unsigned long vector[50][4], unsigned long vectorb[50][4], int cantidad, int cantidadb)
{
	
	unsigned long area;
	int i;
	
	area=0;

	for(i=0;i<cantidad && i<10;i++){
		area+=vector[i][0];
	}

	for(i=0;i<cantidadb;i++){
		area+=vectorb[i][0];
	}

	return area;
}





//----------------------------------------------------------------------------
unsigned long axisExtractor02::conecarea(unsigned long vector[50][4], int cantidad)
{
	
	unsigned long area;
	int i;
	
	area=0;

	for(i=0;i<cantidad && i<10;i++){
		area+=vector[i][0];
	}

	return area;
}






//----------------------------------------------------------------------------
int axisExtractor02::bruled(int candit[10][3], int cantidad, vtkImageData *data4)
{
	int i;
	int burned;
	
	burned=0;
	for(i=0;i<cantidad && i<10;i++){
		if(envolumen(candit[i], data4)!=0){
			burned++;
		}
	}

	return burned;
}




//----------------------------------------------------------------------------
double axisExtractor02::correction(int candit[10][3], int cantidad, vtkImageData *data, int indicecorregido[3], double puntocorregido[3], int indiceanterior[3], double radioanterior, int indicepre[3], double radiopre)
{
	
	
	int rap;
	int punto[3];
	double puntof[3];
	double puntof2[3];
	double centrof[3];
	double radiof;
	double maxcent;
	double mincent;
	double distcorr;
	double distcorr2;
	double distcorr3;
	double dist;
	double prop;
	double rest;
	double *ptr;
	int ext[6];
	int i, j, k, radio;
	int centro[3];
	
	data->GetExtent(ext);
	
	radio=(ext[1]-ext[0])/2;
	radiof=((double)ext[1]-(double)ext[0])/2;
	centro[0]=(ext[1]+ext[0])/2;
	centro[1]=(ext[3]+ext[2])/2;
	centro[2]=(ext[5]+ext[4])/2;

	indicecorregido[0]=centro[0];
	indicecorregido[1]=centro[1];
	indicecorregido[2]=centro[2];
	indextoreal(indicecorregido, puntocorregido);
	indextoreal(indiceanterior, centrof);
	indextoreal(indicepre, puntof2);
		
	rap= (int)( (double)radio/2 );
	maxcent=0;
	mincent=64000;

	for(i=ext[0];i<=ext[1];i++){
		for(j=ext[2];j<=ext[3];j++){
			for(k=ext[4];k<=ext[5];k++){
				ptr=(double *)	data->GetScalarPointer(i,j,k);
				punto[0]=i;
				punto[1]=j;
				punto[2]=k;
				indextoreal(punto, puntof );

				distcorr3=distancia(indiceanterior, indicepre);
				distcorr2=distancia(punto, indicepre);
				distcorr=distancia(punto, indiceanterior);

				dist=distanciaejepunto(puntof, centrof, puntof2);
				prop=proporcioejepunto(puntof, centrof, puntof2);

				if(prop>=0 && prop<=1){
					rest=(radiopre-radioanterior)*prop+radioanterior;
					rest=rest*espprin[0];
				}
				else{
					rest=dist-1;
				}
							
				if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && rest>dist){
				//if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && rest>dist){
				//if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && distcorr2<=distcorr3-radioanterior){
					if( maxcent<*ptr){
							mincent=distcorr;
							maxcent=*ptr;
							indicecorregido[0]=i;
							indicecorregido[1]=j;
							indicecorregido[2]=k;
							indextoreal(indicecorregido, puntocorregido);
					}
					else if( maxcent==*ptr){
						if( mincent>distcorr){
							mincent=distcorr;
							maxcent=*ptr;
							indicecorregido[0]=i;
							indicecorregido[1]=j;
							indicecorregido[2]=k;
							indextoreal(indicecorregido, puntocorregido);
						}
					}
				}
			}
		}
	}

	ptr=(double *)	data->GetScalarPointer(centro[0],centro[1],centro[2]);

	return sqrt(*ptr);
}




//----------------------------------------------------------------------------
double axisExtractor02::correction2(int candit[10][3], int cantidad, vtkImageData *data, int indicecorregido[3], double puntocorregido[3], int indiceanterior[3], double radioanterior)
{
	
	
	int rap;
	int punto[3];
	double maxcent;
	double mincent;
	double distcorr;
	double distcorr2;
	double *ptr;
	int ext[6];
	int i, j, k, radio;
	int centro[3];
	
	data->GetExtent(ext);
	
	radio=(ext[1]-ext[0])/2;
	centro[0]=(ext[1]+ext[0])/2;
	centro[1]=(ext[3]+ext[2])/2;
	centro[2]=(ext[5]+ext[4])/2;

	indicecorregido[0]=centro[0];
	indicecorregido[1]=centro[1];
	indicecorregido[2]=centro[2];
	indextoreal(indicecorregido, puntocorregido);
		
	rap= (int)((double)radio/2);
	maxcent=0;
	mincent=64000;

	for(i=centro[0]-rap;i<=centro[0]+rap;i++){
		for(j=centro[1]-rap;j<=centro[1]+rap;j++){
			for(k=centro[2]-rap;k<=centro[2]+rap;k++){
				ptr=(double *)	data->GetScalarPointer(i,j,k);
				punto[0]=i;
				punto[1]=j;
				punto[2]=k;

				distcorr2=distancia(punto, centro);
				distcorr=distancia(punto, indiceanterior);
				
				if( rap>distcorr2 ){
					if( maxcent<*ptr){
							mincent=distcorr;
							maxcent=*ptr;
							indicecorregido[0]=i;
							indicecorregido[1]=j;
							indicecorregido[2]=k;
							indextoreal(indicecorregido, puntocorregido);
					}
					else if( maxcent==*ptr){
						if( mincent>distcorr){
							mincent=distcorr;
							maxcent=*ptr;
							indicecorregido[0]=i;
							indicecorregido[1]=j;
							indicecorregido[2]=k;
							indextoreal(indicecorregido, puntocorregido);
						}
						
					}
				}
				
			}
		}
	}

	return sqrt(maxcent);
}





//----------------------------------------------------------------------------
void axisExtractor02::avanzar()
{
	
	double puntoactual[3];
	double puntocorregido[3];	
	double puntoanterior[3];	
	double puntosiguiente[3];
	//double puntointer[3];
	double puntopre[3];

	int indiceactual[3];
	int indiceanterior[3];
	int indicecorregido[3];
	//int indiceinter[3];
	int indicepre[3];
	double radioactual;
	double dirant[3];
		
	unsigned char cantidad;
	unsigned char cantidadb;

	unsigned long vector[50][4];
	unsigned long vectorb[50][4];

	int vectortemp[3];
	double tempc;
	
	int extint[6];
	
	int extintreal[6];
		
	int indexp;
	
	double provvc[3];
	
	int radiotrabajo;
	
	int flag =0;
	int flag2 =0;
	int flag3 =0;
	int flag4 =0;
	int flag5 =0;
	int flag6 =0;
	int flag7 =0;
	int flag8 =0;
	int flag9 =0;
	int flag10 =0;
	int flag11 =0;
	int flag12 =0;
	int flag13 =0;
	int flag14 =0;
	int flag15 =0;
	int flag16 =0;
	int flag17 =0;
	int flag18 =0;
	int flag19 =0;
	int flag20 =0;
	int flag21 =0;
	int flag22 =0;

	int i, j;
	
	double radiotrabajof;

	double proprad=2;
	
	double radiopred;
	double radioanterior;
	double radioprinc;

	int burned;
	
	int indant;
	double canditreal[3];
	int indmaxarea;
	unsigned long totalsurf;
	unsigned long conecsurf;

	indmaxarea=0;
	unsigned long caf[4];

	double rel1;
	double rel2;
	double rel3;
	double rel4;
	double rel5;
	double rel6;

	if(!m_Stack0.empty()){

	
	/*	indexp=m_Stack.top();

		puntoactual[0]=m_Stack0.top();
		puntoactual[1]=m_Stack1.top();
		puntoactual[2]=m_Stack2.top();

		puntoanterior[0]=m_Stack3.top();
		puntoanterior[1]=m_Stack4.top();
		puntoanterior[2]=m_Stack5.top();

		puntopre[0]=m_Stack6.top();
		puntopre[1]=m_Stack7.top();
		puntopre[2]=m_Stack8.top();

		

		radiopred=m_Stackr.top();
		radioanterior=m_Stackra.top();
		radioprinc=m_Stackrp.top();*/

		indexp=m_Stack.front();

		puntoactual[0]=m_Stack0.front();
		puntoactual[1]=m_Stack1.front();
		puntoactual[2]=m_Stack2.front();

		puntoanterior[0]=m_Stack3.front();
		puntoanterior[1]=m_Stack4.front();
		puntoanterior[2]=m_Stack5.front();

		puntopre[0]=m_Stack6.front();
		puntopre[1]=m_Stack7.front();
		puntopre[2]=m_Stack8.front();

		radiopred=m_Stackr.front();
		radioanterior=m_Stackra.front();
		radioprinc=m_Stackrp.front();
				
		radiotrabajo= (int) (proprad*radiopred);
		radiotrabajof=proprad*radiopred;

		if(radiotrabajof-radiotrabajo>0.5){
			radiotrabajo=radiotrabajo+1;
		}

// EED 15 Mai 2007 .NET 
//		radiotrabajo=max(minant,radiotrabajo);
		if (minant>radiotrabajo)
		{
			radiotrabajo=minant;
		}

		/*m_Stack0.pop();
		m_Stack1.pop();
		m_Stack2.pop();
		m_Stack3.pop();
		m_Stack4.pop();
		m_Stack5.pop();
		m_Stack6.pop();
		m_Stack7.pop();
		m_Stack8.pop();
		
		m_Stack.pop();
		m_Stackr.pop();
		m_Stackra.pop();
		m_Stackrp.pop();*/

		m_Stack0.pop_front();
		m_Stack1.pop_front();
		m_Stack2.pop_front();
		m_Stack3.pop_front();
		m_Stack4.pop_front();
		m_Stack5.pop_front();
		m_Stack6.pop_front();
		m_Stack7.pop_front();
		m_Stack8.pop_front();
		
		m_Stack.pop_front();
		m_Stackr.pop_front();
		m_Stackra.pop_front();
		m_Stackrp.pop_front();

		dirant[0]=puntoanterior[0]-puntoactual[0];
		dirant[1]=puntoanterior[1]-puntoactual[1];
		dirant[2]=puntoanterior[2]-puntoactual[2];

		
		realtoindex(puntoactual, indiceactual);
		realtoindex(puntoanterior, indiceanterior);
		realtoindex(puntopre, indicepre);

		radioactual=radiopred;
		
		if(radiotrabajo<=maxant){
			
			extint[0]=indiceactual[0]-radiotrabajo;
			extint[1]=indiceactual[0]+radiotrabajo;
			extint[2]=indiceactual[1]-radiotrabajo;
			extint[3]=indiceactual[1]+radiotrabajo;
			extint[4]=indiceactual[2]-radiotrabajo;
			extint[5]=indiceactual[2]+radiotrabajo;

			extrac->SetVOI(extint);
			extrac->UpdateWholeExtent();
			extrac->Update();
			extrac->GetOutput()->GetExtent(extintreal);

			if(extint[0]!=extintreal[0] || extint[1]!=extintreal[1] || extint[2]!=extintreal[2] || extint[3]!=extintreal[3] || extint[4]!=extintreal[4] || extint[5]!=extintreal[5]){
				fprintf(stream, "salio\t");
				flag4 =1;	
			}
			else{
				
				data1->Delete();
				data1=vtkImageData::New();
				data1->SetScalarTypeToUnsignedChar();
				data1->SetExtent(extrac->GetOutput()->GetExtent());
				data1->SetSpacing(extrac->GetOutput()->GetSpacing());

				cilindro(data1, dirant );
				
				optim(extrac->GetOutput(), data1 );

				connect->SetInput(data1);
				connect->RemoveAllSeeds();
				connect->AddSeed(indiceactual[0],indiceactual[1],indiceactual[2]);
				connect->UpdateWholeExtent();
				connect->Update();
							
				data2->Delete();
				data2=vtkImageData::New();
				data2->SetScalarTypeToUnsignedChar();
				data2->SetExtent(connect->GetOutput()->GetExtent());
				data2->SetSpacing(connect->GetOutput()->GetSpacing());

				blanquear(data2);
				
				cantidad=find_components(connect->GetOutput(), data2, 0, vector);

				data3->Delete();
				data3=vtkImageData::New();
				data3->SetScalarTypeToUnsignedChar();
				data3->SetExtent(connect->GetOutput()->GetExtent());
				data3->SetSpacing(connect->GetOutput()->GetSpacing());
				
				blanquear(data3);

				cantidadb=find_componentsb(connect->GetOutput(), data3, 0, vectorb);

				redondear3(connect->GetOutput());

				distance->Update();

				redondear2(distance->GetOutput());

				redondear(connect->GetOutput());
			
				costominimo(distance->GetOutput(), data2);
				
				
				if(buenos==0){
					radioactual=correction2(candit, cantidad, distance->GetOutput(),  indicecorregido,  puntocorregido, indiceanterior, radioanterior);
				}
				else{
					radioactual=correction(candit, cantidad, distance->GetOutput(),  indicecorregido,  puntocorregido, indiceanterior, radioanterior, indicepre, radioprinc);
				}


				//inpeq(extrac->GetOutput(), indicecorregido, radioactual);
								
				indant=mincandit(candit, cantidad, puntoanterior);

				indmaxarea=maxareacandit(vector, cantidad);
				totalsurf=totalarea(vector, vectorb, cantidad, cantidadb);
				conecsurf=conecarea(vector, cantidad);
				
				indextoreal(candit[indant], canditreal);

				burned=bruled(candit, cantidad, data4);
								
				data6->Delete();
				data6=vtkImageData::New();
				data6->SetScalarTypeToUnsignedChar();
				data6->SetExtent(extrac->GetOutput()->GetExtent());
				data6->SetSpacing(extrac->GetOutput()->GetSpacing());

				modelo(data6, cantidad,  vector, candit, radioactual, minis);
				comparacion(connect->GetOutput(), data6, caf);

				rel1=((double)caf[0]+(double)caf[1])/((double)caf[0]+(double)caf[1]+(double)caf[2]+(double)caf[3]);
				rel2=((double)caf[1])/((double)caf[1]+(double)caf[3]);
				rel3=((double)caf[1])/((double)caf[1]+(double)caf[2]);
				rel4=(double)conecsurf/(double)totalsurf;
				rel5=(double)vector[indmaxarea][0]/(double)totalsurf;
				rel6=(double)cant/((double)cant+(double)cant2);
										
					
				if(rel1<0.5 || rel2<0.5 || rel3<0.5 || (rel1<0.75 && rel2<0.75) || (rel2<0.75 && rel3<0.75)|| (rel1<0.75 && rel3<0.75)){
					flag5=1;
				}									
				
				if(indicecorregido[0]!=indiceactual[0] || indicecorregido[1]!=indiceactual[1] || indicecorregido[2]!=indiceactual[2]){
					flag7=1;
				}

				if(burned>=cantidad ){
					flag15=1;
				}

				if(burned==0 ){
					flag16=1;
				}

				if(cantidad<2 ){
					flag17=1;
				}

				if(flag7==1 && flagg>4){
					flag7=0;
				}

				if(buenos==0){
					flag19=1;
				}

				if(mejordst!=0){
					flag18=1;
				}

				if(flag7==1 && flag18==1){
					for(i=0;i<flagg;i++){
						if(visited[i][0]==indiceactual[0] && visited[i][1]==indiceactual[1] && visited[i][2]==indiceactual[2] && visitedrad[i]==radiopred){
							flag2=1;
						}
					}
				}
				
				if((flag19==1 && flagg2<5) || (flag19==1 && flagg2<10 && ( rel4>0.25 || rel5>0.13 || rel6>0.4 || flag17==1  || flag15==1 || flag5==1) ) || (flagg2<10 && ( ((rel4>0.25 || rel5>0.25 || rel6>0.4) && cantidad<3) || flag17==1  || flag15==1 || flag5==1) )){
						flag8=1;
				}

				if(centi<centi2){
					flag3=1;
				}

				if(angulo(dirant[0], dirant[1], dirant[2], canditreal[0]-puntoactual[0], canditreal[1]-puntoactual[1], canditreal[2]-puntoactual[2] )>65){
					flag11=1;
				}

				if(angulo(Vp[0][ejeminp], Vp[1][ejeminp], Vp[2][ejeminp], Vi[0][ejemini], Vi[1][ejemini], Vi[2][ejemini] )>60){
					flag6=1;
				}

				if((flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0)  ) && (flag7==0 && flag8==0)){
					flag9=1;
					//fprintf(stream, "malo\t");
				}

				if(flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0)  ){
					flag20=1;
				}
			
				if((mejordst<radioactual || (mejordst==radioactual && mejorrad>radiopred) || (mejordst==radioactual && mejorrad<=radiopred && mejorcant<cantidad )) && flag17==0 && flag20==0 && flag4==0 && flag12==0 && flag15==0 ){
					mejordst=radioactual;
					mejor[0]=puntoactual[0];
					mejor[1]=puntoactual[1];
					mejor[2]=puntoactual[2];
					mejorrad=radiopred;
					mejorcant=cantidad;
				}
			}
		}
		else{
			fprintf(stream, "grande\t");
			flag12 =1;	
		}

		if(mejordst!=0){
			flag18=1;
		}

		if(flagg>4){
			flag13=1;
		}

		if(flagg2>4){
			flag14=1;
		}

		if( flag9==1 || flag4==1 || flag12==1 || ((cantidad<2 || burned==cantidad) && (flag7==0 && flag8==0)) ){
			frama=1;
		}
		else{
			frama=0;
		}

		if( cantidad>2 && (flag7==0 && flag8==0) ){
			fseg=1;
		}
		else{
			fseg=0;
		}
		
		if((flag7==0 && flag8==0 && (mejordst>radioactual || (mejordst==radioactual && mejorrad<radiopred)) && flag19==0) || flag2==1){
			
			flagg=10;
			flagg2=10;
			mejordst=0;
			
			fprintf(stream, "mejor\t");
															
			m_Stack0.push_front(mejor[0]);
			m_Stack1.push_front(mejor[1]);
			m_Stack2.push_front(mejor[2]);
			m_Stack3.push_front(puntoanterior[0]);
			m_Stack4.push_front(puntoanterior[1]);
			m_Stack5.push_front(puntoanterior[2]);
			m_Stack6.push_front(puntopre[0]);
			m_Stack7.push_front(puntopre[1]);
			m_Stack8.push_front(puntopre[2]);
			m_Stack.push_front(indexp);
			m_Stackr.push_front(mejorrad);
			m_Stackra.push_front(radioanterior);
			m_Stackrp.push_front(radioprinc);

			mejorrad=0;
			mejorcant=0;

		}
		
		else if((flag17==1 || flag9==1 || flag4==1 || flag12==1 || flag15==1) && flag18==1){
			flagg=10;
			flagg2=10;
			mejordst=0;
			
			fprintf(stream, "descorrigio\t");
															
			m_Stack0.push_front(mejor[0]);
			m_Stack1.push_front(mejor[1]);
			m_Stack2.push_front(mejor[2]);
			m_Stack3.push_front(puntoanterior[0]);
			m_Stack4.push_front(puntoanterior[1]);
			m_Stack5.push_front(puntoanterior[2]);
			m_Stack6.push_front(puntopre[0]);
			m_Stack7.push_front(puntopre[1]);
			m_Stack8.push_front(puntopre[2]);
			m_Stack.push_front(indexp);
			m_Stackr.push_front(mejorrad);
			m_Stackra.push_front(radioanterior);
			m_Stackrp.push_front(radioprinc);

			mejorrad=0;
			mejorcant=0;

		}

		else if(flag12==0 && flag4==0 && flag9==0){
			if(flag7==0 && flag8==0){

				if(flag17==0 && flag15==0){
					flagg=0;
					flagg2=0;
					mejordst=0;
					mejorrad=0;
					mejorcant=0;

					fprintf(stream, "inserto\t");
				
					if(flag19==0){
						/*costominimo2(distance->GetOutput(), connect->GetOutput(), indiceactual, candit[indant], indiceinter);
						indextoreal(indiceinter, puntointer);
						realtoreal2(puntointer, provvc);
						points->InsertPoint(buenos,provvc);

						lineas->InsertNextCell(2);
						lineas->InsertCellPoint(indexp);
						lineas->InsertCellPoint(buenos);
						
						buenos++;

						realtoreal2(puntoactual, provvc);
						points->InsertPoint(buenos,provvc);
			
						lineas->InsertNextCell(2);
						lineas->InsertCellPoint(buenos-1);
						lineas->InsertCellPoint(buenos);
						
						buenos++;*/

						realtoreal2(puntoactual, provvc);
						points->InsertPoint(buenos,provvc);
			
						lineas->InsertNextCell(2);
						lineas->InsertCellPoint(indexp);
						lineas->InsertCellPoint(buenos);

						buenos++;
					
					}
					else{
						realtoreal2(puntoactual, provvc);
						points->InsertPoint(buenos,provvc);
										
						buenos++;
					}

					
				}
				else{
					flagg=0;
					flagg2=0;
					mejordst=0;
					mejorrad=0;
					mejorcant=0;
					fprintf(stream, "para\t");
				}
			}
			else if(flag7==0 && flag8==1){
				flagg2++;
				
				fprintf(stream, "agrando\t");

				m_Stack0.push_front(puntoactual[0]);
				m_Stack1.push_front(puntoactual[1]);
				m_Stack2.push_front(puntoactual[2]);
				m_Stack3.push_front(puntoanterior[0]);
				m_Stack4.push_front(puntoanterior[1]);
				m_Stack5.push_front(puntoanterior[2]);
				m_Stack6.push_front(puntopre[0]);
				m_Stack7.push_front(puntopre[1]);
				m_Stack8.push_front(puntopre[2]);

				m_Stack.push_front(indexp);
				m_Stackr.push_front(((double)radiotrabajo+1)/2);
				m_Stackra.push_front(radioanterior);
				m_Stackrp.push_front(radioprinc);
		
			}
			else if(flag7==1 && flag8==1){
				visited[flagg][0]=indiceactual[0];
				visited[flagg][1]=indiceactual[1];
				visited[flagg][2]=indiceactual[2];
				visitedrad[flagg]=radiopred;
				flagg++;
				flagg2++;
				
				fprintf(stream, "agrando y corrigiendo\t");

				m_Stack0.push_front(puntocorregido[0]);
				m_Stack1.push_front(puntocorregido[1]);
				m_Stack2.push_front(puntocorregido[2]);
				m_Stack3.push_front(puntoanterior[0]);
				m_Stack4.push_front(puntoanterior[1]);
				m_Stack5.push_front(puntoanterior[2]);
				m_Stack6.push_front(puntopre[0]);
				m_Stack7.push_front(puntopre[1]);
				m_Stack8.push_front(puntopre[2]);

				m_Stack.push_front(indexp);
				m_Stackr.push_front(((double)radiotrabajo+1)/2);
				m_Stackra.push_front(radioanterior);
				m_Stackrp.push_front(radioprinc);
			}
			else{
				visited[flagg][0]=indiceactual[0];
				visited[flagg][1]=indiceactual[1];
				visited[flagg][2]=indiceactual[2];
				visitedrad[flagg]=radiopred;
				flagg++;
				
				if(flag19==1){
					flagg=10;
				}

				fprintf(stream, "corrigio\t");
															
				m_Stack0.push_front(puntocorregido[0]);
				m_Stack1.push_front(puntocorregido[1]);
				m_Stack2.push_front(puntocorregido[2]);
				m_Stack3.push_front(puntoanterior[0]);
				m_Stack4.push_front(puntoanterior[1]);
				m_Stack5.push_front(puntoanterior[2]);

				m_Stack6.push_front(puntopre[0]);
				m_Stack7.push_front(puntopre[1]);
				m_Stack8.push_front(puntopre[2]);

				m_Stack.push_front(indexp);
				m_Stackr.push_front(radiopred);
				m_Stackra.push_front(radioanterior);
				m_Stackrp.push_front(radioprinc);

				
			}
				
			if(flag7==0 && flag8==0){
					
				if(flag17==0 && flag15==0){
				
					for(i=0;i<cantidad && i<10;i++){
						for(j=i;j<cantidad && i<10;j++){
							if(minis[j]>minis[i]){
								vectortemp[0]=candit[i][0];
								vectortemp[1]=candit[i][1];
								vectortemp[2]=candit[i][2];
								tempc=minis[i];
								
								candit[i][0]=candit[j][0];
								candit[i][1]=candit[j][1];
								candit[i][2]=candit[j][2];
								minis[i]=minis[j];
								
								candit[j][0]=vectortemp[0];
								candit[j][1]=vectortemp[1];
								candit[j][2]=vectortemp[2];
								minis[j]=tempc;
							}
						}
					}
					
	
					if(flag16==1){
						for(i=0;i<cantidad && i<10;i++){
							if(i!=indant || buenos<2){
								indextoreal(candit[i], puntosiguiente);
										
								m_Stack0.push_front(puntosiguiente[0]);
								m_Stack1.push_front(puntosiguiente[1]);
								m_Stack2.push_front(puntosiguiente[2]);
								m_Stack3.push_front(puntoactual[0]);
								m_Stack4.push_front(puntoactual[1]);
								m_Stack5.push_front(puntoactual[2]);
								
								m_Stack6.push_front(puntosiguiente[0]);
								m_Stack7.push_front(puntosiguiente[1]);
								m_Stack8.push_front(puntosiguiente[2]);
				
								double prov1a=sqrt(minis[i]);
								
								m_Stack.push_front(buenos-1);
								m_Stackr.push_front(prov1a);
								m_Stackra.push_front(radioactual);
								m_Stackrp.push_front(prov1a);
								
							}
						}
					}
					else{
						for(i=0;i<cantidad && i<10;i++){
							indextoreal(candit[i], puntosiguiente);
				
							if(envolumen(candit[i], data4)==0){
								m_Stack0.push_front(puntosiguiente[0]);
								m_Stack1.push_front(puntosiguiente[1]);
								m_Stack2.push_front(puntosiguiente[2]);
								m_Stack3.push_front(puntoactual[0]);
								m_Stack4.push_front(puntoactual[1]);
								m_Stack5.push_front(puntoactual[2]);
								
								m_Stack6.push_front(puntosiguiente[0]);
								m_Stack7.push_front(puntosiguiente[1]);
								m_Stack8.push_front(puntosiguiente[2]);
							
								double prov1b=sqrt(minis[i]);
								
								m_Stack.push_front(buenos-1);
								m_Stackr.push_front(prov1b);
								m_Stackra.push_front(radioactual);
								m_Stackrp.push_front(prov1b);
							}
						}
					}
				
					copiar(connect->GetOutput(),  data4 );
					//copiar2(distance->GetOutput(),  data5 );
				}						
			}
		}
		else{
		//	fprintf(stream, "algo");
			mejordst=0;
			mejorrad=0;
			mejorcant=0;
			flagg=0;
			flagg2=0;
				
		}

		
		
		
		/*
		printf("wi[0]: %f\n", wi[0]);
		printf("wi[1]: %f\n", wi[1]);
		printf("wi[2]: %f\n", wi[2]);
		
		printf("inerciari: %f\n", inerciari);
		printf("inerciariy: %f\n", inerciariy);
		printf("inerciariz: %f\n", inerciariz);

		printf("ri1: %f\n", inerciari/inerciariy);
		printf("ri2: %f\n", inerciariy/inerciariz);
		printf("ri3: %f\n", inerciari/inerciariz);

		printf("inerciarp: %f\n", inerciarp);
		printf("inerciarpy: %f\n", inerciarpy);
		printf("inerciarpz: %f\n", inerciarpz);

		printf("rp1: %f\n", inerciarp/inerciarpy);
		printf("rp2: %f\n", inerciarpy/inerciarpz);
		printf("rp3: %f\n", inerciarp/inerciarpz);

		printf("angulo: %f\n",angulo(dirant[0], dirant[1], dirant[2], canditreal[0]-puntoactual[0], canditreal[1]-puntoactual[1], canditreal[2]-puntoactual[2] ));
		printf("angulo2: %f\n", angulo(Vp[0][ejeminp], Vp[1][ejeminp], Vp[2][ejeminp], Vi[0][ejemini], Vi[1][ejemini], Vi[2][ejemini]));

		printf("centi: %f\n", centi);
		printf("centi2: %f\n", centi2);
		*/
		


		fprintf(stream, "%d\t", indexp);
		fprintf(stream, "%d\t", buenos);
		fprintf(stream, "%d\t", radiotrabajo);
		fprintf(stream, "%f\t", radioactual);
		fprintf(stream, "%f\t", radiopred);
		fprintf(stream, "%f\t", radioanterior);
		fprintf(stream, "%f\t", mejordst);
		fprintf(stream, "%f\t", mejorrad);	
		fprintf(stream, "%d\t", mejorcant);	
		
	
		/*
		fprintf(stream, "%d\t", cant);
		fprintf(stream, "%d\t", cant2);
		fprintf(stream, "%f\t", (double)cant/((double)cant+(double)cant2));
		*/	
		
		fprintf(stream, "%d\t", flag2);
		fprintf(stream, "%d\t", flag3);
		fprintf(stream, "%d\t", flag4);
		fprintf(stream, "%d\t", flag5);
		fprintf(stream, "%d\t", flag6);
		fprintf(stream, "%d\t", flag7);
		fprintf(stream, "%d\t", flag8);
		fprintf(stream, "%d\t", flag9);
		fprintf(stream, "%d\t", flag10);
		fprintf(stream, "%d\t", flag11);
		fprintf(stream, "%d\t", flag12);
		fprintf(stream, "%d\t", flag13);
		fprintf(stream, "%d\t", flag14);
		fprintf(stream, "%d\t", flag15);
		fprintf(stream, "%d\t", flag16);
		fprintf(stream, "%d\t", flag17);
		fprintf(stream, "%d\t", flag18);
		fprintf(stream, "%d\t", flag19);
		fprintf(stream, "%d\t", flag20);
		fprintf(stream, "%d\t", flagg);
		fprintf(stream, "%d\t", flagg2);
		fprintf(stream, "%d\t", cantidad);
		fprintf(stream, "%d\t", cantidadb);
		fprintf(stream, "%d\t", burned);

	
		/*
		fprintf(stream, "caf[0] %d\t", caf[0]);
		fprintf(stream, "caf[1] %d\t", caf[1]);
		fprintf(stream, "caf[2] %d\t", caf[2]);
		fprintf(stream, "caf[3] %d\t", caf[3]);
		fprintf(stream, "REL%f\t", ((double)caf[0]+(double)caf[1])/((double)caf[0]+(double)caf[1]+(double)caf[2]+(double)caf[3]));
		fprintf(stream, "REL2%f\t", ((double)caf[1])/((double)caf[3]+(double)caf[1]));
		fprintf(stream, "REL3%f\t", ((double)caf[1])/((double)caf[1]+(double)caf[2]));

		fprintf(stream, "maxareacandit%d\t", vector[indmaxarea][0]);
		fprintf(stream, "totalsurf%d\t", totalsurf);
		fprintf(stream, "conecsurf%d\t", conecsurf);
		fprintf(stream, "ratio%f\t", (double)conecsurf/(double)totalsurf);
		fprintf(stream, "ratio2%f\t", (double)vector[indmaxarea][0]/(double)totalsurf);
		*/
			

		fprintf(stream, "[%f, %f, %f]\t", puntoactual[0], puntoactual[1], puntoactual[2]);
		fprintf(stream, "[%f, %f, %f]\t", puntoanterior[0], puntoanterior[1], puntoanterior[2]);
		fprintf(stream, "[%f, %f, %f]\t", puntopre[0], puntopre[1], puntopre[2]);
		fprintf(stream, "[%f, %f, %f]\t", mejor[0], mejor[1], mejor[2]);

		fprintf(stream, "\n");

	}
	
}




//----------------------------------------------------------------------------
void axisExtractor02::todo()
{
	int i=0;

	while( !m_Stack0.empty() && i<5000){
		avanzar();
		i++;
	}
}



//----------------------------------------------------------------------------
void axisExtractor02::paso()
{

	int i=0;
	int inicio=buenos;

	while( !m_Stack0.empty() && inicio==buenos && i<5000){
		avanzar();
		i++;
	}
}




//----------------------------------------------------------------------------
void axisExtractor02::rama()
{

	int i=0;
	frama=0;
	
	while( !m_Stack0.empty() && frama==0 && i<5000){
		avanzar();
		i++;
	}
}




//----------------------------------------------------------------------------
void axisExtractor02::segmento()
{

	int i=0;
	fseg=0;
	frama=0;
	
	while( !m_Stack0.empty() && frama==0 && fseg==0 && i<5000){
		avanzar();
		i++;
	}
}