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