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

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