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


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


#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 "vtkMath.h"


#include "axisExtractor.h"

#define min(a,b) (((a)<(b))?(a):(b)) 


vtkStandardNewMacro(axisExtractor);


axisExtractor::axisExtractor() {
	
	this->NumberOfRequiredInputs = 1;

	this->humbral=0.45 ;
	this->maxpropradio=100;
	this->maxpropmasa=10;
	this->minpropmasa=0;
	
	this->resample= vtkImageResample::New();
	this->resample->SetDimensionality (3);
	
	this->extrac = vtkExtractVOI::New();
	this->extrac->SetInput(this->resample->GetOutput());
	this->extrac->SetSampleRate(1, 1, 1);
	
	this->thresh = vtkImageThreshold::New();
	this->thresh->SetInput(this->extrac->GetOutput());
	this->thresh->SetInValue(255);
	this->thresh->SetOutValue(0);
	this->thresh->ReleaseDataFlagOff();

	this->cast = vtkImageCast::New();
	this->cast->SetInput(this->thresh->GetOutput());
	this->cast->SetOutputScalarTypeToUnsignedChar();

	this->connect = vtkImageSeedConnectivity::New();
	this->connect->SetInput(this->cast->GetOutput());
	this->connect->SetInputConnectValue(255);
	this->connect->SetOutputConnectedValue(255);
	this->connect->SetOutputUnconnectedValue(0);

	this->dataprov=vtkImageData::New();
	this->dataprov->SetScalarTypeToChar();

	this->datatotal=vtkImageData::New();
	this->datatotal->SetScalarTypeToChar();

	
			

}

//----------------------------------------------------------------------------
// Specify the input data or filter.
void axisExtractor::SetInput(vtkImageData *input)
{
	double minspac;
	double espprin[3];
	
	
	this->vtkProcessObject::SetNthInput(0, input);
	this->resample->SetInput(input);
	input->GetSpacing(espprin);
	

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

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

	
}


//----------------------------------------------------------------------------
// Specify the input data or filter.

void axisExtractor::Execute()
{
	this->points = vtkPoints::New();
	this->lineas = vtkCellArray::New();
	this->buenos=0;
	this->iter=0;
	this->flagg=0;

	this->datatotal->Delete();
	this->datatotal=vtkImageData::New();
	this->datatotal->SetScalarTypeToUnsignedChar();

	this->datatotal->SetExtent(this->resample->GetOutput()->GetExtent());
	this->datatotal->SetSpacing(this->resample->GetOutput()->GetSpacing());

	this->blanquear2();

	
	while( !m_Stack0.empty()){
		this->avanzar();
		this->iter++;
	}


	this->GetOutput()->SetPoints (this->points);
	this->GetOutput()->SetLines(this->lineas);
	this->GetOutput()->Modified();
			
	// this->GetOutput()->GetPointData()->SetVectors(this->salidas);
}








//----------------------------------------------------------------------------
// Specify the input data or filter.
vtkImageData *axisExtractor::GetInput()
{
	if (this->NumberOfInputs < 1){
		return NULL;
    }
  
	return (vtkImageData *)(this->Inputs[0]);
}







//----------------------------------------------------------------------------
// Specify the input data or filter.
vtkImageData *axisExtractor::GetVolumen()
{

	return this->datatotal;
}






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








//----------------------------------------------------------------------------
void axisExtractor::SetMaxPropRadio(double value)
{
	this->maxpropradio=value;
}



//----------------------------------------------------------------------------
double axisExtractor::GetMaxPropRadio()
{
	return this->maxpropradio;
}


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

void axisExtractor::SetHumbral(double value)
{
	this->humbral=value;
}

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

double axisExtractor::GetHumbral()
{
	return this->humbral;
}


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

void axisExtractor::SetMaxPropMasa(double value)
{
	this->maxpropmasa=value;
}

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

double axisExtractor::GetMaxPropMasa()
{
	return this->maxpropmasa;
}


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

void axisExtractor::SetMinPropMasa(double value)
{
	this->minpropmasa=value;
}

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

double axisExtractor::GetMinPropMasa()
{
	return this->minpropmasa;
}







	//----------------------------------------------------------------------------
void axisExtractor::SetPoint(double puntoactualprov[3] )
{

	realtoreal(puntoactualprov,puntoactualprov);

		
	m_Stack0.push(puntoactualprov[0]);
	m_Stack1.push(puntoactualprov[1]);
	m_Stack2.push(puntoactualprov[2]);
	m_Stack3.push(puntoactualprov[0]);
	m_Stack4.push(puntoactualprov[1]);
	m_Stack5.push(puntoactualprov[2]);
	m_Stack.push(0);

}









//----------------------------------------------------------------------------
void axisExtractor::realtoreal(double a[3], double b[3] )
{
		
	double espprin[3];
	int extprin[6];				
	
	this->resample->GetOutput()->GetSpacing(espprin);
	this->resample->GetOutput()->GetExtent(extprin);

	b[0]=a[0]+(espprin[0]*extprin[0]);
	b[1]=a[1]+(espprin[1]*extprin[2]);
	b[2]=a[2];

	

}




//----------------------------------------------------------------------------
void axisExtractor::realtoreal2(double a[3], double b[3] )
{
		
	double espprin[3];
	int extprin[6];				
	
	this->resample->GetOutput()->GetSpacing(espprin);
	this->resample->GetOutput()->GetExtent(extprin);

	b[0]=a[0]-(espprin[0]*extprin[0]);
	b[1]=a[1]-(espprin[1]*extprin[2]);
	b[2]=a[2];

	

}





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

	double minspac;

	double espprin[3];
	int extprin[6];				
	
	this->resample->GetOutput()->GetSpacing(espprin);
	this->resample->GetOutput()->GetExtent(extprin);

	minspac=min(espprin[0],min(espprin[1],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 axisExtractor::indextoreal(int a[3], double b[3] )
{	
	double minspac;

	double espprin[3];
	int extprin[6];				
	
	this->resample->GetOutput()->GetSpacing(espprin);
	this->resample->GetOutput()->GetExtent(extprin);

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





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

	double espprin[3];
	int extprin[6];				
	
	this->resample->GetOutput()->GetSpacing(espprin);
	this->resample->GetOutput()->GetExtent(extprin);

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




//----------------------------------------------------------------------------
void axisExtractor::searc(int i, int j, int k)
{

	unsigned char *ptr;
	unsigned char	ptri;
	char *ptr2;
	
	int radio;

	int i2, j2, k2;
	int i3, j3, k3;
		
	ptr=(unsigned char *)	this->connect->GetOutput()->GetScalarPointer(i,j,k);
	ptr2=(char *)	this->dataprov->GetScalarPointer(i,j,k);

	ptri=*ptr;

	int ext[6];
	
	this->connect->GetOutput()->GetExtent(ext);

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

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


	if(ptri==255){
		*ptr2=this->label;
		this->vector[this->label-1][0]+=1;
		this->vector[this->label-1][1]+=i;
		this->vector[this->label-1][2]+=j;
		this->vector[this->label-1][3]+=k;
	}
	else if(ptri==0){
		*ptr2=-this->label2;
		this->vectorb[this->label2-1][0]+=1;
		this->vectorb[this->label2-1][1]+=i;
		this->vectorb[this->label2-1][2]+=j;
		this->vectorb[this->label2-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 *)	this->connect->GetOutput()->GetScalarPointer(i+i3,j+j3,k+k3);
						ptr2=(char *)	this->dataprov->GetScalarPointer(i+i3,j+j3,k+k3);
						if(ptri==255 && *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))){
							this->searc(i+i3, j+j3, k+k3);
						}
						else if(ptri==0 && *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))){
							this->searc(i+i3, j+j3, k+k3);
						}
					}
				}
			}
		}


	

	
}


//----------------------------------------------------------------------------
void axisExtractor::find_components( )
{
	int ext[6];
	int i, j, k, radio;
	int i2, j2, k2;

	
	unsigned char *ptr;
	char *ptr2;

	this->connect->GetOutput()->GetExtent(ext);


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

	this->label=0;
	this->label2=0;
		


	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 *)	this->connect->GetOutput()->GetScalarPointer(i,j,k);
					ptr2=(char *)	this->dataprov->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)  ){
						this->label=this->label+1;
						this->vector[this->label-1][0]=0;
						this->vector[this->label-1][1]=0;
						this->vector[this->label-1][2]=0;
						this->vector[this->label-1][3]=0;
						this->searc(i, j, k);
						
					}
					else if(*ptr==0 && *ptr2==0 && (radio*radio)>=(i2*i2)+(j2*j2)+(k2*k2) && ((radio-1)*(radio-1))<=(i2*i2)+(j2*j2)+(k2*k2)  ){
						this->label2=this->label2+1;
						this->vectorb[this->label2-1][0]=0;
						this->vectorb[this->label2-1][1]=0;
						this->vectorb[this->label2-1][2]=0;
						this->vectorb[this->label2-1][3]=0;
						this->searc(i, j, k);
						
					}
				}
			}
		}



	

	
	
}




//----------------------------------------------------------------------------
unsigned short axisExtractor::maximo( )
{
	int ext[6];
	int i, j, k;
	unsigned short max=0;

	this->extrac->GetOutput()->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 *)this->extrac->GetOutput()->GetScalarPointer(i,j,k);
					if(*ptr>max){
						max=*ptr;
						
					}
				}
			}
		}

		return max;


}




//----------------------------------------------------------------------------
void axisExtractor::blanquear()
{
	int ext[6];
	int i, j, k;
	
	this->dataprov->GetExtent(ext);

	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=(char *)	this->dataprov->GetScalarPointer(i,j,k);
					*ptr=0;
				}
			}
		}
}






//----------------------------------------------------------------------------
void axisExtractor::blanquear2()
{
	int ext[6];
	int i, j, k;
	
	this->datatotal->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 *)	this->datatotal->GetScalarPointer(i,j,k);
					*ptr=0;
				}
			}
		}
}




//----------------------------------------------------------------------------
void axisExtractor::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();
}





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








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


	unsigned char *ptr;

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

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







//----------------------------------------------------------------------------
void axisExtractor::corte(double punto1[3], double punto2[3], double punto3[3], double centro[3],  double radio )
{
	
	double m1, m2, m3;
	double b1, b2, b3;
	double c0, c1, c2;
	double* roots;
	double root;
	

	m1=punto2[0]-punto1[0];
	m2=punto2[1]-punto1[1];
	m3=punto2[2]-punto1[2];

	b1=punto1[0]-centro[0];
	b2=punto1[1]-centro[1];
	b3=punto1[2]-centro[2];

	c0=m1*m1+m2*m2+m3*m3;
	c1=2*m1*b1+2*m2*b2+2*m3*b3;
	c2=b1*b1+b2*b2+b3*b3-radio*radio;

	roots=vtkMath::SolveQuadratic  (c0, c1, c2);

	if(roots[1]>=0 && roots[1]<=1){
		root=roots[1];
	}
	else{
		root=roots[2];
	}

	punto3[0]=punto1[0]+root*(punto2[0]-punto1[0]);
	punto3[1]=punto1[1]+root*(punto2[1]-punto1[1]);
	punto3[2]=punto1[2]+root*(punto2[2]-punto1[2]);
		
}












//----------------------------------------------------------------------------
void axisExtractor::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;
						
					}
				}
			}
		}
}













//----------------------------------------------------------------------------
void axisExtractor::avanzar()
{



	
	double puntoactual[3];
	double puntoanterior[3];
	double puntoactualdis[3];
	double puntoactualcorre[3];
	int puntoactualarr[3];
	double radioactual;
	double dirant[3];
	double vector2[50][3];
	double vector2b[3];
	unsigned long vectortemp[4];
	double vector2temp[3];
	double maxmasa;
	int extint[6];
	int extintreal[6];
	int indesxdis[3];
	int indexp;
	double norvec;
	double provvc[3];
	double proprov;
	double puntoactualcorre2[3]; 
	int radiotrabajo=1;
	int radiocorte=0;
	int radiocorte2=0;
	int radiocorte3=0;
	int radiocorte4=0;
	int radiocorte5=0;
	int radiocorte6=0;
	int radiocorte7=0;
	int radiocorte8=0;
	int radiocorte9=0;
	int radiocorte10=0;
	int radiocorte11=0;
	int cantidadanterior=0;
	int cantidadanteriorb=0;
	int flag =0;
	int flag2 =0;
	int flag3 =0;
	int flag4 =0;
	int flag5 =0;
	int flag6 =0;
	int flag7 =0;
	int i, j;

				
		
	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();

		m_Stack0.pop();
		m_Stack1.pop();
		m_Stack2.pop();
		m_Stack3.pop();
		m_Stack4.pop();
		m_Stack5.pop();
		m_Stack.pop();

		dirant[0]=puntoanterior[0]-puntoactual[0];
		dirant[1]=puntoanterior[1]-puntoactual[1];
		dirant[2]=puntoanterior[2]-puntoactual[2];
		
		radioactual=distancia(puntoactual, puntoanterior);

		realtoindex(puntoactual, puntoactualarr);
		
		for(this->label=1, this->label2=0;flag4==0  &&(this->label>0) && (this->label==1 && ( radiocorte==0 || radiotrabajo<=radiocorte4 ) ) ||  (this->label>1 && radiotrabajo<=radiocorte4) ;radiotrabajo++){
					
			extint[0]=puntoactualarr[0]-radiotrabajo;
			extint[1]=puntoactualarr[0]+radiotrabajo;
			extint[2]=puntoactualarr[1]-radiotrabajo;
			extint[3]=puntoactualarr[1]+radiotrabajo;
			extint[4]=puntoactualarr[2]-radiotrabajo;
			extint[5]=puntoactualarr[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]){
				flag4 =1;
				radiotrabajo++;
				break;
			}
			
			if(puntoactualarr[0]>=extintreal[0] && puntoactualarr[0]<=extintreal[1] && puntoactualarr[1]>=extintreal[2] && puntoactualarr[1]<=extintreal[3] && puntoactualarr[2]>=extintreal[4] && puntoactualarr[2]<=extintreal[5]  ){

				thresh->ThresholdByUpper(this->maximo()*this->humbral);
				thresh->UpdateWholeExtent();
				thresh->Update();

				redondear(thresh->GetOutput());

				cast->UpdateWholeExtent();
				cast->Update();

				connect->RemoveAllSeeds();
				connect->AddSeed(puntoactualarr[0],puntoactualarr[1],puntoactualarr[2]);
				connect->UpdateWholeExtent();
				connect->Update();

				this->dataprov->Delete();

				this->dataprov=vtkImageData::New();
				this->dataprov->SetScalarTypeToChar();

				this->dataprov->SetExtent(this->connect->GetOutput()->GetExtent());
				this->dataprov->SetSpacing(this->connect->GetOutput()->GetSpacing());
			
				this->blanquear();

				this->find_components();

				if(this->label2>0 && flag==0){
					radiocorte=radiotrabajo;
					vector2b[0]=((double)vectorb[0][1]/(double)vectorb[0][0])-(double)puntoactualarr[0];
					vector2b[1]=((double)vectorb[0][2]/(double)vectorb[0][0])-(double)puntoactualarr[1];
					vector2b[2]=((double)vectorb[0][3]/(double)vectorb[0][0])-(double)puntoactualarr[2];
					flag =1;
				}
				if(this->label2>1 && flag5==0){
					radiocorte5=radiotrabajo;
					flag5 =1;
				}

				if( this->label>1 && flag2==0){
					radiocorte2=radiotrabajo;
					flag2 =1;
				}

				if(radiocorte5==0){
					radiocorte6=radiocorte2;
				}
				else{
					radiocorte6=min(radiocorte2, radiocorte5);
				}

				if((radiocorte6-radiocorte)>1 && flagg<=4){
					flag4 =1;
					radiotrabajo++;
					flag7 =0;
					break;
				}
				else{
					flag7 =1;
				}
				
				if(this->label>2 && flag3==0){
					radiocorte3=radiotrabajo;
					flag3 =1;
				}
				if( this->label==2 && cantidadanterior!=this->label){
					radiocorte7=radiotrabajo;

				}
				if( this->label==3 && cantidadanterior!=this->label){
					radiocorte8=radiotrabajo;
				}
				if( cantidadanterior!=this->label || cantidadanteriorb!=this->label2){
					if(radiocorte2==0){
						radiocorte4=radiotrabajo+10*(radiotrabajo-radiocorte9);
					}
					else{
						radiocorte4=radiotrabajo+(radiotrabajo-radiocorte9);
					}
					radiocorte9=radiotrabajo;
				}
				
				cantidadanterior=this->label;
				cantidadanteriorb=this->label2;
				
			}
			
		}
		radiotrabajo--;
		radiocorte10=radiotrabajo;

		if(flag7==1){
			if(this->label>=2 && radiocorte10!=radiocorte9+1){
				radiotrabajo=radiocorte9+1;
				
				extint[0]=puntoactualarr[0]-radiotrabajo;
				extint[1]=puntoactualarr[0]+radiotrabajo;
				extint[2]=puntoactualarr[1]-radiotrabajo;
				extint[3]=puntoactualarr[1]+radiotrabajo;
				extint[4]=puntoactualarr[2]-radiotrabajo;
				extint[5]=puntoactualarr[2]+radiotrabajo;

				extrac->SetVOI(extint);
				extrac->UpdateWholeExtent();
				extrac->Update();
				extrac->GetOutput()->GetExtent(extintreal);
	
				thresh->ThresholdByUpper(this->maximo()*this->humbral);
				thresh->UpdateWholeExtent();
				thresh->Update();

				redondear(thresh->GetOutput());

				cast->UpdateWholeExtent();
				cast->Update();

				connect->RemoveAllSeeds();
				connect->AddSeed(puntoactualarr[0],puntoactualarr[1],puntoactualarr[2]);
				connect->UpdateWholeExtent();
				connect->Update();

				this->dataprov->Delete();

				this->dataprov=vtkImageData::New();
				this->dataprov->SetScalarTypeToChar();

				this->dataprov->SetExtent(this->connect->GetOutput()->GetExtent());
				this->dataprov->SetSpacing(this->connect->GetOutput()->GetSpacing());

			
				this->blanquear();

				this->find_components();
			}
		}

		if(flag7==1){
			flagg=0;
			if(this->label>1){
				//printf("no corrigio\n");
				
				realtoreal2(puntoactual, provvc);
				points->InsertPoint(buenos,provvc);
										
				if(buenos>0){
					lineas->InsertNextCell(2);
					lineas->InsertCellPoint(indexp);
					lineas->InsertCellPoint(buenos);
				}

				buenos++;
			}
		}
		else{
			flagg++;
			//printf("corrigio\n");
			
			provvc[0]=0;
			provvc[1]=0;
			provvc[2]=0;

			norvec=distancia(provvc, vector2b);

			proprov=(radiocorte6-radiocorte)/(2*norvec);
				
			puntoactualcorre[0]=(puntoactualarr[0]-(vector2b[0]*proprov));
			puntoactualcorre[1]=(puntoactualarr[1]-(vector2b[1]*proprov));
			puntoactualcorre[2]=(puntoactualarr[2]-(vector2b[2]*proprov));

			indextoreal(puntoactualcorre, puntoactualcorre2);
						
			m_Stack0.push(puntoactualcorre2[0]);
			m_Stack1.push(puntoactualcorre2[1]);
			m_Stack2.push(puntoactualcorre2[2]);
			m_Stack3.push(puntoanterior[0]);
			m_Stack4.push(puntoanterior[1]);
			m_Stack5.push(puntoanterior[2]);
			m_Stack.push(indexp);

			
		}

		if(flag7==1){
				
			if(this->label>1){
				
				maxmasa=0;
				for(i=0;i<this->label;i++){
					vector2[i][0]=((double)vector[i][1]/(double)vector[i][0])-(double)puntoactualarr[0];
					vector2[i][1]=((double)vector[i][2]/(double)vector[i][0])-(double)puntoactualarr[1];
					vector2[i][2]=((double)vector[i][3]/(double)vector[i][0])-(double)puntoactualarr[2];
					if(maxmasa<vector[i][0]){
						maxmasa=vector[i][0];
					}	
				}

				for(i=0;i<this->label;i++){
					for(j=i;j<this->label;j++){
						if(vector[j][0]<vector[i][0]){
							vectortemp[0]=vector[i][0];
							vectortemp[1]=vector[i][1];
							vectortemp[2]=vector[i][2];
							vectortemp[3]=vector[i][3];
							vector2temp[0]=vector2[i][0];
							vector2temp[1]=vector2[i][1];
							vector2temp[2]=vector2[i][2];
							
							vector[i][0]=vector[j][0];
							vector[i][1]=vector[j][1];
							vector[i][2]=vector[j][2];
							vector[i][3]=vector[j][3];
							vector2[i][0]=vector2[j][0];
							vector2[i][1]=vector2[j][1];
							vector2[i][2]=vector2[j][2];
						
							vector[j][0]=vectortemp[0];
							vector[j][1]=vectortemp[1];
							vector[j][2]=vectortemp[2];
							vector[j][3]=vectortemp[3];
							vector2[j][0]=vector2temp[0];
							vector2[j][1]=vector2temp[1];
							vector2[j][2]=vector2temp[2];
						}
					}
					
				}

				
				for(i=0;i<this->label;i++){
					if(maxmasa*this->minpropmasa<vector[i][0]){

						indesxdis[0]=(int) (puntoactualarr[0]+vector2[i][0]);
						indesxdis[1]=(int) (puntoactualarr[1]+vector2[i][1]);
						indesxdis[2]=(int) (puntoactualarr[2]+vector2[i][2]);
	
						indextoreal(indesxdis, puntoactualdis);
		
						if(envolumen(indesxdis, datatotal)==0){
							m_Stack0.push(puntoactualdis[0]);
							m_Stack1.push(puntoactualdis[1]);
							m_Stack2.push(puntoactualdis[2]);
							m_Stack3.push(puntoactual[0]);
							m_Stack4.push(puntoactual[1]);
							m_Stack5.push(puntoactual[2]);
							m_Stack.push(buenos-1);
						}
					
					}
					
				}

				if(this->label>1 ){
					copiar(this->connect->GetOutput(),  datatotal );
				}

			}
		}
		
	}
	

}