lib/maracasVisuLib/src/kernel/carotidaBifurcacion.cxx

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