Lib/SpatioTemporalMeanShift/itkSTMS_BlurringSTMS_Spine.txx

   1 /*
   2  #
   3  #  File        : itkSTMS_BlurringSTMS_Spine.txx
   4  #                ( C++ header file - STMS )
   5  #
   6  #  Description : STMS lib that implements the STMS filter and clustering.
   7  #                This file is a part of the STMS Library project.
   8  #                ( https://www.creatis.insa-lyon.fr/site7/fr/realisations )
   9  #
  10  #  [1] S. Mure, Grenier, T., Meier, S., Guttmann, R. G., et Benoit-Cattin, H.,
  11  #       « Unsupervised spatio-temporal filtering of image sequences. A mean-shift specification »,
  12  #       Pattern Recognition Letters, vol. 68, Part 1, p. 48 - 55, 2015.
  13  #
  14  #  Copyright   : Thomas GRENIER - Simon MURE
  15  #                ( https://www.creatis.insa-lyon.fr/~grenier/ )
  16  #
  17  #  License     : CeCILL C
  18  #                ( http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.txt )
  19  #
  20  #  This software is governed by the CeCILL  license under French law and
  21  #  abiding by the rules of distribution of free software.  You can  use,
  22  #  modify and/ or redistribute the software under the terms of the CeCILL
  23  #  license as circulated by CEA, CNRS and INRIA at the following URL
  24  #  "http://www.cecill.info".
  25  #
  26  #  As a counterpart to the access to the source code and  rights to copy,
  27  #  modify and redistribute granted by the license, users are provided only
  28  #  with a limited warranty  and the software's author,  the holder of the
  29  #  economic rights,  and the successive licensors  have only  limited
  30  #  liability.
  31  #
  32  #  In this respect, the user's attention is drawn to the risks associated
  33  #  with loading,  using,  modifying and/or developing or reproducing the
  34  #  software by the user in light of its specific status of free software,
  35  #  that may mean  that it is complicated to manipulate,  and  that  also
  36  #  therefore means  that it is reserved for developers  and  experienced
  37  #  professionals having in-depth computer knowledge. Users are therefore
  38  #  encouraged to load and test the software's suitability as regards their
  39  #  requirements in conditions enabling the security of their systems and/or
  40  #  data to be ensured and,  more generally, to use and operate it in the
  41  #  same conditions as regards security.
  42  #
  43  #  The fact that you are presently reading this means that you have had
  44  #  knowledge of the CeCILL license and that you accept its terms.
  45  #
  46 */
  47 /* Please don't forget to cite our work :
  48   @article {MURE-15a,
  49     title = {Unsupervised spatio-temporal filtering of image sequences. A mean-shift specification},
  50     journal = {Pattern Recognition Letters},
  51     volume = {68, Part 1},
  52     year = {2015},
  53     pages = {48 - 55},
  54     issn = {0167-8655},
  55     doi = {http://dx.doi.org/10.1016/j.patrec.2015.07.021},
  56     url = {http://www.sciencedirect.com/science/article/pii/S0167865515002305},
  57     author = {S. Mure and T Grenier and Meier, S. and Guttmann, R.G. and H. Benoit-Cattin}
  58 }
  59 */
  60 #ifndef itkSTMS_BlurringSTMS_TXX
  61 #define itkSTMS_BlurringSTMS_TXX
  62
  63 #include <random>
  64 #include <algorithm>
  65 #include <iterator>
  66 #include "itkSTMS_BlurringSTMS_Spine.h"
  67
  68
  69 double gettime()
  70 {
  71     struct timespec timestamp;
  72
  73     clock_gettime(CLOCK_REALTIME, &timestamp);
  74     return timestamp.tv_sec * 1000.0 + timestamp.tv_nsec * 1.0e-6;
  75 }
  76
  77
  78 namespace itkSTMS_Spine
  79 {
  80 template < class IndexType, class SpatialType, class PixelType, class ImageType >
  81 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
  82 ::itkSTMS_BlurringSTMS( IndexSampleSetType* idx, IndexSampleSetType* cla, IndexSampleSetType* wei,
  83                         SpatialSampleSetType* sp, RangeSampleSetType* ra, ParametersType* params )
  84 {
  85     this->indexSet   = idx;
  86     this->classSet   = cla;
  87     this->weightsSet = wei;
  88     this->spatialSet = sp;
  89     this->rangeSet   = ra;
  90
  91     this->classSetMemory   = new IndexSampleSetType  (  *this->classSet  );
  92     this->spatialSetMemory = new SpatialSampleSetType( *this->spatialSet );
  93
  94     this->stmsParams     = params;
  95
  96     this->mergeFactor = 1000;
  97 }
  98
  99
 100 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 101 void
 102 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 103 ::GenerateData()
 104 {
 105     double dtime;
 106
 107     if( !stmsParams->merge )
 108     {
 109         std::cout<< std::endl<<"No merge Filtering!"<<std::endl<<"numSamples: "<< indexSet->size() <<std::endl<< std::endl;
 110         NoMergeSTMSFiltering();
 111
 112         dtime=gettime();
 113         ClassificationNoMergeSTMSFiltering();
 114         dtime = gettime()-dtime;
 115         std::cout<< std::endl<<"Classif: " << dtime/1000 << " s" <<std::endl;
 116         FinalMerging();
 117     }
 118     else
 119     {
 120         std::cout<<std::endl<<"Merge Filtering!"<<std::endl<<"numSamples: "<< indexSet->size() <<std::endl<< std::endl;
 121         MergeSTMSFiltering();
 122         FinalMerging();
 123         //        ClassificationNoMergeSTMSFiltering();
 124     }
 125 }
 126
 127 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 128 void
 129 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 130 ::ClassificationNoMergeSTMSFiltering()
 131 {
 132
 133     IndexSampleSetType*   newClassSet
 134             = new IndexSampleSetType();
 135     newClassSet->reserve  ( indexSet->size() );
 136
 137     IndexSampleSetType*   newWeightsSet
 138             = new IndexSampleSetType();
 139     newWeightsSet->reserve( indexSet->size() );
 140
 141     IndexSampleSetType*   newIndexSet
 142             = new IndexSampleSetType();
 143     newIndexSet->reserve  ( indexSet->size() );
 144
 145     SpatialSampleSetType* newSpatialSet
 146             = new SpatialSampleSetType();
 147
 148     newSpatialSet->reserve( indexSet->size() );
 149
 150     RangeSampleSetType*   newRangeSet
 151             = new RangeSampleSetType();
 152     newRangeSet->reserve  ( indexSet->size() );
 153
 154
 155     newClassSet->push_back  ( 1 );
 156     newWeightsSet->push_back( 1 );
 157     newIndexSet->push_back  ( indexSet->at(0)   );
 158     newSpatialSet->push_back( spatialSet->at(0) );
 159     newRangeSet->push_back  ( rangeSet->at(0)   );
 160     classSetMemory->at( indexSet->at(0) ) = newClassSet->at(0);
 161
 162
 163     unsigned int k;
 164     bool newC;
 165     float spDist;
 166     bool  raNorm;
 167
 168     for(unsigned int i=1 ; i<indexSet->size() ; ++i)
 169     {
 170         newC = true;
 171         k=0;
 172
 173         while( (k<newClassSet->size()) && newC )
 174         {
 175             VectorDistance( spDist, newSpatialSet->at(k), spatialSet->at(i) );
 176             if( spDist<=1 )
 177             {
 178                 InfiniteNorm( raNorm, newRangeSet->at(k), rangeSet->at(i) );
 179                 if( raNorm )
 180                 {
 181                     VectorWeightedMean( newSpatialSet->at(k), newWeightsSet->at(k), spatialSet->at(i), weightsSet->at(i) );
 182                     VectorWeightedMean( newRangeSet->at(k)  , newWeightsSet->at(k), rangeSet->at(i)  , weightsSet->at(i) );
 183
 184                     classSetMemory->at( indexSet->at(i) ) = newClassSet->at(k);
 185                     newC = false;
 186                 }
 187                 else
 188                     ++k;
 189             }
 190             else
 191                 ++k;
 192         }
 193
 194         if( newC )
 195         {
 196             newClassSet->push_back( (IndexType)newClassSet->back()+1 );
 197             newWeightsSet->push_back( 1 );
 198             newIndexSet->push_back  ( indexSet->at(i)   );
 199             newSpatialSet->push_back( spatialSet->at(i) );
 200             newRangeSet->push_back  ( rangeSet->at(i)   );
 201             classSetMemory->at( indexSet->at(i) ) = newClassSet->back();
 202         }
 203     }
 204
 205     classSet->swap  ( *newClassSet   );
 206     weightsSet->swap( *newWeightsSet );
 207     indexSet->swap  ( *newIndexSet   );
 208     spatialSet->swap( *newSpatialSet );
 209     rangeSet->swap  ( *newRangeSet   );
 210
 211     delete newClassSet;
 212     delete newWeightsSet;
 213     delete newIndexSet;
 214     delete newSpatialSet;
 215     delete newRangeSet;
 216 }
 217
 218
 219 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 220 void
 221 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 222 ::NoMergeSTMSFiltering()
 223 {
 224
 225     SpatialSampleSetType* newSpatialSet
 226             = new SpatialSampleSetType( indexSet->size(), SpatialVectorType(stmsParams->dim, 0) );
 227     RangeSampleSetType*   newRangeSet
 228             = new RangeSampleSetType( indexSet->size(), RangeVectorType(stmsParams->numTimePoints, 0) );
 229
 230     float globalEvolution = INFINITY;
 231     float epsilon = (stmsParams->epsilon)*(stmsParams->epsilon);
 232     float spDist, raDist;
 233     bool  raNorm;
 234     unsigned int iter = 0;
 235
 236     double dtime;
 237     while( (globalEvolution>epsilon) & (iter<stmsParams->maxIt) )
 238     {
 239         dtime=gettime();
 240         globalEvolution = 0.0;
 241
 242         #pragma omp parallel for reduction(+:globalEvolution) private(spDist, raDist, raNorm)
 243         for(unsigned int i=0 ; i<indexSet->size() ; ++i)
 244         {
 245             unsigned int count = 0;
 246             SpatialVectorType tmpSpatial( stmsParams->dim, 0 );
 247             RangeVectorType   tmpRange  ( stmsParams->numTimePoints, 0 );
 248
 249             for(unsigned int j=0 ; j<indexSet->size() ; ++j)
 250             {
 251                 VectorDistance( spDist, spatialSet->at(i), spatialSet->at(j));
 252
 253                 if( spDist<=1 )
 254                 {
 255                     InfiniteNorm(raNorm, rangeSet->at(i), rangeSet->at(j));
 256
 257                     if( raNorm )
 258                     {
 259                         VectorAcc( tmpSpatial, spatialSet->at(j) );
 260                         VectorAcc( tmpRange  , rangeSet->at(j)   );
 261
 262                         ++count;
 263                     }
 264                 }
 265             }
 266
 267             if( count>1 )
 268             {
 269                 VectorDiv( tmpSpatial, count );
 270                 VectorDiv( tmpRange  , count );
 271             }
 272
 273             VectorDistance(spDist, spatialSet->at(i), tmpSpatial);
 274             VectorDistance(raDist, rangeSet->at(i), tmpRange);
 275
 276             globalEvolution += (spDist+raDist);
 277
 278             #pragma omp critical
 279             {
 280                 newSpatialSet->at(i) = tmpSpatial;
 281                 newRangeSet->at(i)   = tmpRange;
 282             }
 283         }
 284
 285         spatialSet->swap( *newSpatialSet );
 286         rangeSet->swap  (  *newRangeSet  );
 287
 288         ++iter;
 289         dtime = gettime()-dtime;
 290         std::cout<< "Iter: " << iter <<"   "<< dtime/1000 << " s" <<"    GE: "<< globalEvolution <<std::endl;
 291     }
 292
 293     delete newSpatialSet;
 294     delete newRangeSet;
 295 }
 296
 297
 298
 299
 300 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 301 void
 302 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 303 ::MergeSTMSFiltering()
 304 {
 305
 306     float globalEvolution = INFINITY;
 307     float epsilon = (stmsParams->epsilon)*(stmsParams->epsilon);
 308     float spDist, raDist;
 309     bool  raNorm, mergeNorm;
 310     unsigned int iter = 0;
 311     unsigned int merging;
 312
 313     double dtime;
 314     while( (globalEvolution>epsilon) & (iter<stmsParams->maxIt) )
 315     {
 316         SpatialSampleSetType* newSpatialSet
 317                 = new SpatialSampleSetType( *spatialSet );
 318         RangeSampleSetType*   newRangeSet
 319                 = new RangeSampleSetType  (  *rangeSet  );
 320
 321         dtime=gettime();
 322         globalEvolution = 0.0;
 323         merging = 0;
 324
 325         #pragma omp parallel for reduction(+:globalEvolution, merging) private(spDist, raDist, raNorm, mergeNorm)
 326         for(unsigned int i=0 ; i<indexSet->size() ; ++i)
 327         {
 328             unsigned int weight = 0;
 329             bool merge = false;
 330             SpatialVectorType tmpSpatial( stmsParams->dim, 0 );
 331             RangeVectorType   tmpRange  ( stmsParams->numTimePoints, 0 );
 332
 333             for(unsigned int j=0 ; j<indexSet->size() ; ++j)
 334             {
 335                 VectorDistance( spDist, spatialSet->at(i), spatialSet->at(j));
 336
 337                 if( spDist<=1 )
 338                 {
 339                     InfiniteNorm(raNorm, rangeSet->at(i), rangeSet->at(j));
 340
 341                     if( raNorm )
 342                     {
 343                         VectorWeightedAcc(tmpSpatial, spatialSet->at(j), weightsSet->at(j));
 344                         VectorWeightedAcc(tmpRange  , rangeSet->at(j)  , weightsSet->at(j));
 345
 346                         weight += weightsSet->at(j);
 347
 348                         MergeInfiniteNorm(mergeNorm, rangeSet->at(i), rangeSet->at(j));
 349                         if( (spDist <= (1/mergeFactor)) && mergeNorm )
 350                             merge = true;
 351                     }
 352                 }
 353             }
 354
 355             if(weight > weightsSet->at(i))
 356             {
 357                 if( merge )
 358                     merging += 1;
 359                 else
 360                     merge = 0;
 361             }
 362             else
 363                 merge = 0;
 364
 365
 366             VectorDiv( tmpSpatial, weight );
 367             VectorDiv( tmpRange  , weight );
 368
 369             #pragma omp critical
 370             {
 371                 newSpatialSet->at(i) = tmpSpatial;
 372                 newRangeSet->at(i)   = tmpRange;
 373             }
 374
 375             VectorDistance(spDist, spatialSet->at(i), tmpSpatial);
 376             VectorDistance(raDist, rangeSet->at(i), tmpRange);
 377
 378             globalEvolution += (spDist+raDist);
 379         }
 380
 381         spatialSet->swap( *newSpatialSet );
 382         rangeSet->swap  (  *newRangeSet  );
 383
 384         delete newSpatialSet;
 385         delete newRangeSet;
 386
 387         if(merging > 0)
 388             MergeSamples();
 389
 390         ++iter;
 391         dtime = gettime()-dtime;
 392         std::cout<< "Iter: " << iter <<"   "<< dtime/1000 << " s" <<"    GE: "<< globalEvolution << "      numSamples: "<< classSet->size() <<std::endl;
 393     }
 394 }
 395
 396
 397 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 398 void
 399 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 400 ::MergeSamples()
 401 {
 402
 403     IndexSampleSetType*   newClassSet
 404             = new IndexSampleSetType();
 405     newClassSet->reserve  ( indexSet->size() );
 406
 407     IndexSampleSetType*   newWeightsSet
 408             = new IndexSampleSetType();
 409     newWeightsSet->reserve( indexSet->size() );
 410
 411     IndexSampleSetType*   newIndexSet
 412             = new IndexSampleSetType();
 413     newIndexSet->reserve  ( indexSet->size() );
 414
 415     SpatialSampleSetType* newSpatialSet
 416             = new SpatialSampleSetType();
 417
 418     newSpatialSet->reserve( indexSet->size() );
 419
 420     RangeSampleSetType*   newRangeSet
 421             = new RangeSampleSetType();
 422     newRangeSet->reserve  ( indexSet->size() );
 423
 424     IndexSampleSetType* newClassSetMemory
 425             = new IndexSampleSetType( *classSetMemory );
 426
 427     IndexSampleSetType*   indexes
 428             = new IndexSampleSetType( *classSet );
 429
 430     std::random_device rd;
 431     std::mt19937 g(rd());
 432
 433     std::shuffle(indexes->begin(), indexes->end(), g);
 434
 435
 436     newClassSet->push_back  ( 1 );
 437     newWeightsSet->push_back( weightsSet->at( indexes->at(0)-1 ) );
 438     newIndexSet->push_back  ( indexSet->at  ( indexes->at(0)-1 ) );
 439     newSpatialSet->push_back( spatialSet->at( indexes->at(0)-1 ) );
 440     newRangeSet->push_back  ( rangeSet->at  ( indexes->at(0)-1 ) );
 441
 442     for(unsigned int l=0 ; l< classSetMemory->size() ; ++l)
 443     {
 444         if(classSetMemory->at(l) == classSet->at(indexes->at(0)-1))
 445             newClassSetMemory->at(l) = 1;
 446     }
 447
 448     unsigned int k;
 449     bool newC;
 450     float spDist;
 451     bool  raNorm;
 452     unsigned int i;
 453
 454     for(unsigned int m=1 ; m<indexes->size() ; ++m)
 455     {
 456         i = indexes->at(m)-1;
 457
 458         newC = true;
 459         k=0;
 460
 461         while( (k<newClassSet->size()) && newC )
 462         {
 463             VectorDistance( spDist, newSpatialSet->at(k), spatialSet->at(i) );
 464             MergeInfiniteNorm( raNorm, newRangeSet->at(k), rangeSet->at(i) );
 465             if( (spDist<=1/mergeFactor)&&raNorm )
 466             {
 467                 VectorWeightedAcc( newSpatialSet->at(k), newWeightsSet->at(k), spatialSet->at(i), weightsSet->at(i) );
 468                 VectorWeightedAcc( newRangeSet->at(k)  , newWeightsSet->at(k), rangeSet->at(i)  , weightsSet->at(i) );
 469
 470                 newWeightsSet->at(k) += weightsSet->at(i);
 471
 472                 VectorDiv(newSpatialSet->at(k), newWeightsSet->at(k));
 473                 VectorDiv(newRangeSet->at(k), newWeightsSet->at(k));
 474
 475                 for(unsigned int l=0 ; l<classSetMemory->size() ; ++l)
 476                 {
 477                     if(classSetMemory->at(l) == classSet->at(i))
 478                         newClassSetMemory->at(l) = newClassSet->at(k);
 479                 }
 480
 481                 newC = false;
 482             }
 483             else
 484                 ++k;
 485         }
 486
 487         if( newC )
 488         {
 489             newClassSet->push_back( newClassSet->back()+1 );
 490             newWeightsSet->push_back( weightsSet->at(i) );
 491             newIndexSet->push_back  ( indexSet->at(i)   );
 492             newSpatialSet->push_back( spatialSet->at(i) );
 493             newRangeSet->push_back  ( rangeSet->at(i)   );
 494
 495             for(unsigned int l=0 ; l< classSetMemory->size() ; ++l)
 496             {
 497                 if(classSetMemory->at(l) == classSet->at(i))
 498                     newClassSetMemory->at(l) = newClassSet->back();
 499             }
 500         }
 501     }
 502
 503     classSet->swap  ( *newClassSet   );
 504     weightsSet->swap( *newWeightsSet );
 505     indexSet->swap  ( *newIndexSet   );
 506     spatialSet->swap( *newSpatialSet );
 507     rangeSet->swap  ( *newRangeSet   );
 508     classSetMemory->swap  ( *newClassSetMemory );
 509
 510     delete newClassSet;
 511     delete newWeightsSet;
 512     delete newIndexSet;
 513     delete newSpatialSet;
 514     delete newRangeSet;
 515     delete newClassSetMemory;
 516     delete indexes;
 517 }
 518
 519
 520 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 521 void
 522 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 523 ::FinalMerging()
 524 {
 525     unsigned int size = INFINITY;
 526     while(size > classSet->size()){
 527         IndexSampleSetType*   newClassSet
 528                 = new IndexSampleSetType();
 529         newClassSet->reserve  ( indexSet->size() );
 530
 531         IndexSampleSetType*   newWeightsSet
 532                 = new IndexSampleSetType();
 533         newWeightsSet->reserve( indexSet->size() );
 534
 535         IndexSampleSetType*   newIndexSet
 536                 = new IndexSampleSetType();
 537         newIndexSet->reserve  ( indexSet->size() );
 538
 539         SpatialSampleSetType* newSpatialSet
 540                 = new SpatialSampleSetType();
 541
 542         newSpatialSet->reserve( indexSet->size() );
 543
 544         RangeSampleSetType*   newRangeSet
 545                 = new RangeSampleSetType();
 546         newRangeSet->reserve  ( indexSet->size() );
 547
 548         IndexSampleSetType* newClassSetMemory
 549                 = new IndexSampleSetType( *classSetMemory );
 550
 551         newClassSet->push_back  ( 1 );
 552         newWeightsSet->push_back( 1 );
 553         newIndexSet->push_back  ( indexSet->at(0)   );
 554         newSpatialSet->push_back( spatialSet->at(0) );
 555         newRangeSet->push_back  ( rangeSet->at(0)   );
 556
 557         unsigned int k;
 558         bool newC;
 559         bool  raNorm;
 560
 561         for(unsigned int i=1 ; i<classSet->size() ; ++i)
 562         {
 563             newC = true;
 564             k=0;
 565
 566             while( (k<newClassSet->size()) && newC )
 567             {
 568                 InfiniteNorm( raNorm, newRangeSet->at(k), rangeSet->at(i) );
 569                 if( raNorm ){
 570
 571                     ReaderPointer reader = ReaderType::New();
 572                     ImagePointer image = ImageType::New();
 573
 574                     reader->SetFileName( *stmsParams->images.begin() );
 575                     reader->Update();
 576                     image = reader->GetOutput();
 577                     image->FillBuffer( 0 );
 578
 579                     ImageIndexSetType* refClass  = new ImageIndexSetType();
 580                     ImageIndexSetType* candClass = new ImageIndexSetType();
 581                     refClass->reserve( classSet->size()/2 );
 582                     candClass->reserve( classSet->size()/2 );
 583
 584
 585                     for(unsigned int m=0 ; m<newClassSetMemory->size() ; ++m)
 586                     {
 587                         ImageIndexType idx;
 588                         if(newClassSetMemory->at(m) == classSet->at(i)){
 589
 590                             #pragma omp simd
 591                             for(unsigned int n=0 ; n<stmsParams->dim ; ++n)
 592                                 idx[n] = spatialSetMemory->at(m)[n]*stmsParams->spScales[n];
 593
 594                             candClass->push_back( idx );
 595                             image->SetPixel(candClass->back(), 1);
 596                         }
 597                         else
 598                         {
 599                             if(newClassSetMemory->at(m) == newClassSet->at(k)){
 600
 601                                 #pragma omp simd
 602                                 for(unsigned int n=0 ; n<stmsParams->dim ; ++n)
 603                                     idx[n] = spatialSetMemory->at(m)[n]*stmsParams->spScales[n];
 604
 605                                 refClass->push_back( idx );
 606                                 image->SetPixel(refClass->back(), 2);
 607                             }
 608                         }
 609                     }
 610
 611                     bool connex = false;
 612                     if(stmsParams->dim == 2)
 613                     {
 614                         if(refClass->size() < candClass->size())
 615                         {
 616                             unsigned int m=0;
 617                             ImageIndexType idx;
 618                             while((m<refClass->size()) && !connex)
 619                             {
 620                                 for(int x=-1 ; x<=1 ; ++x)
 621                                 {
 622                                     for(int y=-1 ; y<=1 ; ++y)
 623                                     {
 624                                         idx[0] = refClass->at(m)[0]+x;
 625                                         idx[1] = refClass->at(m)[1]+y;
 626
 627                                         if(idx[0]<image->GetBufferedRegion().GetSize()[0] && idx[0]>0 && idx[1]<image->GetBufferedRegion().GetSize()[1] && idx[1]>0)
 628                                         {
 629                                             if(image->GetPixel(idx) == 2)
 630                                             {
 631                                                 connex = true;
 632                                                 x = 2;
 633                                                 y = 2;
 634                                             }
 635                                         }
 636                                     }
 637                                 }
 638
 639                                 ++m;
 640                             }
 641                         }
 642                         else
 643                         {
 644                             unsigned int m=0;
 645                             ImageIndexType idx;
 646                             while((m<candClass->size()) && !connex)
 647                             {
 648                                 for(int x=-1 ; x<=1 ; ++x)
 649                                 {
 650                                     for(int y=-1 ; y<=1 ; ++y)
 651                                     {
 652                                         idx[0] = candClass->at(m)[0]+x;
 653                                         idx[1] = candClass->at(m)[1]+y;
 654
 655                                         if(idx[0]<image->GetBufferedRegion().GetSize()[0] && idx[0]>0 && idx[1]<image->GetBufferedRegion().GetSize()[1] && idx[1]>0)
 656                                         {
 657                                             if(image->GetPixel(idx) == 1)
 658                                             {
 659                                                 connex = true;
 660                                                 x = 2;
 661                                                 y = 2;
 662                                             }
 663                                         }
 664                                     }
 665                                 }
 666
 667                                 ++m;
 668                             }
 669                         }
 670
 671                     }
 672                     else
 673                     {
 674                         if(refClass->size() < candClass->size())
 675                         {
 676                             unsigned int m=0;
 677                             ImageIndexType idx;
 678                             while((m<refClass->size()) && !connex)
 679                             {
 680                                 for(int x=-1 ; x<=1 ; ++x)
 681                                 {
 682                                     for(int y=-1 ; y<=1 ; ++y)
 683                                     {
 684                                         for(int z=-1 ; z<=1 ; ++z)
 685                                         {
 686                                             idx[0] = refClass->at(m)[0]+x;
 687                                             idx[1] = refClass->at(m)[1]+y;
 688                                             idx[2] = refClass->at(m)[2]+z;
 689
 690                                             if(idx[0]<image->GetBufferedRegion().GetSize()[0] && idx[0]>0 && idx[1]<image->GetBufferedRegion().GetSize()[1] && idx[1]>0 && idx[2]<image->GetBufferedRegion().GetSize()[2] && idx[2]>0)
 691                                             {
 692                                                 if(image->GetPixel(idx) == 2)
 693                                                 {
 694                                                     connex = true;
 695                                                     x = 2;
 696                                                     y = 2;
 697                                                     z = 2;
 698                                                 }
 699                                             }
 700                                         }
 701                                     }
 702                                 }
 703
 704                                 ++m;
 705                             }
 706                         }
 707                         else
 708                         {
 709                             unsigned int m=0;
 710                             ImageIndexType idx;
 711                             while((m<candClass->size()) && !connex)
 712                             {
 713                                 for(int x=-1 ; x<=1 ; ++x)
 714                                 {
 715                                     for(int y=-1 ; y<=1 ; ++y)
 716                                     {
 717                                         for(int z=-1 ; z<=1 ; ++z)
 718                                         {
 719                                             idx[0] = candClass->at(m)[0]+x;
 720                                             idx[1] = candClass->at(m)[1]+y;
 721                                             idx[2] = candClass->at(m)[2]+z;
 722
 723                                             if(idx[0]<image->GetBufferedRegion().GetSize()[0] && idx[0]>0 && idx[1]<image->GetBufferedRegion().GetSize()[1] && idx[1]>0 && idx[2]<image->GetBufferedRegion().GetSize()[2] && idx[2]>0)
 724                                             {
 725                                                 if(image->GetPixel(idx) == 1)
 726                                                 {
 727                                                     connex = true;
 728                                                     x = 2;
 729                                                     y = 2;
 730                                                     z = 2;
 731                                                 }
 732                                             }
 733                                         }
 734                                     }
 735                                 }
 736
 737                                 ++m;
 738                             }
 739                         }
 740                     }
 741
 742                     if( connex )
 743                     {
 744
 745                         VectorWeightedAcc( newSpatialSet->at(k), newWeightsSet->at(k), spatialSet->at(i), weightsSet->at(i) );
 746                         VectorWeightedAcc( newRangeSet->at(k)  , newWeightsSet->at(k), rangeSet->at(i)  , weightsSet->at(i) );
 747
 748                         newWeightsSet->at(k) += weightsSet->at(i);
 749
 750                         VectorDiv(newSpatialSet->at(k), newWeightsSet->at(k));
 751                         VectorDiv(newRangeSet->at(k), newWeightsSet->at(k));
 752
 753                         for(unsigned int l=0 ; l< classSetMemory->size() ; ++l)
 754                         {
 755
 756                             if(classSetMemory->at(l) == classSet->at(i))
 757                                 newClassSetMemory->at(l) = newClassSet->at(k);
 758                         }
 759
 760                         newC = false;
 761                     }
 762                     else
 763                         ++k;
 764
 765                     delete refClass;
 766                     delete candClass;
 767                 }
 768                 else
 769                     ++k;
 770             }
 771
 772             if( newC )
 773             {
 774                 newClassSet->push_back( newClassSet->back()+1 );
 775                 newWeightsSet->push_back( weightsSet->at(i) );
 776                 newIndexSet->push_back  ( indexSet->at(i)   );
 777                 newSpatialSet->push_back( spatialSet->at(i) );
 778                 newRangeSet->push_back  ( rangeSet->at(i)   );
 779
 780                 for(unsigned int l=0 ; l< classSetMemory->size() ; ++l)
 781                 {
 782                     if(classSetMemory->at(l) == classSet->at(i))
 783                         newClassSetMemory->at(l) = newClassSet->back();
 784                 }
 785             }
 786         }
 787
 788         size = classSet->size();
 789
 790         classSet->swap  ( *newClassSet   );
 791         weightsSet->swap( *newWeightsSet );
 792         indexSet->swap  ( *newIndexSet   );
 793         spatialSet->swap( *newSpatialSet );
 794         rangeSet->swap  ( *newRangeSet   );
 795         classSetMemory->swap  ( *newClassSetMemory );
 796
 797         delete newClassSet;
 798         delete newWeightsSet;
 799         delete newIndexSet;
 800         delete newSpatialSet;
 801         delete newRangeSet;
 802         delete newClassSetMemory;
 803     }
 804 }
 805
 806
 807
 808 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 809 template< class T >
 810 void
 811 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 812 ::VectorDistance(float &dist, std::vector<T> &a, std::vector<T> &b)
 813 {
 814     dist = 0.0;
 815
 816     #pragma omp simd
 817     for(unsigned int i=0 ; i<a.size() ; ++i)
 818         dist += ( a[i]-b[i] )*( a[i]-b[i] );
 819 }
 820
 821
 822 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 823 void
 824 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 825 ::InfiniteNorm(bool &dist, RangeVectorType &a, RangeVectorType &b)
 826 {
 827     dist = true;
 828
 829     for(unsigned int i=0 ; i<a.size() ; ++i)
 830     {
 831         if(((a[i]-b[i])*(a[i]-b[i])) > 1)
 832         {
 833             dist = false;
 834             i = a.size();
 835         }
 836     }
 837 }
 838
 839
 840 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 841 void
 842 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 843 ::MergeInfiniteNorm(bool &dist, RangeVectorType &a, RangeVectorType &b)
 844 {
 845     dist = true;
 846
 847     for(unsigned int i=0 ; i<a.size() ; ++i)
 848     {
 849         if(((a[i]-b[i])*(a[i]-b[i])) > (1/mergeFactor))
 850         {
 851             dist = false;
 852             i = a.size();
 853         }
 854     }
 855 }
 856
 857
 858 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 859 template< class T >
 860 void
 861 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 862 ::VectorAcc(std::vector<T> &a, std::vector<T> &b)
 863 {
 864     #pragma omp simd
 865     for(unsigned int i=0 ; i<a.size() ; ++i)
 866         a[i] += b[i];
 867 }
 868
 869 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 870 template< class T >
 871 void
 872 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 873 ::VectorWeightedAcc(std::vector<T> &a, std::vector<T> &b, unsigned int &b_w)
 874 {
 875     #pragma omp simd
 876     for(unsigned int i=0 ; i<a.size() ; ++i)
 877         a[i] += b_w*b[i];
 878 }
 879
 880 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 881 template< class T >
 882 void
 883 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 884 ::VectorWeightedAcc(std::vector<T> &a, unsigned int &a_w, std::vector<T> &b, unsigned int &b_w)
 885 {
 886     #pragma omp simd
 887     for(unsigned int i=0 ; i<a.size() ; ++i){
 888         a[i] *= a_w;
 889         a[i] += b_w*b[i];
 890     }
 891 }
 892
 893 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 894 template< class T >
 895 void
 896 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 897 ::VectorWeightedMean(std::vector<T> &a, unsigned int &a_w, std::vector<T> &b, unsigned int &b_w)
 898 {
 899     #pragma omp simd
 900     for(unsigned int i=0 ; i<a.size() ; ++i)
 901         a[i] = ( (a[i]*a_w )+(b[i]*b_w) )/( a_w+b_w );
 902
 903     a_w += b_w;
 904 }
 905
 906
 907 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 908 template< class T >
 909 void
 910 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 911 ::VectorMul(std::vector<T> &a, unsigned int &coeff)
 912 {
 913     #pragma omp simd
 914     for(unsigned int i=0 ; i<a.size() ; ++i)
 915         a[i] *= coeff;
 916 }
 917
 918
 919 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 920 template< class T >
 921 void
 922 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 923 ::VectorDiv(std::vector<T> &a, unsigned int &coeff)
 924 {
 925     #pragma omp simd
 926     for(unsigned int i=0 ; i<a.size() ; ++i)
 927         a[i] /= coeff;
 928 }
 929
 930 } // end of namespace itkSTMS
 931
 932 #endif  // itkmsSTMS_BlurringSTMS_Spine_TXX