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

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