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

itkRankHistogram.h

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkRankHistogram.h,v $
00005   Language:  C++
00006   Date:      $Date: 2009-05-13 21:52:25 $
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 __itkRankHistogram_h
00020 #define __itkRankHistogram_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 #include <sstream>
00034 
00035 template <class TInputPixel>
00036 class RankHistogram
00037 {
00038 public:
00039   RankHistogram() 
00040     {
00041     m_Rank = 0.5;
00042     }
00043   virtual ~RankHistogram(){}
00044 
00045   virtual RankHistogram * Clone() const 
00046     {
00047     return 0;
00048     }
00049    
00050   virtual void AddPixel(const TInputPixel &itkNotUsed(p)){}
00051 
00052   virtual void RemovePixel(const TInputPixel &itkNotUsed(p)){}
00053  
00054   void AddBoundary(){}
00055 
00056   void RemoveBoundary(){}
00057  
00058   virtual TInputPixel GetValue( const TInputPixel & ){return 0;}
00059 
00060   void SetRank(float rank)
00061     {
00062     m_Rank = rank;
00063     }
00064 
00065 protected:
00066   float m_Rank;
00067 };
00068 
00069 template <class TInputPixel, class TCompare>
00070 class RankHistogramMap : public RankHistogram<TInputPixel>
00071 {
00072 public:
00073 
00074   typedef RankHistogram<TInputPixel>       Superclass;
00075 
00076 public:
00077   RankHistogramMap() 
00078     {
00079     m_Below = m_Entries = 0;
00080     // can't set m_RankIt until something has been put in the histogram
00081     m_Initialized = false;
00082     if( m_Compare( NumericTraits< TInputPixel >::max(), 
00083                    NumericTraits< TInputPixel >::NonpositiveMin() ) )
00084       {
00085       m_InitVal = NumericTraits< TInputPixel >::NonpositiveMin();
00086       }
00087     else
00088       {
00089       m_InitVal = NumericTraits< TInputPixel >::max();
00090       }
00091     m_RankValue = m_InitVal;
00092     m_RankIt = m_Map.begin();  // equivalent to setting to the intial value
00093     }
00094   ~RankHistogramMap()
00095     {
00096     }
00097 
00098   void AddPixel(const TInputPixel &p)
00099     {
00100     m_Map[ p ]++;
00101     ++m_Entries;
00102     if (!m_Initialized)
00103       {
00104       m_Initialized = true;
00105       m_RankIt = m_Map.begin();
00106       m_RankValue = p;
00107       }
00108     if (m_Compare(p, m_RankValue) || p == m_RankValue)
00109       {
00110       ++m_Below;
00111       }
00112 
00113     }
00114 
00115   void RemovePixel(const TInputPixel &p)
00116     {
00117     m_Map[ p ]--; 
00118     if (m_Compare(p, m_RankValue) || p == m_RankValue)
00119       {
00120       --m_Below;
00121       }
00122     --m_Entries;
00123     }
00124  
00125   void Initialize()
00126     {
00127     m_RankIt = m_Map.begin();
00128     }
00129 
00130   TInputPixel GetValue( const TInputPixel & )
00131     {
00132     unsigned long target = (int)(this->m_Rank * (m_Entries-1)) + 1;
00133     unsigned long total = m_Below;
00134     unsigned long ThisBin;
00135     bool eraseFlag = false;
00136 
00137     // an itkAssertOrThrowMacro is better than a log message in that case
00138     itkAssertOrThrowMacro( m_Initialized, "Not Initialized" );
00139 
00140     if (total < target)
00141       {
00142       typename MapType::iterator searchIt = m_RankIt;
00143       typename MapType::iterator eraseIt;
00144       
00145       while (searchIt != m_Map.end())
00146         {
00147         // cleaning up the map - probably a better way of organising
00148         // the loop. Currently makes sure that the search iterator is
00149         // incremented before deleting
00150         ++searchIt;
00151         ThisBin = searchIt->second;
00152         total += ThisBin;
00153         if (eraseFlag)
00154           {
00155           m_Map.erase(eraseIt);
00156           eraseFlag = false;
00157           }
00158         if (ThisBin <= 0)
00159           {
00160           eraseFlag = true;
00161           eraseIt = searchIt;
00162           }
00163         if (total >= target)
00164           break;
00165         }
00166       m_RankValue = searchIt->first;
00167       m_RankIt = searchIt;
00168       }
00169     else
00170       {
00171       typename MapType::iterator searchIt = m_RankIt;
00172       typename MapType::iterator eraseIt;
00173 
00174       while(searchIt != m_Map.begin())
00175         {
00176         ThisBin = searchIt->second;
00177         unsigned int tbelow = total - ThisBin;
00178         if (tbelow < target) // we've overshot
00179           break;
00180         if (eraseFlag)
00181           {
00182           m_Map.erase(eraseIt);
00183           eraseFlag = false;
00184           }
00185         if (ThisBin <= 0)
00186           {
00187           eraseIt = searchIt;
00188           eraseFlag = true;
00189           }
00190         total = tbelow;
00191 //         std::cout << searchIt->first << std::endl;
00192 
00193         --searchIt;
00194         }
00195       m_RankValue = searchIt->first;
00196       m_RankIt = searchIt;
00197       }
00198 
00199     m_Below = total;
00200 
00201     return(m_RankValue);
00202 
00203     }
00204 
00205   Superclass * Clone() const
00206     {
00207     RankHistogramMap *result = new RankHistogramMap();
00208     result->m_Map = this->m_Map;
00209     result->m_Rank = this->m_Rank;
00210     result->m_Below = this->m_Below;
00211     result->m_Entries = this->m_Entries;
00212     result->m_InitVal = this->m_InitVal;
00213     result->m_RankValue = this->m_RankValue;
00214     result->m_Initialized = this->m_Initialized;
00215     if (result->m_Initialized)
00216       result->m_RankIt = result->m_Map.find(this->m_RankValue);
00217     return(result);
00218     }
00219 
00220 private:
00221   typedef typename std::map< TInputPixel, unsigned long, TCompare > MapType;
00222   
00223   MapType       m_Map;
00224   unsigned long m_Below;
00225   unsigned long m_Entries;
00226   TInputPixel   m_RankValue;
00227   TInputPixel   m_InitVal;
00228   TCompare      m_Compare;
00229   bool          m_Initialized;
00230   // This iterator will point at the desired rank value
00231   typename MapType::iterator m_RankIt;
00232 };
00233 
00234 template <class TInputPixel, class TCompare>
00235 class RankHistogramVec : public RankHistogram<TInputPixel>
00236 {
00237 public:
00238   typedef RankHistogram<TInputPixel>       Superclass;
00239 
00240 public:
00241   RankHistogramVec() 
00242     {
00243     m_Size = static_cast<unsigned int>( NumericTraits< TInputPixel >::max() - 
00244                                         NumericTraits< TInputPixel >::NonpositiveMin() + 1 );
00245     m_Vec.resize(m_Size, 0 );
00246     if( m_Compare( NumericTraits< TInputPixel >::max(), 
00247                    NumericTraits< TInputPixel >::NonpositiveMin() ) )
00248       {
00249       m_InitVal = NumericTraits< TInputPixel >::NonpositiveMin();
00250       }
00251     else
00252       {
00253       m_InitVal = NumericTraits< TInputPixel >::max();
00254       }
00255     m_Entries = m_Below = 0;
00256     m_RankValue = m_InitVal  - NumericTraits< TInputPixel >::NonpositiveMin();
00257     }
00258 
00259   RankHistogramVec(bool NoInit) 
00260     {
00261     m_Size = static_cast<unsigned int>( NumericTraits< TInputPixel >::max() - 
00262                                         NumericTraits< TInputPixel >::NonpositiveMin() + 1 );
00263     if (!NoInit) m_Vec.resize(m_Size, 0 );
00264     if( m_Compare( NumericTraits< TInputPixel >::max(), 
00265                    NumericTraits< TInputPixel >::NonpositiveMin() ) )
00266       {
00267       m_InitVal = NumericTraits< TInputPixel >::NonpositiveMin();
00268       }
00269     else
00270       {
00271       m_InitVal = NumericTraits< TInputPixel >::max();
00272       }
00273     m_Entries = m_Below = 0;
00274     m_RankValue = m_InitVal  - NumericTraits< TInputPixel >::NonpositiveMin();
00275     }
00276 
00277 
00278   ~RankHistogramVec()
00279     {
00280     }
00281 
00282   TInputPixel GetValue( const TInputPixel & )
00283     {
00284     unsigned long target = (int)(this->m_Rank * (m_Entries-1)) + 1;
00285     unsigned long total = m_Below;
00286     long unsigned int pos = (long unsigned int)(m_RankValue - NumericTraits< TInputPixel >::NonpositiveMin()); 
00287 
00288     if (total < target)
00289       {
00290       while (pos < m_Size)
00291         {
00292         ++pos;
00293         total += m_Vec[pos];
00294         if (total >= target)
00295           break;
00296         }
00297       }
00298     else
00299       {
00300       while(pos > 0)
00301         {
00302         unsigned int tbelow = total - m_Vec[pos];
00303         if (tbelow < target) // we've overshot
00304           break;
00305         total = tbelow;
00306         --pos;
00307         }
00308       }
00309 
00310     m_RankValue = (TInputPixel)(pos + NumericTraits< TInputPixel >::NonpositiveMin());
00311     m_Below = total;
00312     return(m_RankValue);
00313     }
00314 
00315   void AddPixel(const TInputPixel &p)
00316     {
00317     long unsigned int idx = (long unsigned int)(p - NumericTraits< TInputPixel >::NonpositiveMin());
00318     m_Vec[ idx  ]++; 
00319     if (m_Compare(p, m_RankValue) || p == m_RankValue)
00320       {
00321       ++m_Below;
00322       }
00323     ++m_Entries;
00324     }
00325 
00326   void RemovePixel(const TInputPixel &p)
00327     {
00328     itkAssertOrThrowMacro( ( p - NumericTraits< TInputPixel >::NonpositiveMin() >= 0 ),
00329       "pixel value too close to zero");
00330     itkAssertOrThrowMacro( ( p - NumericTraits< TInputPixel >::NonpositiveMin() < (int)m_Vec.size()),
00331       "pixel value outside the range of m_Vec.size()");
00332     itkAssertOrThrowMacro( ( m_Entries >= 1), "Not enough entries");
00333     m_Vec[ (long unsigned int)(p - NumericTraits< TInputPixel >::NonpositiveMin())  ]--; 
00334     --m_Entries;
00335 
00336     if (m_Compare(p, m_RankValue) || p == m_RankValue)
00337       {
00338       --m_Below;
00339       }
00340     }
00341  
00342   Superclass * Clone() const
00343     {
00344     RankHistogramVec *result = new RankHistogramVec(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 private:
00357   typedef typename std::vector<unsigned long> VecType;
00358   
00359   VecType      m_Vec;
00360   unsigned int m_Size;
00361   TCompare     m_Compare;
00362   //unsigned int m_CurrentValue;
00363   //TInputPixel m_CurrentValue;
00364   TInputPixel  m_RankValue;
00365   TInputPixel  m_InitVal;
00366   int          m_Below;
00367   int          m_Entries;
00368 };
00369 
00370 } // end namespace itk
00371 #endif
00372 

Generated at Tue Sep 15 04:30:37 2009 for ITK by doxygen 1.5.8 written by Dimitri van Heesch, © 1997-2000