From 048755bab6b77c1da0050ed091e5007a60564adf Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Tue, 8 Mar 2016 15:41:38 +0100
Subject: Use CompositeGeometryManager for FDK

Also fix a number of scaling/weighting issues in FDK, and
switch to standard cone_bp with FDKWeighting for the BP step.
---
 cuda/2d/fft.cu | 45 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)

(limited to 'cuda/2d/fft.cu')

diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu
index 2bfd493..a84b2cc 100644
--- a/cuda/2d/fft.cu
+++ b/cuda/2d/fft.cu
@@ -35,6 +35,7 @@ $Id$
 #include <fstream>
 
 #include "../../include/astra/Logging.h"
+#include "../../include/astra/Fourier.h"
 
 using namespace astra;
 
@@ -303,6 +304,8 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
 	float * pfFilt = new float[_iFFTFourierDetectorCount];
 	float * pfW = new float[_iFFTFourierDetectorCount];
 
+#if 1
+
 	for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
 	{
 		float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount;
@@ -314,6 +317,48 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
 		// w = 2*pi*(0:size(filt,2)-1)/order
 		pfW[iDetectorIndex] = 3.1415f * 2.0f * fRelIndex;
 	}
+#else
+
+	float *pfRealIn = new float[_iFFTRealDetectorCount];
+	float *pfImagIn = new float[_iFFTRealDetectorCount];
+	float *pfRealOut = new float[_iFFTRealDetectorCount];
+	float *pfImagOut = new float[_iFFTRealDetectorCount];
+
+	for (int i = 0; i < _iFFTRealDetectorCount; ++i) {
+		pfImagIn[i] = 0.0f;
+		pfRealOut[i] = 0.0f;
+		pfImagOut[i] = 0.0f;
+
+		if (i & 1) {
+			int j = i;
+			if (2*j > _iFFTRealDetectorCount)
+				j = _iFFTRealDetectorCount - j;
+			float f = M_PI * j;
+			pfRealIn[i] = -1 / (f*f);
+		} else {
+			pfRealIn[i] = 0.0f;
+		}
+	}
+
+	pfRealIn[0] = 0.25f;
+
+	fastTwoPowerFourierTransform1D(_iFFTRealDetectorCount, pfRealIn, pfImagIn,
+	                               pfRealOut, pfImagOut, 1, 1, false);
+
+	for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+	{
+		float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount;
+
+		pfFilt[iDetectorIndex] = 2.0f * pfRealOut[iDetectorIndex];
+		pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex;
+	}
+
+	delete[] pfRealIn;
+	delete[] pfImagIn;
+	delete[] pfRealOut;
+	delete[] pfImagOut;
+
+#endif
 
 	switch(_eFilter)
 	{
-- 
cgit v1.2.3


From 558820f1421e79b32e542c133c83683e58df6d5e Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 24 Jun 2016 15:28:47 +0200
Subject: Compute FBP filter in spatial domain

---
 cuda/2d/fft.cu | 46 +++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 39 insertions(+), 7 deletions(-)

(limited to 'cuda/2d/fft.cu')

diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu
index 2bfd493..2d259a9 100644
--- a/cuda/2d/fft.cu
+++ b/cuda/2d/fft.cu
@@ -35,6 +35,7 @@ $Id$
 #include <fstream>
 
 #include "../../include/astra/Logging.h"
+#include "astra/Fourier.h"
 
 using namespace astra;
 
@@ -303,16 +304,48 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
 	float * pfFilt = new float[_iFFTFourierDetectorCount];
 	float * pfW = new float[_iFFTFourierDetectorCount];
 
+	// We cache one Fourier transform for repeated FBP's of the same size
+	static float *pfData = 0;
+	static int iFilterCacheSize = 0;
+
+	if (!pfData || iFilterCacheSize != _iFFTRealDetectorCount) {
+		// Compute filter in spatial domain
+
+		delete[] pfData;
+		pfData = new float[2*_iFFTRealDetectorCount];
+		int *ip = new int[int(2+sqrt(_iFFTRealDetectorCount)+1)];
+		ip[0] = 0;
+		float32 *w = new float32[_iFFTRealDetectorCount/2];
+
+		for (int i = 0; i < _iFFTRealDetectorCount; ++i) {
+			pfData[2*i+1] = 0.0f;
+
+			if (i & 1) {
+				int j = i;
+				if (2*j > _iFFTRealDetectorCount)
+					j = _iFFTRealDetectorCount - j;
+				float f = M_PI * j;
+				pfData[2*i] = -1 / (f*f);
+			} else {
+				pfData[2*i] = 0.0f;
+			}
+		}
+
+		pfData[0] = 0.25f;
+
+		cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w);
+		delete[] ip;
+		delete[] w;
+
+		iFilterCacheSize = _iFFTRealDetectorCount;
+	}
+
 	for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
 	{
 		float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount;
 
-		// filt = 2*( 0:(order/2) )./order;
-		pfFilt[iDetectorIndex] = 2.0f * fRelIndex;
-		//pfFilt[iDetectorIndex] = 1.0f;
-
-		// w = 2*pi*(0:size(filt,2)-1)/order
-		pfW[iDetectorIndex] = 3.1415f * 2.0f * fRelIndex;
+		pfFilt[iDetectorIndex] = 2.0f * pfData[2*iDetectorIndex];
+		pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex;
 	}
 
 	switch(_eFilter)
@@ -866,5 +899,4 @@ void downloadDebugFilterReal(float * _pfHostSinogram, int _iProjectionCount,
 	free(pfHostFilter);
 }
 
-
 #endif
-- 
cgit v1.2.3