Convolution of two real vectors (in 1D or 2D) using Fast Fourier Transform. More...

Collaboration diagram for Convolve:

Classes
class	Workspace
	Workspace for Fourier convolution. More...

Public Types
typedef std::vector< double >	double1d_t
	definition of 1d vector of double More...

typedef std::vector< double1d_t >	double2d_t
	definition of 2d vector of double More...

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 modes use LINEAR_SAME or CIRCULAR_SAME_SHIFTED for maximum performance More...

Public Member Functions
	Convolve ()

void	fftconvolve (const double1d_t &source, const double1d_t &kernel, double1d_t &result)
	convolution in 1D More...

void	fftconvolve (const double2d_t &source, const double2d_t &kernel, double2d_t &result)
	convolution in 2D More...

void	init (int h_src, int w_src, int h_kernel, int w_kernel)
	prepare arrays for 2D convolution of given vectors More...

void	setMode (EConvolutionMode mode)
	Sets convolution mode. More...

Private Member Functions
void	fftw_circular_convolution (const double2d_t &source, const double2d_t &kernel)
	compute circual convolution of source and kernel using fast Fourier transformation More...

int	find_closest_factor (int n)
	find closest number X>n that can be factorised according to fftw3 favorite factorisation More...

bool	is_optimal (int n)
	if a number can be factorised using only favorite fftw3 factors More...

Private Attributes
std::vector< size_t >	m_implemented_factors

EConvolutionMode	m_mode
	convolution mode More...

Workspace	ws
	input and output data for fftw3 More...

Detailed Description

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

Definition at line 38 of file Convolve.h.

Member Typedef Documentation

◆ double1d_t

typedef std::vector<double> Convolve::double1d_t

definition of 1d vector of double

Definition at line 41 of file Convolve.h.

◆ double2d_t

typedef std::vector<double1d_t> Convolve::double2d_t

definition of 2d vector of double

Definition at line 44 of file Convolve.h.

Member Enumeration Documentation

◆ EConvolutionMode

enum Convolve::EConvolutionMode

convolution modes use LINEAR_SAME or CIRCULAR_SAME_SHIFTED for maximum performance

Enumerator
FFTW_LINEAR_FULL
FFTW_LINEAR_SAME_UNPADDED
FFTW_LINEAR_SAME
FFTW_LINEAR_VALID
FFTW_CIRCULAR_SAME
FFTW_CIRCULAR_SAME_SHIFTED
FFTW_UNDEFINED

Definition at line 50 of file Convolve.h.

                           {
         FFTW_LINEAR_FULL,
         FFTW_LINEAR_SAME_UNPADDED,
         FFTW_LINEAR_SAME,
         FFTW_LINEAR_VALID,
         FFTW_CIRCULAR_SAME,
         FFTW_CIRCULAR_SAME_SHIFTED,
         FFTW_UNDEFINED
     };

Constructor & Destructor Documentation

◆ Convolve()

Convolve::Convolve ( )

Definition at line 20 of file Convolve.cpp.

                    : m_mode(FFTW_UNDEFINED)
 {
     // storing favorite fftw3 prime factors
     const size_t FFTW_FACTORS[] = {13, 11, 7, 5, 3, 2};
     m_implemented_factors.assign(FFTW_FACTORS,
                                  FFTW_FACTORS + sizeof(FFTW_FACTORS) / sizeof(FFTW_FACTORS[0]));
 }

References m_implemented_factors.

Member Function Documentation

◆ fftconvolve() [1/2]

void Convolve::fftconvolve	(	const double1d_t &	source,
		const double1d_t &	kernel,
		double1d_t &	result
	)

convolution in 1D

Definition at line 138 of file Convolve.cpp.

 {
     // we simply create 2d arrays with length of first dimension equal to 1, and call 2d convolution
     double2d_t source2d, kernel2d;
     source2d.push_back(source);
     kernel2d.push_back(kernel);
  
     double2d_t result2d;
     fftconvolve(source2d, kernel2d, result2d);
     if (result2d.size() != 1)
         throw std::runtime_error("Convolve::fftconvolve -> Panic in 1d");
     result = result2d[0];
 }

Referenced by ConvolutionDetectorResolution::apply1dConvolution(), and ConvolutionDetectorResolution::apply2dConvolution().

◆ fftconvolve() [2/2]

void Convolve::fftconvolve	(	const double2d_t &	source,
		const double2d_t &	kernel,
		double2d_t &	result
	)

convolution in 2D

Definition at line 102 of file Convolve.cpp.

 {
     // set default convolution mode, if not defined
     if (m_mode == FFTW_UNDEFINED)
         setMode(FFTW_LINEAR_SAME);
  
     int h_src = (int)source.size();
     int w_src = (int)(source.size() ? source[0].size() : 0);
     int h_kernel = (int)kernel.size();
     int w_kernel = (kernel.size() ? (int)kernel[0].size() : 0);
  
     // initialisation
     init(h_src, w_src, h_kernel, w_kernel);
     // Compute the circular convolution
     fftw_circular_convolution(source, kernel);
  
     // results
     result.clear();
     result.resize(ws.h_dst);
     for (int i = 0; i < ws.h_dst; i++) {
         result[i].resize(ws.w_dst, 0);
         for (int j = 0; j < ws.w_dst; j++) {
             // result[i][j]=ws.dst[i*ws.w_dst+j];
             if (m_mode == FFTW_CIRCULAR_SAME_SHIFTED) {
                 result[i][j] = ws.dst_fft[((i + int(ws.h_kernel / 2.0)) % ws.h_fftw) * ws.w_fftw
                                           + (j + int(ws.w_kernel / 2.0)) % ws.w_fftw];
             } else {
                 result[i][j] = ws.dst_fft[(i + ws.h_offset) * ws.w_fftw + j + ws.w_offset];
             }
         }
     }
 }

References Convolve::Workspace::dst_fft, fftw_circular_convolution(), FFTW_CIRCULAR_SAME_SHIFTED, FFTW_LINEAR_SAME, FFTW_UNDEFINED, Convolve::Workspace::h_dst, Convolve::Workspace::h_fftw, Convolve::Workspace::h_kernel, Convolve::Workspace::h_offset, init(), m_mode, setMode(), Convolve::Workspace::w_dst, Convolve::Workspace::w_fftw, Convolve::Workspace::w_kernel, Convolve::Workspace::w_offset, and ws.

Here is the call graph for this function:

◆ fftw_circular_convolution()

void Convolve::fftw_circular_convolution	(	const double2d_t &	source,
		const double2d_t &	kernel
	)

private

compute circual convolution of source and kernel using fast Fourier transformation

Definition at line 273 of file Convolve.cpp.

 {
     if (ws.h_fftw <= 0 || ws.w_fftw <= 0)
         throw std::runtime_error("Convolve::fftw_convolve() -> Panic! Initialisation is missed.");
  
     double *ptr, *ptr_end, *ptr2;
  
     // Reset the content of ws.in_src
     for (ptr = ws.in_src, ptr_end = ws.in_src + ws.h_fftw * ws.w_fftw; ptr != ptr_end; ++ptr)
         *ptr = 0.0;
     for (ptr = ws.in_kernel, ptr_end = ws.in_kernel + ws.h_fftw * ws.w_fftw; ptr != ptr_end; ++ptr)
         *ptr = 0.0;
  
     // Then we build our periodic signals
     // ptr = src;
     for (int i = 0; i < ws.h_src; ++i)
         for (int j = 0; j < ws.w_src; ++j, ++ptr)
             ws.in_src[(i % ws.h_fftw) * ws.w_fftw + (j % ws.w_fftw)] += src[i][j];
  
     // ptr = kernel;
     for (int i = 0; i < ws.h_kernel; ++i)
         for (int j = 0; j < ws.w_kernel; ++j, ++ptr)
             ws.in_kernel[(i % ws.h_fftw) * ws.w_fftw + (j % ws.w_fftw)] += kernel[i][j];
  
     // And we compute their packed FFT
     fftw_execute(ws.p_forw_src);
     fftw_execute(ws.p_forw_kernel);
  
     // Compute the element-wise product on the packed terms
     // Let's put the element wise products in ws.in_kernel
     double re_s, im_s, re_k, im_k;
     for (ptr = ws.out_src, ptr2 = ws.out_kernel,
         ptr_end = ws.out_src + 2 * ws.h_fftw * (ws.w_fftw / 2 + 1);
          ptr != ptr_end; ++ptr, ++ptr2) {
         re_s = *ptr;
         im_s = *(++ptr);
         re_k = *ptr2;
         im_k = *(++ptr2);
         *(ptr2 - 1) = re_s * re_k - im_s * im_k;
         *ptr2 = re_s * im_k + im_s * re_k;
     }
  
     // Compute the backward FFT
     // Carefull, The backward FFT does not preserve the output
     fftw_execute(ws.p_back);
     // Scale the transform
     for (ptr = ws.dst_fft, ptr_end = ws.dst_fft + ws.w_fftw * ws.h_fftw; ptr != ptr_end; ++ptr)
         *ptr /= double(ws.h_fftw * ws.w_fftw);
 }

References Convolve::Workspace::dst_fft, Convolve::Workspace::h_fftw, Convolve::Workspace::h_kernel, Convolve::Workspace::h_src, Convolve::Workspace::in_kernel, Convolve::Workspace::in_src, Convolve::Workspace::out_kernel, Convolve::Workspace::out_src, Convolve::Workspace::p_back, Convolve::Workspace::p_forw_kernel, Convolve::Workspace::p_forw_src, Convolve::Workspace::w_fftw, Convolve::Workspace::w_kernel, Convolve::Workspace::w_src, and ws.

Referenced by fftconvolve().

◆ find_closest_factor()

int Convolve::find_closest_factor ( int n )

private

find closest number X>n that can be factorised according to fftw3 favorite factorisation

Definition at line 328 of file Convolve.cpp.

 {
     if (is_optimal(n)) {
         return n;
     } else {
         int j = n + 1;
         while (!is_optimal(j))
             ++j;
         return j;
     }
 }

References is_optimal().

Referenced by init().

Here is the call graph for this function:

◆ init()

void Convolve::init	(	int	h_src,
		int	w_src,
		int	h_kernel,
		int	w_kernel
	)

prepare arrays for 2D convolution of given vectors

Definition at line 155 of file Convolve.cpp.

 {
     if (!h_src || !w_src || !h_kernel || !w_kernel) {
         std::ostringstream os;
         os << "Convolve::init() -> Panic! Wrong dimensions " << h_src << " " << w_src << " "
            << h_kernel << " " << w_kernel << std::endl;
         throw std::runtime_error(os.str());
     }
  
     ws.clear();
     ws.h_src = h_src;
     ws.w_src = w_src;
     ws.h_kernel = h_kernel;
     ws.w_kernel = w_kernel;
     switch (m_mode) {
     case FFTW_LINEAR_FULL:
         // Full Linear convolution
         ws.h_fftw = find_closest_factor(h_src + h_kernel - 1);
         ws.w_fftw = find_closest_factor(w_src + w_kernel - 1);
         ws.h_dst = h_src + h_kernel - 1;
         ws.w_dst = w_src + w_kernel - 1;
         // Here we just keep the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
         ws.h_offset = 0;
         ws.w_offset = 0;
         break;
     case FFTW_LINEAR_SAME_UNPADDED:
         // Same Linear convolution
         ws.h_fftw = h_src + int(h_kernel / 2.0);
         ws.w_fftw = w_src + int(w_kernel / 2.0);
         ws.h_dst = h_src;
         ws.w_dst = w_src;
         // Here we just keep the first [h_filt/2:h_filt/2+h_dst-1 ; w_filt/2:w_filt/2+w_dst-1]
         // real part elements of out_src
         ws.h_offset = int(ws.h_kernel / 2.0);
         ws.w_offset = int(ws.w_kernel / 2.0);
         break;
     case FFTW_LINEAR_SAME:
         // Same Linear convolution
         ws.h_fftw = find_closest_factor(h_src + int(h_kernel / 2.0));
         ws.w_fftw = find_closest_factor(w_src + int(w_kernel / 2.0));
         ws.h_dst = h_src;
         ws.w_dst = w_src;
         // Here we just keep the first [h_filt/2:h_filt/2+h_dst-1 ; w_filt/2:w_filt/2+w_dst-1]
         // real part elements of out_src
         ws.h_offset = int(ws.h_kernel / 2.0);
         ws.w_offset = int(ws.w_kernel / 2.0);
         break;
     case FFTW_LINEAR_VALID:
         // Valid Linear convolution
         if (ws.h_kernel > ws.h_src || ws.w_kernel > ws.w_src) {
             ws.h_fftw = 0;
             ws.w_fftw = 0;
             ws.h_dst = 0;
             ws.w_dst = 0;
             std::cout << "The 'valid' convolution results in an empty matrix" << std::endl;
             throw std::runtime_error("The 'valid' convolution results in an empty matrix");
         } else {
             ws.h_fftw = find_closest_factor(h_src);
             ws.w_fftw = find_closest_factor(w_src);
             ws.h_dst = h_src - h_kernel + 1;
             ws.w_dst = w_src - w_kernel + 1;
         }
         // Here we just take [h_dst x w_dst] elements starting at [h_kernel-1;w_kernel-1]
         ws.h_offset = ws.h_kernel - 1;
         ws.w_offset = ws.w_kernel - 1;
         break;
     case FFTW_CIRCULAR_SAME:
     case FFTW_CIRCULAR_SAME_SHIFTED:
         // Circular convolution, modulo N
         // We don't padd with zeros because if we do, we need to padd with at least h_kernel/2;
         // w_kernel/2 elements plus the wrapp around, which in facts leads to too much
         // computations compared to the gain obtained with the optimal size
         ws.h_fftw = h_src;
         ws.w_fftw = w_src;
         ws.h_dst = h_src;
         ws.w_dst = w_src;
         // We copy the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
         ws.h_offset = 0;
         ws.w_offset = 0;
         break;
     default:
         std::cout
             << "Unrecognized convolution mode, possible modes are "
             << "FFTW_LINEAR_FULL, FFTW_LINEAR_SAME, FFTW_LINEAR_SAME_UNPADDED, FFTW_LINEAR_VALID, "
             << "FFTW_CIRCULAR_SAME, FFTW_CIRCULAR_SHIFTED " << std::endl;
         break;
     }
  
     ws.in_src = new double[ws.h_fftw * ws.w_fftw];
     ws.out_src = (double*)fftw_malloc(sizeof(fftw_complex) * ws.h_fftw * (ws.w_fftw / 2 + 1));
     ws.in_kernel = new double[ws.h_fftw * ws.w_fftw];
     ws.out_kernel = (double*)fftw_malloc(sizeof(fftw_complex) * ws.h_fftw * (ws.w_fftw / 2 + 1));
  
     ws.dst_fft = new double[ws.h_fftw * ws.w_fftw];
     // ws.dst = new double[ws.h_dst * ws.w_dst];
  
     // Initialization of the plans
     ws.p_forw_src = fftw_plan_dft_r2c_2d(ws.h_fftw, ws.w_fftw, ws.in_src, (fftw_complex*)ws.out_src,
                                          FFTW_ESTIMATE);
     if (ws.p_forw_src == nullptr)
         throw std::runtime_error("Convolve::init() -> Error! Can't initialise p_forw_src plan.");
  
     ws.p_forw_kernel = fftw_plan_dft_r2c_2d(ws.h_fftw, ws.w_fftw, ws.in_kernel,
                                             (fftw_complex*)ws.out_kernel, FFTW_ESTIMATE);
     if (ws.p_forw_kernel == nullptr)
         throw std::runtime_error("Convolve::init() -> Error! Can't initialise p_forw_kernel plan.");
  
     // The backward FFT takes ws.out_kernel as input
     ws.p_back = fftw_plan_dft_c2r_2d(ws.h_fftw, ws.w_fftw, (fftw_complex*)ws.out_kernel, ws.dst_fft,
                                      FFTW_ESTIMATE);
     if (ws.p_back == nullptr)
         throw std::runtime_error("Convolve::init() -> Error! Can't initialise p_back plan.");
 }

Referenced by fftconvolve().

Here is the call graph for this function:

◆ is_optimal()

bool Convolve::is_optimal ( int n )

private

if a number can be factorised using only favorite fftw3 factors

Definition at line 343 of file Convolve.cpp.

 {
     if (n == 1)
         return false;
     size_t ntest = n;
     for (size_t i = 0; i < m_implemented_factors.size(); i++) {
         while ((ntest % m_implemented_factors[i]) == 0) {
             ntest = ntest / m_implemented_factors[i];
         }
     }
     return ntest == 1;
 }

References m_implemented_factors.

Referenced by find_closest_factor().

◆ setMode()

void Convolve::setMode ( EConvolutionMode mode )

inline

Sets convolution mode.

Definition at line 70 of file Convolve.h.

70 { m_mode = mode; }

References m_mode.

Referenced by fftconvolve().

Member Data Documentation

◆ m_implemented_factors

std::vector<size_t> Convolve::m_implemented_factors

private

Definition at line 122 of file Convolve.h.

Referenced by Convolve(), and is_optimal().

◆ m_mode

EConvolutionMode Convolve::m_mode

private

convolution mode

Definition at line 121 of file Convolve.h.

Referenced by fftconvolve(), init(), and setMode().

◆ ws

Workspace Convolve::ws

private

input and output data for fftw3

Definition at line 119 of file Convolve.h.

Referenced by fftconvolve(), fftw_circular_convolution(), and init().

The documentation for this class was generated from the following files:

/home/www/ba/Device/Resolution/Convolve.h
/home/www/ba/Device/Resolution/Convolve.cpp

Classes

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ double1d_t

◆ double2d_t

Member Enumeration Documentation

◆ EConvolutionMode

Constructor & Destructor Documentation

◆ Convolve()

Member Function Documentation

◆ fftconvolve() [1/2]

◆ fftconvolve() [2/2]

◆ fftw_circular_convolution()

◆ find_closest_factor()

◆ init()

◆ is_optimal()

◆ setMode()

Member Data Documentation

◆ m_implemented_factors

◆ m_mode

◆ ws