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<(int)(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<(int)(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<(int)(_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<(int)(_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<(int)(_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); // JPRx
 730                 fscanf(fd," %s %d",firstline,&size);
 731                 for(int i=0; i<size; i++)
 732                 {
 733                         fscanf(fd,"%lf %lf %d",&x,&y,&z);
 734                         tempX.push_back(x);
 735                         tempY.push_back(y);
 736                         tempZ.push_back(z);
 737                 }
 738                 SetKeyContours(&tempX,&tempY,&tempZ);
 739                 tempX.clear();
 740                 tempY.clear();
 741                 tempZ.clear();
 742         }
 743 }
 744 //---------------------------------------------------------------------------------------------------------
 745 int PropContour::VectorDirection(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
 746 {
 747         int dir,i;
 748         double SumX = 0,SumY = 0;
 749         double ax,ay,bx,by,axb;
 750         int size = InX->size();
 751         for(i=0; i<size; i++)
 752         {
 753                 SumX = SumX + (*InX)[i];
 754                 SumY = SumY + (*InY)[i];
 755         }
 756         SumX = SumX/size;               //Mass Center: X coord
 757         SumY = SumY/size;               //Mass Center: Y coord
 758
 759         int positive = 0;
 760         int negative = 0;
 761         for(i=0; i<size; i++)
 762         {
 763                 ax = (*InX)[i]-SumX;
 764                 ay = (*InY)[i]-SumY;
 765                 bx = (*InX)[(i+1)%size]-SumX;
 766                 by = (*InY)[(i+1)%size]-SumY;
 767                 axb = (ax*by) - (bx*ay);
 768                 if(axb > 0)
 769                 {
 770                         positive++;
 771                 }
 772                 if(axb < 0)
 773                 {
 774                         negative++;
 775                 }
 776         }
 777         if(positive >= negative)
 778         {
 779                 dir = 1;
 780         }
 781         else
 782         {
 783                 dir = -1;
 784         }
 785
 786         return dir;
 787 }
 788 //----------------------------------------------------------------------------------------------------
 789 void PropContour::VectorOrder(int dir, int posini, std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
 790 {
 791         int size = InX->size();
 792         int i;
 793         std::vector<double> tempX;
 794         std::vector<double> tempY;
 795         std::vector<double> tempZ;
 796
 797         tempX.clear();
 798         tempY.clear();
 799         tempZ.clear();
 800
 801         for(i=0; i<size; i++)
 802         {
 803                 if(dir == 1)
 804                 {
 805                         tempX.push_back((*InX)[posini]);
 806                         tempY.push_back((*InY)[posini]);
 807                         tempZ.push_back((*InZ)[posini]);
 808                         posini++;
 809                         if(posini == size)
 810                         {
 811                                 posini = 0;
 812                         }
 813                 }
 814                 if(dir == -1)
 815                 {
 816                         tempX.push_back((*InX)[posini]);
 817                         tempY.push_back((*InY)[posini]);
 818                         tempZ.push_back((*InZ)[posini]);
 819                         posini--;
 820                         if(posini < 0)
 821                         {
 822                                 posini = size-1;
 823                         }
 824                 }
 825         }
 826         InX->clear();
 827         InY->clear();
 828         InZ->clear();
 829         for(i=0; i<size; i++)
 830         {
 831                 InX->push_back(tempX[i]);
 832                 InY->push_back(tempY[i]);
 833                 InZ->push_back(tempZ[i]);
 834         }
 835
 836 }
 837 //----------------------------------------------------------------------------------------------------
 838 void PropContour::PreparePointsForSpline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ,
 839                                                                                  std::vector<int>*Sizes)
 840 {
 841         int sizeS = Sizes->size();
 842         //int sizeV = InX->size(); // JPRx
 843         int i,j,mem,posinic,dir,cont;
 844         double leX;
 845
 846         std::vector<double> tempX;
 847         std::vector<double> tempY;
 848         std::vector<double> tempZ;
 849         std::vector<double> lstX;
 850         std::vector<double> lstY;
 851         std::vector<double> lstZ;
 852
 853         lstX.clear();
 854         lstY.clear();
 855         lstZ.clear();
 856
 857         mem = 0;
 858         cont = 0;
 859         for(i=0; i<sizeS; i++)
 860         {
 861                 leX=1000;
 862                 tempX.clear();
 863                 tempY.clear();
 864                 tempZ.clear();
 865                 for(j=0; j<(*Sizes)[i]; j++)
 866                 {
 867                         tempX.push_back((*InX)[j+mem]);
 868                         tempY.push_back((*InY)[j+mem]);
 869                         tempZ.push_back((*InZ)[j+mem]);
 870                         if( (*InX)[j] < leX )
 871                         {
 872                                 posinic = j;
 873                                 leX = (*InX)[j];
 874                         }
 875                 }
 876                 mem = mem + (*Sizes)[i];
 877                 dir = VectorDirection(&tempX,&tempY,&tempZ);
 878                 VectorOrder(dir,posinic,&tempX,&tempY,&tempZ);
 879
 880                 for(j=0; j<(*Sizes)[i]; j++)
 881                 {
 882                         lstX.push_back(tempX[j]);
 883                         lstY.push_back(tempY[j]);
 884                         lstZ.push_back((*InZ)[j+cont]);
 885                 }
 886                 cont = cont + (*Sizes)[i];
 887         }
 888
 889 //Fill the Finally lst in X,Y,Z ---------------
 890         int sizetemp = lstX.size();
 891         //printf("\nJSTG-PropContour::PreparePointsForSpline");
 892         InX->clear();
 893         InY->clear();
 894         InZ->clear();
 895         for(i=0; i<sizetemp; i++)
 896         {
 897                 //printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
 898                 InX->push_back(lstX[i]);
 899                 InY->push_back(lstY[i]);
 900                 InZ->push_back(lstZ[i]);
 901         }
 902
 903
 904 }
 905 //----------------------------------------------------------------------------------------------------
 906 vtkImageData* PropContour::method_Spline(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ, std::vector<int>*Sizes)
 907 {
 908         //long interval = wxGetElapsedTime(TRUE); // JPRx
 909
 910         int i,j,k,sizeX,sizeS,sizeInS;
 911         int numspline;
 912         double x,y,z;
 913 //EED   double spc[3]={1,1,1};
 914         std::vector<double> lstX;
 915         std::vector<double> lstY;
 916         std::vector<double> lstZ;
 917         std::vector<double> tempX;
 918         std::vector<double> tempY;
 919         std::vector<double> tempZ;
 920         std::vector<double> tempS;
 921         _mContourModel = new manualContourModel();
 922         _mContourModel->SetNumberOfPointsSpline(_interpnumber);
 923
 924         imagedataValue=NULL;
 925 //EED
 926 //      _dimImage[0] = 200;                     // X axis
 927 //      _dimImage[1] = 200;                     // Y axis
 928 //      _dimImage[2] = 200;                     // Z axis
 929 //
 930 //      unsigned short *pValue;
 931 //      imagedataValue = vtkImageData::New();
 932 //      imagedataValue->SetScalarTypeToUnsignedShort();
 933 //      imagedataValue->SetSpacing(spc);
 934 //      imagedataValue->SetDimensions(_dimImage);
 935 //      imagedataValue->SetExtent(0,_dimImage[0]-1,0,_dimImage[1]-1,0,_dimImage[2]-1);
 936 //      imagedataValue->AllocateScalars();
 937 //      imagedataValue->Update();
 938
 939 //      lstX.clear();
 940 //      lstY.clear();
 941 //      lstZ.clear();
 942
 943         sizeX = InX->size();
 944         for(i=0; i<sizeX; i++)
 945         {
 946                 lstX.push_back((*InX)[i]);
 947                 lstY.push_back((*InY)[i]);
 948                 lstZ.push_back((*InZ)[i]);
 949         }
 950
 951         PreparePointsForSpline(&lstX,&lstY,&lstZ,Sizes);
 952         /*int sizetemp = lstX.size();
 953         printf("\nJSTG-PropContour::method_Spline");
 954         for(i=0; i<sizetemp; i++)
 955         {
 956                 printf("\nlst: Z = %f, X = %f, Y = %f, ",lstZ[i],lstX[i],lstY[i]);
 957         }*/
 958
 959         sizeS = Sizes->size();
 960         int cont = 0;
 961         for(i=0; i<sizeS; i++)
 962         {
 963                 _mContourModel->DeleteAllPoints();
 964                 _mContourModel->SetCloseContour(true);
 965                 sizeInS = (*Sizes)[i];
 966                 for(j=0; j<sizeInS; j++)
 967                 {
 968                         _mContourModel->AddPoint(lstX[cont],lstY[cont],lstZ[cont]);
 969                         cont ++;
 970                 }
 971                 _mContourModel->UpdateSpline();
 972                 numspline = _mContourModel->GetNumberOfPointsSpline();
 973                 for(k=0; k<numspline; k++)
 974                 {
 975                         _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
 976 //EED                   pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
 977 //                      *pValue = 128;
 978                         tempX.push_back(x);
 979                         tempY.push_back(y);
 980                         tempZ.push_back(z);
 981                 }
 982         }
 983
 984         int tam = numspline;
 985         std::vector<Vector>::iterator it;
 986
 987         ResetPlaneVector();
 988
 989         for(i=0; i<numspline; i++)
 990         {
 991                 Vector *vec = new Vector();
 992                 _planevector.push_back(*vec);
 993         }// for i
 994
 995         for(j=0; j<tam; j++)
 996         {
 997                 _mContourModel->DeleteAllPoints();
 998                 _mContourModel->SetCloseContour(false);
 999                 cont = 0;
1000                 for(i=0; i<sizeS; i++)
1001                 {
1002                         //double hh=tempZ[cont+j]; // JPRx
1003                         _mContourModel->AddPoint(tempX[cont+j],tempY[cont+j],tempZ[cont+j]);
1004                         cont = cont + tam;
1005                 } // for i
1006                 _mContourModel->UpdateSpline();
1007                 numspline = _mContourModel->GetNumberOfPointsSpline();
1008                 for(k=0; k<numspline; k++)
1009                 {
1010 //EED002
1011                         _mContourModel->GetSpline_i_Point(k,&x,&y,&z);
1012                         _planevector[k].set_x(x);
1013                         _planevector[k].set_y(y);
1014                         _planevector[k].set_z(z);
1015                         _planevector[k].set_plane(k);
1016 //EED                   pValue = (unsigned short *)imagedataValue->GetScalarPointer(x,y,z);
1017 //EED                   *pValue = 128;
1018                 } // for k
1019         } // for j
1020
1021 /*
1022         int tempsize = _planevector.size();
1023         int sizexx;
1024         for(i=0; i<tempsize; i++)
1025         {
1026                 sizexx = _planevector[i].getsize_x();
1027         }
1028 */
1029 /*
1030         interval = wxGetElapsedTime(FALSE);
1031         int zplanes = abs(tempZ[0]-tempZ[tempZ.size()-1]);
1032         long interPlane = interval/zplanes;
1033         int pointsSize = 0;
1034         for(i=0; i<sizeS; i++)
1035         {
1036                 pointsSize = pointsSize + (*Sizes)[i];
1037         }
1038         for(i=0; i<pointsSize; i++)
1039         {
1040                  pValue = (unsigned short *)imagedataValue->GetScalarPointer( (int)(*InX)[i] , (int)(*InY)[i] ,(int)(*InZ)[i]);
1041                 *pValue=255;
1042         }
1043
1044         printf("\n\n JSTG - PropContour::method_Spline ------------------------");
1045         printf("\n TIME FOR: IMAGE 3D WITHOUT THE CONTROL POINTS.... %d (ms)",interval);
1046         //printf("\n TIME FOR: IMAGE 3D WITH THE CONTROL POINTS....... %lld (ms)",intervalPC);
1047         printf("\n TIME AVERAGE FOR: EVERY PLANE.................... %d (ms)",interPlane);
1048         printf("\n NUMBER OF PLANES................................. %d",zplanes);
1049         printf("\n TOTAL NUMBER OF CONTROL POINTS................... %d",pointsSize);
1050         printf("\n ------------------------------------------------------------");
1051 */
1052         return imagedataValue;
1053 }
1054
1055 //---------------------------------------------------------------------------------------------------
1056 void PropContour::ResetPlaneVector()
1057 {
1058 //      Vector *vec;
1059 //      int ii,iiSize = _planevector.size();
1060 //      for (ii=0 ; ii<iiSize ; ii++)
1061 //      {
1062 //              vec = &(_planevector[ii]);
1063 //              delete vec;
1064 //      }
1065         _planevector.clear();
1066 }
1067
1068
1069 //---------------------------------------------------------------------------------------------------
1070 void PropContour::ResetKeyContours()
1071 {
1072         _KeyContourSizes.clear();
1073         _KeyContourX.clear();
1074         _KeyContourY.clear();
1075         _KeyContourZ.clear();
1076 }
1077 //---------------------------------------------------------------------------------------------------
1078 void PropContour::SetKeyContours(std::vector<double>*InX, std::vector<double>*InY, std::vector<double>*InZ)
1079 {
1080
1081         int idKeyContour                = 0;
1082         int idKeyContourSizes   = 0;
1083         int tmpIdKeyContSizes   = 0;
1084         bool okFind                             = false;
1085         int i;
1086         int sizeKeyContour,Z=(int)(*InZ)[0];
1087         sizeKeyContour = _KeyContourZ.size();
1088         for (i=0; i<sizeKeyContour; i++)
1089         {
1090                 if (i>0)
1091                 {
1092                         if ( (_KeyContourZ[i-1]<Z) && (_KeyContourZ[i]>=Z) )
1093                         {
1094                                 idKeyContour            = i;
1095                                 idKeyContourSizes       = tmpIdKeyContSizes;
1096                                 okFind=true;
1097                                 i=sizeKeyContour;
1098                         }
1099                         if ( (i<sizeKeyContour) && (_KeyContourZ[i-1] != _KeyContourZ[i]) )
1100                         {
1101                                 tmpIdKeyContSizes++;
1102                         }
1103                 } else {
1104                         if  (_KeyContourZ[0]>Z)
1105                         {
1106                                 idKeyContour            = 0;
1107                                 idKeyContourSizes       = 0;
1108                                 okFind                          = true;
1109                                 i                                       = sizeKeyContour;
1110                         }
1111                 } // if >O
1112         } // for
1113
1114         if (okFind==false)
1115         {
1116                 idKeyContour            = _KeyContourX.size();
1117                 idKeyContourSizes       = _KeyContourSizes.size();
1118                 okFind=true;
1119         }
1120
1121         _KeyContourSizes.insert( _KeyContourSizes.begin() + idKeyContourSizes , InX->size() );
1122         for(i=0; i<(int)(InX->size()); i++)
1123         {
1124                 _KeyContourX.insert( _KeyContourX.begin() + idKeyContour, (*InX)[i] );
1125                 _KeyContourY.insert( _KeyContourY.begin() + idKeyContour, (*InY)[i] );
1126                 _KeyContourZ.insert( _KeyContourZ.begin() + idKeyContour, (*InZ)[i] );
1127         }
1128
1129
1130 //EED
1131 //      _KeyContourSizes.push_back( InX->size() );
1132 //      for(i=0; i<InX->size(); i++)
1133 //      {
1134 //              _KeyContourX.push_back( (*InX)[i] );
1135 //              _KeyContourY.push_back( (*InY)[i] );
1136 //              _KeyContourZ.push_back( (*InZ)[i] );
1137 //      }
1138
1139 }
1140 //---------------------------------------------------------------------------------------------------
1141 vtkImageData* PropContour::CalculeSplinePropagation()
1142 {
1143         if(_KeyContourSizes.size() <= 0)
1144         {
1145                 printf("\n There would be at last 1 contour");
1146                 return NULL;
1147         }
1148         if(_KeyContourSizes.size() == 1)
1149         {
1150                 return NULL;
1151         }
1152         if(_KeyContourSizes.size() >= 2)
1153         {
1154                 return method_Spline(&_KeyContourX,&_KeyContourY,&_KeyContourZ,&_KeyContourSizes);
1155         }
1156         return NULL;
1157 }
1158 //---------------------------------------------------------------------------------------------------
1159 void PropContour::GetKeyContours(std::vector<double>*KeyX, std::vector<double>*KeyY, std::vector<double>*KeyZ, std::vector<int>*KeyS)
1160 {
1161         int i;
1162         KeyX->clear();
1163         KeyY->clear();
1164         KeyZ->clear();
1165         KeyS->clear();
1166
1167         for(i=0; i<(int)(_KeyContourSizes.size()); i++)
1168         {
1169                 KeyS->push_back( _KeyContourSizes[i] );
1170         }
1171         for(i=0; i<(int)(_KeyContourX.size()); i++)
1172         {
1173                 KeyX->push_back( _KeyContourX[i] );
1174                 KeyY->push_back( _KeyContourY[i] );
1175                 KeyZ->push_back( _KeyContourZ[i] );
1176         }
1177 }
1178 //---------------------------------------------------------------------------------------------------
1179 void PropContour::GetPropagatedContours( std::vector<Vector>*planevec )
1180 {
1181         int i,j;
1182         planevec->clear();
1183         for(i=0; i<(int)(_planevector.size()); i++)
1184         {
1185                 Vector *temp = new Vector();
1186                 temp->set_plane( _planevector[i].get_plane() );
1187                 for(j=0; j<_planevector[i].getsize_x(); j++)
1188                 {
1189                         temp->set_x( _planevector[i].get_x(j) );
1190                         temp->set_y( _planevector[i].get_y(j) );
1191                         temp->set_z( _planevector[i].get_z(j) );
1192                 }
1193                 planevec->push_back(*temp);
1194                 delete temp;
1195         }
1196 }
1197 //---------------------------------------------------------------------------------------------------
1198 void PropContour::SetInterpNumber(int val)
1199 {
1200         _interpnumber = val;
1201 }
1202
1203 int  PropContour::FindIdWithZ(double z)
1204 {
1205         int result_ID=0;
1206         int k,size=_planevector.size();
1207         double dist,minDist=99999999;
1208         for ( k=0 ; k<size ; k++)
1209         {
1210                 dist=fabs( z - _planevector[k].get_z(0) );
1211                 if (dist<minDist)
1212                 {
1213                         minDist=dist;
1214                         result_ID = k;
1215                 }
1216         }// for i
1217         return result_ID;
1218 }
1219
1220 //---------------------------------------------------------------------------------------------------
1221 void PropContour::GetIdContour(int id, std::vector<double>*vecX, std::vector<double>*vecY, std::vector<double>*vecZ)
1222 {
1223         int i;
1224         vecX->clear();
1225         vecY->clear();
1226         vecZ->clear();
1227         //int sizeplane = _planevector[id].getsize_x();
1228         double tempx;
1229         for(i=0; i<_planevector[id].getsize_x(); i++)
1230         {
1231                 vecX->push_back( _planevector[id].get_x(i) );
1232                 tempx = _planevector[id].get_x(i);
1233                 vecY->push_back( _planevector[id].get_y(i) );
1234                 vecZ->push_back( _planevector[id].get_z(i) );
1235         }
1236 }
1237 //---------------------------------------------------------------------------------------------------
1238 //---------------------------------------------------------------------------------------------------
1239 //---------------------------------------------------------------------------------------------------