Re-implemented median filtering as a function.

bob/ip/base/auxiliary.cpp

@@ -1,5 +1,5 @@
 /**
- * @author Manuel Guenther <manuel.guenthr@idiap.ch>
+ * @author Manuel Guenther <manuel.guenther@idiap.ch>
  * @date Wed Jun 25 18:28:03 CEST 2014
  *
  * @brief Binds auxiliary functions of bob::ip::base class to python

bob/ip/base/cpp/Median.h

@@ -0,0 +1,69 @@
+/**
+ * @date Thu Jul  3 12:37:19 CEST 2014
+ * @author Manuel Guenther <manuel.guenther@idiap.ch>
+ *
+ * @brief This file provides a function to perform median filtering
+ *
+ * Copyright (C) Idiap Research Institute, Martigny, Switzerland
+ */
+
+#ifndef BOB_IP_BASE_MEDIAN_H
+#define BOB_IP_BASE_MEDIAN_H
+
+#include <vector>
+#include <algorithm>
+#include "bob/core/assert.h"
+#include "bob/core/cast.h"
+
+namespace bob { namespace ip { namespace base {
+
+  template <typename T>
+  void medianFilter(
+    const blitz::Array<T,2>& src,
+    blitz::Array<T,2>& dst,
+    const blitz::TinyVector<int,2>& radius
+  ){
+    // Checks
+    bob::core::array::assertZeroBase(src);
+    bob::core::array::assertZeroBase(dst);
+    blitz::TinyVector<int,2> dst_size(src.extent(0) - 2 * radius[0], src.extent(1) - 2 * radius[1]);
+    bob::core::array::assertSameShape(dst, dst_size);
+
+    // compute centeral pixel
+    int center = (2*radius[0]+1)*(2*radius[1]+1)/2;
+    // we only sort the first half of the sequence (this is all we need)
+    std::vector<T> _temp(center+1);
+    // iterate over the destination array
+    for (int y = 0; y < dst_size[0]; ++y)
+      for (int x = 0; x < dst_size[1]; ++x){
+        // get a slice from the src array
+        const blitz::Array<T,2> slice(src(blitz::Range(y, y + 2 * radius[0]), blitz::Range(x, x + 2 * radius[1])));
+        // compute the median
+        // we only sort the first half of the sequence
+        std::partial_sort_copy(slice.begin(), slice.end(), _temp.begin(), _temp.end());
+        // get the central element
+        dst(y,x) = _temp[center];
+    }
+  }
+
+
+  template <typename T>
+  void medianFilter(
+    const blitz::Array<T,3>& src,
+    blitz::Array<T,3>& dst,
+    const blitz::TinyVector<int,2>& radius
+  ){
+    // iterate over the color layers
+    for (int p = 0; p < dst.extent(0); ++p){
+      const blitz::Array<T,2> src_slice = src(p, blitz::Range::all(), blitz::Range::all());
+      blitz::Array<T,2> dst_slice = dst(p, blitz::Range::all(), blitz::Range::all());
+
+      // Apply median filter to the plane
+      medianFilter(src_slice, dst_slice, radius);
+    }
+  }
+
+} } } // namespaces
+
+#endif // BOB_IP_BASE_MEDIAN_H
+

bob/ip/base/filter.cpp

@@ -0,0 +1,83 @@
+/**
+ * @author Manuel Guenther <manuel.guenther@idiap.ch>
+ * @date Thu Jul  3 13:30:38 CEST 2014
+ *
+ * @brief Binds image filter functions of bob::ip::base class to python
+ *
+ * Copyright (C) 2011-2014 Idiap Research Institute, Martigny, Switzerland
+ */
+
+
+#include "main.h"
+#include "cpp/Median.h"
+
+static inline bool f(PyObject* o){return o != 0 && PyObject_IsTrue(o) > 0;}  /* converts PyObject to bool and returns false if object is NULL */
+
+
+bob::extension::FunctionDoc s_median = bob::extension::FunctionDoc(
+  "median",
+  "Performs a median filtering of the input image with the given radius",
+  "This function computes a histogram of the given input image, in several ways.\n\n"
+  "* (version 1 and 2, only valid for uint8 and uint16 types -- and uint32 and uint64 when ``bin_count`` is specified or ``hist`` is given as parameter): For each pixel value of the ``src`` image, a histogram bin is computed, using a fast implementation. "
+  "The number of bins can be limited, and there will be a check that the source image pixels are actually in the desired range ``(0, bin_count-1)``\n\n"
+  "* (version 3 and 4, valid for many data types): The histogram is computed by defining regular bins between the provided minimum and maximum values."
+)
+.add_prototype("src, radius, [dst]", "dst")
+.add_parameter("src", "array_like (2D or 3D)", "The source image to filter, might be a gray level image or a color image")
+.add_parameter("radius", "(int, int)", "The radius of the median filter; the final filter will have the size ``(2*radius[0]+1, 2*radius[1]+1)``")
+.add_parameter("dst", "array_like (2D or 3D)", "The median-filtered image to write; need to be of size ``src.shape - 2*radius``; if not specified, it will be created")
+.add_return("dst", "array_like (2D or 3D)", "The median-filtered image; the same as the ``dst`` parameter, if specified")
+;
+
+template <typename T, int D> PyObject* inner_median(PyBlitzArrayObject* src, PyBlitzArrayObject* dst, const blitz::TinyVector<int,2>& radius) {
+  bob::ip::base::medianFilter(*PyBlitzArrayCxx_AsBlitz<T, D>(src), *PyBlitzArrayCxx_AsBlitz<T, D>(dst), radius);
+  Py_INCREF(dst);
+  return PyBlitzArray_AsNumpyArray(dst, 0);
+}
+
+PyObject* PyBobIpBase_median(PyObject*, PyObject* args, PyObject* kwargs) {
+  TRY
+
+  static char* kwlist[] = {c("src"), c("radius"), c("dst"), NULL};
+
+  PyBlitzArrayObject* src,* dst = 0;
+  blitz::TinyVector<int,2> radius;
+
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&(ii)|O&", kwlist, &PyBlitzArray_Converter, &src, &radius[0], &radius[1], &PyBlitzArray_OutputConverter, &dst)) return 0;
+
+  auto src_ = make_safe(src), dst_ = make_xsafe(dst);
+
+  // allocate the output, if needed
+  if (!dst){
+    if (src->ndim == 2){
+      Py_ssize_t n[] = {src->shape[0] - 2*radius[0], src->shape[1] - 2*radius[1]};
+      dst = reinterpret_cast<PyBlitzArrayObject*>(PyBlitzArray_SimpleNew(src->type_num, 2, n));
+    } else if (src->ndim == 3){
+      Py_ssize_t n[] = {src->shape[0], src->shape[1] - 2*radius[0], src->shape[2] - 2*radius[1]};
+      dst = reinterpret_cast<PyBlitzArrayObject*>(PyBlitzArray_SimpleNew(src->type_num, 3, n));
+    } else {
+      PyErr_Format(PyExc_TypeError, "'median' : only 2D or 3D arrays are supported.");
+      return 0;
+    }
+    dst_ = make_safe(dst);
+  } else {
+    if (dst->type_num != src->type_num || dst->ndim != src->ndim){
+      PyErr_Format(PyExc_TypeError, "'median' : 'src' and 'dst' images must have the same type and number of dimensions, but %s != %s or %d != %d.", PyBlitzArray_TypenumAsString(src->type_num), PyBlitzArray_TypenumAsString(dst->type_num), (int)src->ndim, (int)dst->ndim);
+      return 0;
+    }
+  }
+
+  // compute the median
+  switch (src->type_num){
+    case NPY_UINT8:   if (src->ndim == 2) return inner_median<uint8_t,2>(src, dst, radius);  else return inner_median<uint8_t,3>(src, dst, radius);
+    case NPY_UINT16:  if (src->ndim == 2) return inner_median<uint16_t,2>(src, dst, radius); else return inner_median<uint16_t,3>(src, dst, radius);
+    case NPY_FLOAT16: if (src->ndim == 2) return inner_median<double,2>(src, dst, radius);   else return inner_median<double,3>(src, dst, radius);
+    default:
+      PyErr_Format(PyExc_ValueError, "'median' of %s arrays is currently not supported, only uint8, uint16 or float64 arrays are", PyBlitzArray_TypenumAsString(src->type_num));
+      return 0;
+  }
+
+  CATCH_("in median", 0)
+}
+
+

bob/ip/base/main.cpp

@@ -101,6 +101,12 @@ static PyMethodDef module_methods[] = {
     METH_VARARGS|METH_KEYWORDS,
     s_zigzag.doc()
   },
+  {
+    s_median.name(),
+    (PyCFunction)PyBobIpBase_median,
+    METH_VARARGS|METH_KEYWORDS,
+    s_median.doc()
+  },
   {0}  // Sentinel
 };
 

bob/ip/base/main.h

@@ -170,5 +170,8 @@ extern bob::extension::FunctionDoc s_histogramEqualization;
 PyObject* PyBobIpBase_zigzag(PyObject*, PyObject*, PyObject*);
 extern bob::extension::FunctionDoc s_zigzag;
 
+// filtering
+PyObject* PyBobIpBase_median(PyObject*, PyObject*, PyObject*);
+extern bob::extension::FunctionDoc s_median;
 
 #endif // BOB_IP_BASE_MAIN_H

bob/ip/base/old/main.cc

@@ -23,7 +23,6 @@ void bind_ip_gaussian_scale_space();
 void bind_ip_wgaussian();
 void bind_ip_msr();
 void bind_ip_sqi();
-void bind_ip_median();
 void bind_ip_sobel();
 void bind_ip_hog();
 void bind_ip_glcm_uint8();
@@ -55,7 +54,6 @@ BOOST_PYTHON_MODULE(_old_library) {
   bind_ip_wgaussian();
   bind_ip_msr();
   bind_ip_sqi();
-  bind_ip_median();
   bind_ip_sobel();
   bind_ip_hog();
   bind_ip_glcm_uint8();

bob/ip/base/test/test_filter.py

@@ -0,0 +1,35 @@
+#!/usr/bin/env python
+# vim: set fileencoding=utf-8 :
+# Manuel Guenther <manuel.guenther@idiap.ch>
+# Thu Jul  3 14:31:48 CEST 2014
+#
+# Copyright (C) 2011-2014 Idiap Research Institute, Martigny, Switzerland
+
+"""Tests filtertering"""
+
+import numpy
+import nose.tools
+import bob.ip.base
+import bob.sp
+
+import bob.io.base
+import bob.io.base.test_utils
+import bob.io.image
+
+def test_median():
+  # tests median filtering
+  src = numpy.array([
+      [1, 2, 3, 4, 5],
+      [6, 7, 8, 9, 10],
+      [11, 12, 13, 14, 15],
+      [16, 17, 18, 19, 20]],
+      dtype = numpy.uint16
+  )
+  ref = numpy.array([
+    [7, 8, 9],
+    [12, 13, 14]],
+    dtype = numpy.uint16)
+
+  dst = bob.ip.base.median(src, (1,1))
+  assert numpy.allclose(ref, dst)
+

setup.py

@@ -61,7 +61,6 @@
           "bob/ip/base/old/GLCM.cc",
           "bob/ip/base/old/GLCMProp.cc",
           "bob/ip/base/old/HOG.cc",
-          "bob/ip/base/old/Median.cc",
           "bob/ip/base/old/MultiscaleRetinex.cc",
           "bob/ip/base/old/SelfQuotientImage.cc",
           "bob/ip/base/old/shear.cc",
@@ -105,6 +104,7 @@
           "bob/ip/base/lbp_top.cpp",
           "bob/ip/base/dct_features.cpp",
           "bob/ip/base/tan_triggs.cpp",
+          "bob/ip/base/filter.cpp",
           "bob/ip/base/utils.cpp",
           "bob/ip/base/main.cpp",
           ],