From 223315bcfefd3078cb0fd5140271cc4d12539a7d Mon Sep 17 00:00:00 2001 From: Manuel Guenther Date: Thu, 3 Jul 2014 14:41:55 +0200 Subject: [PATCH] Re-implemented median filtering as a function. --- bob/ip/base/auxiliary.cpp | 2 +- bob/ip/base/cpp/Median.h | 69 +++++++++++++++++++++++++++ bob/ip/base/filter.cpp | 83 +++++++++++++++++++++++++++++++++ bob/ip/base/main.cpp | 6 +++ bob/ip/base/main.h | 3 ++ bob/ip/base/old/main.cc | 2 - bob/ip/base/test/test_filter.py | 35 ++++++++++++++ setup.py | 2 +- 8 files changed, 198 insertions(+), 4 deletions(-) create mode 100644 bob/ip/base/cpp/Median.h create mode 100644 bob/ip/base/filter.cpp create mode 100644 bob/ip/base/test/test_filter.py diff --git a/bob/ip/base/auxiliary.cpp b/bob/ip/base/auxiliary.cpp index e41000e..3c3ed68 100644 --- a/bob/ip/base/auxiliary.cpp +++ b/bob/ip/base/auxiliary.cpp @@ -1,5 +1,5 @@ /** - * @author Manuel Guenther + * @author Manuel Guenther * @date Wed Jun 25 18:28:03 CEST 2014 * * @brief Binds auxiliary functions of bob::ip::base class to python diff --git a/bob/ip/base/cpp/Median.h b/bob/ip/base/cpp/Median.h new file mode 100644 index 0000000..8a24cb1 --- /dev/null +++ b/bob/ip/base/cpp/Median.h @@ -0,0 +1,69 @@ +/** + * @date Thu Jul 3 12:37:19 CEST 2014 + * @author Manuel Guenther + * + * @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 +#include +#include "bob/core/assert.h" +#include "bob/core/cast.h" + +namespace bob { namespace ip { namespace base { + + template + void medianFilter( + const blitz::Array& src, + blitz::Array& dst, + const blitz::TinyVector& radius + ){ + // Checks + bob::core::array::assertZeroBase(src); + bob::core::array::assertZeroBase(dst); + blitz::TinyVector 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 _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 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 + void medianFilter( + const blitz::Array& src, + blitz::Array& dst, + const blitz::TinyVector& radius + ){ + // iterate over the color layers + for (int p = 0; p < dst.extent(0); ++p){ + const blitz::Array src_slice = src(p, blitz::Range::all(), blitz::Range::all()); + blitz::Array 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 + diff --git a/bob/ip/base/filter.cpp b/bob/ip/base/filter.cpp new file mode 100644 index 0000000..b6f9f0b --- /dev/null +++ b/bob/ip/base/filter.cpp @@ -0,0 +1,83 @@ +/** + * @author Manuel Guenther + * @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 PyObject* inner_median(PyBlitzArrayObject* src, PyBlitzArrayObject* dst, const blitz::TinyVector& radius) { + bob::ip::base::medianFilter(*PyBlitzArrayCxx_AsBlitz(src), *PyBlitzArrayCxx_AsBlitz(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 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(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(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(src, dst, radius); else return inner_median(src, dst, radius); + case NPY_UINT16: if (src->ndim == 2) return inner_median(src, dst, radius); else return inner_median(src, dst, radius); + case NPY_FLOAT16: if (src->ndim == 2) return inner_median(src, dst, radius); else return inner_median(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) +} + + diff --git a/bob/ip/base/main.cpp b/bob/ip/base/main.cpp index d82dfb6..b0e55ec 100644 --- a/bob/ip/base/main.cpp +++ b/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 }; diff --git a/bob/ip/base/main.h b/bob/ip/base/main.h index 1d7b093..41bd83a 100644 --- a/bob/ip/base/main.h +++ b/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 diff --git a/bob/ip/base/old/main.cc b/bob/ip/base/old/main.cc index 69a5071..b1b0039 100644 --- a/bob/ip/base/old/main.cc +++ b/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(); diff --git a/bob/ip/base/test/test_filter.py b/bob/ip/base/test/test_filter.py new file mode 100644 index 0000000..c2f368f --- /dev/null +++ b/bob/ip/base/test/test_filter.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python +# vim: set fileencoding=utf-8 : +# Manuel Guenther +# 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) + diff --git a/setup.py b/setup.py index 73a051c..ba7cfe6 100644 --- a/setup.py +++ b/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", ],