Convolve.h

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Device/Resolution/Convolve.h
//! @brief     Defines class Math::Convolve.
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************
 
#ifdef SWIG
#error no need to expose this header to Swig
#endif
 
#ifndef USER_API
#ifndef BORNAGAIN_DEVICE_RESOLUTION_CONVOLVE_H
#define BORNAGAIN_DEVICE_RESOLUTION_CONVOLVE_H
 
#include <fftw3.h>
#include <vector>
 
//! Convolution of two real vectors (in 1D or 2D) using Fast Fourier Transform.
//!
//! Usage:
//! std::vector<double> signal, kernel, result;
//! Convolve cv;
//! cv.fftconvolve(signal, kernel, result)
//!
//! Given code rely on code from Jeremy Fix page, http://jeremy.fix.free.fr/spip.php?article15,
//! see also "Efficient convolution using the Fast Fourier Transform, Application in C++"
//! by Jeremy Fix, May 30, 2011
class Convolve {
public:
    //! definition of 1d vector of double
    using double1d_t = std::vector<double>;
 
    //! definition of 2d vector of double
    using double2d_t = std::vector<double1d_t>;
 
    Convolve();
 
    //! convolution  modes
    //! use LINEAR_SAME or CIRCULAR_SAME_SHIFTED for maximum performance
    enum EConvolutionMode {
        FFTW_LINEAR_FULL,
        FFTW_LINEAR_SAME_UNPADDED,
        FFTW_LINEAR_SAME,
        FFTW_LINEAR_VALID,
        FFTW_CIRCULAR_SAME,
        FFTW_CIRCULAR_SAME_SHIFTED,
        FFTW_UNDEFINED
    };
 
    //! convolution in 1D
    void fftconvolve(const double1d_t& source, const double1d_t& kernel, double1d_t& result);
 
    //! convolution in 2D
    void fftconvolve(const double2d_t& source, const double2d_t& kernel, double2d_t& result);
 
    //! prepare arrays for 2D convolution of given vectors
    void init(int h_src, int w_src, int h_kernel, int w_kernel);
 
    //! Sets convolution mode
    void setMode(EConvolutionMode mode) { m_mode = mode; }
 
private:
    //! compute circual convolution of source and kernel using fast Fourier transformation
    void fftw_circular_convolution(const double2d_t& src, const double2d_t& kernel);
 
    //! find closest number X>n that can be factorized according to fftw3 favorite factorization
    int find_closest_factor(int n);
 
    //! if a number can be factorized using only favorite fftw3 factors
    bool is_optimal(int n);
 
    //! Workspace for Fourier convolution.
 
    //! Workspace contains input (source and kernel), intermediate and output
    //! arrays to run convolution via fft; 'source' it is our signal, 'kernel'
    //! it is our resolution function.
    //! Sizes of input arrays are adjusted; output arrays are alocated via
    //! fftw3 allocation for maximum performance.
    class Workspace {
    public:
        Workspace() = default;
        ~Workspace();
        void clear();
        friend class Convolve;
 
    private:
        int h_src{0}, w_src{0};       // size of original 'source' array in 2 dimensions
        int h_kernel{0}, w_kernel{0}; // size of original 'kernel' array in 2 dimensions
        // size of adjusted source and kernel arrays (in_src, out_src, in_kernel, out_kernel)
        int w_fftw{0}, h_fftw{0};
        //! adjusted input 'source' array
        double* in_src{nullptr};
        //! result of Fourier transformation of source
        double* out_src{nullptr};
        //! adjusted input 'kernel' array
        double* in_kernel{nullptr};
        //! result of Fourier transformation of kernel
        double* out_kernel{nullptr};
        //! result of product of FFT(source) and FFT(kernel)
        double* dst_fft{nullptr};
        int h_dst{0}, w_dst{0};       // size of resulting array
        int h_offset{0}, w_offset{0}; // offsets to copy result into output arrays
        fftw_plan p_forw_src{nullptr};
        fftw_plan p_forw_kernel{nullptr};
        fftw_plan p_back{nullptr};
    };
 
    //! input and output data for fftw3
    Workspace ws;
    //! convolution mode
    EConvolutionMode m_mode;
    std::vector<size_t> m_implemented_factors; // favorite factorization terms of fftw3
};
 
#endif // BORNAGAIN_DEVICE_RESOLUTION_CONVOLVE_H
#endif // USER_API