Commit bed3a7cf authored by Manuel Günther's avatar Manuel Günther

Moved Wiener filter from bob.learn.misc.WienerMachine to bob.ip.base.Wiener

parent 2bc9bac4
/**
* @date Fri Sep 30 16:56:06 2011 +0200
* @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
*
* @brief Implements a Wiener filter
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#include <bob.core/array_copy.h>
#include <bob.core/cast.h>
#include <bob.ip.base/Wiener.h>
#include <complex>
bob::ip::base::Wiener::Wiener(const blitz::Array<double,2>& Ps, const double Pn, const double variance_threshold)
: m_Ps(bob::core::array::ccopy(Ps)),
m_variance_threshold(variance_threshold),
m_Pn(Pn),
m_W(m_Ps.extent(0),m_Ps.extent(1)),
m_fft(m_Ps.extent(0),m_Ps.extent(1)),
m_ifft(m_Ps.extent(0),m_Ps.extent(1)),
m_buffer1(m_Ps.extent(0),m_Ps.extent(1)),
m_buffer2(m_Ps.extent(0),m_Ps.extent(1))
{
computeW();
}
bob::ip::base::Wiener::Wiener(const blitz::TinyVector<int,2>& size, const double Pn, const double variance_threshold)
: m_Ps(size),
m_variance_threshold(variance_threshold),
m_Pn(Pn),
m_W(size),
m_fft(size[0], size[1]),
m_ifft(size[0], size[1]),
m_buffer1(0,0), m_buffer2(0,0)
{
m_Ps = 1.;
computeW();
}
bob::ip::base::Wiener::Wiener(const blitz::Array<double,3>& data, const double variance_threshold)
: m_variance_threshold(variance_threshold)
{
train(data);
}
bob::ip::base::Wiener::Wiener(const bob::ip::base::Wiener& other)
: m_Ps(bob::core::array::ccopy(other.m_Ps)),
m_variance_threshold(other.m_variance_threshold),
m_Pn(other.m_Pn),
m_W(bob::core::array::ccopy(other.m_W)),
m_fft(other.m_fft),
m_ifft(other.m_ifft),
m_buffer1(m_Ps.extent(0),m_Ps.extent(1)),
m_buffer2(m_Ps.extent(0),m_Ps.extent(1))
{
}
bob::ip::base::Wiener::Wiener(bob::io::base::HDF5File& config){
load(config);
}
bob::ip::base::Wiener& bob::ip::base::Wiener::operator=(const bob::ip::base::Wiener& other){
if (this != &other)
{
m_Ps.reference(bob::core::array::ccopy(other.m_Ps));
m_Pn = other.m_Pn;
m_variance_threshold = other.m_variance_threshold;
m_W.reference(bob::core::array::ccopy(other.m_W));
m_fft.setShape(m_Ps.extent(0),m_Ps.extent(1));
m_ifft.setShape(m_Ps.extent(0),m_Ps.extent(1));
m_buffer1.resize(m_Ps.shape());
m_buffer2.resize(m_Ps.shape());
}
return *this;
}
bool bob::ip::base::Wiener::operator==(const bob::ip::base::Wiener& b) const{
return bob::core::array::isEqual(m_Ps, b.m_Ps) &&
m_variance_threshold == b.m_variance_threshold &&
m_Pn == b.m_Pn &&
bob::core::array::isEqual(m_W, b.m_W);
}
bool bob::ip::base::Wiener::operator!=(const bob::ip::base::Wiener& b) const{
return !(this->operator==(b));
}
bool bob::ip::base::Wiener::is_similar_to(const bob::ip::base::Wiener& b, const double r_epsilon, const double a_epsilon) const{
return bob::core::array::isClose(m_Ps, b.m_Ps, r_epsilon, a_epsilon) &&
bob::core::isClose(m_variance_threshold, b.m_variance_threshold, r_epsilon, a_epsilon) &&
bob::core::isClose(m_Pn, b.m_Pn, r_epsilon, a_epsilon) &&
bob::core::array::isClose(m_W, b.m_W, r_epsilon, a_epsilon);
}
void bob::ip::base::Wiener::load(bob::io::base::HDF5File& config){
//reads all data directly into the member variables
m_Ps.reference(config.readArray<double,2>("Ps"));
m_Pn = config.read<double>("Pn");
m_variance_threshold = config.read<double>("variance_threshold");
m_W.reference(config.readArray<double,2>("W"));
m_fft.setShape(m_Ps.extent(0),m_Ps.extent(1));
m_ifft.setShape(m_Ps.extent(0),m_Ps.extent(1));
m_buffer1.resize(m_Ps.shape());
m_buffer2.resize(m_Ps.shape());
}
void bob::ip::base::Wiener::resize(const blitz::TinyVector<int,2>& size){
m_Ps.resizeAndPreserve(size);
m_W.resizeAndPreserve(size);
m_fft.setShape(size[0], size[1]);
m_ifft.setShape(size[0], size[1]);
m_buffer1.resizeAndPreserve(size);
m_buffer2.resizeAndPreserve(size);
}
void bob::ip::base::Wiener::save(bob::io::base::HDF5File& config) const{
config.setArray("Ps", m_Ps);
config.set("Pn", m_Pn);
config.set("variance_threshold", m_variance_threshold);
config.setArray("W", m_W);
}
void bob::ip::base::Wiener::train(const blitz::Array<double,3>& ar){
// Data is checked now and conforms, just proceed w/o any further checks.
const size_t n_samples = ar.extent(0);
const size_t height = ar.extent(1);
const size_t width = ar.extent(2);
// resize with the new dimensions
resize(blitz::TinyVector<int,2>(height, width));
// Loads the data
blitz::Array<double,3> data(height, width, n_samples);
blitz::Array<std::complex<double>,2> sample_fft(height, width);
blitz::Range all = blitz::Range::all();
for (size_t i=0; i<n_samples; ++i) {
blitz::Array<double,2> sample = ar(i,all,all);
blitz::Array<std::complex<double>,2> sample_c = bob::core::array::cast<std::complex<double> >(sample);
m_fft(sample_c, sample_fft);
data(all,all,i) = blitz::abs(sample_fft);
}
// Computes the mean of the training data
blitz::Array<double,2> tmp(height,width);
blitz::thirdIndex k;
tmp = blitz::mean(data,k);
// Removes the mean from the data
for (size_t i=0; i<n_samples; ++i) {
data(all,all,i) -= tmp;
}
// Computes power of 2 values
data *= data;
// Sums to get the variance
tmp = blitz::sum(data,k) / n_samples;
// sets the Wiener filter with the results:
setPs(tmp);
}
void bob::ip::base::Wiener::computeW(){
// W = 1 / (1 + Pn / Ps_thresholded)
m_W = 1. / (1. + m_Pn / m_Ps);
}
void bob::ip::base::Wiener::filter_(const blitz::Array<double,2>& input, blitz::Array<double,2>& output) const{
m_fft(bob::core::array::cast<std::complex<double> >(input), m_buffer1);
m_buffer1 *= m_W;
m_ifft(m_buffer1, m_buffer2);
output = blitz::abs(m_buffer2);
}
void bob::ip::base::Wiener::filter(const blitz::Array<double,2>& input, blitz::Array<double,2>& output) const{
if (m_W.extent(0) != input.extent(0)) { //checks input
boost::format m("number of input rows (%d) is not compatible with internal weight matrix (%d)");
m % input.extent(0) % m_W.extent(0);
throw std::runtime_error(m.str());
}
if (m_W.extent(1) != input.extent(1)) { //checks input
boost::format m("number of input columns (%d) is not compatible with internal weight matrix (%d)");
m % input.extent(1) % m_W.extent(1);
throw std::runtime_error(m.str());
}
if (m_W.extent(0) != output.extent(0)) { //checks output
boost::format m("number of output rows (%d) is not compatible with internal weight matrix (%d)");
m % output.extent(0) % m_W.extent(0);
throw std::runtime_error(m.str());
}
if (m_W.extent(1) != output.extent(1)) { //checks output
boost::format m("number of output columns (%d) is not compatible with internal weight matrix (%d)");
m % output.extent(1) % m_W.extent(1);
throw std::runtime_error(m.str());
}
filter_(input, output);
}
void bob::ip::base::Wiener::setVarianceThreshold(const double variance_threshold){
m_variance_threshold = variance_threshold;
applyVarianceThreshold();
computeW();
}
void bob::ip::base::Wiener::setPs(const blitz::Array<double,2>& Ps){
if (m_Ps.extent(0) != Ps.extent(0)) {
boost::format m("number of rows (%d) for input `Ps' does not match the expected (internal) size (%d)");
m % Ps.extent(0) % m_Ps.extent(0);
throw std::runtime_error(m.str());
}
if (m_Ps.extent(1) != Ps.extent(1)) {
boost::format m("number of columns (%d) for input `Ps' does not match the expected (internal) size (%d)");
m % Ps.extent(1) % m_Ps.extent(1);
throw std::runtime_error(m.str());
}
m_Ps = Ps;
computeW();
}
void bob::ip::base::Wiener::applyVarianceThreshold(){
m_Ps = blitz::where(m_Ps < m_variance_threshold, m_variance_threshold, m_Ps);
}
/**
* @date Fri Sep 30 16:56:06 2011 +0200
* @author Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
*
* Copyright (C) Idiap Research Institute, Martigny, Switzerland
*/
#ifndef BOB_IP_BASE_WIENER_H
#define BOB_IP_BASE_WIENER_H
#include <blitz/array.h>
#include <complex>
#include <bob.io.base/HDF5File.h>
#include <bob.sp/FFT2D.h>
namespace bob { namespace ip { namespace base {
/**
* @brief A Wiener filter, which can be used to denoise a signal,
* by comparing with a statistical model of the noiseless signal.\n
*
* Reference:\n
* Computer Vision: Algorithms and Applications, Richard Szeliski
* (Part 3.4.3)
*/
class Wiener
{
public: //api
/**
* @brief Constructor, builds a new Wiener filter. Wiener filter is
* initialized with the given size, Ps being sets to the variance
* threshold.
*/
Wiener(const blitz::TinyVector<int,2>& size, const double Pn, const double variance_threshold=1e-8);
/**
* @brief Builds a new filter with the given variance estimate Ps and
* noise level Pn.
*/
Wiener(const blitz::Array<double,2>& Ps, const double Pn, const double variance_threshold=1e-8);
/**
* @brief Trains a new Wiener filter with the given data
*/
Wiener(const blitz::Array<double,3>& data, const double variance_threshold=1e-8);
/**
* @brief Copy constructor
*/
Wiener(const Wiener& other);
/**
* @brief Loads a Wiener filter from file
*/
Wiener(bob::io::base::HDF5File& config);
/**
* @brief Assignment operator
*/
Wiener& operator=(const Wiener& other);
/**
* @brief Equal to
*/
bool operator==(const Wiener& other) const;
/**
* @brief Not equal to
*/
bool operator!=(const Wiener& other) const;
/**
* @brief Is similar to
*/
bool is_similar_to(const Wiener& b, const double r_epsilon=1e-5, const double a_epsilon=1e-8) const;
/**
* @brief Loads filter from file
*/
void load(bob::io::base::HDF5File& config);
/**
* @brief Saves filter to file
*/
void save(bob::io::base::HDF5File& config) const;
/**
* @brief Trains the Wiener filter to perform the filtering.
*/
void train(const blitz::Array<double,3>& data);
/**
* @brief filters the given input image
*
* The input and output are NOT checked for compatibility each time. It
* is your responsibility to do it.
*/
void filter_(const blitz::Array<double,2>& input, blitz::Array<double,2>& output) const;
/**
* @brief filters the given input image
*
* The input and output are checked for compatibility each time the
* filtering method is applied.
*/
void filter(const blitz::Array<double,2>& input, blitz::Array<double,2>& output) const;
/**
* @brief Resizes the filter and preserves the data
*/
void resize(const blitz::TinyVector<int,2>& size);
/**
* @brief Returns the current variance Ps estimated at each frequency
*/
const blitz::Array<double, 2>& getPs() const {return m_Ps;}
/**
* @brief Returns the current variance threshold applied to Ps
*/
double getVarianceThreshold() const {return m_variance_threshold;}
/**
* @brief Returns the current noise level Pn
*/
double getPn() const {return m_Pn;}
/**
* @brief Returns the size of the filter/input
*/
blitz::TinyVector<int,2> getSize() const {return m_W.shape();}
/**
* @brief Returns the current Wiener filter (in the frequency domain).
*/
const blitz::Array<double, 2>& getW() const {return m_W;}
/**
* @brief Sets the current variance Ps estimated at each frequency.
* This will also update the Wiener filter, using thresholded values.
*/
void setPs(const blitz::Array<double,2>& Ps);
/**
* @brief Sets the current variance threshold to be used.
* This will also update the Wiener filter
*/
void setVarianceThreshold(const double variance_threshold);
/**
* @brief Sets the current noise level Pn to be considered.
* This will update the Wiener filter
*/
void setPn(const double Pn) {m_Pn = Pn; computeW();}
private: //representation
void computeW(); /// Compute the Wiener filter using Pn, Ps, etc.
void applyVarianceThreshold(); /// Apply variance flooring threshold
blitz::Array<double, 2> m_Ps; ///< variance at each frequency estimated empirically
double m_variance_threshold; ///< Threshold on Ps values when computing the Wiener filter
/// (to avoid division by zero)
double m_Pn; ///< variance of the noise
blitz::Array<double, 2> m_W; ///< Wiener filter in the frequency domain (W=1/(1+Pn/Ps))
bob::sp::FFT2D m_fft;
bob::sp::IFFT2D m_ifft;
mutable blitz::Array<std::complex<double>, 2> m_buffer1; ///< a buffer for speed
mutable blitz::Array<std::complex<double>, 2> m_buffer2; ///< a buffer for speed
};
} } } // namespaces
#endif /* BOB_IP_BASE_WIENER_H */
......@@ -158,6 +158,7 @@ static PyObject* create_module (void) {
if (!init_BobIpBaseSIFT(module)) return 0;
if (!init_BobIpBaseHOG(module)) return 0;
if (!init_BobIpBaseGLCM(module)) return 0;
if (!init_BobIpBaseWiener(module)) return 0;
#ifdef HAVE_VLFEAT
if (!init_BobIpBaseVLFEAT(module)) return 0;
......
......@@ -31,6 +31,7 @@
#include <bob.ip.base/GeomNorm.h>
#include <bob.ip.base/FaceEyesNorm.h>
#include <bob.ip.base/GLCM.h>
#include <bob.ip.base/Wiener.h>
/// inserts the given key, value pair into the given dictionaries
......@@ -264,6 +265,18 @@ int PyBobIpBaseGLCM_Check(PyObject* o);
bool init_BobIpBaseGLCM(PyObject* module);
// .. Wiener filter
typedef struct {
PyObject_HEAD
boost::shared_ptr<bob::ip::base::Wiener> cxx;
} PyBobIpBaseWienerObject;
extern PyTypeObject PyBobIpBaseWiener_Type;
int PyBobIpBaseWiener_Check(PyObject* o);
bool init_BobIpBaseWiener(PyObject* module);
// block
PyObject* PyBobIpBase_block(PyObject*, PyObject*, PyObject*);
extern bob::extension::FunctionDoc s_block;
......
#!/usr/bin/env python
# vim: set fileencoding=utf-8 :
# Laurent El Shafey <Laurent.El-Shafey@idiap.ch>
#
# Copyright (C) 2011-2014 Idiap Research Institute, Martigny, Switzerland
"""Tests the bob.ip.base.Wiener filter
"""
import os
import numpy
import tempfile
import nose.tools
import bob.sp
import bob.io.base
import bob.ip.base
def test_initialization():
# Getters/Setters
m = bob.ip.base.Wiener((5,4),0.5)
nose.tools.eq_(m.size, (5,4))
m.size = (5,6)
nose.tools.eq_(m.size, (5,6))
ps1 = 0.2 + numpy.fabs(numpy.random.randn(5,6))
ps2 = 0.2 + numpy.fabs(numpy.random.randn(5,6))
m.Ps = ps1
assert numpy.allclose(m.Ps, ps1)
m.Ps = ps2
assert numpy.allclose(m.Ps, ps2)
pn1 = 0.5
m.Pn = pn1
assert abs(m.Pn - pn1) < 1e-5
var_thd = 1e-5
m.variance_threshold = var_thd
assert abs(m.variance_threshold - var_thd) < 1e-5
# Comparison operators
m2 = bob.ip.base.Wiener(m)
assert m == m2
assert (m != m2) is False
m3 = bob.ip.base.Wiener(ps2, pn1)
m3.variance_threshold = var_thd
assert m == m3
assert (m != m3 ) is False
# Computation of the Wiener filter W
w_py = 1 / (1. + m.Pn / m.Ps)
assert numpy.allclose(m.w, w_py)
def test_load_save():
m = bob.ip.base.Wiener((5,4),0.5)
# Save and read from file
filename = str(tempfile.mkstemp(".hdf5")[1])
m.save(bob.io.base.HDF5File(filename, 'w'))
m_loaded = bob.ip.base.Wiener(bob.io.base.HDF5File(filename))
assert m == m_loaded
assert (m != m_loaded) is False
assert m.is_similar_to(m_loaded)
# Make them different
m_loaded.variance_threshold = 0.001
assert (m == m_loaded) is False
assert m != m_loaded
# Clean-up
os.unlink(filename)
def test_filter():
ps = 0.2 + numpy.fabs(numpy.random.randn(5,6))
pn = 0.5
m = bob.ip.base.Wiener(ps,pn)
# Python way
sample = numpy.random.randn(5,6)
sample_fft = bob.sp.fft(sample.astype(numpy.complex128))
w = m.w
sample_fft_filtered = sample_fft * m.w
sample_filtered_py = numpy.absolute(bob.sp.ifft(sample_fft_filtered))
# Bob c++ way
sample_filtered0 = m.filter(sample)
sample_filtered1 = m(sample)
sample_filtered2 = numpy.zeros((5,6),numpy.float64)
m.filter(sample, sample_filtered2)
sample_filtered3 = numpy.zeros((5,6),numpy.float64)
m.filter(sample, sample_filtered3)
sample_filtered4 = numpy.zeros((5,6),numpy.float64)
m(sample, sample_filtered4)
assert numpy.allclose(sample_filtered0, sample_filtered_py)
assert numpy.allclose(sample_filtered1, sample_filtered_py)
assert numpy.allclose(sample_filtered2, sample_filtered_py)
assert numpy.allclose(sample_filtered3, sample_filtered_py)
assert numpy.allclose(sample_filtered4, sample_filtered_py)
def test_train():
def train_wiener_ps(training_set):
# Python implementation
n_samples = training_set.shape[0]
height = training_set.shape[1]
width = training_set.shape[2]
training_fftabs = numpy.zeros((n_samples, height, width), dtype=numpy.float64)
for n in range(n_samples):
sample = (training_set[n,:,:]).astype(numpy.complex128)
training_fftabs[n,:,:] = numpy.absolute(bob.sp.fft(sample))
mean = numpy.mean(training_fftabs, axis=0)
for n in range(n_samples):
training_fftabs[n,:,:] -= mean
training_fftabs = training_fftabs * training_fftabs
var_ps = numpy.mean(training_fftabs, axis=0)
return var_ps
n_samples = 20
height = 5
width = 6
training_set = 0.2 + numpy.fabs(numpy.random.randn(n_samples, height, width))
# Python implementation
var_ps = train_wiener_ps(training_set)
# Bob C++ implementation (variant 1) + comparison against python one
m1 = bob.ip.base.Wiener(training_set)
assert numpy.allclose(var_ps, m1.Ps)
# Bob C++ implementation (variant 2) + comparison against python one
m2 = bob.ip.base.Wiener(training_set, 0.)
assert numpy.allclose(var_ps, m2.Ps)
This diff is collapsed.
......@@ -35,6 +35,7 @@ References
.. [Lowe2004] *D. Lowe*. **Distinctive Image Features from Scale-Invariant Keypoints,** In International Journal of Computer Vision, 2004.
.. [Dalal2005] *N. Dalal, B. Triggs*. **Histograms of Oriented Gradients for Human Detection,** In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, 2005.
.. [Haralick1973] *R. M. Haralick, K. Shanmugam, I. Dinstein*. **Textural Features for Image Classification,** In IEEE Transactions on Systems, Man and Cybernetics, vol. SMC-3, No. 6, p. 610-621, 1973.
.. [Szeliski2010] *Richard Szeliski*. **Computer Vision: Algorithms and Applications** (1st ed.). Springer-Verlag New York, USA, 2010.
Indices and tables
------------------
......
......@@ -21,6 +21,7 @@ Classes
bob.ip.base.TanTriggs
bob.ip.base.Gaussian
bob.ip.base.Wiener
bob.ip.base.MultiscaleRetinex
bob.ip.base.WeightedGaussian
bob.ip.base.SelfQuotientImage
......
......@@ -131,6 +131,7 @@ setup(
"bob/ip/base/cpp/SIFT.cpp",
"bob/ip/base/cpp/HOG.cpp",
"bob/ip/base/cpp/GLCM.cpp",
"bob/ip/base/cpp/Wiener.cpp",
],
packages = packages,
boost_modules = boost_modules,
......@@ -163,6 +164,7 @@ setup(
"bob/ip/base/glcm.cpp",
"bob/ip/base/filter.cpp",
"bob/ip/base/utils.cpp",
"bob/ip/base/wiener.cpp",
"bob/ip/base/main.cpp",
],
packages = packages,
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment