BornAgain  1.18.0
Simulate and fit neutron and x-ray scattering at grazing incidence
Convolve.cpp
Go to the documentation of this file.
1 // ************************************************************************** //
2 //
3 // BornAgain: simulate and fit scattering at grazing incidence
4 //
5 //! @file Device/Resolution/Convolve.cpp
6 //! @brief Implements class Convolve.
7 //!
8 //! @homepage http://www.bornagainproject.org
9 //! @license GNU General Public License v3 or higher (see COPYING)
10 //! @copyright Forschungszentrum Jülich GmbH 2018
11 //! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
12 //
13 // ************************************************************************** //
14 
15 #include "Device/Resolution/Convolve.h"
16 #include "Base/Types/Exceptions.h"
17 #include <iostream>
18 #include <sstream>
19 #include <stdexcept> // need overlooked by g++ 5.4
20 
21 Convolve::Convolve() : m_mode(FFTW_UNDEFINED)
22 {
23  // storing favorite fftw3 prime factors
24  const size_t FFTW_FACTORS[] = {13, 11, 7, 5, 3, 2};
25  m_implemented_factors.assign(FFTW_FACTORS,
26  FFTW_FACTORS + sizeof(FFTW_FACTORS) / sizeof(FFTW_FACTORS[0]));
27 }
28 
29 Convolve::Workspace::Workspace()
30  : h_src(0), w_src(0), h_kernel(0), w_kernel(0), w_fftw(0), h_fftw(0), in_src(0), out_src(0),
31  in_kernel(0), out_kernel(0), dst_fft(0), h_dst(0), w_dst(0),
32  // dst(0),
33  h_offset(0), w_offset(0), p_forw_src(nullptr), p_forw_kernel(nullptr), p_back(nullptr)
34 {
35 }
36 
37 Convolve::Workspace::~Workspace()
38 {
39  clear();
40 }
41 
42 void Convolve::Workspace::clear()
43 {
44  h_src = 0;
45  w_src = 0;
46  h_kernel = 0;
47  w_kernel = 0;
48  if (in_src)
49  delete[] in_src;
50  in_src = 0;
51 
52  if (out_src)
53  fftw_free((fftw_complex*)out_src);
54  out_src = 0;
55 
56  if (in_kernel)
57  delete[] in_kernel;
58  in_kernel = 0;
59 
60  if (out_kernel)
61  fftw_free((fftw_complex*)out_kernel);
62  out_kernel = 0;
63 
64  if (dst_fft)
65  delete[] dst_fft;
66  dst_fft = 0;
67 
68  // if(dst) delete[] dst;
69  // dst=0;
70 
71  h_offset = 0;
72  w_offset = 0;
73 
74  if (p_forw_src != nullptr)
75  fftw_destroy_plan(p_forw_src);
76  if (p_forw_kernel != nullptr)
77  fftw_destroy_plan(p_forw_kernel);
78  if (p_back != nullptr)
79  fftw_destroy_plan(p_back);
80 
81  fftw_cleanup();
82 }
83 
84 /* ************************************************************************* */
85 // convolution in 2d
86 /* ************************************************************************* */
87 void Convolve::fftconvolve(const double2d_t& source, const double2d_t& kernel, double2d_t& result)
88 {
89  // set default convolution mode, if not defined
90  if (m_mode == FFTW_UNDEFINED)
91  setMode(FFTW_LINEAR_SAME);
92 
93  int h_src = (int)source.size();
94  int w_src = (int)(source.size() ? source[0].size() : 0);
95  int h_kernel = (int)kernel.size();
96  int w_kernel = (kernel.size() ? (int)kernel[0].size() : 0);
97 
98  // initialisation
99  init(h_src, w_src, h_kernel, w_kernel);
100  // Compute the circular convolution
101  fftw_circular_convolution(source, kernel);
102 
103  // results
104  result.clear();
105  result.resize(ws.h_dst);
106  for (int i = 0; i < ws.h_dst; i++) {
107  result[i].resize(ws.w_dst, 0);
108  for (int j = 0; j < ws.w_dst; j++) {
109  // result[i][j]=ws.dst[i*ws.w_dst+j];
110  if (m_mode == FFTW_CIRCULAR_SAME_SHIFTED) {
111  result[i][j] = ws.dst_fft[((i + int(ws.h_kernel / 2.0)) % ws.h_fftw) * ws.w_fftw
112  + (j + int(ws.w_kernel / 2.0)) % ws.w_fftw];
113  } else {
114  result[i][j] = ws.dst_fft[(i + ws.h_offset) * ws.w_fftw + j + ws.w_offset];
115  }
116  }
117  }
118 }
119 
120 /* ************************************************************************* */
121 // convolution in 1d
122 /* ************************************************************************* */
123 void Convolve::fftconvolve(const double1d_t& source, const double1d_t& kernel, double1d_t& result)
124 {
125  // we simply create 2d arrays with length of first dimension equal to 1, and call 2d convolution
126  double2d_t source2d, kernel2d;
127  source2d.push_back(source);
128  kernel2d.push_back(kernel);
129 
130  double2d_t result2d;
131  fftconvolve(source2d, kernel2d, result2d);
132  if (result2d.size() != 1)
133  throw Exceptions::RuntimeErrorException("Convolve::fftconvolve -> Panic in 1d");
134  result = result2d[0];
135 }
136 
137 /* ************************************************************************* */
138 // initialise input and output arrays for fast Fourier transformation
139 /* ************************************************************************* */
140 void Convolve::init(int h_src, int w_src, int h_kernel, int w_kernel)
141 {
142  if (!h_src || !w_src || !h_kernel || !w_kernel) {
143  std::ostringstream os;
144  os << "Convolve::init() -> Panic! Wrong dimensions " << h_src << " " << w_src << " "
145  << h_kernel << " " << w_kernel << std::endl;
146  throw Exceptions::RuntimeErrorException(os.str());
147  }
148 
149  ws.clear();
150  ws.h_src = h_src;
151  ws.w_src = w_src;
152  ws.h_kernel = h_kernel;
153  ws.w_kernel = w_kernel;
154  switch (m_mode) {
155  case FFTW_LINEAR_FULL:
156  // Full Linear convolution
157  ws.h_fftw = find_closest_factor(h_src + h_kernel - 1);
158  ws.w_fftw = find_closest_factor(w_src + w_kernel - 1);
159  ws.h_dst = h_src + h_kernel - 1;
160  ws.w_dst = w_src + w_kernel - 1;
161  // Here we just keep the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
162  ws.h_offset = 0;
163  ws.w_offset = 0;
164  break;
165  case FFTW_LINEAR_SAME_UNPADDED:
166  // Same Linear convolution
167  ws.h_fftw = h_src + int(h_kernel / 2.0);
168  ws.w_fftw = w_src + int(w_kernel / 2.0);
169  ws.h_dst = h_src;
170  ws.w_dst = w_src;
171  // Here we just keep the first [h_filt/2:h_filt/2+h_dst-1 ; w_filt/2:w_filt/2+w_dst-1]
172  // real part elements of out_src
173  ws.h_offset = int(ws.h_kernel / 2.0);
174  ws.w_offset = int(ws.w_kernel / 2.0);
175  break;
176  case FFTW_LINEAR_SAME:
177  // Same Linear convolution
178  ws.h_fftw = find_closest_factor(h_src + int(h_kernel / 2.0));
179  ws.w_fftw = find_closest_factor(w_src + int(w_kernel / 2.0));
180  ws.h_dst = h_src;
181  ws.w_dst = w_src;
182  // Here we just keep the first [h_filt/2:h_filt/2+h_dst-1 ; w_filt/2:w_filt/2+w_dst-1]
183  // real part elements of out_src
184  ws.h_offset = int(ws.h_kernel / 2.0);
185  ws.w_offset = int(ws.w_kernel / 2.0);
186  break;
187  case FFTW_LINEAR_VALID:
188  // Valid Linear convolution
189  if (ws.h_kernel > ws.h_src || ws.w_kernel > ws.w_src) {
190  ws.h_fftw = 0;
191  ws.w_fftw = 0;
192  ws.h_dst = 0;
193  ws.w_dst = 0;
194  std::cout << "The 'valid' convolution results in an empty matrix" << std::endl;
195  throw std::runtime_error("The 'valid' convolution results in an empty matrix");
196  } else {
197  ws.h_fftw = find_closest_factor(h_src);
198  ws.w_fftw = find_closest_factor(w_src);
199  ws.h_dst = h_src - h_kernel + 1;
200  ws.w_dst = w_src - w_kernel + 1;
201  }
202  // Here we just take [h_dst x w_dst] elements starting at [h_kernel-1;w_kernel-1]
203  ws.h_offset = ws.h_kernel - 1;
204  ws.w_offset = ws.w_kernel - 1;
205  break;
206  case FFTW_CIRCULAR_SAME:
207  case FFTW_CIRCULAR_SAME_SHIFTED:
208  // Circular convolution, modulo N
209  // We don't padd with zeros because if we do, we need to padd with at least h_kernel/2;
210  // w_kernel/2 elements plus the wrapp around, which in facts leads to too much
211  // computations compared to the gain obtained with the optimal size
212  ws.h_fftw = h_src;
213  ws.w_fftw = w_src;
214  ws.h_dst = h_src;
215  ws.w_dst = w_src;
216  // We copy the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
217  ws.h_offset = 0;
218  ws.w_offset = 0;
219  break;
220  default:
221  std::cout
222  << "Unrecognized convolution mode, possible modes are "
223  << "FFTW_LINEAR_FULL, FFTW_LINEAR_SAME, FFTW_LINEAR_SAME_UNPADDED, FFTW_LINEAR_VALID, "
224  << "FFTW_CIRCULAR_SAME, FFTW_CIRCULAR_SHIFTED " << std::endl;
225  break;
226  }
227 
228  ws.in_src = new double[ws.h_fftw * ws.w_fftw];
229  ws.out_src = (double*)fftw_malloc(sizeof(fftw_complex) * ws.h_fftw * (ws.w_fftw / 2 + 1));
230  ws.in_kernel = new double[ws.h_fftw * ws.w_fftw];
231  ws.out_kernel = (double*)fftw_malloc(sizeof(fftw_complex) * ws.h_fftw * (ws.w_fftw / 2 + 1));
232 
233  ws.dst_fft = new double[ws.h_fftw * ws.w_fftw];
234  // ws.dst = new double[ws.h_dst * ws.w_dst];
235 
236  // Initialization of the plans
237  ws.p_forw_src = fftw_plan_dft_r2c_2d(ws.h_fftw, ws.w_fftw, ws.in_src, (fftw_complex*)ws.out_src,
238  FFTW_ESTIMATE);
239  if (ws.p_forw_src == nullptr)
240  throw Exceptions::RuntimeErrorException(
241  "Convolve::init() -> Error! Can't initialise p_forw_src plan.");
242 
243  ws.p_forw_kernel = fftw_plan_dft_r2c_2d(ws.h_fftw, ws.w_fftw, ws.in_kernel,
244  (fftw_complex*)ws.out_kernel, FFTW_ESTIMATE);
245  if (ws.p_forw_kernel == nullptr)
246  throw Exceptions::RuntimeErrorException(
247  "Convolve::init() -> Error! Can't initialise p_forw_kernel plan.");
248 
249  // The backward FFT takes ws.out_kernel as input
250  ws.p_back = fftw_plan_dft_c2r_2d(ws.h_fftw, ws.w_fftw, (fftw_complex*)ws.out_kernel, ws.dst_fft,
251  FFTW_ESTIMATE);
252  if (ws.p_back == nullptr)
253  throw Exceptions::RuntimeErrorException(
254  "Convolve::init() -> Error! Can't initialise p_back plan.");
255 }
256 
257 /* ************************************************************************* */
258 // initialise input and output arrays for fast Fourier transformation
259 /* ************************************************************************* */
260 
261 void Convolve::fftw_circular_convolution(const double2d_t& src, const double2d_t& kernel)
262 {
263  if (ws.h_fftw <= 0 || ws.w_fftw <= 0)
264  throw Exceptions::RuntimeErrorException(
265  "Convolve::fftw_convolve() -> Panic! Initialisation is missed.");
266 
267  double *ptr, *ptr_end, *ptr2;
268 
269  // Reset the content of ws.in_src
270  for (ptr = ws.in_src, ptr_end = ws.in_src + ws.h_fftw * ws.w_fftw; ptr != ptr_end; ++ptr)
271  *ptr = 0.0;
272  for (ptr = ws.in_kernel, ptr_end = ws.in_kernel + ws.h_fftw * ws.w_fftw; ptr != ptr_end; ++ptr)
273  *ptr = 0.0;
274 
275  // Then we build our periodic signals
276  // ptr = src;
277  for (int i = 0; i < ws.h_src; ++i)
278  for (int j = 0; j < ws.w_src; ++j, ++ptr)
279  ws.in_src[(i % ws.h_fftw) * ws.w_fftw + (j % ws.w_fftw)] += src[i][j];
280 
281  // ptr = kernel;
282  for (int i = 0; i < ws.h_kernel; ++i)
283  for (int j = 0; j < ws.w_kernel; ++j, ++ptr)
284  ws.in_kernel[(i % ws.h_fftw) * ws.w_fftw + (j % ws.w_fftw)] += kernel[i][j];
285 
286  // And we compute their packed FFT
287  fftw_execute(ws.p_forw_src);
288  fftw_execute(ws.p_forw_kernel);
289 
290  // Compute the element-wise product on the packed terms
291  // Let's put the element wise products in ws.in_kernel
292  double re_s, im_s, re_k, im_k;
293  for (ptr = ws.out_src, ptr2 = ws.out_kernel,
294  ptr_end = ws.out_src + 2 * ws.h_fftw * (ws.w_fftw / 2 + 1);
295  ptr != ptr_end; ++ptr, ++ptr2) {
296  re_s = *ptr;
297  im_s = *(++ptr);
298  re_k = *ptr2;
299  im_k = *(++ptr2);
300  *(ptr2 - 1) = re_s * re_k - im_s * im_k;
301  *ptr2 = re_s * im_k + im_s * re_k;
302  }
303 
304  // Compute the backward FFT
305  // Carefull, The backward FFT does not preserve the output
306  fftw_execute(ws.p_back);
307  // Scale the transform
308  for (ptr = ws.dst_fft, ptr_end = ws.dst_fft + ws.w_fftw * ws.h_fftw; ptr != ptr_end; ++ptr)
309  *ptr /= double(ws.h_fftw * ws.w_fftw);
310 }
311 
312 /* ************************************************************************* */
313 // find a number closest to the given one, which can be factorised according
314 // to fftw3 favorite factorisation
315 /* ************************************************************************* */
316 
317 int Convolve::find_closest_factor(int n)
318 {
319  if (is_optimal(n)) {
320  return n;
321  } else {
322  int j = n + 1;
323  while (!is_optimal(j))
324  ++j;
325  return j;
326  }
327 }
328 
329 /* ************************************************************************* */
330 // if a number can be factorised using only favorite fftw3 factors
331 /* ************************************************************************* */
332 bool Convolve::is_optimal(int n)
333 {
334  if (n == 1)
335  return false;
336  size_t ntest = n;
337  for (size_t i = 0; i < m_implemented_factors.size(); i++) {
338  while ((ntest % m_implemented_factors[i]) == 0) {
339  ntest = ntest / m_implemented_factors[i];
340  }
341  }
342  return ntest == 1;
343 }
Convolve.h
Defines class MathFunctions::Convolve.
Exceptions.h
Defines many exception classes in namespace Exceptionss.
Convolve::double2d_t
std::vector< double1d_t > double2d_t
definition of 2d vector of double
Definition: Convolve.h:40
Convolve::double1d_t
std::vector< double > double1d_t
definition of 1d vector of double
Definition: Convolve.h:37
Convolve::init
void init(int h_src, int w_src, int h_kernel, int w_kernel)
prepare arrays for 2D convolution of given vectors
Definition: Convolve.cpp:140
Convolve::setMode
void setMode(EConvolutionMode mode)
Sets convolution mode.
Definition: Convolve.h:66
Convolve::fftconvolve
void fftconvolve(const double1d_t &source, const double1d_t &kernel, double1d_t &result)
convolution in 1D
Definition: Convolve.cpp:123