From 559d3e599b7306e2de64f2a584d72bc5c98b692b Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 4 Feb 2016 13:56:06 +0100
Subject: Refactor CUDA projector params into struct
---
cuda/2d/astra.cu | 1 -
1 file changed, 1 deletion(-)
(limited to 'cuda/2d')
diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu
index 4c69628..7317d69 100644
--- a/cuda/2d/astra.cu
+++ b/cuda/2d/astra.cu
@@ -341,7 +341,6 @@ bool AstraFBP::run()
dims3d.iProjAngles = pData->dims.iProjAngles;
dims3d.iProjU = pData->dims.iProjDets;
dims3d.iProjV = 1;
- dims3d.iRaysPerDetDim = dims3d.iRaysPerVoxelDim = 1;
astraCUDA3d::FDK_PreWeight(tmp, pData->fOriginSourceDistance,
pData->fOriginDetectorDistance, 0.0f, 0.0f,
--
cgit v1.2.3
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/astra.cu | 2 +-
cuda/2d/fft.cu | 45 +++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 46 insertions(+), 1 deletion(-)
(limited to 'cuda/2d')
diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu
index 7317d69..b56ddf9 100644
--- a/cuda/2d/astra.cu
+++ b/cuda/2d/astra.cu
@@ -343,7 +343,7 @@ bool AstraFBP::run()
dims3d.iProjV = 1;
astraCUDA3d::FDK_PreWeight(tmp, pData->fOriginSourceDistance,
- pData->fOriginDetectorDistance, 0.0f, 0.0f,
+ pData->fOriginDetectorDistance, 0.0f,
pData->dims.fDetScale, 1.0f, // TODO: Are these correct?
pData->bShortScan, dims3d, pData->angles);
}
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')
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