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