Add option to specify custom filter for FDK

6 files changed, 48 insertions, 7 deletions

diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu
index 5670873..35e3cd4 100644
--- a/cuda/3d/astra3d.cu
+++ b/cuda/3d/astra3d.cu
@@ -1311,7 +1311,7 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections,
                   const CVolumeGeometry3D* pVolGeom,
                   const CConeProjectionGeometry3D* pProjGeom,
                   bool bShortScan,
-                  int iGPUIndex, int iVoxelSuperSampling)
+                  int iGPUIndex, int iVoxelSuperSampling, const float* filter)
 {
 	SDimensions3D dims;
 
@@ -1369,7 +1369,7 @@ bool astraCudaFDK(float* pfVolume, const float* pfProjections,
 	// TODO: Offer interface for SrcZ, DetZ
 	ok &= FDK(D_volumeData, D_projData, fOriginSourceDistance,
 	          fOriginDetectorDistance, 0, 0, fDetUSize, fDetVSize,
-	          dims, pfAngles, bShortScan);
+	          dims, pfAngles, bShortScan, filter);
 
 	ok &= copyVolumeFromDevice(pfVolume, D_volumeData, dims, dims.iVolX);
 
diff --git a/cuda/3d/astra3d.h b/cuda/3d/astra3d.h
index 2137587..dde1347 100644
--- a/cuda/3d/astra3d.h
+++ b/cuda/3d/astra3d.h
@@ -314,7 +314,7 @@ _AstraExport bool astraCudaFDK(float* pfVolume, const float* pfProjections,
                   const CVolumeGeometry3D* pVolGeom,
                   const CConeProjectionGeometry3D* pProjGeom,
                   bool bShortScan,
-                  int iGPUIndex, int iVoxelSuperSampling);
+                  int iGPUIndex, int iVoxelSuperSampling, const float* filter);
 
 
 }
diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 0e13be1..4899ad1 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -394,7 +394,8 @@ bool FDK(cudaPitchedPtr D_volumeData,
          cudaPitchedPtr D_projData,
          float fSrcOrigin, float fDetOrigin,
          float fSrcZ, float fDetZ, float fDetUSize, float fDetVSize,
-         const SDimensions3D& dims, const float* angles, bool bShortScan)
+         const SDimensions3D& dims, const float* angles, bool bShortScan,
+	     const float* filter)
 {
 	bool ok;
 	// Generate filter
@@ -412,7 +413,14 @@ bool FDK(cudaPitchedPtr D_volumeData,
 	cufftComplex *pHostFilter = new cufftComplex[dims.iProjAngles * iHalfFFTSize];
 	memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iHalfFFTSize);
 
-	genFilter(FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize);
+	if (filter==NULL){
+		genFilter(FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize);
+	}else{
+		for(int i=0;i<dims.iProjAngles * iHalfFFTSize;i++){
+			pHostFilter[i].x = filter[i];
+			pHostFilter[i].y = 0;
+		}
+	}
 
 
 	allocateComplexOnDevice(dims.iProjAngles, iHalfFFTSize, &D_filter);
diff --git a/cuda/3d/fdk.h b/cuda/3d/fdk.h
index de7fe53..7a9d318 100644
--- a/cuda/3d/fdk.h
+++ b/cuda/3d/fdk.h
@@ -43,7 +43,8 @@ bool FDK(cudaPitchedPtr D_volumeData,
          cudaPitchedPtr D_projData,
          float fSrcOrigin, float fDetOrigin,
          float fSrcZ, float fDetZ, float fDetUSize, float fDetVSize,
-         const SDimensions3D& dims, const float* angles, bool bShortScan);
+         const SDimensions3D& dims, const float* angles, bool bShortScan,
+         const float* filter);
 
 
 }
diff --git a/include/astra/CudaFDKAlgorithm3D.h b/include/astra/CudaFDKAlgorithm3D.h
index 63f07fd..477bf34 100644
--- a/include/astra/CudaFDKAlgorithm3D.h
+++ b/include/astra/CudaFDKAlgorithm3D.h
@@ -151,6 +151,7 @@ protected:
 
 	int m_iGPUIndex;
 	int m_iVoxelSuperSampling;
+	int m_iFilterDataId;
 	bool m_bShortScan;
 
 	void initializeFromProjector();
diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp
index b5ce545..bc15a3d 100644
--- a/src/CudaFDKAlgorithm3D.cpp
+++ b/src/CudaFDKAlgorithm3D.cpp
@@ -36,8 +36,11 @@ $Id$
 #include "astra/Logging.h"
 
 #include "../cuda/3d/astra3d.h"
+#include "../cuda/2d/fft.h"
+#include "../cuda/3d/util3d.h"
 
 using namespace std;
+using namespace astraCUDA3d;
 
 namespace astra {
 
@@ -129,6 +132,28 @@ bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg)
 	CC.markOptionParsed("GPUIndex");
 	if (!_cfg.self.hasOption("GPUIndex"))
 		CC.markOptionParsed("GPUindex");
+	
+	// filter
+	if (_cfg.self.hasOption("FilterSinogramId")){
+		m_iFilterDataId = (int)_cfg.self.getOptionInt("FilterSinogramId");
+		const CFloat32ProjectionData2D * pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(m_iFilterDataId));
+		if (!pFilterData){
+			ASTRA_ERROR("Incorrect FilterSinogramId");
+			return false;
+		}
+		const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry();
+		const CProjectionGeometry2D* filtgeom = pFilterData->getGeometry();
+		int iPaddedDetCount = calcNextPowerOfTwo(2 * projgeom->getDetectorColCount());
+		int iHalfFFTSize = calcFFTFourSize(iPaddedDetCount);
+		if(filtgeom->getDetectorCount()!=iHalfFFTSize || filtgeom->getProjectionAngleCount()!=projgeom->getProjectionCount()){
+			ASTRA_ERROR("Filter size does not match required size (%i angles, %i detectors)",projgeom->getProjectionCount(),iHalfFFTSize);
+			return false;
+		}
+	}else
+	{
+		m_iFilterDataId = -1;
+	}
+	CC.markOptionParsed("FilterSinogramId");
 
 
 
@@ -196,10 +221,16 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations)
 
 
 	bool ok = true;
+	
+	const float *filter = NULL;
+	if(m_iFilterDataId!=-1){
+		const CFloat32ProjectionData2D * pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(m_iFilterDataId));
+		filter = pFilterData->getDataConst();
+	}
 
 	ok = astraCudaFDK(pReconMem->getData(), pSinoMem->getDataConst(),
 	                  &volgeom, conegeom,
-	                  m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling);
+	                  m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling, filter);
 
 	ASTRA_ASSERT(ok);