ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkRankHistogram.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 // histogram from the moving histogram operations
00019 #ifndef __itkRankHistogram_h
00020 #define __itkRankHistogram_h
00021 
00022 #include "itkIntTypes.h"
00023 #include "itkNumericTraits.h"
00024 
00025 #include <map>
00026 #include <vector>
00027 
00028 namespace itk
00029 {
00030 namespace Function
00031 {
00032 // a simple histogram class hierarchy. One subclass will be maps, the
00033 // other vectors.
00034 // This version is intended for keeping track of arbitary ranks. It is
00035 // based on the code from consolidatedMorphology.
00036 //
00037 // Support for different TCompare hasn't been tested, and shouldn't be
00038 // necessary for the rank filters.
00039 //
00040 // This is a modified version for use with masks. Need to allow for
00041 // the situation in which the map is empty
00042 
00043 /*
00044  *
00045  * This code was contributed in the Insight Journal paper:
00046  * "Efficient implementation of kernel filtering"
00047  * by Beare R., Lehmann G
00048  * http://hdl.handle.net/1926/555
00049  * http://www.insight-journal.org/browse/publication/160
00050  *
00051  */
00052 
00053 
00054 template< class TInputPixel >
00055 class RankHistogram
00056 {
00057 public:
00058 
00059   typedef std::less< TInputPixel > TCompare;
00060 
00061   RankHistogram()
00062   {
00063     m_Rank = 0.5;
00064     m_Below = m_Entries = 0;
00065     // can't set m_RankIt until something has been put in the histogram
00066     m_Initialized = false;
00067     if ( m_Compare( NumericTraits< TInputPixel >::max(),
00068                     NumericTraits< TInputPixel >::NonpositiveMin() ) )
00069       {
00070       m_InitVal = NumericTraits< TInputPixel >::max();
00071       }
00072     else
00073       {
00074       m_InitVal = NumericTraits< TInputPixel >::NonpositiveMin();
00075       }
00076     m_RankValue = m_InitVal;
00077     m_RankIt = m_Map.begin();  // equivalent to setting to the intial value
00078   }
00079 
00080   ~RankHistogram()
00081   {}
00082 
00083   RankHistogram & operator=( const RankHistogram & hist )
00084   {
00085     m_Map = hist.m_Map;
00086     m_Rank = hist.m_Rank;
00087     m_Below = hist.m_Below;
00088     m_Entries = hist.m_Entries;
00089     m_InitVal = hist.m_InitVal;
00090     m_RankValue = hist.m_RankValue;
00091     m_Initialized = hist.m_Initialized;
00092     if ( m_Initialized )
00093       {
00094       m_RankIt = m_Map.find(m_RankValue);
00095       }
00096     return *this;
00097   }
00098 
00099   void AddPixel(const TInputPixel & p)
00100   {
00101     m_Map[p]++;
00102     if ( !m_Initialized )
00103       {
00104       m_Initialized = true;
00105       m_RankIt = m_Map.begin();
00106       m_Entries = m_Below = 0;
00107       m_RankValue = p;
00108       }
00109     if ( m_Compare(p, m_RankValue) || p == m_RankValue )
00110       {
00111       ++m_Below;
00112       }
00113     ++m_Entries;
00114   }
00115 
00116   void RemovePixel(const TInputPixel & p)
00117   {
00118     m_Map[p]--;
00119     if ( m_Compare(p, m_RankValue) || p == m_RankValue )
00120       {
00121       --m_Below;
00122       }
00123     --m_Entries;
00124     // this is the change that makes this version less efficient. The
00125     // simplest approach I can think of with maps, though
00126     if ( m_Entries <= 0 )
00127       {
00128       m_Initialized = false;
00129       m_Below = 0;
00130       m_Map.clear();
00131       }
00132   }
00133 
00134   bool IsValid()
00135   {
00136     return m_Initialized;
00137   }
00138 
00139   TInputPixel GetValueBruteForce()
00140   {
00141     SizeValueType count = 0;
00142     SizeValueType target = (int)( m_Rank * ( m_Entries - 1 ) ) + 1;
00143     for( typename MapType::iterator it=m_Map.begin(); it != m_Map.end(); it++ )
00144       {
00145       count += it->second;
00146       if( count >= target )
00147         {
00148         return it->first;
00149         }
00150       }
00151     return NumericTraits< TInputPixel >::max();
00152   }
00153 
00154   TInputPixel GetValue(const TInputPixel &)
00155   {
00156     SizeValueType target = (SizeValueType)( m_Rank * ( m_Entries - 1 ) ) + 1;
00157     SizeValueType total = m_Below;
00158     SizeValueType ThisBin;
00159     bool          eraseFlag = false;
00160 
00161     if ( total < target )
00162       {
00163       typename MapType::iterator searchIt = m_RankIt;
00164       typename MapType::iterator eraseIt;
00165 
00166       while ( searchIt != m_Map.end() )
00167         {
00168         // cleaning up the map - probably a better way of organising
00169         // the loop. Currently makes sure that the search iterator is
00170         // incremented before deleting
00171         ++searchIt;
00172         ThisBin = searchIt->second;
00173         total += ThisBin;
00174         if ( eraseFlag )
00175           {
00176           m_Map.erase(eraseIt);
00177           eraseFlag = false;
00178           }
00179         if ( ThisBin <= 0 )
00180           {
00181           eraseFlag = true;
00182           eraseIt = searchIt;
00183           }
00184         if ( total >= target )
00185           {
00186           break;
00187           }
00188         }
00189       m_RankValue = searchIt->first;
00190       m_RankIt = searchIt;
00191       }
00192     else
00193       {
00194       typename MapType::iterator searchIt = m_RankIt;
00195       typename MapType::iterator eraseIt;
00196 
00197       while ( searchIt != m_Map.begin() )
00198         {
00199         ThisBin = searchIt->second;
00200         unsigned int tbelow = total - ThisBin;
00201         if ( tbelow < target ) // we've overshot
00202           {
00203           break;
00204           }
00205         if ( eraseFlag )
00206           {
00207           m_Map.erase(eraseIt);
00208           eraseFlag = false;
00209           }
00210         if ( ThisBin <= 0 )
00211           {
00212           eraseIt = searchIt;
00213           eraseFlag = true;
00214           }
00215         total = tbelow;
00216 
00217         --searchIt;
00218         }
00219       m_RankValue = searchIt->first;
00220       m_RankIt = searchIt;
00221       }
00222 
00223     m_Below = total;
00224     itkAssertInDebugAndIgnoreInReleaseMacro( m_RankValue == GetValueBruteForce() );
00225     return ( m_RankValue );
00226   }
00227 
00228   void SetRank(float rank)
00229   {
00230     m_Rank = rank;
00231   }
00232 
00233   void AddBoundary(){}
00234 
00235   void RemoveBoundary(){}
00236 
00237   static bool UseVectorBasedAlgorithm()
00238   {
00239     return false;
00240   }
00241 
00242 protected:
00243   float m_Rank;
00244 
00245 private:
00246   typedef typename std::map< TInputPixel, SizeValueType, TCompare > MapType;
00247 
00248   MapType       m_Map;
00249   SizeValueType m_Below;
00250   SizeValueType m_Entries;
00251   TInputPixel   m_RankValue;
00252   TInputPixel   m_InitVal;
00253   TCompare      m_Compare;
00254   bool          m_Initialized;
00255   // This iterator will point at the desired rank value
00256   typename MapType::iterator m_RankIt;
00257 };
00258 
00259 
00260 template< class TInputPixel >
00261 class VectorRankHistogram
00262 {
00263 public:
00264   typedef std::less< TInputPixel > TCompare;
00265 
00266   VectorRankHistogram()
00267   {
00268     m_Size = (OffsetValueType)NumericTraits< TInputPixel >::max() - (OffsetValueType)NumericTraits< TInputPixel >::NonpositiveMin() + 1;
00269     m_Vec.resize(m_Size, 0);
00270     if ( m_Compare( NumericTraits< TInputPixel >::max(),
00271                     NumericTraits< TInputPixel >::NonpositiveMin() ) )
00272       {
00273       m_InitVal = NumericTraits< TInputPixel >::NonpositiveMin();
00274       }
00275     else
00276       {
00277       m_InitVal = NumericTraits< TInputPixel >::max();
00278       }
00279     m_Entries = m_Below = 0;
00280     m_RankValue = m_InitVal  - NumericTraits< TInputPixel >::NonpositiveMin();
00281     m_Rank = 0.5;
00282   }
00283 
00284   ~VectorRankHistogram() {}
00285 
00286   bool IsValid()
00287   {
00288     return m_Entries > 0;
00289   }
00290 
00291   TInputPixel GetValueBruteForce()
00292   {
00293     SizeValueType count = 0;
00294     SizeValueType target = (SizeValueType)( m_Rank * ( m_Entries - 1 ) ) + 1;
00295     for( SizeValueType i=0; i<m_Size; i++ )
00296       {
00297       count += m_Vec[i];
00298       if( count >= target )
00299         {
00300         return i + NumericTraits< TInputPixel >::NonpositiveMin();
00301         }
00302       }
00303     return NumericTraits< TInputPixel >::max();
00304   }
00305 
00306   TInputPixel GetValue(const TInputPixel &)
00307   {
00308     return GetValueBruteForce();
00309   }
00310 
00311   void AddPixel(const TInputPixel & p)
00312   {
00313     OffsetValueType q = (OffsetValueType)p - NumericTraits< TInputPixel >::NonpositiveMin();
00314 
00315     m_Vec[q]++;
00316     if ( m_Compare(p, m_RankValue) || p == m_RankValue )
00317       {
00318       ++m_Below;
00319       }
00320     ++m_Entries;
00321   }
00322 
00323   void RemovePixel(const TInputPixel & p)
00324   {
00325     const OffsetValueType q = (OffsetValueType)p - NumericTraits< TInputPixel >::NonpositiveMin();
00326 
00327     itkAssertInDebugAndIgnoreInReleaseMacro( q >= 0 );
00328     itkAssertInDebugAndIgnoreInReleaseMacro( q < (int)m_Vec.size() );
00329     itkAssertInDebugAndIgnoreInReleaseMacro( m_Entries >= 1 );
00330     itkAssertInDebugAndIgnoreInReleaseMacro( m_Vec[q] > 0 );
00331 
00332     m_Vec[q]--;
00333     --m_Entries;
00334 
00335     if ( m_Compare(p, m_RankValue) || p == m_RankValue )
00336       {
00337       --m_Below;
00338       }
00339   }
00340 
00341   void SetRank(float rank)
00342   {
00343     m_Rank = rank;
00344   }
00345 
00346   void AddBoundary(){}
00347 
00348   void RemoveBoundary(){}
00349 
00350   static bool UseVectorBasedAlgorithm()
00351   {
00352     return true;
00353   }
00354 
00355 protected:
00356   float m_Rank;
00357 
00358 private:
00359   typedef typename std::vector< SizeValueType > VecType;
00360 
00361   VecType       m_Vec;
00362   SizeValueType m_Size;
00363   TCompare      m_Compare;
00364   TInputPixel   m_RankValue;
00365   TInputPixel   m_InitVal;
00366   int           m_Below;
00367   int           m_Entries;
00368 };
00369 
00370 // now create MorphologicalGradientHistogram specilizations using the VectorMorphologicalGradientHistogram
00371 // as base class
00372 
00375 template<>
00376 class RankHistogram<unsigned char>:
00377   public VectorRankHistogram<unsigned char>
00378 {
00379 };
00380 
00381 template<>
00382 class RankHistogram<signed char>:
00383   public VectorRankHistogram<signed char>
00384 {
00385 };
00386 
00387 template<>
00388 class RankHistogram<bool>:
00389   public VectorRankHistogram<bool>
00390 {
00391 };
00392 
00395 } // end namespace Function
00396 } // end namespace itk
00397 #endif
00398