itk/clitkExplosionControlledThresholdConnectedImageFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================*/
  18 #ifndef __clitkExplosionControlledThresholdConnectedImageFilter_txx
  19 #define __clitkExplosionControlledThresholdConnectedImageFilter_txx
  20
  21 #include "clitkExplosionControlledThresholdConnectedImageFilter.h"
  22 #include "itkBinaryThresholdImageFunction.h"
  23 #include "itkFloodFilledImageFunctionConditionalIterator.h"
  24 #include "itkShapedFloodFilledImageFunctionConditionalIterator.h"
  25 #include "itkProgressReporter.h"
  26
  27 namespace clitk
  28 {
  29
  30   //--------------------------------------------------------------------
  31   /**
  32    * Constructor
  33    */
  34   template <class TInputImage, class TOutputImage>
  35   ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
  36   ::ExplosionControlledThresholdConnectedImageFilter()
  37   {
  38     m_Lower = itk::NumericTraits<InputImagePixelType>::NonpositiveMin();
  39     m_Upper = itk::NumericTraits<InputImagePixelType>::max();
  40     m_FinalLower = itk::NumericTraits<InputImagePixelType>::NonpositiveMin();
  41     m_FinalUpper = itk::NumericTraits<InputImagePixelType>::max();
  42     m_AdaptLowerBorder=false;
  43     m_AdaptUpperBorder=false;
  44     m_MaximumUpperThreshold=itk::NumericTraits<InputImagePixelType>::max();
  45     m_MinimumLowerThreshold=itk::NumericTraits<InputImagePixelType>::NonpositiveMin();
  46     m_Multiplier=2.0;
  47     m_ReplaceValue = itk::NumericTraits<OutputImagePixelType>::One;
  48     m_MinimumThresholdStepSize=itk::NumericTraits<InputImagePixelType>::One;
  49     m_ThresholdStepSize=64;
  50     m_Verbose=false;
  51     m_FullyConnected=false;
  52     m_MinimumSize=5000;
  53   }
  54   //--------------------------------------------------------------------
  55
  56
  57   //--------------------------------------------------------------------
  58   template <class TInputImage, class TOutputImage>
  59   void
  60   ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
  61   ::ClearSeeds()
  62   {
  63     if( this->m_Seeds.size() > 0 )
  64       {
  65         this->m_Seeds.clear();
  66         this->Modified();
  67       }
  68   }
  69   //--------------------------------------------------------------------
  70
  71
  72   //--------------------------------------------------------------------
  73   template <class TInputImage, class TOutputImage>
  74   void
  75   ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
  76   ::SetSeed(const IndexType & seed)
  77   {
  78     this->ClearSeeds();
  79     this->AddSeed ( seed );
  80   }
  81   //--------------------------------------------------------------------
  82
  83
  84   //--------------------------------------------------------------------
  85   template <class TInputImage, class TOutputImage>
  86   void
  87   ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
  88   ::AddSeed ( const IndexType & seed )
  89   {
  90     this->m_Seeds.push_back ( seed );
  91     this->Modified();
  92   }
  93   //--------------------------------------------------------------------
  94
  95
  96   //--------------------------------------------------------------------
  97   /**
  98    * Standard PrintSelf method.
  99    */
 100   template <class TInputImage, class TOutputImage>
 101   void
 102   ExplosionControlledThresholdConnectedImageFilter<TInputImage, TOutputImage>
 103   ::PrintSelf(std::ostream& os, itk::Indent indent) const
 104   {
 105     this->Superclass::PrintSelf(os, indent);
 106     os << indent << "Upper: "
 107        << static_cast<typename itk::NumericTraits<InputImagePixelType>::PrintType>(m_Upper)
 108        << std::endl;
 109     os << indent << "Lower: "
 110        << static_cast<typename itk::NumericTraits<InputImagePixelType>::PrintType>(m_Lower)
 111        << std::endl;
 112     os << indent << "ReplaceValue: "
 113        << static_cast<double>(m_Multiplier)
 114        << std::endl;
 115     os << indent << "ReplaceValue: "
 116        << static_cast<typename itk::NumericTraits<OutputImagePixelType>::PrintType>(m_ReplaceValue)
 117        << std::endl;
 118   }
 119   //--------------------------------------------------------------------
 120
 121
 122   //--------------------------------------------------------------------
 123   template <class TInputImage, class TOutputImage>
 124   void
 125   ExplosionControlledThresholdConnectedImageFilter<TInputImage,TOutputImage>
 126   ::GenerateInputRequestedRegion()
 127   {
 128     Superclass::GenerateInputRequestedRegion();
 129     if ( this->GetInput() )
 130       {
 131         InputImagePointer image =
 132           const_cast< InputImageType * >( this->GetInput() );
 133         image->SetRequestedRegionToLargestPossibleRegion();
 134       }
 135   }
 136   //--------------------------------------------------------------------
 137
 138
 139   //--------------------------------------------------------------------
 140   template <class TInputImage, class TOutputImage>
 141   void
 142   ExplosionControlledThresholdConnectedImageFilter<TInputImage,TOutputImage>
 143   ::EnlargeOutputRequestedRegion(itk::DataObject *output)
 144   {
 145     Superclass::EnlargeOutputRequestedRegion(output);
 146     output->SetRequestedRegionToLargestPossibleRegion();
 147   }
 148   //--------------------------------------------------------------------
 149
 150
 151   //--------------------------------------------------------------------
 152   template <class TInputImage, class TOutputImage>
 153   void
 154   ExplosionControlledThresholdConnectedImageFilter<TInputImage,TOutputImage>
 155   ::GenerateData()
 156   {
 157     typename Superclass::InputImageConstPointer inputImage  = this->GetInput();
 158     typename Superclass::OutputImagePointer     outputImage = this->GetOutput();
 159
 160     // Zero the output
 161     outputImage->SetBufferedRegion( outputImage->GetRequestedRegion() );
 162     outputImage->Allocate();
 163     outputImage->FillBuffer ( itk::NumericTraits<OutputImagePixelType>::Zero );
 164
 165
 166     //--------------------------------------------------
 167     // Get initial region size
 168     //--------------------------------------------------
 169     typedef itk::BinaryThresholdImageFunction<InputImageType> FunctionType;
 170     //typedef itk::ShapedFloodFilledImageFunctionConditionalIterator<OutputImageType, FunctionType> IteratorType;
 171     typedef itk::FloodFilledImageFunctionConditionalIterator<OutputImageType, FunctionType> IteratorType;
 172     typename FunctionType::Pointer function = FunctionType::New();
 173     function->SetInputImage ( inputImage );
 174     function->ThresholdBetween ( m_Lower, m_Upper );
 175     IteratorType it = IteratorType ( outputImage, function, m_Seeds );
 176     //it.SetFullyConnected(m_FullyConnected);
 177     unsigned int initialSize=0;
 178
 179     it.GoToBegin();
 180     while( !it.IsAtEnd())
 181       {
 182         ++initialSize;
 183         ++it;
 184       }
 185     if (m_Verbose)  std::cout<<"Initial region size using thresholds ["<<m_Lower<<", "<<m_Upper<<"], is "<<initialSize<<"..."<<std::endl;
 186
 187
 188     //--------------------------------------------------
 189     // Decrease lower threshold
 190     //--------------------------------------------------
 191     m_FinalLower=m_Lower;
 192     if (m_AdaptLowerBorder)
 193       {
 194         // Init
 195         InputImagePixelType currentLower=m_Lower;
 196         unsigned int currentSize=initialSize;
 197         unsigned int previousSize=initialSize;
 198
 199         // Lower the threshold till explosion
 200         while ( (previousSize<m_MinimumSize) || (  (currentLower> m_MinimumLowerThreshold) && ( (unsigned int) currentSize <= (unsigned int )m_Multiplier* previousSize)  ) )
 201           {
 202             // Update threshold
 203             currentLower-=m_ThresholdStepSize;
 204             function->ThresholdBetween(currentLower, m_Upper);
 205
 206             // Get current region size
 207             previousSize=currentSize;
 208             currentSize=0;
 209             IteratorType it = IteratorType ( outputImage, function, m_Seeds );
 210             it.GoToBegin();
 211             while( !it.IsAtEnd())
 212               {
 213                 ++currentSize;
 214                 ++it;
 215               }
 216             if (m_Verbose)  std::cout<<"Decreasing lower threshold to "<<currentLower<<", size is "<<currentSize  <<"(previously "<<previousSize<<") ..."<<std::endl;
 217           }
 218
 219
 220         // Explosion occured, increase lower theshold
 221         if ( (double)currentSize > m_Multiplier*(double) previousSize)
 222           {
 223             // Raise lower threshold
 224             if (m_Verbose)  std::cout<<"Explosion detected, adapting lower threshold ..."<<std::endl;
 225             currentLower+=m_ThresholdStepSize;
 226             InputImagePixelType currentStepSize=m_ThresholdStepSize;
 227
 228             while (currentStepSize>m_MinimumThresholdStepSize)
 229               {
 230                 currentStepSize/=2;
 231                 currentLower-=currentStepSize;
 232
 233                 // Get current region size
 234                 currentSize=0;
 235                 IteratorType it = IteratorType ( outputImage, function, m_Seeds );
 236                 it.GoToBegin();
 237                 while( !it.IsAtEnd())
 238                   {
 239                     ++currentSize;
 240                     ++it;
 241                   }
 242
 243                 if (m_Verbose)  std::cout<<"Adapting lower threshold to "<<currentLower<<", size is "<<currentSize  <<"(previously "<<previousSize<<") ..."<<std::endl;
 244
 245                 // Explosion: go back
 246                 if ((double)currentSize > m_Multiplier* (double) previousSize)
 247                   currentLower+=currentStepSize;
 248
 249                 // No explosion: adapt previous
 250                 else previousSize=currentSize;
 251               }
 252
 253             // Converged
 254             m_FinalLower=currentLower;
 255             if (m_Verbose)  std::cout<<"Final lower threshold  (precision="<<m_MinimumThresholdStepSize<<") is set to "<<m_FinalLower<<"..."<<std::endl;
 256           }
 257         else {
 258           if (m_Verbose) std::cout<<"No explosion before minimum lower threshold reached!"<<std::endl;
 259         }
 260
 261       } // end update lower
 262
 263
 264     //--------------------------------------------------
 265     // Increase upper threshold
 266     //--------------------------------------------------
 267     m_FinalUpper=m_Upper;
 268     if (m_AdaptUpperBorder)
 269       {
 270         // Init
 271         InputImagePixelType currentUpper=m_Upper;
 272         unsigned int currentSize=initialSize;
 273         unsigned int previousSize=initialSize;
 274
 275         // Upper the threshold till explosion
 276         while ((previousSize<m_MinimumSize) || ( (currentUpper< m_MaximumUpperThreshold) && ((double)currentSize <= m_Multiplier* (double) previousSize) )  )
 277           {
 278             // Update threshold
 279             currentUpper+=m_ThresholdStepSize;
 280             function->ThresholdBetween(m_Lower,currentUpper);
 281
 282             // Get current region size
 283             previousSize=currentSize;
 284             currentSize=0;
 285             IteratorType it = IteratorType ( outputImage, function, m_Seeds );
 286             it.GoToBegin();
 287             while( !it.IsAtEnd())
 288               {
 289                 ++currentSize;
 290                 ++it;
 291               }
 292             if (m_Verbose)  std::cout<<"Increasing upper threshold to "<<currentUpper<<", size is "<<currentSize  <<"(previously "<<previousSize<<") ..."<<std::endl;
 293
 294           }
 295
 296         // Explosion occured, decrease upper theshold
 297         if ((double)currentSize > m_Multiplier* (double) previousSize)
 298           {
 299             // Lower upper threshold
 300             if (m_Verbose)  std::cout<<"Explosion detected, adapting upper threshold ..."<<std::endl;
 301             currentUpper-=m_ThresholdStepSize;
 302             InputImagePixelType currentStepSize=m_ThresholdStepSize;
 303
 304             while (currentStepSize>m_MinimumThresholdStepSize)
 305               {
 306                 currentStepSize/=2;
 307                 currentUpper+=currentStepSize;
 308
 309                 // Get current region size
 310                 currentSize=0;
 311                 function->ThresholdBetween(m_Lower,currentUpper);
 312                 IteratorType it = IteratorType ( outputImage, function, m_Seeds );
 313                 it.GoToBegin();
 314                 while( !it.IsAtEnd())
 315                   {
 316                     ++currentSize;
 317                     ++it;
 318                   }
 319                 if (m_Verbose)  std::cout<<"Adapting upper threshold to "<<currentUpper<<", size is "<<currentSize  <<"(previously "<<previousSize<<") ..."<<std::endl;
 320
 321                 // Explosion: go back
 322                 if ((double)currentSize > m_Multiplier* (double) previousSize)
 323                   currentUpper-=currentStepSize;
 324
 325                 // No explosion: adapt previous
 326                 else previousSize=currentSize;
 327               }
 328
 329             // Converged
 330             m_FinalUpper=currentUpper;
 331             if (m_Verbose)  std::cout<<"Final upper threshold  (precision="<<m_MinimumThresholdStepSize<<") is set to "<<m_FinalUpper<<"..."<<std::endl;
 332
 333           }
 334         else {  if (m_Verbose)  std::cout<<"No explosion before maximum upper threshold reached!"<<std::endl; }
 335
 336       } // end update upper
 337
 338
 339
 340         //--------------------------------------------------
 341         // Update the output with final thresholds
 342         //--------------------------------------------------
 343     function->ThresholdBetween(m_FinalLower, m_FinalUpper);
 344     IteratorType it2 = IteratorType ( outputImage, function, m_Seeds );
 345     it2.GoToBegin();
 346     while( !it2.IsAtEnd())
 347       {
 348         it2.Set(m_ReplaceValue);
 349         ++it2;
 350       }
 351
 352   }
 353   //--------------------------------------------------------------------
 354
 355
 356 } // end namespace clitk
 357
 358 #endif