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

   1 /*=========================================================================
   2
   3
   4 =========================================================================*/
   5 #include "axisExtractor02.h"
   6 #include "vtkPolyData.h"
   7 #include "vtkObjectFactory.h"
   8 #include "vtkImageThreshold.h"
   9 #include "vtkImageCast.h"
  10 #include "vtkImageSeedConnectivity.h"
  11 #include "vtkImageData.h"
  12 #include "vtkMarchingCubes.h"
  13 #include "vtkDoubleArray.h"
  14 #include "vtkPointData.h"
  15 #include "vtkExtractVOI.h"
  16 #include "vtkImageConstantPad.h"
  17 #include "vtkImageTranslateExtent.h"
  18 #include "vtkMath.h"
  19 #include "vtkTransformPolyDataFilter.h"
  20 #include "vtkTransform.h"
  21 #include <time.h>
  22
  23
  24
  25
  26
  27 vtkStandardNewMacro(axisExtractor02);
  28
  29
  30 //----------------------------------------------------------------------------
  31 void axisExtractor02::SetParam(double value)
  32 {
  33         this->param=value;
  34 }
  35
  36
  37
  38 //----------------------------------------------------------------------------
  39 double axisExtractor02::GetParam()
  40 {
  41         return this->param;
  42 }
  43
  44
  45
  46
  47 //----------------------------------------------------------------------------
  48 void axisExtractor02::SetParam2(double value)
  49 {
  50         this->param2=value;
  51 }
  52
  53
  54
  55 //----------------------------------------------------------------------------
  56 double axisExtractor02::GetParam2()
  57 {
  58         return this->param2;
  59 }
  60
  61
  62
  63
  64 //----------------------------------------------------------------------------
  65 void axisExtractor02::SetParam3(double value)
  66 {
  67         this->param3=value;
  68 }
  69
  70
  71
  72 //----------------------------------------------------------------------------
  73 double axisExtractor02::GetParam3()
  74 {
  75         return this->param3;
  76 }
  77
  78
  79
  80 //----------------------------------------------------------------------------
  81 void axisExtractor02::SetMaxant(int value)
  82 {
  83         this->maxant=value;
  84 }
  85
  86
  87
  88 //----------------------------------------------------------------------------
  89 int axisExtractor02::GetMaxant()
  90 {
  91         return this->maxant;
  92 }
  93
  94
  95
  96 //----------------------------------------------------------------------------
  97 void axisExtractor02::SetMinant(int value)
  98 {
  99         this->minant=value;
 100 }
 101
 102
 103
 104 //----------------------------------------------------------------------------
 105 int axisExtractor02::GetMinant()
 106 {
 107         return this->minant;
 108 }
 109
 110
 111 //----------------------------------------------------------------------------
 112 void axisExtractor02::SetPoint(double value[3])
 113 {
 114
 115
 116         double spa[3];
 117
 118         this->GetInput()->GetSpacing (spa);
 119
 120         value[0]=value[0]+maxant*spa[0];
 121         value[1]=value[1]+maxant*spa[1];
 122         value[2]=value[2]+maxant*spa[2];
 123
 124         this->m_Stack0.push_front(value[0]);
 125         this->m_Stack1.push_front(value[1]);
 126         this->m_Stack2.push_front(value[2]);
 127         this->m_Stack3.push_front(value[0]);
 128         this->m_Stack4.push_front(value[1]);
 129         this->m_Stack5.push_front(value[2]);
 130         this->m_Stack6.push_front(value[0]);
 131         this->m_Stack7.push_front(value[1]);
 132         this->m_Stack8.push_front(value[2]);
 133         this->m_Stack.push_front(0);
 134         this->m_Stackr.push_front(0);
 135         this->m_Stackra.push_front(0);
 136         this->m_Stackrp.push_front(0);
 137 }
 138
 139
 140
 141 //----------------------------------------------------------------------------
 142 vtkImageData *axisExtractor02::GetVolumen()
 143 {
 144
 145         vtkImageTranslateExtent *trans;
 146
 147         trans = vtkImageTranslateExtent::New();
 148
 149         trans->SetInput(this->data4);
 150
 151         trans->SetTranslation (maxant, maxant, maxant);
 152
 153         trans->Update();
 154
 155         return trans->GetOutput();
 156 }
 157
 158
 159
 160
 161
 162 //----------------------------------------------------------------------------
 163 vtkPolyData  *axisExtractor02::GetOutput ()
 164 {
 165
 166
 167
 168
 169         vtkTransform *transL1;
 170         vtkTransformPolyDataFilter *trans;
 171
 172         transL1 = vtkTransform::New();
 173         trans = vtkTransformPolyDataFilter::New();
 174
 175 // EED 30 Oct 2006
 176         double spa[3];
 177         this->GetInput()->GetSpacing (spa);
 178         transL1->Translate(-maxant*spa[0], -maxant*spa[1], -maxant*spa[2]);
 179 //      transL1->Translate(-maxant, -maxant, -maxant);
 180
 181
 182         trans->SetInput((vtkPolyData *)(this->Outputs[0]));
 183
 184         trans->SetTransform (transL1);
 185
 186         trans->Update();
 187
 188         return trans->GetOutput();
 189 }
 190
 191
 192
 193
 194
 195
 196
 197 //----------------------------------------------------------------------------
 198 axisExtractor02::axisExtractor02() {
 199
 200
 201         this->NumberOfRequiredInputs = 1;
 202         this->param=1;
 203         this->param2=1;
 204         this->param3=0.5;
 205         this->param4=0;
 206
 207         /*this->resample= vtkImageResample::New();
 208         this->resample->SetDimensionality (3);*/
 209
 210         this->data4=vtkImageData::New();
 211
 212         this->data1=vtkImageData::New();
 213         this->data2=vtkImageData::New();
 214         this->data3=vtkImageData::New();
 215         this->data6=vtkImageData::New();
 216
 217         this->extrac = vtkExtractVOI::New();
 218 //      this->extrac->SetInput(resample->GetOutput());
 219         this->extrac->SetSampleRate(1, 1, 1);
 220
 221         this->connect = vtkImageSeedConnectivity::New();
 222         this->connect->SetInput(this->data1);
 223         this->connect->SetInputConnectValue(255);
 224         this->connect->SetOutputConnectedValue(255);
 225         this->connect->SetOutputUnconnectedValue(0);
 226
 227         this->distance= vtkImageEuclideanDistance::New();
 228         this->distance->SetInput(this->connect->GetOutput());
 229         this->distance->InitializeOn();
 230         this->distance->ConsiderAnisotropyOff();
 231
 232
 233 }
 234
 235
 236
 237 //----------------------------------------------------------------------------
 238 void axisExtractor02::SetInput(vtkImageData *input)
 239 {
 240
 241         vtkImageConstantPad *pad;
 242         vtkImageTranslateExtent *trans;
 243
 244         pad = vtkImageConstantPad::New();
 245         trans = vtkImageTranslateExtent::New();
 246
 247         pad->SetInput(input);
 248         trans->SetInput(pad->GetOutput());
 249
 250
 251         pad->SetConstant(0);
 252
 253         pad->SetInput(input);
 254
 255         int ext[6];
 256
 257         input->GetExtent(ext);
 258
 259         ext[0]=ext[0]-maxant;
 260         ext[1]=ext[1]+maxant;
 261         ext[2]=ext[2]-maxant;
 262         ext[3]=ext[3]+maxant;
 263         ext[4]=ext[4]-maxant;
 264         ext[5]=ext[5]+maxant;
 265
 266         pad->SetOutputWholeExtent(ext);
 267
 268         trans->SetTranslation (maxant, maxant, maxant);
 269
 270         trans->Update();
 271
 272
 273         //this->vtkProcessObject::SetNthInput(0, input);
 274         this->vtkProcessObject::SetNthInput(0, trans->GetOutput());
 275
 276         //this->GetInput()->SetOrigin(0,0,0);
 277         //this->resample->SetInput(this->GetInput());
 278         this->extrac->SetInput(this->GetInput());
 279
 280 }
 281
 282
 283
 284 //----------------------------------------------------------------------------
 285 void axisExtractor02::distanciaejes(vtkPolyData *eje1,  vtkPolyData *eje2)
 286 {
 287         vtkIdType totalpuntos1=eje1->GetNumberOfPoints();
 288         vtkIdType totalpuntos2=eje2->GetNumberOfPoints();
 289         vtkIdType totallineas=eje2->GetNumberOfLines();
 290
 291         vtkCell * lineas;
 292
 293         int i, j;
 294
 295         double point[3], point2[3], point3[3], point4[3];
 296
 297         double v[3], v2[3];
 298
 299         double x, dist, min, y;
 300         FILE *ff;
 301
 302         ff=fopen("comparacion.txt","w");
 303
 304
 305
 306         for(i=0;i<totalpuntos1;i++){
 307                 eje1->GetPoint(i, point);
 308                 min =5000;
 309                 for(j=0;j<totallineas;j++){
 310                         lineas=eje2->GetCell(j);
 311
 312                                         lineas->GetPoints()->GetPoint(0, point3);
 313                                         lineas->GetPoints()->GetPoint(1, point2);
 314                                         v[0]=point3[0]-point2[0];
 315                                         v[1]=point3[1]-point2[1];
 316                                         v[2]=point3[2]-point2[2];
 317
 318                                         x=((v[0]*point[0])+(v[1]*point[1])+(v[2]*point[2])-(v[0]*point2[0])-(v[1]*point2[1])-(v[2]*point2[2]))/((v[0]*v[0])+(v[1]*v[1])+(v[2]*v[2]));
 319
 320
 321                                         point4[0]=point2[0]+(x*v[0]);
 322                                         point4[1]=point2[1]+(x*v[1]);
 323                                         point4[2]=point2[2]+(x*v[2]);
 324
 325                                         v2[0]=point4[0]-point[0];
 326                                         v2[1]=point4[1]-point[1];
 327                                         v2[2]=point4[2]-point[2];
 328
 329                                         y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);
 330
 331                                         if(x>=0 && x<=1){
 332                                                 dist=sqrt(((point4[0]-point[0])*(point4[0]-point[0]))+((point4[1]-point[1])*(point4[1]-point[1]))+((point4[2]-point[2])*(point4[2]-point[2])));
 333                                                 if(dist<min){
 334                                                         min=dist;
 335
 336
 337                                                         /*printf("punto: %f %f %f\n", point[0], point[1], point[2]);
 338                                                         printf("punto2: %f %f %f\n", point2[0], point2[1], point2[2]);
 339                                                         printf("punto3: %f %f %f\n", point3[0], point3[1], point3[2]);
 340                                                         printf("punto4: %f %f %f\n", point4[0], point4[1], point4[2]);
 341                                                         printf("v: %f %f %f\n", v[0], v[1], v[2]);
 342                                                         printf("v2: %f %f %f\n", v2[0], v2[1], v2[2]);
 343                                                         printf("%f\n", min);
 344                                                         printf("%f\n", x);
 345                                                         printf("%f\n", y);*/
 346
 347
 348
 349                                                 }
 350                                         }
 351
 352                         /*}*/
 353                 }
 354                 if(min<5){
 355                         fprintf(ff,"%f\n", min);
 356                 }
 357
 358
 359
 360         }
 361
 362         fclose(ff);
 363
 364
 365
 366 }
 367
 368
 369
 370
 371 //----------------------------------------------------------------------------
 372 void axisExtractor02::Execute()
 373 {
 374         flagg=0;
 375         flagg2=0;
 376         frama=0;
 377         fseg=0;
 378         mejordst=0;
 379         mejorrad=0;
 380         mejorcant=0;
 381
 382         //double minspac;
 383
 384         time_t start,end;
 385         double dif;
 386 //      FILE *ff;
 387
 388         this->GetInput()->GetExtent(extprin0);
 389         this->GetInput()->GetExtent(extprin);
 390         this->GetInput()->GetSpacing(espprin);
 391
 392 // EED
 393 //      stream  = fopen( "resultado.txt", "w" );
 394 //      fprintf(stream, "accion\tindexp\tbuenos\tradiotrabajo\tradioactual\tradiopred\tradioanterior\tmejordst\tmejorrad\tmejorcant\tflag2\tflag3\tflag4\tflag5\tflag6\tflag7\tflag8\tflag9\tflag10\tflag11\tflag12\tflag13\tflag14\tflag15\tflag16\tflag17\tflag18\tflag19\tflag20\tflagg\tflagg2\tcantidad\tcantidadb\tburned\tpuntoactual\tpuntoanterior\tpuntopre\tmejor\n");
 395
 396         /*minspac=min(espprin[0],min(espprin[1],espprin[2]));
 397
 398         this->resample->SetAxisOutputSpacing( 0, minspac);
 399         this->resample->SetAxisOutputSpacing( 1, minspac);
 400         this->resample->SetAxisOutputSpacing( 2, minspac);
 401         resample->Update();*/
 402
 403
 404         this->data4->SetScalarTypeToUnsignedChar();
 405         this->data4->SetExtent(this->GetInput()->GetExtent());
 406         this->data4->SetSpacing(this->GetInput()->GetSpacing());
 407
 408         this->blanquear(this->data4);
 409
 410         this->points = vtkPoints::New();
 411         this->buenos=0;
 412         this->lineas = vtkCellArray::New();
 413
 414
 415         time (&start);
 416         this->todo();
 417
 418 // EED
 419 //      fclose( stream );
 420
 421         ((vtkPolyData *)this->Outputs[0])->SetPoints (this->points);
 422         ((vtkPolyData *)this->Outputs[0])->SetLines(this->lineas);
 423         time (&end);
 424         dif = difftime (end,start);
 425
 426 //      ff=fopen("time.txt","w");
 427 //      fprintf(ff,"%d \n",this->buenos);
 428 //      fprintf(ff,"%.2lf \n",dif);
 429 //      fprintf(ff,"%.2lf \n",(double)this->buenos/dif);
 430 //      fclose(ff);
 431
 432 }
 433
 434
 435
 436
 437
 438
 439 //----------------------------------------------------------------------------
 440 vtkImageData *axisExtractor02::GetInput()
 441 {
 442         if (this->NumberOfInputs < 1){
 443                 return NULL;
 444     }
 445
 446         return (vtkImageData *)(this->Inputs[0]);
 447 }
 448
 449
 450
 451
 452
 453
 454
 455 //----------------------------------------------------------------------------
 456 void axisExtractor02::PrintSelf(ostream& os, vtkIndent indent)
 457 {
 458         this->Superclass::PrintSelf(os,indent);
 459 }
 460
 461
 462
 463
 464 //----------------------------------------------------------------------------
 465 void axisExtractor02::realtoreal(double a[3], double b[3] )
 466 {
 467
 468         b[0]=a[0]+(espprin[0]*extprin[0]);
 469         b[1]=a[1]+(espprin[1]*extprin[2]);
 470         b[2]=a[2];
 471
 472
 473
 474 }
 475
 476
 477
 478 //----------------------------------------------------------------------------
 479 void axisExtractor02::realtoreal2(double a[3], double b[3] )
 480 {
 481
 482         b[0]=a[0]-(espprin[0]*extprin[0]);
 483         b[1]=a[1]-(espprin[1]*extprin[2]);
 484         b[2]=a[2];
 485
 486 }
 487
 488
 489 //----------------------------------------------------------------------------
 490 void axisExtractor02::realtoindex(double a[3], int b[3] )
 491 {
 492         double c[3];
 493
 494         double minspac;
 495
 496 // EED 15 Mai 2007 .NET
 497 //      minspac=min(espprin[0],min(espprin[1],espprin[2]));
 498
 499         minspac=espprin[0];
 500         if (espprin[1]<minspac) { minspac=espprin[1]; }
 501         if (espprin[2]<minspac) { minspac=espprin[2]; }
 502
 503
 504         c[0]=(a[0]/minspac);
 505         c[1]=(a[1]/minspac);
 506         c[2]=(a[2]/minspac);
 507
 508         b[0]=(int)(a[0]/minspac);
 509         b[1]=(int)(a[1]/minspac);
 510         b[2]=(int)(a[2]/minspac);
 511
 512         if(c[0]-b[0]>0.5){
 513                 b[0]+=1;
 514         }
 515         if(c[1]-b[1]>0.5){
 516                 b[1]+=1;
 517         }
 518         if(c[2]-b[2]>0.5){
 519                 b[2]+=1;
 520         }
 521
 522
 523 }
 524
 525
 526
 527
 528
 529
 530
 531 //----------------------------------------------------------------------------
 532 void axisExtractor02::indextoreal(int a[3], double b[3] )
 533 {
 534         double minspac;
 535
 536 // EED 15 Mai 2007 .NET
 537 //      minspac=min(espprin[0],min(espprin[1],espprin[2]));
 538         minspac=espprin[0];
 539         if (espprin[1]<minspac) { minspac=espprin[1]; }
 540         if (espprin[2]<minspac) { minspac=espprin[2]; }
 541
 542         b[0]=(a[0])*minspac;
 543         b[1]=(a[1])*minspac;
 544         b[2]=a[2]*minspac;
 545 }
 546
 547
 548
 549 //----------------------------------------------------------------------------
 550 void axisExtractor02::indextoreal(double a[3], double b[3] )
 551 {
 552         double minspac;
 553
 554 // EED 15 Mai 2007 .NET
 555 //      minspac=min(espprin[0],min(espprin[1],espprin[2]));
 556
 557         minspac=espprin[0];
 558         if (espprin[1]<minspac) { minspac=espprin[1]; }
 559         if (espprin[2]<minspac) { minspac=espprin[2]; }
 560
 561         b[0]=(a[0])*minspac;
 562         b[1]=(a[1])*minspac;
 563         b[2]=a[2]*minspac;
 564 }
 565
 566
 567
 568 //----------------------------------------------------------------------------
 569 double axisExtractor02::distanciaejepunto(double point[3], double point2[3], double point3[3])
 570 {
 571
 572
 573         double point4[3];
 574
 575         double v[3], v2[3], v3[3];
 576
 577         double x, dist=0,  y;
 578
 579
 580         v[0]=point3[0]-point2[0];
 581         v[1]=point3[1]-point2[1];
 582         v[2]=point3[2]-point2[2];
 583
 584         v3[0]=point[0]-point2[0];
 585         v3[1]=point[1]-point2[1];
 586         v3[2]=point[2]-point2[2];
 587
 588         x=((v[0]*v3[0])+(v[1]*v3[1])+(v[2]*v3[2]))/((v[0]*v[0])+(v[1]*v[1])+(v[2]*v[2]));
 589
 590
 591         point4[0]=point2[0]+(x*v[0]);
 592         point4[1]=point2[1]+(x*v[1]);
 593         point4[2]=point2[2]+(x*v[2]);
 594
 595         v2[0]=point4[0]-point[0];
 596         v2[1]=point4[1]-point[1];
 597         v2[2]=point4[2]-point[2];
 598
 599         y=(v[0]*v2[0])+(v[1]*v2[1])+(v[2]*v2[2]);
 600
 601         dist=sqrt((v2[0]*v2[0])+(v2[1]*v2[1])+(v2[2]*v2[2]));
 602
 603         return dist;
 604
 605 }
 606
 607
 608
 609 //----------------------------------------------------------------------------
 610 double axisExtractor02::proporcioejepunto(double point[3], double point2[3], double point3[3])
 611 {
 612
 613
 614         double v[3], v3[3];
 615
 616         double x;
 617
 618         v[0]=point3[0]-point2[0];
 619         v[1]=point3[1]-point2[1];
 620         v[2]=point3[2]-point2[2];
 621
 622         v3[0]=point[0]-point2[0];
 623         v3[1]=point[1]-point2[1];
 624         v3[2]=point[2]-point2[2];
 625
 626         x=((v[0]*v3[0])+(v[1]*v3[1])+(v[2]*v3[2]))/((v[0]*v[0])+(v[1]*v[1])+(v[2]*v[2]));
 627
 628         return x;
 629
 630
 631
 632
 633 }
 634
 635
 636
 637
 638
 639
 640
 641
 642 //----------------------------------------------------------------------------
 643 void axisExtractor02::searc(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
 644 {
 645
 646         unsigned char *ptr;
 647         unsigned char *ptr2;
 648
 649         int radio;
 650
 651         int i2, j2, k2;
 652         int i3, j3, k3;
 653
 654         ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
 655         ptr2=(unsigned char *)  data2->GetScalarPointer(i,j,k);
 656
 657         int ext[6];
 658
 659         data->GetExtent(ext);
 660
 661         radio=(ext[1]-ext[0])/2;
 662
 663         i2=i-ext[0]-radio;
 664         j2=j-ext[2]-radio;
 665         k2=k-ext[4]-radio;
 666
 667         *ptr2=label;
 668         vector[label-1][0]+=1;
 669         vector[label-1][1]+=i;
 670         vector[label-1][2]+=j;
 671         vector[label-1][3]+=k;
 672
 673         for(i3=-1;i3<=1;i3++){
 674                         for(j3=-1;j3<=1;j3++){
 675                                 for(k3=-1;k3<=1;k3++){
 676                                         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] ){
 677                                                 ptr=(unsigned char *)   data->GetScalarPointer(i+i3,j+j3,k+k3);
 678                                                 ptr2=(unsigned char *)  data2->GetScalarPointer(i+i3,j+j3,k+k3);
 679                                                 if(*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)))){
 680                                                         searc(i+i3, j+j3, k+k3, data, data2, label, vector);
 681                                                 }
 682                                         }
 683                                 }
 684                         }
 685                 }
 686
 687 }
 688
 689
 690
 691
 692 //----------------------------------------------------------------------------
 693 void axisExtractor02::searcb(int i, int j, int k, vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
 694 {
 695
 696         unsigned char *ptr;
 697         unsigned char *ptr2;
 698
 699         int radio;
 700
 701         int i2, j2, k2;
 702         int i3, j3, k3;
 703
 704         ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
 705         ptr2=(unsigned char *)  data2->GetScalarPointer(i,j,k);
 706
 707         int ext[6];
 708
 709         data->GetExtent(ext);
 710
 711         radio=(ext[1]-ext[0])/2;
 712
 713         i2=i-ext[0]-radio;
 714         j2=j-ext[2]-radio;
 715         k2=k-ext[4]-radio;
 716
 717         *ptr2=label;
 718         vector[label-1][0]+=1;
 719         vector[label-1][1]+=i;
 720         vector[label-1][2]+=j;
 721         vector[label-1][3]+=k;
 722
 723         for(i3=-1;i3<=1;i3++){
 724                         for(j3=-1;j3<=1;j3++){
 725                                 for(k3=-1;k3<=1;k3++){
 726                                         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] ){
 727                                                 ptr=(unsigned char *)   data->GetScalarPointer(i+i3,j+j3,k+k3);
 728                                                 ptr2=(unsigned char *)  data2->GetScalarPointer(i+i3,j+j3,k+k3);
 729                                                 if(*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)))){
 730                                                         searcb(i+i3, j+j3, k+k3, data, data2, label, vector);
 731                                                 }
 732                                         }
 733                                 }
 734                         }
 735                 }
 736
 737 }
 738
 739
 740
 741
 742 //----------------------------------------------------------------------------
 743 unsigned char axisExtractor02::find_components(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
 744 {
 745         int ext[6];
 746         int i, j, k, radio;
 747         int i2, j2, k2;
 748
 749         data->GetExtent(ext);
 750         double *ff=  data->GetOrigin();
 751
 752
 753         unsigned char *ptr;
 754         unsigned char *ptr2;
 755
 756         radio=(ext[1]-ext[0])/2;
 757
 758         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
 759                         for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
 760                                 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
 761                                         ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
 762                                         ptr2=(unsigned char *)  data2->GetScalarPointer(i,j,k);
 763
 764                                         if(*ptr==255 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2))  ){
 765                                                 label=label+1;
 766                                                 vector[label-1][0]=0;
 767                                                 vector[label-1][1]=0;
 768                                                 vector[label-1][2]=0;
 769                                                 vector[label-1][3]=0;
 770                                                 searc(i, j, k, data, data2, label, vector);
 771
 772                                         }
 773                                 }
 774                         }
 775                 }
 776
 777                 return label;
 778 }
 779
 780
 781
 782
 783 //----------------------------------------------------------------------------
 784 unsigned char axisExtractor02::find_componentsb(vtkImageData *data, vtkImageData *data2, unsigned char label, unsigned long vector[50][4] )
 785 {
 786         int ext[6];
 787         int i, j, k, radio;
 788         int i2, j2, k2;
 789
 790         data->GetExtent(ext);
 791
 792         unsigned char *ptr;
 793         unsigned char *ptr2;
 794
 795         radio=(ext[1]-ext[0])/2;
 796
 797         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
 798                         for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
 799                                 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
 800                                         ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
 801                                         ptr2=(unsigned char *)  data2->GetScalarPointer(i,j,k);
 802
 803                                         if(*ptr==0 && *ptr2==0 && (radio*radio)>=((i2*i2)+(j2*j2)+(k2*k2)) && ((radio-1)*(radio-1))<=((i2*i2)+(j2*j2)+(k2*k2))  ){
 804                                                 label=label+1;
 805                                                 vector[label-1][0]=0;
 806                                                 vector[label-1][1]=0;
 807                                                 vector[label-1][2]=0;
 808                                                 vector[label-1][3]=0;
 809                                                 searcb(i, j, k, data, data2, label, vector);
 810
 811                                         }
 812                                 }
 813                         }
 814                 }
 815
 816                 return label;
 817 }
 818
 819
 820
 821
 822
 823 //----------------------------------------------------------------------------
 824 int axisExtractor02::proporcion(vtkImageData *data )
 825 {
 826         int ext[6];
 827         int i, j, k;
 828         int max=0;
 829         int total=0;
 830
 831         data->GetExtent(ext);
 832
 833         unsigned short *ptr;
 834
 835         for(i=ext[0];i<=ext[1];i++){
 836                 for(j=ext[2];j<=ext[3];j++){
 837                         for(k=ext[4];k<=ext[5];k++){
 838                                 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
 839                                 if(*ptr!=0){
 840                                         max++;
 841                                 }
 842                                 total++;
 843                         }
 844                 }
 845         }
 846
 847         return (max*100)/total;
 848 }
 849
 850
 851
 852
 853
 854
 855
 856
 857 //----------------------------------------------------------------------------
 858 bool axisExtractor02::border(vtkImageData *data, int p1[3] )
 859 {
 860         int ext[6];
 861         int i, j, k;
 862         short val;
 863         int total1, total2;
 864         int centro[3];
 865
 866
 867         int p2[3];
 868
 869         data->GetExtent(ext);
 870
 871         centro[0]=(ext[1]+ext[0])/2;
 872         centro[1]=(ext[3]+ext[2])/2;
 873         centro[2]=(ext[5]+ext[4])/2;
 874
 875         unsigned char *ptr;
 876
 877         ptr=(unsigned char *)   data->GetScalarPointer(p1);
 878
 879         val=*ptr;
 880
 881         bool res=false;
 882         total1=0;
 883         total2=0;
 884
 885         for(i=-1;i<=1;i++){
 886                 for(j=-1;j<=1;j++){
 887                         for(k=-1;k<=1;k++){
 888                                 p2[0]=p1[0]+i;
 889                                 p2[1]=p1[1]+j;
 890                                 p2[2]=p1[2]+k;
 891                                 if(ext[0]<=p2[0] && ext[1]>=p2[0] && ext[2]<=p2[1] && ext[3]>=p2[1] && ext[4]<=p2[2] && ext[5]>=p2[2] ){
 892                                         ptr=(unsigned char *)   data->GetScalarPointer(p2);
 893                                         if(*ptr!=val){
 894                                                 total1++;
 895                                         }
 896                                         else{
 897                                                 total2++;
 898                                         }
 899                                 }
 900                         }
 901                 }
 902         }
 903
 904         if(total1>0){
 905                 res=true;
 906         }
 907         /*if(p1[0]==centro[0] && p1[1]==centro[1] && p1[2]==centro[2]){
 908                 res=false;
 909         }*/
 910
 911         return res;
 912 }
 913
 914
 915
 916
 917
 918
 919 //----------------------------------------------------------------------------
 920 void axisExtractor02::optim(vtkImageData *data, vtkImageData *data2 )
 921 {
 922         int ext[6];
 923         int i, j, k, radio;
 924         int centro[3];
 925
 926         double w0, w1;
 927
 928         double w0p, w1p;
 929
 930         double inerciaxx, inerciayy, inerciazz, inerciaxy, inerciayz, inerciazx;
 931         double sumx, sumy, sumz;
 932         double sumxx, sumyy, sumzz, sumxy, sumyz, sumzx;
 933
 934         double sumix, sumiy, sumiz;
 935         double sumixx, sumiyy, sumizz, sumixy, sumiyz, sumizx;
 936
 937         double inerciaixx, inerciaiyy, inerciaizz, inerciaixy, inerciaiyz, inerciaizx;
 938
 939         double inerciaxx2, inerciayy2, inerciazz2, inerciaxy2, inerciayz2, inerciazx2;
 940         double sumx2, sumy2, sumz2;
 941         double sumxx2, sumyy2, sumzz2, sumxy2, sumyz2, sumzx2;
 942
 943         double sumix2, sumiy2, sumiz2;
 944         double sumixx2, sumiyy2, sumizz2, sumixy2, sumiyz2, sumizx2;
 945
 946         //double sumixp, sumiyp, sumizp;
 947         //double sumix2p, sumiy2p, sumiz2p;
 948
 949         double inerciaxxp, inerciayyp, inerciazzp, inerciaxyp, inerciayzp, inerciazxp;
 950         double sumxp, sumyp, sumzp;
 951         double sumxxp, sumyyp, sumzzp, sumxyp, sumyzp, sumzxp;
 952
 953         double inerciaxx2p, inerciayy2p, inerciazz2p, inerciaxy2p, inerciayz2p, inerciazx2p;
 954         double sumx2p, sumy2p, sumz2p;
 955         double sumxx2p, sumyy2p, sumzz2p, sumxy2p, sumyz2p, sumzx2p;
 956
 957         //double inerciaixx2p, inerciaiyy2p, inerciaizz2p, inerciaixy2p, inerciaiyz2p, inerciaizx2p;
 958
 959         unsigned long cantp;
 960         double sumip;
 961         double sumiip;
 962         double inerciaiip;
 963         double centip;
 964
 965         unsigned long cant2p;
 966         double sumi2p;
 967         double sumii2p;
 968         double inerciaii2p;
 969         double centi2p;
 970
 971         int point[3];
 972
 973
 974         double sumitotal=0;
 975         double distcent=0;
 976
 977         double im, jm, km, lm;
 978
 979
 980         max2=0;
 981         cant=0;
 982         sumi=0;
 983         sumii=0;
 984         des2=0;
 985         min2=3200;
 986         centi=0;
 987
 988         max3=0;
 989         cant2=0;
 990         sumi2=0;
 991         sumii2=0;
 992         des3=0;
 993         min3=3200;
 994         centi2=0;
 995
 996         sumx=0;
 997         sumy=0;
 998         sumz=0;
 999
1000         sumxx=0;
1001         sumyy=0;
1002         sumzz=0;
1003
1004         sumxy=0;
1005         sumyz=0;
1006         sumzx=0;
1007
1008         sumix=0;
1009         sumiy=0;
1010         sumiz=0;
1011
1012         sumixx=0;
1013         sumiyy=0;
1014         sumizz=0;
1015
1016         sumixy=0;
1017         sumiyz=0;
1018         sumizx=0;
1019
1020         sumx2=0;
1021         sumy2=0;
1022         sumz2=0;
1023
1024         sumxx2=0;
1025         sumyy2=0;
1026         sumzz2=0;
1027
1028         sumxy2=0;
1029         sumyz2=0;
1030         sumzx2=0;
1031
1032         sumix2=0;
1033         sumiy2=0;
1034         sumiz2=0;
1035
1036         sumixx2=0;
1037         sumiyy2=0;
1038         sumizz2=0;
1039
1040         sumixy2=0;
1041         sumiyz2=0;
1042         sumizx2=0;
1043
1044         max=0;
1045         min=3200;
1046
1047
1048         data->GetExtent(ext);
1049
1050         unsigned short *ptr;
1051         unsigned char *ptr2;
1052
1053         radio=(ext[1]-ext[0])/2;
1054
1055         centro[0]=(ext[1]+ext[0])/2;
1056         centro[1]=(ext[3]+ext[2])/2;
1057         centro[2]=(ext[5]+ext[4])/2;
1058
1059         double tmpsqrt;
1060
1061         for(i=ext[0];i<=ext[1];i++){
1062                 for(j=ext[2];j<=ext[3];j++){
1063                         for(k=ext[4];k<=ext[5];k++){
1064                                 tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
1065                                 if( radio>=sqrt( tmpsqrt) ){
1066                                         ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1067                                         if(*ptr>max){
1068                                                 max=*ptr;
1069
1070                                         }
1071                                         if(*ptr<min){
1072                                                 min=*ptr;
1073
1074                                         }
1075                                 }
1076                         }
1077                 }
1078         }
1079
1080         for(i=ext[0];i<=ext[1];i++){
1081                 for(j=ext[2];j<=ext[3];j++){
1082                         for(k=ext[4];k<=ext[5];k++){
1083                                 tmpsqrt=((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
1084                                 if( radio>=sqrt(tmpsqrt) ){
1085                                         ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1086                                         ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
1087                                         point[0]=i;
1088                                         point[1]=j;
1089                                         point[2]=k;
1090
1091                                         im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
1092                                         jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
1093                                         km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
1094                                         lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
1095
1096                                         if(*ptr2!=0){
1097                                                 sumi+=lm;
1098                                                 sumii+=lm*lm;
1099                                                 cant+=1;
1100
1101                                                 sumx+=im;
1102                                                 sumy+=jm;
1103                                                 sumz+=km;
1104
1105                                                 sumxx+=im*im;
1106                                                 sumyy+=jm*jm;
1107                                                 sumzz+=km*km;
1108                                                 sumxy+=im*jm;
1109                                                 sumyz+=jm*km;
1110                                                 sumzx+=km*im;
1111
1112                                                 sumix+=lm*im;
1113                                                 sumiy+=lm*jm;
1114                                                 sumiz+=lm*km;
1115
1116                                                 sumixx+=lm*im*im;
1117                                                 sumiyy+=lm*jm*jm;
1118                                                 sumizz+=lm*km*km;
1119                                                 sumixy+=lm*im*jm;
1120                                                 sumiyz+=lm*jm*km;
1121                                                 sumizx+=lm*km*im;
1122                                         }
1123                                         else{
1124                                                 sumi2+=lm;
1125                                                 sumii2+=lm*lm;
1126                                                 cant2+=1;
1127
1128                                                 sumx2+=im;
1129                                                 sumy2+=jm;
1130                                                 sumz2+=km;
1131
1132                                                 sumxx2+=im*im;
1133                                                 sumyy2+=jm*jm;
1134                                                 sumzz2+=km*km;
1135                                                 sumxy2+=im*jm;
1136                                                 sumyz2+=jm*km;
1137                                                 sumzx2+=km*im;
1138
1139                                                 sumix2+=lm*im;
1140                                                 sumiy2+=lm*jm;
1141                                                 sumiz2+=lm*km;
1142
1143                                                 sumixx2+=lm*im*im;
1144                                                 sumiyy2+=lm*jm*jm;
1145                                                 sumizz2+=lm*km*km;
1146                                                 sumixy2+=lm*im*jm;
1147                                                 sumiyz2+=lm*jm*km;
1148                                                 sumizx2+=lm*km*im;
1149                                         }
1150                                 }
1151                         }
1152                 }
1153         }
1154
1155         centi=(double)sumi/(double)cant;
1156         inerciaii= ((double)sumii/(double)cant)-(centi*centi);
1157
1158         centi2=(double)sumi2/(double)cant2;
1159         inerciaii2= ((double)sumii2/(double)cant2)-(centi2*centi2);
1160
1161         centx=(double)sumx/(double)cant;
1162         centy=(double)sumy/(double)cant;
1163         centz=(double)sumz/(double)cant;
1164
1165         centx2=(double)sumx2/(double)cant2;
1166         centy2=(double)sumy2/(double)cant2;
1167         centz2=(double)sumz2/(double)cant2;
1168
1169         inerciaxx= ((double)sumxx/(double)cant)-centx*centx;
1170         inerciayy= ((double)sumyy/(double)cant)-centy*centy;
1171         inerciazz= ((double)sumzz/(double)cant)-centz*centz;
1172
1173         inerciaxx2= ((double)sumxx2/(double)cant2)-centx2*centx2;
1174         inerciayy2= ((double)sumyy2/(double)cant2)-centy2*centy2;
1175         inerciazz2= ((double)sumzz2/(double)cant2)-centz2*centz2;
1176
1177         inerciaxy= ((double)sumxy/(double)cant)-centx*centy;
1178         inerciayz= ((double)sumyz/(double)cant)-centy*centz;
1179         inerciazx= ((double)sumzx/(double)cant)-centz*centx;
1180
1181         inerciaxy2= ((double)sumxy2/(double)cant2)-centx2*centy2;
1182         inerciayz2= ((double)sumyz2/(double)cant2)-centy2*centz2;
1183         inerciazx2= ((double)sumzx2/(double)cant2)-centz2*centx2;
1184
1185
1186         centix=((double)sumix+(double)sumix2)/((double)sumi+(double)sumi2);
1187         centiy=((double)sumiy+(double)sumiy2)/((double)sumi+(double)sumi2);
1188         centiz=((double)sumiz+(double)sumiz2)/((double)sumi+(double)sumi2);
1189
1190         inerciaixx= ( ((double)sumixx+(double)sumixx2) / ((double)sumi+(double)sumi2)  )-(centix*centix);
1191         inerciaiyy= ( ((double)sumiyy+(double)sumiyy2) / ((double)sumi+(double)sumi2)  )-(centiy*centiy);
1192         inerciaizz= ( ((double)sumizz+(double)sumizz2) / ((double)sumi+(double)sumi2)  )-(centiz*centiz);
1193
1194
1195         inerciaixy= ( ((double)sumixy+(double)sumixy2) / ((double)sumi+(double)sumi2)  )-(centix*centiy);
1196         inerciaiyz= ( ((double)sumiyz+(double)sumiyz2) / ((double)sumi+(double)sumi2)  )-(centiy*centiz);
1197         inerciaizx= ( ((double)sumizx+(double)sumizx2) / ((double)sumi+(double)sumi2)  )-(centiz*centix);
1198
1199
1200         A[0][0]=inerciaxx;
1201         A[1][1]=inerciayy;
1202         A[2][2]=inerciazz;
1203         A[0][1]=A[1][0]=inerciaxy;
1204         A[0][2]=A[2][0]=inerciazx;
1205         A[1][2]=A[2][1]=inerciayz;
1206
1207         A2[0][0]=inerciaxx2;
1208         A2[1][1]=inerciayy2;
1209         A2[2][2]=inerciazz2;
1210         A2[0][1]=A2[1][0]=inerciaxy2;
1211         A2[0][2]=A2[2][0]=inerciazx2;
1212         A2[1][2]=A2[2][1]=inerciayz2;
1213
1214         Ai[0][0]=inerciaixx;
1215         Ai[1][1]=inerciaiyy;
1216         Ai[2][2]=inerciaizz;
1217         Ai[0][1]=Ai[1][0]=inerciaixy;
1218         Ai[0][2]=Ai[2][0]=inerciaizx;
1219         Ai[1][2]=Ai[2][1]=inerciaiyz;
1220
1221         vtkMath::Diagonalize3x3  (A, w, V);
1222         vtkMath::Diagonalize3x3  (A2, w2, V2);
1223         vtkMath::Diagonalize3x3  (Ai, wi, Vi);
1224
1225
1226         if(w[0]>w[1]){
1227                 if(w[0]>w[2]){
1228                         ejemin=0;
1229                         inerciar=w[1]+w[2];
1230                         if(w[1]>w[2]){
1231                                 inerciary=w[0]+w[2];
1232                                 inerciarz=w[0]+w[1];
1233                         }
1234                         else{
1235                                 inerciary=w[0]+w[1];
1236                                 inerciarz=w[0]+w[2];
1237                         }
1238                 }
1239                 else{
1240                         ejemin=2;
1241                         inerciar=w[0]+w[1];
1242                         if(w[1]>w[0]){
1243                                 inerciary=w[0]+w[2];
1244                                 inerciarz=w[1]+w[2];
1245                         }
1246                         else{
1247                                 inerciary=w[1]+w[2];
1248                                 inerciarz=w[0]+w[2];
1249                         }
1250                 }
1251         }
1252         else{
1253                 if(w[1]>w[2]){
1254                         ejemin=1;
1255                         inerciar=w[2]+w[0];
1256                         if(w[2]>w[0]){
1257                                 inerciary=w[1]+w[0];
1258                                 inerciarz=w[1]+w[2];
1259                         }
1260                         else{
1261                                 inerciary=w[1]+w[2];
1262                                 inerciarz=w[1]+w[0];
1263                         }
1264                 }
1265                 else{
1266                         ejemin=2;
1267                         inerciar=w[1]+w[0];
1268                         if(w[1]>w[0]){
1269                                 inerciary=w[0]+w[2];
1270                                 inerciarz=w[1]+w[2];
1271                         }
1272                         else{
1273                                 inerciary=w[1]+w[2];
1274                                 inerciarz=w[0]+w[2];
1275                         }
1276                 }
1277         }
1278
1279
1280         if(w2[0]>w2[1]){
1281                 if(w2[0]>w2[2]){
1282                         ejemin2=0;
1283                         inercia2r=w2[1]+w2[2];
1284                 }
1285                 else{
1286                         ejemin2=2;
1287                         inercia2r=w2[1]+w2[0];
1288                 }
1289         }
1290         else{
1291                 if(w2[1]>w2[2]){
1292                         ejemin2=1;
1293                         inercia2r=w2[2]+w2[0];
1294                 }
1295                 else{
1296                         ejemin2=2;
1297                         inercia2r=w2[1]+w2[0];
1298                 }
1299         }
1300
1301
1302         if(wi[0]>wi[1]){
1303                 if(wi[0]>wi[2]){
1304                         ejemini=0;
1305                         inerciari=wi[1]+wi[2];
1306                         if(wi[1]>wi[2]){
1307                                 inerciariy=wi[0]+wi[2];
1308                                 inerciariz=wi[0]+wi[1];
1309                         }
1310                         else{
1311                                 inerciariy=wi[0]+wi[1];
1312                                 inerciariz=wi[0]+wi[2];
1313                         }
1314                 }
1315                 else{
1316                         ejemini=2;
1317                         inerciari=wi[1]+wi[0];
1318                         if(wi[1]>wi[0]){
1319                                 inerciariy=wi[0]+wi[2];
1320                                 inerciariz=wi[1]+wi[2];
1321                         }
1322                         else{
1323                                 inerciariy=wi[1]+wi[2];
1324                                 inerciariz=wi[0]+wi[2];
1325                         }
1326
1327                 }
1328         }
1329         else{
1330                 if(wi[1]>wi[2]){
1331                         ejemini=1;
1332                         inerciari=wi[2]+wi[0];
1333                         if(wi[2]>wi[0]){
1334                                 inerciariy=wi[1]+wi[0];
1335                                 inerciariz=wi[1]+wi[2];
1336                         }
1337                         else{
1338                                 inerciariy=wi[1]+wi[2];
1339                                 inerciariz=wi[1]+wi[0];
1340                         }
1341                 }
1342                 else{
1343                         ejemini=2;
1344                         inerciari=wi[1]+wi[0];
1345                         if(wi[1]>wi[0]){
1346                                 inerciariy=wi[0]+wi[2];
1347                                 inerciariz=wi[1]+wi[2];
1348                         }
1349                         else{
1350                                 inerciariy=wi[1]+wi[2];
1351                                 inerciariz=wi[0]+wi[2];
1352                         }
1353                 }
1354         }
1355
1356
1357         /*
1358         for(i=ext[0];i<=ext[1];i++){
1359                 for(j=ext[2];j<=ext[3];j++){
1360                         for(k=ext[4];k<=ext[5];k++){
1361                                 if( radio>=sqrt(((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]))) ){
1362                                         ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1363                                         ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
1364
1365
1366                                         im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
1367                                         jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
1368                                         km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
1369
1370                                         point1[0]=im;
1371                                         point1[1]=jm;
1372                                         point1[2]=km;
1373
1374                                         point2[0]=centx+V[0][ejemin];
1375                                         point2[1]=centy+V[1][ejemin];
1376                                         point2[2]=centz+V[2][ejemin];
1377
1378                                         point3[0]=centx-V[0][ejemin];
1379                                         point3[1]=centy-V[1][ejemin];
1380                                         point3[2]=centz-V[2][ejemin];
1381
1382                                         //distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
1383                                         distcent=1;
1384
1385                                         if(*ptr2==0){
1386
1387                                                 sumitotal+=(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
1388
1389                                         }
1390                                 }
1391                         }
1392                 }
1393         }
1394
1395         inercia2r=sumitotal/cant2;
1396         */
1397
1398
1399
1400
1401         /*printf("centx: %f\n", centx);
1402         printf("centy: %f\n", centy);
1403         printf("centz: %f\n", centz);
1404
1405         printf("w[0]: %f\n", w[0]);
1406         printf("w[1]: %f\n", w[1]);
1407         printf("w[2]: %f\n", w[2]);
1408         printf("inerciar: %f\n", inerciar);
1409
1410         printf("w[ejemin]: %f\n", w[ejemin]);
1411
1412         printf("w[ejemin]/inerciar: %f\n", w[ejemin]/inerciar);
1413
1414         printf("V[0][0]: %f\n", V[0][0]);
1415         printf("V[1][0]: %f\n", V[1][0]);
1416         printf("V[2][0]: %f\n", V[2][0]);
1417
1418         printf("V[0][1]: %f\n", V[0][1]);
1419         printf("V[1][1]: %f\n", V[1][1]);
1420         printf("V[2][1]: %f\n", V[2][1]);
1421
1422         printf("V[0][2]: %f\n", V[0][2]);
1423         printf("V[1][2]: %f\n", V[1][2]);
1424         printf("V[2][2]: %f\n", V[2][2]);*/
1425
1426         w0=(double)cant/((double)cant+(double)cant2);
1427         w1=(double)cant2/((double)cant+(double)cant2);
1428
1429
1430         /*
1431         ut=w0*centi+w1*centi2;
1432         intervar=w0*inerciaii+w1*inerciaii2;
1433         */
1434
1435         //w0=1;
1436         //w1=1;
1437
1438
1439
1440         //costo=w0*inerciaii+w1*inerciaii2+param*inerciar;
1441         //costo=param*(w0*param3*inerciaii+w1*(1-param3)*inerciaii2)+(1-param)*(w0*param2*inerciar+w1*(1-param2)*inercia2r);
1442         costo=param*w0*inerciaii+w1*param2*inerciaii2+w0*param3*inerciar+w1*param4*inercia2r;
1443
1444         int changes=4;
1445
1446         int total=0;
1447
1448         while(changes>3 && total<500){
1449                 changes=0;
1450                 total++;
1451                 for(i=ext[0];i<=ext[1];i++){
1452                         for(j=ext[2];j<=ext[3];j++){
1453                                 for(k=ext[4];k<=ext[5];k++){
1454                                         tmpsqrt = ((i-centro[0])*(i-centro[0]))+((j-centro[1])*(j-centro[1]))+((k-centro[2])*(k-centro[2]));
1455                                         if( radio>=sqrt(tmpsqrt) ){
1456                                                 ptr=(unsigned short *)data->GetScalarPointer(i,j,k);
1457                                                 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
1458                                                 point[0]=i;
1459                                                 point[1]=j;
1460                                                 point[2]=k;
1461
1462                                                 im=2*((((double)i-(double)ext[0])/((double)ext[1]-(double)ext[0]))-0.5);
1463                                                 jm=2*((((double)j-(double)ext[2])/((double)ext[3]-(double)ext[2]))-0.5);
1464                                                 km=2*((((double)k-(double)ext[4])/((double)ext[5]-(double)ext[4]))-0.5);
1465                                                 lm=2*(((double)*ptr-(double)min)/((double)max-(double)min));
1466
1467                                                 if(border(data2, point)){
1468                                                         if(*ptr2!=0){
1469
1470                                                                 sumip=sumi-lm;
1471                                                                 sumiip=sumii-lm*lm;
1472
1473                                                                 cantp=cant-1;
1474
1475                                                                 sumi2p=sumi2+lm;
1476                                                                 sumii2p=sumii2+lm*lm;
1477
1478                                                                 cant2p=cant2+1;
1479
1480                                                                 sumxp=sumx-im;
1481                                                                 sumyp=sumy-jm;
1482                                                                 sumzp=sumz-km;
1483                                                                 sumxxp=sumxx-im*im;
1484                                                                 sumyyp=sumyy-jm*jm;
1485                                                                 sumzzp=sumzz-km*km;
1486                                                                 sumxyp=sumxy-im*jm;
1487                                                                 sumyzp=sumyz-jm*km;
1488                                                                 sumzxp=sumzx-km*im;
1489
1490                                                                 sumx2p=sumx2+im;
1491                                                                 sumy2p=sumy2+jm;
1492                                                                 sumz2p=sumz2+km;
1493                                                                 sumxx2p=sumxx2+im*im;
1494                                                                 sumyy2p=sumyy2+jm*jm;
1495                                                                 sumzz2p=sumzz2+km*km;
1496                                                                 sumxy2p=sumxy2+im*jm;
1497                                                                 sumyz2p=sumyz2+jm*km;
1498                                                                 sumzx2p=sumzx2+km*im;
1499
1500                                                                 centip=(double)sumip/(double)cantp;
1501                                                                 inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
1502
1503                                                                 centi2p=(double)sumi2p/(double)cant2p;
1504                                                                 inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
1505
1506                                                                 centxp=(double)sumxp/(double)cantp;
1507                                                                 centyp=(double)sumyp/(double)cantp;
1508                                                                 centzp=(double)sumzp/(double)cantp;
1509
1510                                                                 centx2p=(double)sumx2p/(double)cant2p;
1511                                                                 centy2p=(double)sumy2p/(double)cant2p;
1512                                                                 centz2p=(double)sumz2p/(double)cant2p;
1513
1514                                                                 inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
1515                                                                 inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
1516                                                                 inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;
1517
1518                                                                 inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
1519                                                                 inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
1520                                                                 inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;
1521
1522                                                                 inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
1523                                                                 inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
1524                                                                 inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;
1525
1526                                                                 inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
1527                                                                 inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
1528                                                                 inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;
1529
1530                                                                 Ap[0][0]=inerciaxxp;
1531                                                                 Ap[1][1]=inerciayyp;
1532                                                                 Ap[2][2]=inerciazzp;
1533                                                                 Ap[0][1]=Ap[1][0]=inerciaxyp;
1534                                                                 Ap[0][2]=Ap[2][0]=inerciazxp;
1535                                                                 Ap[1][2]=Ap[2][1]=inerciayzp;
1536
1537
1538                                                                 A2p[0][0]=inerciaxx2p;
1539                                                                 A2p[1][1]=inerciayy2p;
1540                                                                 A2p[2][2]=inerciazz2p;
1541                                                                 A2p[0][1]=A2p[1][0]=inerciaxy2p;
1542                                                                 A2p[0][2]=A2p[2][0]=inerciazx2p;
1543                                                                 A2p[1][2]=A2p[2][1]=inerciayz2p;
1544
1545
1546
1547                                                                 w0p=(double)cantp/((double)cantp+(double)cant2p);
1548                                                                 w1p=(double)cant2p/((double)cantp+(double)cant2p);
1549
1550                                                                 //w0p=1;
1551                                                                 //w1p=1;
1552
1553                                                                 vtkMath::Diagonalize3x3  (Ap, wp, Vp);
1554                                                                 vtkMath::Diagonalize3x3  (A2p, w2p, V2p);
1555
1556
1557                                                                 if(wp[0]>wp[1]){
1558                                                                         if(wp[0]>wp[2]){
1559                                                                                 ejeminp=0;
1560                                                                                 inerciarp=wp[2]+wp[1];
1561                                                                                 if(wp[1]>wp[2]){
1562                                                                                         inerciarpy=wp[0]+wp[2];
1563                                                                                         inerciarpz=wp[0]+wp[1];
1564                                                                                 }
1565                                                                                 else{
1566                                                                                         inerciarpy=wp[0]+wp[1];
1567                                                                                         inerciarpz=wp[0]+wp[2];
1568                                                                                 }
1569                                                                         }
1570                                                                         else{
1571                                                                                 ejeminp=2;
1572                                                                                 inerciarp=wp[1]+wp[0];
1573                                                                                 if(wp[1]>wp[0]){
1574                                                                                         inerciarpy=wp[0]+wp[2];
1575                                                                                         inerciarpz=wp[1]+wp[2];
1576                                                                                 }
1577                                                                                 else{
1578                                                                                         inerciarpy=wp[1]+wp[2];
1579                                                                                         inerciarpz=wp[0]+wp[2];
1580                                                                                 }
1581                                                                         }
1582                                                                 }
1583                                                                 else{
1584                                                                         if(wp[1]>wp[2]){
1585                                                                                 ejeminp=1;
1586                                                                                 inerciarp=wp[2]+wp[0];
1587                                                                                 if(wp[2]>wp[0]){
1588                                                                                         inerciarpy=wp[1]+wp[0];
1589                                                                                         inerciarpz=wp[1]+wp[2];
1590                                                                                 }
1591                                                                                 else{
1592                                                                                         inerciarpy=wp[1]+wp[2];
1593                                                                                         inerciarpz=wp[1]+wp[0];
1594                                                                                 }
1595                                                                         }
1596                                                                         else{
1597                                                                                 ejeminp=2;
1598                                                                                 inerciarp=wp[1]+wp[0];
1599                                                                                 if(wp[1]>wp[0]){
1600                                                                                         inerciarpy=wp[0]+wp[2];
1601                                                                                         inerciarpz=wp[1]+wp[2];
1602                                                                                 }
1603                                                                                 else{
1604                                                                                         inerciarpy=wp[1]+wp[2];
1605                                                                                         inerciarpz=wp[0]+wp[2];
1606                                                                                 }
1607                                                                         }
1608                                                                 }
1609
1610
1611                                                                 if(w2p[0]>w2p[1]){
1612                                                                         if(w2p[0]>w2p[2]){
1613                                                                                 ejemin2p=0;
1614                                                                                 inercia2rp=w2p[2]+w2p[1];
1615                                                                         }
1616                                                                         else{
1617                                                                                 ejemin2p=2;
1618                                                                                 inercia2rp=w2p[1]+w2p[0];
1619                                                                         }
1620                                                                 }
1621                                                                 else{
1622                                                                         if(w2p[1]>w2p[2]){
1623                                                                                 ejemin2p=1;
1624                                                                                 inercia2rp=w2p[2]+w2p[0];
1625                                                                         }
1626                                                                         else{
1627                                                                                 ejemin2p=2;
1628                                                                                 inercia2rp=w2p[1]+w2p[0];
1629                                                                         }
1630                                                                 }
1631
1632                                                         /*      point1[0]=im;
1633                                                                 point1[1]=jm;
1634                                                                 point1[2]=km;
1635
1636                                                                 point2[0]=centxp+Vp[0][ejeminp];
1637                                                                 point2[1]=centyp+Vp[1][ejeminp];
1638                                                                 point2[2]=centzp+Vp[2][ejeminp];
1639
1640                                                                 point3[0]=centxp-Vp[0][ejeminp];
1641                                                                 point3[1]=centyp-Vp[1][ejeminp];
1642                                                                 point3[2]=centzp-Vp[2][ejeminp];*/
1643
1644                                                                 //distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
1645                                                         /*      distcent=1;
1646
1647                                                                 sumitotalp=sumitotal+(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
1648
1649                                                                 inercia2rp=sumitotalp/cant2p;*/
1650
1651                                                                 //costop=w0p*inerciaiip+w1p*inerciaii2p+param*inerciarp;
1652
1653                                                                 //costop=param*(w0p*param3*inerciaiip+w1p*(1-param3)*inerciaii2p)+(1-param)*(w0p*param2*inerciarp+w1p*(1-param2)*inercia2rp);
1654
1655
1656                                                                 costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;
1657
1658                                                                 if(costop<costo  && centip>centi2p){
1659
1660                                                                         *ptr2=0;
1661                                                                         changes++;
1662
1663                                                                         sumi=sumip;
1664                                                                         sumii=sumiip;
1665                                                                         cant=cantp;
1666
1667                                                                         sumi2=sumi2p;
1668                                                                         sumii2=sumii2p;
1669                                                                         cant2=cant2p;
1670
1671                                                                         sumx=sumxp;
1672                                                                         sumy=sumyp;
1673                                                                         sumz=sumzp;
1674                                                                         sumxx=sumxxp;
1675                                                                         sumyy=sumyyp;
1676                                                                         sumzz=sumzzp;
1677                                                                         sumxy=sumxyp;
1678                                                                         sumyz=sumyzp;
1679                                                                         sumzx=sumzxp;
1680
1681                                                                         sumx2=sumx2p;
1682                                                                         sumy2=sumy2p;
1683                                                                         sumz2=sumz2p;
1684                                                                         sumxx2=sumxx2p;
1685                                                                         sumyy2=sumyy2p;
1686                                                                         sumzz2=sumzz2p;
1687                                                                         sumxy2=sumxy2p;
1688                                                                         sumyz2=sumyz2p;
1689                                                                         sumzx2=sumzx2p;
1690
1691                                                                         centi=centip;
1692                                                                         inerciaii= inerciaiip;
1693
1694                                                                         centi2=centi2p;
1695                                                                         inerciaii2p= inerciaii2p;
1696
1697                                                                         centx=centxp;
1698                                                                         centy=centyp;
1699                                                                         centz=centzp;
1700
1701                                                                         centx2=centx2p;
1702                                                                         centy2=centy2p;
1703                                                                         centz2=centz2p;
1704
1705                                                                         inerciaxx= inerciaxxp;
1706                                                                         inerciayy= inerciayyp;
1707                                                                         inerciazz= inerciazzp;
1708
1709                                                                         inerciaxx2= inerciaxx2p;
1710                                                                         inerciayy2= inerciayy2p;
1711                                                                         inerciazz2= inerciazz2p;
1712
1713                                                                         inerciaxy= inerciaxyp;
1714                                                                         inerciayz= inerciayzp;
1715                                                                         inerciazx= inerciazxp;
1716
1717                                                                         inerciaxy2= inerciaxy2p;
1718                                                                         inerciayz2= inerciayz2p;
1719                                                                         inerciazx2= inerciazx2p;
1720
1721                                                                 //      sumitotal=sumitotalp;
1722
1723                                                                         costo=costop;
1724
1725                                                                 }
1726                                                         }
1727                                                         else{
1728                                                                 sumip=sumi+lm;
1729                                                                 sumiip=sumii+lm*lm;
1730
1731                                                                 cantp=cant+1;
1732
1733                                                                 sumi2p=sumi2-lm;
1734                                                                 sumii2p=sumii2-lm*lm;
1735
1736                                                                 cant2p=cant2-1;
1737
1738                                                                 sumxp=sumx+im;
1739                                                                 sumyp=sumy+jm;
1740                                                                 sumzp=sumz+km;
1741                                                                 sumxxp=sumxx+im*im;
1742                                                                 sumyyp=sumyy+jm*jm;
1743                                                                 sumzzp=sumzz+km*km;
1744                                                                 sumxyp=sumxy+im*jm;
1745                                                                 sumyzp=sumyz+jm*km;
1746                                                                 sumzxp=sumzx+km*im;
1747
1748                                                                 sumx2p=sumx2-im;
1749                                                                 sumy2p=sumy2-jm;
1750                                                                 sumz2p=sumz2-km;
1751                                                                 sumxx2p=sumxx2-im*im;
1752                                                                 sumyy2p=sumyy2-jm*jm;
1753                                                                 sumzz2p=sumzz2-km*km;
1754                                                                 sumxy2p=sumxy2-im*jm;
1755                                                                 sumyz2p=sumyz2-jm*km;
1756                                                                 sumzx2p=sumzx2-km*im;
1757
1758                                                                 centip=(double)sumip/(double)cantp;
1759                                                                 inerciaiip= ((double)sumiip/(double)cantp)-(centip*centip);
1760
1761                                                                 centi2p=(double)sumi2p/(double)cant2p;
1762                                                                 inerciaii2p= ((double)sumii2p/(double)cant2p)-(centi2p*centi2p);
1763
1764                                                                 centxp=(double)sumxp/(double)cantp;
1765                                                                 centyp=(double)sumyp/(double)cantp;
1766                                                                 centzp=(double)sumzp/(double)cantp;
1767
1768                                                                 centx2p=(double)sumx2p/(double)cant2p;
1769                                                                 centy2p=(double)sumy2p/(double)cant2p;
1770                                                                 centz2p=(double)sumz2p/(double)cant2p;
1771
1772
1773                                                                 inerciaxxp= ((double)sumxxp/(double)cantp)-centxp*centxp;
1774                                                                 inerciayyp= ((double)sumyyp/(double)cantp)-centyp*centyp;
1775                                                                 inerciazzp= ((double)sumzzp/(double)cantp)-centzp*centzp;
1776
1777                                                                 inerciaxx2p= ((double)sumxx2p/(double)cant2p)-centx2p*centx2p;
1778                                                                 inerciayy2p= ((double)sumyy2p/(double)cant2p)-centy2p*centy2p;
1779                                                                 inerciazz2p= ((double)sumzz2p/(double)cant2p)-centz2p*centz2p;
1780
1781                                                                 inerciaxyp= ((double)sumxyp/(double)cantp)-centxp*centyp;
1782                                                                 inerciayzp= ((double)sumyzp/(double)cantp)-centyp*centzp;
1783                                                                 inerciazxp= ((double)sumzxp/(double)cantp)-centzp*centxp;
1784
1785                                                                 inerciaxy2p= ((double)sumxy2p/(double)cant2p)-centx2p*centy2p;
1786                                                                 inerciayz2p= ((double)sumyz2p/(double)cant2p)-centy2p*centz2p;
1787                                                                 inerciazx2p= ((double)sumzx2p/(double)cant2p)-centz2p*centx2p;
1788
1789                                                                 Ap[0][0]=inerciaxxp;
1790                                                                 Ap[1][1]=inerciayyp;
1791                                                                 Ap[2][2]=inerciazzp;
1792                                                                 Ap[0][1]=Ap[1][0]=inerciaxyp;
1793                                                                 Ap[0][2]=Ap[2][0]=inerciazxp;
1794                                                                 Ap[1][2]=Ap[2][1]=inerciayzp;
1795
1796                                                                 A2p[0][0]=inerciaxx2p;
1797                                                                 A2p[1][1]=inerciayy2p;
1798                                                                 A2p[2][2]=inerciazz2p;
1799                                                                 A2p[0][1]=A2p[1][0]=inerciaxy2p;
1800                                                                 A2p[0][2]=A2p[2][0]=inerciazx2p;
1801                                                                 A2p[1][2]=A2p[2][1]=inerciayz2p;
1802
1803                                                                 w0p=(double)cantp/((double)cantp+(double)cant2p);
1804                                                                 w1p=(double)cant2p/((double)cantp+(double)cant2p);
1805
1806                                                                 //w0p=1;
1807                                                                 //w1p=1;
1808
1809                                                                 vtkMath::Diagonalize3x3  (Ap, wp, Vp);
1810                                                                 vtkMath::Diagonalize3x3  (A2p, w2p, V2p);
1811
1812                                                                 if(wp[0]>wp[1]){
1813                                                                         if(wp[0]>wp[2]){
1814                                                                                 ejeminp=0;
1815                                                                                 inerciarp=wp[2]+wp[1];
1816                                                                                 if(wp[1]>wp[2]){
1817                                                                                         inerciarpy=wp[0]+wp[2];
1818                                                                                         inerciarpz=wp[0]+wp[1];
1819                                                                                 }
1820                                                                                 else{
1821                                                                                         inerciarpy=wp[0]+wp[1];
1822                                                                                         inerciarpz=wp[0]+wp[2];
1823                                                                                 }
1824                                                                         }
1825                                                                         else{
1826                                                                                 ejeminp=2;
1827                                                                                 inerciarp=wp[1]+wp[0];
1828                                                                                 if(wp[1]>wp[0]){
1829                                                                                         inerciarpy=wp[0]+wp[2];
1830                                                                                         inerciarpz=wp[1]+wp[2];
1831                                                                                 }
1832                                                                                 else{
1833                                                                                         inerciarpy=wp[1]+wp[2];
1834                                                                                         inerciarpz=wp[0]+wp[2];
1835                                                                                 }
1836                                                                         }
1837                                                                 }
1838                                                                 else{
1839                                                                         if(wp[1]>wp[2]){
1840                                                                                 ejeminp=1;
1841                                                                                 inerciarp=wp[2]+wp[0];
1842                                                                                 if(wp[2]>wp[0]){
1843                                                                                         inerciarpy=wp[1]+wp[0];
1844                                                                                         inerciarpz=wp[1]+wp[2];
1845                                                                                 }
1846                                                                                 else{
1847                                                                                         inerciarpy=wp[1]+wp[2];
1848                                                                                         inerciarpz=wp[1]+wp[0];
1849                                                                                 }
1850                                                                         }
1851                                                                         else{
1852                                                                                 ejeminp=2;
1853                                                                                 inerciarp=wp[1]+wp[0];
1854                                                                                 if(wp[1]>wp[0]){
1855                                                                                         inerciarpy=wp[0]+wp[2];
1856                                                                                         inerciarpz=wp[1]+wp[2];
1857                                                                                 }
1858                                                                                 else{
1859                                                                                         inerciarpy=wp[1]+wp[2];
1860                                                                                         inerciarpz=wp[0]+wp[2];
1861                                                                                 }
1862                                                                         }
1863                                                                 }
1864
1865
1866                                                                 if(w2p[0]>w2p[1]){
1867                                                                         if(w2p[0]>w2p[2]){
1868                                                                                 ejemin2p=0;
1869                                                                                 inercia2rp=w2p[2]+w2p[1];
1870                                                                         }
1871                                                                         else{
1872                                                                                 ejemin2p=2;
1873                                                                                 inercia2rp=w2p[1]+w2p[0];
1874                                                                         }
1875                                                                 }
1876                                                                 else{
1877                                                                         if(w2p[1]>w2p[2]){
1878                                                                                 ejemin2p=1;
1879                                                                                 inercia2rp=w2p[2]+w2p[0];
1880                                                                         }
1881                                                                         else{
1882                                                                                 ejemin2p=2;
1883                                                                                 inercia2rp=w2p[1]+w2p[0];
1884                                                                         }
1885                                                                 }
1886
1887
1888                                                         /*      point1[0]=im;
1889                                                                 point1[1]=jm;
1890                                                                 point1[2]=km;
1891
1892                                                                 point2[0]=centxp+Vp[0][ejeminp];
1893                                                                 point2[1]=centyp+Vp[1][ejeminp];
1894                                                                 point2[2]=centzp+Vp[2][ejeminp];
1895
1896                                                                 point3[0]=centxp-Vp[0][ejeminp];
1897                                                                 point3[1]=centyp-Vp[1][ejeminp];
1898                                                                 point3[2]=centzp-Vp[2][ejeminp];*/
1899
1900                                                                 //distcent=1+sqrt(centx*centx+centy*centy+centz*centz);
1901                                                         /*      distcent=1;
1902
1903                                                                 sumitotalp=sumitotal-(distcent-fabs(distanciaejepunto(point1, point2, point3)))*(distcent-fabs(distanciaejepunto(point1, point2, point3)));
1904
1905                                                                 inercia2rp=sumitotalp/cant2p;*/
1906
1907
1908                                                                 //costop=w0p*inerciaiip+w1p*inerciaii2p+param*inerciarp;
1909                                                                 //costop=param*(w0p*param3*inerciaiip+w1p*(1-param3)*inerciaii2p)+(1-param)*(w0p*param2*inerciarp+w1p*(1-param2)*inercia2rp);
1910
1911
1912                                                                 costop=param*w0p*inerciaiip+w1p*param2*inerciaii2p+w0p*param3*inerciarp+w1p*param4*inercia2rp;
1913
1914                                                                 if(costop<costo && centip>centi2p){
1915
1916                                                                         *ptr2=255;
1917                                                                         changes++;
1918
1919                                                                         sumi=sumip;
1920                                                                         sumii=sumiip;
1921                                                                         cant=cantp;
1922
1923                                                                         sumi2=sumi2p;
1924                                                                         sumii2=sumii2p;
1925                                                                         cant2=cant2p;
1926
1927                                                                         sumx=sumxp;
1928                                                                         sumy=sumyp;
1929                                                                         sumz=sumzp;
1930                                                                         sumxx=sumxxp;
1931                                                                         sumyy=sumyyp;
1932                                                                         sumzz=sumzzp;
1933                                                                         sumxy=sumxyp;
1934                                                                         sumyz=sumyzp;
1935                                                                         sumzx=sumzxp;
1936
1937                                                                         sumx2=sumx2p;
1938                                                                         sumy2=sumy2p;
1939                                                                         sumz2=sumz2p;
1940                                                                         sumxx2=sumxx2p;
1941                                                                         sumyy2=sumyy2p;
1942                                                                         sumzz2=sumzz2p;
1943                                                                         sumxy2=sumxy2p;
1944                                                                         sumyz2=sumyz2p;
1945                                                                         sumzx2=sumzx2p;
1946
1947                                                                         centi=centip;
1948                                                                         inerciaii= inerciaiip;
1949
1950                                                                         centi2=centi2p;
1951                                                                         inerciaii2p= inerciaii2p;
1952
1953                                                                         centx=centxp;
1954                                                                         centy=centyp;
1955                                                                         centz=centzp;
1956
1957                                                                         centx2=centx2p;
1958                                                                         centy2=centy2p;
1959                                                                         centz2=centz2p;
1960
1961                                                                         inerciaxx= inerciaxxp;
1962                                                                         inerciayy= inerciayyp;
1963                                                                         inerciazz= inerciazzp;
1964
1965                                                                         inerciaxx2= inerciaxx2p;
1966                                                                         inerciayy2= inerciayy2p;
1967                                                                         inerciazz2= inerciazz2p;
1968
1969                                                                         inerciaxy= inerciaxyp;
1970                                                                         inerciayz= inerciayzp;
1971                                                                         inerciazx= inerciazxp;
1972
1973                                                                         inerciaxy2= inerciaxy2p;
1974                                                                         inerciayz2= inerciayz2p;
1975                                                                         inerciazx2= inerciazx2p;
1976
1977                                                                 //      sumitotal=sumitotalp;
1978
1979                                                                         costo=costop;
1980
1981                                                                 }
1982
1983                                                         }
1984                                                 }
1985                                         }
1986                                 }
1987                         }
1988                 }
1989
1990         //      printf("changes: %d\n", changes);
1991
1992
1993         }
1994
1995         //printf("total: %d\n", total);
1996
1997
1998 }
1999
2000
2001
2002
2003
2004 //----------------------------------------------------------------------------
2005 void axisExtractor02::costominimo(vtkImageData *data,  vtkImageData *data2 )
2006 {
2007         int ext[6];
2008         int i, j, k, radio;
2009         int centro[3];
2010
2011
2012         data->GetExtent(ext);
2013
2014         double *ptr;
2015         unsigned char *ptr2;
2016
2017         int ind;
2018
2019         radio=(ext[1]-ext[0])/2;
2020
2021         centro[0]=(ext[1]+ext[0])/2;
2022         centro[1]=(ext[3]+ext[2])/2;
2023         centro[2]=(ext[5]+ext[4])/2;
2024
2025
2026         for(i=0;i<=10;i++){
2027                 minis[i]=0;
2028         }
2029
2030         for(i=ext[0];i<=ext[1];i++){
2031                 for(j=ext[2];j<=ext[3];j++){
2032                         for(k=ext[4];k<=ext[5];k++){
2033                                 ptr=(double *)data->GetScalarPointer(i,j,k);
2034                                 ptr2=(unsigned char *)data2->GetScalarPointer(i,j,k);
2035
2036                                 ind=(int)*ptr2;
2037
2038                                 if(*ptr2!=0){
2039                                         if(minis[ind-1]<=*ptr && *ptr!=0 && ind<11){
2040                                                 minis[ind-1]=*ptr;
2041                                                 candit[ind-1][0]=i;
2042                                                 candit[ind-1][1]=j;
2043                                                 candit[ind-1][2]=k;
2044                                         }
2045                                 }
2046                         }
2047                 }
2048         }
2049 }
2050
2051
2052
2053
2054
2055 //----------------------------------------------------------------------------
2056 void axisExtractor02::costominimo2(vtkImageData *data,  vtkImageData *data3, int p1[3], int p2[3], int p3[3])
2057 {
2058         int ext[6];
2059         int i, j, k, radio;
2060         int punto[3];
2061         int puntob[3];
2062         int punto2[3];
2063         int puntob2[3];
2064         int puntomin[3];
2065         int puntobmin[3];
2066
2067         double *ptr;
2068         double *ptr2;
2069         unsigned char *ptr3;
2070         unsigned char *ptr4;
2071         double *ptr5;
2072         unsigned char *ptr6;
2073         double *ptr7;
2074         double *ptr8;
2075
2076         vtkImageData *data4;
2077
2078         data4=vtkImageData::New();
2079         data4->SetScalarTypeToUnsignedChar();
2080         data4->SetExtent(data->GetExtent());
2081         data4->SetSpacing(data->GetSpacing());
2082
2083         vtkImageData *data6;
2084
2085         data6=vtkImageData::New();
2086         data6->SetScalarTypeToUnsignedChar();
2087         data6->SetExtent(data->GetExtent());
2088         data6->SetSpacing(data->GetSpacing());
2089
2090
2091         vtkImageData *data2;
2092
2093         data2=vtkImageData::New();
2094         data2->SetScalarTypeToDouble();
2095         data2->SetExtent(data->GetExtent());
2096         data2->SetSpacing(data->GetSpacing());
2097
2098         vtkImageData *data8;
2099
2100         data8=vtkImageData::New();
2101         data8->SetScalarTypeToDouble();
2102         data8->SetExtent(data->GetExtent());
2103         data8->SetSpacing(data->GetSpacing());
2104
2105
2106         double mincst, mincstb;
2107
2108         bool flag=true;
2109
2110         double costoactual;
2111
2112         data->GetExtent(ext);
2113
2114         radio=(ext[1]-ext[0])/2;
2115
2116         for(i=ext[0];i<=ext[1];i++){
2117                 for(j=ext[2];j<=ext[3];j++){
2118                         for(k=ext[4];k<=ext[5];k++){
2119                                 ptr2=(double *)data2->GetScalarPointer(i,j,k);
2120                                 ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
2121                                 ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
2122                                 ptr8=(double *)data8->GetScalarPointer(i,j,k);
2123
2124                                 *ptr2 = 1000000000;
2125                                 *ptr4 = 0;
2126                                 *ptr6 = 0;
2127                                 *ptr8 = 1000000000;
2128                         }
2129                 }
2130         }
2131
2132
2133
2134         punto[0]=p1[0];
2135         punto[1]=p1[1];
2136         punto[2]=p1[2];
2137         puntob[0]=p2[0];
2138         puntob[1]=p2[1];
2139         puntob[2]=p2[2];
2140
2141
2142         ptr2=(double *)data2->GetScalarPointer(p1[0],p1[1],p1[2]);
2143         *ptr2=0;
2144
2145         ptr4=(unsigned char *)data4->GetScalarPointer(p1[0],p1[1],p1[2]);
2146         *ptr4=1;
2147
2148         ptr8=(double *)data8->GetScalarPointer(p2[0],p2[1],p2[2]);
2149         *ptr8=0;
2150
2151         ptr6=(unsigned char *)data6->GetScalarPointer(p2[0],p2[1],p2[2]);
2152         *ptr6=1;
2153
2154
2155
2156
2157
2158
2159         while(flag){
2160
2161                 ptr5=(double *)data2->GetScalarPointer(punto[0],punto[1],punto[2]);
2162                 ptr7=(double *)data8->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
2163
2164                         for(i=-1;i<=1;i++){
2165                                 for(j=-1;j<=1;j++){
2166                                         for(k=-1;k<=1;k++){
2167                                                 if(i!=0 || j!=0 || k!=0){
2168                                                         punto2[0]=punto[0]+i;
2169                                                         punto2[1]=punto[1]+j;
2170                                                         punto2[2]=punto[2]+k;
2171
2172                                                         puntob2[0]=puntob[0]+i;
2173                                                         puntob2[1]=puntob[1]+j;
2174                                                         puntob2[2]=puntob[2]+k;
2175
2176                                                         if(punto2[0]>=ext[0] && punto2[0]<=ext[1] && punto2[1]>=ext[2] && punto2[1]<=ext[3] && punto2[2]>=ext[4] &&  punto2[2]<=ext[5] ){
2177
2178                                                                 ptr=(double *)data->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2179                                                                 ptr2=(double *)data2->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2180                                                                 ptr3=(unsigned char *)data3->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2181                                                                 ptr4=(unsigned char *)data4->GetScalarPointer(punto2[0],punto2[1],punto2[2]);
2182
2183
2184                                                                 if(*ptr3!=0 && *ptr4!=2){
2185                                                                         costoactual=*ptr5+(1/(*ptr));
2186                                                                         if(*ptr2>costoactual){
2187                                                                                 *ptr2=costoactual;
2188                                                                                 *ptr4=1;
2189                                                                         }
2190                                                                 }
2191                                                         }
2192
2193                                                         if(puntob2[0]>=ext[0] && puntob2[0]<=ext[1] && puntob2[1]>=ext[2] && puntob2[1]<=ext[3] && puntob2[2]>=ext[4] &&  puntob2[2]<=ext[5] ){
2194
2195                                                                 ptr=(double *)data->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2196                                                                 ptr8=(double *)data8->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2197                                                                 ptr3=(unsigned char *)data3->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2198                                                                 ptr6=(unsigned char *)data6->GetScalarPointer(puntob2[0],puntob2[1],puntob2[2]);
2199
2200                                                                 if(*ptr3!=0 && *ptr6!=2){
2201                                                                         costoactual=*ptr7+(1/(*ptr));
2202                                                                         if(*ptr8>costoactual){
2203                                                                                 *ptr8=costoactual;
2204                                                                                 *ptr6=1;
2205                                                                         }
2206                                                                 }
2207                                                         }
2208                                                 }
2209                                         }
2210                                 }
2211                         }
2212
2213                         ptr4=(unsigned char *)data4->GetScalarPointer(punto[0],punto[1],punto[2]);
2214                         *ptr4=2;
2215
2216                         ptr6=(unsigned char *)data6->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
2217                         *ptr6=2;
2218
2219                         mincst=1000000000;
2220                         mincstb=1000000000;
2221
2222                         for(i=ext[0];i<=ext[1];i++){
2223                                 for(j=ext[2];j<=ext[3];j++){
2224                                         for(k=ext[4];k<=ext[5];k++){
2225                                                 ptr2=(double *)data2->GetScalarPointer(i,j,k);
2226                                                 ptr8=(double *)data8->GetScalarPointer(i,j,k);
2227                                                 ptr4=(unsigned char *)data4->GetScalarPointer(i,j,k);
2228                                                 ptr6=(unsigned char *)data6->GetScalarPointer(i,j,k);
2229
2230                                                 if(*ptr4==1){
2231                                                         if(mincst>=*ptr2){
2232                                                                 mincst=*ptr2;
2233                                                                 puntomin[0]=i;
2234                                                                 puntomin[1]=j;
2235                                                                 puntomin[2]=k;
2236                                                         }
2237                                                 }
2238
2239                                                 if(*ptr6==1){
2240                                                         if(mincstb>=*ptr8){
2241                                                                 mincstb=*ptr8;
2242                                                                 puntobmin[0]=i;
2243                                                                 puntobmin[1]=j;
2244                                                                 puntobmin[2]=k;
2245                                                         }
2246                                                 }
2247                                         }
2248                                 }
2249                         }
2250
2251                         punto[0]=puntomin[0];
2252                         punto[1]=puntomin[1];
2253                         punto[2]=puntomin[2];
2254
2255                         puntob[0]=puntobmin[0];
2256                         puntob[1]=puntobmin[1];
2257                         puntob[2]=puntobmin[2];
2258
2259                         ptr4=(unsigned char *)data4->GetScalarPointer(puntob[0],puntob[1],puntob[2]);
2260                         ptr6=(unsigned char *)data6->GetScalarPointer(punto[0],punto[1],punto[2]);
2261
2262
2263                         if(*ptr4==2){
2264                                 p3[0]=puntob[0];
2265                                 p3[1]=puntob[1];
2266                                 p3[2]=puntob[2];
2267                                 flag=false;
2268                         }
2269
2270                         if(*ptr6==2){
2271                                 p3[0]=punto[0];
2272                                 p3[1]=punto[1];
2273                                 p3[2]=punto[2];
2274                                 flag=false;
2275                         }
2276                 }
2277
2278                 data2->Delete();
2279                 data4->Delete();
2280                 data6->Delete();
2281                 data8->Delete();
2282 }
2283
2284
2285
2286
2287
2288 //----------------------------------------------------------------------------
2289 void axisExtractor02::invertir(vtkImageData *data )
2290 {
2291         int ext[6];
2292         int i, j, k, radio;
2293         int i2, j2, k2;
2294
2295         data->GetExtent(ext);
2296
2297         unsigned char *ptr;
2298
2299         radio=(ext[1]-ext[0])/2;
2300
2301         double tmpsqrt;
2302
2303         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2304                 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2305                         for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2306                                 ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2307                                 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2308                                 if( radio<sqrt(tmpsqrt) ){
2309                                         if(*ptr==0){
2310                                                 *ptr=255;
2311                                         }
2312                                         else{
2313                                                 *ptr=0;
2314                                         }
2315                                 }
2316                         }
2317                 }
2318         }
2319 }
2320
2321
2322
2323 //----------------------------------------------------------------------------
2324 void axisExtractor02::redondear(vtkImageData *data )
2325 {
2326         int ext[6];
2327         int i, j, k, radio;
2328         int i2, j2, k2;
2329
2330         data->GetExtent(ext);
2331
2332         unsigned char *ptr;
2333
2334         radio=(ext[1]-ext[0])/2;
2335
2336         double tmpsqrt;
2337
2338         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2339                 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2340                         for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2341                                 ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2342                                 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2343                                 if( radio<sqrt(tmpsqrt) ){
2344                                         *ptr=0;
2345                                 }
2346                                 /*if( i2==0 && j2==0 && k2==0 ){
2347                                         *ptr=255;
2348                                 }*/
2349                         }
2350                 }
2351         }
2352 }
2353
2354
2355
2356
2357
2358 //----------------------------------------------------------------------------
2359 void axisExtractor02::redondear2(vtkImageData *data )
2360 {
2361         int ext[6];
2362         int i, j, k, radio;
2363         int i2, j2, k2;
2364
2365         data->GetExtent(ext);
2366
2367         double *ptr;
2368
2369         radio=(ext[1]-ext[0])/2;
2370
2371         double tmpsqrt;
2372
2373         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2374                 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2375                         for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2376                                 ptr=(double *)  data->GetScalarPointer(i,j,k);
2377                                 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2378                                 if( radio<sqrt(tmpsqrt) ){
2379                                         *ptr=0;
2380
2381                                 }
2382
2383                         }
2384                 }
2385         }
2386 }
2387
2388
2389
2390
2391 //----------------------------------------------------------------------------
2392 void axisExtractor02::redondear3(vtkImageData *data )
2393 {
2394         int ext[6];
2395         int i, j, k, radio;
2396         int i2, j2, k2;
2397
2398         data->GetExtent(ext);
2399
2400         unsigned char *ptr;
2401
2402         radio=(ext[1]-ext[0])/2;
2403
2404         double tmpsqrt;
2405
2406         for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2407                 for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2408                         for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2409                                 ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2410                                 tmpsqrt=(i2*i2)+(j2*j2)+(k2*k2);
2411                                 if( radio<sqrt(tmpsqrt) ){
2412                                         *ptr=255;
2413
2414                                 }
2415
2416                         }
2417                 }
2418         }
2419 }
2420
2421
2422
2423
2424 //----------------------------------------------------------------------------
2425 double axisExtractor02::distancia(double a[3], double b[3] )
2426 {
2427         double tmpsqrt=((a[0]-b[0])*(a[0]-b[0]))+((a[1]-b[1])*(a[1]-b[1]))+((a[2]-b[2])*(a[2]-b[2]));
2428         return sqrt(tmpsqrt);
2429 }
2430
2431
2432
2433
2434
2435 //----------------------------------------------------------------------------
2436 double axisExtractor02::distancia(int a[3], int b[3] )
2437 {
2438         double tmpsqrt=((a[0]-b[0])*(a[0]-b[0]))+((a[1]-b[1])*(a[1]-b[1]))+((a[2]-b[2])*(a[2]-b[2]));
2439         return sqrt(tmpsqrt);
2440 }
2441
2442
2443
2444
2445
2446
2447
2448 //----------------------------------------------------------------------------
2449 void axisExtractor02::blanquear(vtkImageData *data )
2450 {
2451         int ext[6];
2452         int i, j, k;
2453
2454         data->GetExtent(ext);
2455
2456         unsigned char *ptr;
2457
2458         for(i=ext[0];i<=ext[1];i++){
2459                 for(j=ext[2];j<=ext[3];j++){
2460                         for(k=ext[4];k<=ext[5];k++){
2461                                 ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2462                                 *ptr=0;
2463                         }
2464                 }
2465         }
2466 }
2467
2468
2469
2470 //----------------------------------------------------------------------------
2471 void axisExtractor02::blanquear3(vtkImageData *data )
2472 {
2473         int ext[6];
2474         int i, j, k;
2475
2476         data->GetExtent(ext);
2477
2478         double *ptr;
2479
2480         for(i=ext[0];i<=ext[1];i++){
2481                 for(j=ext[2];j<=ext[3];j++){
2482                         for(k=ext[4];k<=ext[5];k++){
2483                                 ptr=(double *)  data->GetScalarPointer(i,j,k);
2484                                 *ptr=0;
2485                         }
2486                 }
2487         }
2488 }
2489
2490
2491
2492
2493 //----------------------------------------------------------------------------
2494 void axisExtractor02::blanquear2(vtkImageData *data )
2495 {
2496         int ext[6];
2497         int centro[3];
2498         int i, j, k;
2499
2500         data->GetExtent(ext);
2501
2502         centro[0]=(ext[1]+ext[0])/2;
2503         centro[1]=(ext[3]+ext[2])/2;
2504         centro[2]=(ext[5]+ext[4])/2;
2505
2506         unsigned char *ptr;
2507
2508         for(i=ext[0];i<=ext[1];i++){
2509                 for(j=ext[2];j<=ext[3];j++){
2510                         for(k=ext[4];k<=ext[5];k++){
2511                                 ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2512                                 if(i==centro[0] && j==centro[1] && k==centro[2]){
2513                                         *ptr=255;
2514                                 }
2515                                 else{
2516                                         *ptr=0;
2517                                 }
2518                         }
2519                 }
2520         }
2521 }
2522
2523
2524
2525
2526 //----------------------------------------------------------------------------
2527 void axisExtractor02::cilindro(vtkImageData *data, double vector[3] )
2528 {
2529         int ext[6];
2530         int i, j, k, radio;
2531         int i2, j2, k2;
2532         int centro[3];
2533         int punto[3];
2534         double centrof[3];
2535         double puntof[3];
2536         double radiof;
2537         double puntof2[3];
2538
2539
2540         data->GetExtent(ext);
2541
2542         unsigned char *ptr;
2543
2544         radio=(ext[1]-ext[0])/2;
2545         radiof=espprin[0]*(((double)ext[1]-(double)ext[0])/4);
2546         centro[0]=(ext[1]+ext[0])/2;
2547         centro[1]=(ext[3]+ext[2])/2;
2548         centro[2]=(ext[5]+ext[4])/2;
2549
2550         indextoreal(centro, centrof );
2551
2552         if(vector[0]==0 && vector[1]==0 && vector[2]==0){
2553                 blanquear2(data);
2554         }
2555         else{
2556                 puntof2[0]=centrof[0]+vector[0];
2557                 puntof2[1]=centrof[1]+vector[1];
2558                 puntof2[2]=centrof[2]+vector[2];
2559
2560                 double tmpsqrt;
2561                 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2562                         for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2563                                 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2564                                         ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2565                                         punto[0]=i;
2566                                         punto[1]=j;
2567                                         punto[2]=k;
2568                                         indextoreal(punto, puntof );
2569                                         tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
2570                                         if( radio<sqrt(tmpsqrt) ){
2571                                                 *ptr=0;
2572
2573                                         }
2574                                         else if(radiof<distanciaejepunto(puntof, centrof, puntof2)){
2575                                                 *ptr=0;
2576                                         }
2577                                         else{
2578                                                 *ptr=255;
2579                                         }
2580                                 }
2581                         }
2582                 }
2583         }
2584 }
2585
2586
2587
2588
2589
2590 //----------------------------------------------------------------------------
2591 void axisExtractor02::modelo(vtkImageData *data, unsigned char cantidad,  unsigned long vector[50][4], int candit[10][3], double radioactual, double minis[10])
2592 {
2593         int ext[6];
2594         int i, j, k, radio, radio2;
2595         int i2, j2, k2;
2596         int t;
2597         int centro[3];
2598         int punto[3];
2599         int punto2[3];
2600         double centrof[3];
2601         double puntof[3];
2602         double radiof;
2603         double puntof2[3];
2604
2605         double dist, prop, rest;
2606
2607
2608         data->GetExtent(ext);
2609
2610         unsigned char *ptr;
2611
2612         radio=(ext[1]-ext[0])/2;
2613         radio2=(int)radioactual;
2614
2615         centro[0]=(ext[1]+ext[0])/2;
2616         centro[1]=(ext[3]+ext[2])/2;
2617         centro[2]=(ext[5]+ext[4])/2;
2618
2619         indextoreal(centro, centrof );
2620
2621         double tmpsqrt;
2622
2623         for(i=ext[0];i<=ext[1];i++){
2624                 for(j=ext[2];j<=ext[3];j++){
2625                         for(k=ext[4];k<=ext[5];k++){
2626                                 ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2627                                 punto[0]=i;
2628                                 punto[1]=j;
2629                                 punto[2]=k;
2630
2631                                 if( radioactual<distancia(punto, centro) ){
2632                                         *ptr=0;
2633                                 }
2634                                 else{
2635                                         *ptr=255;
2636
2637                                 }
2638                         }
2639                 }
2640         }
2641
2642         for(t=0;t<cantidad;t++){
2643                 radiof=sqrt(minis[t]);
2644
2645                 punto2[0]=candit[t][0];
2646                 punto2[1]=candit[t][1];
2647                 punto2[2]=candit[t][2];
2648
2649                 indextoreal(punto2, puntof2 );
2650
2651
2652                 for(i=ext[0], i2=-radio;i<=ext[1];i++, i2++){
2653                         for(j=ext[2], j2=-radio;j<=ext[3];j++, j2++){
2654                                 for(k=ext[4], k2=-radio;k<=ext[5];k++, k2++){
2655                                         ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2656                                         punto[0]=i;
2657                                         punto[1]=j;
2658                                         punto[2]=k;
2659                                         indextoreal(punto, puntof );
2660
2661                                         dist=distanciaejepunto(puntof, centrof, puntof2);
2662                                         prop=proporcioejepunto(puntof, centrof, puntof2);
2663
2664                                         if(prop>=0 && prop<=1){
2665                                                 rest=(radiof-radioactual)*prop+radioactual;
2666                                                 rest=rest+0.25;
2667                                                 rest=rest*espprin[0];
2668                                         }
2669                                         else{
2670                                                 rest=dist-1;
2671                                         }
2672                                         tmpsqrt = (i2*i2)+(j2*j2)+(k2*k2);
2673                                         if( radio>sqrt(tmpsqrt) && rest>dist){
2674                                                 *ptr=255;
2675                                         }
2676                                 }
2677                         }
2678                 }
2679         }
2680 }
2681
2682
2683
2684
2685 //----------------------------------------------------------------------------
2686 void axisExtractor02::comparacion(vtkImageData *data, vtkImageData *data2, unsigned long minis[4])
2687 {
2688         int ext[6];
2689         int i, j, k;
2690         double radio;
2691
2692         data->GetExtent(ext);
2693
2694         unsigned char *ptr;
2695         unsigned char *ptr2;
2696
2697
2698         minis[0]=0;
2699         minis[1]=0;
2700         minis[2]=0;
2701         minis[3]=0;
2702
2703         int centro[3];
2704         int punto[3];
2705
2706         radio=((double)ext[1]-(double)ext[0])/2;
2707
2708         centro[0]=(ext[1]+ext[0])/2;
2709         centro[1]=(ext[3]+ext[2])/2;
2710         centro[2]=(ext[5]+ext[4])/2;
2711
2712         for(i=ext[0];i<=ext[1];i++){
2713                 for(j=ext[2];j<=ext[3];j++){
2714                         for(k=ext[4];k<=ext[5];k++){
2715                                 ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2716                                 ptr2=(unsigned char *)  data2->GetScalarPointer(i,j,k);
2717
2718                                 punto[0]=i;
2719                                 punto[1]=j;
2720                                 punto[2]=k;
2721
2722                                 if( radio>distancia(punto, centro) ){
2723
2724                                         if( *ptr==0 && *ptr2==0 ){
2725                                                 minis[0]++;
2726                                         }
2727                                         else if( *ptr!=0 && *ptr2!=0 ){
2728                                                 minis[1]++;
2729                                         }
2730                                         else if(*ptr==0 && *ptr2!=0){
2731                                                 minis[2]++;
2732                                         }
2733                                         else{
2734                                                 minis[3]++;
2735                                         }
2736                                 }
2737                         }
2738                 }
2739         }
2740 }
2741
2742
2743
2744 //----------------------------------------------------------------------------
2745 void axisExtractor02::copiar(vtkImageData *data, vtkImageData *data2 )
2746 {
2747         int ext[6];
2748         int i, j, k;
2749
2750         data->GetExtent(ext);
2751
2752         unsigned char *ptr, *ptr2;
2753
2754         for(i=ext[0];i<=ext[1];i++){
2755                 for(j=ext[2];j<=ext[3];j++){
2756                         for(k=ext[4];k<=ext[5];k++){
2757                                 ptr=(unsigned char *)   data->GetScalarPointer(i,j,k);
2758                                 ptr2=(unsigned char *)  data2->GetScalarPointer(i,j,k);
2759                                 if(*ptr!=0 && *ptr2==0){
2760                                         *ptr2=*ptr;
2761                                 }
2762                         }
2763                 }
2764         }
2765
2766         data2->Modified();
2767 }
2768
2769
2770
2771
2772 //----------------------------------------------------------------------------
2773 void axisExtractor02::copiar2(vtkImageData *data, vtkImageData *data2 )
2774 {
2775         int ext[6];
2776         int i, j, k;
2777
2778         data->GetExtent(ext);
2779
2780         double *ptr, *ptr2;
2781
2782         for(i=ext[0];i<=ext[1];i++){
2783                 for(j=ext[2];j<=ext[3];j++){
2784                         for(k=ext[4];k<=ext[5];k++){
2785                                 ptr=(double *)  data->GetScalarPointer(i,j,k);
2786                                 ptr2=(double *) data2->GetScalarPointer(i,j,k);
2787                                 if(*ptr!=0 && *ptr2==0){
2788                                         *ptr2=*ptr;
2789                                 }
2790                         }
2791                 }
2792         }
2793
2794         data2->Modified();
2795 }
2796
2797
2798
2799
2800
2801 //----------------------------------------------------------------------------
2802 double axisExtractor02::angulo(double a[3], double b[3] )
2803 {
2804         double m1=sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]));
2805         double m2=sqrt((b[0]*b[0])+(b[1]*b[1])+(b[2]*b[2]));
2806         double res=(acos(((a[0]*b[0])+(a[1]*b[1])+(a[2]*b[2]))/(m1*m2))*180)/3.1415;
2807
2808         if(res<0){
2809                 res=-res;
2810         }
2811
2812         if(res>=0 && res<90){
2813                 res=res;
2814         }
2815         else if(res>=90 && res<180){
2816                 res=180-res;
2817         }
2818         else if(res>=180 && res<270){
2819                 res=res-180;
2820         }
2821         else{
2822                 res=360-res;
2823         }
2824
2825         return res;
2826 }
2827
2828
2829
2830
2831 //----------------------------------------------------------------------------
2832 double axisExtractor02::angulo(double i1, double j1, double k1, double i2, double j2, double k2 )
2833 {
2834         double m1=sqrt((i1*i1)+(j1*j1)+(k1*k1));
2835         double m2=sqrt((i2*i2)+(j2*j2)+(k2*k2));
2836         double res=(acos(((i1*i2)+(j1*j2)+(k1*k2))/(m1*m2))*180)/3.1415;
2837
2838         if(res<0){
2839                 res=-res;
2840         }
2841
2842         if(res>=0 && res<90){
2843                 res=res;
2844         }
2845         else if(res>=90 && res<180){
2846                 res=180-res;
2847         }
2848         else if(res>=180 && res<270){
2849                 res=res-180;
2850         }
2851         else{
2852                 res=360-res;
2853         }
2854
2855         return res;
2856 }
2857
2858
2859
2860
2861
2862 //----------------------------------------------------------------------------
2863 int axisExtractor02::envolumen(int a[3], vtkImageData *datae )
2864 {
2865
2866         int res;
2867
2868         unsigned short *ptr;
2869
2870         ptr=(unsigned short *)  datae->GetScalarPointer(a);
2871
2872         if(*ptr==0){
2873                 res=0;
2874         }
2875         else{
2876                 res=1;
2877         }
2878
2879         return res;
2880 }
2881
2882
2883
2884
2885 //----------------------------------------------------------------------------
2886 int axisExtractor02::mincandit(int candit[10][3], int cantidad, double puntoanterior[3])
2887 {
2888         double distentpu;
2889         double distmientrp;
2890         int     inminenpu=0;
2891
2892         double canditreal[3];
2893         int i;
2894
2895
2896         distmientrp=64000;
2897         inminenpu=0;
2898
2899         for(i=0;i<cantidad && i<10;i++){
2900                 indextoreal(candit[i], canditreal);
2901                 distentpu=distancia(canditreal, puntoanterior);
2902                 if(distentpu<distmientrp){
2903                         distmientrp=distentpu;
2904                         inminenpu=i;
2905                 }
2906         }
2907
2908         return inminenpu;
2909
2910 }
2911
2912
2913
2914 //----------------------------------------------------------------------------
2915 int axisExtractor02::maxareacandit(unsigned long vector[50][4], int cantidad)
2916 {
2917
2918         unsigned long max;
2919         int i;
2920         int indmax=0;
2921
2922
2923         max=0;
2924
2925         for(i=0;i<cantidad && i<10;i++){
2926                 if(max<vector[i][0]){
2927                         max=vector[i][0];
2928                         indmax=i;
2929                 }
2930         }
2931
2932         return indmax;
2933
2934 }
2935
2936
2937
2938
2939 //----------------------------------------------------------------------------
2940 unsigned long axisExtractor02::totalarea(unsigned long vector[50][4], unsigned long vectorb[50][4], int cantidad, int cantidadb)
2941 {
2942
2943         unsigned long area;
2944         int i;
2945
2946         area=0;
2947
2948         for(i=0;i<cantidad && i<10;i++){
2949                 area+=vector[i][0];
2950         }
2951
2952         for(i=0;i<cantidadb;i++){
2953                 area+=vectorb[i][0];
2954         }
2955
2956         return area;
2957 }
2958
2959
2960
2961
2962
2963 //----------------------------------------------------------------------------
2964 unsigned long axisExtractor02::conecarea(unsigned long vector[50][4], int cantidad)
2965 {
2966
2967         unsigned long area;
2968         int i;
2969
2970         area=0;
2971
2972         for(i=0;i<cantidad && i<10;i++){
2973                 area+=vector[i][0];
2974         }
2975
2976         return area;
2977 }
2978
2979
2980
2981
2982
2983
2984 //----------------------------------------------------------------------------
2985 int axisExtractor02::bruled(int candit[10][3], int cantidad, vtkImageData *data4)
2986 {
2987         int i;
2988         int burned;
2989
2990         burned=0;
2991         for(i=0;i<cantidad && i<10;i++){
2992                 if(envolumen(candit[i], data4)!=0){
2993                         burned++;
2994                 }
2995         }
2996
2997         return burned;
2998 }
2999
3000
3001
3002
3003 //----------------------------------------------------------------------------
3004 double axisExtractor02::correction(int candit[10][3], int cantidad, vtkImageData *data, int indicecorregido[3], double puntocorregido[3], int indiceanterior[3], double radioanterior, int indicepre[3], double radiopre)
3005 {
3006
3007
3008         int rap;
3009         int punto[3];
3010         double puntof[3];
3011         double puntof2[3];
3012         double centrof[3];
3013         double radiof;
3014         double maxcent;
3015         double mincent;
3016         double distcorr;
3017         double distcorr2;
3018         double distcorr3;
3019         double dist;
3020         double prop;
3021         double rest;
3022         double *ptr;
3023         int ext[6];
3024         int i, j, k, radio;
3025         int centro[3];
3026
3027         data->GetExtent(ext);
3028
3029         radio=(ext[1]-ext[0])/2;
3030         radiof=((double)ext[1]-(double)ext[0])/2;
3031         centro[0]=(ext[1]+ext[0])/2;
3032         centro[1]=(ext[3]+ext[2])/2;
3033         centro[2]=(ext[5]+ext[4])/2;
3034
3035         indicecorregido[0]=centro[0];
3036         indicecorregido[1]=centro[1];
3037         indicecorregido[2]=centro[2];
3038         indextoreal(indicecorregido, puntocorregido);
3039         indextoreal(indiceanterior, centrof);
3040         indextoreal(indicepre, puntof2);
3041
3042         rap= (int)( (double)radio/2 );
3043         maxcent=0;
3044         mincent=64000;
3045
3046         for(i=ext[0];i<=ext[1];i++){
3047                 for(j=ext[2];j<=ext[3];j++){
3048                         for(k=ext[4];k<=ext[5];k++){
3049                                 ptr=(double *)  data->GetScalarPointer(i,j,k);
3050                                 punto[0]=i;
3051                                 punto[1]=j;
3052                                 punto[2]=k;
3053                                 indextoreal(punto, puntof );
3054
3055                                 distcorr3=distancia(indiceanterior, indicepre);
3056                                 distcorr2=distancia(punto, indicepre);
3057                                 distcorr=distancia(punto, indiceanterior);
3058
3059                                 dist=distanciaejepunto(puntof, centrof, puntof2);
3060                                 prop=proporcioejepunto(puntof, centrof, puntof2);
3061
3062                                 if(prop>=0 && prop<=1){
3063                                         rest=(radiopre-radioanterior)*prop+radioanterior;
3064                                         rest=rest*espprin[0];
3065                                 }
3066                                 else{
3067                                         rest=dist-1;
3068                                 }
3069
3070                                 if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && rest>dist){
3071                                 //if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && rest>dist){
3072                                 //if( radioanterior<distcorr && distcorr3>=distcorr && distcorr2<distcorr3 && distcorr+radio>distcorr3+1 && distcorr-radio<distcorr3-1 && distcorr2<=distcorr3-radioanterior){
3073                                         if( maxcent<*ptr){
3074                                                         mincent=distcorr;
3075                                                         maxcent=*ptr;
3076                                                         indicecorregido[0]=i;
3077                                                         indicecorregido[1]=j;
3078                                                         indicecorregido[2]=k;
3079                                                         indextoreal(indicecorregido, puntocorregido);
3080                                         }
3081                                         else if( maxcent==*ptr){
3082                                                 if( mincent>distcorr){
3083                                                         mincent=distcorr;
3084                                                         maxcent=*ptr;
3085                                                         indicecorregido[0]=i;
3086                                                         indicecorregido[1]=j;
3087                                                         indicecorregido[2]=k;
3088                                                         indextoreal(indicecorregido, puntocorregido);
3089                                                 }
3090                                         }
3091                                 }
3092                         }
3093                 }
3094         }
3095
3096         ptr=(double *)  data->GetScalarPointer(centro[0],centro[1],centro[2]);
3097
3098         return sqrt(*ptr);
3099 }
3100
3101
3102
3103
3104 //----------------------------------------------------------------------------
3105 double axisExtractor02::correction2(int candit[10][3], int cantidad, vtkImageData *data, int indicecorregido[3], double puntocorregido[3], int indiceanterior[3], double radioanterior)
3106 {
3107
3108
3109         int rap;
3110         int punto[3];
3111         double maxcent;
3112         double mincent;
3113         double distcorr;
3114         double distcorr2;
3115         double *ptr;
3116         int ext[6];
3117         int i, j, k, radio;
3118         int centro[3];
3119
3120         data->GetExtent(ext);
3121
3122         radio=(ext[1]-ext[0])/2;
3123         centro[0]=(ext[1]+ext[0])/2;
3124         centro[1]=(ext[3]+ext[2])/2;
3125         centro[2]=(ext[5]+ext[4])/2;
3126
3127         indicecorregido[0]=centro[0];
3128         indicecorregido[1]=centro[1];
3129         indicecorregido[2]=centro[2];
3130         indextoreal(indicecorregido, puntocorregido);
3131
3132         rap= (int)((double)radio/2);
3133         maxcent=0;
3134         mincent=64000;
3135
3136         for(i=centro[0]-rap;i<=centro[0]+rap;i++){
3137                 for(j=centro[1]-rap;j<=centro[1]+rap;j++){
3138                         for(k=centro[2]-rap;k<=centro[2]+rap;k++){
3139                                 ptr=(double *)  data->GetScalarPointer(i,j,k);
3140                                 punto[0]=i;
3141                                 punto[1]=j;
3142                                 punto[2]=k;
3143
3144                                 distcorr2=distancia(punto, centro);
3145                                 distcorr=distancia(punto, indiceanterior);
3146
3147                                 if( rap>distcorr2 ){
3148                                         if( maxcent<*ptr){
3149                                                         mincent=distcorr;
3150                                                         maxcent=*ptr;
3151                                                         indicecorregido[0]=i;
3152                                                         indicecorregido[1]=j;
3153                                                         indicecorregido[2]=k;
3154                                                         indextoreal(indicecorregido, puntocorregido);
3155                                         }
3156                                         else if( maxcent==*ptr){
3157                                                 if( mincent>distcorr){
3158                                                         mincent=distcorr;
3159                                                         maxcent=*ptr;
3160                                                         indicecorregido[0]=i;
3161                                                         indicecorregido[1]=j;
3162                                                         indicecorregido[2]=k;
3163                                                         indextoreal(indicecorregido, puntocorregido);
3164                                                 }
3165
3166                                         }
3167                                 }
3168
3169                         }
3170                 }
3171         }
3172
3173         return sqrt(maxcent);
3174 }
3175
3176
3177
3178
3179
3180 //----------------------------------------------------------------------------
3181 void axisExtractor02::avanzar()
3182 {
3183
3184         double puntoactual[3];
3185         double puntocorregido[3];
3186         double puntoanterior[3];
3187         double puntosiguiente[3];
3188         //double puntointer[3];
3189         double puntopre[3];
3190
3191         int indiceactual[3];
3192         int indiceanterior[3];
3193         int indicecorregido[3];
3194         //int indiceinter[3];
3195         int indicepre[3];
3196         double radioactual;
3197         double dirant[3];
3198
3199         unsigned char cantidad;
3200         unsigned char cantidadb;
3201
3202         unsigned long vector[50][4];
3203         unsigned long vectorb[50][4];
3204
3205         int vectortemp[3];
3206         double tempc;
3207
3208         int extint[6];
3209
3210         int extintreal[6];
3211
3212         int indexp;
3213
3214         double provvc[3];
3215
3216         int radiotrabajo;
3217
3218         int flag =0;
3219         int flag2 =0;
3220         int flag3 =0;
3221         int flag4 =0;
3222         int flag5 =0;
3223         int flag6 =0;
3224         int flag7 =0;
3225         int flag8 =0;
3226         int flag9 =0;
3227         int flag10 =0;
3228         int flag11 =0;
3229         int flag12 =0;
3230         int flag13 =0;
3231         int flag14 =0;
3232         int flag15 =0;
3233         int flag16 =0;
3234         int flag17 =0;
3235         int flag18 =0;
3236         int flag19 =0;
3237         int flag20 =0;
3238         int flag21 =0;
3239         int flag22 =0;
3240
3241         int i, j;
3242
3243         double radiotrabajof;
3244
3245         double proprad=2;
3246
3247         double radiopred;
3248         double radioanterior;
3249         double radioprinc;
3250
3251         int burned;
3252
3253         int indant;
3254         double canditreal[3];
3255         int indmaxarea;
3256         unsigned long totalsurf;
3257         unsigned long conecsurf;
3258
3259         indmaxarea=0;
3260         unsigned long caf[4];
3261
3262         double rel1;
3263         double rel2;
3264         double rel3;
3265         double rel4;
3266         double rel5;
3267         double rel6;
3268
3269         if(!m_Stack0.empty()){
3270
3271
3272         /*      indexp=m_Stack.top();
3273
3274                 puntoactual[0]=m_Stack0.top();
3275                 puntoactual[1]=m_Stack1.top();
3276                 puntoactual[2]=m_Stack2.top();
3277
3278                 puntoanterior[0]=m_Stack3.top();
3279                 puntoanterior[1]=m_Stack4.top();
3280                 puntoanterior[2]=m_Stack5.top();
3281
3282                 puntopre[0]=m_Stack6.top();
3283                 puntopre[1]=m_Stack7.top();
3284                 puntopre[2]=m_Stack8.top();
3285
3286
3287
3288                 radiopred=m_Stackr.top();
3289                 radioanterior=m_Stackra.top();
3290                 radioprinc=m_Stackrp.top();*/
3291
3292                 indexp=m_Stack.front();
3293
3294                 puntoactual[0]=m_Stack0.front();
3295                 puntoactual[1]=m_Stack1.front();
3296                 puntoactual[2]=m_Stack2.front();
3297
3298                 puntoanterior[0]=m_Stack3.front();
3299                 puntoanterior[1]=m_Stack4.front();
3300                 puntoanterior[2]=m_Stack5.front();
3301
3302                 puntopre[0]=m_Stack6.front();
3303                 puntopre[1]=m_Stack7.front();
3304                 puntopre[2]=m_Stack8.front();
3305
3306                 radiopred=m_Stackr.front();
3307                 radioanterior=m_Stackra.front();
3308                 radioprinc=m_Stackrp.front();
3309
3310                 radiotrabajo= (int) (proprad*radiopred);
3311                 radiotrabajof=proprad*radiopred;
3312
3313                 if(radiotrabajof-radiotrabajo>0.5){
3314                         radiotrabajo=radiotrabajo+1;
3315                 }
3316
3317 // EED 15 Mai 2007 .NET
3318 //              radiotrabajo=max(minant,radiotrabajo);
3319                 if (minant>radiotrabajo)
3320                 {
3321                         radiotrabajo=minant;
3322                 }
3323
3324                 /*m_Stack0.pop();
3325                 m_Stack1.pop();
3326                 m_Stack2.pop();
3327                 m_Stack3.pop();
3328                 m_Stack4.pop();
3329                 m_Stack5.pop();
3330                 m_Stack6.pop();
3331                 m_Stack7.pop();
3332                 m_Stack8.pop();
3333
3334                 m_Stack.pop();
3335                 m_Stackr.pop();
3336                 m_Stackra.pop();
3337                 m_Stackrp.pop();*/
3338
3339                 m_Stack0.pop_front();
3340                 m_Stack1.pop_front();
3341                 m_Stack2.pop_front();
3342                 m_Stack3.pop_front();
3343                 m_Stack4.pop_front();
3344                 m_Stack5.pop_front();
3345                 m_Stack6.pop_front();
3346                 m_Stack7.pop_front();
3347                 m_Stack8.pop_front();
3348
3349                 m_Stack.pop_front();
3350                 m_Stackr.pop_front();
3351                 m_Stackra.pop_front();
3352                 m_Stackrp.pop_front();
3353
3354                 dirant[0]=puntoanterior[0]-puntoactual[0];
3355                 dirant[1]=puntoanterior[1]-puntoactual[1];
3356                 dirant[2]=puntoanterior[2]-puntoactual[2];
3357
3358
3359                 realtoindex(puntoactual, indiceactual);
3360                 realtoindex(puntoanterior, indiceanterior);
3361                 realtoindex(puntopre, indicepre);
3362
3363                 radioactual=radiopred;
3364
3365                 if(radiotrabajo<=maxant){
3366
3367                         extint[0]=indiceactual[0]-radiotrabajo;
3368                         extint[1]=indiceactual[0]+radiotrabajo;
3369                         extint[2]=indiceactual[1]-radiotrabajo;
3370                         extint[3]=indiceactual[1]+radiotrabajo;
3371                         extint[4]=indiceactual[2]-radiotrabajo;
3372                         extint[5]=indiceactual[2]+radiotrabajo;
3373
3374                         extrac->SetVOI(extint);
3375                         extrac->UpdateWholeExtent();
3376                         extrac->Update();
3377                         extrac->GetOutput()->GetExtent(extintreal);
3378
3379                         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]){
3380                                 fprintf(stream, "salio\t");
3381                                 flag4 =1;
3382                         }
3383                         else{
3384
3385                                 data1->Delete();
3386                                 data1=vtkImageData::New();
3387                                 data1->SetScalarTypeToUnsignedChar();
3388                                 data1->SetExtent(extrac->GetOutput()->GetExtent());
3389                                 data1->SetSpacing(extrac->GetOutput()->GetSpacing());
3390
3391                                 cilindro(data1, dirant );
3392
3393                                 optim(extrac->GetOutput(), data1 );
3394
3395                                 connect->SetInput(data1);
3396                                 connect->RemoveAllSeeds();
3397                                 connect->AddSeed(indiceactual[0],indiceactual[1],indiceactual[2]);
3398                                 connect->UpdateWholeExtent();
3399                                 connect->Update();
3400
3401                                 data2->Delete();
3402                                 data2=vtkImageData::New();
3403                                 data2->SetScalarTypeToUnsignedChar();
3404                                 data2->SetExtent(connect->GetOutput()->GetExtent());
3405                                 data2->SetSpacing(connect->GetOutput()->GetSpacing());
3406
3407                                 blanquear(data2);
3408
3409                                 cantidad=find_components(connect->GetOutput(), data2, 0, vector);
3410
3411                                 data3->Delete();
3412                                 data3=vtkImageData::New();
3413                                 data3->SetScalarTypeToUnsignedChar();
3414                                 data3->SetExtent(connect->GetOutput()->GetExtent());
3415                                 data3->SetSpacing(connect->GetOutput()->GetSpacing());
3416
3417                                 blanquear(data3);
3418
3419                                 cantidadb=find_componentsb(connect->GetOutput(), data3, 0, vectorb);
3420
3421                                 redondear3(connect->GetOutput());
3422
3423                                 distance->Update();
3424
3425                                 redondear2(distance->GetOutput());
3426
3427                                 redondear(connect->GetOutput());
3428
3429                                 costominimo(distance->GetOutput(), data2);
3430
3431
3432                                 if(buenos==0){
3433                                         radioactual=correction2(candit, cantidad, distance->GetOutput(),  indicecorregido,  puntocorregido, indiceanterior, radioanterior);
3434                                 }
3435                                 else{
3436                                         radioactual=correction(candit, cantidad, distance->GetOutput(),  indicecorregido,  puntocorregido, indiceanterior, radioanterior, indicepre, radioprinc);
3437                                 }
3438
3439
3440                                 //inpeq(extrac->GetOutput(), indicecorregido, radioactual);
3441
3442                                 indant=mincandit(candit, cantidad, puntoanterior);
3443
3444                                 indmaxarea=maxareacandit(vector, cantidad);
3445                                 totalsurf=totalarea(vector, vectorb, cantidad, cantidadb);
3446                                 conecsurf=conecarea(vector, cantidad);
3447
3448                                 indextoreal(candit[indant], canditreal);
3449
3450                                 burned=bruled(candit, cantidad, data4);
3451
3452                                 data6->Delete();
3453                                 data6=vtkImageData::New();
3454                                 data6->SetScalarTypeToUnsignedChar();
3455                                 data6->SetExtent(extrac->GetOutput()->GetExtent());
3456                                 data6->SetSpacing(extrac->GetOutput()->GetSpacing());
3457
3458                                 modelo(data6, cantidad,  vector, candit, radioactual, minis);
3459                                 comparacion(connect->GetOutput(), data6, caf);
3460
3461                                 rel1=((double)caf[0]+(double)caf[1])/((double)caf[0]+(double)caf[1]+(double)caf[2]+(double)caf[3]);
3462                                 rel2=((double)caf[1])/((double)caf[1]+(double)caf[3]);
3463                                 rel3=((double)caf[1])/((double)caf[1]+(double)caf[2]);
3464                                 rel4=(double)conecsurf/(double)totalsurf;
3465                                 rel5=(double)vector[indmaxarea][0]/(double)totalsurf;
3466                                 rel6=(double)cant/((double)cant+(double)cant2);
3467
3468
3469                                 if(rel1<0.5 || rel2<0.5 || rel3<0.5 || (rel1<0.75 && rel2<0.75) || (rel2<0.75 && rel3<0.75)|| (rel1<0.75 && rel3<0.75)){
3470                                         flag5=1;
3471                                 }
3472
3473                                 if(indicecorregido[0]!=indiceactual[0] || indicecorregido[1]!=indiceactual[1] || indicecorregido[2]!=indiceactual[2]){
3474                                         flag7=1;
3475                                 }
3476
3477                                 if(burned>=cantidad ){
3478                                         flag15=1;
3479                                 }
3480
3481                                 if(burned==0 ){
3482                                         flag16=1;
3483                                 }
3484
3485                                 if(cantidad<2 ){
3486                                         flag17=1;
3487                                 }
3488
3489                                 if(flag7==1 && flagg>4){
3490                                         flag7=0;
3491                                 }
3492
3493                                 if(buenos==0){
3494                                         flag19=1;
3495                                 }
3496
3497                                 if(mejordst!=0){
3498                                         flag18=1;
3499                                 }
3500
3501                                 if(flag7==1 && flag18==1){
3502                                         for(i=0;i<flagg;i++){
3503                                                 if(visited[i][0]==indiceactual[0] && visited[i][1]==indiceactual[1] && visited[i][2]==indiceactual[2] && visitedrad[i]==radiopred){
3504                                                         flag2=1;
3505                                                 }
3506                                         }
3507                                 }
3508
3509                                 if((flag19==1 && flagg2<5) || (flag19==1 && flagg2<10 && ( rel4>0.25 || rel5>0.13 || rel6>0.4 || flag17==1  || flag15==1 || flag5==1) ) || (flagg2<10 && ( ((rel4>0.25 || rel5>0.25 || rel6>0.4) && cantidad<3) || flag17==1  || flag15==1 || flag5==1) )){
3510                                                 flag8=1;
3511                                 }
3512
3513                                 if(centi<centi2){
3514                                         flag3=1;
3515                                 }
3516
3517                                 if(angulo(dirant[0], dirant[1], dirant[2], canditreal[0]-puntoactual[0], canditreal[1]-puntoactual[1], canditreal[2]-puntoactual[2] )>65){
3518                                         flag11=1;
3519                                 }
3520
3521                                 if(angulo(Vp[0][ejeminp], Vp[1][ejeminp], Vp[2][ejeminp], Vi[0][ejemini], Vi[1][ejemini], Vi[2][ejemini] )>60){
3522                                         flag6=1;
3523                                 }
3524
3525                                 if((flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0)  ) && (flag7==0 && flag8==0)){
3526                                         flag9=1;
3527                                         //fprintf(stream, "malo\t");
3528                                 }
3529
3530                                 if(flag3==1 || flag5==1 || (flag16==1 && flag11==1 && flag19==0)  ){
3531                                         flag20=1;
3532                                 }
3533
3534                                 if((mejordst<radioactual || (mejordst==radioactual && mejorrad>radiopred) || (mejordst==radioactual && mejorrad<=radiopred && mejorcant<cantidad )) && flag17==0 && flag20==0 && flag4==0 && flag12==0 && flag15==0 ){
3535                                         mejordst=radioactual;
3536                                         mejor[0]=puntoactual[0];
3537                                         mejor[1]=puntoactual[1];
3538                                         mejor[2]=puntoactual[2];
3539                                         mejorrad=radiopred;
3540                                         mejorcant=cantidad;
3541                                 }
3542                         }
3543                 }
3544                 else{
3545                         fprintf(stream, "grande\t");
3546                         flag12 =1;
3547                 }
3548
3549                 if(mejordst!=0){
3550                         flag18=1;
3551                 }
3552
3553                 if(flagg>4){
3554                         flag13=1;
3555                 }
3556
3557                 if(flagg2>4){
3558                         flag14=1;
3559                 }
3560
3561                 if( flag9==1 || flag4==1 || flag12==1 || ((cantidad<2 || burned==cantidad) && (flag7==0 && flag8==0)) ){
3562                         frama=1;
3563                 }
3564                 else{
3565                         frama=0;
3566                 }
3567
3568                 if( cantidad>2 && (flag7==0 && flag8==0) ){
3569                         fseg=1;
3570                 }
3571                 else{
3572                         fseg=0;
3573                 }
3574
3575                 if((flag7==0 && flag8==0 && (mejordst>radioactual || (mejordst==radioactual && mejorrad<radiopred)) && flag19==0) || flag2==1){
3576
3577                         flagg=10;
3578                         flagg2=10;
3579                         mejordst=0;
3580
3581                         fprintf(stream, "mejor\t");
3582
3583                         m_Stack0.push_front(mejor[0]);
3584                         m_Stack1.push_front(mejor[1]);
3585                         m_Stack2.push_front(mejor[2]);
3586                         m_Stack3.push_front(puntoanterior[0]);
3587                         m_Stack4.push_front(puntoanterior[1]);
3588                         m_Stack5.push_front(puntoanterior[2]);
3589                         m_Stack6.push_front(puntopre[0]);
3590                         m_Stack7.push_front(puntopre[1]);
3591                         m_Stack8.push_front(puntopre[2]);
3592                         m_Stack.push_front(indexp);
3593                         m_Stackr.push_front(mejorrad);
3594                         m_Stackra.push_front(radioanterior);
3595                         m_Stackrp.push_front(radioprinc);
3596
3597                         mejorrad=0;
3598                         mejorcant=0;
3599
3600                 }
3601
3602                 else if((flag17==1 || flag9==1 || flag4==1 || flag12==1 || flag15==1) && flag18==1){
3603                         flagg=10;
3604                         flagg2=10;
3605                         mejordst=0;
3606
3607                         fprintf(stream, "descorrigio\t");
3608
3609                         m_Stack0.push_front(mejor[0]);
3610                         m_Stack1.push_front(mejor[1]);
3611                         m_Stack2.push_front(mejor[2]);
3612                         m_Stack3.push_front(puntoanterior[0]);
3613                         m_Stack4.push_front(puntoanterior[1]);
3614                         m_Stack5.push_front(puntoanterior[2]);
3615                         m_Stack6.push_front(puntopre[0]);
3616                         m_Stack7.push_front(puntopre[1]);
3617                         m_Stack8.push_front(puntopre[2]);
3618                         m_Stack.push_front(indexp);
3619                         m_Stackr.push_front(mejorrad);
3620                         m_Stackra.push_front(radioanterior);
3621                         m_Stackrp.push_front(radioprinc);
3622
3623                         mejorrad=0;
3624                         mejorcant=0;
3625
3626                 }
3627
3628                 else if(flag12==0 && flag4==0 && flag9==0){
3629                         if(flag7==0 && flag8==0){
3630
3631                                 if(flag17==0 && flag15==0){
3632                                         flagg=0;
3633                                         flagg2=0;
3634                                         mejordst=0;
3635                                         mejorrad=0;
3636                                         mejorcant=0;
3637
3638                                         fprintf(stream, "inserto\t");
3639
3640                                         if(flag19==0){
3641                                                 /*costominimo2(distance->GetOutput(), connect->GetOutput(), indiceactual, candit[indant], indiceinter);
3642                                                 indextoreal(indiceinter, puntointer);
3643                                                 realtoreal2(puntointer, provvc);
3644                                                 points->InsertPoint(buenos,provvc);
3645
3646                                                 lineas->InsertNextCell(2);
3647                                                 lineas->InsertCellPoint(indexp);
3648                                                 lineas->InsertCellPoint(buenos);
3649
3650                                                 buenos++;
3651
3652                                                 realtoreal2(puntoactual, provvc);
3653                                                 points->InsertPoint(buenos,provvc);
3654
3655                                                 lineas->InsertNextCell(2);
3656                                                 lineas->InsertCellPoint(buenos-1);
3657                                                 lineas->InsertCellPoint(buenos);
3658
3659                                                 buenos++;*/
3660
3661                                                 realtoreal2(puntoactual, provvc);
3662                                                 points->InsertPoint(buenos,provvc);
3663
3664                                                 lineas->InsertNextCell(2);
3665                                                 lineas->InsertCellPoint(indexp);
3666                                                 lineas->InsertCellPoint(buenos);
3667
3668                                                 buenos++;
3669
3670                                         }
3671                                         else{
3672                                                 realtoreal2(puntoactual, provvc);
3673                                                 points->InsertPoint(buenos,provvc);
3674
3675                                                 buenos++;
3676                                         }
3677
3678
3679                                 }
3680                                 else{
3681                                         flagg=0;
3682                                         flagg2=0;
3683                                         mejordst=0;
3684                                         mejorrad=0;
3685                                         mejorcant=0;
3686                                         fprintf(stream, "para\t");
3687                                 }
3688                         }
3689                         else if(flag7==0 && flag8==1){
3690                                 flagg2++;
3691
3692                                 fprintf(stream, "agrando\t");
3693
3694                                 m_Stack0.push_front(puntoactual[0]);
3695                                 m_Stack1.push_front(puntoactual[1]);
3696                                 m_Stack2.push_front(puntoactual[2]);
3697                                 m_Stack3.push_front(puntoanterior[0]);
3698                                 m_Stack4.push_front(puntoanterior[1]);
3699                                 m_Stack5.push_front(puntoanterior[2]);
3700                                 m_Stack6.push_front(puntopre[0]);
3701                                 m_Stack7.push_front(puntopre[1]);
3702                                 m_Stack8.push_front(puntopre[2]);
3703
3704                                 m_Stack.push_front(indexp);
3705                                 m_Stackr.push_front(((double)radiotrabajo+1)/2);
3706                                 m_Stackra.push_front(radioanterior);
3707                                 m_Stackrp.push_front(radioprinc);
3708
3709                         }
3710                         else if(flag7==1 && flag8==1){
3711                                 visited[flagg][0]=indiceactual[0];
3712                                 visited[flagg][1]=indiceactual[1];
3713                                 visited[flagg][2]=indiceactual[2];
3714                                 visitedrad[flagg]=radiopred;
3715                                 flagg++;
3716                                 flagg2++;
3717
3718                                 fprintf(stream, "agrando y corrigiendo\t");
3719
3720                                 m_Stack0.push_front(puntocorregido[0]);
3721                                 m_Stack1.push_front(puntocorregido[1]);
3722                                 m_Stack2.push_front(puntocorregido[2]);
3723                                 m_Stack3.push_front(puntoanterior[0]);
3724                                 m_Stack4.push_front(puntoanterior[1]);
3725                                 m_Stack5.push_front(puntoanterior[2]);
3726                                 m_Stack6.push_front(puntopre[0]);
3727                                 m_Stack7.push_front(puntopre[1]);
3728                                 m_Stack8.push_front(puntopre[2]);
3729
3730                                 m_Stack.push_front(indexp);
3731                                 m_Stackr.push_front(((double)radiotrabajo+1)/2);
3732                                 m_Stackra.push_front(radioanterior);
3733                                 m_Stackrp.push_front(radioprinc);
3734                         }
3735                         else{
3736                                 visited[flagg][0]=indiceactual[0];
3737                                 visited[flagg][1]=indiceactual[1];
3738                                 visited[flagg][2]=indiceactual[2];
3739                                 visitedrad[flagg]=radiopred;
3740                                 flagg++;
3741
3742                                 if(flag19==1){
3743                                         flagg=10;
3744                                 }
3745
3746                                 fprintf(stream, "corrigio\t");
3747
3748                                 m_Stack0.push_front(puntocorregido[0]);
3749                                 m_Stack1.push_front(puntocorregido[1]);
3750                                 m_Stack2.push_front(puntocorregido[2]);
3751                                 m_Stack3.push_front(puntoanterior[0]);
3752                                 m_Stack4.push_front(puntoanterior[1]);
3753                                 m_Stack5.push_front(puntoanterior[2]);
3754
3755                                 m_Stack6.push_front(puntopre[0]);
3756                                 m_Stack7.push_front(puntopre[1]);
3757                                 m_Stack8.push_front(puntopre[2]);
3758
3759                                 m_Stack.push_front(indexp);
3760                                 m_Stackr.push_front(radiopred);
3761                                 m_Stackra.push_front(radioanterior);
3762                                 m_Stackrp.push_front(radioprinc);
3763
3764
3765                         }
3766
3767                         if(flag7==0 && flag8==0){
3768
3769                                 if(flag17==0 && flag15==0){
3770
3771                                         for(i=0;i<cantidad && i<10;i++){
3772                                                 for(j=i;j<cantidad && i<10;j++){
3773                                                         if(minis[j]>minis[i]){
3774                                                                 vectortemp[0]=candit[i][0];
3775                                                                 vectortemp[1]=candit[i][1];
3776                                                                 vectortemp[2]=candit[i][2];
3777                                                                 tempc=minis[i];
3778
3779                                                                 candit[i][0]=candit[j][0];
3780                                                                 candit[i][1]=candit[j][1];
3781                                                                 candit[i][2]=candit[j][2];
3782                                                                 minis[i]=minis[j];
3783
3784                                                                 candit[j][0]=vectortemp[0];
3785                                                                 candit[j][1]=vectortemp[1];
3786                                                                 candit[j][2]=vectortemp[2];
3787                                                                 minis[j]=tempc;
3788                                                         }
3789                                                 }
3790                                         }
3791
3792
3793                                         if(flag16==1){
3794                                                 for(i=0;i<cantidad && i<10;i++){
3795                                                         if(i!=indant || buenos<2){
3796                                                                 indextoreal(candit[i], puntosiguiente);
3797
3798                                                                 m_Stack0.push_front(puntosiguiente[0]);
3799                                                                 m_Stack1.push_front(puntosiguiente[1]);
3800                                                                 m_Stack2.push_front(puntosiguiente[2]);
3801                                                                 m_Stack3.push_front(puntoactual[0]);
3802                                                                 m_Stack4.push_front(puntoactual[1]);
3803                                                                 m_Stack5.push_front(puntoactual[2]);
3804
3805                                                                 m_Stack6.push_front(puntosiguiente[0]);
3806                                                                 m_Stack7.push_front(puntosiguiente[1]);
3807                                                                 m_Stack8.push_front(puntosiguiente[2]);
3808
3809                                                                 double prov1a=sqrt(minis[i]);
3810
3811                                                                 m_Stack.push_front(buenos-1);
3812                                                                 m_Stackr.push_front(prov1a);
3813                                                                 m_Stackra.push_front(radioactual);
3814                                                                 m_Stackrp.push_front(prov1a);
3815
3816                                                         }
3817                                                 }
3818                                         }
3819                                         else{
3820                                                 for(i=0;i<cantidad && i<10;i++){
3821                                                         indextoreal(candit[i], puntosiguiente);
3822
3823                                                         if(envolumen(candit[i], data4)==0){
3824                                                                 m_Stack0.push_front(puntosiguiente[0]);
3825                                                                 m_Stack1.push_front(puntosiguiente[1]);
3826                                                                 m_Stack2.push_front(puntosiguiente[2]);
3827                                                                 m_Stack3.push_front(puntoactual[0]);
3828                                                                 m_Stack4.push_front(puntoactual[1]);
3829                                                                 m_Stack5.push_front(puntoactual[2]);
3830
3831                                                                 m_Stack6.push_front(puntosiguiente[0]);
3832                                                                 m_Stack7.push_front(puntosiguiente[1]);
3833                                                                 m_Stack8.push_front(puntosiguiente[2]);
3834
3835                                                                 double prov1b=sqrt(minis[i]);
3836
3837                                                                 m_Stack.push_front(buenos-1);
3838                                                                 m_Stackr.push_front(prov1b);
3839                                                                 m_Stackra.push_front(radioactual);
3840                                                                 m_Stackrp.push_front(prov1b);
3841                                                         }
3842                                                 }
3843                                         }
3844
3845                                         copiar(connect->GetOutput(),  data4 );
3846                                         //copiar2(distance->GetOutput(),  data5 );
3847                                 }
3848                         }
3849                 }
3850                 else{
3851                 //      fprintf(stream, "algo");
3852                         mejordst=0;
3853                         mejorrad=0;
3854                         mejorcant=0;
3855                         flagg=0;
3856                         flagg2=0;
3857
3858                 }
3859
3860
3861
3862
3863                 /*
3864                 printf("wi[0]: %f\n", wi[0]);
3865                 printf("wi[1]: %f\n", wi[1]);
3866                 printf("wi[2]: %f\n", wi[2]);
3867
3868                 printf("inerciari: %f\n", inerciari);
3869                 printf("inerciariy: %f\n", inerciariy);
3870                 printf("inerciariz: %f\n", inerciariz);
3871
3872                 printf("ri1: %f\n", inerciari/inerciariy);
3873                 printf("ri2: %f\n", inerciariy/inerciariz);
3874                 printf("ri3: %f\n", inerciari/inerciariz);
3875
3876                 printf("inerciarp: %f\n", inerciarp);
3877                 printf("inerciarpy: %f\n", inerciarpy);
3878                 printf("inerciarpz: %f\n", inerciarpz);
3879
3880                 printf("rp1: %f\n", inerciarp/inerciarpy);
3881                 printf("rp2: %f\n", inerciarpy/inerciarpz);
3882                 printf("rp3: %f\n", inerciarp/inerciarpz);
3883
3884                 printf("angulo: %f\n",angulo(dirant[0], dirant[1], dirant[2], canditreal[0]-puntoactual[0], canditreal[1]-puntoactual[1], canditreal[2]-puntoactual[2] ));
3885                 printf("angulo2: %f\n", angulo(Vp[0][ejeminp], Vp[1][ejeminp], Vp[2][ejeminp], Vi[0][ejemini], Vi[1][ejemini], Vi[2][ejemini]));
3886
3887                 printf("centi: %f\n", centi);
3888                 printf("centi2: %f\n", centi2);
3889                 */
3890
3891
3892
3893                 fprintf(stream, "%d\t", indexp);
3894                 fprintf(stream, "%d\t", buenos);
3895                 fprintf(stream, "%d\t", radiotrabajo);
3896                 fprintf(stream, "%f\t", radioactual);
3897                 fprintf(stream, "%f\t", radiopred);
3898                 fprintf(stream, "%f\t", radioanterior);
3899                 fprintf(stream, "%f\t", mejordst);
3900                 fprintf(stream, "%f\t", mejorrad);
3901                 fprintf(stream, "%d\t", mejorcant);
3902
3903
3904                 /*
3905                 fprintf(stream, "%d\t", cant);
3906                 fprintf(stream, "%d\t", cant2);
3907                 fprintf(stream, "%f\t", (double)cant/((double)cant+(double)cant2));
3908                 */
3909
3910                 fprintf(stream, "%d\t", flag2);
3911                 fprintf(stream, "%d\t", flag3);
3912                 fprintf(stream, "%d\t", flag4);
3913                 fprintf(stream, "%d\t", flag5);
3914                 fprintf(stream, "%d\t", flag6);
3915                 fprintf(stream, "%d\t", flag7);
3916                 fprintf(stream, "%d\t", flag8);
3917                 fprintf(stream, "%d\t", flag9);
3918                 fprintf(stream, "%d\t", flag10);
3919                 fprintf(stream, "%d\t", flag11);
3920                 fprintf(stream, "%d\t", flag12);
3921                 fprintf(stream, "%d\t", flag13);
3922                 fprintf(stream, "%d\t", flag14);
3923                 fprintf(stream, "%d\t", flag15);
3924                 fprintf(stream, "%d\t", flag16);
3925                 fprintf(stream, "%d\t", flag17);
3926                 fprintf(stream, "%d\t", flag18);
3927                 fprintf(stream, "%d\t", flag19);
3928                 fprintf(stream, "%d\t", flag20);
3929                 fprintf(stream, "%d\t", flagg);
3930                 fprintf(stream, "%d\t", flagg2);
3931                 fprintf(stream, "%d\t", cantidad);
3932                 fprintf(stream, "%d\t", cantidadb);
3933                 fprintf(stream, "%d\t", burned);
3934
3935
3936                 /*
3937                 fprintf(stream, "caf[0] %d\t", caf[0]);
3938                 fprintf(stream, "caf[1] %d\t", caf[1]);
3939                 fprintf(stream, "caf[2] %d\t", caf[2]);
3940                 fprintf(stream, "caf[3] %d\t", caf[3]);
3941                 fprintf(stream, "REL%f\t", ((double)caf[0]+(double)caf[1])/((double)caf[0]+(double)caf[1]+(double)caf[2]+(double)caf[3]));
3942                 fprintf(stream, "REL2%f\t", ((double)caf[1])/((double)caf[3]+(double)caf[1]));
3943                 fprintf(stream, "REL3%f\t", ((double)caf[1])/((double)caf[1]+(double)caf[2]));
3944
3945                 fprintf(stream, "maxareacandit%d\t", vector[indmaxarea][0]);
3946                 fprintf(stream, "totalsurf%d\t", totalsurf);
3947                 fprintf(stream, "conecsurf%d\t", conecsurf);
3948                 fprintf(stream, "ratio%f\t", (double)conecsurf/(double)totalsurf);
3949                 fprintf(stream, "ratio2%f\t", (double)vector[indmaxarea][0]/(double)totalsurf);
3950                 */
3951
3952
3953                 fprintf(stream, "[%f, %f, %f]\t", puntoactual[0], puntoactual[1], puntoactual[2]);
3954                 fprintf(stream, "[%f, %f, %f]\t", puntoanterior[0], puntoanterior[1], puntoanterior[2]);
3955                 fprintf(stream, "[%f, %f, %f]\t", puntopre[0], puntopre[1], puntopre[2]);
3956                 fprintf(stream, "[%f, %f, %f]\t", mejor[0], mejor[1], mejor[2]);
3957
3958                 fprintf(stream, "\n");
3959
3960         }
3961
3962 }
3963
3964
3965
3966
3967 //----------------------------------------------------------------------------
3968 void axisExtractor02::todo()
3969 {
3970         int i=0;
3971
3972         while( !m_Stack0.empty() && i<5000){
3973                 avanzar();
3974                 i++;
3975         }
3976 }
3977
3978
3979
3980 //----------------------------------------------------------------------------
3981 void axisExtractor02::paso()
3982 {
3983
3984         int i=0;
3985         int inicio=buenos;
3986
3987         while( !m_Stack0.empty() && inicio==buenos && i<5000){
3988                 avanzar();
3989                 i++;
3990         }
3991 }
3992
3993
3994
3995
3996 //----------------------------------------------------------------------------
3997 void axisExtractor02::rama()
3998 {
3999
4000         int i=0;
4001         frama=0;
4002
4003         while( !m_Stack0.empty() && frama==0 && i<5000){
4004                 avanzar();
4005                 i++;
4006         }
4007 }
4008
4009
4010
4011
4012 //----------------------------------------------------------------------------
4013 void axisExtractor02::segmento()
4014 {
4015
4016         int i=0;
4017         fseg=0;
4018         frama=0;
4019
4020         while( !m_Stack0.empty() && frama==0 && fseg==0 && i<5000){
4021                 avanzar();
4022                 i++;
4023         }
4024 }
4025
4026