itk/clitkExplosionControlledThresholdConnectedImageFilter.txx

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