axisExtractor.cxx

Go to the documentation of this file.

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


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


#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 );
                                }

                        }
                }
                
        }
        

}