ConvolutionEngine1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_CONVOLUTIONENGINE1D_HH_
#define NPSTAT_CONVOLUTIONENGINE1D_HH_
 
/*!
// \file ConvolutionEngine1D.hh
//
// \brief Fast one-dimensional convolutions via Fourier
//        transforms (FFTW interface)
//
// Author: I. Volobouev
//
// November 2009
*/
 
#include <map>
 
#include "npstat/nm/FourierImage.hh"
#include "npstat/nm/ArrayND.hh"
 
#ifdef SWIG
#include "npstat/wrap/arrayNDToNumpy.hh"
#include "npstat/wrap/numpyArrayUtils.hh"
#endif // SWIG
 
namespace npstat {
    /**
    // Class for performing convolutions by FFT using FFTW library.
    // Internally, transforms are calculated in double precision.
    */
    class ConvolutionEngine1D
    {
    public:
        /**
        // Constructor takes data length and FFTW optimization flag
        // as parameters. Make sure the data length is sufficient
        // to avoid the convolution wrap-around.
        */
        explicit ConvolutionEngine1D(unsigned dataLen,
                                     unsigned optimization = FFTW_ESTIMATE);
        ~ConvolutionEngine1D();
 
        /**
        // Provide the filter for subsequent convolutions.
        //
        // The "shift" argument allows you to rotate the filter data
        // buffer so that the element with index "shift" becomes new
        // element with index 0 (essentially, the buffer is rotated
        // left by "shift" elements).
        */
        template <typename Real>
        void setFilter(const Real* data, unsigned dataLen, unsigned shift=0);
 
        /**
        // Provide the filter for subsequent convolutions using its
        // Fourier space representation
        */
        void setFilterImage(const std::complex<double>* image, unsigned imLen);
 
        /**
        // Convolve provided data with the filter previously established
        // by setFilter.
        //
        // The "shift" argument allows you to rotate the input data
        // buffer so that the element with index "shift" becomes new
        // element with index 0 (essentially, the buffer is rotated
        // left by "shift" elements).
        */
        template <typename Real1, typename Real2>
        void convolveWithFilter(const Real1* in, Real2* out, unsigned dataLen,
                                unsigned shift=0);
 
        /** Deconvolve the input using the previously established filter */
        template <typename Real1, typename Real2>
        void deconvolveFilter(const Real1* in, Real2* out, unsigned dataLen,
                              unsigned shift=0);
 
        /**
        // Provide a filter for subsequent convolutions.
        //
        // The "shift" argument allows you to rotate the filter data
        // buffer so that the element with index "shift" becomes new
        // element with index 0 (essentially, the buffer is rotated
        // left by "shift" elements).
        //
        // "id" is the user id for the filter which can be later reused
        // with the method "convolveWithDeposit".
        */
        template <typename Real>
        void depositFilter(unsigned long id,
                           const Real* data, unsigned dataLen,
                           unsigned shift=0);
 
        /**
        // Provide a filter for subsequent convolutions using its
        // Fourier space representation
        //
        // "id" is the user id for the filter which can be later reused
        // with the method "convolveWithDeposit".
        */
        void depositFilterImage(unsigned long id,
                            const std::complex<double>* image, unsigned imLen);
 
        /**
        // Convolve provided data with one of the filters previously
        // established by the depositFilter method. The "id" argument
        // is used to identify a particular filter.
        //
        // The "shift" argument allows you to rotate the input data
        // buffer so that the element with index "shift" becomes new
        // element with index 0 (essentially, the buffer is rotated
        // left by "shift" elements).
        */
        template <typename Real1, typename Real2>
        void convolveWithDeposit(unsigned long id,
                                 const Real1* in, Real2* out, unsigned dataLen,
                                 unsigned shift=0);
 
        /** Deconvolve the data with a previously deposited filter */
        template <typename Real1, typename Real2>
        void deconvolveDeposit(unsigned long id,
                               const Real1* in, Real2* out, unsigned dataLen,
                               unsigned shift=0);
 
        /**
        // Convolve the filter established by the "setFilter" method
        // with a previously deposited filter
        */
        template <typename Real2>
        void convolveFilterAndDeposit(unsigned long id,
                                      Real2* out, unsigned dataLen);
 
        /**
        // Deconvolve the filter established by the "setFilter" method
        // with a previously deposited filter
        */
        template <typename Real2>
        void deconvolveFilterAndDeposit(unsigned long id,
                                        Real2* out, unsigned dataLen);
 
        /**
        // Discard a filter previously deposited by "depositFilter" method.
        // This method returns "true" if the filter with the given id
        // was indeed discarded and "false" if such a filter was not found.
        */
        bool discardFilter(unsigned long id);
 
        /** Discard all filters deposited by "depositFilter" method */
        void clearFilterBank();
 
        /** Expected data length */
        inline unsigned dataLength() const {return dataLen_;}
 
        /** Expected image length */
        inline unsigned imageLength() const {return dataLen_/2 + 1;}
 
        // Some accessors
        inline bool isFilterSet() const {return validFilter_;}
        inline unsigned filterBankSize() const {return filterMap_.size();}
 
    private:
        typedef std::map<unsigned long,FourierImage*> FilterMap;
 
        ConvolutionEngine1D();
        ConvolutionEngine1D(const ConvolutionEngine1D&);
        ConvolutionEngine1D& operator=(const ConvolutionEngine1D&);
 
        template <typename Real1, typename Real2>
        void convolveWithImage(const fftw_complex* image,
                               const Real1* in, Real2* out, unsigned dataLen,
                               unsigned shift);
 
        template <typename Real1, typename Real2>
        void deconvolveImage(const fftw_complex* image,
                             const Real1* in, Real2* out, unsigned dataLen,
                             unsigned shift);
        fftw_plan pf;
        fftw_plan pb;
 
        double* in;
        fftw_complex* out;
        fftw_complex* filterImage;
 
        unsigned dataLen_;
        bool validFilter_;
 
        FilterMap filterMap_;
 
#ifdef SWIG
    public:
        inline void setFilter2(PyObject* input, unsigned shift=0)
        {
            if (!isNumpyShapeCompatible(input, &dataLen_, 1U)) throw std::invalid_argument(
                "In npstat::ConvolutionEngine1D::setFilter2: incompatible array shape");
            const int typenum = numpyArrayType(input);
            if (typenum == NumpyTypecode<double>::code)
                setFilter((double*)PyArray_DATA((PyArrayObject*)input), dataLen_, shift);
            else if (typenum == NumpyTypecode<float>::code)
                setFilter((float*)PyArray_DATA((PyArrayObject*)input), dataLen_, shift);
            else if (typenum == NumpyTypecode<long long>::code)
                setFilter((long long*)PyArray_DATA((PyArrayObject*)input), dataLen_, shift);
            else if (typenum == NumpyTypecode<long int>::code)
                setFilter((long int*)PyArray_DATA((PyArrayObject*)input), dataLen_, shift);
            else if (typenum == NumpyTypecode<unsigned char>::code)
                setFilter((unsigned char*)PyArray_DATA((PyArrayObject*)input), dataLen_, shift);
            else if (typenum == NumpyTypecode<int>::code)
                setFilter((int*)PyArray_DATA((PyArrayObject*)input), dataLen_, shift);
            else
                throw std::invalid_argument("In npstat::ConvolutionEngine1D::setFilter2: "
                                            "unsupported input array type");
        }
 
        inline void depositFilter2(unsigned long id,
                                   const double* in, unsigned dataLen,
                                   unsigned shift=0)
            {depositFilter(id, in, dataLen, shift);}
 
        inline PyObject* convolveWithFilter2(const double* in, unsigned dataLen,
                                             const unsigned shift = 0,
                                             const bool makeNonNegative = false)
        {
            const int typenum = NumpyTypecode<double>::code;
            npy_intp sh = dataLen;
            PyObject* xarr = PyArray_SimpleNew(1, &sh, typenum);
            if (xarr)
            {
                double* out = (double*)PyArray_DATA((PyArrayObject*)xarr);
                try
                {
                    convolveWithFilter(in, out, dataLen, shift);
                }
                catch (const std::exception& e)
                {
                    Py_DECREF(xarr);
                    throw e;
                }
                if (makeNonNegative)
                    for (unsigned i=0; i<dataLen; ++i)
                        if (out[i] < 0.0)
                            out[i] = 0.0;
            }
            return xarr;
        }
 
        inline PyObject* convolveWithDeposit2(const unsigned long id,
                                              const double* in, unsigned dataLen,
                                              const unsigned shift = 0,
                                              const bool makeNonNegative = false)
        {
            const int typenum = NumpyTypecode<double>::code;
            npy_intp sh = dataLen;
            PyObject* xarr = PyArray_SimpleNew(1, &sh, typenum);
            if (xarr)
            {
                double* out = (double*)PyArray_DATA((PyArrayObject*)xarr);
                try
                {
                    convolveWithDeposit(id, in, out, dataLen, shift);
                }
                catch (const std::exception& e)
                {
                    Py_DECREF(xarr);
                    throw e;
                }
                if (makeNonNegative)
                    for (unsigned i=0; i<dataLen; ++i)
                        if (out[i] < 0.0)
                            out[i] = 0.0;
            }
            return xarr;
        }
 
        inline PyObject* convolveFilterAndDeposit2(const unsigned long id,
                                                   const bool makeNonNegative = false)
        {
            const int typenum = NumpyTypecode<double>::code;
            npy_intp sh = dataLen_;
            PyObject* xarr = PyArray_SimpleNew(1, &sh, typenum);
            if (xarr)
            {
                double* out = (double*)PyArray_DATA((PyArrayObject*)xarr);
                try
                {
                    convolveFilterAndDeposit(id, out, dataLen_);
                }
                catch (const std::exception& e)
                {
                    Py_DECREF(xarr);
                    throw e;
                }
                if (makeNonNegative)
                    for (unsigned i=0; i<dataLen_; ++i)
                        if (out[i] < 0.0)
                            out[i] = 0.0;
            }
            return xarr;
        }
#endif // SWIG
    };
}
 
#include "npstat/nm/ConvolutionEngine1D.icc"
 
#endif // NPSTAT_CONVOLUTIONENGINE1D_HH_