lib/maracasVisuLib/src/interface/wxWindows/Contour/Propagation.cxx

   1 #include "Propagation.h"
   2
   3 //------------------------------------------------------------------------------------------------------------------------------------------
   4 //CLASS: Vector ------------------------------------------------------------------------------------------------------------------------
   5 //------------------------------------------------------------------------------------------------------------------------------------------
   6 //Constructor
   7 Vector::Vector()
   8 {
   9         _plane = -1;
  10         _var = -1;
  11 }
  12 //Destructor
  13 Vector::~Vector()
  14 {
  15         _vec.clear();
  16         _vecX.clear();
  17         _vecY.clear();
  18         _vecZ.clear();
  19 }
  20
  21 //------------------------------------------------------------------------------------------------------------------------------------------
  22 void Vector::set_vec(double val)
  23 {
  24         _vec.push_back(val);
  25 }
  26 //------------------------------------------------------------------------------------------------------------------------------------------
  27 void Vector::set_var(double val)
  28 {
  29         _var = val;
  30 }
  31 //------------------------------------------------------------------------------------------------------------------------------------------
  32 double Vector::get_vec(int id)
  33 {
  34         if(_vec.size() != 0)
  35         {
  36                 return _vec[id];
  37         }
  38         return -1;
  39 }
  40 //------------------------------------------------------------------------------------------------------------------------------------------
  41 double Vector::get_var()
  42 {
  43         return _var;
  44 }
  45 //------------------------------------------------------------------------------------------------------------------------------------------
  46 int Vector::getsize_vec()
  47 {
  48         return _vec.size();
  49 }
  50 //------------------------------------------------------------------------------------------------------------------------------------------
  51 void Vector::copyVector( std::vector<Vector>*vec1, std::vector<Vector>*vec2 )
  52 {
  53         int i,j;
  54         if( vec1->size() != 0 )
  55         {
  56                 vec2->clear();
  57                 for(i=0; i<vec1->size(); i++)
  58                 {
  59                         Vector *temp = new Vector();
  60                         temp->set_var( (*vec1)[i].get_var() );
  61                         for(j=0; j<(*vec1)[i].getsize_vec(); j++)
  62                         {
  63                                 temp->set_vec( (*vec1)[i].get_vec(j) );
  64                         }
  65                         vec2->push_back(*temp);
  66                         delete temp;
  67                 }
  68         }
  69 }
  70 //------------------------------------------------------------------------------------------------------------------------------------------
  71 void Vector::printVector(std::vector<Vector>*vec1)
  72 {
  73         int i,j;
  74         for(i=0; i<vec1->size(); i++)
  75         {
  76                 printf("\n Pos (%d) => var = %f",i,(*vec1)[i].get_var());
  77                 for(j=0; j<(*vec1)[i].getsize_vec(); j++)
  78                 {
  79                         printf("\n vec(%d) = %f",j,(*vec1)[i].get_vec(j));
  80                 }
  81         }
  82 }
  83 //------------------------------------------------------------------------------------------------------------------------------------------
  84 void Vector::set_plane(int val)
  85 {
  86         _plane = val;
  87 }
  88 //------------------------------------------------------------------------------------------------------------------------------------------
  89 void Vector::set_x(double val)
  90 {
  91         _vecX.push_back(val);
  92 }
  93 //------------------------------------------------------------------------------------------------------------------------------------------
  94 //------------------------------------------------------------------------------------------------------------------------------------------
  95 void Vector::set_y(double val)
  96 {
  97         _vecY.push_back(val);
  98 }
  99 //------------------------------------------------------------------------------------------------------------------------------------------
 100 //------------------------------------------------------------------------------------------------------------------------------------------
 101 void Vector::set_z(double val)
 102 {
 103         _vecZ.push_back(val);
 104 }
 105 //------------------------------------------------------------------------------------------------------------------------------------------
 106 int Vector::get_plane()
 107 {
 108         return _plane;
 109 }
 110 //------------------------------------------------------------------------------------------------------------------------------------------
 111 double Vector::get_x(int id)
 112 {
 113         if( (-1<id) && (id<_vecX.size()) )
 114         {
 115                 return _vecX[id];
 116         }
 117         return -1;
 118 }
 119 //------------------------------------------------------------------------------------------------------------------------------------------
 120 int Vector::getsize_x()
 121 {
 122         if(_vecX.size() != 0)
 123         {
 124                 return _vecX.size();
 125         }
 126         return -1;
 127 }
 128 //------------------------------------------------------------------------------------------------------------------------------------------
 129 double Vector::get_y(int id)
 130 {
 131         if( (-1<id) && (id<_vecY.size()) )
 132         {
 133                 return _vecY[id];
 134         }
 135         return -1;
 136 }
 137 //------------------------------------------------------------------------------------------------------------------------------------------
 138 int Vector::getsize_y()
 139 {
 140         if(_vecY.size() != 0)
 141         {
 142                 return _vecY.size();
 143         }
 144         return -1;
 145 }
 146 //------------------------------------------------------------------------------------------------------------------------------------------
 147 int Vector::getsize_z()
 148 {
 149         if(_vecZ.size() != 0)
 150         {
 151                 return _vecZ.size();
 152         }
 153         return -1;
 154 }
 155 //------------------------------------------------------------------------------------------------------------------------------------------
 156 double Vector::get_z(int id)
 157 {
 158         if( (-1<id) && (id<_vecZ.size()) )
 159         {
 160                 return _vecZ[id];
 161         }
 162         return -1;
 163 }
 164 //------------------------------------------------------------------------------------------------------------------------------------------
 165 std::vector<double> Vector::getVec()
 166 {
 167         return _vec;
 168 }
 169 //------------------------------------------------------------------------------------------------------------------------------------------
 170 void Vector::resetVec()
 171 {
 172         _vec.clear();
 173 }
 174 //------------------------------------------------------------------------------------------------------------------------------------------
 175 //------------------------------------------------------------------------------------------------------------------------------------------
 176 //------------------------------------------------------------------------------------------------------------------------------------------
 177
 178
 179 //Constructor
 180 PropContour::PropContour()
 181 {
 182         _interpnumber = 100;
 183 }
 184
 185 //Destructor
 186 PropContour::~PropContour()
 187 {
 188         ResetPlaneVector();
 189         ResetKeyContours();
 190 }
 191
 192 /* EED 03/07/2008
 193 //------------------------------------------------------------------------------------
 194 double PropContour::RBF_WendLand(double norm, double m_rad)
 195 {
 196         double y;
 197           norm = norm / m_rad;
 198           y = pow( 1-norm,4 ) * ( (4*norm) + 1 );
 199           if(norm >= 1)
 200           {
 201                 y = 0;
 202           }
 203         return y;
 204 }
 205 //------------------------------------------------------------------------------------
 206 double PropContour::RBF_ThinPlate(double norm)
 207 {
 208         double y;
 209         if(norm == 0)
 210         {
 211                 y = 0;
 212         }
 213         else
 214         {
 215                 y = pow(norm,2)*log(norm);
 216         }
 217         return y;
 218 }
 219 //------------------------------------------------------------------------------------
 220 vtkImageData* PropContour::method_RBF ( double rad, std::vector<double>*CoordX, std::vector<double>*CoordY,
 221                                                                                 std::vector<double>*CoordZ )
 222 {
 223         _dimImage[0] = 100;                     // X Axis
 224         _dimImage[1] = 100;                     // Y Axis
 225         _dimImage[2] = 1;                       // Z Axis
 226         int pointsSize = CoordX->size();
 227         double spc[3]={1,1,1};
 228         double norm = 0, val = 0;
 229
 230         vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
 231         vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
 232         vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
 233         vnl_vector<double> H( pointsSize, 1 );
 234         vnl_vector<double> D( pointsSize, 0.0 );
 235         vnl_vector<double> Q( 2, 0.0 );
 236         vnl_matrix<double> Impl( _dimImage[0],_dimImage[1], 0.0 );
 237
 238         unsigned short *pValue;
 239         imagedataValue = vtkImageData::New();
 240         imagedataValue->SetScalarTypeToUnsignedShort();
 241         imagedataValue->SetSpacing(spc);
 242         imagedataValue->SetDimensions(_dimImage);
 243         imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
 244         imagedataValue->AllocateScalars();
 245         imagedataValue->Update();
 246
 247     int i,j,h;
 248         for(i=0; i<pointsSize; i++)
 249         {
 250                 for(j=0; j<pointsSize; j++)
 251                 {
 252                         norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) );
 253                         A(i,j) = RBF_WendLand(norm,rad);
 254                         //A(i,j) = RBF_ThinPlate(norm);
 255                 }
 256         }
 257
 258         Ainv = vnl_matrix_inverse<double>(A);
 259         D = Ainv* H;
 260
 261         for(i=0; i < _dimImage[1]; i++)
 262         {
 263                 for(j=0; j < _dimImage[0]; j++)
 264                 {
 265                         Q(0) = i;
 266                         Q(1) = j;
 267                         val = 0;
 268                         for(h=0; h<pointsSize; h++)
 269                         {
 270                                 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) );
 271                                 //val = val + ( D(h) * RBF_ThinPlate(norm) );
 272                                 val = val + ( D(h) * RBF_WendLand(norm,rad) );
 273                         }
 274                         Impl(i,j) = val - 1.0;
 275
 276 //                      if ( (Impl(i,j)>=-0.001) && (Impl(i,j) <=0.001)  )
 277 //                      {Impl(i,j)=128;
 278 //                      }
 279 //                      else
 280 //                      {Impl(i,j)=0;}
 281
 282
 283 //                      pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0);
 284 //                      vvalue = (Impl(i,j)*256+256);
 285 //                      if (vvalue) < 0)
 286 //                      {
 287 //                              *pValue = 0;
 288 //                      } else {
 289 //                              *pValue = vvalue;
 290 //                      }
 291
 292                         pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,0);
 293                         if ( Impl(i,j) >=0 )
 294                         {
 295                                 *pValue = 128;
 296                         }
 297                         else
 298                         {
 299                                 *pValue = 0;
 300                         }
 301                 }
 302         }
 303
 304         for(i=0; i<pointsSize; i=i+5)
 305         {
 306          pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,0);
 307                  *pValue=255;
 308         }
 309
 310         return imagedataValue;
 311 }
 312 //----------------------------------------------------------------------------------------------------
 313 double PropContour::RBF_ThinPlate_3D(double norm)
 314 {
 315         return norm = pow( norm,3 );
 316 }
 317 //----------------------------------------------------------------------------------------------------
 318 vtkImageData* PropContour::method_RBF_3D (  double rad, std::vector<double>*CoordX, std::vector<double>*CoordY,
 319                                                                                         std::vector<double>*CoordZ )
 320 {
 321         long interval = wxGetElapsedTime(TRUE);
 322
 323         int i,j,k,h;
 324         int pointsSize = CoordX->size();
 325         double max = -1, min = 1000, minz = 100000, maxz = -10000;
 326
 327         for( i=0; i<pointsSize; i++ )
 328         {
 329                 if( (*CoordX)[i] > max )
 330                 {
 331                         max = (*CoordX)[i];
 332                 }
 333                 if( (*CoordY)[i] > max )
 334                 {
 335                         max = (*CoordY)[i];
 336                 }
 337                 if( (*CoordX)[i] < min )
 338                 {
 339                         min = (*CoordX)[i];
 340                 }
 341                 if( (*CoordY)[i] < min )
 342                 {
 343                         min = (*CoordY)[i];
 344                 }
 345                 if( (*CoordZ)[i] < minz )
 346                 {
 347                         minz = (*CoordZ)[i];
 348                 }
 349                 if( (*CoordZ)[i] > maxz )
 350                 {
 351                         maxz = (*CoordZ)[i];
 352                 }
 353         }
 354
 355         _dimImage[0] = 200;                     // X axis
 356         _dimImage[1] = 200;                     // Y axis
 357         _dimImage[2] = 200;                     // Z axis
 358
 359         double spc[3]={1,1,1};
 360         double norm = 0, val, vvalue;
 361
 362         vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
 363         vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
 364         vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
 365         vnl_vector<double> H( pointsSize, 1 );
 366         vnl_vector<double> D( pointsSize, 0.0 );
 367         vnl_vector<double> Q( 3, 0.0 );
 368
 369         unsigned short *pValue;
 370         imagedataValue = vtkImageData::New();
 371         imagedataValue->SetScalarTypeToUnsignedShort();
 372         imagedataValue->SetSpacing(spc);
 373         imagedataValue->SetDimensions(_dimImage);
 374         imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
 375         imagedataValue->AllocateScalars();
 376         imagedataValue->Update();
 377
 378         for(i=0; i<pointsSize; i++)
 379         {
 380                 for(j=0; j<pointsSize; j++)
 381                 {
 382                         norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) );
 383                         A(i,j) = RBF_WendLand(norm,rad);
 384                 }
 385         }
 386         D = vnl_matrix_inverse<double>(A)* H;
 387
 388 //Inicialization
 389 //      for(k=0; k<_dimImage[2]; k++)
 390 //      {
 391 //              for(j=0; j<_dimImage[1]; j++)
 392 //              {
 393 //                      for(i=0; i<_dimImage[0]; i++)
 394 //                      {
 395 //                              pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
 396 //                              *pValue = 0;
 397 //                      }
 398 //              }
 399 //      }
 400
 401 //Filling
 402 //      pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0);
 403 //      i = 0;
 404 //      j = 0;
 405 //      k = 0;
 406 //      for(h=0; h<pointsSize; h++)
 407 //      {
 408 //              norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
 409 //              //val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
 410 //              val = val + ( D(h) * RBF_WendLand(norm,rad) );
 411 //              if( h == pointsSize-1 )
 412 //              {
 413 //                      val = val - 1;
 414 //                      pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
 415 //                      if ( val >=0 )
 416 //                      {
 417 //                              *pValue = 128;
 418 //                      }
 419 //                      else
 420 //                      {
 421 //                              *pValue = 0;
 422 //                      }
 423 //                      //pValue++;
 424 //                      j++;
 425 //                      h = -1;
 426 //                      val = 0;
 427 //                      if( j == _dimImage[0] )
 428 //                      {
 429 //                              i++;
 430 //                              j = 0;
 431 //                              h = -1;
 432 //                              if( i == _dimImage[1])
 433 //                              {
 434 //                                      k++;
 435 //                                      i = 0;
 436 //                                      j = 0;
 437 //                                      h = -1;
 438 //                                      if(k == _dimImage[2])
 439 //                                      {
 440 //                                              h = pointsSize;
 441 //                                      }
 442 //                              }
 443 //                      }
 444 //              }
 445 //      }
 446
 447         i = (int)min-10;
 448         j = (int)min-10;
 449         k = (int)minz-10;
 450         val = 0;
 451
 452         for(h=0; h<pointsSize; h++)
 453         {
 454                 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
 455                 val = val + ( D(h) * RBF_WendLand(norm,rad) );
 456                 if( h == pointsSize-1 )
 457                 {
 458                         val = val - 1;
 459                         pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
 460                         if ( val >=0 )
 461                         {
 462                                 *pValue = 128;
 463                         }
 464                         else
 465                         {
 466                                 *pValue = 0;
 467                         }
 468                         //pValue++;
 469                         j++;
 470                         h = -1;
 471                         val = 0;
 472                         if( j == (int)max+10)
 473                         {
 474                                 i++;
 475                                 j = min-10;
 476                                 h = -1;
 477                                 if( i == (int)max+10)
 478                                 {
 479                                         k++;
 480                                         i = min-10;
 481                                         j = min-10;
 482                                         h = -1;
 483                                         if(k == maxz+10)
 484                                         {
 485                                                 h = pointsSize;
 486                                         }
 487                                 }
 488                         }
 489                 }
 490         }
 491
 492         interval = wxGetElapsedTime(FALSE);
 493         long interPlane = interval/_dimImage[2];
 494
 495         for(i=0; i<pointsSize; i++)
 496         {
 497                  pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]);
 498                 *pValue=255;
 499         }
 500
 501         long intervalPC = wxGetElapsedTime();
 502
 503         printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------");
 504         printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval);
 505         printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
 506         printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane);
 507         printf("\n NUMBER OF PLANES................................. %d",k);
 508         printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
 509         printf("\n ------------------------------------------------------------");
 510
 511         return imagedataValue;
 512 }
 513 //----------------------------------------------------------------------------------------------------
 514
 515 vtkImageData* PropContour::method_RBF_3D_ThinPlate (  double rad, std::vector<double>*CoordX, std::vector<double>*CoordY,
 516                                                                                                           std::vector<double>*CoordZ )
 517 {
 518         long interval = wxGetElapsedTime(TRUE);
 519
 520         int i,j,k,h;
 521         int pointsSize = CoordX->size();
 522         double max = -1, min = 1000, minz = 100000, maxz = -10000;
 523
 524         for( i=0; i<pointsSize; i++ )
 525         {
 526                 if( (*CoordX)[i] > max )
 527                 {
 528                         max = (*CoordX)[i];
 529                 }
 530                 if( (*CoordY)[i] > max )
 531                 {
 532                         max = (*CoordY)[i];
 533                 }
 534                 if( (*CoordX)[i] < min )
 535                 {
 536                         min = (*CoordX)[i];
 537                 }
 538                 if( (*CoordY)[i] < min )
 539                 {
 540                         min = (*CoordY)[i];
 541                 }
 542                 if( (*CoordZ)[i] < minz )
 543                 {
 544                         minz = (*CoordZ)[i];
 545                 }
 546                 if( (*CoordZ)[i] > maxz )
 547                 {
 548                         maxz = (*CoordZ)[i];
 549                 }
 550         }
 551
 552         _dimImage[0] = 190;                     // X axis
 553         _dimImage[1] = 190;                     // Y axis
 554         _dimImage[2] = 190;                     // Z axis
 555
 556         double spc[3]={1,1,1};
 557         double norm = 0, val, vvalue;
 558
 559         vnl_matrix<double> A( pointsSize,pointsSize, 0.0 );
 560         vnl_matrix<double> Ainv( pointsSize,pointsSize, 0.0 );
 561         vnl_matrix<double> I( pointsSize,pointsSize, 0.0 );
 562         vnl_vector<double> H( pointsSize, 1 );
 563         vnl_vector<double> D( pointsSize, 0.0 );
 564         vnl_vector<double> Q( 3, 0.0 );
 565
 566         unsigned short *pValue;
 567         imagedataValue = vtkImageData::New();
 568         imagedataValue->SetScalarTypeToUnsignedShort();
 569         imagedataValue->SetSpacing(spc);
 570         imagedataValue->SetDimensions(_dimImage);
 571         imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,-10,_dimImage[2]-1);
 572         imagedataValue->AllocateScalars();
 573         imagedataValue->Update();
 574
 575         for(i=0; i<pointsSize; i++)
 576         {
 577                 for(j=0; j<pointsSize; j++)
 578                 {
 579                         norm = sqrt( pow((*CoordX)[i]-(*CoordX)[j],2) + pow((*CoordY)[i]-(*CoordY)[j],2) + pow((*CoordZ)[i]-(*CoordZ)[j],2) );
 580                         A(i,j) = RBF_ThinPlate_3D(norm);
 581                 }
 582         }
 583         D = vnl_matrix_inverse<double>(A)* H;
 584
 585 //Inicialization
 586 //      for(k=0; k<_dimImage[2]; k++)
 587 //      {
 588 //              for(j=0; j<_dimImage[1]; j++)
 589 //              {
 590 //                      for(i=0; i<_dimImage[0]; i++)
 591 //                      {
 592 //                              pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
 593 //                              *pValue = 0;
 594 //                      }
 595 //              }
 596 //      }
 597
 598 //Filling
 599
 600 //      pValue = (unsigned short *)imagedataValue->GetScalarPointer(0,0,0);
 601 //      i = 0;
 602 //      j = 0;
 603 //      k = 0;
 604 //      for(h=0; h<pointsSize; h++)
 605 //      {
 606 //              norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
 607 //              //val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
 608 //              val = val + ( D(h) * RBF_WendLand(norm,rad) );
 609 //              if( h == pointsSize-1 )
 610 //              {
 611 //                      val = val - 1;
 612 //                      pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
 613 //                      if ( val >=0 )
 614 //                      {
 615 //                              *pValue = 128;
 616 //                      }
 617 //                      else
 618 //                      {
 619 //                              *pValue = 0;
 620 //                      }
 621 //                      //pValue++;
 622 //                      j++;
 623 //                      h = -1;
 624 //                      val = 0;
 625 //                      if( j == _dimImage[0] )
 626 //                      {
 627 //                              i++;
 628 //                              j = 0;
 629 //                              h = -1;
 630 //                              if( i == _dimImage[1])
 631 //                              {
 632 //                                      k++;
 633 //                                      i = 0;
 634 //                                      j = 0;
 635 //                                      h = -1;
 636 //                                      if(k == _dimImage[2])
 637 //                                      {
 638 //                                              h = pointsSize;
 639 //                                      }
 640 //                              }
 641 //                      }
 642 //              }
 643 //      }
 644
 645         i = (int)min-10;
 646         j = (int)min-10;
 647         k = (int)minz-10;
 648         val = 0;
 649
 650         for(h=0; h<pointsSize; h++)
 651         {
 652                 norm = sqrt( pow( i-(*CoordX)[h],2 ) + pow( j-(*CoordY)[h],2 ) + pow( k-(*CoordZ)[h],2 ));
 653                 val = val + ( D(h) * RBF_ThinPlate_3D(norm) );
 654                 if( h == pointsSize-1 )
 655                 {
 656                         val = val - 1;
 657                         pValue = (unsigned short *)imagedataValue->GetScalarPointer(i,j,k);
 658                         if ( val >=0 )
 659                         {
 660                                 *pValue = 128;
 661                         }
 662                         else
 663                         {
 664                                 *pValue = 0;
 665                         }
 666                         //pValue++;
 667                         j++;
 668                         h = -1;
 669                         val = 0;
 670                         if( j == (int)max+10)
 671                         {
 672                                 i++;
 673                                 j = min-10;
 674                                 h = -1;
 675                                 if( i == (int)max+10)
 676                                 {
 677                                         k++;
 678                                         i = min-10;
 679                                         j = min-10;
 680                                         h = -1;
 681                                         if(k == maxz+10)
 682                                         {
 683                                                 h = pointsSize;
 684                                         }
 685                                 }
 686                         }
 687                 }
 688         }
 689
 690         interval = wxGetElapsedTime(FALSE);
 691         long interPlane = interval/_dimImage[2];
 692
 693         for(i=0; i<pointsSize; i++)
 694         {
 695                  pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*CoordX)[i] , (int)(*CoordY)[i] ,(int)(*CoordZ)[i]);
 696                 *pValue=255;
 697         }
 698
 699         long intervalPC = wxGetElapsedTime();
 700
 701         printf("\n\n JSTG - PropContour::method_RBF_3D ------------------------");
 702         printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %lld (ms)",interval);
 703         printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
 704         printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %lld (ms)",interPlane);
 705         printf("\n NUMBER OF PLANES................................. %d",k);
 706         printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
 707         printf("\n ------------------------------------------------------------");
 708
 709         return imagedataValue;
 710 }
 711
 712 */
 713
 714 //---------------------------------------------------------------------------------------------------------
 715 void PropContour::ReadKeyContour(FILE* fd)
 716 {
 717         char firstline[30];
 718         int     size;
 719         double x,y;
 720         int z;
 721         std::vector<double> tempX;
 722         std::vector<double> tempY;
 723         std::vector<double> tempZ;
 724         tempX.clear();
 725         tempY.clear();
 726         tempZ.clear();
 727         while(!feof(fd))
 728         {
 729                 fscanf(fd," %s %d",&firstline,&size);
 730                 for(int i=0; i<size; i++)
 731                 {
 732                         fscanf(fd,"%lf %lf %d",&x,&y,&z);
 733                         tempX.push_back(x);
 734                         tempY.push_back(y);
 735                         tempZ.push_back(z);
 736                 }
 737                 SetKeyContours(&tempX,&tempY,&tempZ);
 738                 tempX.clear();
 739                 tempY.clear();
 740                 tempZ.clear();
 741         }
 742 }
 743 //---------------------------------------------------------------------------------------------------------
 744 int PropContour::VectorDirection(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
 745 {
 746         int dir,i;
 747         double SumX = 0,SumY = 0;
 748         double ax,ay,bx,by,axb;
 749         int size = InX->size();
 750         for(i=0; i<size; i++)
 751         {
 752                 SumX = SumX + (*InX)[i];
 753                 SumY = SumY + (*InY)[i];
 754         }
 755         SumX = SumX/size;               //Mass Center: X coord
 756         SumY = SumY/size;               //Mass Center: Y coord
 757
 758         int positive = 0;
 759         int negative = 0;
 760         for(i=0; i<size; i++)
 761         {
 762                 ax = (*InX)[i]-SumX;
 763                 ay = (*InY)[i]-SumY;
 764                 bx = (*InX)[i+1]-SumX;
 765                 by = (*InY)[i+1]-SumY;
 766                 axb = (ax*by) - (bx*ay);
 767                 if(axb > 0)
 768                 {
 769                         positive++;
 770                 }
 771                 if(axb < 0)
 772                 {
 773                         negative++;
 774                 }
 775         }
 776         if(positive >= negative)
 777         {
 778                 dir = 1;
 779         }
 780         else
 781         {
 782                 dir = -1;
 783         }
 784
 785         return dir;
 786 }
 787 //----------------------------------------------------------------------------------------------------
 788 void PropContour::VectorOrder(int dir, int posini, std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
 789 {
 790         int size = InX->size();
 791         int i;
 792         std::vector<double> tempX;
 793         std::vector<double> tempY;
 794         std::vector<double> tempZ;
 795
 796         tempX.clear();
 797         tempY.clear();
 798         tempZ.clear();
 799
 800         for(i=0; i<size; i++)
 801         {
 802                 if(dir == 1)
 803                 {
 804                         tempX.push_back((*InX)[posini]);
 805                         tempY.push_back((*InY)[posini]);
 806                         tempZ.push_back((*InZ)[posini]);
 807                         posini++;
 808                         if(posini == size)
 809                         {
 810                                 posini = 0;
 811                         }
 812                 }
 813                 if(dir == -1)
 814                 {
 815                         tempX.push_back((*InX)[posini]);
 816                         tempY.push_back((*InY)[posini]);
 817                         tempZ.push_back((*InZ)[posini]);
 818                         posini--;
 819                         if(posini < 0)
 820                         {
 821                                 posini = size-1;
 822                         }
 823                 }
 824         }
 825         InX->clear();
 826         InY->clear();
 827         InZ->clear();
 828         for(i=0; i<size; i++)
 829         {
 830                 InX->push_back(tempX[i]);
 831                 InY->push_back(tempY[i]);
 832                 InZ->push_back(tempZ[i]);
 833         }
 834
 835 }
 836 //----------------------------------------------------------------------------------------------------
 837 void PropContour::PreparePointsForSpline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ,
 838                                                                                  std::vector<int>*Sizes)
 839 {
 840         int sizeS = Sizes->size();
 841         int sizeV = InX->size();
 842         int i,j,mem,posinic,dir,cont;
 843         double leX;
 844
 845         std::vector<double> tempX;
 846         std::vector<double> tempY;
 847         std::vector<double> tempZ;
 848         std::vector<double> lstX;
 849         std::vector<double> lstY;
 850         std::vector<double> lstZ;
 851
 852         lstX.clear();
 853         lstY.clear();
 854         lstZ.clear();
 855
 856         mem = 0;
 857         cont = 0;
 858         for(i=0; i<sizeS; i++)
 859         {
 860                 leX=1000;
 861                 tempX.clear();
 862                 tempY.clear();
 863                 tempZ.clear();
 864                 for(j=0; j<(*Sizes)[i]; j++)
 865                 {
 866                         tempX.push_back((*InX)[j+mem]);
 867                         tempY.push_back((*InY)[j+mem]);
 868                         tempZ.push_back((*InZ)[j+mem]);
 869                         if( (*InX)[j] < leX )
 870                         {
 871                                 posinic = j;
 872                                 leX = (*InX)[j];
 873                         }
 874                 }
 875                 mem = mem + (*Sizes)[i];
 876                 dir = VectorDirection(&tempX,&tempY,&tempZ);
 877                 VectorOrder(dir,posinic,&tempX,&tempY,&tempZ);
 878
 879                 for(j=0; j<(*Sizes)[i]; j++)
 880                 {
 881                         lstX.push_back(tempX[j]);
 882                         lstY.push_back(tempY[j]);
 883                         lstZ.push_back((*InZ)[j+cont]);
 884                 }
 885                 cont = cont + (*Sizes)[i];
 886         }
 887
 888 //Fill the Finally lst in X,Y,Z ---------------
 889         int sizetemp = lstX.size();
 890         //printf("\nJSTG-PropContour::PreparePointsForSpline");
 891         InX->clear();
 892         InY->clear();
 893         InZ->clear();
 894         for(i=0; i<sizetemp; i++)
 895         {
 896                 //printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
 897                 InX->push_back(lstX[i]);
 898                 InY->push_back(lstY[i]);
 899                 InZ->push_back(lstZ[i]);
 900         }
 901
 902
 903 }
 904 //----------------------------------------------------------------------------------------------------
 905 vtkImageData* PropContour::method_Spline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ, std::vector<int>*Sizes)
 906 {
 907         long interval = wxGetElapsedTime(TRUE);
 908
 909         int i,j,k,sizeX,sizeS,sizeInS;
 910         int numspline;
 911         double x,y,z;
 912 //EED   double spc[3]={1,1,1};
 913         std::vector<double> lstX;
 914         std::vector<double> lstY;
 915         std::vector<double> lstZ;
 916         std::vector<double> tempX;
 917         std::vector<double> tempY;
 918         std::vector<double> tempZ;
 919         std::vector<double> tempS;
 920         _mContourModel = new manualContourModel();
 921         _mContourModel->SetNumberOfPointsSpline(_interpnumber);
 922
 923         imagedataValue=NULL;
 924 //EED
 925 //      _dimImage[0] = 200;                     // X axis
 926 //      _dimImage[1] = 200;                     // Y axis
 927 //      _dimImage[2] = 200;                     // Z axis
 928 //
 929 //      unsigned short *pValue;
 930 //      imagedataValue = vtkImageData::New();
 931 //      imagedataValue->SetScalarTypeToUnsignedShort();
 932 //      imagedataValue->SetSpacing(spc);
 933 //      imagedataValue->SetDimensions(_dimImage);
 934 //      imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
 935 //      imagedataValue->AllocateScalars();
 936 //      imagedataValue->Update();
 937
 938 //      lstX.clear();
 939 //      lstY.clear();
 940 //      lstZ.clear();
 941
 942         sizeX = InX->size();
 943         for(i=0; i<sizeX; i++)
 944         {
 945                 lstX.push_back((*InX)[i]);
 946                 lstY.push_back((*InY)[i]);
 947                 lstZ.push_back((*InZ)[i]);
 948         }
 949
 950         PreparePointsForSpline(&lstX,&lstY,&lstZ,Sizes);
 951         /*int sizetemp = lstX.size();
 952         printf("\nJSTG-PropContour::method_Spline");
 953         for(i=0; i<sizetemp; i++)
 954         {
 955                 printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
 956         }*/
 957
 958         sizeS = Sizes->size();
 959         int cont = 0;
 960         for(i=0; i<sizeS; i++)
 961         {
 962                 _mContourModel->DeleteAllPoints();
 963                 _mContourModel->SetCloseContour(true);
 964                 sizeInS = (*Sizes)[i];
 965                 for(j=0; j<sizeInS; j++)
 966                 {
 967                         _mContourModel->AddPoint(lstX[cont],lstY[cont],lstZ[cont]);
 968                         cont ++;
 969                 }
 970                 _mContourModel->UpdateSpline();
 971                 numspline = _mContourModel->GetNumberOfPointsSpline();
 972                 for(k=0; k<numspline; k++)
 973                 {
 974                         _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
 975 //EED                   pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
 976 //                      *pValue = 128;
 977                         tempX.push_back(x);
 978                         tempY.push_back(y);
 979                         tempZ.push_back(z);
 980                 }
 981         }
 982
 983         int tam = numspline;
 984         std::vector<Vector>::iterator it;
 985
 986         ResetPlaneVector();
 987
 988         for(i=0; i<numspline; i++)
 989         {
 990                 Vector *vec = new Vector();
 991                 _planevector.push_back(*vec);
 992         }// for i
 993
 994         for(j=0; j<tam; j++)
 995         {
 996                 _mContourModel->DeleteAllPoints();
 997                 _mContourModel->SetCloseContour(false);
 998                 cont = 0;
 999                 for(i=0; i<sizeS; i++)
1000                 {
1001                         double hh=tempZ[cont+j];
1002                         _mContourModel->AddPoint(tempX[cont+j],tempY[cont+j],tempZ[cont+j]);
1003                         cont = cont + tam;
1004                 } // for i
1005                 _mContourModel->UpdateSpline();
1006                 numspline = _mContourModel->GetNumberOfPointsSpline();
1007                 for(k=0; k<numspline; k++)
1008                 {
1009 //EED002
1010                         _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
1011                         _planevector[k].set_x(x);
1012                         _planevector[k].set_y(y);
1013                         _planevector[k].set_z(z);
1014                         _planevector[k].set_plane(k);
1015 //EED                   pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
1016 //EED                   *pValue = 128;
1017                 } // for k
1018         } // for j
1019
1020 /*
1021         int tempsize = _planevector.size();
1022         int sizexx;
1023         for(i=0; i<tempsize; i++)
1024         {
1025                 sizexx = _planevector[i].getsize_x();
1026         }
1027 */
1028 /*
1029         interval = wxGetElapsedTime(FALSE);
1030         int zplanes = abs(tempZ[0]-tempZ[tempZ.size()-1]);
1031         long interPlane = interval/zplanes;
1032         int pointsSize = 0;
1033         for(i=0; i<sizeS; i++)
1034         {
1035                 pointsSize = pointsSize + (*Sizes)[i];
1036         }
1037         for(i=0; i<pointsSize; i++)
1038         {
1039                  pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*InX)[i] , (int)(*InY)[i] ,(int)(*InZ)[i]);
1040                 *pValue=255;
1041         }
1042
1043         printf("\n\n JSTG - PropContour::method_Spline ------------------------");
1044         printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %d (ms)",interval);
1045         //printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
1046         printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %d (ms)",interPlane);
1047         printf("\n NUMBER OF PLANES................................. %d",zplanes);
1048         printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
1049         printf("\n ------------------------------------------------------------");
1050 */
1051         return imagedataValue;
1052 }
1053
1054 //---------------------------------------------------------------------------------------------------
1055 void PropContour::ResetPlaneVector()
1056 {
1057 //      Vector *vec;
1058 //      int ii,iiSize = _planevector.size();
1059 //      for (ii=0 ; ii<iiSize ; ii++)
1060 //      {
1061 //              vec = &(_planevector[ii]);
1062 //              delete vec;
1063 //      }
1064         _planevector.clear();
1065 }
1066
1067
1068 //---------------------------------------------------------------------------------------------------
1069 void PropContour::ResetKeyContours()
1070 {
1071         _KeyContourSizes.clear();
1072         _KeyContourX.clear();
1073         _KeyContourY.clear();
1074         _KeyContourZ.clear();
1075 }
1076 //---------------------------------------------------------------------------------------------------
1077 void PropContour::SetKeyContours(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
1078 {
1079
1080         int idKeyContour                = 0;
1081         int idKeyContourSizes   = 0;
1082         int tmpIdKeyContSizes   = 0;
1083         bool okFind                             = false;
1084         int i;
1085         int sizeKeyContour,Z=(*InZ)[0];
1086         sizeKeyContour = _KeyContourZ.size();
1087         for (i=0; i<sizeKeyContour; i++)
1088         {
1089                 if (i>0)
1090                 {
1091                         if ( (_KeyContourZ[i-1]<Z) && (_KeyContourZ[i]>=Z) )
1092                         {
1093                                 idKeyContour            = i;
1094                                 idKeyContourSizes       = tmpIdKeyContSizes;
1095                                 okFind=true;
1096                                 i=sizeKeyContour;
1097                         }
1098                         if ( (_KeyContourZ[i-1] != _KeyContourZ[i]) )
1099                         {
1100                                 tmpIdKeyContSizes++;
1101                         }
1102                 } else {
1103                         if  (_KeyContourZ[0]>Z)
1104                         {
1105                                 idKeyContour            = 0;
1106                                 idKeyContourSizes       = 0;
1107                                 okFind                          = true;
1108                                 i                                       = sizeKeyContour;
1109                         }
1110                 } // if >O
1111         } // for
1112
1113         if (okFind==false)
1114         {
1115                 idKeyContour            = _KeyContourX.size();
1116                 idKeyContourSizes       = _KeyContourSizes.size();
1117                 okFind=true;
1118         }
1119
1120         _KeyContourSizes.insert( _KeyContourSizes.begin() + idKeyContourSizes , InX->size() );
1121         for(i=0; i<InX->size(); i++)
1122         {
1123                 _KeyContourX.insert( _KeyContourX.begin() + idKeyContour, (*InX)[i] );
1124                 _KeyContourY.insert( _KeyContourY.begin() + idKeyContour, (*InY)[i] );
1125                 _KeyContourZ.insert( _KeyContourZ.begin() + idKeyContour, (*InZ)[i] );
1126         }
1127
1128
1129 //EED
1130 //      _KeyContourSizes.push_back( InX->size() );
1131 //      for(i=0; i<InX->size(); i++)
1132 //      {
1133 //              _KeyContourX.push_back( (*InX)[i] );
1134 //              _KeyContourY.push_back( (*InY)[i] );
1135 //              _KeyContourZ.push_back( (*InZ)[i] );
1136 //      }
1137
1138 }
1139 //---------------------------------------------------------------------------------------------------
1140 vtkImageData* PropContour::CalculeSplinePropagation()
1141 {
1142         if(_KeyContourSizes.size() <= 0)
1143         {
1144                 printf("\n There would be at last 1 contour");
1145                 return NULL;
1146         }
1147         if(_KeyContourSizes.size() == 1)
1148         {
1149                 return NULL;
1150         }
1151         if(_KeyContourSizes.size() >= 2)
1152         {
1153                 return method_Spline(&_KeyContourX,&_KeyContourY,&_KeyContourZ,&_KeyContourSizes);
1154         }
1155 }
1156 //---------------------------------------------------------------------------------------------------
1157 void PropContour::GetKeyContours(std::vector<double>*KeyX, std::vector<double>*KeyY, std::vector<double>*KeyZ, std::vector<int>*KeyS)
1158 {
1159         int i;
1160         KeyX->clear();
1161         KeyY->clear();
1162         KeyZ->clear();
1163         KeyS->clear();
1164
1165         for(i=0; i<_KeyContourSizes.size(); i++)
1166         {
1167                 KeyS->push_back( _KeyContourSizes[i] );
1168         }
1169         for(i=0; i<_KeyContourX.size(); i++)
1170         {
1171                 KeyX->push_back( _KeyContourX[i] );
1172                 KeyY->push_back( _KeyContourY[i] );
1173                 KeyZ->push_back( _KeyContourZ[i] );
1174         }
1175 }
1176 //---------------------------------------------------------------------------------------------------
1177 void PropContour::GetPropagatedContours( std::vector<Vector>*planevec )
1178 {
1179         int i,j;
1180         planevec->clear();
1181         for(i=0; i<_planevector.size(); i++)
1182         {
1183                 Vector *temp = new Vector();
1184                 temp->set_plane( _planevector[i].get_plane() );
1185                 for(j=0; j<_planevector[i].getsize_x(); j++)
1186                 {
1187                         temp->set_x( _planevector[i].get_x(j) );
1188                         temp->set_y( _planevector[i].get_y(j) );
1189                         temp->set_z( _planevector[i].get_z(j) );
1190                 }
1191                 planevec->push_back(*temp);
1192                 delete temp;
1193         }
1194 }
1195 //---------------------------------------------------------------------------------------------------
1196 void PropContour::SetInterpNumber(int val)
1197 {
1198         _interpnumber = val;
1199 }
1200
1201 int  PropContour::FindIdWithZ(double z)
1202 {
1203         int result_ID=0;
1204         int k,size=_planevector.size();
1205         double dist,minDist=99999999;
1206         for ( k=0 ; k<size ; k++)
1207         {
1208                 dist=fabs( z - _planevector[k].get_z(0) );
1209                 if (dist<minDist)
1210                 {
1211                         minDist=dist;
1212                         result_ID = k;
1213                 }
1214         }// for i
1215         return result_ID;
1216 }
1217
1218 //---------------------------------------------------------------------------------------------------
1219 void PropContour::GetIdContour(int id, std::vector<double>*vecX, std::vector<double>*vecY, std::vector<double>*vecZ)
1220 {
1221         int i;
1222         vecX->clear();
1223         vecY->clear();
1224         vecZ->clear();
1225         int sizeplane = _planevector[id].getsize_x();
1226         double tempx;
1227         for(i=0; i<_planevector[id].getsize_x(); i++)
1228         {
1229                 vecX->push_back( _planevector[id].get_x(i) );
1230                 tempx = _planevector[id].get_x(i);
1231                 vecY->push_back( _planevector[id].get_y(i) );
1232                 vecZ->push_back( _planevector[id].get_z(i) );
1233         }
1234 }
1235 //---------------------------------------------------------------------------------------------------
1236 //---------------------------------------------------------------------------------------------------
1237 //---------------------------------------------------------------------------------------------------