Lib/SpatioTemporalMeanShift/itkSTMS_BlurringSTMS.txx

   1 /*
   2  #
   3  #  File        : itkSTMS_BlurringSTMS.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.h"
  67 #include "itkSTMS_XMLFileParser.h"
  68
  69
  70 double gettime()
  71 {
  72     struct timespec timestamp;
  73
  74     clock_gettime(CLOCK_REALTIME, &timestamp);
  75     return timestamp.tv_sec * 1000.0 + timestamp.tv_nsec * 1.0e-6;
  76 }
  77
  78
  79 namespace itkSTMS
  80 {
  81 template < class IndexType, class SpatialType, class PixelType, class ImageType >
  82 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
  83 ::itkSTMS_BlurringSTMS( IndexSampleSetType* idx, IndexSampleSetType* cla, IndexSampleSetType* wei,
  84                         SpatialSampleSetType* sp, RangeSampleSetType* ra, ParametersType* params, ParserOutputType* desc)
  85 {
  86     this->indexSet   = idx;
  87     this->classSet   = cla;
  88     this->weightsSet = wei;
  89     this->spatialSet = sp;
  90     this->rangeSet   = ra;
  91
  92     this->classSetMemory   = new IndexSampleSetType  (  *this->classSet  );
  93     this->spatialSetMemory = new SpatialSampleSetType( *this->spatialSet );
  94
  95     this->stmsParams     = params;
  96     this->expDescription = desc;
  97
  98     this->mergeFactor = 1000;
  99 }
 100
 101
 102 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 103 void
 104 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 105 ::GenerateData()
 106 {
 107     double dtime;
 108
 109     if( !stmsParams->merge )
 110     {
 111         std::cout<< std::endl<<"No merge Filtering!"<<std::endl<<"numSamples: "<< indexSet->size() <<std::endl<< std::endl;
 112         NoMergeSTMSFiltering();
 113
 114         dtime=gettime();
 115         ClassificationNoMergeSTMSFiltering();
 116         dtime = gettime()-dtime;
 117         std::cout<< std::endl<<"Classif: " << dtime/1000 << " s" <<std::endl;
 118         FinalMerging();
 119     }
 120     else
 121     {
 122         std::cout<<std::endl<<"Merge Filtering!"<<std::endl<<"numSamples: "<< indexSet->size() <<std::endl<< std::endl;
 123         MergeSTMSFiltering();
 124         FinalMerging();
 125         //        ClassificationNoMergeSTMSFiltering();
 126     }
 127 }
 128
 129 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 130 void
 131 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 132 ::ClassificationNoMergeSTMSFiltering()
 133 {
 134
 135     IndexSampleSetType*   newClassSet
 136             = new IndexSampleSetType();
 137     newClassSet->reserve  ( indexSet->size() );
 138
 139     IndexSampleSetType*   newWeightsSet
 140             = new IndexSampleSetType();
 141     newWeightsSet->reserve( indexSet->size() );
 142
 143     IndexSampleSetType*   newIndexSet
 144             = new IndexSampleSetType();
 145     newIndexSet->reserve  ( indexSet->size() );
 146
 147     SpatialSampleSetType* newSpatialSet
 148             = new SpatialSampleSetType();
 149
 150     newSpatialSet->reserve( indexSet->size() );
 151
 152     RangeSampleSetType*   newRangeSet
 153             = new RangeSampleSetType();
 154     newRangeSet->reserve  ( indexSet->size() );
 155
 156
 157     newClassSet->push_back  ( 1 );
 158     newWeightsSet->push_back( 1 );
 159     newIndexSet->push_back  ( indexSet->at(0)   );
 160     newSpatialSet->push_back( spatialSet->at(0) );
 161     newRangeSet->push_back  ( rangeSet->at(0)   );
 162     classSetMemory->at( indexSet->at(0) ) = newClassSet->at(0);
 163
 164
 165     unsigned int k;
 166     bool newC;
 167     float spDist;
 168     bool  raNorm;
 169
 170     for(unsigned int i=1 ; i<indexSet->size() ; ++i)
 171     {
 172         newC = true;
 173         k=0;
 174
 175         while( (k<newClassSet->size()) && newC )
 176         {
 177             VectorDistance( spDist, newSpatialSet->at(k), spatialSet->at(i) );
 178             if( spDist<=1 )
 179             {
 180                 InfiniteNorm( raNorm, newRangeSet->at(k), rangeSet->at(i) );
 181                 if( raNorm )
 182                 {
 183                     VectorWeightedMean( newSpatialSet->at(k), newWeightsSet->at(k), spatialSet->at(i), weightsSet->at(i) );
 184                     VectorWeightedMean( newRangeSet->at(k)  , newWeightsSet->at(k), rangeSet->at(i)  , weightsSet->at(i) );
 185
 186                     classSetMemory->at( indexSet->at(i) ) = newClassSet->at(k);
 187                     newC = false;
 188                 }
 189                 else
 190                     ++k;
 191             }
 192             else
 193                 ++k;
 194         }
 195
 196         if( newC )
 197         {
 198             newClassSet->push_back( (IndexType)newClassSet->back()+1 );
 199             newWeightsSet->push_back( 1 );
 200             newIndexSet->push_back  ( indexSet->at(i)   );
 201             newSpatialSet->push_back( spatialSet->at(i) );
 202             newRangeSet->push_back  ( rangeSet->at(i)   );
 203             classSetMemory->at( indexSet->at(i) ) = newClassSet->back();
 204         }
 205     }
 206
 207     classSet->swap  ( *newClassSet   );
 208     weightsSet->swap( *newWeightsSet );
 209     indexSet->swap  ( *newIndexSet   );
 210     spatialSet->swap( *newSpatialSet );
 211     rangeSet->swap  ( *newRangeSet   );
 212
 213     delete newClassSet;
 214     delete newWeightsSet;
 215     delete newIndexSet;
 216     delete newSpatialSet;
 217     delete newRangeSet;
 218 }
 219
 220
 221 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 222 void
 223 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 224 ::NoMergeSTMSFiltering()
 225 {
 226
 227     SpatialSampleSetType* newSpatialSet
 228             = new SpatialSampleSetType( indexSet->size(), SpatialVectorType(stmsParams->dim, 0) );
 229     RangeSampleSetType*   newRangeSet
 230             = new RangeSampleSetType( indexSet->size(), RangeVectorType(stmsParams->numTimePoints, 0) );
 231
 232     float globalEvolution = INFINITY;
 233     float epsilon = (stmsParams->epsilon)*(stmsParams->epsilon);
 234     float spDist, raDist;
 235     bool  raNorm;
 236     unsigned int iter = 0;
 237
 238     double dtime;
 239     while( (globalEvolution>epsilon) & (iter<stmsParams->maxIt) )
 240     {
 241         dtime=gettime();
 242         globalEvolution = 0.0;
 243
 244         #pragma omp parallel for reduction(+:globalEvolution) private(spDist, raDist, raNorm)
 245         for(unsigned int i=0 ; i<indexSet->size() ; ++i)
 246         {
 247             unsigned int count = 0;
 248             SpatialVectorType tmpSpatial( stmsParams->dim, 0 );
 249             RangeVectorType   tmpRange  ( stmsParams->numTimePoints, 0 );
 250
 251             for(unsigned int j=0 ; j<indexSet->size() ; ++j)
 252             {
 253                 VectorDistance( spDist, spatialSet->at(i), spatialSet->at(j));
 254
 255                 if( spDist<=1 )
 256                 {
 257                     InfiniteNorm(raNorm, rangeSet->at(i), rangeSet->at(j));
 258
 259                     if( raNorm )
 260                     {
 261                         VectorAcc( tmpSpatial, spatialSet->at(j) );
 262                         VectorAcc( tmpRange  , rangeSet->at(j)   );
 263
 264                         ++count;
 265                     }
 266                 }
 267             }
 268
 269             if( count>1 )
 270             {
 271                 VectorDiv( tmpSpatial, count );
 272                 VectorDiv( tmpRange  , count );
 273             }
 274
 275             VectorDistance(spDist, spatialSet->at(i), tmpSpatial);
 276             VectorDistance(raDist, rangeSet->at(i), tmpRange);
 277
 278             globalEvolution += (spDist+raDist);
 279
 280             #pragma omp critical
 281             {
 282                 newSpatialSet->at(i) = tmpSpatial;
 283                 newRangeSet->at(i)   = tmpRange;
 284             }
 285         }
 286
 287         spatialSet->swap( *newSpatialSet );
 288         rangeSet->swap  (  *newRangeSet  );
 289
 290         ++iter;
 291         dtime = gettime()-dtime;
 292         std::cout<< "Iter: " << iter <<"   "<< dtime/1000 << " s" <<"    GE: "<< globalEvolution <<std::endl;
 293     }
 294
 295     delete newSpatialSet;
 296     delete newRangeSet;
 297 }
 298
 299
 300
 301
 302 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 303 void
 304 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 305 ::MergeSTMSFiltering()
 306 {
 307
 308     float globalEvolution = INFINITY;
 309     float epsilon = (stmsParams->epsilon)*(stmsParams->epsilon);
 310     float spDist, raDist;
 311     bool  raNorm, mergeNorm;
 312     unsigned int iter = 0;
 313     unsigned int merging;
 314
 315     double dtime;
 316     while( (globalEvolution>epsilon) & (iter<stmsParams->maxIt) )
 317     {
 318         SpatialSampleSetType* newSpatialSet
 319                 = new SpatialSampleSetType( *spatialSet );
 320         RangeSampleSetType*   newRangeSet
 321                 = new RangeSampleSetType  (  *rangeSet  );
 322
 323         dtime=gettime();
 324         globalEvolution = 0.0;
 325         merging = 0;
 326
 327         #pragma omp parallel for reduction(+:globalEvolution, merging) private(spDist, raDist, raNorm, mergeNorm)
 328         for(unsigned int i=0 ; i<indexSet->size() ; ++i)
 329         {
 330             unsigned int weight = 0;
 331             bool merge = false;
 332             SpatialVectorType tmpSpatial( stmsParams->dim, 0 );
 333             RangeVectorType   tmpRange  ( stmsParams->numTimePoints, 0 );
 334
 335             for(unsigned int j=0 ; j<indexSet->size() ; ++j)
 336             {
 337                 VectorDistance( spDist, spatialSet->at(i), spatialSet->at(j));
 338
 339                 if( spDist<=1 )
 340                 {
 341                     InfiniteNorm(raNorm, rangeSet->at(i), rangeSet->at(j));
 342
 343                     if( raNorm )
 344                     {
 345                         VectorWeightedAcc(tmpSpatial, spatialSet->at(j), weightsSet->at(j));
 346                         VectorWeightedAcc(tmpRange  , rangeSet->at(j)  , weightsSet->at(j));
 347
 348                         weight += weightsSet->at(j);
 349
 350                         MergeInfiniteNorm(mergeNorm, rangeSet->at(i), rangeSet->at(j));
 351                         if( (spDist <= (1/mergeFactor)) && mergeNorm )
 352                             merge = true;
 353                     }
 354                 }
 355             }
 356
 357             if(weight > weightsSet->at(i))
 358             {
 359                 if( merge )
 360                     merging += 1;
 361                 else
 362                     merge = 0;
 363             }
 364             else
 365                 merge = 0;
 366
 367
 368             VectorDiv( tmpSpatial, weight );
 369             VectorDiv( tmpRange  , weight );
 370
 371             #pragma omp critical
 372             {
 373                 newSpatialSet->at(i) = tmpSpatial;
 374                 newRangeSet->at(i)   = tmpRange;
 375             }
 376
 377             VectorDistance(spDist, spatialSet->at(i), tmpSpatial);
 378             VectorDistance(raDist, rangeSet->at(i), tmpRange);
 379
 380             globalEvolution += (spDist+raDist);
 381         }
 382
 383         spatialSet->swap( *newSpatialSet );
 384         rangeSet->swap  (  *newRangeSet  );
 385
 386         delete newSpatialSet;
 387         delete newRangeSet;
 388
 389         if(merging > 0)
 390             MergeSamples();
 391
 392         ++iter;
 393         dtime = gettime()-dtime;
 394         std::cout<< "Iter: " << iter <<"   "<< dtime/1000 << " s" <<"    GE: "<< globalEvolution << "      numSamples: "<< classSet->size() <<std::endl;
 395     }
 396 }
 397
 398
 399 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 400 void
 401 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 402 ::MergeSamples()
 403 {
 404
 405     IndexSampleSetType*   newClassSet
 406             = new IndexSampleSetType();
 407     newClassSet->reserve  ( indexSet->size() );
 408
 409     IndexSampleSetType*   newWeightsSet
 410             = new IndexSampleSetType();
 411     newWeightsSet->reserve( indexSet->size() );
 412
 413     IndexSampleSetType*   newIndexSet
 414             = new IndexSampleSetType();
 415     newIndexSet->reserve  ( indexSet->size() );
 416
 417     SpatialSampleSetType* newSpatialSet
 418             = new SpatialSampleSetType();
 419
 420     newSpatialSet->reserve( indexSet->size() );
 421
 422     RangeSampleSetType*   newRangeSet
 423             = new RangeSampleSetType();
 424     newRangeSet->reserve  ( indexSet->size() );
 425
 426     IndexSampleSetType* newClassSetMemory
 427             = new IndexSampleSetType( *classSetMemory );
 428
 429     IndexSampleSetType*   indexes
 430             = new IndexSampleSetType( *classSet );
 431
 432     std::random_device rd;
 433     std::mt19937 g(rd());
 434
 435     std::shuffle(indexes->begin(), indexes->end(), g);
 436
 437
 438     newClassSet->push_back  ( 1 );
 439     newWeightsSet->push_back( weightsSet->at( indexes->at(0)-1 ) );
 440     newIndexSet->push_back  ( indexSet->at  ( indexes->at(0)-1 ) );
 441     newSpatialSet->push_back( spatialSet->at( indexes->at(0)-1 ) );
 442     newRangeSet->push_back  ( rangeSet->at  ( indexes->at(0)-1 ) );
 443
 444     for(unsigned int l=0 ; l< classSetMemory->size() ; ++l)
 445     {
 446         if(classSetMemory->at(l) == classSet->at(indexes->at(0)-1))
 447             newClassSetMemory->at(l) = 1;
 448     }
 449
 450     unsigned int k;
 451     bool newC;
 452     float spDist;
 453     bool  raNorm;
 454     unsigned int i;
 455
 456     for(unsigned int m=1 ; m<indexes->size() ; ++m)
 457     {
 458         i = indexes->at(m)-1;
 459
 460         newC = true;
 461         k=0;
 462
 463         while( (k<newClassSet->size()) && newC )
 464         {
 465             VectorDistance( spDist, newSpatialSet->at(k), spatialSet->at(i) );
 466             MergeInfiniteNorm( raNorm, newRangeSet->at(k), rangeSet->at(i) );
 467             if( (spDist<=1/mergeFactor)&&raNorm )
 468             {
 469                 VectorWeightedAcc( newSpatialSet->at(k), newWeightsSet->at(k), spatialSet->at(i), weightsSet->at(i) );
 470                 VectorWeightedAcc( newRangeSet->at(k)  , newWeightsSet->at(k), rangeSet->at(i)  , weightsSet->at(i) );
 471
 472                 newWeightsSet->at(k) += weightsSet->at(i);
 473
 474                 VectorDiv(newSpatialSet->at(k), newWeightsSet->at(k));
 475                 VectorDiv(newRangeSet->at(k), newWeightsSet->at(k));
 476
 477                 for(unsigned int l=0 ; l<classSetMemory->size() ; ++l)
 478                 {
 479                     if(classSetMemory->at(l) == classSet->at(i))
 480                         newClassSetMemory->at(l) = newClassSet->at(k);
 481                 }
 482
 483                 newC = false;
 484             }
 485             else
 486                 ++k;
 487         }
 488
 489         if( newC )
 490         {
 491             newClassSet->push_back( newClassSet->back()+1 );
 492             newWeightsSet->push_back( weightsSet->at(i) );
 493             newIndexSet->push_back  ( indexSet->at(i)   );
 494             newSpatialSet->push_back( spatialSet->at(i) );
 495             newRangeSet->push_back  ( rangeSet->at(i)   );
 496
 497             for(unsigned int l=0 ; l< classSetMemory->size() ; ++l)
 498             {
 499                 if(classSetMemory->at(l) == classSet->at(i))
 500                     newClassSetMemory->at(l) = newClassSet->back();
 501             }
 502         }
 503     }
 504
 505     classSet->swap  ( *newClassSet   );
 506     weightsSet->swap( *newWeightsSet );
 507     indexSet->swap  ( *newIndexSet   );
 508     spatialSet->swap( *newSpatialSet );
 509     rangeSet->swap  ( *newRangeSet   );
 510     classSetMemory->swap  ( *newClassSetMemory );
 511
 512     delete newClassSet;
 513     delete newWeightsSet;
 514     delete newIndexSet;
 515     delete newSpatialSet;
 516     delete newRangeSet;
 517     delete newClassSetMemory;
 518     delete indexes;
 519 }
 520
 521
 522 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 523 void
 524 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 525 ::FinalMerging()
 526 {
 527     unsigned int size = INFINITY;
 528     std::string imagePath = expDescription->experimentPath+expDescription->inputFolder+expDescription->inputCommonRoot+STMS_NUMBERING_FORM_ONE+expDescription->imageExtension;
 529
 530     while(size > classSet->size()){
 531         IndexSampleSetType*   newClassSet
 532                 = new IndexSampleSetType();
 533         newClassSet->reserve  ( indexSet->size() );
 534
 535         IndexSampleSetType*   newWeightsSet
 536                 = new IndexSampleSetType();
 537         newWeightsSet->reserve( indexSet->size() );
 538
 539         IndexSampleSetType*   newIndexSet
 540                 = new IndexSampleSetType();
 541         newIndexSet->reserve  ( indexSet->size() );
 542
 543         SpatialSampleSetType* newSpatialSet
 544                 = new SpatialSampleSetType();
 545
 546         newSpatialSet->reserve( indexSet->size() );
 547
 548         RangeSampleSetType*   newRangeSet
 549                 = new RangeSampleSetType();
 550         newRangeSet->reserve  ( indexSet->size() );
 551
 552         IndexSampleSetType* newClassSetMemory
 553                 = new IndexSampleSetType( *classSetMemory );
 554
 555         newClassSet->push_back  ( 1 );
 556         newWeightsSet->push_back( 1 );
 557         newIndexSet->push_back  ( indexSet->at(0)   );
 558         newSpatialSet->push_back( spatialSet->at(0) );
 559         newRangeSet->push_back  ( rangeSet->at(0)   );
 560
 561         unsigned int k;
 562         bool newC;
 563         bool  raNorm;
 564
 565         for(unsigned int i=1 ; i<classSet->size() ; ++i)
 566         {
 567             newC = true;
 568             k=0;
 569
 570             while( (k<newClassSet->size()) && newC )
 571             {
 572                 InfiniteNorm( raNorm, newRangeSet->at(k), rangeSet->at(i) );
 573                 if( raNorm ){
 574
 575                     ReaderPointer reader = ReaderType::New();
 576                     ImagePointer image = ImageType::New();
 577
 578                     reader->SetFileName( imagePath );
 579                     reader->Update();
 580                     image = reader->GetOutput();
 581                     image->FillBuffer( 0 );
 582
 583                     ImageIndexSetType* refClass  = new ImageIndexSetType();
 584                     ImageIndexSetType* candClass = new ImageIndexSetType();
 585                     refClass->reserve( classSet->size()/2 );
 586                     candClass->reserve( classSet->size()/2 );
 587
 588
 589                     for(unsigned int m=0 ; m<newClassSetMemory->size() ; ++m)
 590                     {
 591                         ImageIndexType idx;
 592                         if(newClassSetMemory->at(m) == classSet->at(i)){
 593
 594                             #pragma omp simd
 595                             for(unsigned int n=0 ; n<stmsParams->dim ; ++n)
 596                                 idx[n] = spatialSetMemory->at(m)[n]*stmsParams->spScales[n];
 597
 598                             candClass->push_back( idx );
 599                             image->SetPixel(candClass->back(), 1);
 600                         }
 601                         else
 602                         {
 603                             if(newClassSetMemory->at(m) == newClassSet->at(k)){
 604
 605                                 #pragma omp simd
 606                                 for(unsigned int n=0 ; n<stmsParams->dim ; ++n)
 607                                     idx[n] = spatialSetMemory->at(m)[n]*stmsParams->spScales[n];
 608
 609                                 refClass->push_back( idx );
 610                                 image->SetPixel(refClass->back(), 2);
 611                             }
 612                         }
 613                     }
 614
 615                     bool connex = false;
 616                     if(stmsParams->dim == 2)
 617                     {
 618                         if(refClass->size() < candClass->size())
 619                         {
 620                             unsigned int m=0;
 621                             ImageIndexType idx;
 622                             while((m<refClass->size()) && !connex)
 623                             {
 624                                 for(int x=-1 ; x<=1 ; ++x)
 625                                 {
 626                                     for(int y=-1 ; y<=1 ; ++y)
 627                                     {
 628                                         idx[0] = refClass->at(m)[0]+x;
 629                                         idx[1] = refClass->at(m)[1]+y;
 630
 631                                         if(idx[0]<image->GetBufferedRegion().GetSize()[0] && idx[0]>0 && idx[1]<image->GetBufferedRegion().GetSize()[1] && idx[1]>0)
 632                                         {
 633                                             if(image->GetPixel(idx) == 2)
 634                                             {
 635                                                 connex = true;
 636                                                 x = 2;
 637                                                 y = 2;
 638                                             }
 639                                         }
 640                                     }
 641                                 }
 642
 643                                 ++m;
 644                             }
 645                         }
 646                         else
 647                         {
 648                             unsigned int m=0;
 649                             ImageIndexType idx;
 650                             while((m<candClass->size()) && !connex)
 651                             {
 652                                 for(int x=-1 ; x<=1 ; ++x)
 653                                 {
 654                                     for(int y=-1 ; y<=1 ; ++y)
 655                                     {
 656                                         idx[0] = candClass->at(m)[0]+x;
 657                                         idx[1] = candClass->at(m)[1]+y;
 658
 659                                         if(idx[0]<image->GetBufferedRegion().GetSize()[0] && idx[0]>0 && idx[1]<image->GetBufferedRegion().GetSize()[1] && idx[1]>0)
 660                                         {
 661                                             if(image->GetPixel(idx) == 1)
 662                                             {
 663                                                 connex = true;
 664                                                 x = 2;
 665                                                 y = 2;
 666                                             }
 667                                         }
 668                                     }
 669                                 }
 670
 671                                 ++m;
 672                             }
 673                         }
 674
 675                     }
 676                     else
 677                     {
 678                         if(refClass->size() < candClass->size())
 679                         {
 680                             unsigned int m=0;
 681                             ImageIndexType idx;
 682                             while((m<refClass->size()) && !connex)
 683                             {
 684                                 for(int x=-1 ; x<=1 ; ++x)
 685                                 {
 686                                     for(int y=-1 ; y<=1 ; ++y)
 687                                     {
 688                                         for(int z=-1 ; z<=1 ; ++z)
 689                                         {
 690                                             idx[0] = refClass->at(m)[0]+x;
 691                                             idx[1] = refClass->at(m)[1]+y;
 692                                             idx[2] = refClass->at(m)[2]+z;
 693
 694                                             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)
 695                                             {
 696                                                 if(image->GetPixel(idx) == 2)
 697                                                 {
 698                                                     connex = true;
 699                                                     x = 2;
 700                                                     y = 2;
 701                                                     z = 2;
 702                                                 }
 703                                             }
 704                                         }
 705                                     }
 706                                 }
 707
 708                                 ++m;
 709                             }
 710                         }
 711                         else
 712                         {
 713                             unsigned int m=0;
 714                             ImageIndexType idx;
 715                             while((m<candClass->size()) && !connex)
 716                             {
 717                                 for(int x=-1 ; x<=1 ; ++x)
 718                                 {
 719                                     for(int y=-1 ; y<=1 ; ++y)
 720                                     {
 721                                         for(int z=-1 ; z<=1 ; ++z)
 722                                         {
 723                                             idx[0] = candClass->at(m)[0]+x;
 724                                             idx[1] = candClass->at(m)[1]+y;
 725                                             idx[2] = candClass->at(m)[2]+z;
 726
 727                                             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)
 728                                             {
 729                                                 if(image->GetPixel(idx) == 1)
 730                                                 {
 731                                                     connex = true;
 732                                                     x = 2;
 733                                                     y = 2;
 734                                                     z = 2;
 735                                                 }
 736                                             }
 737                                         }
 738                                     }
 739                                 }
 740
 741                                 ++m;
 742                             }
 743                         }
 744                     }
 745
 746                     if( connex )
 747                     {
 748
 749                         VectorWeightedAcc( newSpatialSet->at(k), newWeightsSet->at(k), spatialSet->at(i), weightsSet->at(i) );
 750                         VectorWeightedAcc( newRangeSet->at(k)  , newWeightsSet->at(k), rangeSet->at(i)  , weightsSet->at(i) );
 751
 752                         newWeightsSet->at(k) += weightsSet->at(i);
 753
 754                         VectorDiv(newSpatialSet->at(k), newWeightsSet->at(k));
 755                         VectorDiv(newRangeSet->at(k), newWeightsSet->at(k));
 756
 757                         for(unsigned int l=0 ; l< classSetMemory->size() ; ++l)
 758                         {
 759
 760                             if(classSetMemory->at(l) == classSet->at(i))
 761                                 newClassSetMemory->at(l) = newClassSet->at(k);
 762                         }
 763
 764                         newC = false;
 765                     }
 766                     else
 767                         ++k;
 768
 769                     delete refClass;
 770                     delete candClass;
 771                 }
 772                 else
 773                     ++k;
 774             }
 775
 776             if( newC )
 777             {
 778                 newClassSet->push_back( newClassSet->back()+1 );
 779                 newWeightsSet->push_back( weightsSet->at(i) );
 780                 newIndexSet->push_back  ( indexSet->at(i)   );
 781                 newSpatialSet->push_back( spatialSet->at(i) );
 782                 newRangeSet->push_back  ( rangeSet->at(i)   );
 783
 784                 for(unsigned int l=0 ; l< classSetMemory->size() ; ++l)
 785                 {
 786                     if(classSetMemory->at(l) == classSet->at(i))
 787                         newClassSetMemory->at(l) = newClassSet->back();
 788                 }
 789             }
 790         }
 791
 792         size = classSet->size();
 793
 794         classSet->swap  ( *newClassSet   );
 795         weightsSet->swap( *newWeightsSet );
 796         indexSet->swap  ( *newIndexSet   );
 797         spatialSet->swap( *newSpatialSet );
 798         rangeSet->swap  ( *newRangeSet   );
 799         classSetMemory->swap  ( *newClassSetMemory );
 800
 801         delete newClassSet;
 802         delete newWeightsSet;
 803         delete newIndexSet;
 804         delete newSpatialSet;
 805         delete newRangeSet;
 806         delete newClassSetMemory;
 807     }
 808 }
 809
 810
 811
 812 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 813 template< class T >
 814 void
 815 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 816 ::VectorDistance(float &dist, std::vector<T> &a, std::vector<T> &b)
 817 {
 818     dist = 0.0;
 819
 820     #pragma omp simd
 821     for(unsigned int i=0 ; i<a.size() ; ++i)
 822         dist += ( a[i]-b[i] )*( a[i]-b[i] );
 823 }
 824
 825
 826 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 827 void
 828 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 829 ::InfiniteNorm(bool &dist, RangeVectorType &a, RangeVectorType &b)
 830 {
 831     dist = true;
 832
 833     for(unsigned int i=0 ; i<a.size() ; ++i)
 834     {
 835         if(((a[i]-b[i])*(a[i]-b[i])) > 1)
 836         {
 837             dist = false;
 838             i = a.size();
 839         }
 840     }
 841 }
 842
 843
 844 template < class IndexType, class SpatialType, class PixelType, class ImageType >
 845 void
 846 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 847 ::MergeInfiniteNorm(bool &dist, RangeVectorType &a, RangeVectorType &b)
 848 {
 849     dist = true;
 850
 851     for(unsigned int i=0 ; i<a.size() ; ++i)
 852     {
 853         if(((a[i]-b[i])*(a[i]-b[i])) > (1/mergeFactor))
 854         {
 855             dist = false;
 856             i = a.size();
 857         }
 858     }
 859 }
 860
 861
 862 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 863 template< class T >
 864 void
 865 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 866 ::VectorAcc(std::vector<T> &a, std::vector<T> &b)
 867 {
 868     #pragma omp simd
 869     for(unsigned int i=0 ; i<a.size() ; ++i)
 870         a[i] += b[i];
 871 }
 872
 873 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 874 template< class T >
 875 void
 876 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 877 ::VectorWeightedAcc(std::vector<T> &a, std::vector<T> &b, unsigned int &b_w)
 878 {
 879     #pragma omp simd
 880     for(unsigned int i=0 ; i<a.size() ; ++i)
 881         a[i] += b_w*b[i];
 882 }
 883
 884 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 885 template< class T >
 886 void
 887 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 888 ::VectorWeightedAcc(std::vector<T> &a, unsigned int &a_w, std::vector<T> &b, unsigned int &b_w)
 889 {
 890     #pragma omp simd
 891     for(unsigned int i=0 ; i<a.size() ; ++i){
 892         a[i] *= a_w;
 893         a[i] += b_w*b[i];
 894     }
 895 }
 896
 897 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 898 template< class T >
 899 void
 900 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 901 ::VectorWeightedMean(std::vector<T> &a, unsigned int &a_w, std::vector<T> &b, unsigned int &b_w)
 902 {
 903     #pragma omp simd
 904     for(unsigned int i=0 ; i<a.size() ; ++i)
 905         a[i] = ( (a[i]*a_w )+(b[i]*b_w) )/( a_w+b_w );
 906
 907     a_w += b_w;
 908 }
 909
 910
 911 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 912 template< class T >
 913 void
 914 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 915 ::VectorMul(std::vector<T> &a, unsigned int &coeff)
 916 {
 917     #pragma omp simd
 918     for(unsigned int i=0 ; i<a.size() ; ++i)
 919         a[i] *= coeff;
 920 }
 921
 922
 923 template< class IndexType, class SpatialType, class PixelType, class ImageType >
 924 template< class T >
 925 void
 926 itkSTMS_BlurringSTMS< IndexType, SpatialType, PixelType, ImageType >
 927 ::VectorDiv(std::vector<T> &a, unsigned int &coeff)
 928 {
 929     #pragma omp simd
 930     for(unsigned int i=0 ; i<a.size() ; ++i)
 931         a[i] /= coeff;
 932 }
 933
 934 } // end of namespace itkSTMS
 935
 936 #endif  // itkmsSTMS_BlurringSTMS_TXX