1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
|
// ************************************************************************************************
//
// 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 // SWIG
#ifndef BORNAGAIN_DEVICE_RESOLUTION_CONVOLVE_H
#define BORNAGAIN_DEVICE_RESOLUTION_CONVOLVE_H
#include "Base/Type/Field2D.h"
#include <fftw3.h>
//! 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)
//!
//! Inspired by code from Jeremy Fix.
//! Original provenance http://jeremy.fix.free.fr/spip.php?article15 no longer online.
//! Code repository now at https://github.com/jeremyfix/FFTConvolution.
//! Paper "Efficient convolution using the Fast Fourier Transform, Application in C++"
//! by Jeremy Fix, May 30, 2011, preserved in ba-inter/lit/JFix_FFTConvolution.pdf
class Convolve {
public:
Convolve();
//! convolution modes
//! use LINEAR_SAME or CIRCULAR_SAME_SHIFTED for maximum performance
enum ConvolutionMode {
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 fftconvolve1D(const std::vector<double>& source, const std::vector<double>& kernel,
std::vector<double>& result);
//! convolution in 2D
void fftconvolve2D(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);
private:
//! compute circual convolution of source and kernel using fast Fourier transformation
void fftw_circular_convolution(const double2d_t& src, const double2d_t& kernel);
//! 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
ConvolutionMode m_mode;
};
#endif // BORNAGAIN_DEVICE_RESOLUTION_CONVOLVE_H
|
