Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkMaskedRankHistogram.h

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkMaskedRankHistogram.h,v $
00005   Language:  C++
00006   Date:      $Date: 2009-05-13 22:17:41 $
00007   Version:   $Revision: 1.5 $
00008 
00009   Copyright (c) Insight Software Consortium. All rights reserved.
00010   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
00011 
00012     This software is distributed WITHOUT ANY WARRANTY; without even 
00013     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
00014     PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 
00018 // histogram from the moving histogram operations
00019 #ifndef __itkMaskedRankHistogram_h
00020 #define __itkMaskedRankHistogram_h
00021 #include "itkNumericTraits.h"
00022 
00023 namespace itk {
00024 
00025 // a simple histogram class hierarchy. One subclass will be maps, the
00026 // other vectors.
00027 // This version is intended for keeping track of arbitary ranks. It is
00028 // based on the code from consolidatedMorphology.
00029 //
00030 // Support for different TCompare hasn't been tested, and shouldn't be
00031 // necessary for the rank filters.
00032 //
00033 // This is a modified version for use with masks. Need to allow for
00034 // the situation in which the map is empty
00035 template <class TInputPixel>
00036 class MaskedRankHistogram
00037 {
00038 public:
00039   MaskedRankHistogram() 
00040     {
00041     m_Rank = 0.5;
00042     }
00043   virtual ~MaskedRankHistogram(){}
00044 
00045   virtual MaskedRankHistogram * Clone() const { return 0; }
00046   
00047   virtual void Reset(){}
00048     
00049   virtual void AddPixel(const TInputPixel & itkNotUsed(p)){}
00050 
00051   virtual void RemovePixel(const TInputPixel & itkNotUsed(p)){}
00052  
00053   // For the map based version - to be called after there is some data
00054   // included. Meant to be an optimization so that the rank value
00055   // iterator is properly set.
00056   // virtual void Initialize(){};
00057 
00058   virtual bool IsValid()
00059     {
00060     return false;
00061     }
00062 
00063   virtual TInputPixel GetValue( const TInputPixel & )
00064     {
00065     return NumericTraits<TInputPixel>::Zero;
00066     }
00067 
00068   void SetRank(float rank)
00069     {
00070     m_Rank = rank;
00071     }
00072 
00073   void AddBoundary(){}
00074 
00075   void RemoveBoundary(){}
00076  
00077 protected:
00078   float m_Rank;
00079 
00080 };
00081 
00082 template <class TInputPixel, class TCompare>
00083 class MaskedRankHistogramMap : public MaskedRankHistogram<TInputPixel>
00084 {
00085 public:
00086 
00087   typedef MaskedRankHistogram<TInputPixel>  Superclass;
00088 
00089 public:
00090   MaskedRankHistogramMap() 
00091     {
00092     m_Below = m_Entries = 0;
00093     // can't set m_RankIt until something has been put in the histogram
00094     m_Initialized = false;
00095     if( m_Compare( NumericTraits< TInputPixel >::max(), 
00096                    NumericTraits< TInputPixel >::NonpositiveMin() ) )
00097       {
00098       m_InitVal = NumericTraits< TInputPixel >::max();
00099       }
00100     else
00101       {
00102       m_InitVal = NumericTraits< TInputPixel >::NonpositiveMin();
00103       }
00104     m_RankValue = m_InitVal;
00105     m_RankIt = m_Map.begin();  // equivalent to setting to the intial value
00106      }
00107   ~MaskedRankHistogramMap()
00108     {
00109     }
00110 
00111   Superclass * Clone() const
00112     {
00113     MaskedRankHistogramMap *result = new MaskedRankHistogramMap();
00114     result->m_Map = this->m_Map;
00115     result->m_Rank = this->m_Rank;
00116     result->m_Below = this->m_Below;
00117     result->m_Entries = this->m_Entries;
00118     result->m_InitVal = this->m_InitVal;
00119     result->m_RankValue = this->m_RankValue;
00120     result->m_Initialized = this->m_Initialized;
00121     if (result->m_Initialized)
00122       result->m_RankIt = result->m_Map.find(this->m_RankValue);
00123     return(result);
00124     }
00125   void Reset()
00126     {
00127     m_Map.clear();
00128     m_RankValue = m_InitVal;
00129     m_Entries = m_Below = 0;
00130     }
00131     
00132   void AddPixel(const TInputPixel &p)
00133     {
00134     m_Map[ p ]++;
00135     if (!m_Initialized)
00136       {
00137       m_Initialized = true;
00138       m_RankIt = m_Map.begin();
00139       m_Entries = m_Below = 0;
00140       m_RankValue = p;
00141       }
00142     if (m_Compare(p, m_RankValue) || p == m_RankValue)
00143       {
00144       ++m_Below;
00145       }
00146     ++m_Entries;
00147     }
00148 
00149   void RemovePixel(const TInputPixel &p)
00150     {
00151     m_Map[ p ]--; 
00152     if (m_Compare(p, m_RankValue) || p == m_RankValue)
00153       {
00154       --m_Below;
00155       }
00156     --m_Entries;
00157     // this is the change that makes this version less efficient. The
00158     // simplest approach I can think of with maps, though
00159     if (m_Entries <= 0)
00160       {
00161       m_Initialized = false;
00162       m_Below = 0;
00163       m_Map.clear();
00164       }
00165     }
00166  
00167 //   void Initialize()
00168 //   {
00169 //     m_RankIt = m_Map.begin();
00170 //   }
00171 
00172   virtual bool IsValid()
00173     {
00174     return m_Initialized;
00175     }
00176 
00177 
00178   TInputPixel GetValue( const TInputPixel & )
00179     {
00180     unsigned long target = (int)(this->m_Rank * (m_Entries-1)) + 1;
00181     unsigned long total = m_Below;
00182     unsigned long ThisBin;
00183     bool eraseFlag = false;
00184 
00185     if (total < target)
00186       {
00187       typename MapType::iterator searchIt = m_RankIt;
00188       typename MapType::iterator eraseIt;
00189       
00190       while (searchIt != m_Map.end())
00191         {
00192         // cleaning up the map - probably a better way of organising
00193         // the loop. Currently makes sure that the search iterator is
00194         // incremented before deleting
00195         ++searchIt;
00196         ThisBin = searchIt->second;
00197         total += ThisBin;
00198         if (eraseFlag)
00199           {
00200           m_Map.erase(eraseIt);
00201           eraseFlag = false;
00202           }
00203         if (ThisBin <= 0)
00204           {
00205           eraseFlag = true;
00206           eraseIt = searchIt;
00207           }
00208         if (total >= target)
00209           break;
00210         }
00211       m_RankValue = searchIt->first;
00212       m_RankIt = searchIt;
00213       }
00214     else
00215       {
00216       typename MapType::iterator searchIt = m_RankIt;
00217       typename MapType::iterator eraseIt;
00218 
00219       while(searchIt != m_Map.begin())
00220         {
00221         ThisBin = searchIt->second;
00222         unsigned int tbelow = total - ThisBin;
00223         if (tbelow < target) // we've overshot
00224           break;
00225         if (eraseFlag)
00226           {
00227           m_Map.erase(eraseIt);
00228           eraseFlag = false;
00229           }
00230         if (ThisBin <= 0)
00231           {
00232           eraseIt = searchIt;
00233           eraseFlag = true;
00234           }
00235         total = tbelow;
00236         
00237         --searchIt;
00238         }
00239       m_RankValue = searchIt->first;
00240       m_RankIt = searchIt;
00241       }
00242 
00243     m_Below = total;
00244     return(m_RankValue);
00245 
00246     }
00247 
00248 #if 0
00249   TInputPixel GetValue( const TInputPixel & )
00250     {    // clean the map
00251     typename MapType::iterator mapIt = m_Map.begin();
00252     while( mapIt != m_Map.end() )
00253       {
00254       if( mapIt->second <= 0 )
00255         {
00256         // this value must be removed from the histogram
00257         // The value must be stored and the iterator updated before removing the value
00258         // or the iterator is invalidated.
00259         TInputPixel toErase = mapIt->first;
00260         mapIt++;
00261         m_Map.erase( toErase );
00262         }
00263       else
00264         {
00265         mapIt++;
00266         // don't remove all the zero value found, just remove the one before the current maximum value
00267         // the histogram may become quite big on real type image, but it's an important increase of performances
00268         break;
00269         }
00270       }
00271     
00272     // and return the value
00273     return m_Map.begin()->first;
00274     }
00275 #endif
00276 
00277 private:
00278   typedef typename std::map< TInputPixel, unsigned long, TCompare > MapType;
00279   
00280   MapType       m_Map;
00281   unsigned long m_Below;
00282   unsigned long m_Entries;
00283   TInputPixel   m_RankValue;
00284   TInputPixel   m_InitVal;
00285   TCompare      m_Compare;
00286   bool          m_Initialized;
00287   // This iterator will point at the desired rank value
00288   typename MapType::iterator m_RankIt;
00289 
00290 };
00291 
00292 template <class TInputPixel, class TCompare>
00293 class MaskedRankHistogramVec : public MaskedRankHistogram<TInputPixel>
00294 {
00295 public:
00296 
00297   typedef MaskedRankHistogram<TInputPixel>  Superclass;
00298 
00299 public:
00300 
00301   MaskedRankHistogramVec() 
00302     {
00303     m_Size = static_cast<unsigned int>( NumericTraits< TInputPixel >::max() - 
00304                                         NumericTraits< TInputPixel >::NonpositiveMin() + 1 );
00305     m_Vec.resize(m_Size, 0 );
00306     if( m_Compare( NumericTraits< TInputPixel >::max(), 
00307                    NumericTraits< TInputPixel >::NonpositiveMin() ) )
00308       {
00309       m_InitVal = NumericTraits< TInputPixel >::NonpositiveMin();
00310       }
00311     else
00312       {
00313       m_InitVal = NumericTraits< TInputPixel >::max();
00314       }
00315     m_Entries = m_Below = 0;
00316     m_RankValue = m_InitVal  - NumericTraits< TInputPixel >::NonpositiveMin();
00317     }
00318 
00319   MaskedRankHistogramVec(bool NoInit) 
00320     {
00321     m_Size = static_cast<unsigned int>( NumericTraits< TInputPixel >::max() - 
00322                                         NumericTraits< TInputPixel >::NonpositiveMin() + 1 );
00323     if (!NoInit) m_Vec.resize(m_Size, 0 );
00324     if( m_Compare( NumericTraits< TInputPixel >::max(), 
00325                    NumericTraits< TInputPixel >::NonpositiveMin() ) )
00326       {
00327       m_InitVal = NumericTraits< TInputPixel >::NonpositiveMin();
00328       }
00329     else
00330       {
00331       m_InitVal = NumericTraits< TInputPixel >::max();
00332       }
00333     m_Entries = m_Below = 0;
00334     m_RankValue = m_InitVal  - NumericTraits< TInputPixel >::NonpositiveMin();
00335     }
00336 
00337 
00338   ~MaskedRankHistogramVec()
00339     {
00340     }
00341 
00342   Superclass * Clone() const
00343     {
00344     MaskedRankHistogramVec *result = new MaskedRankHistogramVec(true);
00345     result->m_Vec = this->m_Vec;
00346     result->m_Size = this->m_Size;
00347     //result->m_CurrentValue = this->m_CurrentValue;
00348     result->m_InitVal = this->m_InitVal;
00349     result->m_Entries = this->m_Entries;
00350     result->m_Below = this->m_Below;
00351     result->m_Rank = this->m_Rank;
00352     result->m_RankValue = this->m_RankValue;
00353     return(result);
00354     }
00355 
00356   virtual bool IsValid()
00357     {
00358     return m_Entries > 0;
00359     }
00360 
00361 
00362   TInputPixel GetValue( const TInputPixel & )
00363     {
00364     unsigned long target = (int)(this->m_Rank * (m_Entries-1)) + 1;
00365     unsigned long total = m_Below;
00366     long unsigned int pos = (long unsigned int)(m_RankValue - NumericTraits< TInputPixel >::NonpositiveMin()); 
00367 
00368     if (total < target)
00369       {
00370       while (pos < m_Size)
00371         {
00372         ++pos;
00373         total += m_Vec[pos];
00374         if (total >= target)
00375           break;
00376         }
00377       }
00378     else
00379       {
00380       while(pos > 0)
00381         {
00382         unsigned int tbelow = total - m_Vec[pos];
00383         if (tbelow < target) // we've overshot
00384           break;
00385         total = tbelow;
00386         --pos;
00387         }
00388       }
00389 
00390     m_RankValue = (TInputPixel)(pos + NumericTraits< TInputPixel >::NonpositiveMin());
00391     m_Below = total;
00392     return(m_RankValue);
00393     }
00394 
00395   void Reset()
00396     {
00397     std::fill(&(m_Vec[0]), &(m_Vec[m_Size-1]),0);
00398     m_RankValue = m_InitVal;
00399     m_Entries = m_Below = 0;
00400     }
00401     
00402   void AddPixel(const TInputPixel &p)
00403     {
00404     long unsigned int idx = (long unsigned int)(p - NumericTraits< TInputPixel >::NonpositiveMin());
00405     m_Vec[ idx  ]++; 
00406     if (m_Compare(p, m_RankValue) || p == m_RankValue)
00407       {
00408       ++m_Below;
00409       }
00410     ++m_Entries;
00411     }
00412 
00413   void RemovePixel(const TInputPixel &p)
00414     {
00415     const long int q = p - NumericTraits< TInputPixel >::NonpositiveMin();
00416     itkAssertOrThrowMacro( ( q >= 0 ), "Input pixel value is out of range");
00417     itkAssertOrThrowMacro( ( q < (int)m_Vec.size() ), "Input pixel value is out of range");
00418     itkAssertOrThrowMacro( (m_Entries >= 1), "Insufficient entries");
00419 
00420     m_Vec[ static_cast<long unsigned int>(q)  ]--; 
00421     --m_Entries;
00422 
00423     if (m_Compare(p, m_RankValue) || p == m_RankValue)
00424       {
00425       --m_Below;
00426       }
00427     }
00428 
00429 private:
00430   typedef typename std::vector<unsigned long> VecType;
00431   
00432   VecType      m_Vec;
00433   unsigned int m_Size;
00434   TCompare     m_Compare;
00435   TInputPixel  m_RankValue;
00436   TInputPixel  m_InitVal;
00437   int          m_Below;
00438   int          m_Entries;
00439 };
00440 
00441 } // end namespace itk
00442 #endif
00443 

Generated at Thu May 28 10:42:53 2009 for ITK by doxygen 1.5.5 written by Dimitri van Heesch, © 1997-2000