From 1570eb840f238b61914d8ac05d61826b4a1b54c5 Mon Sep 17 00:00:00 2001
From: pospelov <pospelov@fz-juelich.de>
Date: Fri, 8 Jun 2012 18:05:16 +0200
Subject: [PATCH] new Convolve class for fast fourier transform convolution in
2d
---
App/App.pro | 6 +-
App/inc/TestInstrument.h | 45 ++
App/src/DrawHelper.cpp | 7 +-
App/src/TestDiffuseReflection.cpp | 11 +-
App/src/TestFactory.cpp | 2 +
App/src/TestFresnelCoeff.cpp | 4 +-
App/src/TestInstrument.cpp | 364 +++++++++++++++
Core/Core.pro | 6 +-
Core/inc/Convolve.h | 162 +++++++
Core/inc/Exceptions.h | 6 +
Core/inc/MathFunctions.h | 2 +
Core/src/Convolve.cpp | 750 ++++++++++++++++++++++++++++++
Core/src/Exceptions.cpp | 6 +
Core/src/IFunctionalTest.cpp | 2 +-
Core/src/MathFunctions.cpp | 20 +
GISASFW.pro | 3 +-
16 files changed, 1380 insertions(+), 16 deletions(-)
create mode 100644 App/inc/TestInstrument.h
create mode 100644 App/src/TestInstrument.cpp
create mode 100644 Core/inc/Convolve.h
create mode 100644 Core/src/Convolve.cpp
diff --git a/App/App.pro b/App/App.pro
index e3aac81a099..5e1188a8c84 100644
--- a/App/App.pro
+++ b/App/App.pro
@@ -16,7 +16,8 @@ SOURCES += \
src/StandardSamples.cpp \
src/CommandLine.cpp \
src/TestFactory.cpp \
- src/TestDiffuseReflection.cpp
+ src/TestDiffuseReflection.cpp \
+ src/TestInstrument.cpp
HEADERS += \
inc/DrawHelper.h \
@@ -30,7 +31,8 @@ HEADERS += \
inc/StandardSamples.h \
inc/CommandLine.h \
inc/TestFactory.h \
- inc/TestDiffuseReflection.h
+ inc/TestDiffuseReflection.h \
+ inc/TestInstrument.h
INCLUDEPATH += ./inc
DEPENDPATH += ./inc
diff --git a/App/inc/TestInstrument.h b/App/inc/TestInstrument.h
new file mode 100644
index 00000000000..79073497c70
--- /dev/null
+++ b/App/inc/TestInstrument.h
@@ -0,0 +1,45 @@
+#ifndef TESTINSTRUMENT_H
+#define TESTINSTRUMENT_H
+// ********************************************************************
+// * The BornAgain project *
+// * Simulation of neutron and x-ray scattering at grazing incidence *
+// * *
+// * LICENSE AND DISCLAIMER *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit. Mauris *
+// * eget quam orci. Quisque porta varius dui, quis posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor. *
+// ********************************************************************
+//! @file TestConvolution.h
+//! @brief Definition of TestConvolution class for testing of Convolve class
+//! @author Scientific Computing Group at FRM II
+//! @date 02.05.2012
+
+#include "IFunctionalTest.h"
+#include "Convolve.h"
+#include <string>
+#include <vector>
+
+
+//- -------------------------------------------------------------------
+//! @class TestConvolution
+//! @brief Testing Convolve class for instrumental effects studies
+//- -------------------------------------------------------------------
+class TestInstrument : public IFunctionalTest
+{
+public:
+ TestInstrument();
+
+ void execute();
+
+ //! testing convolution in 1d
+ void test_convolve1d();
+
+ //! testing convolution in 2d
+ void test_convolve2d();
+
+private:
+ typedef std::pair<std::string, MathFunctions::Convolve::Mode> mode_pair_t;
+ std::vector<mode_pair_t> m_modes;
+};
+
+#endif // TESTINSTRUMENT_H
diff --git a/App/src/DrawHelper.cpp b/App/src/DrawHelper.cpp
index ffbfaba2b22..132e1b8fe0c 100644
--- a/App/src/DrawHelper.cpp
+++ b/App/src/DrawHelper.cpp
@@ -162,12 +162,13 @@ void DrawHelper::SetStyle()
scattStyle->SetMarkerSize(0.2);
scattStyle->SetHistLineWidth(2.);
scattStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes
+ scattStyle->SetFuncWidth(1.);
// get rid of X error bars and y error bar caps
//scattStyle->SetErrorX(0.001);
// do not display any of the standard histogram decorations
- scattStyle->SetOptTitle(0);
+ scattStyle->SetOptTitle(1);
//scattStyle->SetOptStat(1111);
scattStyle->SetOptStat(1);
//scattStyle->SetOptFit(1111);
@@ -266,7 +267,7 @@ void DrawHelper::DrawMultilayer(const MultiLayer *sample)
double z2 = sample->getLayerBottomZ(nlayers-1);
double size = std::abs(z2 - z1); // typical size of sample
- double margin = size*0.02; // safety margin to avoid overlapping of layers
+ double margin = size*0.01; // safety margin to avoid overlapping of layers
double fake_thickness = size*0.3; // thickness for layers without defined thickness (like ambience&substrate)
// frame to draw
@@ -277,7 +278,7 @@ void DrawHelper::DrawMultilayer(const MultiLayer *sample)
// drawing properties for interfaces, material and their names
TLine interface;
- interface.SetLineWidth(3);
+ interface.SetLineWidth(1);
interface.SetLineStyle(3);
interface.SetLineColor(kRed);
TBox bulk;
diff --git a/App/src/TestDiffuseReflection.cpp b/App/src/TestDiffuseReflection.cpp
index c96900e5042..ddb74fcbb40 100644
--- a/App/src/TestDiffuseReflection.cpp
+++ b/App/src/TestDiffuseReflection.cpp
@@ -35,7 +35,7 @@ void TestDiffuseReflection::execute()
DWBADiffuseReflection calc;
std::vector<std::string > snames;
- //snames.push_back("MultilayerOffspecTestcase1a");
+ snames.push_back("MultilayerOffspecTestcase1a");
snames.push_back("MultilayerOffspecTestcase1b");
//snames.push_back("MultilayerOffspecTestcase2a");
//snames.push_back("MultilayerOffspecTestcase2b");
@@ -162,11 +162,12 @@ void TestDiffuseReflection::draw()
TCanvas *c1 = new TCanvas(cname.c_str(),"Diffuse reflection",1024,768);
DrawHelper::instance().SetMagnifier(c1);
- c1->Divide(2,2);
- c1->cd(1);
- gr->Draw("apl");
+// c1->Divide(2,2);
+// c1->cd(1);
+// gr->Draw("apl");
- c1->cd(2);
+// c1->cd(2);
+ c1->cd();
gPad->SetRightMargin(0.2);
gStyle->SetPalette(1);
gPad->SetLogz();
diff --git a/App/src/TestFactory.cpp b/App/src/TestFactory.cpp
index c4cf7dad00a..97b6cf7705f 100644
--- a/App/src/TestFactory.cpp
+++ b/App/src/TestFactory.cpp
@@ -4,6 +4,7 @@
#include "TestFormFactor.h"
#include "TestDWBAFormFactor.h"
#include "TestDiffuseReflection.h"
+#include "TestInstrument.h"
#include "TBenchmark.h"
@@ -19,6 +20,7 @@ TestFactory::TestFactory() : m_benchmark(0)
registerItem("formfactor", IFactoryCreateFunction<TestFormFactor, IFunctionalTest> );
registerItem("dwba", IFactoryCreateFunction<TestDWBAFormFactor, IFunctionalTest> );
registerItem("diffuse", IFactoryCreateFunction<TestDiffuseReflection, IFunctionalTest> );
+ registerItem("instrument", IFactoryCreateFunction<TestInstrument, IFunctionalTest> );
m_benchmark = new TBenchmark();
}
diff --git a/App/src/TestFresnelCoeff.cpp b/App/src/TestFresnelCoeff.cpp
index efe08ef1435..58778877c4d 100644
--- a/App/src/TestFresnelCoeff.cpp
+++ b/App/src/TestFresnelCoeff.cpp
@@ -35,10 +35,10 @@ void TestFresnelCoeff::execute()
std::cout << "TestFresnelCoeff::execute() -> Info." << std::endl;
std::vector<std::string > snames;
-// snames.push_back("AirOnSubstrate");
+ snames.push_back("AirOnSubstrate");
// snames.push_back("SubstrateOnSubstrate");
// snames.push_back("SimpleMultilayer");
- snames.push_back("MultilayerOffspecTestcase1a");
+// snames.push_back("MultilayerOffspecTestcase1a");
// loop over standard samples defined in SampleFactory and StandardSamples
for(size_t i_sample=0; i_sample<snames.size(); i_sample++){
diff --git a/App/src/TestInstrument.cpp b/App/src/TestInstrument.cpp
new file mode 100644
index 00000000000..98923c866ad
--- /dev/null
+++ b/App/src/TestInstrument.cpp
@@ -0,0 +1,364 @@
+#include "TestInstrument.h"
+#include "MathFunctions.h"
+#include "Convolve.h"
+#include "DrawHelper.h"
+
+#include <iostream>
+
+#include "TMath.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TF1.h"
+#include "TF2.h"
+#include "TRandom.h"
+#include "TCanvas.h"
+#include "TProfile.h"
+#include "TGraph.h"
+#include "TBenchmark.h"
+#include "TStyle.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()
+{
+ std::cout << "TestInstrument::TestInstrument() -> 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));
+ m_modes.push_back( mode_pair_t("LINEAR_SAME", MathFunctions::Convolve::FFTW_LINEAR_SAME));
+ m_modes.push_back( mode_pair_t("CIRCULAR_SAME", MathFunctions::Convolve::FFTW_CIRCULAR_SAME));
+ m_modes.push_back( mode_pair_t("CIRCULAR_SAME_SHIFTED", MathFunctions::Convolve::FFTW_CIRCULAR_SAME_SHIFTED));
+}
+
+
+void TestInstrument::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
+
+}
+
+
+
+void TestInstrument::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;
+
+ 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];
+ 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();
+ }
+ 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);
+ 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;
+ }
+
+ // building responce function
+ std::vector<double> xresp;
+ xresp.resize(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];
+// }
+
+ 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;
+ }
+
+ 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);
+ href->DrawCopy();
+ gr_resp->Draw("pl same");
+
+
+ TBenchmark benchmark;
+ for(size_t i_mode=0; i_mode<m_modes.size(); i_mode++) {
+ std::string sname = m_modes[i_mode].first;
+ MathFunctions::Convolve::Mode mode = m_modes[i_mode].second;
+
+ // result
+ TGraph *gr_result = new TGraph(npoints);
+ benchmark.Start(sname.c_str());
+
+ MathFunctions::Convolve cv;
+ cv.setMode(mode);
+ std::vector<double> xresult;
+ 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);
+ }
+ benchmark.Stop(sname.c_str());
+ benchmark.Show(sname.c_str());
+ for(size_t i=0; i<xresult.size(); i++) {
+ double x = xmin + i*dx;
+ gr_result->SetPoint(i,x,xresult[i]);
+ }
+ c1->cd(4+i_mode);
+ href->SetMaximum(20.0);
+ href->DrawCopy();
+ gr_result->Draw("pl same");
+ }
+ float_t rp, cp;
+ std::cout << "--------------" << std::endl;
+ benchmark.Summary(rp, cp);
+
+}
+
+
+/* ************************************************************************* */
+//
+/* ************************************************************************* */
+void TestInstrument::test_convolve2d()
+{
+ double xmin(0.0), xmax(100.);
+ size_t npx = 100;
+ double dx = (xmax-xmin)/(npx);
+
+ double ymin(0.0), ymax(40.);
+ size_t npy = 40;
+ double dy = (ymax-ymin)/(npy);
+
+ double sigmax=(xmax-xmin)*0.01;
+ double sigmay=(ymax-ymin)*0.01;
+
+ // filling signal pattern - spikes on top of flat background
+ std::vector<std::vector<double> > source;
+ source.resize(npx);
+ for(size_t ix=0; ix<npx; ix++) {
+ source[ix].resize(npy, 1);
+ for(size_t iy=0; iy<npy; iy++){
+ // background
+ source[ix][iy] = 0.5*double(ix)/double(npx) + 0.5*double(iy)/double(npy);
+ // 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
+ TH2F *h2_signal = new TH2F("h2_signal","h2_signal", npx, xmin, xmax, npy, ymin, ymax);
+ h2_signal->SetStats(0);
+ h2_signal->SetMinimum(0);
+ h2_signal->SetMaximum(3.0);
+ for(size_t ix=0; ix<npx; ix++){
+ double x = xmin + dx*ix;
+ for(size_t iy=0; iy<npy; iy++){
+ double y = ymin + dy*iy;
+ h2_signal->Fill(x,y,source[ix][iy]);
+ }
+ }
+
+ // resolution function (two dimensional gaus)
+ TF2 *fres = new TF2("fres","[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])", xmin, xmax, ymin, ymax);
+ fres->SetParameters(1,xmin+(xmax-xmin)/2.,sigmax, ymin + (ymax-ymin)/2., sigmay);
+ // 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);
+ kernel.resize(npx);
+ for(size_t ix=0; ix<npx; ix++) {
+ kernel[ix].resize(npy);
+ double x = xmin + dx*ix;
+ for(size_t iy=0; iy<npy; iy++){
+ double y = ymin + dy*iy;
+ kernel[ix][iy] = fres->Eval(x,y);
+ h2_kernel->Fill(x,y, kernel[ix][iy]);
+ }
+ }
+
+ // simulation of measurements with smearing measurement with resolution function
+ TH2F *h2_measured = new TH2F("h2_measured","h2_measured", npx, xmin, xmax, npy, ymin, ymax);
+ TRandom mr;
+ for(int i_meas=0; i_meas<1000000; i_meas++) {
+ // generating distribution according to our signal
+ double x,y;
+ h2_signal->GetRandom2(x,y);
+ // smearing signal to get "measured" signal
+ h2_measured->Fill(x + mr.Gaus(0,sigmax), y + mr.Gaus(0,sigmay));
+ }
+
+ // drawing
+ TCanvas *c1 = new TCanvas("c1_test_convolve2d","c1_test_convolve2d",1024, 768);
+ DrawHelper::instance().SetMagnifier(c1);
+ gStyle->SetPalette(1);
+
+ c1->Divide(3,3);
+ c1->cd(1); gPad->SetRightMargin(0.15);
+ h2_signal->Draw("colz");
+
+ c1->cd(2); gPad->SetRightMargin(0.15);
+ h2_signal->Draw("lego");
+
+ c1->cd(3); gPad->SetRightMargin(0.15);
+ h2_kernel->Draw("lego");
+
+ c1->cd(4); gPad->SetRightMargin(0.15);
+ h2_measured->Draw("colz");
+
+ // calculation of convolution
+ std::vector<std::vector<double> > result;
+ TBenchmark benchmark;
+ for(size_t i_mode=0; i_mode<m_modes.size(); i_mode++) {
+ std::string sname = m_modes[i_mode].first;
+ MathFunctions::Convolve::Mode mode = m_modes[i_mode].second;
+
+ // running convolution several times to get statistics for benchmarking
+ benchmark.Start(sname.c_str());
+ MathFunctions::Convolve cv;
+ cv.setMode(mode);
+ for(int i=0; i<1000; i++) {
+ cv.fftconvolve(source, kernel, result);
+ }
+ benchmark.Stop(sname.c_str());
+
+ // drawing results of last convolution
+ c1->cd(5+i_mode);
+ 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++) {
+ double x = xmin + dx*ix;
+ for(size_t iy=0; iy<npy; iy++){
+ double y = ymin + dy*iy;
+ h2_result.Fill(x,y, result[ix][iy]);
+ }
+ }
+ h2_result.SetTitle(sname.c_str());
+ h2_result.DrawCopy("colz");
+ }
+
+ // benchmark summary
+ float_t rp, cp;
+ std::cout << "--------------" << std::endl;
+ benchmark.Summary(rp, cp);
+
+}
+
+
+
+
+std::vector<double> Convolve_MyFFT1D_B(int mode, const std::vector<double> &data1, const std::vector<double> &data2)
+{
+ size_t Ns = data1.size();
+ double * src = new double[Ns];
+ for(size_t i=0; i<Ns; i++) {
+ src[i] = data1[i];
+ }
+
+ 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;
+}
+
+
+
+
+
diff --git a/Core/Core.pro b/Core/Core.pro
index a447ecc215f..d82f519cbda 100644
--- a/Core/Core.pro
+++ b/Core/Core.pro
@@ -36,7 +36,8 @@ SOURCES += \
src/ISingleton.cpp \
src/IFunctionalTest.cpp \
src/IFactory.cpp \
- src/DWBADiffuseReflection.cpp
+ src/DWBADiffuseReflection.cpp \
+ src/Convolve.cpp
HEADERS += \
inc/ISample.h \
@@ -72,7 +73,8 @@ HEADERS += \
inc/IFunctionalTest.h \
inc/IFactory.h \
inc/Numeric.h \
- inc/DWBADiffuseReflection.h
+ inc/DWBADiffuseReflection.h \
+ inc/Convolve.h
INCLUDEPATH += ./inc
DEPENDPATH += ./inc
diff --git a/Core/inc/Convolve.h b/Core/inc/Convolve.h
new file mode 100644
index 00000000000..2ba520f64b6
--- /dev/null
+++ b/Core/inc/Convolve.h
@@ -0,0 +1,162 @@
+#ifndef CONVOLVE_H
+#define CONVOLVE_H
+// ********************************************************************
+// * The BornAgain project *
+// * Simulation of neutron and x-ray scattering at grazing incidence *
+// * *
+// * LICENSE AND DISCLAIMER *
+// * Lorem ipsum dolor sit amet, consectetur adipiscing elit. Mauris *
+// * eget quam orci. Quisque porta varius dui, quis posuere nibh *
+// * mollis quis. Mauris commodo rhoncus porttitor. *
+// ********************************************************************
+//! @file Convolve.h
+//! @brief Definition of Convolve class
+//! @author Scientific Computing Group at FRM II
+//! @date 30.05.2012
+
+
+#include <fftw3.h>
+#include <vector>
+
+
+
+namespace MathFunctions
+{
+
+
+//- -------------------------------------------------------------------
+//! @class Convolve
+//! @brief Convolution of two real vectors (1D or 2D ) using Fast Fourier Transformation.
+//!
+//! 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
+//- -------------------------------------------------------------------
+class Convolve
+{
+public:
+ Convolve();
+ ~Convolve();
+
+ //! definition of 1d vector of double
+ typedef std::vector<double> double1d_t;
+
+ //! definition of 2d vector of double
+ typedef std::vector<double1d_t> double2d_t;
+
+ //! convolution modes
+ 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
+ void fftconvolve(const double1d_t &source, const double1d_t &kernel, double1d_t &result);
+
+ //! convolution in 2D
+ void fftconvolve(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);
+
+ //! set convolution mode
+ 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
+ int find_closest_factor(int n);
+
+ //! if a number can be factorised using only favorite fftw3 factors
+ 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
+ //! Sizes of input arrays are adjusted; output arrays are alocated via fftw3 allocation for maximum performance
+ class Workspace
+ {
+ public:
+ Workspace();
+ ~Workspace();
+ void clear();
+ friend class Convolve;
+ private:
+ int h_src, w_src; // size of original 'source' array in 2 dimensions
+ int h_kernel, w_kernel; // size of original 'kernel' array in 2 dimensions
+ int w_fftw, h_fftw; // size of adjusted source and kernel arrays (in_src, out_src, in_kernel, out_kernel)
+ double *in_src; // adjusted input 'source' array
+ double *out_src; // result of fourier transformation of source
+ double *in_kernel; // adjusted input 'kernel' array
+ double *out_kernel; // result of fourier transformation of kernel
+ 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
+ fftw_plan p_forw_src;
+ fftw_plan p_forw_kernel;
+ fftw_plan p_back;
+ };
+
+ Workspace ws; // input and output data for fftw3
+ Mode m_mode; // convolution mode
+ std::vector<size_t > m_implemented_factors; // favorite factorization terms of fftw3
+};
+
+
+} // 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/inc/Exceptions.h b/Core/inc/Exceptions.h
index d35aba238c5..fd740afb836 100644
--- a/Core/inc/Exceptions.h
+++ b/Core/inc/Exceptions.h
@@ -72,6 +72,12 @@ public:
LogicErrorException(const std::string& message);
};
+class RuntimeErrorException : public std::runtime_error
+{
+public:
+ RuntimeErrorException(const std::string& message);
+};
+
class DivisionByZeroException : public std::runtime_error
{
public:
diff --git a/Core/inc/MathFunctions.h b/Core/inc/MathFunctions.h
index 6a205a81bd5..93f0003eff9 100644
--- a/Core/inc/MathFunctions.h
+++ b/Core/inc/MathFunctions.h
@@ -42,6 +42,8 @@ std::vector<complex_t > FastFourierTransform(const std::vector<complex_t > &data
std::vector<complex_t > FastFourierTransform(const std::vector<double > &data, TransformCase tcase);
+std::vector<complex_t> ConvolveFFT(const std::vector<double> &signal, const std::vector<double> &resfunc);
+
} // Namespace MathFunctions
inline double MathFunctions::GenerateNormalRandom(double average, double std_dev)
diff --git a/Core/src/Convolve.cpp b/Core/src/Convolve.cpp
new file mode 100644
index 00000000000..fa2e7c70bc7
--- /dev/null
+++ b/Core/src/Convolve.cpp
@@ -0,0 +1,750 @@
+#include "Convolve.h"
+#include <iostream>
+#include <stdexcept>
+#include "Exceptions.h"
+
+MathFunctions::Convolve::Convolve() : 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]));
+}
+
+
+
+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),
+ p_forw_src(NULL), p_forw_kernel(NULL), p_back(NULL)
+{
+
+}
+
+
+MathFunctions::Convolve::Workspace::~Workspace()
+{
+ clear();
+}
+
+
+void MathFunctions::Convolve::Workspace::clear()
+{
+ h_src=0;
+ w_src=0;
+ h_kernel=0;
+ w_kernel = 0;
+ if(in_src) delete[] in_src;
+ in_src = 0;
+
+ if(out_src) fftw_free((fftw_complex*)out_src);
+ out_src = 0;
+
+ if(in_kernel) delete[] in_kernel;
+ in_kernel=0;
+
+ if(out_kernel) fftw_free((fftw_complex*)out_kernel);
+ out_kernel=0;
+
+ if(dst_fft) delete[] dst_fft;
+ dst_fft=0;
+
+ if(dst) delete[] dst;
+ dst=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);
+
+ // this returns fftw3 into completely initial state but is dramatically slow
+ //fftw_cleanup();
+}
+
+
+/* ************************************************************************* */
+// convolution in dd
+/* ************************************************************************* */
+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;
+ setMode(FFTW_LINEAR_SAME);
+ }
+
+ size_t h_src = source.size();
+ size_t w_src = (source.size() ? source[0].size() : 0);
+ 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");
+ }
+
+ 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;
+
+ // 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];
+ }
+ }
+
+}
+
+
+/* ************************************************************************* */
+// convolution in 1d
+/* ************************************************************************* */
+void MathFunctions::Convolve::fftconvolve(const double1d_t &source, const double1d_t &kernel, double1d_t &result)
+{
+ // 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) {
+ std::cout << "MathFunctions::Convolve::fftconvolve -> Panic in 1d" << std::endl;
+ throw std::runtime_error("Panic!");
+ }
+ result = result2d[0];
+}
+
+
+/* ************************************************************************* */
+// initialise input and output arrays for fast fourier transformation
+/* ************************************************************************* */
+void MathFunctions::Convolve::init(int h_src, int w_src, int h_kernel, int w_kernel)
+{
+ 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;
+ 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));
+ ws.w_fftw = find_closest_factor(w_src + int(w_kernel/2.0));
+ 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)
+ {
+ 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;
+ }
+ 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;
+ 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;
+ }
+
+ 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 == NULL ) throw RuntimeErrorException("MathFunctions::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 == NULL ) throw RuntimeErrorException("MathFunctions::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 == NULL ) throw RuntimeErrorException("MathFunctions::Convolve::init() -> Error! Can't initialise p_back plan.");
+
+
+}
+
+
+/* ************************************************************************* */
+// 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)
+{
+ 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];
+ //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][j];
+ //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);
+
+}
+
+
+
+/* ************************************************************************* */
+// 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;
+ }
+}
+
+
+
+/* ************************************************************************* */
+// find a number closest to the given one, which can be factorised according
+// to fftw3 favorite factorisation
+/* ************************************************************************* */
+int MathFunctions::Convolve::find_closest_factor(int n)
+{
+ if(is_optimal(n) ) {
+ return n;
+ } else {
+ int j = n+1;
+ while( !is_optimal(j) ) ++j;
+ return j;
+ }
+}
+
+
+/* ************************************************************************* */
+// if a number can be factorised using only favorite fftw3 factors
+/* ************************************************************************* */
+bool MathFunctions::Convolve::is_optimal(int n)
+{
+ if(n==1) return false;
+ int 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];
+ }
+ }
+ if(ntest==1) return true;
+ return false;
+}
+
+
+
+
+
+/////////////////
+///
+///
+
+
+// ******************** 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");
+ }
+}
+
diff --git a/Core/src/Exceptions.cpp b/Core/src/Exceptions.cpp
index 72b7e5bb4a3..70350420c1c 100644
--- a/Core/src/Exceptions.cpp
+++ b/Core/src/Exceptions.cpp
@@ -60,6 +60,12 @@ LogicErrorException::LogicErrorException(const std::string &message)
LogExceptionMessage(message);
}
+RuntimeErrorException::RuntimeErrorException(const std::string &message)
+ : std::runtime_error(message)
+{
+ LogExceptionMessage(message);
+}
+
DivisionByZeroException::DivisionByZeroException(const std::string &message)
: std::runtime_error(message)
{
diff --git a/Core/src/IFunctionalTest.cpp b/Core/src/IFunctionalTest.cpp
index 8da8babfa99..8946990022e 100644
--- a/Core/src/IFunctionalTest.cpp
+++ b/Core/src/IFunctionalTest.cpp
@@ -8,6 +8,6 @@ IFunctionalTest::IFunctionalTest()
void IFunctionalTest::execute()
{
- throw NotImplementedException("This test can't run.");
+ throw NotImplementedException("This test can't run. Undefined execute method.");
}
diff --git a/Core/src/MathFunctions.cpp b/Core/src/MathFunctions.cpp
index f9961a4a986..f0cb5848bd0 100644
--- a/Core/src/MathFunctions.cpp
+++ b/Core/src/MathFunctions.cpp
@@ -88,3 +88,23 @@ std::vector<complex_t > MathFunctions::FastFourierTransform(const std::vector<do
}
+/* ************************************************************************* */
+// convolution of two real vectors of equal size
+/* ************************************************************************* */
+std::vector<complex_t> MathFunctions::ConvolveFFT(const std::vector<double> &signal, const std::vector<double> &resfunc)
+{
+ if(signal.size() != resfunc.size() ) {
+ throw std::runtime_error("MathFunctions::ConvolveFFT() -> This convolution works only for two vectors of equal size. Use Convolve class instead.");
+ }
+ std::vector<complex_t > fft_signal = MathFunctions::FastFourierTransform(signal, MathFunctions::ForwardFFT);
+ std::vector<complex_t > fft_resfunc = MathFunctions::FastFourierTransform(resfunc, MathFunctions::ForwardFFT);
+
+ std::vector<complex_t > fft_prod;
+ fft_prod.resize(fft_signal.size());
+ for(size_t i=0; i<fft_signal.size(); i++){
+ fft_prod[i] = fft_signal[i] * fft_resfunc[i];
+ }
+
+ std::vector<complex_t > result = MathFunctions::FastFourierTransform(fft_prod, MathFunctions::BackwardFFT);
+ return result;
+}
diff --git a/GISASFW.pro b/GISASFW.pro
index 9a8faa9233f..c2414637001 100644
--- a/GISASFW.pro
+++ b/GISASFW.pro
@@ -3,10 +3,11 @@ TEMPLATE = subdirs
SUBDIRS += \
Core \
ThirdParty/gtest \
+# ThirdParty/fftw++ \
App \
UnitTests/TestCore
-TestCore.depends = ThirdParty/gtest
+TestCore.depends = ThirdParty/gtest ThirdParty/fftw++
# means that compilation will be in the listed order
CONFIG += ordered
--
GitLab