lib/maracasVisuLib/src/kernel/axisExtractor.cxx

   1 /*=========================================================================
   2
   3
   4 =========================================================================*/
   5
   6
   7 #include "vtkPolyData.h"
   8 #include "vtkObjectFactory.h"
   9
  10 //#include "vtkimagethreshold.h"
  11 //#include "vtkImageCast.h"
  12 //#include "vtkImageSeedConnectivity.h"
  13 //#include "vtkImageData.h"
  14 //#include "vtkMarchingCubes.h"
  15 //#include "vtkDoubleArray.h"
  16 //#include "vtkPointData.h"
  17 //#include "vtkextractvoi.h"
  18
  19 #include "vtkMath.h"
  20
  21
  22 #include "axisExtractor.h"
  23
  24 #define min(a,b) (((a)<(b))?(a):(b))
  25
  26
  27 vtkStandardNewMacro(axisExtractor);
  28
  29
  30 axisExtractor::axisExtractor() {
  31
  32         this->NumberOfRequiredInputs = 1;
  33
  34         this->humbral=0.45 ;
  35         this->maxpropradio=100;
  36         this->maxpropmasa=10;
  37         this->minpropmasa=0;
  38
  39         this->resample= vtkImageResample::New();
  40         this->resample->SetDimensionality (3);
  41
  42         this->extrac = vtkExtractVOI::New();
  43         this->extrac->SetInput(this->resample->GetOutput());
  44         this->extrac->SetSampleRate(1, 1, 1);
  45
  46         this->thresh = vtkImageThreshold::New();
  47         this->thresh->SetInput(this->extrac->GetOutput());
  48         this->thresh->SetInValue(255);
  49         this->thresh->SetOutValue(0);
  50         this->thresh->ReleaseDataFlagOff();
  51
  52         this->cast = vtkImageCast::New();
  53         this->cast->SetInput(this->thresh->GetOutput());
  54         this->cast->SetOutputScalarTypeToUnsignedChar();
  55
  56         this->connect = vtkImageSeedConnectivity::New();
  57         this->connect->SetInput(this->cast->GetOutput());
  58         this->connect->SetInputConnectValue(255);
  59         this->connect->SetOutputConnectedValue(255);
  60         this->connect->SetOutputUnconnectedValue(0);
  61
  62         this->dataprov=vtkImageData::New();
  63         this->dataprov->SetScalarTypeToChar();
  64
  65         this->datatotal=vtkImageData::New();
  66         this->datatotal->SetScalarTypeToChar();
  67
  68
  69
  70
  71 }
  72
  73 //----------------------------------------------------------------------------
  74 // Specify the input data or filter.
  75 void axisExtractor::SetInput(vtkImageData *input)
  76 {
  77         double minspac;
  78         double espprin[3];
  79
  80
  81         this->vtkProcessObject::SetNthInput(0, input);
  82         this->resample->SetInput(input);
  83         input->GetSpacing(espprin);
  84
  85
  86         minspac=min(espprin[0],min(espprin[1],espprin[2]));
  87
  88
  89         this->resample->SetAxisOutputSpacing( 0, minspac);
  90         this->resample->SetAxisOutputSpacing( 1, minspac);
  91         this->resample->SetAxisOutputSpacing( 2, minspac);
  92         this->resample->Update();
  93
  94
  95 }
  96
  97
  98 //----------------------------------------------------------------------------
  99 // Specify the input data or filter.
 100
 101 void axisExtractor::Execute()
 102 {
 103         this->points = vtkPoints::New();
 104         this->lineas = vtkCellArray::New();
 105         this->buenos=0;
 106         this->iter=0;
 107         this->flagg=0;
 108
 109         this->datatotal->Delete();
 110         this->datatotal=vtkImageData::New();
 111         this->datatotal->SetScalarTypeToUnsignedChar();
 112
 113         this->datatotal->SetExtent(this->resample->GetOutput()->GetExtent());
 114         this->datatotal->SetSpacing(this->resample->GetOutput()->GetSpacing());
 115
 116         this->blanquear2();
 117
 118
 119         while( !m_Stack0.empty()){
 120                 this->avanzar();
 121                 this->iter++;
 122         }
 123
 124
 125         this->GetOutput()->SetPoints (this->points);
 126         this->GetOutput()->SetLines(this->lineas);
 127         this->GetOutput()->Modified();
 128
 129         // this->GetOutput()->GetPointData()->SetVectors(this->salidas);
 130 }
 131
 132
 133
 134
 135
 136
 137
 138
 139 //----------------------------------------------------------------------------
 140 // Specify the input data or filter.
 141 vtkImageData *axisExtractor::GetInput()
 142 {
 143         if (this->NumberOfInputs < 1){
 144                 return NULL;
 145     }
 146
 147         return (vtkImageData *)(this->Inputs[0]);
 148 }
 149
 150
 151
 152
 153
 154
 155
 156 //----------------------------------------------------------------------------
 157 // Specify the input data or filter.
 158 vtkImageData *axisExtractor::GetVolumen()
 159 {
 160
 161         return this->datatotal;
 162 }
 163
 164
 165
 166
 167
 168
 169 //----------------------------------------------------------------------------
 170 void axisExtractor::PrintSelf(ostream& os, vtkIndent indent)
 171 {
 172         this->Superclass::PrintSelf(os,indent);
 173 }
 174
 175
 176
 177
 178
 179
 180
 181
 182 //----------------------------------------------------------------------------
 183 void axisExtractor::SetMaxPropRadio(double value)
 184 {
 185         this->maxpropradio=value;
 186 }
 187
 188
 189
 190 //----------------------------------------------------------------------------
 191 double axisExtractor::GetMaxPropRadio()
 192 {
 193         return this->maxpropradio;
 194 }
 195
 196
 197
 198 //----------------------------------------------------------------------------
 199
 200 void axisExtractor::SetHumbral(double value)
 201 {
 202         this->humbral=value;
 203 }
 204
 205
 206
 207   //----------------------------------------------------------------------------
 208
 209 double axisExtractor::GetHumbral()
 210 {
 211         return this->humbral;
 212 }
 213
 214
 215 //----------------------------------------------------------------------------
 216
 217 void axisExtractor::SetMaxPropMasa(double value)
 218 {
 219         this->maxpropmasa=value;
 220 }
 221
 222
 223
 224   //----------------------------------------------------------------------------
 225
 226 double axisExtractor::GetMaxPropMasa()
 227 {
 228         return this->maxpropmasa;
 229 }
 230
 231
 232 //----------------------------------------------------------------------------
 233
 234 void axisExtractor::SetMinPropMasa(double value)
 235 {
 236         this->minpropmasa=value;
 237 }
 238
 239
 240
 241   //----------------------------------------------------------------------------
 242
 243 double axisExtractor::GetMinPropMasa()
 244 {
 245         return this->minpropmasa;
 246 }
 247
 248
 249
 250
 251
 252
 253
 254         //----------------------------------------------------------------------------
 255 void axisExtractor::SetPoint(double puntoactualprov[3] )
 256 {
 257
 258         realtoreal(puntoactualprov,puntoactualprov);
 259
 260
 261         m_Stack0.push(puntoactualprov[0]);
 262         m_Stack1.push(puntoactualprov[1]);
 263         m_Stack2.push(puntoactualprov[2]);
 264         m_Stack3.push(puntoactualprov[0]);
 265         m_Stack4.push(puntoactualprov[1]);
 266         m_Stack5.push(puntoactualprov[2]);
 267         m_Stack.push(0);
 268
 269 }
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279 //----------------------------------------------------------------------------
 280 void axisExtractor::realtoreal(double a[3], double b[3] )
 281 {
 282
 283         double espprin[3];
 284         int extprin[6];
 285
 286         this->resample->GetOutput()->GetSpacing(espprin);
 287         this->resample->GetOutput()->GetExtent(extprin);
 288
 289         b[0]=a[0]+(espprin[0]*extprin[0]);
 290         b[1]=a[1]+(espprin[1]*extprin[2]);
 291         b[2]=a[2];
 292
 293
 294
 295 }
 296
 297
 298
 299
 300 //----------------------------------------------------------------------------
 301 void axisExtractor::realtoreal2(double a[3], double b[3] )
 302 {
 303
 304         double espprin[3];
 305         int extprin[6];
 306
 307         this->resample->GetOutput()->GetSpacing(espprin);
 308         this->resample->GetOutput()->GetExtent(extprin);
 309
 310         b[0]=a[0]-(espprin[0]*extprin[0]);
 311         b[1]=a[1]-(espprin[1]*extprin[2]);
 312         b[2]=a[2];
 313
 314
 315
 316 }
 317
 318
 319
 320
 321
 322 //----------------------------------------------------------------------------
 323 void axisExtractor::realtoindex(double a[3], int b[3] )
 324 {
 325         double c[3];
 326
 327         double minspac;
 328
 329         double espprin[3];
 330         int extprin[6];
 331
 332         this->resample->GetOutput()->GetSpacing(espprin);
 333         this->resample->GetOutput()->GetExtent(extprin);
 334
 335         minspac=min(espprin[0],min(espprin[1],espprin[2]));
 336
 337         c[0]=(a[0]/minspac);
 338         c[1]=(a[1]/minspac);
 339         c[2]=(a[2]/minspac);
 340
 341         b[0]=(int)(a[0]/minspac);
 342         b[1]=(int)(a[1]/minspac);
 343         b[2]=(int)(a[2]/minspac);
 344
 345         if(c[0]-b[0]>0.5){
 346                 b[0]+=1;
 347         }
 348         if(c[1]-b[1]>0.5){
 349                 b[1]+=1;
 350         }
 351         if(c[2]-b[2]>0.5){
 352                 b[2]+=1;
 353         }
 354
 355
 356 }
 357
 358
 359
 360
 361
 362
 363 //----------------------------------------------------------------------------
 364 void axisExtractor::indextoreal(int a[3], double b[3] )
 365 {
 366         double minspac;
 367
 368         double espprin[3];
 369         int extprin[6];
 370
 371         this->resample->GetOutput()->GetSpacing(espprin);
 372         this->resample->GetOutput()->GetExtent(extprin);
 373
 374         minspac=min(espprin[0],min(espprin[1],espprin[2]));
 375
 376         b[0]=(a[0])*minspac;
 377         b[1]=(a[1])*minspac;
 378         b[2]=a[2]*minspac;
 379 }
 380
 381
 382
 383
 384
 385 //----------------------------------------------------------------------------
 386 void axisExtractor::indextoreal(double a[3], double b[3] )
 387 {
 388         double minspac;
 389
 390         double espprin[3];
 391         int extprin[6];
 392
 393         this->resample->GetOutput()->GetSpacing(espprin);
 394         this->resample->GetOutput()->GetExtent(extprin);
 395
 396         minspac=min(espprin[0],min(espprin[1],espprin[2]));
 397
 398         b[0]=(a[0])*minspac;
 399         b[1]=(a[1])*minspac;
 400         b[2]=a[2]*minspac;
 401 }
 402
 403
 404
 405
 406 //----------------------------------------------------------------------------
 407 void axisExtractor::searc(int i, int j, int k)
 408 {
 409
 410         unsigned char *ptr;
 411         unsigned char   ptri;
 412         char *ptr2;
 413
 414         int radio;
 415
 416         int i2, j2, k2;
 417         int i3, j3, k3;
 418
 419         ptr=(unsigned char *)   this->connect->GetOutput()->GetScalarPointer(i,j,k);
 420         ptr2=(char *)   this->dataprov->GetScalarPointer(i,j,k);
 421
 422         ptri=*ptr;
 423
 424         int ext[6];
 425
 426         this->connect->GetOutput()->GetExtent(ext);
 427
 428         radio=(ext[1]-ext[0])/2;
 429
 430         i2=i-ext[0]-radio;
 431         j2=j-ext[2]-radio;
 432         k2=k-ext[4]-radio;
 433
 434
 435         if(ptri==255){
 436                 *ptr2=this->label;
 437                 this->vector[this->label-1][0]+=1;
 438                 this->vector[this->label-1][1]+=i;
 439                 this->vector[this->label-1][2]+=j;
 440                 this->vector[this->label-1][3]+=k;
 441         }
 442         else if(ptri==0){
 443                 *ptr2=-this->label2;
 444                 this->vectorb[this->label2-1][0]+=1;
 445                 this->vectorb[this->label2-1][1]+=i;
 446                 this->vectorb[this->label2-1][2]+=j;
 447                 this->vectorb[this->label2-1][3]+=k;
 448         }
 449
 450
 451         for(i3=-1;i3<=1;i3++){
 452                         for(j3=-1;j3<=1;j3++){
 453                                 for(k3=-1;k3<=1;k3++){
 454                                         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] ){
 455                                                 ptr=(unsigned char *)   this->connect->GetOutput()->GetScalarPointer(i+i3,j+j3,k+k3);
 456                                                 ptr2=(char *)   this->dataprov->GetScalarPointer(i+i3,j+j3,k+k3);
 457                                                 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))){
 458                                                         this->searc(i+i3, j+j3, k+k3);
 459                                                 }
 460                                                 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))){
 461                                                         this->searc(i+i3, j+j3, k+k3);
 462                                                 }
 463                                         }
 464                                 }
 465                         }
 466                 }
 467
 468
 469
 470
 471
 472 }
 473
 474
 475 //----------------------------------------------------------------------------
 476 void axisExtractor::find_components( )
 477 {
 478         int ext[6];
 479         int i, j, k, radio;
 480         int i2, j2, k2;
 481
 482
 483         unsigned char *ptr;
 484         char *ptr2;
 485
 486         this->connect->GetOutput()->GetExtent(ext);
 487
 488
 489         radio=(ext[1]-ext[0])/2;
 490
 491         this->label=0;
 492         this->label2=0;
 493
 494
 495
 496         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
 497                         for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
 498                                 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
 499                                         ptr=(unsigned char *)   this->connect->GetOutput()->GetScalarPointer(i,j,k);
 500                                         ptr2=(char *)   this->dataprov->GetScalarPointer(i,j,k);
 501
 502                                         if(*ptr==255 && *ptr2==0 && (radio*radio)>=(i2*i2)+(j2*j2)+(k2*k2) && ((radio-1)*(radio-1))<=(i2*i2)+(j2*j2)+(k2*k2)  ){
 503                                                 this->label=this->label+1;
 504                                                 this->vector[this->label-1][0]=0;
 505                                                 this->vector[this->label-1][1]=0;
 506                                                 this->vector[this->label-1][2]=0;
 507                                                 this->vector[this->label-1][3]=0;
 508                                                 this->searc(i, j, k);
 509
 510                                         }
 511                                         else if(*ptr==0 && *ptr2==0 && (radio*radio)>=(i2*i2)+(j2*j2)+(k2*k2) && ((radio-1)*(radio-1))<=(i2*i2)+(j2*j2)+(k2*k2)  ){
 512                                                 this->label2=this->label2+1;
 513                                                 this->vectorb[this->label2-1][0]=0;
 514                                                 this->vectorb[this->label2-1][1]=0;
 515                                                 this->vectorb[this->label2-1][2]=0;
 516                                                 this->vectorb[this->label2-1][3]=0;
 517                                                 this->searc(i, j, k);
 518
 519                                         }
 520                                 }
 521                         }
 522                 }
 523
 524
 525
 526
 527
 528
 529
 530 }
 531
 532
 533
 534
 535 //----------------------------------------------------------------------------
 536 unsigned short axisExtractor::maximo( )
 537 {
 538         int ext[6];
 539         int i, j, k;
 540         unsigned short max=0;
 541
 542         this->extrac->GetOutput()->GetExtent(ext);
 543
 544         unsigned short *ptr;
 545
 546
 547         for(i=ext[0];i<=ext[1];i++){
 548                         for(j=ext[2];j<=ext[3];j++){
 549                                 for(k=ext[4];k<=ext[5];k++){
 550                                         ptr=(unsigned short *)this->extrac->GetOutput()->GetScalarPointer(i,j,k);
 551                                         if(*ptr>max){
 552                                                 max=*ptr;
 553
 554                                         }
 555                                 }
 556                         }
 557                 }
 558
 559                 return max;
 560
 561
 562 }
 563
 564
 565
 566
 567 //----------------------------------------------------------------------------
 568 void axisExtractor::blanquear()
 569 {
 570         int ext[6];
 571         int i, j, k;
 572
 573         this->dataprov->GetExtent(ext);
 574
 575         char *ptr;
 576
 577
 578
 579         for(i=ext[0];i<=ext[1];i++){
 580                         for(j=ext[2];j<=ext[3];j++){
 581                                 for(k=ext[4];k<=ext[5];k++){
 582                                         ptr=(char *)    this->dataprov->GetScalarPointer(i,j,k);
 583                                         *ptr=0;
 584                                 }
 585                         }
 586                 }
 587 }
 588
 589
 590
 591
 592
 593
 594 //----------------------------------------------------------------------------
 595 void axisExtractor::blanquear2()
 596 {
 597         int ext[6];
 598         int i, j, k;
 599
 600         this->datatotal->GetExtent(ext);
 601
 602         unsigned char *ptr;
 603
 604
 605
 606         for(i=ext[0];i<=ext[1];i++){
 607                         for(j=ext[2];j<=ext[3];j++){
 608                                 for(k=ext[4];k<=ext[5];k++){
 609                                         ptr=(unsigned char *)   this->datatotal->GetScalarPointer(i,j,k);
 610                                         *ptr=0;
 611                                 }
 612                         }
 613                 }
 614 }
 615
 616
 617
 618
 619 //----------------------------------------------------------------------------
 620 void axisExtractor::copiar(vtkImageData *data, vtkImageData *data2 )
 621 {
 622         int ext[6];
 623         int i, j, k;
 624
 625         data->GetExtent(ext);
 626
 627         unsigned char *ptr, *ptr2;
 628
 629         for(i=ext[0];i<=ext[1];i++){
 630                         for(j=ext[2];j<=ext[3];j++){
 631                                 for(k=ext[4];k<=ext[5];k++){
 632                                         ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
 633                                         ptr2=(unsigned char *)  data2->GetScalarPointer(i,j,k);
 634                                         if(*ptr!=0 && *ptr2==0){
 635                                                 *ptr2=*ptr;
 636                                         }
 637                                 }
 638                         }
 639                 }
 640
 641         data2->Modified();
 642 }
 643
 644
 645
 646
 647
 648 //----------------------------------------------------------------------------
 649 double axisExtractor::distancia(double a[3], double b[3] )
 650 {
 651
 652
 653         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])));
 654 }
 655
 656
 657
 658
 659
 660
 661
 662
 663 //----------------------------------------------------------------------------
 664 int axisExtractor::envolumen(int a[3], vtkImageData *datae )
 665 {
 666
 667         int res;
 668
 669
 670         unsigned char *ptr;
 671
 672
 673         ptr=(unsigned char *)   datae->GetScalarPointer(a);
 674
 675         if(*ptr==0){
 676                 res=0;
 677         }
 678         else{
 679                 res=1;
 680         }
 681
 682
 683         return res;
 684 }
 685
 686
 687
 688
 689
 690
 691
 692 //----------------------------------------------------------------------------
 693 void axisExtractor::corte(double punto1[3], double punto2[3], double punto3[3], double centro[3],  double radio )
 694 {
 695
 696         double m1, m2, m3;
 697         double b1, b2, b3;
 698         double c0, c1, c2;
 699         double* roots;
 700         double root;
 701
 702
 703         m1=punto2[0]-punto1[0];
 704         m2=punto2[1]-punto1[1];
 705         m3=punto2[2]-punto1[2];
 706
 707         b1=punto1[0]-centro[0];
 708         b2=punto1[1]-centro[1];
 709         b3=punto1[2]-centro[2];
 710
 711         c0=m1*m1+m2*m2+m3*m3;
 712         c1=2*m1*b1+2*m2*b2+2*m3*b3;
 713         c2=b1*b1+b2*b2+b3*b3-radio*radio;
 714
 715         roots=vtkMath::SolveQuadratic  (c0, c1, c2);
 716
 717         if(roots[1]>=0 && roots[1]<=1){
 718                 root=roots[1];
 719         }
 720         else{
 721                 root=roots[2];
 722         }
 723
 724         punto3[0]=punto1[0]+root*(punto2[0]-punto1[0]);
 725         punto3[1]=punto1[1]+root*(punto2[1]-punto1[1]);
 726         punto3[2]=punto1[2]+root*(punto2[2]-punto1[2]);
 727
 728 }
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741 //----------------------------------------------------------------------------
 742 void axisExtractor::redondear(vtkImageData *data )
 743 {
 744         int ext[6];
 745         int i, j, k, radio;
 746         int i2, j2, k2;
 747
 748         data->GetExtent(ext);
 749
 750         unsigned char *ptr;
 751
 752         radio=(ext[1]-ext[0])/2;
 753
 754         double tmpsqrt;
 755         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
 756                         for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
 757                                 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
 758                                         ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
 759                                         tmpsqrt =       (i2*i2)+(j2*j2)+(k2*k2);
 760                                         if( radio<sqrt(tmpsqrt) ){
 761                                                 *ptr=0;
 762
 763                                         }
 764                                 }
 765                         }
 766                 }
 767 }
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781 //----------------------------------------------------------------------------
 782 void axisExtractor::avanzar()
 783 {
 784
 785
 786
 787
 788         double puntoactual[3];
 789         double puntoanterior[3];
 790         double puntoactualdis[3];
 791         double puntoactualcorre[3];
 792         int puntoactualarr[3];
 793         double radioactual;
 794         double dirant[3];
 795         double vector2[50][3];
 796         double vector2b[3];
 797         unsigned long vectortemp[4];
 798         double vector2temp[3];
 799         double maxmasa;
 800         int extint[6];
 801         int extintreal[6];
 802         int indesxdis[3];
 803         int indexp;
 804         double norvec;
 805         double provvc[3];
 806         double proprov;
 807         double puntoactualcorre2[3];
 808         int radiotrabajo=1;
 809         int radiocorte=0;
 810         int radiocorte2=0;
 811         int radiocorte3=0;
 812         int radiocorte4=0;
 813         int radiocorte5=0;
 814         int radiocorte6=0;
 815         int radiocorte7=0;
 816         int radiocorte8=0;
 817         int radiocorte9=0;
 818         int radiocorte10=0;
 819         int radiocorte11=0;
 820         int cantidadanterior=0;
 821         int cantidadanteriorb=0;
 822         int flag =0;
 823         int flag2 =0;
 824         int flag3 =0;
 825         int flag4 =0;
 826         int flag5 =0;
 827         int flag6 =0;
 828         int flag7 =0;
 829         int i, j;
 830
 831
 832
 833         if(!m_Stack0.empty()){
 834
 835                 indexp=m_Stack.top();
 836
 837                 puntoactual[0]=m_Stack0.top();
 838                 puntoactual[1]=m_Stack1.top();
 839                 puntoactual[2]=m_Stack2.top();
 840
 841                 puntoanterior[0]=m_Stack3.top();
 842                 puntoanterior[1]=m_Stack4.top();
 843                 puntoanterior[2]=m_Stack5.top();
 844
 845                 m_Stack0.pop();
 846                 m_Stack1.pop();
 847                 m_Stack2.pop();
 848                 m_Stack3.pop();
 849                 m_Stack4.pop();
 850                 m_Stack5.pop();
 851                 m_Stack.pop();
 852
 853                 dirant[0]=puntoanterior[0]-puntoactual[0];
 854                 dirant[1]=puntoanterior[1]-puntoactual[1];
 855                 dirant[2]=puntoanterior[2]-puntoactual[2];
 856
 857                 radioactual=distancia(puntoactual, puntoanterior);
 858
 859                 realtoindex(puntoactual, puntoactualarr);
 860
 861                 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++){
 862
 863                         extint[0]=puntoactualarr[0]-radiotrabajo;
 864                         extint[1]=puntoactualarr[0]+radiotrabajo;
 865                         extint[2]=puntoactualarr[1]-radiotrabajo;
 866                         extint[3]=puntoactualarr[1]+radiotrabajo;
 867                         extint[4]=puntoactualarr[2]-radiotrabajo;
 868                         extint[5]=puntoactualarr[2]+radiotrabajo;
 869
 870                         extrac->SetVOI(extint);
 871                         extrac->UpdateWholeExtent();
 872                         extrac->Update();
 873                         extrac->GetOutput()->GetExtent(extintreal);
 874
 875                         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]){
 876                                 flag4 =1;
 877                                 radiotrabajo++;
 878                                 break;
 879                         }
 880
 881                         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]  ){
 882
 883                                 thresh->ThresholdByUpper(this->maximo()*this->humbral);
 884                                 thresh->UpdateWholeExtent();
 885                                 thresh->Update();
 886
 887                                 redondear(thresh->GetOutput());
 888
 889                                 cast->UpdateWholeExtent();
 890                                 cast->Update();
 891
 892                                 connect->RemoveAllSeeds();
 893                                 connect->AddSeed(puntoactualarr[0],puntoactualarr[1],puntoactualarr[2]);
 894                                 connect->UpdateWholeExtent();
 895                                 connect->Update();
 896
 897                                 this->dataprov->Delete();
 898
 899                                 this->dataprov=vtkImageData::New();
 900                                 this->dataprov->SetScalarTypeToChar();
 901
 902                                 this->dataprov->SetExtent(this->connect->GetOutput()->GetExtent());
 903                                 this->dataprov->SetSpacing(this->connect->GetOutput()->GetSpacing());
 904
 905                                 this->blanquear();
 906
 907                                 this->find_components();
 908
 909                                 if(this->label2>0 && flag==0){
 910                                         radiocorte=radiotrabajo;
 911                                         vector2b[0]=((double)vectorb[0][1]/(double)vectorb[0][0])-(double)puntoactualarr[0];
 912                                         vector2b[1]=((double)vectorb[0][2]/(double)vectorb[0][0])-(double)puntoactualarr[1];
 913                                         vector2b[2]=((double)vectorb[0][3]/(double)vectorb[0][0])-(double)puntoactualarr[2];
 914                                         flag =1;
 915                                 }
 916                                 if(this->label2>1 && flag5==0){
 917                                         radiocorte5=radiotrabajo;
 918                                         flag5 =1;
 919                                 }
 920
 921                                 if( this->label>1 && flag2==0){
 922                                         radiocorte2=radiotrabajo;
 923                                         flag2 =1;
 924                                 }
 925
 926                                 if(radiocorte5==0){
 927                                         radiocorte6=radiocorte2;
 928                                 }
 929                                 else{
 930                                         radiocorte6=min(radiocorte2, radiocorte5);
 931                                 }
 932
 933                                 if((radiocorte6-radiocorte)>1 && flagg<=4){
 934                                         flag4 =1;
 935                                         radiotrabajo++;
 936                                         flag7 =0;
 937                                         break;
 938                                 }
 939                                 else{
 940                                         flag7 =1;
 941                                 }
 942
 943                                 if(this->label>2 && flag3==0){
 944                                         radiocorte3=radiotrabajo;
 945                                         flag3 =1;
 946                                 }
 947                                 if( this->label==2 && cantidadanterior!=this->label){
 948                                         radiocorte7=radiotrabajo;
 949
 950                                 }
 951                                 if( this->label==3 && cantidadanterior!=this->label){
 952                                         radiocorte8=radiotrabajo;
 953                                 }
 954                                 if( cantidadanterior!=this->label || cantidadanteriorb!=this->label2){
 955                                         if(radiocorte2==0){
 956                                                 radiocorte4=radiotrabajo+10*(radiotrabajo-radiocorte9);
 957                                         }
 958                                         else{
 959                                                 radiocorte4=radiotrabajo+(radiotrabajo-radiocorte9);
 960                                         }
 961                                         radiocorte9=radiotrabajo;
 962                                 }
 963
 964                                 cantidadanterior=this->label;
 965                                 cantidadanteriorb=this->label2;
 966
 967                         }
 968
 969                 }
 970                 radiotrabajo--;
 971                 radiocorte10=radiotrabajo;
 972
 973                 if(flag7==1){
 974                         if(this->label>=2 && radiocorte10!=radiocorte9+1){
 975                                 radiotrabajo=radiocorte9+1;
 976
 977                                 extint[0]=puntoactualarr[0]-radiotrabajo;
 978                                 extint[1]=puntoactualarr[0]+radiotrabajo;
 979                                 extint[2]=puntoactualarr[1]-radiotrabajo;
 980                                 extint[3]=puntoactualarr[1]+radiotrabajo;
 981                                 extint[4]=puntoactualarr[2]-radiotrabajo;
 982                                 extint[5]=puntoactualarr[2]+radiotrabajo;
 983
 984                                 extrac->SetVOI(extint);
 985                                 extrac->UpdateWholeExtent();
 986                                 extrac->Update();
 987                                 extrac->GetOutput()->GetExtent(extintreal);
 988
 989                                 thresh->ThresholdByUpper(this->maximo()*this->humbral);
 990                                 thresh->UpdateWholeExtent();
 991                                 thresh->Update();
 992
 993                                 redondear(thresh->GetOutput());
 994
 995                                 cast->UpdateWholeExtent();
 996                                 cast->Update();
 997
 998                                 connect->RemoveAllSeeds();
 999                                 connect->AddSeed(puntoactualarr[0],puntoactualarr[1],puntoactualarr[2]);
1000                                 connect->UpdateWholeExtent();
1001                                 connect->Update();
1002
1003                                 this->dataprov->Delete();
1004
1005                                 this->dataprov=vtkImageData::New();
1006                                 this->dataprov->SetScalarTypeToChar();
1007
1008                                 this->dataprov->SetExtent(this->connect->GetOutput()->GetExtent());
1009                                 this->dataprov->SetSpacing(this->connect->GetOutput()->GetSpacing());
1010
1011
1012                                 this->blanquear();
1013
1014                                 this->find_components();
1015                         }
1016                 }
1017
1018                 if(flag7==1){
1019                         flagg=0;
1020                         if(this->label>1){
1021                                 //printf("no corrigio\n");
1022
1023                                 realtoreal2(puntoactual, provvc);
1024                                 points->InsertPoint(buenos,provvc);
1025
1026                                 if(buenos>0){
1027                                         lineas->InsertNextCell(2);
1028                                         lineas->InsertCellPoint(indexp);
1029                                         lineas->InsertCellPoint(buenos);
1030                                 }
1031
1032                                 buenos++;
1033                         }
1034                 }
1035                 else{
1036                         flagg++;
1037                         //printf("corrigio\n");
1038
1039                         provvc[0]=0;
1040                         provvc[1]=0;
1041                         provvc[2]=0;
1042
1043                         norvec=distancia(provvc, vector2b);
1044
1045                         proprov=(radiocorte6-radiocorte)/(2*norvec);
1046
1047                         puntoactualcorre[0]=(puntoactualarr[0]-(vector2b[0]*proprov));
1048                         puntoactualcorre[1]=(puntoactualarr[1]-(vector2b[1]*proprov));
1049                         puntoactualcorre[2]=(puntoactualarr[2]-(vector2b[2]*proprov));
1050
1051                         indextoreal(puntoactualcorre, puntoactualcorre2);
1052
1053                         m_Stack0.push(puntoactualcorre2[0]);
1054                         m_Stack1.push(puntoactualcorre2[1]);
1055                         m_Stack2.push(puntoactualcorre2[2]);
1056                         m_Stack3.push(puntoanterior[0]);
1057                         m_Stack4.push(puntoanterior[1]);
1058                         m_Stack5.push(puntoanterior[2]);
1059                         m_Stack.push(indexp);
1060
1061
1062                 }
1063
1064                 if(flag7==1){
1065
1066                         if(this->label>1){
1067
1068                                 maxmasa=0;
1069                                 for(i=0;i<this->label;i++){
1070                                         vector2[i][0]=((double)vector[i][1]/(double)vector[i][0])-(double)puntoactualarr[0];
1071                                         vector2[i][1]=((double)vector[i][2]/(double)vector[i][0])-(double)puntoactualarr[1];
1072                                         vector2[i][2]=((double)vector[i][3]/(double)vector[i][0])-(double)puntoactualarr[2];
1073                                         if(maxmasa<vector[i][0]){
1074                                                 maxmasa=vector[i][0];
1075                                         }
1076                                 }
1077
1078                                 for(i=0;i<this->label;i++){
1079                                         for(j=i;j<this->label;j++){
1080                                                 if(vector[j][0]<vector[i][0]){
1081                                                         vectortemp[0]=vector[i][0];
1082                                                         vectortemp[1]=vector[i][1];
1083                                                         vectortemp[2]=vector[i][2];
1084                                                         vectortemp[3]=vector[i][3];
1085                                                         vector2temp[0]=vector2[i][0];
1086                                                         vector2temp[1]=vector2[i][1];
1087                                                         vector2temp[2]=vector2[i][2];
1088
1089                                                         vector[i][0]=vector[j][0];
1090                                                         vector[i][1]=vector[j][1];
1091                                                         vector[i][2]=vector[j][2];
1092                                                         vector[i][3]=vector[j][3];
1093                                                         vector2[i][0]=vector2[j][0];
1094                                                         vector2[i][1]=vector2[j][1];
1095                                                         vector2[i][2]=vector2[j][2];
1096
1097                                                         vector[j][0]=vectortemp[0];
1098                                                         vector[j][1]=vectortemp[1];
1099                                                         vector[j][2]=vectortemp[2];
1100                                                         vector[j][3]=vectortemp[3];
1101                                                         vector2[j][0]=vector2temp[0];
1102                                                         vector2[j][1]=vector2temp[1];
1103                                                         vector2[j][2]=vector2temp[2];
1104                                                 }
1105                                         }
1106
1107                                 }
1108
1109
1110                                 for(i=0;i<this->label;i++){
1111                                         if(maxmasa*this->minpropmasa<vector[i][0]){
1112
1113                                                 indesxdis[0]=(int) (puntoactualarr[0]+vector2[i][0]);
1114                                                 indesxdis[1]=(int) (puntoactualarr[1]+vector2[i][1]);
1115                                                 indesxdis[2]=(int) (puntoactualarr[2]+vector2[i][2]);
1116
1117                                                 indextoreal(indesxdis, puntoactualdis);
1118
1119                                                 if(envolumen(indesxdis, datatotal)==0){
1120                                                         m_Stack0.push(puntoactualdis[0]);
1121                                                         m_Stack1.push(puntoactualdis[1]);
1122                                                         m_Stack2.push(puntoactualdis[2]);
1123                                                         m_Stack3.push(puntoactual[0]);
1124                                                         m_Stack4.push(puntoactual[1]);
1125                                                         m_Stack5.push(puntoactual[2]);
1126                                                         m_Stack.push(buenos-1);
1127                                                 }
1128
1129                                         }
1130
1131                                 }
1132
1133                                 if(this->label>1 ){
1134                                         copiar(this->connect->GetOutput(),  datatotal );
1135                                 }
1136
1137                         }
1138                 }
1139
1140         }
1141
1142
1143 }