SourceXtractorPlusPlus  0.16
Please provide a description of the project.
DFT.h
Go to the documentation of this file.
1 
17 /*
18  * @file SEFramework/Convolution/DFT.h
19  * @date 17/09/18
20  * @author aalvarez
21  */
22 
23 #ifndef _SEFRAMEWORK_CONVOLUTION_DFT_H
24 #define _SEFRAMEWORK_CONVOLUTION_DFT_H
25 
27 #include "SEFramework/FFT/FFT.h"
33 
34 #include <fftw3.h>
35 
36 
37 namespace SourceXtractor {
38 
46 template<typename T = SeFloat, class TPadding = PaddedImage<T, Reflect101Coordinates>>
48 public:
49  typedef T real_t;
50  typedef typename FFT<T>::complex_t complex_t;
51 
58  ConvolutionContext() = default;
59 
60  private:
64 
65  friend class DFTConvolution<T, TPadding>;
66  };
67 
74  : m_kernel{std::move(img)} {
75  }
76 
80  virtual ~DFTConvolution() = default;
81 
87  return m_kernel->getWidth();
88  }
89 
95  return m_kernel->getHeight();
96  }
97 
107  auto context = Euclid::make_unique<ConvolutionContext>();
108 
109  // Dimension of the working padded images
110  context->m_padded_width = model_ptr->getWidth() + m_kernel->getWidth() - 1;
111  context->m_padded_height = model_ptr->getHeight() + m_kernel->getHeight() - 1;
112 
113  // For performance, use a size that is convenient for FFTW
114  context->m_padded_width = fftRoundDimension(context->m_padded_width);
115  context->m_padded_height = fftRoundDimension(context->m_padded_height);
116 
117  // Need to add extra padding for storing the complex part too
118  // (See FFTW documentation)
119  context->m_transform_padding = 2 * (context->m_padded_width / 2 + 1) - context->m_padded_width;
120  int work_area_size = context->m_padded_height * (context->m_padded_width / 2 + 1) * 2;
121 
122  // Pre-allocate buffers for the transformations
123  context->m_kernel_transform.resize(work_area_size);
124  context->m_work_area.resize(work_area_size);
125 
126  // Since we already have the buffers, get the plans too
127  context->m_fwd_plan = FFT<T>::createForwardPlan(context->m_padded_width, context->m_padded_height,
128  context->m_work_area);
129  context->m_inv_plan = FFT<T>::createInversePlan(context->m_padded_width, context->m_padded_height,
130  context->m_work_area);
131 
132  // Transform here the kernel into frequency space
133  padKernel(*context);
134  FFT<T>::executeForward(context->m_fwd_plan, context->m_kernel_transform);
135 
136  return context;
137  }
138 
151  template <typename ...Args>
154  Args... padding_args) const {
155  assert(image_ptr->getWidth() <= context->m_padded_width);
156  assert(image_ptr->getHeight() <= context->m_padded_height);
157 
158  // Padded image
159  auto padded = TPadding::create(image_ptr,
160  context->m_padded_width, context->m_padded_height,
161  std::forward<Args>(padding_args)...);
162 
163  // Create a matrix with the padded image
164  dumpImage(padded, context->m_work_area);
165 
166  // Transform the image
167  FFT<T>::executeForward(context->m_fwd_plan, context->m_work_area);
168 
169  // Multiply the two DFT
170  complex_t* kernel_complex = reinterpret_cast<complex_t*>(context->m_kernel_transform.data());
171  complex_t* img_complex = reinterpret_cast<complex_t*>(context->m_work_area.data());
172  size_t ncomplex = (context->m_padded_width / 2 + 1) * context->m_padded_height;
173  for (size_t i = 0; i < ncomplex; ++i) {
174  const auto& a = img_complex[i];
175  const auto& b = kernel_complex[i];
176  float re = a[0] * b[0] - a[1] * b[1];
177  float im = a[0] * b[1] + a[1] * b[0];
178 
179  img_complex[i][0] = re;
180  img_complex[i][1] = im;
181  }
182 
183  // Inverse DFT
184  FFT<T>::executeInverse(context->m_inv_plan, context->m_work_area);
185 
186  // Copy to the output, removing the pad
187  auto wpad = ::div(context->m_padded_width - image_ptr->getWidth(), 2);
188  auto lpad = wpad.quot;
189  auto rpad = wpad.quot + wpad.rem;
190  auto hpad = ::div(context->m_padded_height - image_ptr->getHeight(), 2);
191  auto tpad = hpad.quot;
192  auto bpad = hpad.quot + hpad.rem;
193  copyFFTWorkAreaToImage(context->m_work_area, *image_ptr, rpad, lpad, tpad, bpad, true);
194  }
195 
208  template <typename ...Args>
209  void convolve(std::shared_ptr<WriteableImage<T>> image_ptr, Args... padding_args) const {
210  auto context = prepare(image_ptr);
211  convolve(image_ptr, context, std::forward(padding_args)...);
212  }
213 
219  return m_kernel;
220  }
221 
222 protected:
223  void padKernel(ConvolutionContext& context) const {
224  auto padded = PaddedImage<T>::create(m_kernel, context.m_padded_width, context.m_padded_height);
225  auto center = PixelCoordinate{context.m_padded_width / 2, context.m_padded_height / 2};
226  if (context.m_padded_width % 2 == 0) center.m_x--;
227  if (context.m_padded_height % 2 == 0) center.m_y--;
228  auto recenter = RecenterImage<T>::create(padded, center);
229 
230  dumpImage(recenter, context.m_kernel_transform);
231  }
232 
233  void dumpImage(const std::shared_ptr<const Image<T>> &img, std::vector<T>& work_area) const {
234  const auto chunk = img->getChunk(0, 0, img->getWidth(), img->getHeight());
235  copyImageToFFTWorkArea(*chunk, work_area);
236  }
237 
238 private:
240 };
241 
242 } // end SourceXtractor
243 
244 #endif // _SEFRAMEWORK_CONVOLUTION_DFT_H
std::size_t getHeight() const
Definition: DFT.h:94
std::size_t getWidth() const
Definition: DFT.h:86
void dumpImage(const std::shared_ptr< const Image< T >> &img, std::vector< T > &work_area) const
Definition: DFT.h:233
std::shared_ptr< const Image< T > > m_kernel
Definition: DFT.h:239
FFT< T >::complex_t complex_t
Definition: DFT.h:50
std::unique_ptr< ConvolutionContext > prepare(const std::shared_ptr< const Image< T >> &model_ptr) const
Definition: DFT.h:106
virtual ~DFTConvolution()=default
std::shared_ptr< const Image< T > > getKernel() const
Definition: DFT.h:218
DFTConvolution(std::shared_ptr< const Image< T >> img)
Definition: DFT.h:73
void convolve(std::shared_ptr< WriteableImage< T >> image_ptr, std::unique_ptr< ConvolutionContext > &context, Args... padding_args) const
Definition: DFT.h:152
void padKernel(ConvolutionContext &context) const
Definition: DFT.h:223
void convolve(std::shared_ptr< WriteableImage< T >> image_ptr, Args... padding_args) const
Definition: DFT.h:209
Interface representing an image.
Definition: Image.h:43
static std::shared_ptr< PaddedImage< T, CoordinateInterpolation > > create(Args &&... args)
Definition: PaddedImage.h:89
static std::shared_ptr< RecenterImage< T > > create(Args &&... args)
Definition: RecenterImage.h:44
T forward(T... args)
T move(T... args)
static void copyFFTWorkAreaToImage(std::vector< T > &buffer, Img< T > &dest, int rpad=0, int lpad=0, int tpad=0, int bpad=0, bool normalize=true)
Definition: FFTHelper.h:87
static void copyImageToFFTWorkArea(Img< T > &origin, std::vector< T > &buffer)
Definition: FFTHelper.h:44
int fftRoundDimension(int size)
Definition: FFT.cpp:49
STL namespace.
std::vector< real_t > m_kernel_transform
Definition: DFT.h:62
static void executeInverse(plan_ptr_t &plan, std::vector< T > &inout)
Definition: FFT.cpp:196
static plan_ptr_t createForwardPlan(int width, int height, std::vector< T > &inout)
Definition: FFT.cpp:111
fftw_traits::complex_t complex_t
Definition: FFT.h:92
static plan_ptr_t createInversePlan(int width, int height, std::vector< T > &inout)
Definition: FFT.cpp:151
static void executeForward(plan_ptr_t &plan, std::vector< T > &inout)
Definition: FFT.cpp:191
A pixel coordinate made of two integers m_x and m_y.