/******************************************************************** ** Image Component Library (ICL) ** ** ** ** Copyright (C) 2006-2013 CITEC, University of Bielefeld ** ** Neuroinformatics Group ** ** Website: www.iclcv.org and ** ** http://opensource.cit-ec.de/projects/icl ** ** ** ** File : ICLFilter/src/ICLFilter/MedianOp.h ** ** Module : ICLFilter ** ** Authors: Christof Elbrechter, Robert Haschke, Andre Justus, ** ** Sergius Gaulik ** ** ** ** ** ** GNU LESSER GENERAL PUBLIC LICENSE ** ** This file may be used under the terms of the GNU Lesser General ** ** Public License version 3.0 as published by the ** ** ** ** Free Software Foundation and appearing in the file LICENSE.LGPL ** ** included in the packaging of this file. Please review the ** ** following information to ensure the license requirements will ** ** be met: http://www.gnu.org/licenses/lgpl-3.0.txt ** ** ** ** The development of this software was supported by the ** ** Excellence Cluster EXC 277 Cognitive Interaction Technology. ** ** The Excellence Cluster EXC 277 is a grant of the Deutsche ** ** Forschungsgemeinschaft (DFG) in the context of the German ** ** Excellence Initiative. ** ** ** ********************************************************************/ #pragma once #include #include namespace icl { namespace filter{ /// Class that provides median filter abilities \ingroup UNARY \ingroup NBH /** The median class provides the ability for arbitrary mask sizes. Although the functionality of the IPP-optimized implementation and the fallback C++-implementation is identical, the performances are of different orders of magnitude. Fallback-functions for masks with dimensions 3x3 and 5x5 are optimizied with SSE-instructions, but the computation time for images of types Img8u and Img16s is still greater than the one of IPP-implementations. However SSE-implementations show a better performance for the type icl32f. Here the algorithm just takes the median using min and max functions. Images of the type Img8u and Img16s with other mask sizes are processed by the algorithm, that was introduced by Thomas S. Huang. With the help of a histogram the algorithm runs in O(n). For all other types a trivial implementation was designed to calculate the median values. It uses the naive algorithm of sorting all N pixel values inside the median mask, and setting the destination pixel value to the mid-element of the sorted pixel list. This algorithm runs in O(w*h*N*log(N)) where (w,h) is the size of source images ROI, and N=n*n is the mask size used. The following code extract explains the operation of the fallback algorithm in Img-style notation for a single channel image, and Img8u type (the real implementation uses some special optimizations, that are not mentioned further): 
        void channel_median_8u(Img8u &src, Img8u &dst, int w, int h, int c)
        {
            std::vector list;
            
            for(Img8u::iterator s=src.begin(c), d=dst.begin(c); s.inRegion() ; s++, d++ )
            {
                for(Img8u::iterator sR(s,w,h); sR.inRegion(); sR++)
                {
                   list.push_back(*sR);
                }
                std::sort(list.begin(),list.end());
                *d = list[w*h/2];
                list.clear();
            }
        }
The IPP implementation uses a fast median algorithm that estimates the median not by sorting a list, but by working on a histogram. Look at the IPPI manual for more detail.

IPP for floats

Currently the IPP implementation of 2D-median filtering for Ipp32f type is slower than the C++-implementation, so the C++-fallback is always used then.

Mask-Sizes

Although the fallback C++ implementation can work with arbitrary mask sizes, the Median will internally use odd mask dimension like 3x3 or 5x7. If an even width or height parameter is given to the Median constructor, the next higher odd value is used.

Benchmarks

table

3x3 5x5 51x51 51x3 3x51
icl8u, ipp ~1.8ms ~2.6ms ~292ms ~32.4ms ~328ms
icl8u, c++ ~2.6ms ~11.7ms ~323ms ~350ms ~47ms
icl16s, ipp ~2.7ms ~4ms ~21100ms ~338ms ~420ms
icl16s, c++ ~4.2ms ~16.2ms ~370ms ~384ms ~90ms
float, ipp ~25ms ~257ms ~?ms ~52000ms ~49000ms
float, c++ ~11.6ms ~39ms ~217000ms ~8700ms ~8700ms

Details

Test A 1000x1000 icl8u-image with IPP

- mask size 3x3 ~1.8ms (highly optimized) - mask size 5x5 ~2.6ms - mask size 51x51 ~292ms (still usable!) - mask size 51x3 ~32.4ms - mask size 3x51 ~328ms (mask height specifies time usage)

Test B 1000x1000 icl8u-image no IPP

- mask size 3x3 ~2.6ms (less than 1.5 times slower) - mask size 5x5 ~11.7ms (about 4 times slower) - mask size 51x51 ~323ms (IPP is not much faster) - mask size 51x3 ~350ms (mask width specifies time usage) - mask size 3x51 ~47ms

Test A 1000x1000 icl8u-image with IPP

- mask size 3x3 ~2.7ms (highly optimized) - mask size 5x5 ~4ms - mask size 51x51 ~21100ms (unusable!) - mask size 51x3 ~338ms - mask size 3x51 ~420ms

Test B 1000x1000 icl16s-image no IPP

- mask size 3x3 ~4.2ms (1.5 times slower) - mask size 5x5 ~16.2ms ( 4 times slower) - mask size 51x51 ~370ms (IPP is not much faster) - mask size 51x3 ~384ms (mask width specifies time usage) - mask size 3x51 ~90ms

Test C 1000x1000 icl32f-image with IPP

- mask size 3x3 ~25ms - mask size 5x5 ~257ms - mask size 51x51 ~?ms (calculation time is to great) - mask size 51x3 ~52000ms (unusable!) - mask size 3x51 ~49000ms (unusable!)

Test C 1000x1000 icl32f-image no IPP

- mask size 3x3 ~11.6ms (sse optimization) - mask size 5x5 ~39ms (sse optimization) - mask size 51x51 ~217000ms (unusable!) - mask size 51x3 ~8700ms - mask size 3x51 ~8700ms

Example

Here is an examle, how to use the Median object. 
      // create source and destination image
      Img8u src(640,480,3), *poDst=0;
    
      // acquire some image data
      ...
    
      // create the median object
      Median m(Size(5,5));
  
      // apply the median on the images - first call (slow)
      // destination image is renewed to 640x480x3 (memory allocation)
      m.apply(&src,&poDst);
  
      // enter iteration loop 
      while(1)
      {
          // aquire some new image data
          ...       
  
          // apply the median filter (now fast, as no memory
          // allocation must be performed)
          m.apply(&src,&poDst);
  
          // further processing steps
          ...
      }
*/ class ICLFilter_API MedianOp : public NeighborhoodOp { public: /// Constructor that creates a median filter object, with specified mask size /** @param maskSize of odd width and height Even width or height is increased to next higher odd value. **/ MedianOp (const utils::Size &maskSize):NeighborhoodOp(adaptSize(maskSize)){} /// applies the median operation on poSrc and stores the result in poDst /** The depth, channel count and size of poDst is adapted to poSrc' ROI: @param poSrc source image @param ppoDst pointer to destination image **/ void apply(const core::ImgBase *poSrc, core::ImgBase **ppoDst); /// Import unaryOps apply function without destination image using NeighborhoodOp::apply; /// ensures that mask width and height are odd /** This is a workaround, necessary because of an ipp Bug that allows no even mask sizes here! @param size size to adjust @return adjusted size (width and height are rounded up to the next higher odd value **/ virtual utils::Size adaptSize(const utils::Size &size){ return utils::Size(1+ 2*(size.width/2),1+ 2*(size.height/2)); } }; } // namespace filter }