From 559d3e599b7306e2de64f2a584d72bc5c98b692b Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn 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 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 #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 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 #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