From feb171c696fde6174d3627d669bfcfed5bd13f94 Mon Sep 17 00:00:00 2001
From: pospelov <pospelov@fz-juelich.de>
Date: Mon, 11 Jun 2012 13:21:17 +0200
Subject: [PATCH] Functional test for Convolve class, code cleanup.
---
App/App.pro | 4 +-
.../{TestInstrument.h => TestConvolution.h} | 19 +-
...TestInstrument.cpp => TestConvolution.cpp} | 276 ++++-----
App/src/TestFactory.cpp | 14 +-
Core/inc/Convolve.h | 52 +-
Core/src/Convolve.cpp | 556 ++++--------------
6 files changed, 262 insertions(+), 659 deletions(-)
rename App/inc/{TestInstrument.h => TestConvolution.h} (75%)
rename App/src/{TestInstrument.cpp => TestConvolution.cpp} (57%)
diff --git a/App/App.pro b/App/App.pro
index 5e1188a8c84..b3af846a5cd 100644
--- a/App/App.pro
+++ b/App/App.pro
@@ -17,7 +17,7 @@ SOURCES += \
src/CommandLine.cpp \
src/TestFactory.cpp \
src/TestDiffuseReflection.cpp \
- src/TestInstrument.cpp
+ src/TestConvolution.cpp
HEADERS += \
inc/DrawHelper.h \
@@ -32,7 +32,7 @@ HEADERS += \
inc/CommandLine.h \
inc/TestFactory.h \
inc/TestDiffuseReflection.h \
- inc/TestInstrument.h
+ inc/TestConvolution.h
INCLUDEPATH += ./inc
DEPENDPATH += ./inc
diff --git a/App/inc/TestInstrument.h b/App/inc/TestConvolution.h
similarity index 75%
rename from App/inc/TestInstrument.h
rename to App/inc/TestConvolution.h
index 79073497c70..b3b375894c3 100644
--- a/App/inc/TestInstrument.h
+++ b/App/inc/TestConvolution.h
@@ -1,5 +1,5 @@
-#ifndef TESTINSTRUMENT_H
-#define TESTINSTRUMENT_H
+#ifndef TESTCONVOLUTION_H
+#define TESTCONVOLUTION_H
// ********************************************************************
// * The BornAgain project *
// * Simulation of neutron and x-ray scattering at grazing incidence *
@@ -24,10 +24,10 @@
//! @class TestConvolution
//! @brief Testing Convolve class for instrumental effects studies
//- -------------------------------------------------------------------
-class TestInstrument : public IFunctionalTest
+class TestConvolution : public IFunctionalTest
{
public:
- TestInstrument();
+ TestConvolution();
void execute();
@@ -38,8 +38,17 @@ public:
void test_convolve2d();
private:
+ //! test function with many gaus'es on top of flat background for convolution studies
+ double fpeaks(double *x, double *par);
+
+ //! typedef of pair for the description of convolution mode
typedef std::pair<std::string, MathFunctions::Convolve::Mode> mode_pair_t;
+
+ //! vector of pairs defined above
std::vector<mode_pair_t> m_modes;
+
+ //! number of peaks for convolution test function
+ int m_npeaks;
};
-#endif // TESTINSTRUMENT_H
+#endif // TESTCONVOLUTION_H
diff --git a/App/src/TestInstrument.cpp b/App/src/TestConvolution.cpp
similarity index 57%
rename from App/src/TestInstrument.cpp
rename to App/src/TestConvolution.cpp
index 98923c866ad..c48eebadf31 100644
--- a/App/src/TestInstrument.cpp
+++ b/App/src/TestConvolution.cpp
@@ -1,4 +1,4 @@
-#include "TestInstrument.h"
+#include "TestConvolution.h"
#include "MathFunctions.h"
#include "Convolve.h"
#include "DrawHelper.h"
@@ -16,29 +16,14 @@
#include "TGraph.h"
#include "TBenchmark.h"
#include "TStyle.h"
+#include "TText.h"
-std::vector<double > Convolve_MyFFT1D_B(int mode, const std::vector<double> &data1, const std::vector<double> &data2);
-int npeaks(30);
-Double_t fpeaks(Double_t *x, Double_t *par) {
- Double_t result = par[0] + par[1]*x[0];
- for (Int_t p=0;p<npeaks;p++) {
- Double_t norm = par[3*p+2];
- Double_t mean = par[3*p+3];
- Double_t sigma = par[3*p+4];
- result += norm*TMath::Gaus(x[0],mean,sigma);
- }
- return result;
-}
-
-
-
-
-TestInstrument::TestInstrument()
+TestConvolution::TestConvolution()
{
- std::cout << "TestInstrument::TestInstrument() -> Info." << std::endl;
+ std::cout << "TestConvolution::TestConvolution() -> Info." << std::endl;
m_modes.push_back( mode_pair_t("LINEAR_FULL", MathFunctions::Convolve::FFTW_LINEAR_FULL));
m_modes.push_back( mode_pair_t("LINEAR_SAME_UNPADDED", MathFunctions::Convolve::FFTW_LINEAR_SAME_UNPADDED));
@@ -48,107 +33,120 @@ TestInstrument::TestInstrument()
}
-void TestInstrument::execute()
+void TestConvolution::execute()
{
test_convolve1d();
- // results are about
- // LINEAR_FULL: Real Time = 0.37 seconds Cpu Time = 0.36 seconds
- // LINEAR_SAME_UNPADDED: Real Time = 1.03 seconds Cpu Time = 1.03 seconds
- // LINEAR_SAME: Real Time = 0.23 seconds Cpu Time = 0.23 seconds
- // CIRCULAR_SAME: Real Time = 0.50 seconds Cpu Time = 0.50 seconds
- // CIRCULAR_SAME_SHIFTED: Real Time = 0.51 seconds Cpu Time = 0.51 seconds
-
- //test_convolve2d();
- // LINEAR_FULL: Real Time = 1.19 seconds Cpu Time = 1.18 seconds
- // LINEAR_SAME_UNPADDED: Real Time = 0.72 seconds Cpu Time = 0.72 seconds
- // LINEAR_SAME: Real Time = 0.72 seconds Cpu Time = 0.73 seconds
- // CIRCULAR_SAME: Real Time = 0.53 seconds Cpu Time = 0.54 seconds
- // CIRCULAR_SAME_SHIFTED: Real Time = 0.61 seconds Cpu Time = 0.61 seconds
-}
+ test_convolve2d();
+
+ /*
+ --- TestConvolution::test_convolve_1d() -> Peformance.
+ LINEAR_FULL: Real Time = 0.37 seconds Cpu Time = 0.37 seconds
+ LINEAR_SAME_UNPADDED: Real Time = 0.98 seconds Cpu Time = 0.98 seconds
+ LINEAR_SAME: Real Time = 0.22 seconds Cpu Time = 0.22 seconds
+ CIRCULAR_SAME: Real Time = 0.51 seconds Cpu Time = 0.50 seconds
+ CIRCULAR_SAME_SHIFTED: Real Time = 0.51 seconds Cpu Time = 0.51 seconds
+ --- TestConvolution::test_convolve_2d() -> Peformance.
+ LINEAR_FULL: Real Time = 1.26 seconds Cpu Time = 1.26 seconds
+ LINEAR_SAME_UNPADDED: Real Time = 0.76 seconds Cpu Time = 0.76 seconds
+ LINEAR_SAME: Real Time = 0.73 seconds Cpu Time = 0.73 seconds
+ CIRCULAR_SAME: Real Time = 0.57 seconds Cpu Time = 0.57 seconds
+ CIRCULAR_SAME_SHIFTED: Real Time = 0.62 seconds Cpu Time = 0.63 seconds
+ */
+
+}
-void TestInstrument::test_convolve1d()
+/* ************************************************************************* */
+// test of convolution of two arrays in 1d
+/* ************************************************************************* */
+void TestConvolution::test_convolve1d()
{
- npeaks = TMath::Abs(10);
double xmin(0.0);
double xmax(1000.);
size_t npoints = 501;
- double dx = (xmax-xmin)/(npoints-1);
-
- double sigma0=(xmax-xmin)*0.001;
+ double dx = (xmax-xmin)/double(npoints-1);
+ double sigmax=(xmax-xmin)*0.01; // precision of the measurement of x-value
- TH1F *href = new TH1F("h","test", npoints-1, xmin, xmax);
- href->SetMinimum(0.0);
- href->SetMaximum(2.0);
-
- //generate n peaks at random
- Double_t par[3000];
+ // generate signal (n gaus peaks at random on top of flat background)
+ std::vector<double > signal;
+ m_npeaks = 10;
+ double *par = new double[m_npeaks*3+2];
par[0] = 1.0;
par[1] = -0.6/(xmax-xmin);
- Int_t p;
- for (p=0;p<npeaks;p++) {
- par[3*p+2] = 1;
- par[3*p+3] = (xmax-xmin)*0.01+gRandom->Rndm()*(xmax-xmin);
- par[3*p+4] = sigma0+sigma0*5*gRandom->Rndm();
+ for (int i_par=0; i_par<m_npeaks;i_par++) {
+ par[3*i_par+2] = 1;
+ par[3*i_par+3] = (xmax-xmin)*0.01+gRandom->Rndm()*(xmax-xmin);
+ par[3*i_par+4] = sigmax/5.+sigmax*gRandom->Rndm();
}
- TF1 *fspect = new TF1("f",fpeaks,0,1000,2+3*npeaks);
- fspect->SetNpx(10000);
- fspect->SetParameters(par);
- TCanvas *c1 = new TCanvas("c1_test_convolve1d","c1_test_convolve1d",1024, 768);
- c1->Divide(3,3);
- c1->cd(1);
- //href->FillRandom("f",20000);
- href->DrawCopy();
- fspect->Draw("pl same");
-
- c1->cd(2);
- double sigma_measure = sigma0*5;
- TH1F *hmeasured = new TH1F("hp","hp",500, xmin, xmax);
-// double sigma_measure=sigma0*3;
- for(int i=0; i<1000000; i++) {
- hmeasured->Fill(fspect->GetRandom(xmin, xmax) + gRandom->Gaus(0.0,sigma_measure));
- }
- hmeasured->Scale(1./hmeasured->GetEntries());
- hmeasured->SetLineColor(kRed);
- hmeasured->Draw();
-
- std::vector<double> xspect;
- xspect.resize(npoints);
+ TF1 *fsignal = new TF1("f",this, &TestConvolution::fpeaks, xmin, xmax,2+3*m_npeaks, "TestConvolution","fpeaks");
+ fsignal->SetNpx(10000);
+ fsignal->SetParameters(par);
+ signal.resize(npoints);
for(size_t i=0; i<npoints; i++) {
double x = xmin + i*dx;
- xspect[i] = fspect->Eval(x);
- std::cout << i << " " << x << " " << xspect[i] << std::endl;
+ signal[i] = fsignal->Eval(x);
}
- // building responce function
- std::vector<double> xresp;
- xresp.resize(npoints);
+ // building kernel
+ std::vector<double> kernel;
+ kernel.resize(npoints);
+ TGraph *gr_kernel = new TGraph(npoints);
// for(size_t i=0; i<npoints/2; i++) {
// double x = i*dx;
// double y = std::exp(-x*x/sigma_measure/sigma_measure);
-// xresp[i] = y;
-// xresp[npoints-1-i] = xresp[i];
+// kernel[i] = y;
+// kernel[npoints-1-i] = kernel[i];
// }
-
for(size_t i=0; i<npoints; i++) {
double x = -(xmax-xmin)/2.+i*dx;
- double y = std::exp(-x*x/sigma_measure/sigma_measure);
- xresp[i] = y;
+ double y = std::exp(-x*x/sigmax/sigmax);
+ kernel[i] = y;
+ gr_kernel->SetPoint(i,x,kernel[i]);
}
- TGraph *gr_resp = new TGraph(npoints);
- for(size_t i=0; i<xresp.size(); i++){
- double x = xmin + i*dx;
- gr_resp->SetPoint(i,x,xresp[i]);
- }
- c1->cd(3);
+
+ // --------------
+ // drawing
+ // --------------
+ TCanvas *c1 = new TCanvas("c1_test_convolve1d","c1_test_convolve1d",1024, 768);
+ DrawHelper::instance().SetMagnifier(c1);
+ c1->Divide(3,3);
+
+ // drawing signal
+ c1->cd(1);
+ gPad->SetTopMargin(0.1);
+ TH1F *href = new TH1F("h","test", npoints-1, xmin, xmax);
+ href->SetMinimum(0.0);
+ href->SetMaximum(2.0);
+ href->SetStats(0);
href->DrawCopy();
- gr_resp->Draw("pl same");
+ href->SetBit(TH1::kNoTitle, false);
+ fsignal->SetTitle("signal");
+ fsignal->Draw("pl same");
+
+ // drawing measured signal
+ c1->cd(2);
+ gPad->SetTopMargin(0.1);
+ TH1F *hmeasured = new TH1F("hmeasured","hmeasured",500, xmin, xmax);
+ for(int i=0; i<1000000; i++) {
+ hmeasured->Fill(fsignal->GetRandom(xmin, xmax) + gRandom->Gaus(0.0,sigmax));
+ }
+ hmeasured->Scale(1./hmeasured->GetEntries());
+ hmeasured->SetLineColor(kRed);
+ hmeasured->SetStats(0);
+ hmeasured->SetBit(TH1::kNoTitle, false);
+ hmeasured->Draw();
+ // drawing kernel
+ c1->cd(3);
+ gPad->SetTopMargin(0.1);
+ gr_kernel->SetTitle("kernel");
+ gr_kernel->Draw("apl");
+ // calculation of convolution
TBenchmark benchmark;
for(size_t i_mode=0; i_mode<m_modes.size(); i_mode++) {
std::string sname = m_modes[i_mode].first;
@@ -160,34 +158,34 @@ void TestInstrument::test_convolve1d()
MathFunctions::Convolve cv;
cv.setMode(mode);
- std::vector<double> xresult;
+ std::vector<double> result;
for(int i=0; i<1000; i++) {
- //int mm = (int)modes[i_mode];
- // xresult = Convolve_MyFFT1D_B(mm, xspect, xresp);
- cv.fftconvolve(xspect, xresp, xresult);
+ cv.fftconvolve(signal, kernel, result);
}
benchmark.Stop(sname.c_str());
- benchmark.Show(sname.c_str());
- for(size_t i=0; i<xresult.size(); i++) {
+ //benchmark.Show(sname.c_str());
+ for(size_t i=0; i<result.size(); i++) {
double x = xmin + i*dx;
- gr_result->SetPoint(i,x,xresult[i]);
+ gr_result->SetPoint(i,x,result[i]);
}
c1->cd(4+i_mode);
+ gPad->SetTopMargin(0.1);
href->SetMaximum(20.0);
+ href->SetTitle(sname.c_str());
href->DrawCopy();
gr_result->Draw("pl same");
}
float_t rp, cp;
- std::cout << "--------------" << std::endl;
+ std::cout << "--- TestConvolution::test_convolve_1d() -> Peformance." << std::endl;
benchmark.Summary(rp, cp);
}
/* ************************************************************************* */
-//
+// test of convolution of two arrays in 2d
/* ************************************************************************* */
-void TestInstrument::test_convolve2d()
+void TestConvolution::test_convolve2d()
{
double xmin(0.0), xmax(100.);
size_t npx = 100;
@@ -211,13 +209,12 @@ void TestInstrument::test_convolve2d()
// spikes
if(ix>1 && iy!=0) {
if( (ix%20==0 || (ix-1)%20==0) && iy%8==0) {
- std::cout << ix << " " << iy << std::endl;
source[ix][iy]+=5;
}
}
} // iy
} // ix
- // signal histogram
+ // filling signal histogram with values from signal array
TH2F *h2_signal = new TH2F("h2_signal","h2_signal", npx, xmin, xmax, npy, ymin, ymax);
h2_signal->SetStats(0);
h2_signal->SetMinimum(0);
@@ -236,6 +233,7 @@ void TestInstrument::test_convolve2d()
// preparing resolution kernel and histogram
std::vector<std::vector<double> > kernel;
TH2F *h2_kernel = new TH2F("h2_kernel","h2_kernel", npx, xmin, xmax, npy, ymin, ymax);
+ h2_kernel->SetStats(0);
kernel.resize(npx);
for(size_t ix=0; ix<npx; ix++) {
kernel[ix].resize(npy);
@@ -249,6 +247,7 @@ void TestInstrument::test_convolve2d()
// simulation of measurements with smearing measurement with resolution function
TH2F *h2_measured = new TH2F("h2_measured","h2_measured", npx, xmin, xmax, npy, ymin, ymax);
+ h2_measured->SetStats(0);
TRandom mr;
for(int i_meas=0; i_meas<1000000; i_meas++) {
// generating distribution according to our signal
@@ -264,16 +263,24 @@ void TestInstrument::test_convolve2d()
gStyle->SetPalette(1);
c1->Divide(3,3);
- c1->cd(1); gPad->SetRightMargin(0.15);
+ c1->cd(1);
+ gPad->SetRightMargin(0.12);
+ gPad->SetTopMargin(0.1);
h2_signal->Draw("colz");
- c1->cd(2); gPad->SetRightMargin(0.15);
+ c1->cd(2);
+ gPad->SetRightMargin(0.12);
+ gPad->SetTopMargin(0.1);
h2_signal->Draw("lego");
- c1->cd(3); gPad->SetRightMargin(0.15);
+ c1->cd(3);
+ gPad->SetRightMargin(0.12);
+ gPad->SetTopMargin(0.1);
h2_kernel->Draw("lego");
- c1->cd(4); gPad->SetRightMargin(0.15);
+ c1->cd(4);
+ gPad->SetRightMargin(0.12);
+ gPad->SetTopMargin(0.1);
h2_measured->Draw("colz");
// calculation of convolution
@@ -294,6 +301,8 @@ void TestInstrument::test_convolve2d()
// drawing results of last convolution
c1->cd(5+i_mode);
+ gPad->SetRightMargin(0.12);
+ gPad->SetTopMargin(0.1);
TH2F h2_result("h2_result","h2_result", npx, xmin, xmax, npy, ymin, ymax);
h2_result.SetStats(0);
for(size_t ix=0; ix<npx; ix++) {
@@ -304,61 +313,32 @@ void TestInstrument::test_convolve2d()
}
}
h2_result.SetTitle(sname.c_str());
+ h2_result.SetBit(TH1::kNoTitle, false);
h2_result.DrawCopy("colz");
}
// benchmark summary
float_t rp, cp;
- std::cout << "--------------" << std::endl;
+ std::cout << "--- TestConvolution::test_convolve_2d() -> Peformance." << std::endl;
benchmark.Summary(rp, cp);
}
-
-
-std::vector<double> Convolve_MyFFT1D_B(int mode, const std::vector<double> &data1, const std::vector<double> &data2)
+/* ************************************************************************* */
+// test function for convolution
+/* ************************************************************************* */
+double TestConvolution::fpeaks(double *x, double *par)
{
- size_t Ns = data1.size();
- double * src = new double[Ns];
- for(size_t i=0; i<Ns; i++) {
- src[i] = data1[i];
+ Double_t result = par[0] + par[1]*x[0];
+ for (Int_t p=0;p<m_npeaks;p++) {
+ Double_t norm = par[3*p+2];
+ Double_t mean = par[3*p+3];
+ Double_t sigma = par[3*p+4];
+ result += norm*TMath::Gaus(x[0],mean,sigma);
}
-
- size_t Nk = data2.size();
- double * kernel = new double[Nk];
- for(size_t i=0; i<Nk; i++) {
- kernel[i] = data2[i];
- }
-
- FFTW_Workspace ws;
-// init_workspace_fftw(ws, FFTW_LINEAR_FULL, 1,Ns,1,Nk);
-// init_workspace_fftw(ws, FFTW_LINEAR_SAME_UNPADDED, 1,Ns,1,Nk);
-// init_workspace_fftw(ws, FFTW_LINEAR_SAME, 1,Ns,1,Nk);
-// init_workspace_fftw(ws, FFTW_LINEAR_VALID, 1,Ns,1,Nk);
-// init_workspace_fftw(ws, FFTW_CIRCULAR_SAME, 1,Ns,1,Nk);
-// init_workspace_fftw(ws, FFTW_CIRCULAR_SAME_SHIFTED, 1,Ns,1,Nk);
-
- FFTW_Convolution_Mode xm = (FFTW_Convolution_Mode)mode;
- init_workspace_fftw(ws, xm, 1,Ns,1,Nk);
-
-
- fftw_convolve(ws, src, kernel);
-
-// std::cout << "!!! " << ws.w_dst << std::endl;
-// for(int j = 0 ; j < ws.w_dst ; ++j) std::cout << ws.dst[j] << " ";
-
- std::vector<double > fft_result;
- fft_result.resize(ws.w_dst);
- for(size_t i=0; i<data1.size(); i++){
- fft_result[i] = ws.dst[i];
- }
-
- clear_workspace_fftw(ws);
- return fft_result;
+ return result;
}
-
-
diff --git a/App/src/TestFactory.cpp b/App/src/TestFactory.cpp
index 97b6cf7705f..c60096d74f6 100644
--- a/App/src/TestFactory.cpp
+++ b/App/src/TestFactory.cpp
@@ -4,7 +4,7 @@
#include "TestFormFactor.h"
#include "TestDWBAFormFactor.h"
#include "TestDiffuseReflection.h"
-#include "TestInstrument.h"
+#include "TestConvolution.h"
#include "TBenchmark.h"
@@ -15,12 +15,12 @@ TestFactory::TestFactory() : m_benchmark(0)
setStoreObjects(true);
setDeleteObjects(true);
- registerItem("roughness", IFactoryCreateFunction<TestRoughness, IFunctionalTest> );
- registerItem("fresnel", IFactoryCreateFunction<TestFresnelCoeff, IFunctionalTest> );
- registerItem("formfactor", IFactoryCreateFunction<TestFormFactor, IFunctionalTest> );
- registerItem("dwba", IFactoryCreateFunction<TestDWBAFormFactor, IFunctionalTest> );
- registerItem("diffuse", IFactoryCreateFunction<TestDiffuseReflection, IFunctionalTest> );
- registerItem("instrument", IFactoryCreateFunction<TestInstrument, IFunctionalTest> );
+ registerItem("roughness", IFactoryCreateFunction<TestRoughness, IFunctionalTest> );
+ registerItem("fresnel", IFactoryCreateFunction<TestFresnelCoeff, IFunctionalTest> );
+ registerItem("formfactor", IFactoryCreateFunction<TestFormFactor, IFunctionalTest> );
+ registerItem("dwba", IFactoryCreateFunction<TestDWBAFormFactor, IFunctionalTest> );
+ registerItem("diffuse", IFactoryCreateFunction<TestDiffuseReflection, IFunctionalTest> );
+ registerItem("convolution", IFactoryCreateFunction<TestConvolution, IFunctionalTest> );
m_benchmark = new TBenchmark();
}
diff --git a/Core/inc/Convolve.h b/Core/inc/Convolve.h
index 2ba520f64b6..7b0671a8cf8 100644
--- a/Core/inc/Convolve.h
+++ b/Core/inc/Convolve.h
@@ -52,6 +52,7 @@ public:
typedef std::vector<double1d_t> double2d_t;
//! convolution modes
+ //! use LINEAR_SAME or CIRCULAR_SAME_SHIFTED for maximum performance
enum Mode { 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
@@ -67,12 +68,7 @@ public:
void setMode(Mode mode) { m_mode = mode; }
private:
- //! compute convolution of source and kernel using fast fourier transformation
- void fftw_convolve(const double2d_t &source, const double2d_t &kernel);
- //void fftw_convolve(double * src,double * kernel);
-
//! compute circual convolution of source and kernel using fast fourier transformation
- //void fftw_circular_convolution(double * src, double * kernel);
void fftw_circular_convolution(const double2d_t &source, const double2d_t &kernel);
//! find closest number X>n, which can be factorised according to fftw3 favorite factorisation
@@ -82,7 +78,7 @@ private:
bool is_optimal(int n);
//! 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 (or delta-responce) function
+ //! 'source' it is our signal, 'kernel' it is our resolution (also known as delta-responce) function
//! Sizes of input arrays are adjusted; output arrays are alocated via fftw3 allocation for maximum performance
class Workspace
{
@@ -102,6 +98,7 @@ private:
double *dst_fft; // result of production of FFT(source) and FFT(kernel)
int h_dst, w_dst; // size of resulting array
double *dst; // The array containing the result
+ int h_offset, w_offset; // offsets to copy result into output arrays
fftw_plan p_forw_src;
fftw_plan p_forw_kernel;
fftw_plan p_back;
@@ -116,47 +113,4 @@ private:
} // namespace MathFunctions
-typedef enum
-{
- FFTW_LINEAR_FULL,
- FFTW_LINEAR_SAME_UNPADDED,
- FFTW_LINEAR_SAME,
- FFTW_LINEAR_VALID,
- FFTW_CIRCULAR_SAME,
- FFTW_CIRCULAR_SAME_SHIFTED,
- FFTW_CIRCULAR_CUSTOM
-} FFTW_Convolution_Mode;
-
-typedef struct FFTW_Workspace
-{
- //fftw_complex * in_src, *out_src;
- double * in_src, *out_src, *in_kernel, *out_kernel;
- int h_src, w_src, h_kernel, w_kernel;
- int w_fftw, h_fftw;
- FFTW_Convolution_Mode mode;
- double * dst_fft;
- double * dst; // The array containing the result
- int h_dst, w_dst; // its size ; This is automatically set by init_workspace, despite you use FFTW_CIRCULAR_CUSTOM
- fftw_plan p_forw_src;
- fftw_plan p_forw_kernel;
- fftw_plan p_back;
-
-} FFTW_Workspace;
-
-extern void factorize (const int n,
- int *n_factors,
- int factors[],
- int * implemented_factors);
-
-extern bool is_optimal(int n, int * implemented_factors);
-extern int find_closest_factor(int n, int * implemented_factor);
-extern void init_workspace_fftw(FFTW_Workspace & ws, FFTW_Convolution_Mode mode, int h_src, int w_src, int h_kernel, int w_kernel);
-extern void clear_workspace_fftw(FFTW_Workspace & ws);
-extern void fftw_convolve(FFTW_Workspace &ws, double * src,double * kernel);
-extern void fftw_circular_convolution(FFTW_Workspace &ws, double * src, double * kernel);
-
-
-
-
-
#endif // CONVOLVE_H
diff --git a/Core/src/Convolve.cpp b/Core/src/Convolve.cpp
index fa2e7c70bc7..86a37ec4958 100644
--- a/Core/src/Convolve.cpp
+++ b/Core/src/Convolve.cpp
@@ -1,6 +1,7 @@
#include "Convolve.h"
#include <iostream>
#include <stdexcept>
+#include <sstream>
#include "Exceptions.h"
MathFunctions::Convolve::Convolve() : m_mode(FFTW_UNDEFINED)
@@ -21,6 +22,7 @@ MathFunctions::Convolve::~Convolve()
MathFunctions::Convolve::Workspace::Workspace() :
h_src(0), w_src(0), h_kernel(0), w_kernel(0),
w_fftw(0), h_fftw(0), in_src(0), out_src(0), in_kernel(0), out_kernel(0), dst_fft(0), h_dst(0), w_dst(0), dst(0),
+ h_offset(0), w_offset(0),
p_forw_src(NULL), p_forw_kernel(NULL), p_back(NULL)
{
@@ -57,6 +59,9 @@ void MathFunctions::Convolve::Workspace::clear()
if(dst) delete[] dst;
dst=0;
+ h_offset = 0;
+ w_offset = 0;
+
if(p_forw_src != NULL) fftw_destroy_plan(p_forw_src);
if(p_forw_kernel != NULL) fftw_destroy_plan(p_forw_kernel);
if(p_back != NULL) fftw_destroy_plan(p_back);
@@ -67,13 +72,13 @@ void MathFunctions::Convolve::Workspace::clear()
/* ************************************************************************* */
-// convolution in dd
+// convolution in 2d
/* ************************************************************************* */
void MathFunctions::Convolve::fftconvolve(const double2d_t &source, const double2d_t &kernel, double2d_t &result)
{
// set default convolution mode, if not defined
if(m_mode == FFTW_UNDEFINED) {
- std::cout << "setting FFTW_LINEAR_SAME" << std::endl;
+ std::cout << "MathFunctions::Convolve::fftconvolve() -> Info. Setting FFTW_LINEAR_SAME convolution mode." << std::endl;
setMode(FFTW_LINEAR_SAME);
}
@@ -82,32 +87,10 @@ void MathFunctions::Convolve::fftconvolve(const double2d_t &source, const double
size_t h_kernel = kernel.size();
size_t w_kernel = (kernel.size() ? kernel[0].size() : 0);
- if(!h_src || !w_src || !h_kernel || !w_kernel) {
- std::cout << "MathFunctions::Convolve::fftconvolve() -> Panic! Wrong dimensions " << h_src << " " << w_src << " " << h_kernel << " " << w_kernel << std::endl;
- throw std::runtime_error("XXX");
- }
-
+ // initialisation
init(h_src, w_src, h_kernel, w_kernel);
-
-// double *src = new double[h_src*w_src];
-// for(size_t i=0; i<h_src; i++) {
-// for(size_t j=0; j<w_src; j++) {
-// src[i*w_src + j] = source[i][j];
-// }
-// }
-
-// double *krn = new double[h_kernel*w_kernel];
-// for(size_t i=0; i<h_kernel; i++) {
-// for(size_t j=0; j<w_kernel; j++) {
-// krn[i*w_kernel + j] = kernel[i][j];
-// }
-// }
-
-// fftw_convolve(src, krn);
- fftw_convolve(source, kernel);
-
-// delete [] src;
-// delete [] krn;
+ // Compute the circular convolution
+ fftw_circular_convolution(source, kernel);
// results
result.clear();
@@ -115,7 +98,12 @@ void MathFunctions::Convolve::fftconvolve(const double2d_t &source, const double
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];
+ //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];
+ }
}
}
@@ -147,6 +135,12 @@ void MathFunctions::Convolve::fftconvolve(const double1d_t &source, const double
/* ************************************************************************* */
void MathFunctions::Convolve::init(int h_src, int w_src, int h_kernel, int w_kernel)
{
+ if(!h_src || !w_src || !h_kernel || !w_kernel) {
+ std::ostringstream os;
+ os << "MathFunctions::Convolve::init() -> Panic! Wrong dimensions " << h_src << " " << w_src << " " << h_kernel << " " << w_kernel << std::endl;
+ throw RuntimeErrorException(os.str());
+ }
+
ws.clear();
ws.h_src = h_src;
ws.w_src = w_src;
@@ -160,6 +154,9 @@ void MathFunctions::Convolve::init(int h_src, int w_src, int h_kernel, int w_ker
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
@@ -167,6 +164,9 @@ void MathFunctions::Convolve::init(int h_src, int w_src, int h_kernel, int w_ker
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
@@ -174,6 +174,9 @@ void MathFunctions::Convolve::init(int h_src, int w_src, int h_kernel, int w_ker
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
@@ -191,6 +194,9 @@ void MathFunctions::Convolve::init(int h_src, int w_src, int h_kernel, int w_ker
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:
// Circular convolution, modulo N
@@ -201,6 +207,9 @@ void MathFunctions::Convolve::init(int h_src, int w_src, int h_kernel, int w_ker
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;
case FFTW_CIRCULAR_SAME_SHIFTED:
// Circular convolution, modulo N, shifted by M/2
@@ -211,6 +220,9 @@ void MathFunctions::Convolve::init(int h_src, int w_src, int h_kernel, int w_ker
ws.w_fftw = w_src;
ws.h_dst = h_src;
ws.w_dst = w_src;
+ // We copy the [h_kernel/2:h_kernel/2+h_dst-1 ; w_kernel/2:w_kernel/2+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;
@@ -242,9 +254,13 @@ void MathFunctions::Convolve::init(int h_src, int w_src, int h_kernel, int w_ker
/* ************************************************************************* */
// initialise input and output arrays for fast fourier transformation
/* ************************************************************************* */
-//void MathFunctions::Convolve::fftw_circular_convolution(double * src, double * kernel)
void MathFunctions::Convolve::fftw_circular_convolution(const double2d_t &src, const double2d_t &kernel)
{
+ if(ws.h_fftw <= 0 || ws.w_fftw <= 0)
+ {
+ throw RuntimeErrorException("MathFunctions::Convolve::fftw_convolve() -> Panic! Initialisation is missed.");
+ }
+
double * ptr, *ptr_end, *ptr2;
// Reset the content of ws.in_src
@@ -294,74 +310,69 @@ void MathFunctions::Convolve::fftw_circular_convolution(const double2d_t &src, c
-/* ************************************************************************* */
-// initialise input and output arrays for fast fourier transformation
-/* ************************************************************************* */
-//void MathFunctions::Convolve::fftw_convolve(double * src,double * kernel)
-void MathFunctions::Convolve::fftw_convolve(const double2d_t &src, const double2d_t &kernel)
-{
- if(ws.h_fftw <= 0 || ws.w_fftw <= 0)
- {
- std::cout << "MathFunctions::Convolve::fftw_convolve() -> Panic! Initialisation is missed." << std::endl;
- throw std::runtime_error("MathFunctions::Convolve::fftw_convolve() -> Panic! Initialisation is missed.");
- }
-
- // Compute the circular convolution
- fftw_circular_convolution(src, kernel);
-
- // Depending on the type of convolution one is looking for, we extract the appropriate part of the result from out_src
- int h_offset(0), w_offset(0);
-
- switch(m_mode)
- {
- case FFTW_LINEAR_FULL:
- // Full Linear convolution
- // Here we just keep the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[i*ws.w_fftw], ws.w_dst*sizeof(double));
- break;
- case FFTW_LINEAR_SAME_UNPADDED:
- // Same linear convolution
- // 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
- h_offset = int(ws.h_kernel/2.0);
- w_offset = int(ws.w_kernel/2.0);
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[(i+h_offset)*ws.w_fftw+w_offset], ws.w_dst*sizeof(double));
- break;
- case FFTW_LINEAR_SAME:
- // Same linear convolution
- // 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
- h_offset = int(ws.h_kernel/2.0);
- w_offset = int(ws.w_kernel/2.0);
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[(i+h_offset)*ws.w_fftw+w_offset], ws.w_dst*sizeof(double));
- break;
- case FFTW_LINEAR_VALID:
- // Valid linear convolution
- // Here we just take [h_dst x w_dst] elements starting at [h_kernel-1;w_kernel-1]
- h_offset = ws.h_kernel - 1;
- w_offset = ws.w_kernel - 1;
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[(i+h_offset)*ws.w_fftw+w_offset], ws.w_dst*sizeof(double));
- break;
- case FFTW_CIRCULAR_SAME:
- // Circular convolution
- // We copy the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[i*ws.w_fftw], ws.w_dst*sizeof(double) );
- break;
- case FFTW_CIRCULAR_SAME_SHIFTED:
- // Shifted Circular convolution
- // We copy the [h_kernel/2:h_kernel/2+h_dst-1 ; w_kernel/2:w_kernel/2+w_dst-1] real part elements of out_src
- for(int i = 0 ; i < ws.h_dst ; ++i)
- for(int j = 0 ; j < ws.w_dst ; ++j)
- ws.dst[i*ws.w_dst + 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];
- 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;
- }
-}
-
+//void MathFunctions::Convolve::fftw_convolve(const double2d_t &src, const double2d_t &kernel)
+//{
+// if(ws.h_fftw <= 0 || ws.w_fftw <= 0)
+// {
+// std::cout << "MathFunctions::Convolve::fftw_convolve() -> Panic! Initialisation is missed." << std::endl;
+// throw std::runtime_error("MathFunctions::Convolve::fftw_convolve() -> Panic! Initialisation is missed.");
+// }
+//
+// // Compute the circular convolution
+// fftw_circular_convolution(src, kernel);
+//
+// // Depending on the type of convolution one is looking for, we extract the appropriate part of the result from out_src
+// int h_offset(0), w_offset(0);
+//
+// switch(m_mode)
+// {
+// case FFTW_LINEAR_FULL:
+// // Full Linear convolution
+// // Here we just keep the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
+// for(int i = 0 ; i < ws.h_dst ; ++i)
+// memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[i*ws.w_fftw], ws.w_dst*sizeof(double));
+// break;
+// case FFTW_LINEAR_SAME_UNPADDED:
+// // Same linear convolution
+// // 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
+// h_offset = int(ws.h_kernel/2.0);
+// w_offset = int(ws.w_kernel/2.0);
+// for(int i = 0 ; i < ws.h_dst ; ++i)
+// memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[(i+h_offset)*ws.w_fftw+w_offset], ws.w_dst*sizeof(double));
+// break;
+// case FFTW_LINEAR_SAME:
+// // Same linear convolution
+// // 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
+// h_offset = int(ws.h_kernel/2.0);
+// w_offset = int(ws.w_kernel/2.0);
+// for(int i = 0 ; i < ws.h_dst ; ++i)
+// memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[(i+h_offset)*ws.w_fftw+w_offset], ws.w_dst*sizeof(double));
+// break;
+// case FFTW_LINEAR_VALID:
+// // Valid linear convolution
+// // Here we just take [h_dst x w_dst] elements starting at [h_kernel-1;w_kernel-1]
+// h_offset = ws.h_kernel - 1;
+// w_offset = ws.w_kernel - 1;
+// for(int i = 0 ; i < ws.h_dst ; ++i)
+// memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[(i+h_offset)*ws.w_fftw+w_offset], ws.w_dst*sizeof(double));
+// break;
+// case FFTW_CIRCULAR_SAME:
+// // Circular convolution
+// // We copy the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
+// for(int i = 0 ; i < ws.h_dst ; ++i)
+// memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[i*ws.w_fftw], ws.w_dst*sizeof(double) );
+// break;
+// case FFTW_CIRCULAR_SAME_SHIFTED:
+// // Shifted Circular convolution
+// // We copy the [h_kernel/2:h_kernel/2+h_dst-1 ; w_kernel/2:w_kernel/2+w_dst-1] real part elements of out_src
+// for(int i = 0 ; i < ws.h_dst ; ++i)
+// for(int j = 0 ; j < ws.w_dst ; ++j)
+// ws.dst[i*ws.w_dst + 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];
+// 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;
+// }
+//}
/* ************************************************************************* */
@@ -397,354 +408,3 @@ bool MathFunctions::Convolve::is_optimal(int n)
}
-
-
-
-/////////////////
-///
-///
-
-
-// ******************** Begin of factorization code ***********************//
-// A piece of code to determine if a number "n" can be written as products of only the integers given in implemented_factors
-void factorize (const int n,
- int *n_factors,
- int factors[],
- int * implemented_factors)
-{
- int nf = 0;
- int ntest = n;
- int factor;
- int i = 0;
-
- if (n == 0)
- {
- printf("Length n must be positive integer\n");
- return ;
- }
-
- if (n == 1)
- {
- factors[0] = 1;
- *n_factors = 1;
- return ;
- }
-
- /* deal with the implemented factors */
-
- while (implemented_factors[i] && ntest != 1)
- {
- factor = implemented_factors[i];
- while ((ntest % factor) == 0)
- {
- ntest = ntest / factor;
- factors[nf] = factor;
- nf++;
- }
- i++;
- }
-
- // Ok that's it
- if(ntest != 1)
- {
- factors[nf] = ntest;
- nf++;
- }
-
- /* check that the factorization is correct */
- {
- int product = 1;
-
- for (i = 0; i < nf; i++)
- {
- product *= factors[i];
- }
-
- if (product != n)
- {
- printf("factorization failed");
- }
- }
-
- *n_factors = nf;
-}
-
-bool is_optimal(int n, int * implemented_factors)
-{
- int nf;
- int factors[64];
- int i = 0;
- factorize(n, &nf, factors,implemented_factors);
-
- // We just have to check if the last factor belongs to GSL_FACTORS
- while(implemented_factors[i])
- {
- if(factors[nf-1] == implemented_factors[i])
- return true;
- ++i;
- }
- return false;
-}
-
-int find_closest_factor(int n, int * implemented_factor)
-{
- int j;
- if(is_optimal(n,implemented_factor))
- return n;
- else
- {
- j = n+1;
- while(!is_optimal(j,implemented_factor))
- ++j;
- return j;
- }
-}
-// ******************** End of factorization code ***********************//
-
-int FFTW_FACTORS[7] = {13,11,7,5,3,2,0}; // end with zero to detect the end of the array
-
-
-void init_workspace_fftw(FFTW_Workspace & ws, FFTW_Convolution_Mode mode, int h_src, int w_src, int h_kernel, int w_kernel)
-{
- ws.h_src = h_src;
- ws.w_src = w_src;
- ws.h_kernel = h_kernel;
- ws.w_kernel = w_kernel;
- ws.mode = mode;
-
- switch(mode)
- {
- case FFTW_LINEAR_FULL:
- // Full Linear convolution
- ws.h_fftw = find_closest_factor(h_src + h_kernel - 1,FFTW_FACTORS);
- ws.w_fftw = find_closest_factor(w_src + w_kernel - 1,FFTW_FACTORS);
- ws.h_dst = h_src + h_kernel-1;
- ws.w_dst = w_src + w_kernel-1;
- 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;
- break;
- case FFTW_LINEAR_SAME:
- // Same Linear convolution
- ws.h_fftw = find_closest_factor(h_src + int(h_kernel/2.0),FFTW_FACTORS);
- ws.w_fftw = find_closest_factor(w_src + int(w_kernel/2.0),FFTW_FACTORS);
- ws.h_dst = h_src;
- ws.w_dst = w_src;
- break;
- case FFTW_LINEAR_VALID:
- // Valid Linear convolution
- if(ws.h_kernel > ws.h_src || ws.w_kernel > ws.w_src)
- {
- printf("Warning : The 'valid' convolution results in an empty matrix\n");
- ws.h_fftw = 0;
- ws.w_fftw = 0;
- ws.h_dst = 0;
- ws.w_dst = 0;
- }
- else
- {
- ws.h_fftw = find_closest_factor(h_src, FFTW_FACTORS);
- ws.w_fftw = find_closest_factor(w_src, FFTW_FACTORS);
- ws.h_dst = h_src - h_kernel+1;
- ws.w_dst = w_src - w_kernel+1;
- }
- break;
- case FFTW_CIRCULAR_SAME:
- // 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;
- break;
- case FFTW_CIRCULAR_SAME_SHIFTED:
- // Circular convolution, modulo N, shifted by M/2
- // 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;
- break;
- case FFTW_CIRCULAR_CUSTOM:
- // We here want to compute a circular convolution modulo h_dst, w_dst
- // These two variables must have been set before calling init_workscape !!
- ws.h_fftw = ws.h_dst;
- ws.w_fftw = ws.w_dst;
- break;
- default:
- printf("Unrecognized convolution mode, possible modes are :\n");
- printf(" - FFTW_LINEAR_FULL \n");
- printf(" - FFTW_LINEAR_SAME \n");
- printf(" - FFTW_LINEAR_SAME_UNPADDED\n");
- printf(" - FFTW_LINEAR_VALID \n");
- printf(" - FFTW_CIRCULAR_SAME \n");
- printf(" - FFTW_CIRCULAR_SHIFTED\n");
- printf(" - FFTW_CIRCULAR_CUSTOM\n");
- }
-
- 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);
- 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);
-
- // 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);
-}
-
-void clear_workspace_fftw(FFTW_Workspace & ws)
-{
- delete[] ws.in_src;
- fftw_free((fftw_complex*)ws.out_src);
- delete[] ws.in_kernel;
- fftw_free((fftw_complex*)ws.out_kernel);
-
- delete[] ws.dst_fft;
- delete[] ws.dst;
-
- // Destroy the plans
- fftw_destroy_plan(ws.p_forw_src);
- fftw_destroy_plan(ws.p_forw_kernel);
- fftw_destroy_plan(ws.p_back);
-}
-
-// Compute the circular convolution of src and kernel modulo ws.h_fftw, ws.w_fftw
-// using the Fast Fourier Transform
-// The result is in ws.dst
-void fftw_circular_convolution(FFTW_Workspace &ws, double * src, double * kernel)
-{
- 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*ws.w_src + 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*ws.w_kernel + 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);
-
- // That's it !
-}
-
-void fftw_convolve(FFTW_Workspace &ws, double * src,double * kernel)
-{
- if(ws.h_fftw <= 0 || ws.w_fftw <= 0)
- return;
-
- // Compute the circular convolution
- fftw_circular_convolution(ws, src, kernel);
-
- // Depending on the type of convolution one is looking for, we extract the appropriate part of the result from out_src
- int h_offset, w_offset;
-
-// double * dst_ptr;
-// double * src_ptr;
- switch(ws.mode)
- {
- case FFTW_LINEAR_FULL:
- // Full Linear convolution
- // Here we just keep the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[i*ws.w_fftw], ws.w_dst*sizeof(double));
- break;
- case FFTW_LINEAR_SAME_UNPADDED:
- // Same linear convolution
- // 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
- h_offset = int(ws.h_kernel/2.0);
- w_offset = int(ws.w_kernel/2.0);
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[(i+h_offset)*ws.w_fftw+w_offset], ws.w_dst*sizeof(double));
- //for(int j = 0 ; j < ws.w_dst ; ++j)
- // ws.dst[i*ws.w_dst + j] = ws.out_src[2*((i+h_offset)*ws.w_fftw+j+w_offset)+0]/double(ws.w_fftw * ws.h_fftw);
- break;
- case FFTW_LINEAR_SAME:
- // Same linear convolution
- // 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
- h_offset = int(ws.h_kernel/2.0);
- w_offset = int(ws.w_kernel/2.0);
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[(i+h_offset)*ws.w_fftw+w_offset], ws.w_dst*sizeof(double));
- break;
- case FFTW_LINEAR_VALID:
- // Valid linear convolution
- // Here we just take [h_dst x w_dst] elements starting at [h_kernel-1;w_kernel-1]
- h_offset = ws.h_kernel - 1;
- w_offset = ws.w_kernel - 1;
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[(i+h_offset)*ws.w_fftw+w_offset], ws.w_dst*sizeof(double));
- break;
- case FFTW_CIRCULAR_SAME:
- // Circular convolution
- // We copy the first [0:h_dst-1 ; 0:w_dst-1] real part elements of out_src
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[i*ws.w_fftw], ws.w_dst*sizeof(double) );
- break;
- case FFTW_CIRCULAR_SAME_SHIFTED:
- // Shifted Circular convolution
- // We copy the [h_kernel/2:h_kernel/2+h_dst-1 ; w_kernel/2:w_kernel/2+w_dst-1] real part elements of out_src
- for(int i = 0 ; i < ws.h_dst ; ++i)
- for(int j = 0 ; j < ws.w_dst ; ++j)
- ws.dst[i*ws.w_dst + 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];
- break;
- case FFTW_CIRCULAR_CUSTOM:
- for(int i = 0 ; i < ws.h_dst ; ++i)
- memcpy(&ws.dst[i*ws.w_dst], &ws.dst_fft[i*ws.w_fftw], ws.w_dst*sizeof(double) );
- break;
- default:
- printf("Unrecognized convolution mode, possible modes are :\n");
- printf(" - FFTW_LINEAR_FULL \n");
- printf(" - FFTW_LINEAR_SAME \n");
- printf(" - FFTW_LINEAR_SAME_UNPADDED\n");
- printf(" - FFTW_LINEAR_VALID \n");
- printf(" - FFTW_CIRCULAR_SAME \n");
- printf(" - FFTW_CIRCULAR_SHIFTED\n");
- printf(" - FFTW_CIRCULAR_CUSTOM\n");
- }
-}
-
--
GitLab