lib/maracasVisuLib/src/kernel/include/axisExtractor.cxx

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