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


=========================================================================*/
#include "carotidaBifurcacion.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 "vtkMath.h"



vtkStandardNewMacro(carotidaBifurcacion);


carotidaBifurcacion::carotidaBifurcacion() {
	
	this->NumberOfRequiredInputs = 2;

	this->humbral=0.45 ;
	this->maxpropradio=5;
	this->maxpropmasa=10;
	this->minpropmasa=0.001;
	
	this->extrac = vtkExtractVOI::New();
	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->SetScalarTypeToUnsignedChar();

	
			

}
//----------------------------------------------------------------------------
// Specify the input data or filter.
void carotidaBifurcacion::SetInput(vtkPolyData *input)
{
	this->vtkProcessObject::SetNthInput(0, input);
}


//----------------------------------------------------------------------------
// Specify the input data or filter.
void carotidaBifurcacion::SetInput2(vtkImageData *input)
{
	this->vtkProcessObject::SetNthInput(1, input);
	this->extrac->SetInput(this->GetInput2());
  
}


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

void carotidaBifurcacion::Execute()
{
	vtkIdType totalpuntos=this->GetInput()->GetNumberOfPoints();
	this->points = vtkPoints::New();
	this->points2 = vtkPoints::New();
	this->points3 = vtkPoints::New();
	this->points4 = vtkPoints::New();
	this->buenos2=0;
	this->iter=0;
	/*this->lineas = vtkCellArray::New();
	this->salidas= vtkDoubleArray::New();
	this->salidas->SetNumberOfComponents(3);*/

	for(this->iter=0;this->iter<totalpuntos;this->iter++){
		this->avanzar(this->iter);
	}

	this->limpiar(points, points2, points3, points4 );
   
	this->GetOutput()->SetPoints (this->points);
	/*this->GetOutput()->SetLines(this->lineas);
	this->GetOutput()->GetPointData()->SetVectors(this->salidas);*/
}







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







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







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








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



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


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

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

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

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


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

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

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

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


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

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

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

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







void carotidaBifurcacion::searc(int i, int j, int k)
{

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

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

	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;


	*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;


	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=(unsigned char *)	this->dataprov->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-2)*(radio-2))<((i2+i3)*(i2+i3))+((j2+j3)*(j2+j3))+((k2+k3)*(k2+k3))){
							this->searc(i+i3, j+j3, k+k3);
						}
					}
				}
			}
		}
}



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

	
	unsigned char *ptr;
	unsigned char *ptr2;

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


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

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

					if(*ptr==255 && *ptr2==0 && (radio*radio)>(i2*i2)+(j2*j2)+(k2*k2) && ((radio-2)*(radio-2))<(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);
						
					}
				}
			}
		}
	
}





unsigned short carotidaBifurcacion::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 carotidaBifurcacion::blanquear()
{
	int ext[6];
	int i, j, k;
	
	this->dataprov->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->dataprov->GetScalarPointer(i,j,k);
					*ptr=0;
				}
			}
		}
}





double carotidaBifurcacion::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<180){
		res=res;
	}
	else{
		res=360-res;
	}
	
	return res;
}



void carotidaBifurcacion::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 carotidaBifurcacion::direcciones(vtkPolyData *profile, int iter, double radio, double puntocortea[3], double puntocortes[3] )
{
	double puntoactual[3];
	double puntoant[3];
	double puntosig[3];
	double puntoant2[3];
	double puntosig2[3];
	double dist;
	int i,j;
	
	vtkIdType totalpuntos=profile->GetNumberOfPoints();

	if(iter<totalpuntos){

		profile->GetPoint(iter, puntoactual);
	
		if(iter==0){
			profile->GetPoint(iter+1, puntosig);
			profile->GetPoint(iter, puntosig2);
			i=1;

			dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));

			while(dist< radio){
			
				profile->GetPoint(iter+1+i, puntosig);
				profile->GetPoint(iter+i, puntosig2);
				dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
				i++;

			}

			this->corte(puntosig, puntosig2, puntocortes, puntoactual,  radio );

			puntocortea[0]=2*puntoactual[0]-puntocortes[0];
			puntocortea[1]=2*puntoactual[1]-puntocortes[1];
			puntocortea[2]=2*puntoactual[2]-puntocortes[2];
			
		}
		else if(iter==totalpuntos-1){
			profile->GetPoint(iter-1, puntoant);
			profile->GetPoint(iter, puntoant2);
			j=1;

			dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));

			while(dist< radio){
			
				profile->GetPoint(iter-1-j, puntoant);
				profile->GetPoint(iter-j, puntoant2);
				dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
				j++;

			}

			this->corte(puntoant, puntoant2, puntocortea, puntoactual,  radio );

			puntocortes[0]=2*puntoactual[0]-puntocortea[0];
			puntocortes[1]=2*puntoactual[1]-puntocortea[1];
			puntocortes[2]=2*puntoactual[2]-puntocortea[2];
		}
		else{
			profile->GetPoint(iter+1, puntosig);
			profile->GetPoint(iter, puntosig2);
			i=1;

			dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));

			while(dist< radio && iter+1+i<totalpuntos){
			
				profile->GetPoint(iter+1+i, puntosig);
				profile->GetPoint(iter+i, puntosig2);
				dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntosig ));
				i++;

			}

			profile->GetPoint(iter-1, puntoant);
			profile->GetPoint(iter, puntoant2);
			j=1;

			dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));

			while(dist< radio && iter-1-j>=0){
			
				profile->GetPoint(iter-1-j, puntoant);
				profile->GetPoint(iter-j, puntoant2);
				dist=sqrt(vtkMath::Distance2BetweenPoints(puntoactual, puntoant ));
				j++;
			}


			if(iter+1+i<totalpuntos){
				this->corte(puntosig, puntosig2, puntocortes, puntoactual,  radio );				
			}

			if(iter-1-j>=0){
				this->corte(puntoant, puntoant2, puntocortea, puntoactual,  radio );	
			}

			if(iter+1+i>=totalpuntos){
				puntocortes[0]=2*puntoactual[0]-puntocortea[0];
				puntocortes[1]=2*puntoactual[1]-puntocortea[1];
				puntocortes[2]=2*puntoactual[2]-puntocortea[2];
			}
			
			if(iter-1-j<0){
				puntocortea[0]=2*puntoactual[0]-puntocortes[0];
				puntocortea[1]=2*puntoactual[1]-puntocortes[1];
				puntocortea[2]=2*puntoactual[2]-puntocortes[2];
			}
		}
	}
	
}

int carotidaBifurcacion::igual(double puntoactual[3], double puntoactualdis[3], double puntoantguar[3], double puntoantguardis[3] )
{
	int res=0;
	double mat[3][3];
	double matinv[3][3];
				
	double t[3];
	double b[3];

	double disact;
	double disant;
	double discru;

								

	disact=sqrt(vtkMath::Distance2BetweenPoints (puntoactual,puntoactualdis));
	disant=sqrt(vtkMath::Distance2BetweenPoints (puntoantguar,puntoantguardis));
	discru=sqrt(vtkMath::Distance2BetweenPoints (puntoactual,puntoantguar));

							
	if(discru<disact+disant){
		mat[0][0]=puntoantguardis[0]-puntoantguar[0];
		mat[0][1]=-(puntoactualdis[0]-puntoactual[0]);
		mat[0][2]=0;
										
		mat[1][0]=puntoantguardis[1]-puntoantguar[1];
		mat[1][1]=-(puntoactualdis[1]-puntoactual[1]);
		mat[1][2]=0;

		mat[2][0]=0;
		mat[2][1]=0;
		mat[2][2]=1;

		b[0]=puntoactual[0]-puntoantguar[0];
		b[1]=puntoactual[1]-puntoantguar[1];
		b[2]=0;

									
		vtkMath::Invert3x3(mat,matinv);
		vtkMath::Multiply3x3 (matinv, b, t);


								
		if(t[0]>0 && t[1]>0 ){
			res=1;
		}

									
	}

	return res;
}



int carotidaBifurcacion::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;
}





void carotidaBifurcacion::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 carotidaBifurcacion::limpiar(vtkPoints *profile, vtkPoints *profile2, vtkPoints *profile3, vtkPoints *profile4 )
{
	int i,j,k,l,m, n, p;
	double puntoactual[3];
	double puntoactualdis[3];
	double puntoantguar[3];
	double puntoantguardis[3];
	double promedio[3];
	int cant;
	int total[100][10];
	int pairs[100][2];
	int flag;
	
	vtkIdType totalpuntos=profile2->GetNumberOfPoints();

	k=0;

	for(i=0;i<totalpuntos;i++){
		flag=0;
		profile2->GetPoint(i, puntoactual);
		profile3->GetPoint(i, puntoactualdis);
			
		for(j=i+1;j<totalpuntos;j++){
			profile2->GetPoint(j, puntoantguar);
			profile3->GetPoint(j, puntoantguardis);
			if(puntoactual[0]!=puntoantguar[0] && puntoactual[1]!=puntoantguar[1] && puntoactual[2]!=puntoantguar[2] && 1==this->igual(puntoactual, puntoactualdis, puntoantguar, puntoantguardis)){
				pairs[k][0]=i;
				pairs[k][1]=j;
				k++;
				flag=1;
			}

		}
		for(n=0;n<k;n++){
			if(pairs[n][1]==i){
				flag=1;
			}

		}
		if(flag==0){
			pairs[k][0]=i;
			pairs[k][1]=-1;
			k++;
		}
	}

	for(i=0;i<10;i++){
		for(j=0;j<100;j++){
			total[j][i]=-1;
		}
	}


	l=0;

	for(i=0;i<k;i++){
		m=0;
		if(pairs[i][0]!=-1){
			total[l][m]=pairs[i][0];
			m++;
		}
		if(pairs[i][1]!=-1){
			total[l][m]=pairs[i][1];
			m++;
		}
		if(pairs[i][0]!=-1 || pairs[i][1]!=-1){
			pairs[i][0]=-1;
			pairs[i][1]=-1;
			l++;
			
		
			for(j=i+1;j<k;j++){

				for(n=0;n<m;n++){
				
					if(total[l-1][n]==pairs[j][0]){
						flag=0;
						for(p=0;p<m;p++){
							if(total[l-1][p]==pairs[j][1]){
								flag=1;
							}
						}
						if(flag==0){
							total[l-1][m]=pairs[j][1];
							m++;
						}
						pairs[j][0]=-1;
						pairs[j][1]=-1;
							

					}
					else if(total[l-1][n]==pairs[j][1]){
						flag=0;
						for(p=0;p<m;p++){
							if(total[l-1][p]==pairs[j][0]){
								flag=1;
							}
						}
						if(flag==0){
							total[l-1][m]=pairs[j][0];
							m++;
						}
						pairs[j][0]=-1;
						pairs[j][1]=-1;
					}
					
				}
			}
		}
	
	}

	for(i=0;i<l;i++){
		promedio[0]=0;
		promedio[1]=0;
		promedio[2]=0;
		cant=0;
		flag=0;
		m=0;

		while(flag==0){
			if(total[i][m]!=-1){
				profile4->GetPoint(total[i][m], puntoactual);
				promedio[0]+=puntoactual[0];
				promedio[1]+=puntoactual[1];
				promedio[2]+=puntoactual[2];
				cant++;
			}
			else{
				flag=1;
			}
			m++;
		}
		promedio[0]=promedio[0]/cant;
		promedio[1]=promedio[1]/cant;
		promedio[2]=promedio[2]/cant;
		
		profile->InsertPoint(i,promedio);


	}


}


void carotidaBifurcacion::avanzar(vtkIdType iter)
{


	double espprin[3];
	int extprin[6];
	double puntoactual[3];
	double puntoactualdis[3];
	int puntoactualarr[3];
	double radioactual;
	double dirant[3];
	double dirsig[3];
	
	
	
	double vector2[50][3];
	double angulos[50];

	unsigned long vectortemp[4];
	double vector2temp[3];
	double angulostemp;

	double maxmasa;
		
	int extint[6];
	int extintreal[6];
	int prop;

	int i, j;
	int radiotrabajo;

	
				
	
	this->GetInput2()->GetExtent(extprin);
	this->GetInput2()->GetSpacing(espprin);

	vtkIdType totalpuntos=this->GetInput()->GetNumberOfPoints();

	if(iter<totalpuntos){

		vtkDoubleArray *radio = (vtkDoubleArray *)this->GetInput()->GetPointData()->GetScalars("radio");
		//vtkDoubleArray *dir = (vtkDoubleArray *)this->GetInput()->GetPointData()->GetVectors("norm");
		this->GetInput()->GetPoint(iter, puntoactual);
		radioactual=radio->GetValue(iter);
		radiotrabajo=1;

		for(this->label=0, prop=100;radiotrabajo<radioactual*this->maxpropradio && (this->label<2 || prop>this->maxpropmasa);radiotrabajo++){

			puntoactualarr[0]= (int) ( (puntoactual[0]/espprin[0])+extprin[0] );
			puntoactualarr[1]= (int) ( (puntoactual[1]/espprin[1])+extprin[2] );
			puntoactualarr[2]= (int) ( (puntoactual[2]/espprin[2])            );

			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(puntoactualarr[0]>=extintreal[0] && puntoactualarr[0]<=extintreal[1] && puntoactualarr[1]>=extintreal[2] && puntoactualarr[1]<=extintreal[3] && puntoactualarr[2]>=extintreal[4] && puntoactualarr[2]<=extintreal[5]  ){

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

			prop=this->proporcion(this->thresh->GetOutput());
			this->redondear(this->thresh->GetOutput());


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


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

			this->dataprov->Delete();

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

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

			this->find_components();
			
				

				
				
			}
		}


		this->direcciones(this->GetInput(), this->iter, radiotrabajo-1, dirant, dirsig);

			

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




			
			
			maxmasa=0;
			if(this->label>2){
				for(i=0;i<this->label;i++){
					vector2[i][0]=((double)this->vector[i][1]/(double)this->vector[i][0])-(double)puntoactualarr[0];
					vector2[i][1]=((double)this->vector[i][2]/(double)this->vector[i][0])-(double)puntoactualarr[1];
					vector2[i][2]=((double)this->vector[i][3]/(double)this->vector[i][0])-(double)puntoactualarr[2];
					double tmp_EED_A = this->angulo(dirsig[0], dirsig[1], dirsig[2], vector2[i][0], vector2[i][1], vector2[i][2]);
					double tmp_EED_B = this->angulo(dirant[0], dirant[1], dirant[2], vector2[i][0], vector2[i][1], vector2[i][2]);
//					angulos[i] =min( tmp_EED_A , tmp_EED_B );
					angulos[i] = tmp_EED_A ;
					if (tmp_EED_B<angulos[i] ) 
					{
						angulos[i]=tmp_EED_B;
					}
					if(maxmasa<this->vector[i][0]){
						maxmasa=this->vector[i][0];
					}
					
				}

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

							this->vector[i][0]=this->vector[j][0];
							this->vector[i][1]=this->vector[j][1];
							this->vector[i][2]=this->vector[j][2];
							this->vector[i][3]=this->vector[j][3];
							vector2[i][0]=vector2[j][0];
							vector2[i][1]=vector2[j][1];
							vector2[i][2]=vector2[j][2];
							angulos[i]=angulos[j];

							this->vector[j][0]=vectortemp[0];
							this->vector[j][1]=vectortemp[1];
							this->vector[j][2]=vectortemp[2];
							this->vector[j][3]=vectortemp[3];
							vector2[j][0]=vector2temp[0];
							vector2[j][1]=vector2temp[1];
							vector2[j][2]=vector2temp[2];
							angulos[j]=angulostemp;
						}
					}
					
				}


				for(i=this->label-2;i<this->label;i++){
					this->vector[i][0]=0;
				}



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

							this->vector[i][0]=this->vector[j][0];
							this->vector[i][1]=this->vector[j][1];
							this->vector[i][2]=this->vector[j][2];
							this->vector[i][3]=this->vector[j][3];
							vector2[i][0]=vector2[j][0];
							vector2[i][1]=vector2[j][1];
							vector2[i][2]=vector2[j][2];
							angulos[i]=angulos[j];

							this->vector[j][0]=vectortemp[0];
							this->vector[j][1]=vectortemp[1];
							this->vector[j][2]=vectortemp[2];
							this->vector[j][3]=vectortemp[3];
							vector2[j][0]=vector2temp[0];
							vector2[j][1]=vector2temp[1];
							vector2[j][2]=vector2temp[2];
							angulos[j]=angulostemp;
						}
					}
					
				}

				
				for(i=0;i<this->label;i++){
					if(maxmasa*this->minpropmasa<this->vector[i][0]){
						puntoactualdis[0]=(puntoactualarr[0]+vector2[i][0]-extprin[0])*espprin[0];
						puntoactualdis[1]=(puntoactualarr[1]+vector2[i][1]-extprin[2])*espprin[1];
						puntoactualdis[2]=(puntoactualarr[2]+vector2[i][2])*espprin[2];

						points2->InsertPoint(buenos2,puntoactual);
						points3->InsertPoint(buenos2,puntoactualdis);
						points4->InsertPoint(buenos2,puntoactual[0]+(radioactual/radiotrabajo)*(puntoactualdis[0]-puntoactual[0]), puntoactual[1]+(radioactual/radiotrabajo)*(puntoactualdis[1]-puntoactual[1]), puntoactual[2]+(radioactual/radiotrabajo)*(puntoactualdis[2]-puntoactual[2]));
						buenos2++;
					}
					
				}
					
			}
	
	}


}