m_convolve.h

Go to the documentation of this file.

/* $Id: m_convolve.h 15035 2013-07-29 23:32:59Z predoehl $ */

/* {{{======================================================================= *
   |
   |  Copyright (c) 1994-2013 by members of the Interdisciplinary Visual
   |  Intelligence Laboratory.
   |
   |  Personal and educational use of this code is granted, provided that this
   |  header is kept intact, and that the authorship is not misrepresented, that
   |  its use is acknowledged in publications, and relevant papers are cited.
   |
   |  For other use contact kobus AT sista DOT arizona DOT edu.
   |
   |  Please note that the code in this file has not necessarily been adequately
   |  tested. Naturally, there is no guarantee of performance, support, or
   |  fitness for any particular task. Nonetheless, I am interested in hearing
   |  about problems that you encounter.
   |
 * ====================================================================== }}}*/

// vim: tabstop=4 shiftwidth=4 foldmethod=marker

#ifndef KJB_CPP_M_CPP_M_CONVOLVE_H
#define KJB_CPP_M_CPP_M_CONVOLVE_H

#include <m_cpp/m_matrix.h>

#ifdef KJB_HAVE_FFTW /* Include the header, if we can. */

/*
 * Their header has a built-in extern-C gadget, so we do not need one here.
 * However, we would like to wrap their symbols in a namespace.
 */
namespace FFTW {
#include <fftw3.h>
}

#else /* KJB_HAVE_FFTW not defined */
#warning "This code requires FFTW version 3 to work properly."
const int FFTW_ESTIMATE = 0;
const int FFTW_MEASURE = 0;
const int FFTW_PATIENT = 0;
const int FFTW_EXHAUSTIVE = 0;
const int FFTW_DESTROY_INPUT = 0;
#endif /* KJB_HAVE_FFTW */


#include <utility>
#include <boost/scoped_ptr.hpp>
#include <boost/shared_ptr.hpp>


namespace kjb
{

#ifdef KJB_HAVE_FFTW /* This block defines a speciality vector for FFTW */


// We don't need the end pointer currently, but future changes might demand it.
#define KJB_FFTW_VECTOR_NEEDS_END 0


template <typename T>
class FFTW_vector
{
    T* m_begin;
#if KJB_FFTW_VECTOR_NEEDS_END
    T* m_end;
#endif
    FFTW_vector<T>(const FFTW_vector<T>&); // teaser
    FFTW_vector<T>& operator=(const FFTW_vector<T>&); // teaser

public:
    FFTW_vector(size_t n = 0)
    :   m_begin(0)
    {
        if (n && 0 == (m_begin = (T*) FFTW::fftw_malloc(sizeof(T) * n)))
        {
            KJB_THROW(Resource_exhaustion);
        }
#if KJB_FFTW_VECTOR_NEEDS_END
        m_end = m_begin + n;
#endif
    }

    ~FFTW_vector()
    {
        if (m_begin) FFTW::fftw_free(m_begin);
    }

    T* begin()
    {
        return m_begin;
    }

    const T* begin() const
    {
        return m_begin;
    }

#if KJB_FFTW_VECTOR_NEEDS_END
    T* end()
    {
        return m_end;
    }

    const T* end() const
    {
        return m_end;
    }
#endif

    void swap(FFTW_vector<T>& v)
    {
        using std::swap;
        swap(m_begin, v.m_begin);
#if KJB_FFTW_VECTOR_NEEDS_END
        swap(m_end, v.m_end);
#endif
    }
};


typedef FFTW_vector<FFTW::fftw_complex> FFTW_complex_vector;

typedef FFTW_vector<double> FFTW_real_vector;



#else /* what to do if FFTW library is not available? */
typedef void* FFTW_complex_vector;
typedef void* FFTW_real_vector;
#endif






template <class T, class U> class FFTW_Plan2d {};



#ifdef KJB_HAVE_FFTW /* This block defines FFTW plan wrappers. */

// "forward" FFT plan management
template<>
class FFTW_Plan2d<double, FFTW::fftw_complex>
{
    FFTW::fftw_plan plan;

    FFTW_Plan2d<double, FFTW::fftw_complex>(
            const FFTW_Plan2d<double, FFTW::fftw_complex>&); // teaser
    FFTW_Plan2d<double, FFTW::fftw_complex>& operator=(
            const FFTW_Plan2d<double, FFTW::fftw_complex>&); // teaser
public:
    FFTW_Plan2d<double, FFTW::fftw_complex>(
        int n0,
        int n1,
        double* in,
        FFTW::fftw_complex* out,
        unsigned flags
    )
    :   plan(FFTW::fftw_plan_dft_r2c_2d(n0, n1, in, out, flags))
    {
        if (0 == plan) KJB_THROW(Resource_exhaustion);
    }

    ~FFTW_Plan2d<double, FFTW::fftw_complex>()
    {
        FFTW::fftw_destroy_plan(plan);
    }

    void execute() const
    {
        FFTW::fftw_execute(plan);
    }

    void new_array_exec(double *in, FFTW::fftw_complex *out) const
    {
        FFTW::fftw_execute_dft_r2c(plan, in, out);
    }
};


// inverse-FFT plan management
template<>
class FFTW_Plan2d<FFTW::fftw_complex, double>
{
    FFTW::fftw_plan plan;

    FFTW_Plan2d<FFTW::fftw_complex, double>(
            const FFTW_Plan2d<FFTW::fftw_complex, double>&); // teaser
    FFTW_Plan2d<FFTW::fftw_complex, double>& operator=(
            const FFTW_Plan2d<FFTW::fftw_complex, double>&); // teaser
public:
    FFTW_Plan2d<FFTW::fftw_complex, double>(
        int n0,
        int n1,
        FFTW::fftw_complex* in,
        double* out,
        unsigned flags
    )
    :   plan(FFTW::fftw_plan_dft_c2r_2d(n0, n1, in, out, flags))
    {
        if (0 == plan) KJB_THROW(Resource_exhaustion);
    }

    ~FFTW_Plan2d<FFTW::fftw_complex, double>()
    {
        FFTW::fftw_destroy_plan(plan);
    }

    void execute() const
    {
        FFTW::fftw_execute(plan);
    }

    void new_array_exec(FFTW::fftw_complex *in, double* out) const
    {
        FFTW::fftw_execute_dft_c2r(plan, in, out);
    }
};


typedef FFTW_Plan2d<double, FFTW::fftw_complex> FFTW_plan_r2c;

typedef FFTW_Plan2d<FFTW::fftw_complex, double> FFTW_plan_c2r;

#endif






class Fftw_convolution_2d
{
public:
    typedef std::pair< boost::shared_ptr<FFTW_real_vector>,
                       boost::shared_ptr<FFTW_complex_vector> > Work_buffer;

    struct Sizes
    {
        int data_rows, data_cols, mask_rows, mask_cols, pad_rows, pad_cols;
        Sizes(int, int, int, int);

        bool is_matrix_size_same_as_data_size(const Matrix& m) const
        {
            return m.get_num_rows()==data_rows && m.get_num_cols()==data_cols;
        }
        bool is_matrix_size_within_mask_size(const Matrix& m) const
        {
            return m.get_num_rows()<=mask_rows && m.get_num_cols()<=mask_cols;
        }
        int Nreal() const
        {
            return pad_rows * pad_cols;
        }
        int Ncomplex() const
        {
            return pad_rows * (pad_cols/2 + 1);
        }
    };


    Fftw_convolution_2d(
        int data_num_rows,
        int data_num_cols,
        int mask_max_rows,
        int mask_max_cols,
        int fft_alg_type = FFTW_MEASURE
    );

    void set_mask(const Matrix&);

    void set_gaussian_mask(double sigma);

    void convolve(const Matrix&, Matrix&) const;

    void reflect_and_convolve(const Matrix&, Matrix&) const;

    void execute(const Matrix& i, Matrix& o) const { convolve(i, o); }

    // not thread safe
    Work_buffer allocate_work_buffer() const;


/*  ABOVE THE LINE:  NOT THREAD SAFE
 * -----------------------------------------------------------------
 *  BELOW THE LINE:  THREAD SAFE
 */

    void convolve(const Matrix&, Matrix&, Work_buffer) const;

    void reflect_and_convolve(const Matrix&, Matrix&, Work_buffer) const;

    const Sizes& get_sizes() const
    {
        return sizes_;
    }

    bool is_mask_set() const
    {
        return is_mask_set_;
    }

private:

#ifdef KJB_HAVE_FFTW
    Work_buffer alloc_work_buf_impl_() const;

    void reset_data_plan_() const;

    void convolution_by_dft_() const;

    void convolution_by_dft_(Work_buffer) const;
#endif


    // because this field is const, all its members are automatically const.
    const Sizes sizes_;

    const int fft_algorithm_type_;

    bool is_mask_set_;

    Work_buffer mask_;
    mutable Work_buffer data_;


#ifdef KJB_HAVE_FFTW             /* Define plan members in the real thing. */
    boost::scoped_ptr<FFTW_plan_r2c> mask_plan_;
    mutable boost::scoped_ptr<FFTW_plan_r2c> data_plan_;
    mutable boost::scoped_ptr<FFTW_plan_c2r> data_inv_plan_;
#endif

};


inline bool work_buffer_is_unique(const Fftw_convolution_2d::Work_buffer& b)
{
    return b.first.unique() && b.second.unique();
}


namespace debug
{

int test_reflect_into_input_buf(
    const Matrix&,
    Matrix*,
    const Fftw_convolution_2d::Sizes&
);

} // namespace kjb::debug


} // namespace kjb

#endif