diff options
36 files changed, 1179 insertions, 1715 deletions
diff --git a/astra_vc09.vcproj b/astra_vc09.vcproj index f2cf62a..428bf8d 100644 --- a/astra_vc09.vcproj +++ b/astra_vc09.vcproj @@ -1077,6 +1077,10 @@ > </File> <File + RelativePath=".\include\astra\ParallelVecProjectionGeometry2D.h" + > + </File> + <File RelativePath=".\include\astra\ParallelVecProjectionGeometry3D.h" > </File> @@ -1121,6 +1125,10 @@ > </File> <File + RelativePath=".\src\GeometryUtil2D.cpp" + > + </File> + <File RelativePath=".\src\GeometryUtil3D.cpp" > </File> @@ -1133,6 +1141,10 @@ > </File> <File + RelativePath=".\src\ParallelVecProjectionGeometry2D.cpp" + > + </File> + <File RelativePath=".\src\ParallelVecProjectionGeometry3D.cpp" > </File> @@ -2185,6 +2197,10 @@ > </File> <File + RelativePath=".\cuda\2d\fbp.h" + > + </File> + <File RelativePath=".\cuda\2d\fft.h" > </File> @@ -2557,6 +2573,42 @@ </FileConfiguration> </File> <File + RelativePath=".\cuda\2d\fbp.cu" + > + <FileConfiguration + Name="Debug|Win32" + ExcludedFromBuild="true" + > + <Tool + Name="Cudart Build Rule" + /> + </FileConfiguration> + <FileConfiguration + Name="Debug|x64" + ExcludedFromBuild="true" + > + <Tool + Name="Cudart Build Rule" + /> + </FileConfiguration> + <FileConfiguration + Name="Release|Win32" + ExcludedFromBuild="true" + > + <Tool + Name="Cudart Build Rule" + /> + </FileConfiguration> + <FileConfiguration + Name="Release|x64" + ExcludedFromBuild="true" + > + <Tool + Name="Cudart Build Rule" + /> + </FileConfiguration> + </File> + <File RelativePath=".\cuda\2d\fft.cu" > <FileConfiguration diff --git a/astra_vc11.vcxproj b/astra_vc11.vcxproj index cd73076..d48796b 100644 --- a/astra_vc11.vcxproj +++ b/astra_vc11.vcxproj @@ -529,6 +529,7 @@ <ClCompile Include="src\Float32VolumeData3DMemory.cpp" /> <ClCompile Include="src\ForwardProjectionAlgorithm.cpp" /> <ClCompile Include="src\Fourier.cpp" /> + <ClCompile Include="src\GeometryUtil2D.cpp" /> <ClCompile Include="src\GeometryUtil3D.cpp" /> <ClCompile Include="src\Globals.cpp" /> <ClCompile Include="src\Logging.cpp" /> @@ -569,6 +570,7 @@ <ClInclude Include="cuda\2d\em.h" /> <ClInclude Include="cuda\2d\fan_bp.h" /> <ClInclude Include="cuda\2d\fan_fp.h" /> + <ClInclude Include="cuda\2d\fbp.h" /> <ClInclude Include="cuda\2d\fbp_filters.h" /> <ClInclude Include="cuda\2d\fft.h" /> <ClInclude Include="cuda\2d\par_bp.h" /> @@ -652,8 +654,8 @@ <ClInclude Include="include\astra\ParallelBeamStripKernelProjector2D.h" /> <ClInclude Include="include\astra\ParallelProjectionGeometry2D.h" /> <ClInclude Include="include\astra\ParallelProjectionGeometry3D.h" /> - <ClInclude Include="include\astra\ParallelVecProjectionGeometry3D.h" /> <ClInclude Include="include\astra\ParallelVecProjectionGeometry2D.h" /> + <ClInclude Include="include\astra\ParallelVecProjectionGeometry3D.h" /> <ClInclude Include="include\astra\PlatformDepSystemCode.h" /> <ClInclude Include="include\astra\PluginAlgorithm.h" /> <ClInclude Include="include\astra\ProjectionGeometry2D.h" /> @@ -727,6 +729,12 @@ <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">true</ExcludedFromBuild> <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</ExcludedFromBuild> </CudaCompile> + <CudaCompile Include="cuda\2d\fbp.cu"> + <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">true</ExcludedFromBuild> + <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">true</ExcludedFromBuild> + <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">true</ExcludedFromBuild> + <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</ExcludedFromBuild> + </CudaCompile> <CudaCompile Include="cuda\2d\fft.cu"> <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">true</ExcludedFromBuild> <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">true</ExcludedFromBuild> diff --git a/astra_vc11.vcxproj.filters b/astra_vc11.vcxproj.filters index 4143cd9..d6748f1 100644 --- a/astra_vc11.vcxproj.filters +++ b/astra_vc11.vcxproj.filters @@ -25,6 +25,9 @@ <CudaCompile Include="cuda\2d\fan_fp.cu"> <Filter>CUDA\cuda source</Filter> </CudaCompile> + <CudaCompile Include="cuda\2d\fbp.cu"> + <Filter>CUDA\cuda source</Filter> + </CudaCompile> <CudaCompile Include="cuda\2d\fft.cu"> <Filter>CUDA\cuda source</Filter> </CudaCompile> @@ -198,6 +201,9 @@ <ClCompile Include="src\FanFlatVecProjectionGeometry2D.cpp"> <Filter>Geometries\source</Filter> </ClCompile> + <ClCompile Include="src\GeometryUtil2D.cpp"> + <Filter>Geometries\source</Filter> + </ClCompile> <ClCompile Include="src\GeometryUtil3D.cpp"> <Filter>Geometries\source</Filter> </ClCompile> @@ -207,6 +213,9 @@ <ClCompile Include="src\ParallelProjectionGeometry3D.cpp"> <Filter>Geometries\source</Filter> </ClCompile> + <ClCompile Include="src\ParallelVecProjectionGeometry2D.cpp"> + <Filter>Geometries\source</Filter> + </ClCompile> <ClCompile Include="src\ParallelVecProjectionGeometry3D.cpp"> <Filter>Geometries\source</Filter> </ClCompile> @@ -321,10 +330,6 @@ <ClCompile Include="src\CudaSirtAlgorithm3D.cpp"> <Filter>CUDA\astra source</Filter> </ClCompile> - <ClCompile Include="src\GeometryUtil3D.cpp" /> - <ClCompile Include="src\ParallelVecProjectionGeometry2D.cpp"> - <Filter>Geometries\source</Filter> - </ClCompile> </ItemGroup> <ItemGroup> <ClInclude Include="include\astra\Algorithm.h"> @@ -474,6 +479,9 @@ <ClInclude Include="include\astra\ParallelProjectionGeometry3D.h"> <Filter>Geometries\headers</Filter> </ClInclude> + <ClInclude Include="include\astra\ParallelVecProjectionGeometry2D.h"> + <Filter>Geometries\headers</Filter> + </ClInclude> <ClInclude Include="include\astra\ParallelVecProjectionGeometry3D.h"> <Filter>Geometries\headers</Filter> </ClInclude> @@ -615,6 +623,9 @@ <ClInclude Include="cuda\2d\fbp_filters.h"> <Filter>CUDA\cuda headers</Filter> </ClInclude> + <ClInclude Include="cuda\2d\fbp.h"> + <Filter>CUDA\cuda headers</Filter> + </ClInclude> <ClInclude Include="cuda\2d\fft.h"> <Filter>CUDA\cuda headers</Filter> </ClInclude> @@ -675,11 +686,6 @@ <ClInclude Include="cuda\3d\util3d.h"> <Filter>CUDA\cuda headers</Filter> </ClInclude> - <ClInclude Include="include\astra\GeometryUtil2D.h" /> - <ClInclude Include="include\astra\GeometryUtil3D.h" /> - <ClInclude Include="include\astra\ParallelVecProjectionGeometry2D.h"> - <Filter>Geometries\headers</Filter> - </ClInclude> </ItemGroup> <ItemGroup> <None Include="include\astra\DataProjectorPolicies.inl"> diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index c72e20f..2af705f 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -139,6 +139,7 @@ BASE_OBJECTS=\ src/Float32VolumeData3DMemory.lo \ src/ForwardProjectionAlgorithm.lo \ src/Fourier.lo \ + src/GeometryUtil2D.lo \ src/GeometryUtil3D.lo \ src/Globals.lo \ src/Logging.lo \ @@ -197,6 +198,7 @@ CUDA_OBJECTS=\ cuda/2d/par_bp.lo \ cuda/2d/fan_fp.lo \ cuda/2d/fan_bp.lo \ + cuda/2d/fbp.lo \ cuda/2d/sirt.lo \ cuda/2d/sart.lo \ cuda/2d/cgls.lo \ diff --git a/build/msvc/gen.py b/build/msvc/gen.py index cc69a62..57f4539 100644 --- a/build/msvc/gen.py +++ b/build/msvc/gen.py @@ -156,6 +156,7 @@ P_astra["filters"]["CUDA\\cuda source"] = [ "cuda\\2d\\em.cu", "cuda\\2d\\fan_bp.cu", "cuda\\2d\\fan_fp.cu", +"cuda\\2d\\fbp.cu", "cuda\\2d\\fft.cu", "cuda\\2d\\par_bp.cu", "cuda\\2d\\par_fp.cu", @@ -225,9 +226,11 @@ P_astra["filters"]["Geometries\\source"] = [ "src\\ConeVecProjectionGeometry3D.cpp", "src\\FanFlatProjectionGeometry2D.cpp", "src\\FanFlatVecProjectionGeometry2D.cpp", +"src\\GeometryUtil2D.cpp", "src\\GeometryUtil3D.cpp", "src\\ParallelProjectionGeometry2D.cpp", "src\\ParallelProjectionGeometry3D.cpp", +"src\\ParallelVecProjectionGeometry2D.cpp", "src\\ParallelVecProjectionGeometry3D.cpp", "src\\ProjectionGeometry2D.cpp", "src\\ProjectionGeometry3D.cpp", @@ -285,6 +288,7 @@ P_astra["filters"]["CUDA\\cuda headers"] = [ "cuda\\2d\\fan_bp.h", "cuda\\2d\\fan_fp.h", "cuda\\2d\\fbp_filters.h", +"cuda\\2d\\fbp.h", "cuda\\2d\\fft.h", "cuda\\2d\\par_bp.h", "cuda\\2d\\par_fp.h", @@ -366,6 +370,7 @@ P_astra["filters"]["Geometries\\headers"] = [ "include\\astra\\GeometryUtil3D.h", "include\\astra\\ParallelProjectionGeometry2D.h", "include\\astra\\ParallelProjectionGeometry3D.h", +"include\\astra\\ParallelVecProjectionGeometry2D.h", "include\\astra\\ParallelVecProjectionGeometry3D.h", "include\\astra\\ProjectionGeometry2D.h", "include\\astra\\ProjectionGeometry3D.h", diff --git a/cuda/2d/algo.cu b/cuda/2d/algo.cu index dc74e51..bde4f42 100644 --- a/cuda/2d/algo.cu +++ b/cuda/2d/algo.cu @@ -35,13 +35,13 @@ $Id$ #include "fan_bp.h" #include "util.h" #include "arith.h" +#include "astra.h" namespace astraCUDA { ReconAlgo::ReconAlgo() { - angles = 0; - TOffsets = 0; + parProjs = 0; fanProjs = 0; shouldAbort = false; @@ -66,8 +66,7 @@ ReconAlgo::~ReconAlgo() void ReconAlgo::reset() { - delete[] angles; - delete[] TOffsets; + delete[] parProjs; delete[] fanProjs; if (freeGPUMemory) { @@ -77,8 +76,7 @@ void ReconAlgo::reset() cudaFree(D_volumeData); } - angles = 0; - TOffsets = 0; + parProjs = 0; fanProjs = 0; shouldAbort = false; @@ -124,46 +122,40 @@ bool ReconAlgo::enableSinogramMask() } -bool ReconAlgo::setGeometry(const SDimensions& _dims, const float* _angles) +bool ReconAlgo::setGeometry(const astra::CVolumeGeometry2D* pVolGeom, + const astra::CProjectionGeometry2D* pProjGeom) { - dims = _dims; + bool ok; - angles = new float[dims.iProjAngles]; + ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); - memcpy(angles, _angles, sizeof(angles[0]) * dims.iProjAngles); + if (!ok) + return false; + delete[] parProjs; + parProjs = 0; delete[] fanProjs; fanProjs = 0; - return true; -} - -bool ReconAlgo::setFanGeometry(const SDimensions& _dims, - const SFanProjection* _projs) -{ - dims = _dims; - fanProjs = new SFanProjection[dims.iProjAngles]; - - memcpy(fanProjs, _projs, sizeof(fanProjs[0]) * dims.iProjAngles); - - delete[] angles; - angles = 0; + fOutputScale = 1.0f; + ok = convertAstraGeometry(pVolGeom, pProjGeom, parProjs, fanProjs, fOutputScale); + if (!ok) + return false; return true; } - -bool ReconAlgo::setTOffsets(const float* _TOffsets) +bool ReconAlgo::setSuperSampling(int raysPerDet, int raysPerPixelDim) { - // TODO: determine if they're all zero? - TOffsets = new float[dims.iProjAngles]; - memcpy(TOffsets, _TOffsets, sizeof(angles[0]) * dims.iProjAngles); + if (raysPerDet <= 0 || raysPerPixelDim <= 0) + return false; + + dims.iRaysPerDet = raysPerDet; + dims.iRaysPerPixelDim = raysPerPixelDim; return true; } - - bool ReconAlgo::setVolumeMask(float* _D_maskData, unsigned int _maskPitch) { assert(useVolumeMask); @@ -324,14 +316,14 @@ bool ReconAlgo::callFP(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, float outputScale) { - if (angles) { + if (parProjs) { assert(!fanProjs); return FP(D_volumeData, volumePitch, D_projData, projPitch, - dims, angles, TOffsets, outputScale); + dims, parProjs, fOutputScale * outputScale); } else { assert(fanProjs); return FanFP(D_volumeData, volumePitch, D_projData, projPitch, - dims, fanProjs, outputScale); + dims, fanProjs, fOutputScale * outputScale); } } @@ -339,10 +331,10 @@ bool ReconAlgo::callBP(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, float outputScale) { - if (angles) { + if (parProjs) { assert(!fanProjs); return BP(D_volumeData, volumePitch, D_projData, projPitch, - dims, angles, TOffsets, outputScale); + dims, parProjs, outputScale); } else { assert(fanProjs); return FanBP(D_volumeData, volumePitch, D_projData, projPitch, diff --git a/cuda/2d/algo.h b/cuda/2d/algo.h index 99959c8..0b6a066 100644 --- a/cuda/2d/algo.h +++ b/cuda/2d/algo.h @@ -31,6 +31,17 @@ $Id$ #include "util.h" +namespace astra { + +class CParallelProjectionGeometry2D; +class CParallelVecProjectionGeometry2D; +class CFanFlatProjectionGeometry2D; +class CFanFlatVecProjectionGeometry2D; +class CVolumeGeometry2D; +class CProjectionGeometry2D; + +} + namespace astraCUDA { class _AstraExport ReconAlgo { @@ -40,11 +51,10 @@ public: bool setGPUIndex(int iGPUIndex); - bool setGeometry(const SDimensions& dims, const float* angles); - bool setFanGeometry(const SDimensions& dims, const SFanProjection* projs); + bool setGeometry(const astra::CVolumeGeometry2D* pVolGeom, + const astra::CProjectionGeometry2D* pProjGeom); - // setTOffsets should (optionally) be called after setGeometry - bool setTOffsets(const float* TOffsets); + bool setSuperSampling(int raysPerDet, int raysPerPixelDim); void signalAbort() { shouldAbort = true; } @@ -123,9 +133,9 @@ protected: SDimensions dims; - float* angles; - float* TOffsets; + SParProjection* parProjs; SFanProjection* fanProjs; + float fOutputScale; volatile bool shouldAbort; diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu index c1e6566..b39e0a9 100644 --- a/cuda/2d/astra.cu +++ b/cuda/2d/astra.cu @@ -42,8 +42,10 @@ $Id$ #include <fstream> #include <cuda.h> +#include "../../include/astra/GeometryUtil2D.h" #include "../../include/astra/VolumeGeometry2D.h" #include "../../include/astra/ParallelProjectionGeometry2D.h" +#include "../../include/astra/ParallelVecProjectionGeometry2D.h" #include "../../include/astra/FanFlatProjectionGeometry2D.h" #include "../../include/astra/FanFlatVecProjectionGeometry2D.h" @@ -64,515 +66,6 @@ enum CUDAProjectionType { }; -class AstraFBP_internal { -public: - SDimensions dims; - float* angles; - float* TOffsets; - astraCUDA::SFanProjection* fanProjections; - - float fOriginSourceDistance; - float fOriginDetectorDistance; - - float fPixelSize; - - bool bFanBeam; - bool bShortScan; - - bool initialized; - bool setStartReconstruction; - - float* D_sinoData; - unsigned int sinoPitch; - - float* D_volumeData; - unsigned int volumePitch; - - cufftComplex * m_pDevFilter; -}; - -AstraFBP::AstraFBP() -{ - pData = new AstraFBP_internal(); - - pData->angles = 0; - pData->fanProjections = 0; - pData->TOffsets = 0; - pData->D_sinoData = 0; - pData->D_volumeData = 0; - - pData->dims.iVolWidth = 0; - pData->dims.iProjAngles = 0; - pData->dims.fDetScale = 1.0f; - pData->dims.iRaysPerDet = 1; - pData->dims.iRaysPerPixelDim = 1; - - pData->initialized = false; - pData->setStartReconstruction = false; - - pData->m_pDevFilter = NULL; -} - -AstraFBP::~AstraFBP() -{ - delete[] pData->angles; - pData->angles = 0; - - delete[] pData->TOffsets; - pData->TOffsets = 0; - - delete[] pData->fanProjections; - pData->fanProjections = 0; - - cudaFree(pData->D_sinoData); - pData->D_sinoData = 0; - - cudaFree(pData->D_volumeData); - pData->D_volumeData = 0; - - if(pData->m_pDevFilter != NULL) - { - freeComplexOnDevice(pData->m_pDevFilter); - pData->m_pDevFilter = NULL; - } - - delete pData; - pData = 0; -} - -bool AstraFBP::setReconstructionGeometry(unsigned int iVolWidth, - unsigned int iVolHeight, - float fPixelSize) -{ - if (pData->initialized) - return false; - - pData->dims.iVolWidth = iVolWidth; - pData->dims.iVolHeight = iVolHeight; - - pData->fPixelSize = fPixelSize; - - return (iVolWidth > 0 && iVolHeight > 0 && fPixelSize > 0.0f); -} - -bool AstraFBP::setProjectionGeometry(unsigned int iProjAngles, - unsigned int iProjDets, - const float* pfAngles, - float fDetSize) -{ - if (pData->initialized) - return false; - - pData->dims.iProjAngles = iProjAngles; - pData->dims.iProjDets = iProjDets; - pData->dims.fDetScale = fDetSize / pData->fPixelSize; - - if (iProjAngles == 0 || iProjDets == 0 || pfAngles == 0) - return false; - - pData->angles = new float[iProjAngles]; - memcpy(pData->angles, pfAngles, iProjAngles * sizeof(pfAngles[0])); - - pData->bFanBeam = false; - - return true; -} - -bool AstraFBP::setFanGeometry(unsigned int iProjAngles, - unsigned int iProjDets, - const astraCUDA::SFanProjection *fanProjs, - const float* pfAngles, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetSize, - bool bShortScan) -{ - // Slightly abusing setProjectionGeometry for this... - if (!setProjectionGeometry(iProjAngles, iProjDets, pfAngles, fDetSize)) - return false; - - pData->fOriginSourceDistance = fOriginSourceDistance; - pData->fOriginDetectorDistance = fOriginDetectorDistance; - - pData->fanProjections = new astraCUDA::SFanProjection[iProjAngles]; - memcpy(pData->fanProjections, fanProjs, iProjAngles * sizeof(fanProjs[0])); - - pData->bFanBeam = true; - pData->bShortScan = bShortScan; - - return true; -} - - -bool AstraFBP::setPixelSuperSampling(unsigned int iPixelSuperSampling) -{ - if (pData->initialized) - return false; - - if (iPixelSuperSampling == 0) - return false; - - pData->dims.iRaysPerPixelDim = iPixelSuperSampling; - - return true; -} - - -bool AstraFBP::setTOffsets(const float* pfTOffsets) -{ - if (pData->initialized) - return false; - - if (pfTOffsets == 0) - return false; - - pData->TOffsets = new float[pData->dims.iProjAngles]; - memcpy(pData->TOffsets, pfTOffsets, pData->dims.iProjAngles * sizeof(pfTOffsets[0])); - - return true; -} - -bool AstraFBP::init(int iGPUIndex) -{ - if (pData->initialized) - { - return false; - } - - if (pData->dims.iProjAngles == 0 || pData->dims.iVolWidth == 0) - { - return false; - } - - if (iGPUIndex != -1) { - cudaSetDevice(iGPUIndex); - cudaError_t err = cudaGetLastError(); - - // Ignore errors caused by calling cudaSetDevice multiple times - if (err != cudaSuccess && err != cudaErrorSetOnActiveProcess) - { - return false; - } - } - - bool ok = allocateVolumeData(pData->D_volumeData, pData->volumePitch, pData->dims); - if (!ok) - { - return false; - } - - ok = allocateProjectionData(pData->D_sinoData, pData->sinoPitch, pData->dims); - if (!ok) - { - cudaFree(pData->D_volumeData); - pData->D_volumeData = 0; - return false; - } - - pData->initialized = true; - - return true; -} - -bool AstraFBP::setSinogram(const float* pfSinogram, - unsigned int iSinogramPitch) -{ - if (!pData->initialized) - return false; - if (!pfSinogram) - return false; - - bool ok = copySinogramToDevice(pfSinogram, iSinogramPitch, - pData->dims, - pData->D_sinoData, pData->sinoPitch); - if (!ok) - return false; - - // rescale sinogram to adjust for pixel size - processSino<opMul>(pData->D_sinoData, - 1.0f/(pData->fPixelSize*pData->fPixelSize), - pData->sinoPitch, pData->dims); - - pData->setStartReconstruction = false; - - return true; -} - -static int calcNextPowerOfTwo(int _iValue) -{ - int iOutput = 1; - - while(iOutput < _iValue) - { - iOutput *= 2; - } - - return iOutput; -} - -bool AstraFBP::run() -{ - if (!pData->initialized) - { - return false; - } - - zeroVolumeData(pData->D_volumeData, pData->volumePitch, pData->dims); - - bool ok = false; - - if (pData->bFanBeam) { - // Call FDK_PreWeight to handle fan beam geometry. We treat - // this as a cone beam setup of a single slice: - - // TODO: TOffsets affects this preweighting... - - // We create a fake cudaPitchedPtr - cudaPitchedPtr tmp; - tmp.ptr = pData->D_sinoData; - tmp.pitch = pData->sinoPitch * sizeof(float); - tmp.xsize = pData->dims.iProjDets; - tmp.ysize = pData->dims.iProjAngles; - // and a fake Dimensions3D - astraCUDA3d::SDimensions3D dims3d; - dims3d.iVolX = pData->dims.iVolWidth; - dims3d.iVolY = pData->dims.iVolHeight; - dims3d.iVolZ = 1; - 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, - pData->dims.fDetScale, 1.0f, // TODO: Are these correct? - pData->bShortScan, dims3d, pData->angles); - } - - if (pData->m_pDevFilter) { - - int iFFTRealDetCount = calcNextPowerOfTwo(2 * pData->dims.iProjDets); - int iFFTFourDetCount = calcFFTFourSize(iFFTRealDetCount); - - cufftComplex * pDevComplexSinogram = NULL; - - allocateComplexOnDevice(pData->dims.iProjAngles, iFFTFourDetCount, &pDevComplexSinogram); - - runCudaFFT(pData->dims.iProjAngles, pData->D_sinoData, pData->sinoPitch, pData->dims.iProjDets, iFFTRealDetCount, iFFTFourDetCount, pDevComplexSinogram); - - applyFilter(pData->dims.iProjAngles, iFFTFourDetCount, pDevComplexSinogram, pData->m_pDevFilter); - - runCudaIFFT(pData->dims.iProjAngles, pDevComplexSinogram, pData->D_sinoData, pData->sinoPitch, pData->dims.iProjDets, iFFTRealDetCount, iFFTFourDetCount); - - freeComplexOnDevice(pDevComplexSinogram); - - } - - float fOutputScale = (M_PI / 2.0f) / (float)pData->dims.iProjAngles; - - if (pData->bFanBeam) { - ok = FanBP_FBPWeighted(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->fanProjections, fOutputScale); - - } else { - ok = BP(pData->D_volumeData, pData->volumePitch, pData->D_sinoData, pData->sinoPitch, pData->dims, pData->angles, pData->TOffsets, fOutputScale); - } - if(!ok) - { - return false; - } - - return true; -} - -bool AstraFBP::getReconstruction(float* pfReconstruction, unsigned int iReconstructionPitch) const -{ - if (!pData->initialized) - return false; - - bool ok = copyVolumeFromDevice(pfReconstruction, iReconstructionPitch, - pData->dims, - pData->D_volumeData, pData->volumePitch); - if (!ok) - return false; - - return true; -} - -int AstraFBP::calcFourierFilterSize(int _iDetectorCount) -{ - int iFFTRealDetCount = calcNextPowerOfTwo(2 * _iDetectorCount); - int iFreqBinCount = calcFFTFourSize(iFFTRealDetCount); - - // CHECKME: Matlab makes this at least 64. Do we also need to? - return iFreqBinCount; -} - -bool AstraFBP::setFilter(E_FBPFILTER _eFilter, const float * _pfHostFilter /* = NULL */, int _iFilterWidth /* = 0 */, float _fD /* = 1.0f */, float _fFilterParameter /* = -1.0f */) -{ - if(pData->m_pDevFilter != 0) - { - freeComplexOnDevice(pData->m_pDevFilter); - pData->m_pDevFilter = 0; - } - - if (_eFilter == FILTER_NONE) - return true; // leave pData->m_pDevFilter set to 0 - - - int iFFTRealDetCount = calcNextPowerOfTwo(2 * pData->dims.iProjDets); - int iFreqBinCount = calcFFTFourSize(iFFTRealDetCount); - - cufftComplex * pHostFilter = new cufftComplex[pData->dims.iProjAngles * iFreqBinCount]; - memset(pHostFilter, 0, sizeof(cufftComplex) * pData->dims.iProjAngles * iFreqBinCount); - - allocateComplexOnDevice(pData->dims.iProjAngles, iFreqBinCount, &(pData->m_pDevFilter)); - - switch(_eFilter) - { - case FILTER_NONE: - // handled above - break; - case FILTER_RAMLAK: - case FILTER_SHEPPLOGAN: - case FILTER_COSINE: - case FILTER_HAMMING: - case FILTER_HANN: - case FILTER_TUKEY: - case FILTER_LANCZOS: - case FILTER_TRIANGULAR: - case FILTER_GAUSSIAN: - case FILTER_BARTLETTHANN: - case FILTER_BLACKMAN: - case FILTER_NUTTALL: - case FILTER_BLACKMANHARRIS: - case FILTER_BLACKMANNUTTALL: - case FILTER_FLATTOP: - case FILTER_PARZEN: - { - genFilter(_eFilter, _fD, pData->dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount, _fFilterParameter); - uploadComplexArrayToDevice(pData->dims.iProjAngles, iFreqBinCount, pHostFilter, pData->m_pDevFilter); - - break; - } - case FILTER_PROJECTION: - { - // make sure the offered filter has the correct size - assert(_iFilterWidth == iFreqBinCount); - - for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++) - { - float fValue = _pfHostFilter[iFreqBinIndex]; - - for(int iProjectionIndex = 0; iProjectionIndex < (int)pData->dims.iProjAngles; iProjectionIndex++) - { - pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].x = fValue; - pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].y = 0.0f; - } - } - uploadComplexArrayToDevice(pData->dims.iProjAngles, iFreqBinCount, pHostFilter, pData->m_pDevFilter); - break; - } - case FILTER_SINOGRAM: - { - // make sure the offered filter has the correct size - assert(_iFilterWidth == iFreqBinCount); - - for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++) - { - for(int iProjectionIndex = 0; iProjectionIndex < (int)pData->dims.iProjAngles; iProjectionIndex++) - { - float fValue = _pfHostFilter[iFreqBinIndex + iProjectionIndex * _iFilterWidth]; - - pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].x = fValue; - pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].y = 0.0f; - } - } - uploadComplexArrayToDevice(pData->dims.iProjAngles, iFreqBinCount, pHostFilter, pData->m_pDevFilter); - break; - } - case FILTER_RPROJECTION: - { - int iProjectionCount = pData->dims.iProjAngles; - int iRealFilterElementCount = iProjectionCount * iFFTRealDetCount; - float * pfHostRealFilter = new float[iRealFilterElementCount]; - memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount); - - int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount); - int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2; - int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth; - - int iFilterShiftSize = _iFilterWidth / 2; - - for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++) - { - int iFFTInFilterIndex = (iDetectorIndex + iFFTRealDetCount - iFilterShiftSize) % iFFTRealDetCount; - float fValue = _pfHostFilter[iDetectorIndex]; - - for(int iProjectionIndex = 0; iProjectionIndex < (int)pData->dims.iProjAngles; iProjectionIndex++) - { - pfHostRealFilter[iFFTInFilterIndex + iProjectionIndex * iFFTRealDetCount] = fValue; - } - } - - float* pfDevRealFilter = NULL; - cudaMalloc((void **)&pfDevRealFilter, sizeof(float) * iRealFilterElementCount); // TODO: check for errors - cudaMemcpy(pfDevRealFilter, pfHostRealFilter, sizeof(float) * iRealFilterElementCount, cudaMemcpyHostToDevice); - delete[] pfHostRealFilter; - - runCudaFFT(iProjectionCount, pfDevRealFilter, iFFTRealDetCount, iFFTRealDetCount, iFFTRealDetCount, iFreqBinCount, pData->m_pDevFilter); - - cudaFree(pfDevRealFilter); - - break; - } - case FILTER_RSINOGRAM: - { - int iProjectionCount = pData->dims.iProjAngles; - int iRealFilterElementCount = iProjectionCount * iFFTRealDetCount; - float* pfHostRealFilter = new float[iRealFilterElementCount]; - memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount); - - int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount); - int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2; - int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth; - - int iFilterShiftSize = _iFilterWidth / 2; - - for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++) - { - int iFFTInFilterIndex = (iDetectorIndex + iFFTRealDetCount - iFilterShiftSize) % iFFTRealDetCount; - - for(int iProjectionIndex = 0; iProjectionIndex < (int)pData->dims.iProjAngles; iProjectionIndex++) - { - float fValue = _pfHostFilter[iDetectorIndex + iProjectionIndex * _iFilterWidth]; - pfHostRealFilter[iFFTInFilterIndex + iProjectionIndex * iFFTRealDetCount] = fValue; - } - } - - float* pfDevRealFilter = NULL; - cudaMalloc((void **)&pfDevRealFilter, sizeof(float) * iRealFilterElementCount); // TODO: check for errors - cudaMemcpy(pfDevRealFilter, pfHostRealFilter, sizeof(float) * iRealFilterElementCount, cudaMemcpyHostToDevice); - delete[] pfHostRealFilter; - - runCudaFFT(iProjectionCount, pfDevRealFilter, iFFTRealDetCount, iFFTRealDetCount, iFFTRealDetCount, iFreqBinCount, pData->m_pDevFilter); - - cudaFree(pfDevRealFilter); - - break; - } - default: - { - ASTRA_ERROR("AstraFBP::setFilter: Unknown filter type requested"); - delete [] pHostFilter; - return false; - } - } - - delete [] pHostFilter; - - return true; -} - BPalgo::BPalgo() { @@ -617,18 +110,17 @@ float BPalgo::computeDiffNorm() bool astraCudaFP(const float* pfVolume, float* pfSinogram, unsigned int iVolWidth, unsigned int iVolHeight, unsigned int iProjAngles, unsigned int iProjDets, - const float *pfAngles, const float *pfOffsets, - float fDetSize, unsigned int iDetSuperSampling, + const SParProjection *pAngles, + unsigned int iDetSuperSampling, float fOutputScale, int iGPUIndex) { SDimensions dims; - if (iProjAngles == 0 || iProjDets == 0 || pfAngles == 0) + if (iProjAngles == 0 || iProjDets == 0 || pAngles == 0) return false; dims.iProjAngles = iProjAngles; dims.iProjDets = iProjDets; - dims.fDetScale = fDetSize; if (iDetSuperSampling == 0) return false; @@ -678,7 +170,7 @@ bool astraCudaFP(const float* pfVolume, float* pfSinogram, } zeroProjectionData(D_sinoData, sinoPitch, dims); - ok = FP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, pfAngles, pfOffsets, fOutputScale); + ok = FP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, pAngles, fOutputScale); if (!ok) { cudaFree(D_volumeData); cudaFree(D_sinoData); @@ -713,7 +205,6 @@ bool astraCudaFanFP(const float* pfVolume, float* pfSinogram, dims.iProjAngles = iProjAngles; dims.iProjDets = iProjDets; - dims.fDetScale = 1.0f; // TODO? if (iDetSuperSampling == 0) return false; @@ -789,118 +280,99 @@ bool astraCudaFanFP(const float* pfVolume, float* pfSinogram, } -bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, - const CParallelProjectionGeometry2D* pProjGeom, - float*& detectorOffsets, float*& projectionAngles, - float& detSize, float& outputScale) +// adjust pProjs to normalize volume geometry +template<typename ProjectionT> +static bool convertAstraGeometry_internal(const CVolumeGeometry2D* pVolGeom, + unsigned int iProjectionAngleCount, + ProjectionT*& pProjs, + float& fOutputScale) { - assert(pVolGeom); - assert(pProjGeom); - assert(pProjGeom->getProjectionAngles()); - + // TODO: Make EPS relative const float EPS = 0.00001f; - int nth = pProjGeom->getProjectionAngleCount(); - // Check if pixels are square if (abs(pVolGeom->getPixelLengthX() - pVolGeom->getPixelLengthY()) > EPS) return false; + float dx = -(pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2; + float dy = -(pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) / 2; - // Scale volume pixels to 1x1 - detSize = pProjGeom->getDetectorWidth() / pVolGeom->getPixelLengthX(); + float factor = 1.0f / pVolGeom->getPixelLengthX(); - // Copy angles - float *angles = new float[nth]; - for (int i = 0; i < nth; ++i) - angles[i] = pProjGeom->getProjectionAngles()[i]; - projectionAngles = angles; - - // Check if we need to translate - bool offCenter = false; - if (abs(pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) > EPS || - abs(pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) > EPS) - { - offCenter = true; + for (int i = 0; i < iProjectionAngleCount; ++i) { + // CHECKME: Order of scaling and translation + pProjs[i].translate(dx, dy); + pProjs[i].scale(factor); } + // CHECKME: Check factor + fOutputScale *= pVolGeom->getPixelLengthX() * pVolGeom->getPixelLengthY(); - // If there are existing detector offsets, or if we need to translate, - // we need to return offsets - if (offCenter) - { - float* offset = new float[nth]; + return true; +} - for (int i = 0; i < nth; ++i) - offset[i] = 0.0f; - if (offCenter) { - float dx = (pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2; - float dy = (pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) / 2; - // CHECKME: Is d in pixels or in units? +bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, + const CParallelProjectionGeometry2D* pProjGeom, + SParProjection*& pProjs, + float& fOutputScale) +{ + assert(pVolGeom); + assert(pProjGeom); + assert(pProjGeom->getProjectionAngles()); - for (int i = 0; i < nth; ++i) { - float d = dx * cos(angles[i]) + dy * sin(angles[i]); - offset[i] += d; - } - } + int nth = pProjGeom->getProjectionAngleCount(); - // CHECKME: Order of scaling and translation + pProjs = genParProjections(nth, + pProjGeom->getDetectorCount(), + pProjGeom->getDetectorWidth(), + pProjGeom->getProjectionAngles(), 0); - // Scale volume pixels to 1x1 - for (int i = 0; i < nth; ++i) { - //offset[i] /= pVolGeom->getPixelLengthX(); - //offset[i] *= detSize; - } + bool ok; + fOutputScale = 1.0f; + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); - detectorOffsets = offset; - } else { - detectorOffsets = 0; + if (!ok) { + delete[] pProjs; + pProjs = 0; } - outputScale = pVolGeom->getPixelLengthX(); - outputScale *= outputScale; - - return true; + return ok; } -static void convertAstraGeometry_internal(const CVolumeGeometry2D* pVolGeom, - unsigned int iProjectionAngleCount, - astraCUDA::SFanProjection*& pProjs, - float& outputScale) +bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, + const CParallelVecProjectionGeometry2D* pProjGeom, + SParProjection*& pProjs, + float& fOutputScale) { - // Translate - float dx = (pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2; - float dy = (pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) / 2; + assert(pVolGeom); + assert(pProjGeom); + assert(pProjGeom->getProjectionVectors()); - for (int i = 0; i < iProjectionAngleCount; ++i) { - pProjs[i].fSrcX -= dx; - pProjs[i].fSrcY -= dy; - pProjs[i].fDetSX -= dx; - pProjs[i].fDetSY -= dy; + int nth = pProjGeom->getProjectionAngleCount(); + + pProjs = new SParProjection[nth]; + + for (int i = 0; i < nth; ++i) { + pProjs[i] = pProjGeom->getProjectionVectors()[i]; } - // CHECKME: Order of scaling and translation + bool ok; + fOutputScale = 1.0f; - // Scale - float factor = 1.0f / pVolGeom->getPixelLengthX(); - for (int i = 0; i < iProjectionAngleCount; ++i) { - pProjs[i].fSrcX *= factor; - pProjs[i].fSrcY *= factor; - pProjs[i].fDetSX *= factor; - pProjs[i].fDetSY *= factor; - pProjs[i].fDetUX *= factor; - pProjs[i].fDetUY *= factor; + ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale); + if (!ok) { + delete[] pProjs; + pProjs = 0; } - // CHECKME: Check factor - outputScale = pVolGeom->getPixelLengthX(); -// outputScale *= outputScale; + return ok; } + bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, const CFanFlatProjectionGeometry2D* pProjGeom, astraCUDA::SFanProjection*& pProjs, @@ -910,6 +382,7 @@ bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, assert(pProjGeom); assert(pProjGeom->getProjectionAngles()); + // TODO: Make EPS relative const float EPS = 0.00001f; int nth = pProjGeom->getProjectionAngleCount(); @@ -928,23 +401,9 @@ bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, float fDetSize = pProjGeom->getDetectorWidth(); const float *pfAngles = pProjGeom->getProjectionAngles(); - pProjs = new SFanProjection[nth]; - - float fSrcX0 = 0.0f; - float fSrcY0 = -fOriginSourceDistance; - float fDetUX0 = fDetSize; - float fDetUY0 = 0.0f; - float fDetSX0 = pProjGeom->getDetectorCount() * fDetUX0 / -2.0f; - float fDetSY0 = fOriginDetectorDistance; - -#define ROTATE0(name,i,alpha) do { pProjs[i].f##name##X = f##name##X0 * cos(alpha) - f##name##Y0 * sin(alpha); pProjs[i].f##name##Y = f##name##X0 * sin(alpha) + f##name##Y0 * cos(alpha); } while(0) - for (int i = 0; i < nth; ++i) { - ROTATE0(Src, i, pfAngles[i]); - ROTATE0(DetS, i, pfAngles[i]); - ROTATE0(DetU, i, pfAngles[i]); - } - -#undef ROTATE0 + pProjs = genFanProjections(nth, pProjGeom->getDetectorCount(), + fOriginSourceDistance, fOriginDetectorDistance, + fDetSize, pfAngles); convertAstraGeometry_internal(pVolGeom, nth, pProjs, outputScale); @@ -961,6 +420,7 @@ bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, assert(pProjGeom); assert(pProjGeom->getProjectionVectors()); + // TODO: Make EPS relative const float EPS = 0.00001f; int nx = pVolGeom->getGridColCount(); @@ -983,6 +443,52 @@ bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, } +bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, + const CProjectionGeometry2D* pProjGeom, + astraCUDA::SParProjection*& pParProjs, + astraCUDA::SFanProjection*& pFanProjs, + float& outputScale) +{ + const CParallelProjectionGeometry2D* parProjGeom = dynamic_cast<const CParallelProjectionGeometry2D*>(pProjGeom); + const CParallelVecProjectionGeometry2D* parVecProjGeom = dynamic_cast<const CParallelVecProjectionGeometry2D*>(pProjGeom); + const CFanFlatProjectionGeometry2D* fanProjGeom = dynamic_cast<const CFanFlatProjectionGeometry2D*>(pProjGeom); + const CFanFlatVecProjectionGeometry2D* fanVecProjGeom = dynamic_cast<const CFanFlatVecProjectionGeometry2D*>(pProjGeom); + + bool ok = false; + + if (parProjGeom) { + ok = convertAstraGeometry(pVolGeom, parProjGeom, pParProjs, outputScale); + } else if (parVecProjGeom) { + ok = convertAstraGeometry(pVolGeom, parVecProjGeom, pParProjs, outputScale); + } else if (fanProjGeom) { + ok = convertAstraGeometry(pVolGeom, fanProjGeom, pFanProjs, outputScale); + } else if (fanVecProjGeom) { + ok = convertAstraGeometry(pVolGeom, fanVecProjGeom, pFanProjs, outputScale); + } else { + ok = false; + } + + return ok; +} + +bool convertAstraGeometry_dims(const CVolumeGeometry2D* pVolGeom, + const CProjectionGeometry2D* pProjGeom, + SDimensions& dims) +{ + dims.iVolWidth = pVolGeom->getGridColCount(); + dims.iVolHeight = pVolGeom->getGridRowCount(); + + dims.iProjAngles = pProjGeom->getProjectionAngleCount(); + dims.iProjDets = pProjGeom->getDetectorCount(); + + dims.iRaysPerDet = 1; + dims.iRaysPerPixelDim = 1; + + return true; +} + + + } diff --git a/cuda/2d/astra.h b/cuda/2d/astra.h index 617883b..8e91feb 100644 --- a/cuda/2d/astra.h +++ b/cuda/2d/astra.h @@ -29,7 +29,6 @@ $Id$ #ifndef _CUDA_ASTRA_H #define _CUDA_ASTRA_H -#include "fft.h" #include "fbp_filters.h" #include "dims.h" #include "algo.h" @@ -43,140 +42,12 @@ enum Cuda2DProjectionKernel { }; class CParallelProjectionGeometry2D; +class CParallelVecProjectionGeometry2D; class CFanFlatProjectionGeometry2D; class CFanFlatVecProjectionGeometry2D; class CVolumeGeometry2D; +class CProjectionGeometry2D; -class AstraFBP_internal; - -class _AstraExport AstraFBP { -public: - // Constructor - AstraFBP(); - - // Destructor - ~AstraFBP(); - - // Set the size of the reconstruction rectangle. - // Volume pixels are currently assumed to be 1x1 squares. - bool setReconstructionGeometry(unsigned int iVolWidth, - unsigned int iVolHeight, - float fPixelSize = 1.0f); - - // Set the projection angles and number of detector pixels per angle. - // pfAngles must be a float array of length iProjAngles. - // fDetSize indicates the size of a detector pixel compared to a - // volume pixel edge. - // - // pfAngles will only be read from during this call. - bool setProjectionGeometry(unsigned int iProjAngles, - unsigned int iProjDets, - const float *pfAngles, - float fDetSize = 1.0f); - // Set the projection angles and number of detector pixels per angle. - // pfAngles must be a float array of length iProjAngles. - // fDetSize indicates the size of a detector pixel compared to a - // volume pixel edge. - // - // pfAngles, fanProjs will only be read from during this call. - bool setFanGeometry(unsigned int iProjAngles, - unsigned int iProjDets, - const astraCUDA::SFanProjection *fanProjs, - const float *pfAngles, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetSize = 1.0f, - bool bShortScan = false); - - // Set linear supersampling factor for the BP. - // (The number of rays is the square of this) - // - // This may optionally be called before init(). - bool setPixelSuperSampling(unsigned int iPixelSuperSampling); - - // Set per-detector shifts. - // - // pfTOffsets will only be read from during this call. - bool setTOffsets(const float *pfTOffsets); - - // Returns the required size of a filter in the fourier domain - // when multiplying it with the fft of the projection data. - // Its value is equal to the smallest power of two larger than - // or equal to twice the number of detectors in the spatial domain. - // - // _iDetectorCount is the number of detectors in the spatial domain. - static int calcFourierFilterSize(int _iDetectorCount); - - // Sets the filter type. Some filter types require the user to supply an - // array containing the filter. - // The number of elements in a filter in the fourier domain should be equal - // to the value returned by calcFourierFilterSize(). - // The following types require a filter: - // - // - FILTER_PROJECTION: - // The filter size should be equal to the output of - // calcFourierFilterSize(). The filtered sinogram is - // multiplied with the supplied filter. - // - // - FILTER_SINOGRAM: - // Same as FILTER_PROJECTION, but now the filter should contain a row for - // every projection direction. - // - // - FILTER_RPROJECTION: - // The filter should now contain one kernel (= ifft of filter), with the - // peak in the center. The filter width - // can be any value. If odd, the peak is assumed to be in the center, if - // even, it is assumed to be at floor(filter-width/2). - // - // - FILTER_RSINOGRAM - // Same as FILTER_RPROJECTION, but now the supplied filter should contain a - // row for every projection direction. - // - // A large number of other filters (FILTER_RAMLAK, FILTER_SHEPPLOGAN, - // FILTER_COSINE, FILTER_HAMMING, and FILTER_HANN) - // have a D variable, which gives the cutoff point in the frequency domain. - // Setting this value to 1.0 will include the whole filter - bool setFilter(E_FBPFILTER _eFilter, - const float * _pfHostFilter = NULL, - int _iFilterWidth = 0, float _fD = 1.0f, float _fFilterParameter = -1.0f); - - // Initialize CUDA, allocate GPU buffers and - // precompute geometry-specific data. - // - // CUDA is set up to use GPU number iGPUIndex. - // - // This must be called after calling setReconstructionGeometry() and - // setProjectionGeometry(). - bool init(int iGPUIndex = 0); - - // Setup input sinogram for a slice. - // pfSinogram must be a float array of size iProjAngles*iSinogramPitch. - // NB: iSinogramPitch is measured in floats, not in bytes. - // - // This must be called after init(), and before iterate(). It may be - // called again after iterate()/getReconstruction() to start a new slice. - // - // pfSinogram will only be read from during this call. - bool setSinogram(const float* pfSinogram, unsigned int iSinogramPitch); - - // Runs an FBP reconstruction. - // This must be called after setSinogram(). - // - // run can be called before setFilter, but will then use the default Ram-Lak filter - bool run(); - - // Get the reconstructed slice. - // pfReconstruction must be a float array of size - // iVolHeight*iReconstructionPitch. - // NB: iReconstructionPitch is measured in floats, not in bytes. - // - // This may be called after run(). - bool getReconstruction(float* pfReconstruction, - unsigned int iReconstructionPitch) const; - -private: - AstraFBP_internal* pData; -}; class _AstraExport BPalgo : public astraCUDA::ReconAlgo { public: @@ -199,8 +70,8 @@ public: _AstraExport bool astraCudaFP(const float* pfVolume, float* pfSinogram, unsigned int iVolWidth, unsigned int iVolHeight, unsigned int iProjAngles, unsigned int iProjDets, - const float *pfAngles, const float *pfOffsets, - float fDetSize = 1.0f, unsigned int iDetSuperSampling = 1, + const SParProjection *pAngles, + unsigned int iDetSuperSampling = 1, float fOutputScale = 1.0f, int iGPUIndex = 0); _AstraExport bool astraCudaFanFP(const float* pfVolume, float* pfSinogram, @@ -213,8 +84,14 @@ _AstraExport bool astraCudaFanFP(const float* pfVolume, float* pfSinogram, _AstraExport bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, const CParallelProjectionGeometry2D* pProjGeom, - float*& pfDetectorOffsets, float*& pfProjectionAngles, - float& fDetSize, float& fOutputScale); + astraCUDA::SParProjection*& pProjs, + float& fOutputScale); + +_AstraExport bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, + const CParallelVecProjectionGeometry2D* pProjGeom, + astraCUDA::SParProjection*& pProjs, + float& fOutputScale); + _AstraExport bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, const CFanFlatProjectionGeometry2D* pProjGeom, @@ -226,6 +103,15 @@ _AstraExport bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, astraCUDA::SFanProjection*& pProjs, float& outputScale); +_AstraExport bool convertAstraGeometry_dims(const CVolumeGeometry2D* pVolGeom, + const CProjectionGeometry2D* pProjGeom, + astraCUDA::SDimensions& dims); + +_AstraExport bool convertAstraGeometry(const CVolumeGeometry2D* pVolGeom, + const CProjectionGeometry2D* pProjGeom, + astraCUDA::SParProjection*& pParProjs, + astraCUDA::SFanProjection*& pFanProjs, + float& outputScale); } #endif diff --git a/cuda/2d/cgls.cu b/cuda/2d/cgls.cu index f402914..3f06581 100644 --- a/cuda/2d/cgls.cu +++ b/cuda/2d/cgls.cu @@ -207,42 +207,6 @@ float CGLS::computeDiffNorm() return sqrt(s); } -bool doCGLS(float* D_volumeData, unsigned int volumePitch, - float* D_sinoData, unsigned int sinoPitch, - const SDimensions& dims, /*const SAugmentedData& augs,*/ - const float* angles, const float* TOffsets, unsigned int iterations) -{ - CGLS cgls; - bool ok = true; - - ok &= cgls.setGeometry(dims, angles); -#if 0 - if (D_maskData) - ok &= cgls.enableVolumeMask(); -#endif - if (TOffsets) - ok &= cgls.setTOffsets(TOffsets); - - if (!ok) - return false; - - ok = cgls.init(); - if (!ok) - return false; - -#if 0 - if (D_maskData) - ok &= cgls.setVolumeMask(D_maskData, maskPitch); -#endif - - ok &= cgls.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch); - if (!ok) - return false; - - ok = cgls.iterate(iterations); - - return ok; -} } diff --git a/cuda/2d/dims.h b/cuda/2d/dims.h index e870da5..f74d0e4 100644 --- a/cuda/2d/dims.h +++ b/cuda/2d/dims.h @@ -34,9 +34,11 @@ $Id$ namespace astraCUDA { +using astra::SParProjection; using astra::SFanProjection; + struct SDimensions { // Width, height of reconstruction volume unsigned int iVolWidth; @@ -48,9 +50,6 @@ struct SDimensions { // Number of detector pixels unsigned int iProjDets; - // size of detector compared to volume pixels - float fDetScale; - // in FP, number of rays to trace per detector pixel. // This should usually be set to 1. // If fDetScale > 1, this should be set to an integer of roughly diff --git a/cuda/2d/em.cu b/cuda/2d/em.cu index 8593b08..47ec548 100644 --- a/cuda/2d/em.cu +++ b/cuda/2d/em.cu @@ -170,34 +170,6 @@ float EM::computeDiffNorm() } -bool doEM(float* D_volumeData, unsigned int volumePitch, - float* D_sinoData, unsigned int sinoPitch, - const SDimensions& dims, const float* angles, - const float* TOffsets, unsigned int iterations) -{ - EM em; - bool ok = true; - - ok &= em.setGeometry(dims, angles); - if (TOffsets) - ok &= em.setTOffsets(TOffsets); - - if (!ok) - return false; - - ok = em.init(); - if (!ok) - return false; - - ok &= em.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch); - if (!ok) - return false; - - ok = em.iterate(iterations); - - return ok; -} - } #ifdef STANDALONE diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu new file mode 100644 index 0000000..2feac7d --- /dev/null +++ b/cuda/2d/fbp.cu @@ -0,0 +1,347 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp + 2014-2015, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://sf.net/projects/astra-toolbox + +This file is part of the ASTRA Toolbox. + + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +$Id$ +*/ + +#include "fbp.h" +#include "fft.h" +#include "par_bp.h" +#include "fan_bp.h" + +// For fan-beam preweighting +#include "../3d/fdk.h" + +#include "astra/Logging.h" + +#include <cuda.h> + +namespace astraCUDA { + + + +static int calcNextPowerOfTwo(int n) +{ + int x = 1; + while (x < n) + x *= 2; + + return x; +} + +// static +int FBP::calcFourierFilterSize(int _iDetectorCount) +{ + int iFFTRealDetCount = calcNextPowerOfTwo(2 * _iDetectorCount); + int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount); + + // CHECKME: Matlab makes this at least 64. Do we also need to? + return iFreqBinCount; +} + + + + +FBP::FBP() : ReconAlgo() +{ + D_filter = 0; + +} + +FBP::~FBP() +{ + reset(); +} + +void FBP::reset() +{ + if (D_filter) { + freeComplexOnDevice((cufftComplex *)D_filter); + D_filter = 0; + } +} + +bool FBP::init() +{ + return true; +} + +bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* = NULL */, int _iFilterWidth /* = 0 */, float _fD /* = 1.0f */, float _fFilterParameter /* = -1.0f */) +{ + if (D_filter) + { + freeComplexOnDevice((cufftComplex*)D_filter); + D_filter = 0; + } + + if (_eFilter == astra::FILTER_NONE) + return true; // leave D_filter set to 0 + + + int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets); + int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount); + + cufftComplex * pHostFilter = new cufftComplex[dims.iProjAngles * iFreqBinCount]; + memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iFreqBinCount); + + allocateComplexOnDevice(dims.iProjAngles, iFreqBinCount, (cufftComplex**)&D_filter); + + switch(_eFilter) + { + case astra::FILTER_NONE: + // handled above + break; + case astra::FILTER_RAMLAK: + case astra::FILTER_SHEPPLOGAN: + case astra::FILTER_COSINE: + case astra::FILTER_HAMMING: + case astra::FILTER_HANN: + case astra::FILTER_TUKEY: + case astra::FILTER_LANCZOS: + case astra::FILTER_TRIANGULAR: + case astra::FILTER_GAUSSIAN: + case astra::FILTER_BARTLETTHANN: + case astra::FILTER_BLACKMAN: + case astra::FILTER_NUTTALL: + case astra::FILTER_BLACKMANHARRIS: + case astra::FILTER_BLACKMANNUTTALL: + case astra::FILTER_FLATTOP: + case astra::FILTER_PARZEN: + { + genFilter(_eFilter, _fD, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount, _fFilterParameter); + uploadComplexArrayToDevice(dims.iProjAngles, iFreqBinCount, pHostFilter, (cufftComplex*)D_filter); + + break; + } + case astra::FILTER_PROJECTION: + { + // make sure the offered filter has the correct size + assert(_iFilterWidth == iFreqBinCount); + + for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++) + { + float fValue = _pfHostFilter[iFreqBinIndex]; + + for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++) + { + pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].x = fValue; + pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].y = 0.0f; + } + } + uploadComplexArrayToDevice(dims.iProjAngles, iFreqBinCount, pHostFilter, (cufftComplex*)D_filter); + break; + } + case astra::FILTER_SINOGRAM: + { + // make sure the offered filter has the correct size + assert(_iFilterWidth == iFreqBinCount); + + for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++) + { + for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++) + { + float fValue = _pfHostFilter[iFreqBinIndex + iProjectionIndex * _iFilterWidth]; + + pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].x = fValue; + pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].y = 0.0f; + } + } + uploadComplexArrayToDevice(dims.iProjAngles, iFreqBinCount, pHostFilter, (cufftComplex*)D_filter); + break; + } + case astra::FILTER_RPROJECTION: + { + int iProjectionCount = dims.iProjAngles; + int iRealFilterElementCount = iProjectionCount * iFFTRealDetCount; + float * pfHostRealFilter = new float[iRealFilterElementCount]; + memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount); + + int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount); + int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2; + int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth; + + int iFilterShiftSize = _iFilterWidth / 2; + + for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++) + { + int iFFTInFilterIndex = (iDetectorIndex + iFFTRealDetCount - iFilterShiftSize) % iFFTRealDetCount; + float fValue = _pfHostFilter[iDetectorIndex]; + + for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++) + { + pfHostRealFilter[iFFTInFilterIndex + iProjectionIndex * iFFTRealDetCount] = fValue; + } + } + + float* pfDevRealFilter = NULL; + cudaMalloc((void **)&pfDevRealFilter, sizeof(float) * iRealFilterElementCount); // TODO: check for errors + cudaMemcpy(pfDevRealFilter, pfHostRealFilter, sizeof(float) * iRealFilterElementCount, cudaMemcpyHostToDevice); + delete[] pfHostRealFilter; + + runCudaFFT(iProjectionCount, pfDevRealFilter, iFFTRealDetCount, iFFTRealDetCount, iFFTRealDetCount, iFreqBinCount, (cufftComplex*)D_filter); + + cudaFree(pfDevRealFilter); + + break; + } + case astra::FILTER_RSINOGRAM: + { + int iProjectionCount = dims.iProjAngles; + int iRealFilterElementCount = iProjectionCount * iFFTRealDetCount; + float* pfHostRealFilter = new float[iRealFilterElementCount]; + memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount); + + int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount); + int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2; + int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth; + + int iFilterShiftSize = _iFilterWidth / 2; + + for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++) + { + int iFFTInFilterIndex = (iDetectorIndex + iFFTRealDetCount - iFilterShiftSize) % iFFTRealDetCount; + + for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++) + { + float fValue = _pfHostFilter[iDetectorIndex + iProjectionIndex * _iFilterWidth]; + pfHostRealFilter[iFFTInFilterIndex + iProjectionIndex * iFFTRealDetCount] = fValue; + } + } + + float* pfDevRealFilter = NULL; + cudaMalloc((void **)&pfDevRealFilter, sizeof(float) * iRealFilterElementCount); // TODO: check for errors + cudaMemcpy(pfDevRealFilter, pfHostRealFilter, sizeof(float) * iRealFilterElementCount, cudaMemcpyHostToDevice); + delete[] pfHostRealFilter; + + runCudaFFT(iProjectionCount, pfDevRealFilter, iFFTRealDetCount, iFFTRealDetCount, iFFTRealDetCount, iFreqBinCount, (cufftComplex*)D_filter); + + cudaFree(pfDevRealFilter); + + break; + } + default: + { + ASTRA_ERROR("FBP::setFilter: Unknown filter type requested"); + delete [] pHostFilter; + return false; + } + } + + delete [] pHostFilter; + + return true; +} + +bool FBP::iterate(unsigned int iterations) +{ + zeroVolumeData(D_volumeData, volumePitch, dims); + + bool ok = false; + + if (fanProjs) { + // Call FDK_PreWeight to handle fan beam geometry. We treat + // this as a cone beam setup of a single slice: + + // TODO: TOffsets affects this preweighting... + + // TODO: We take the fan parameters from the last projection here + // without checking if they're the same in all projections + + float *pfAngles = new float[dims.iProjAngles]; + + float fOriginSource, fOriginDetector, fDetSize, fOffset; + for (unsigned int i = 0; i < dims.iProjAngles; ++i) { + bool ok = astra::getFanParameters(fanProjs[i], dims.iProjDets, + pfAngles[i], + fOriginSource, fOriginDetector, + fDetSize, fOffset); + if (!ok) { + ASTRA_ERROR("FBP_CUDA: Failed to extract circular fan beam parameters from fan beam geometry"); + return false; + } + } + + // We create a fake cudaPitchedPtr + cudaPitchedPtr tmp; + tmp.ptr = D_sinoData; + tmp.pitch = sinoPitch * sizeof(float); + tmp.xsize = dims.iProjDets; + tmp.ysize = dims.iProjAngles; + // and a fake Dimensions3D + astraCUDA3d::SDimensions3D dims3d; + dims3d.iVolX = dims.iVolWidth; + dims3d.iVolY = dims.iVolHeight; + dims3d.iVolZ = 1; + dims3d.iProjAngles = dims.iProjAngles; + dims3d.iProjU = dims.iProjDets; + dims3d.iProjV = 1; + dims3d.iRaysPerDetDim = dims3d.iRaysPerVoxelDim = 1; + + astraCUDA3d::FDK_PreWeight(tmp, fOriginSource, + fOriginDetector, 0.0f, 0.0f, + fDetSize, 1.0f, + m_bShortScan, dims3d, pfAngles); + } else { + // TODO: How should different detector pixel size in different + // projections be handled? + } + + if (D_filter) { + + int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets); + int iFFTFourDetCount = calcFFTFourierSize(iFFTRealDetCount); + + cufftComplex * pDevComplexSinogram = NULL; + + allocateComplexOnDevice(dims.iProjAngles, iFFTFourDetCount, &pDevComplexSinogram); + + runCudaFFT(dims.iProjAngles, D_sinoData, sinoPitch, dims.iProjDets, iFFTRealDetCount, iFFTFourDetCount, pDevComplexSinogram); + + applyFilter(dims.iProjAngles, iFFTFourDetCount, pDevComplexSinogram, (cufftComplex*)D_filter); + + runCudaIFFT(dims.iProjAngles, pDevComplexSinogram, D_sinoData, sinoPitch, dims.iProjDets, iFFTRealDetCount, iFFTFourDetCount); + + freeComplexOnDevice(pDevComplexSinogram); + + } + + float fOutputScale = (M_PI / 2.0f) / (float)dims.iProjAngles; + + if (fanProjs) { + ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fOutputScale); + + } else { + ok = BP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, parProjs, fOutputScale); + } + if(!ok) + { + return false; + } + + return true; +} + + +} diff --git a/cuda/2d/fbp.h b/cuda/2d/fbp.h new file mode 100644 index 0000000..8b4d64d --- /dev/null +++ b/cuda/2d/fbp.h @@ -0,0 +1,98 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp + 2014-2015, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://sf.net/projects/astra-toolbox + +This file is part of the ASTRA Toolbox. + + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +$Id$ +*/ + +#include "algo.h" +#include "fbp_filters.h" + +namespace astraCUDA { + +class _AstraExport FBP : public ReconAlgo { +public: + FBP(); + ~FBP(); + + virtual bool useSinogramMask() { return false; } + virtual bool useVolumeMask() { return false; } + + // Returns the required size of a filter in the fourier domain + // when multiplying it with the fft of the projection data. + // Its value is equal to the smallest power of two larger than + // or equal to twice the number of detectors in the spatial domain. + // + // _iDetectorCount is the number of detectors in the spatial domain. + static int calcFourierFilterSize(int _iDetectorCount); + + // Sets the filter type. Some filter types require the user to supply an + // array containing the filter. + // The number of elements in a filter in the fourier domain should be equal + // to the value returned by calcFourierFilterSize(). + // The following types require a filter: + // + // - FILTER_PROJECTION: + // The filter size should be equal to the output of + // calcFourierFilterSize(). The filtered sinogram is + // multiplied with the supplied filter. + // + // - FILTER_SINOGRAM: + // Same as FILTER_PROJECTION, but now the filter should contain a row for + // every projection direction. + // + // - FILTER_RPROJECTION: + // The filter should now contain one kernel (= ifft of filter), with the + // peak in the center. The filter width + // can be any value. If odd, the peak is assumed to be in the center, if + // even, it is assumed to be at floor(filter-width/2). + // + // - FILTER_RSINOGRAM + // Same as FILTER_RPROJECTION, but now the supplied filter should contain a + // row for every projection direction. + // + // A large number of other filters (FILTER_RAMLAK, FILTER_SHEPPLOGAN, + // FILTER_COSINE, FILTER_HAMMING, and FILTER_HANN) + // have a D variable, which gives the cutoff point in the frequency domain. + // Setting this value to 1.0 will include the whole filter + bool setFilter(astra::E_FBPFILTER _eFilter, + const float * _pfHostFilter = NULL, + int _iFilterWidth = 0, float _fD = 1.0f, float _fFilterParameter = -1.0f); + + bool setShortScan(bool ss) { m_bShortScan = ss; return true; } + + virtual bool init(); + + virtual bool iterate(unsigned int iterations); + + virtual float computeDiffNorm() { return 0.0f; } // TODO + +protected: + void reset(); + + void* D_filter; // cufftComplex* + bool m_bShortScan; +}; + +} diff --git a/cuda/2d/fbp_filters.h b/cuda/2d/fbp_filters.h index b206a5d..71c7572 100644 --- a/cuda/2d/fbp_filters.h +++ b/cuda/2d/fbp_filters.h @@ -29,6 +29,8 @@ $Id$ #ifndef FBP_FILTERS_H #define FBP_FILTERS_H +namespace astra { + enum E_FBPFILTER { FILTER_NONE, //< no filter (regular BP) @@ -55,4 +57,6 @@ enum E_FBPFILTER FILTER_RSINOGRAM, //< sinogram filter in real space }; +} + #endif /* FBP_FILTERS_H */ diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index 2d259a9..8a7cc10 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -64,6 +64,8 @@ using namespace astra; } } while (0) +namespace astraCUDA { + __global__ static void applyFilter_kernel(int _iProjectionCount, int _iFreqBinCount, cufftComplex * _pSinogram, @@ -276,11 +278,11 @@ bool runCudaIFFT(int _iProjectionCount, const cufftComplex* _pDevSourceComplex, // Because the input is real, the Fourier transform is symmetric. // CUFFT only outputs the first half (ignoring the redundant second half), // and expects the same as input for the IFFT. -int calcFFTFourSize(int _iFFTRealSize) +int calcFFTFourierSize(int _iFFTRealSize) { - int iFFTFourSize = _iFFTRealSize / 2 + 1; + int iFFTFourierSize = _iFFTRealSize / 2 + 1; - return iFFTFourSize; + return iFFTFourierSize; } void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, @@ -695,6 +697,10 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, delete[] pfW; } + +} + + #ifdef STANDALONE __global__ static void doubleFourierOutput_kernel(int _iProjectionCount, diff --git a/cuda/2d/fft.h b/cuda/2d/fft.h index e985fc6..b29f17a 100644 --- a/cuda/2d/fft.h +++ b/cuda/2d/fft.h @@ -34,6 +34,8 @@ $Id$ #include "fbp_filters.h" +namespace astraCUDA { + bool allocateComplexOnDevice(int _iProjectionCount, int _iDetectorCount, cufftComplex ** _ppDevComplex); @@ -57,13 +59,15 @@ bool runCudaIFFT(int _iProjectionCount, const cufftComplex* _pDevSourceComplex, void applyFilter(int _iProjectionCount, int _iFreqBinCount, cufftComplex * _pSinogram, cufftComplex * _pFilter); -int calcFFTFourSize(int _iFFTRealSize); +int calcFFTFourierSize(int _iFFTRealSize); -void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, +void genFilter(astra::E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount, float _fParameter = -1.0f); void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount); +} + #endif /* FFT_H */ diff --git a/cuda/2d/par_bp.cu b/cuda/2d/par_bp.cu index d9f7325..cf0a684 100644 --- a/cuda/2d/par_bp.cu +++ b/cuda/2d/par_bp.cu @@ -53,8 +53,8 @@ const unsigned int g_blockSlices = 16; const unsigned int g_MaxAngles = 2560; -__constant__ float gC_angle_sin[g_MaxAngles]; -__constant__ float gC_angle_cos[g_MaxAngles]; +__constant__ float gC_angle_scaled_sin[g_MaxAngles]; +__constant__ float gC_angle_scaled_cos[g_MaxAngles]; __constant__ float gC_angle_offset[g_MaxAngles]; static bool bindProjDataTexture(float* data, unsigned int pitch, unsigned int width, unsigned int height, cudaTextureAddressMode mode = cudaAddressModeBorder) @@ -73,7 +73,7 @@ static bool bindProjDataTexture(float* data, unsigned int pitch, unsigned int wi return true; } -__global__ void devBP(float* D_volData, unsigned int volPitch, unsigned int startAngle, bool offsets, const SDimensions dims, float fOutputScale) +__global__ void devBP(float* D_volData, unsigned int volPitch, unsigned int startAngle, const SDimensions dims, float fOutputScale) { const int relX = threadIdx.x; const int relY = threadIdx.y; @@ -87,47 +87,30 @@ __global__ void devBP(float* D_volData, unsigned int volPitch, unsigned int star if (X >= dims.iVolWidth || Y >= dims.iVolHeight) return; - const float fX = ( X - 0.5f*dims.iVolWidth + 0.5f ) / dims.fDetScale; - const float fY = ( Y - 0.5f*dims.iVolHeight + 0.5f ) / dims.fDetScale; + const float fX = ( X - 0.5f*dims.iVolWidth + 0.5f ); + const float fY = ( Y - 0.5f*dims.iVolHeight + 0.5f ); float* volData = (float*)D_volData; float fVal = 0.0f; float fA = startAngle + 0.5f; - const float fT_base = 0.5f*dims.iProjDets - 0.5f + 0.5f; - if (offsets) { - - for (int angle = startAngle; angle < endAngle; ++angle) - { - const float cos_theta = gC_angle_cos[angle]; - const float sin_theta = gC_angle_sin[angle]; - const float TOffset = gC_angle_offset[angle]; - - const float fT = fT_base + fX * cos_theta - fY * sin_theta + TOffset; - fVal += tex2D(gT_projTexture, fT, fA); - fA += 1.0f; - } - - } else { - - for (int angle = startAngle; angle < endAngle; ++angle) - { - const float cos_theta = gC_angle_cos[angle]; - const float sin_theta = gC_angle_sin[angle]; - - const float fT = fT_base + fX * cos_theta - fY * sin_theta; - fVal += tex2D(gT_projTexture, fT, fA); - fA += 1.0f; - } + for (int angle = startAngle; angle < endAngle; ++angle) + { + const float scaled_cos_theta = gC_angle_scaled_cos[angle]; + const float scaled_sin_theta = gC_angle_scaled_sin[angle]; + const float TOffset = gC_angle_offset[angle]; + const float fT = fX * scaled_cos_theta - fY * scaled_sin_theta + TOffset; + fVal += tex2D(gT_projTexture, fT, fA); + fA += 1.0f; } volData[Y*volPitch+X] += fVal * fOutputScale; } // supersampling version -__global__ void devBP_SS(float* D_volData, unsigned int volPitch, unsigned int startAngle, bool offsets, const SDimensions dims, float fOutputScale) +__global__ void devBP_SS(float* D_volData, unsigned int volPitch, unsigned int startAngle, const SDimensions dims, float fOutputScale) { const int relX = threadIdx.x; const int relY = threadIdx.y; @@ -141,61 +124,35 @@ __global__ void devBP_SS(float* D_volData, unsigned int volPitch, unsigned int s if (X >= dims.iVolWidth || Y >= dims.iVolHeight) return; - const float fX = ( X - 0.5f*dims.iVolWidth + 0.5f - 0.5f + 0.5f/dims.iRaysPerPixelDim) / dims.fDetScale; - const float fY = ( Y - 0.5f*dims.iVolHeight + 0.5f - 0.5f + 0.5f/dims.iRaysPerPixelDim) / dims.fDetScale; + const float fX = ( X - 0.5f*dims.iVolWidth + 0.5f - 0.5f + 0.5f/dims.iRaysPerPixelDim); + const float fY = ( Y - 0.5f*dims.iVolHeight + 0.5f - 0.5f + 0.5f/dims.iRaysPerPixelDim); - const float fSubStep = 1.0f/(dims.iRaysPerPixelDim * dims.fDetScale); + const float fSubStep = 1.0f/(dims.iRaysPerPixelDim); // * dims.fDetScale); float* volData = (float*)D_volData; float fVal = 0.0f; float fA = startAngle + 0.5f; - const float fT_base = 0.5f*dims.iProjDets - 0.5f + 0.5f; fOutputScale /= (dims.iRaysPerPixelDim * dims.iRaysPerPixelDim); - if (offsets) { - - for (int angle = startAngle; angle < endAngle; ++angle) - { - const float cos_theta = gC_angle_cos[angle]; - const float sin_theta = gC_angle_sin[angle]; - const float TOffset = gC_angle_offset[angle]; - - float fT = fT_base + fX * cos_theta - fY * sin_theta + TOffset; - - for (int iSubX = 0; iSubX < dims.iRaysPerPixelDim; ++iSubX) { - float fTy = fT; - fT += fSubStep * cos_theta; - for (int iSubY = 0; iSubY < dims.iRaysPerPixelDim; ++iSubY) { - fVal += tex2D(gT_projTexture, fTy, fA); - fTy -= fSubStep * sin_theta; - } - } - fA += 1.0f; - } - - } else { + for (int angle = startAngle; angle < endAngle; ++angle) + { + const float cos_theta = gC_angle_scaled_cos[angle]; + const float sin_theta = gC_angle_scaled_sin[angle]; + const float TOffset = gC_angle_offset[angle]; - for (int angle = startAngle; angle < endAngle; ++angle) - { - const float cos_theta = gC_angle_cos[angle]; - const float sin_theta = gC_angle_sin[angle]; + float fT = fX * cos_theta - fY * sin_theta + TOffset; - float fT = fT_base + fX * cos_theta - fY * sin_theta; - - for (int iSubX = 0; iSubX < dims.iRaysPerPixelDim; ++iSubX) { - float fTy = fT; - fT += fSubStep * cos_theta; - for (int iSubY = 0; iSubY < dims.iRaysPerPixelDim; ++iSubY) { - fVal += tex2D(gT_projTexture, fTy, fA); - fTy -= fSubStep * sin_theta; - } + for (int iSubX = 0; iSubX < dims.iRaysPerPixelDim; ++iSubX) { + float fTy = fT; + fT += fSubStep * cos_theta; + for (int iSubY = 0; iSubY < dims.iRaysPerPixelDim; ++iSubY) { + fVal += tex2D(gT_projTexture, fTy, fA); + fTy -= fSubStep * sin_theta; } - fA += 1.0f; - } - + fA += 1.0f; } volData[Y*volPitch+X] += fVal * fOutputScale; @@ -212,12 +169,10 @@ __global__ void devBP_SART(float* D_volData, unsigned int volPitch, float offset if (X >= dims.iVolWidth || Y >= dims.iVolHeight) return; - const float fX = ( X - 0.5f*dims.iVolWidth + 0.5f ) / dims.fDetScale; - const float fY = ( Y - 0.5f*dims.iVolHeight + 0.5f ) / dims.fDetScale; - - const float fT_base = 0.5f*dims.iProjDets - 0.5f + 0.5f; + const float fX = ( X - 0.5f*dims.iVolWidth + 0.5f ); + const float fY = ( Y - 0.5f*dims.iVolHeight + 0.5f ); - const float fT = fT_base + fX * angle_cos - fY * angle_sin + offset; + const float fT = fX * angle_cos - fY * angle_sin + offset; const float fVal = tex2D(gT_projTexture, fT, 0.5f); D_volData[Y*volPitch+X] += fVal * fOutputScale; @@ -226,32 +181,35 @@ __global__ void devBP_SART(float* D_volData, unsigned int volPitch, float offset bool BP_internal(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, - const SDimensions& dims, const float* angles, const float* TOffsets, float fOutputScale) + const SDimensions& dims, const SParProjection* angles, + float fOutputScale) { - // TODO: process angles block by block assert(dims.iProjAngles <= g_MaxAngles); - float* angle_sin = new float[dims.iProjAngles]; - float* angle_cos = new float[dims.iProjAngles]; + float* angle_scaled_sin = new float[dims.iProjAngles]; + float* angle_scaled_cos = new float[dims.iProjAngles]; + float* angle_offset = new float[dims.iProjAngles]; bindProjDataTexture(D_projData, projPitch, dims.iProjDets, dims.iProjAngles); for (unsigned int i = 0; i < dims.iProjAngles; ++i) { - angle_sin[i] = sinf(angles[i]); - angle_cos[i] = cosf(angles[i]); + double d = angles[i].fDetUX * angles[i].fRayY - angles[i].fDetUY * angles[i].fRayX; + angle_scaled_cos[i] = angles[i].fRayY / d; + angle_scaled_sin[i] = -angles[i].fRayX / d; // TODO: Check signs + angle_offset[i] = (angles[i].fDetSY * angles[i].fRayX - angles[i].fDetSX * angles[i].fRayY) / d; } - cudaError_t e1 = cudaMemcpyToSymbol(gC_angle_sin, angle_sin, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - cudaError_t e2 = cudaMemcpyToSymbol(gC_angle_cos, angle_cos, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); + + cudaError_t e1 = cudaMemcpyToSymbol(gC_angle_scaled_sin, angle_scaled_sin, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); + cudaError_t e2 = cudaMemcpyToSymbol(gC_angle_scaled_cos, angle_scaled_cos, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); + cudaError_t e3 = cudaMemcpyToSymbol(gC_angle_offset, angle_offset, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); assert(e1 == cudaSuccess); assert(e2 == cudaSuccess); + assert(e3 == cudaSuccess); - if (TOffsets) { - cudaError_t e3 = cudaMemcpyToSymbol(gC_angle_offset, TOffsets, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - assert(e3 == cudaSuccess); - } - delete[] angle_sin; - delete[] angle_cos; + delete[] angle_scaled_sin; + delete[] angle_scaled_cos; + delete[] angle_offset; dim3 dimBlock(g_blockSlices, g_blockSliceSize); dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices, @@ -263,9 +221,9 @@ bool BP_internal(float* D_volumeData, unsigned int volumePitch, for (unsigned int i = 0; i < dims.iProjAngles; i += g_anglesPerBlock) { if (dims.iRaysPerPixelDim > 1) - devBP_SS<<<dimGrid, dimBlock, 0, stream>>>(D_volumeData, volumePitch, i, (TOffsets != 0), dims, fOutputScale); + devBP_SS<<<dimGrid, dimBlock, 0, stream>>>(D_volumeData, volumePitch, i, dims, fOutputScale); else - devBP<<<dimGrid, dimBlock, 0, stream>>>(D_volumeData, volumePitch, i, (TOffsets != 0), dims, fOutputScale); + devBP<<<dimGrid, dimBlock, 0, stream>>>(D_volumeData, volumePitch, i, dims, fOutputScale); } cudaThreadSynchronize(); @@ -278,7 +236,7 @@ bool BP_internal(float* D_volumeData, unsigned int volumePitch, bool BP(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, - const SDimensions& dims, const float* angles, const float* TOffsets, float fOutputScale) + const SDimensions& dims, const SParProjection* angles, float fOutputScale) { for (unsigned int iAngle = 0; iAngle < dims.iProjAngles; iAngle += g_MaxAngles) { SDimensions subdims = dims; @@ -290,9 +248,7 @@ bool BP(float* D_volumeData, unsigned int volumePitch, bool ret; ret = BP_internal(D_volumeData, volumePitch, D_projData + iAngle * projPitch, projPitch, - subdims, angles + iAngle, - TOffsets ? TOffsets + iAngle : 0, - fOutputScale); + subdims, angles + iAngle, fOutputScale); if (!ret) return false; } @@ -303,25 +259,23 @@ bool BP(float* D_volumeData, unsigned int volumePitch, bool BP_SART(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, unsigned int angle, const SDimensions& dims, - const float* angles, const float* TOffsets, float fOutputScale) + const SParProjection* angles, float fOutputScale) { // Only one angle. // We need to Clamp to the border pixels instead of to zero, because // SART weights with ray length. bindProjDataTexture(D_projData, projPitch, dims.iProjDets, 1, cudaAddressModeClamp); - float angle_sin = sinf(angles[angle]); - float angle_cos = cosf(angles[angle]); - - float offset = 0.0f; - if (TOffsets) - offset = TOffsets[angle]; + double d = angles[angle].fDetUX * angles[angle].fRayY - angles[angle].fDetUY * angles[angle].fRayX; + float angle_scaled_cos = angles[angle].fRayY / d; + float angle_scaled_sin = -angles[angle].fRayX / d; // TODO: Check signs + float angle_offset = (angles[angle].fDetSY * angles[angle].fRayX - angles[angle].fDetSX * angles[angle].fRayY) / d; dim3 dimBlock(g_blockSlices, g_blockSliceSize); dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices, (dims.iVolHeight+g_blockSliceSize-1)/g_blockSliceSize); - devBP_SART<<<dimGrid, dimBlock>>>(D_volumeData, volumePitch, offset, angle_sin, angle_cos, dims, fOutputScale); + devBP_SART<<<dimGrid, dimBlock>>>(D_volumeData, volumePitch, angle_offset, angle_scaled_sin, angle_scaled_cos, dims, fOutputScale); cudaThreadSynchronize(); cudaTextForceKernelsCompletion(); diff --git a/cuda/2d/par_bp.h b/cuda/2d/par_bp.h index 64bcd34..102cb2b 100644 --- a/cuda/2d/par_bp.h +++ b/cuda/2d/par_bp.h @@ -35,13 +35,13 @@ namespace astraCUDA { _AstraExport bool BP(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, - const SDimensions& dims, const float* angles, - const float* TOffsets, float fOutputScale); + const SDimensions& dims, const SParProjection* angles, + float fOutputScale); _AstraExport bool BP_SART(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, unsigned int angle, const SDimensions& dims, - const float* angles, const float* TOffsets, float fOutputScale); + const SParProjection* angles, float fOutputScale); } diff --git a/cuda/2d/par_fp.cu b/cuda/2d/par_fp.cu index bb8b909..3c010b3 100644 --- a/cuda/2d/par_fp.cu +++ b/cuda/2d/par_fp.cu @@ -30,6 +30,7 @@ $Id$ #include <cassert> #include <iostream> #include <list> +#include <cmath> #include "util.h" #include "arith.h" @@ -51,6 +52,7 @@ namespace astraCUDA { static const unsigned g_MaxAngles = 2560; __constant__ float gC_angle[g_MaxAngles]; __constant__ float gC_angle_offset[g_MaxAngles]; +__constant__ float gC_angle_detsize[g_MaxAngles]; // optimization parameters @@ -63,10 +65,6 @@ static const unsigned int g_blockSlices = 64; #define iPREC_FACTOR 16 -// if necessary, a buffer of zeroes of size g_MaxAngles -static float* g_pfZeroes = 0; - - static bool bindVolumeDataTexture(float* data, cudaArray*& dataArray, unsigned int pitch, unsigned int width, unsigned int height) { cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>(); @@ -89,9 +87,10 @@ static bool bindVolumeDataTexture(float* data, cudaArray*& dataArray, unsigned i return true; } + // projection for angles that are roughly horizontal -// theta between 45 and 135 degrees -__global__ void FPhorizontal(float* D_projData, unsigned int projPitch, unsigned int startSlice, unsigned int startAngle, unsigned int endAngle, int regionOffset, const SDimensions dims, float outputScale) +// (detector roughly vertical) +__global__ void FPhorizontal_simple(float* D_projData, unsigned int projPitch, unsigned int startSlice, unsigned int startAngle, unsigned int endAngle, const SDimensions dims, float outputScale) { const int relDet = threadIdx.x; const int relAngle = threadIdx.y; @@ -106,33 +105,7 @@ __global__ void FPhorizontal(float* D_projData, unsigned int projPitch, unsigned const float sin_theta = __sinf(theta); // compute start detector for this block/angle: - // (The same for all threadIdx.x) - // ------------------------------------- - - const int midSlice = startSlice + g_blockSlices / 2; - - // ASSUMPTION: fDetScale >= 1.0f - // problem: detector regions get skipped because slice blocks aren't large - // enough - const unsigned int g_blockSliceSize = g_detBlockSize; - - // project (midSlice,midRegion) on this thread's detector - - const float fBase = 0.5f*dims.iProjDets - 0.5f + - ( - (midSlice - 0.5f*dims.iVolWidth + 0.5f) * cos_theta - - (g_blockSliceSize/2 - 0.5f*dims.iVolHeight + 0.5f) * sin_theta - + gC_angle_offset[angle] - ) / dims.fDetScale; - int iBase = (int)floorf(fBase * fPREC_FACTOR); - int iInc = (int)floorf(g_blockSliceSize * sin_theta / dims.fDetScale * -fPREC_FACTOR); - - // ASSUMPTION: 16 > regionOffset / fDetScale - const int detRegion = (iBase + (blockIdx.y - regionOffset)*iInc + 16*iPREC_FACTOR*g_detBlockSize) / (iPREC_FACTOR * g_detBlockSize) - 16; - const int detPrevRegion = (iBase + (blockIdx.y - regionOffset - 1)*iInc + 16*iPREC_FACTOR*g_detBlockSize) / (iPREC_FACTOR * g_detBlockSize) - 16; - - if (blockIdx.y > 0 && detRegion == detPrevRegion) - return; + const int detRegion = blockIdx.y; const int detector = detRegion * g_detBlockSize + relDet; @@ -142,7 +115,7 @@ __global__ void FPhorizontal(float* D_projData, unsigned int projPitch, unsigned if (detector < 0 || detector >= dims.iProjDets) return; - const float fDetStep = -dims.fDetScale / sin_theta; + const float fDetStep = -gC_angle_detsize[angle] / sin_theta; float fSliceStep = cos_theta / sin_theta; float fDistCorr; if (sin_theta > 0.0f) @@ -192,9 +165,10 @@ __global__ void FPhorizontal(float* D_projData, unsigned int projPitch, unsigned D_projData[angle*projPitch+detector] += fVal * fDistCorr; } + // projection for angles that are roughly vertical -// theta between 0 and 45, or 135 and 180 degrees -__global__ void FPvertical(float* D_projData, unsigned int projPitch, unsigned int startSlice, unsigned int startAngle, unsigned int endAngle, int regionOffset, const SDimensions dims, float outputScale) +// (detector roughly horizontal) +__global__ void FPvertical_simple(float* D_projData, unsigned int projPitch, unsigned int startSlice, unsigned int startAngle, unsigned int endAngle, const SDimensions dims, float outputScale) { const int relDet = threadIdx.x; const int relAngle = threadIdx.y; @@ -209,33 +183,7 @@ __global__ void FPvertical(float* D_projData, unsigned int projPitch, unsigned i const float sin_theta = __sinf(theta); // compute start detector for this block/angle: - // (The same for all threadIdx.x) - // ------------------------------------- - - const int midSlice = startSlice + g_blockSlices / 2; - - // project (midSlice,midRegion) on this thread's detector - - // ASSUMPTION: fDetScale >= 1.0f - // problem: detector regions get skipped because slice blocks aren't large - // enough - const unsigned int g_blockSliceSize = g_detBlockSize; - - const float fBase = 0.5f*dims.iProjDets - 0.5f + - ( - (g_blockSliceSize/2 - 0.5f*dims.iVolWidth + 0.5f) * cos_theta - - (midSlice - 0.5f*dims.iVolHeight + 0.5f) * sin_theta - + gC_angle_offset[angle] - ) / dims.fDetScale; - int iBase = (int)floorf(fBase * fPREC_FACTOR); - int iInc = (int)floorf(g_blockSliceSize * cos_theta / dims.fDetScale * fPREC_FACTOR); - - // ASSUMPTION: 16 > regionOffset / fDetScale - const int detRegion = (iBase + (blockIdx.y - regionOffset)*iInc + 16*iPREC_FACTOR*g_detBlockSize) / (iPREC_FACTOR * g_detBlockSize) - 16; - const int detPrevRegion = (iBase + (blockIdx.y - regionOffset-1)*iInc + 16*iPREC_FACTOR*g_detBlockSize) / (iPREC_FACTOR * g_detBlockSize) - 16; - - if (blockIdx.y > 0 && detRegion == detPrevRegion) - return; + const int detRegion = blockIdx.y; const int detector = detRegion * g_detBlockSize + relDet; @@ -245,7 +193,7 @@ __global__ void FPvertical(float* D_projData, unsigned int projPitch, unsigned i if (detector < 0 || detector >= dims.iProjDets) return; - const float fDetStep = dims.fDetScale / cos_theta; + const float fDetStep = gC_angle_detsize[angle] / cos_theta; float fSliceStep = sin_theta / cos_theta; float fDistCorr; if (cos_theta < 0.0f) @@ -293,183 +241,67 @@ __global__ void FPvertical(float* D_projData, unsigned int projPitch, unsigned i D_projData[angle*projPitch+detector] += fVal * fDistCorr; } -// projection for angles that are roughly horizontal -// (detector roughly vertical) -__global__ void FPhorizontal_simple(float* D_projData, unsigned int projPitch, unsigned int startSlice, unsigned int startAngle, unsigned int endAngle, const SDimensions dims, float outputScale) -{ - const int relDet = threadIdx.x; - const int relAngle = threadIdx.y; - - int angle = startAngle + blockIdx.x * g_anglesPerBlock + relAngle; - - if (angle >= endAngle) - return; - - const float theta = gC_angle[angle]; - const float cos_theta = __cosf(theta); - const float sin_theta = __sinf(theta); - - // compute start detector for this block/angle: - const int detRegion = blockIdx.y; - - const int detector = detRegion * g_detBlockSize + relDet; - - // Now project the part of the ray to angle,detector through - // slices startSlice to startSlice+g_blockSlices-1 - - if (detector < 0 || detector >= dims.iProjDets) - return; - - const float fDetStep = -dims.fDetScale / sin_theta; - float fSliceStep = cos_theta / sin_theta; - float fDistCorr; - if (sin_theta > 0.0f) - fDistCorr = -fDetStep; - else - fDistCorr = fDetStep; - fDistCorr *= outputScale; - - float fVal = 0.0f; - // project detector on slice - float fP = (detector - 0.5f*dims.iProjDets + 0.5f - gC_angle_offset[angle]) * fDetStep + (startSlice - 0.5f*dims.iVolWidth + 0.5f) * fSliceStep + 0.5f*dims.iVolHeight - 0.5f + 0.5f; - float fS = startSlice + 0.5f; - int endSlice = startSlice + g_blockSlices; - if (endSlice > dims.iVolWidth) - endSlice = dims.iVolWidth; - if (dims.iRaysPerDet > 1) { - fP += (-0.5f*dims.iRaysPerDet + 0.5f)/dims.iRaysPerDet * fDetStep; - const float fSubDetStep = fDetStep / dims.iRaysPerDet; - fDistCorr /= dims.iRaysPerDet; - fSliceStep -= dims.iRaysPerDet * fSubDetStep; +// Coordinates of center of detector pixel number t: +// x = (t - 0.5*nDets + 0.5 - fOffset) * fSize * cos(fAngle) +// y = - (t - 0.5*nDets + 0.5 - fOffset) * fSize * sin(fAngle) - for (int slice = startSlice; slice < endSlice; ++slice) - { - for (int iSubT = 0; iSubT < dims.iRaysPerDet; ++iSubT) { - fVal += tex2D(gT_volumeTexture, fS, fP); - fP += fSubDetStep; - } - fP += fSliceStep; - fS += 1.0f; - } - - } else { - for (int slice = startSlice; slice < endSlice; ++slice) - { - fVal += tex2D(gT_volumeTexture, fS, fP); - fP += fSliceStep; - fS += 1.0f; - } - - - } - - D_projData[angle*projPitch+detector] += fVal * fDistCorr; -} - - -// projection for angles that are roughly vertical -// (detector roughly horizontal) -__global__ void FPvertical_simple(float* D_projData, unsigned int projPitch, unsigned int startSlice, unsigned int startAngle, unsigned int endAngle, const SDimensions dims, float outputScale) +static void convertAndUploadAngles(const SParProjection *projs, unsigned int nth, unsigned int ndets) { - const int relDet = threadIdx.x; - const int relAngle = threadIdx.y; - - int angle = startAngle + blockIdx.x * g_anglesPerBlock + relAngle; - - if (angle >= endAngle) - return; - - const float theta = gC_angle[angle]; - const float cos_theta = __cosf(theta); - const float sin_theta = __sinf(theta); - - // compute start detector for this block/angle: - const int detRegion = blockIdx.y; - - const int detector = detRegion * g_detBlockSize + relDet; + float *angles = new float[nth]; + float *offset = new float[nth]; + float *detsizes = new float[nth]; - // Now project the part of the ray to angle,detector through - // slices startSlice to startSlice+g_blockSlices-1 + for (int i = 0; i < nth; ++i) { - if (detector < 0 || detector >= dims.iProjDets) - return; +#warning TODO: Replace by getParParamaters - const float fDetStep = dims.fDetScale / cos_theta; - float fSliceStep = sin_theta / cos_theta; - float fDistCorr; - if (cos_theta < 0.0f) - fDistCorr = -fDetStep; - else - fDistCorr = fDetStep; - fDistCorr *= outputScale; + // Take part of DetU orthogonal to Ray + double ux = projs[i].fDetUX; + double uy = projs[i].fDetUY; - float fVal = 0.0f; - float fP = (detector - 0.5f*dims.iProjDets + 0.5f - gC_angle_offset[angle]) * fDetStep + (startSlice - 0.5f*dims.iVolHeight + 0.5f) * fSliceStep + 0.5f*dims.iVolWidth - 0.5f + 0.5f; - float fS = startSlice+0.5f; - int endSlice = startSlice + g_blockSlices; - if (endSlice > dims.iVolHeight) - endSlice = dims.iVolHeight; + double t = (ux * projs[i].fRayX + uy * projs[i].fRayY) / (projs[i].fRayX * projs[i].fRayX + projs[i].fRayY * projs[i].fRayY); - if (dims.iRaysPerDet > 1) { + ux -= t * projs[i].fRayX; + uy -= t * projs[i].fRayY; - fP += (-0.5f*dims.iRaysPerDet + 0.5f)/dims.iRaysPerDet * fDetStep; - const float fSubDetStep = fDetStep / dims.iRaysPerDet; - fDistCorr /= dims.iRaysPerDet; + double angle = atan2(uy, ux); - fSliceStep -= dims.iRaysPerDet * fSubDetStep; + angles[i] = angle; - for (int slice = startSlice; slice < endSlice; ++slice) - { - for (int iSubT = 0; iSubT < dims.iRaysPerDet; ++iSubT) { - fVal += tex2D(gT_volumeTexture, fP, fS); - fP += fSubDetStep; - } - fP += fSliceStep; - fS += 1.0f; - } + double norm2 = uy * uy + ux * ux; - } else { - - for (int slice = startSlice; slice < endSlice; ++slice) - { - fVal += tex2D(gT_volumeTexture, fP, fS); - fP += fSliceStep; - fS += 1.0f; - } + detsizes[i] = sqrt(norm2); + // CHECKME: SIGNS? + offset[i] = -0.5*ndets - (projs[i].fDetSY*uy + projs[i].fDetSX*ux) / norm2; } - D_projData[angle*projPitch+detector] += fVal * fDistCorr; + cudaMemcpyToSymbol(gC_angle, angles, nth*sizeof(float), 0, cudaMemcpyHostToDevice); + cudaMemcpyToSymbol(gC_angle_offset, offset, nth*sizeof(float), 0, cudaMemcpyHostToDevice); + cudaMemcpyToSymbol(gC_angle_detsize, detsizes, nth*sizeof(float), 0, cudaMemcpyHostToDevice); } bool FP_simple_internal(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, - const SDimensions& dims, const float* angles, - const float* TOffsets, float outputScale) + const SDimensions& dims, const SParProjection* angles, + float outputScale) { - // TODO: load angles into constant memory in smaller blocks assert(dims.iProjAngles <= g_MaxAngles); + assert(angles); + cudaArray* D_dataArray; bindVolumeDataTexture(D_volumeData, D_dataArray, volumePitch, dims.iVolWidth, dims.iVolHeight); - cudaMemcpyToSymbol(gC_angle, angles, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - if (TOffsets) { - cudaMemcpyToSymbol(gC_angle_offset, TOffsets, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - } else { - if (!g_pfZeroes) { - g_pfZeroes = new float[g_MaxAngles]; - memset(g_pfZeroes, 0, g_MaxAngles * sizeof(float)); - } - cudaMemcpyToSymbol(gC_angle_offset, g_pfZeroes, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - } + convertAndUploadAngles(angles, dims.iProjAngles, dims.iProjDets); + dim3 dimBlock(g_detBlockSize, g_anglesPerBlock); // detector block size, angles @@ -492,7 +324,7 @@ bool FP_simple_internal(float* D_volumeData, unsigned int volumePitch, // Maybe we should detect corner cases and put them in the optimal // group of angles. if (a != dims.iProjAngles) - vertical = (fabsf(sinf(angles[a])) <= fabsf(cosf(angles[a]))); + vertical = (fabsf(angles[a].fRayX) <= fabsf(angles[a].fRayY)); if (a == dims.iProjAngles || vertical != blockVertical) { // block done @@ -536,8 +368,8 @@ bool FP_simple_internal(float* D_volumeData, unsigned int volumePitch, bool FP_simple(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, - const SDimensions& dims, const float* angles, - const float* TOffsets, float outputScale) + const SDimensions& dims, const SParProjection* angles, + float outputScale) { for (unsigned int iAngle = 0; iAngle < dims.iProjAngles; iAngle += g_MaxAngles) { SDimensions subdims = dims; @@ -550,7 +382,7 @@ bool FP_simple(float* D_volumeData, unsigned int volumePitch, ret = FP_simple_internal(D_volumeData, volumePitch, D_projData + iAngle * projPitch, projPitch, subdims, angles + iAngle, - TOffsets ? TOffsets + iAngle : 0, outputScale); + outputScale); if (!ret) return false; } @@ -559,106 +391,12 @@ bool FP_simple(float* D_volumeData, unsigned int volumePitch, bool FP(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, - const SDimensions& dims, const float* angles, - const float* TOffsets, float outputScale) + const SDimensions& dims, const SParProjection* angles, + float outputScale) { return FP_simple(D_volumeData, volumePitch, D_projData, projPitch, - dims, angles, TOffsets, outputScale); - - // TODO: Fix bug in this non-simple FP with large detscale and TOffsets -#if 0 - - // TODO: load angles into constant memory in smaller blocks - assert(dims.iProjAngles <= g_MaxAngles); - - // TODO: compute region size dynamically to resolve these two assumptions - // ASSUMPTION: 16 > regionOffset / fDetScale - const unsigned int g_blockSliceSize = g_detBlockSize; - assert(16 > (g_blockSlices / g_blockSliceSize) / dims.fDetScale); - // ASSUMPTION: fDetScale >= 1.0f - assert(dims.fDetScale > 0.9999f); - - cudaArray* D_dataArray; - bindVolumeDataTexture(D_volumeData, D_dataArray, volumePitch, dims.iVolWidth, dims.iVolHeight); - - cudaMemcpyToSymbol(gC_angle, angles, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - - if (TOffsets) { - cudaMemcpyToSymbol(gC_angle_offset, TOffsets, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - } else { - if (!g_pfZeroes) { - g_pfZeroes = new float[g_MaxAngles]; - memset(g_pfZeroes, 0, g_MaxAngles * sizeof(float)); - } - cudaMemcpyToSymbol(gC_angle_offset, g_pfZeroes, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); - } - - int regionOffset = g_blockSlices / g_blockSliceSize; - - dim3 dimBlock(g_detBlockSize, g_anglesPerBlock); // region size, angles - - std::list<cudaStream_t> streams; - + dims, angles, outputScale); - // Run over all angles, grouping them into groups of the same - // orientation (roughly horizontal vs. roughly vertical). - // Start a stream of grids for each such group. - - // TODO: Check if it's worth it to store this info instead - // of recomputing it every FP. - - unsigned int blockStart = 0; - unsigned int blockEnd = 0; - bool blockVertical = false; - for (unsigned int a = 0; a <= dims.iProjAngles; ++a) { - bool vertical; - // TODO: Having <= instead of < below causes a 5% speedup. - // Maybe we should detect corner cases and put them in the optimal - // group of angles. - if (a != dims.iProjAngles) - vertical = (fabsf(sinf(angles[a])) <= fabsf(cosf(angles[a]))); - if (a == dims.iProjAngles || vertical != blockVertical) { - // block done - - blockEnd = a; - if (blockStart != blockEnd) { - unsigned int length = dims.iVolHeight; - if (blockVertical) - length = dims.iVolWidth; - dim3 dimGrid((blockEnd-blockStart+g_anglesPerBlock-1)/g_anglesPerBlock, - (length+g_blockSliceSize-1)/g_blockSliceSize+2*regionOffset); // angle blocks, regions - // TODO: check if we can't immediately - // destroy the stream after use - cudaStream_t stream; - cudaStreamCreate(&stream); - streams.push_back(stream); - //printf("angle block: %d to %d, %d\n", blockStart, blockEnd, blockVertical); - if (!blockVertical) - for (unsigned int i = 0; i < dims.iVolWidth; i += g_blockSlices) - FPhorizontal<<<dimGrid, dimBlock, 0, stream>>>(D_projData, projPitch, i, blockStart, blockEnd, regionOffset, dims, outputScale); - else - for (unsigned int i = 0; i < dims.iVolHeight; i += g_blockSlices) - FPvertical<<<dimGrid, dimBlock, 0, stream>>>(D_projData, projPitch, i, blockStart, blockEnd, regionOffset, dims, outputScale); - } - blockVertical = vertical; - blockStart = a; - } - } - - for (std::list<cudaStream_t>::iterator iter = streams.begin(); iter != streams.end(); ++iter) - cudaStreamDestroy(*iter); - - streams.clear(); - - cudaThreadSynchronize(); - - cudaTextForceKernelsCompletion(); - - cudaFreeArray(D_dataArray); - - - return true; -#endif } diff --git a/cuda/2d/par_fp.h b/cuda/2d/par_fp.h index 010d064..5fd593c 100644 --- a/cuda/2d/par_fp.h +++ b/cuda/2d/par_fp.h @@ -33,8 +33,8 @@ namespace astraCUDA { _AstraExport bool FP(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, - const SDimensions& dims, const float* angles, - const float* TOffsets, float fOutputScale); + const SDimensions& dims, const SParProjection* angles, + float fOutputScale); } diff --git a/cuda/2d/sart.cu b/cuda/2d/sart.cu index c8608a3..a1085cd 100644 --- a/cuda/2d/sart.cu +++ b/cuda/2d/sart.cu @@ -254,10 +254,10 @@ bool SART::callFP_SART(float* D_volumeData, unsigned int volumePitch, { SDimensions d = dims; d.iProjAngles = 1; - if (angles) { + if (parProjs) { assert(!fanProjs); return FP(D_volumeData, volumePitch, D_projData, projPitch, - d, &angles[angle], TOffsets, outputScale); + d, &parProjs[angle], outputScale); } else { assert(fanProjs); return FanFP(D_volumeData, volumePitch, D_projData, projPitch, @@ -269,10 +269,10 @@ bool SART::callBP_SART(float* D_volumeData, unsigned int volumePitch, float* D_projData, unsigned int projPitch, unsigned int angle, float outputScale) { - if (angles) { + if (parProjs) { assert(!fanProjs); return BP_SART(D_volumeData, volumePitch, D_projData, projPitch, - angle, dims, angles, TOffsets, outputScale); + angle, dims, parProjs, outputScale); } else { assert(fanProjs); return FanBP_SART(D_volumeData, volumePitch, D_projData, projPitch, diff --git a/cuda/2d/sirt.cu b/cuda/2d/sirt.cu index 4baaccb..9dc4f64 100644 --- a/cuda/2d/sirt.cu +++ b/cuda/2d/sirt.cu @@ -152,7 +152,7 @@ bool SIRT::precomputeWeights() bool SIRT::doSlabCorrections() { // This function compensates for effectively infinitely large slab-like - // objects of finite thickness 1. + // objects of finite thickness 1 in a parallel beam geometry. // Each ray through the object has an intersection of length d/cos(alpha). // The length of the ray actually intersecting the reconstruction volume is @@ -170,6 +170,10 @@ bool SIRT::doSlabCorrections() if (useVolumeMask || useSinogramMask) return false; + // Parallel-beam only + if (!parProjs) + return false; + // multiply by line weights processSino<opDiv>(D_sinoData, D_lineWeight, projPitch, dims); @@ -181,7 +185,10 @@ bool SIRT::doSlabCorrections() float bound = cosf(1.3963f); float* t = (float*)D_sinoData; for (int i = 0; i < dims.iProjAngles; ++i) { - float f = fabs(cosf(angles[i])); + // TODO: Checkme + // TODO: Replace by getParParameters + double angle = atan2(parProjs[i].fRayX, -parProjs[i].fRayY); + float f = fabs(cosf(angle)); if (f < bound) f = bound; @@ -189,7 +196,6 @@ bool SIRT::doSlabCorrections() processSino<opMul>(t, f, sinoPitch, subdims); t += sinoPitch; } - return true; } @@ -299,40 +305,6 @@ float SIRT::computeDiffNorm() } -bool doSIRT(float* D_volumeData, unsigned int volumePitch, - float* D_sinoData, unsigned int sinoPitch, - float* D_maskData, unsigned int maskPitch, - const SDimensions& dims, const float* angles, - const float* TOffsets, unsigned int iterations) -{ - SIRT sirt; - bool ok = true; - - ok &= sirt.setGeometry(dims, angles); - if (D_maskData) - ok &= sirt.enableVolumeMask(); - if (TOffsets) - ok &= sirt.setTOffsets(TOffsets); - - if (!ok) - return false; - - ok = sirt.init(); - if (!ok) - return false; - - if (D_maskData) - ok &= sirt.setVolumeMask(D_maskData, maskPitch); - - ok &= sirt.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch); - if (!ok) - return false; - - ok = sirt.iterate(iterations); - - return ok; -} - } #ifdef STANDALONE diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu index 0e13be1..e7418a2 100644 --- a/cuda/3d/fdk.cu +++ b/cuda/3d/fdk.cu @@ -353,7 +353,7 @@ bool FDK_Filter(cudaPitchedPtr D_projData, // The filtering is a regular ramp filter per detector line. int iPaddedDetCount = calcNextPowerOfTwo(2 * dims.iProjU); - int iHalfFFTSize = calcFFTFourSize(iPaddedDetCount); + int iHalfFFTSize = astraCUDA::calcFFTFourierSize(iPaddedDetCount); int projPitch = D_projData.pitch/sizeof(float); @@ -361,22 +361,22 @@ bool FDK_Filter(cudaPitchedPtr D_projData, float* D_sinoData = (float*)D_projData.ptr; cufftComplex * D_sinoFFT = NULL; - allocateComplexOnDevice(dims.iProjAngles, iHalfFFTSize, &D_sinoFFT); + astraCUDA::allocateComplexOnDevice(dims.iProjAngles, iHalfFFTSize, &D_sinoFFT); bool ok = true; for (int v = 0; v < dims.iProjV; ++v) { - ok = runCudaFFT(dims.iProjAngles, D_sinoData, projPitch, + ok = astraCUDA::runCudaFFT(dims.iProjAngles, D_sinoData, projPitch, dims.iProjU, iPaddedDetCount, iHalfFFTSize, D_sinoFFT); if (!ok) break; - applyFilter(dims.iProjAngles, iHalfFFTSize, D_sinoFFT, D_filter); + astraCUDA::applyFilter(dims.iProjAngles, iHalfFFTSize, D_sinoFFT, D_filter); - ok = runCudaIFFT(dims.iProjAngles, D_sinoFFT, D_sinoData, projPitch, + ok = astraCUDA::runCudaIFFT(dims.iProjAngles, D_sinoFFT, D_sinoData, projPitch, dims.iProjU, iPaddedDetCount, iHalfFFTSize); if (!ok) break; @@ -384,7 +384,7 @@ bool FDK_Filter(cudaPitchedPtr D_projData, D_sinoData += (dims.iProjAngles * projPitch); } - freeComplexOnDevice(D_sinoFFT); + astraCUDA::freeComplexOnDevice(D_sinoFFT); return ok; } @@ -401,7 +401,7 @@ bool FDK(cudaPitchedPtr D_volumeData, // TODO: Check errors cufftComplex * D_filter; int iPaddedDetCount = calcNextPowerOfTwo(2 * dims.iProjU); - int iHalfFFTSize = calcFFTFourSize(iPaddedDetCount); + int iHalfFFTSize = astraCUDA::calcFFTFourierSize(iPaddedDetCount); ok = FDK_PreWeight(D_projData, fSrcOrigin, fDetOrigin, fSrcZ, fDetZ, fDetUSize, fDetVSize, @@ -412,11 +412,11 @@ 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); + astraCUDA::genFilter(astra::FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize); - allocateComplexOnDevice(dims.iProjAngles, iHalfFFTSize, &D_filter); - uploadComplexArrayToDevice(dims.iProjAngles, iHalfFFTSize, pHostFilter, D_filter); + astraCUDA::allocateComplexOnDevice(dims.iProjAngles, iHalfFFTSize, &D_filter); + astraCUDA::uploadComplexArrayToDevice(dims.iProjAngles, iHalfFFTSize, pHostFilter, D_filter); delete [] pHostFilter; @@ -430,7 +430,7 @@ bool FDK(cudaPitchedPtr D_volumeData, bShortScan, dims, angles); // Clean up filter - freeComplexOnDevice(D_filter); + astraCUDA::freeComplexOnDevice(D_filter); if (!ok) diff --git a/include/astra/CudaFilteredBackProjectionAlgorithm.h b/include/astra/CudaFilteredBackProjectionAlgorithm.h index cf1f19f..180e001 100644 --- a/include/astra/CudaFilteredBackProjectionAlgorithm.h +++ b/include/astra/CudaFilteredBackProjectionAlgorithm.h @@ -26,28 +26,24 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. $Id$ */ -#ifndef CUDAFILTEREDBACKPROJECTIONALGORITHM2_H -#define CUDAFILTEREDBACKPROJECTIONALGORITHM2_H +#ifndef CUDAFILTEREDBACKPROJECTIONALGORITHM_H +#define CUDAFILTEREDBACKPROJECTIONALGORITHM_H #include <astra/Float32ProjectionData2D.h> #include <astra/Float32VolumeData2D.h> -#include <astra/ReconstructionAlgorithm2D.h> +#include <astra/CudaReconstructionAlgorithm2D.h> #include "../../cuda/2d/astra.h" namespace astra { -class _AstraExport CCudaFilteredBackProjectionAlgorithm : public CReconstructionAlgorithm2D +class _AstraExport CCudaFilteredBackProjectionAlgorithm : public CCudaReconstructionAlgorithm2D { public: static std::string type; private: - CFloat32ProjectionData2D * m_pSinogram; - CFloat32VolumeData2D * m_pReconstruction; - int m_iGPUIndex; - int m_iPixelSuperSampling; E_FBPFILTER m_eFilter; float * m_pfFilter; int m_iFilterWidth; // number of elements per projection direction in filter @@ -64,15 +60,6 @@ public: virtual bool initialize(const Config& _cfg); bool initialize(CFloat32ProjectionData2D * _pSinogram, CFloat32VolumeData2D * _pReconstruction, E_FBPFILTER _eFilter, const float * _pfFilter = NULL, int _iFilterWidth = 0, int _iGPUIndex = -1, float _fFilterParameter = -1.0f); - virtual void run(int _iNrIterations = 0); - - static int calcIdealRealFilterWidth(int _iDetectorCount); - static int calcIdealFourierFilterWidth(int _iDetectorCount); - - //debug - static void testGenFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount); - static int getGPUCount(); - /** Get a description of the class. * * @return description string @@ -82,12 +69,7 @@ public: protected: bool check(); - AstraFBP* m_pFBP; - - bool m_bAstraFBPInit; - - void initializeFromProjector(); - virtual bool requiresProjector() const { return false; } + virtual void initCUDAAlgorithm(); }; // inline functions @@ -95,4 +77,4 @@ inline std::string CCudaFilteredBackProjectionAlgorithm::description() const { r } -#endif /* CUDAFILTEREDBACKPROJECTIONALGORITHM2_H */ +#endif /* CUDAFILTEREDBACKPROJECTIONALGORITHM_H */ diff --git a/include/astra/GeometryUtil2D.h b/include/astra/GeometryUtil2D.h index 680cecd..2f03062 100644 --- a/include/astra/GeometryUtil2D.h +++ b/include/astra/GeometryUtil2D.h @@ -32,27 +32,72 @@ $Id$ namespace astra { struct SParProjection { - // the ray direction - float fRayX, fRayY; + // the ray direction + float fRayX, fRayY; + + // the start of the (linear) detector + float fDetSX, fDetSY; + + // the length of a single detector pixel + float fDetUX, fDetUY; + + + void translate(double dx, double dy) { + fDetSX += dx; + fDetSY += dy; + } + void scale(double factor) { + fRayX *= factor; + fRayY *= factor; + fDetSX *= factor; + fDetSY *= factor; + fDetUX *= factor; + fDetUY *= factor; + } +}; - // the start of the (linear) detector - float fDetSX, fDetSY; - // the length of a single detector pixel - float fDetUX, fDetUY; +struct SFanProjection { + // the source + float fSrcX, fSrcY; + + // the start of the (linear) detector + float fDetSX, fDetSY; + + // the length of a single detector pixel + float fDetUX, fDetUY; + + void translate(double dx, double dy) { + fSrcX += dx; + fSrcY += dy; + fDetSX += dx; + fDetSY += dy; + } + void scale(double factor) { + fSrcX *= factor; + fSrcY *= factor; + fDetSX *= factor; + fDetSY *= factor; + fDetUX *= factor; + fDetUY *= factor; + } }; -struct SFanProjection { - // the source - float fSrcX, fSrcY; - // the start of the (linear) detector - float fDetSX, fDetSY; +SParProjection* genParProjections(unsigned int iProjAngles, + unsigned int iProjDets, + double fDetSize, + const float *pfAngles, + const float *pfExtraOffsets); - // the length of a single detector pixel - float fDetUX, fDetUY; -}; +SFanProjection* genFanProjections(unsigned int iProjAngles, + unsigned int iProjDets, + double fOriginSource, double fOriginDetector, + double fDetSize, + const float *pfAngles); + +bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float &fAngle, float &fOriginSource, float &fOriginDetector, float &fDetSize, float &fOffset); } diff --git a/include/astra/ParallelProjectionGeometry2D.h b/include/astra/ParallelProjectionGeometry2D.h index 36b4b6f..ca67ba7 100644 --- a/include/astra/ParallelProjectionGeometry2D.h +++ b/include/astra/ParallelProjectionGeometry2D.h @@ -84,8 +84,7 @@ public: CParallelProjectionGeometry2D(int _iProjectionAngleCount, int _iDetectorCount, float32 _fDetectorWidth, - const float32* _pfProjectionAngles, - const float32* _pfExtraDetectorOffsets = 0); + const float32* _pfProjectionAngles); /** Copy constructor. */ @@ -117,8 +116,7 @@ public: bool initialize(int _iProjectionAngleCount, int _iDetectorCount, float32 _fDetectorWidth, - const float32* _pfProjectionAngles, - const float32* _pfExtraDetectorOffsets = 0); + const float32* _pfProjectionAngles); /** Create a hard copy. */ diff --git a/include/astra/ProjectionGeometry2D.h b/include/astra/ProjectionGeometry2D.h index d26e6a7..b656d97 100644 --- a/include/astra/ProjectionGeometry2D.h +++ b/include/astra/ProjectionGeometry2D.h @@ -90,8 +90,7 @@ protected: CProjectionGeometry2D(int _iProjectionAngleCount, int _iDetectorCount, float32 _fDetectorWidth, - const float32* _pfProjectionAngles, - const float32* _pfExtraDetectorOffsets = 0); + const float32* _pfProjectionAngles); /** Copy constructor. */ @@ -117,8 +116,7 @@ protected: bool _initialize(int _iProjectionAngleCount, int _iDetectorCount, float32 _fDetectorWidth, - const float32* _pfProjectionAngles, - const float32* _pfExtraDetectorOffsets = 0); + const float32* _pfProjectionAngles); public: diff --git a/samples/matlab/s014_FBP.m b/samples/matlab/s014_FBP.m index b73149c..faf91fb 100644 --- a/samples/matlab/s014_FBP.m +++ b/samples/matlab/s014_FBP.m @@ -9,7 +9,7 @@ % ----------------------------------------------------------------------- vol_geom = astra_create_vol_geom(256, 256); -proj_geom = astra_create_proj_geom('parallel', 1.0, 384, linspace2(0,pi,180)); +proj_geom = astra_create_proj_geom('fanflat', 1.0, 384, linspace2(0,2*pi,1800), 500, 0); % As before, create a sinogram from a phantom P = phantom(256); diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index aa97eec..9d23253 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -33,6 +33,7 @@ $Id$ #include "astra/AstraObjectManager.h" #include "astra/CudaProjector2D.h" #include "../cuda/2d/astra.h" +#include "../cuda/2d/fbp.h" #include "astra/Logging.h" @@ -44,8 +45,7 @@ string CCudaFilteredBackProjectionAlgorithm::type = "FBP_CUDA"; CCudaFilteredBackProjectionAlgorithm::CCudaFilteredBackProjectionAlgorithm() { m_bIsInitialized = false; - CReconstructionAlgorithm2D::_clear(); - m_pFBP = 0; + CCudaReconstructionAlgorithm2D::_clear(); m_pfFilter = NULL; m_fFilterParameter = -1.0f; m_fFilterD = 1.0f; @@ -53,35 +53,8 @@ CCudaFilteredBackProjectionAlgorithm::CCudaFilteredBackProjectionAlgorithm() CCudaFilteredBackProjectionAlgorithm::~CCudaFilteredBackProjectionAlgorithm() { - if(m_pfFilter != NULL) - { - delete [] m_pfFilter; - m_pfFilter = NULL; - } - - if(m_pFBP != NULL) - { - delete m_pFBP; - m_pFBP = NULL; - } -} - -void CCudaFilteredBackProjectionAlgorithm::initializeFromProjector() -{ - m_iPixelSuperSampling = 1; - m_iGPUIndex = -1; - - // Projector - CCudaProjector2D* pCudaProjector = dynamic_cast<CCudaProjector2D*>(m_pProjector); - if (!pCudaProjector) { - if (m_pProjector) { - ASTRA_WARN("non-CUDA Projector2D passed to FBP_CUDA"); - } - } else { - m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling(); - m_iGPUIndex = pCudaProjector->getGPUIndex(); - } - + delete[] m_pfFilter; + m_pfFilter = NULL; } bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) @@ -92,54 +65,28 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) // if already initialized, clear first if (m_bIsInitialized) { +#warning FIXME Necessary? clear(); } - // Projector - XMLNode node = _cfg.self.getSingleNode("ProjectorId"); - CCudaProjector2D* pCudaProjector = 0; - if (node) { - int id = node.getContentInt(); - CProjector2D *projector = CProjector2DManager::getSingleton().get(id); - pCudaProjector = dynamic_cast<CCudaProjector2D*>(projector); - if (!pCudaProjector) { - ASTRA_WARN("non-CUDA Projector2D passed"); - } - } - CC.markNodeParsed("ProjectorId"); - - - // sinogram data - node = _cfg.self.getSingleNode("ProjectionDataId"); - ASTRA_CONFIG_CHECK(node, "CudaFBP", "No ProjectionDataId tag specified."); - int id = node.getContentInt(); - m_pSinogram = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id)); - CC.markNodeParsed("ProjectionDataId"); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_cfg); + if (!m_bIsInitialized) + return false; - // reconstruction data - node = _cfg.self.getSingleNode("ReconstructionDataId"); - ASTRA_CONFIG_CHECK(node, "CudaFBP", "No ReconstructionDataId tag specified."); - id = node.getContentInt(); - m_pReconstruction = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id)); - CC.markNodeParsed("ReconstructionDataId"); // filter type - node = _cfg.self.getSingleNode("FilterType"); + XMLNode node = _cfg.self.getSingleNode("FilterType"); if (node) - { m_eFilter = _convertStringToFilter(node.getContent().c_str()); - } else - { m_eFilter = FILTER_RAMLAK; - } CC.markNodeParsed("FilterType"); // filter node = _cfg.self.getSingleNode("FilterSinogramId"); if (node) { - id = node.getContentInt(); + int id = node.getContentInt(); const CFloat32ProjectionData2D * pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id)); m_iFilterWidth = pFilterData->getGeometry()->getDetectorCount(); int iFilterProjectionCount = pFilterData->getGeometry()->getProjectionAngleCount(); @@ -188,20 +135,9 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) initializeFromProjector(); - // Deprecated options - m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling); - CC.markOptionParsed("PixelSuperSampling"); - // GPU number - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); - CC.markOptionParsed("GPUIndex"); - if (!_cfg.self.hasOption("GPUIndex")) - CC.markOptionParsed("GPUindex"); - - - m_pFBP = new AstraFBP; - m_bAstraFBPInit = false; + m_pAlgo = new astraCUDA::FBP(); + m_bAlgoInit = false; return check(); } @@ -226,9 +162,8 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * // success m_bIsInitialized = true; - m_pFBP = new AstraFBP; - - m_bAstraFBPInit = false; + m_pAlgo = new astraCUDA::FBP(); + m_bAlgoInit = false; if(_pfFilter != NULL) { @@ -256,107 +191,26 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * return check(); } -void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) -{ - // check initialized - ASTRA_ASSERT(m_bIsInitialized); - - if (!m_bAstraFBPInit) { - - const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); - const CParallelProjectionGeometry2D* parprojgeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry()); - const CFanFlatProjectionGeometry2D* fanprojgeom = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry()); - - bool ok = true; - - // TODO: non-square pixels? - ok &= m_pFBP->setReconstructionGeometry(volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getPixelLengthX()); - int iDetectorCount; - if (parprojgeom) { - - float *offsets, *angles, detSize, outputScale; - - ok = convertAstraGeometry(&volgeom, parprojgeom, offsets, angles, detSize, outputScale); - - - ok &= m_pFBP->setProjectionGeometry(parprojgeom->getProjectionAngleCount(), - parprojgeom->getDetectorCount(), - angles, - parprojgeom->getDetectorWidth()); - iDetectorCount = parprojgeom->getDetectorCount(); - // TODO: Are detSize and outputScale handled correctly? - - if (offsets) - ok &= m_pFBP->setTOffsets(offsets); - ASTRA_ASSERT(ok); - - delete[] offsets; - delete[] angles; - - } else if (fanprojgeom) { - - astraCUDA::SFanProjection* projs; - float outputScale; - - // FIXME: Implement this, and clean up the interface to AstraFBP. - if (abs(volgeom.getWindowMinX() + volgeom.getWindowMaxX()) > 0.00001 * volgeom.getPixelLengthX()) { - // Off-center volume geometry isn't supported yet - ASTRA_ASSERT(false); - } - if (abs(volgeom.getWindowMinY() + volgeom.getWindowMaxY()) > 0.00001 * volgeom.getPixelLengthY()) { - // Off-center volume geometry isn't supported yet - ASTRA_ASSERT(false); - } - - ok = convertAstraGeometry(&volgeom, fanprojgeom, projs, outputScale); - - // CHECKME: outputScale? - - ok &= m_pFBP->setFanGeometry(fanprojgeom->getProjectionAngleCount(), - fanprojgeom->getDetectorCount(), - projs, - fanprojgeom->getProjectionAngles(), - fanprojgeom->getOriginSourceDistance(), - fanprojgeom->getOriginDetectorDistance(), - - fanprojgeom->getDetectorWidth(), - m_bShortScan); - - iDetectorCount = fanprojgeom->getDetectorCount(); - - delete[] projs; - } else { - assert(false); - } - - ok &= m_pFBP->setPixelSuperSampling(m_iPixelSuperSampling); - - ASTRA_ASSERT(ok); +void CCudaFilteredBackProjectionAlgorithm::initCUDAAlgorithm() +{ + CCudaReconstructionAlgorithm2D::initCUDAAlgorithm(); - ok &= m_pFBP->init(m_iGPUIndex); - ASTRA_ASSERT(ok); + astraCUDA::FBP* pFBP = dynamic_cast<astraCUDA::FBP*>(m_pAlgo); - ok &= m_pFBP->setSinogram(m_pSinogram->getDataConst(), iDetectorCount); + bool ok = pFBP->setFilter(m_eFilter, m_pfFilter, m_iFilterWidth, m_fFilterD, m_fFilterParameter); + if (!ok) { + ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set filter"); ASTRA_ASSERT(ok); - - ok &= m_pFBP->setFilter(m_eFilter, m_pfFilter, m_iFilterWidth, m_fFilterD, m_fFilterParameter); - ASTRA_ASSERT(ok); - - m_bAstraFBPInit = true; } - bool ok = m_pFBP->run(); - ASTRA_ASSERT(ok); - - const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); - ok &= m_pFBP->getReconstruction(m_pReconstruction->getData(), volgeom.getGridColCount()); - - ASTRA_ASSERT(ok); + ok &= pFBP->setShortScan(m_bShortScan); + if (!ok) { + ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set short-scan mode"); + } } + bool CCudaFilteredBackProjectionAlgorithm::check() { // check pointers @@ -395,16 +249,6 @@ static int calcNextPowerOfTwo(int _iValue) return iOutput; } -int CCudaFilteredBackProjectionAlgorithm::calcIdealRealFilterWidth(int _iDetectorCount) -{ - return calcNextPowerOfTwo(_iDetectorCount); -} - -int CCudaFilteredBackProjectionAlgorithm::calcIdealFourierFilterWidth(int _iDetectorCount) -{ - return (calcNextPowerOfTwo(_iDetectorCount) / 2 + 1); -} - static bool stringCompareLowerCase(const char * _stringA, const char * _stringB) { int iCmpReturn = 0; @@ -518,12 +362,3 @@ E_FBPFILTER CCudaFilteredBackProjectionAlgorithm::_convertStringToFilter(const c return output; } -void CCudaFilteredBackProjectionAlgorithm::testGenFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount) -{ - genFilter(_eFilter, _fD, _iProjectionCount, _pFilter, _iFFTRealDetectorCount, _iFFTFourierDetectorCount); -} - -int CCudaFilteredBackProjectionAlgorithm::getGPUCount() -{ - return 0; -} diff --git a/src/CudaForwardProjectionAlgorithm.cpp b/src/CudaForwardProjectionAlgorithm.cpp index 80f2e02..bdd7dd0 100644 --- a/src/CudaForwardProjectionAlgorithm.cpp +++ b/src/CudaForwardProjectionAlgorithm.cpp @@ -221,56 +221,48 @@ void CCudaForwardProjectionAlgorithm::run(int) // check initialized assert(m_bIsInitialized); - CVolumeGeometry2D* pVolGeom = m_pVolume->getGeometry(); - const CParallelProjectionGeometry2D* parProjGeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry()); - const CFanFlatProjectionGeometry2D* fanProjGeom = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry()); - const CFanFlatVecProjectionGeometry2D* fanVecProjGeom = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pSinogram->getGeometry()); + bool ok; - bool ok = false; - if (parProjGeom) { + const CVolumeGeometry2D* pVolGeom = m_pVolume->getGeometry(); + const CProjectionGeometry2D* pProjGeom = m_pSinogram->getGeometry(); + astraCUDA::SDimensions dims; - float *offsets, *angles, detSize, outputScale; - ok = convertAstraGeometry(pVolGeom, parProjGeom, offsets, angles, detSize, outputScale); + ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims); - ASTRA_ASSERT(ok); // FIXME + if (!ok) + return; - // FIXME: Output scaling + astraCUDA::SParProjection* pParProjs = 0; + astraCUDA::SFanProjection* pFanProjs = 0; + float fOutputScale = 1.0f; - ok = astraCudaFP(m_pVolume->getDataConst(), m_pSinogram->getData(), - pVolGeom->getGridColCount(), pVolGeom->getGridRowCount(), - parProjGeom->getProjectionAngleCount(), - parProjGeom->getDetectorCount(), - angles, offsets, detSize, - m_iDetectorSuperSampling, 1.0f * outputScale, m_iGPUIndex); - - delete[] offsets; - delete[] angles; + ok = convertAstraGeometry(pVolGeom, pProjGeom, pParProjs, pFanProjs, fOutputScale); + if (!ok) + return; - } else if (fanProjGeom || fanVecProjGeom) { + if (pParProjs) { + assert(!pFanProjs); - astraCUDA::SFanProjection* projs; - float outputScale; + ok = astraCudaFP(m_pVolume->getDataConst(), m_pSinogram->getData(), + pVolGeom->getGridColCount(), pVolGeom->getGridRowCount(), + pProjGeom->getProjectionAngleCount(), + pProjGeom->getDetectorCount(), + pParProjs, + m_iDetectorSuperSampling, 1.0f * fOutputScale, m_iGPUIndex); - if (fanProjGeom) { - ok = convertAstraGeometry(pVolGeom, fanProjGeom, projs, outputScale); - } else { - ok = convertAstraGeometry(pVolGeom, fanVecProjGeom, projs, outputScale); - } + delete[] pParProjs; - ASTRA_ASSERT(ok); + } else { + assert(pFanProjs); ok = astraCudaFanFP(m_pVolume->getDataConst(), m_pSinogram->getData(), pVolGeom->getGridColCount(), pVolGeom->getGridRowCount(), - m_pSinogram->getGeometry()->getProjectionAngleCount(), - m_pSinogram->getGeometry()->getDetectorCount(), - projs, - m_iDetectorSuperSampling, outputScale, m_iGPUIndex); - - delete[] projs; - - } else { + pProjGeom->getProjectionAngleCount(), + pProjGeom->getDetectorCount(), + pFanProjs, + m_iDetectorSuperSampling, fOutputScale, m_iGPUIndex); - ASTRA_ASSERT(false); + delete[] pFanProjs; } diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 2798434..326ef1e 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -250,69 +250,14 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry() ok = m_pAlgo->setGPUIndex(m_iGPUIndex); if (!ok) return false; - astraCUDA::SDimensions dims; - const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); + const CProjectionGeometry2D& projgeom = *m_pSinogram->getGeometry(); - // TODO: non-square pixels? - dims.iVolWidth = volgeom.getGridColCount(); - dims.iVolHeight = volgeom.getGridRowCount(); - float fPixelSize = volgeom.getPixelLengthX(); - - dims.iRaysPerDet = m_iDetectorSuperSampling; - dims.iRaysPerPixelDim = m_iPixelSuperSampling; - - - const CParallelProjectionGeometry2D* parProjGeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry()); - const CFanFlatProjectionGeometry2D* fanProjGeom = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry()); - const CFanFlatVecProjectionGeometry2D* fanVecProjGeom = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pSinogram->getGeometry()); - - if (parProjGeom) { - - float *offsets, *angles, detSize, outputScale; - - ok = convertAstraGeometry(&volgeom, parProjGeom, offsets, angles, detSize, outputScale); - - dims.iProjAngles = parProjGeom->getProjectionAngleCount(); - dims.iProjDets = parProjGeom->getDetectorCount(); - dims.fDetScale = parProjGeom->getDetectorWidth() / fPixelSize; - - ok = m_pAlgo->setGeometry(dims, parProjGeom->getProjectionAngles()); - ok &= m_pAlgo->setTOffsets(offsets); - - // CHECKME: outputScale? detSize? - - delete[] offsets; - delete[] angles; - - } else if (fanProjGeom || fanVecProjGeom) { - - astraCUDA::SFanProjection* projs; - float outputScale; - - if (fanProjGeom) { - ok = convertAstraGeometry(&volgeom, fanProjGeom, projs, outputScale); - } else { - ok = convertAstraGeometry(&volgeom, fanVecProjGeom, projs, outputScale); - } - - dims.iProjAngles = m_pSinogram->getGeometry()->getProjectionAngleCount(); - dims.iProjDets = m_pSinogram->getGeometry()->getDetectorCount(); - dims.fDetScale = m_pSinogram->getGeometry()->getDetectorWidth() / fPixelSize; - - ok = m_pAlgo->setFanGeometry(dims, projs); - - // CHECKME: outputScale? - - delete[] projs; - - } else { - - ASTRA_ASSERT(false); - - } + ok = m_pAlgo->setGeometry(&volgeom, &projgeom); if (!ok) return false; + ok = m_pAlgo->setSuperSampling(m_iDetectorSuperSampling, m_iPixelSuperSampling); + if (!ok) return false; if (m_bUseReconstructionMask) ok &= m_pAlgo->enableVolumeMask(); diff --git a/src/FanFlatProjectionGeometry2D.cpp b/src/FanFlatProjectionGeometry2D.cpp index 3aab582..4eec9c4 100644 --- a/src/FanFlatProjectionGeometry2D.cpp +++ b/src/FanFlatProjectionGeometry2D.cpp @@ -28,6 +28,8 @@ $Id$ #include "astra/FanFlatProjectionGeometry2D.h" +#include "astra/GeometryUtil2D.h" + #include <cstring> #include <sstream> @@ -218,16 +220,12 @@ Config* CFanFlatProjectionGeometry2D::getConfiguration() const //---------------------------------------------------------------------------------------- CFanFlatVecProjectionGeometry2D* CFanFlatProjectionGeometry2D::toVectorGeometry() { - SFanProjection* vectors = new SFanProjection[m_iProjectionAngleCount]; - for (int i = 0; i < m_iProjectionAngleCount; ++i) - { - vectors[i].fSrcX = sinf(m_pfProjectionAngles[i]) * m_fOriginSourceDistance; - vectors[i].fSrcY = -cosf(m_pfProjectionAngles[i]) * m_fOriginSourceDistance; - vectors[i].fDetUX = cosf(m_pfProjectionAngles[i]) * m_fDetectorWidth; - vectors[i].fDetUY = sinf(m_pfProjectionAngles[i]) * m_fDetectorWidth; - vectors[i].fDetSX = -sinf(m_pfProjectionAngles[i]) * m_fOriginDetectorDistance - 0.5f * m_iDetectorCount * vectors[i].fDetUX; - vectors[i].fDetSY = cosf(m_pfProjectionAngles[i]) * m_fOriginDetectorDistance - 0.5f * m_iDetectorCount * vectors[i].fDetUY; - } + SFanProjection* vectors = genFanProjections(m_iProjectionAngleCount, + m_iDetectorCount, + m_fOriginSourceDistance, + m_fOriginDetectorDistance, + m_fDetectorWidth, + m_pfProjectionAngles); CFanFlatVecProjectionGeometry2D* vecGeom = new CFanFlatVecProjectionGeometry2D(); vecGeom->initialize(m_iProjectionAngleCount, m_iDetectorCount, vectors); diff --git a/src/GeometryUtil2D.cpp b/src/GeometryUtil2D.cpp new file mode 100644 index 0000000..1bca2bc --- /dev/null +++ b/src/GeometryUtil2D.cpp @@ -0,0 +1,161 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp + 2014-2015, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://sf.net/projects/astra-toolbox + +This file is part of the ASTRA Toolbox. + + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +$Id$ +*/ + +#include "astra/GeometryUtil2D.h" + +#include <cmath> + +namespace astra { + +SParProjection* genParProjections(unsigned int iProjAngles, + unsigned int iProjDets, + double fDetSize, + const float *pfAngles, + const float *pfExtraOffsets) +{ + SParProjection base; + base.fRayX = 0.0f; + base.fRayY = 1.0f; + + base.fDetSX = iProjDets * fDetSize * -0.5f; + base.fDetSY = 0.0f; + + base.fDetUX = fDetSize; + base.fDetUY = 0.0f; + + SParProjection* p = new SParProjection[iProjAngles]; + +#define ROTATE0(name,i,alpha) do { p[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); p[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); } while(0) + + for (unsigned int i = 0; i < iProjAngles; ++i) { + if (pfExtraOffsets) { + // TODO + } + + ROTATE0(Ray, i, pfAngles[i]); + ROTATE0(DetS, i, pfAngles[i]); + ROTATE0(DetU, i, pfAngles[i]); + + if (pfExtraOffsets) { + float d = pfExtraOffsets[i]; + p[i].fDetSX -= d * p[i].fDetUX; + p[i].fDetSY -= d * p[i].fDetUY; + } + } + +#undef ROTATE0 + + return p; +} + + +SFanProjection* genFanProjections(unsigned int iProjAngles, + unsigned int iProjDets, + double fOriginSource, double fOriginDetector, + double fDetSize, + const float *pfAngles) +// const float *pfExtraOffsets) +{ + SFanProjection *pProjs = new SFanProjection[iProjAngles]; + + float fSrcX0 = 0.0f; + float fSrcY0 = -fOriginSource; + float fDetUX0 = fDetSize; + float fDetUY0 = 0.0f; + float fDetSX0 = iProjDets * fDetUX0 / -2.0f; + float fDetSY0 = fOriginDetector; + +#define ROTATE0(name,i,alpha) do { pProjs[i].f##name##X = f##name##X0 * cos(alpha) - f##name##Y0 * sin(alpha); pProjs[i].f##name##Y = f##name##X0 * sin(alpha) + f##name##Y0 * cos(alpha); } while(0) + for (unsigned int i = 0; i < iProjAngles; ++i) { + ROTATE0(Src, i, pfAngles[i]); + ROTATE0(DetS, i, pfAngles[i]); + ROTATE0(DetU, i, pfAngles[i]); + } + +#undef ROTATE0 + + return pProjs; +} + +// Convert a SParProjection back into its set of "standard" circular parallel +// beam parameters. This is always possible. +bool getParParameters(const SParProjection &proj /* , ... */) +{ + // angle + // det size + // offset + + // (see convertAndUploadAngles in par_fp.cu) + + return true; +} + +// Convert a SFanProjection back into its set of "standard" circular fan beam +// parameters. This will return false if it can not be represented in this way. +bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float &fAngle, float &fOriginSource, float &fOriginDetector, float &fDetSize, float &fOffset) +{ + // angle + // det size + // offset + // origin-source + // origin-detector + + // Need to check if line source-origin is orthogonal to vector ux,uy + // (including the case source==origin) + + // (equivalent: source and origin project to same point on detector) + + double dp = proj.fSrcX * proj.fDetUX + proj.fSrcY * proj.fDetUY; + + double rel = (proj.fSrcX*proj.fSrcX + proj.fSrcY*proj.fSrcY) * (proj.fDetUX*proj.fDetUX + proj.fDetUY*proj.fDetUY); + rel = sqrt(rel); + + if (std::abs(dp) > rel * 0.0001) + return false; + + fOriginSource = sqrt(proj.fSrcX*proj.fSrcX + proj.fSrcY*proj.fSrcY); + + fDetSize = sqrt(proj.fDetUX*proj.fDetUX + proj.fDetUY*proj.fDetUY); + + // project origin on detector line ( == project source on detector line) + + double t = (- proj.fDetSX) * proj.fDetUX + (- proj.fDetSY) * proj.fDetUY; + + fOffset = (float)t - 0.5*iProjDets; + + // TODO: CHECKME + fOriginDetector = sqrt((proj.fDetSX + t * proj.fDetUX)*(proj.fDetSX + t * proj.fDetUX) + (proj.fDetSY + t * proj.fDetUY)*(proj.fDetSY + t * proj.fDetUY)); + + //float fAngle = atan2(proj.fDetSX + t * proj.fDetUX - proj.fSrcX, proj.fDetSY + t * proj.fDetUY); // TODO: Fix order + sign + fAngle = atan2(proj.fDetUY, proj.fDetUX); // TODO: Check order + sign + + return true; +} + + +} diff --git a/src/ParallelProjectionGeometry2D.cpp b/src/ParallelProjectionGeometry2D.cpp index 105e24b..e1f93aa 100644 --- a/src/ParallelProjectionGeometry2D.cpp +++ b/src/ParallelProjectionGeometry2D.cpp @@ -28,6 +28,8 @@ $Id$ #include "astra/ParallelProjectionGeometry2D.h" +#include "astra/GeometryUtil2D.h" + #include <cstring> using namespace std; @@ -48,15 +50,13 @@ CParallelProjectionGeometry2D::CParallelProjectionGeometry2D() : CParallelProjectionGeometry2D::CParallelProjectionGeometry2D(int _iProjectionAngleCount, int _iDetectorCount, float32 _fDetectorWidth, - const float32* _pfProjectionAngles, - const float32* _pfExtraDetectorOffsets) + const float32* _pfProjectionAngles) { _clear(); initialize(_iProjectionAngleCount, _iDetectorCount, _fDetectorWidth, - _pfProjectionAngles, - _pfExtraDetectorOffsets); + _pfProjectionAngles); } //---------------------------------------------------------------------------------------- @@ -115,14 +115,12 @@ bool CParallelProjectionGeometry2D::initialize(const Config& _cfg) bool CParallelProjectionGeometry2D::initialize(int _iProjectionAngleCount, int _iDetectorCount, float32 _fDetectorWidth, - const float32* _pfProjectionAngles, - const float32* _pfExtraDetectorOffsets) + const float32* _pfProjectionAngles) { _initialize(_iProjectionAngleCount, _iDetectorCount, _fDetectorWidth, - _pfProjectionAngles, - _pfExtraDetectorOffsets); + _pfProjectionAngles); // success m_bInitialized = _check(); @@ -178,11 +176,6 @@ Config* CParallelProjectionGeometry2D::getConfiguration() const cfg->self.addChildNode("DetectorCount", getDetectorCount()); cfg->self.addChildNode("DetectorWidth", getDetectorWidth()); cfg->self.addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); - if(m_pfExtraDetectorOffset!=NULL){ - XMLNode opt = cfg->self.addChildNode("Option"); - opt.addAttribute("key","ExtraDetectorOffset"); - opt.setContent(m_pfExtraDetectorOffset, m_iProjectionAngleCount); - } return cfg; } @@ -203,17 +196,11 @@ CVector3D CParallelProjectionGeometry2D::getProjectionDirection(int _iProjection //---------------------------------------------------------------------------------------- CParallelVecProjectionGeometry2D* CParallelProjectionGeometry2D::toVectorGeometry() { - SParProjection* vectors = new SParProjection[m_iProjectionAngleCount]; - for (int i = 0; i < m_iProjectionAngleCount; ++i) - { - vectors[i].fRayX = sinf(m_pfProjectionAngles[i]); - vectors[i].fRayY = -cosf(m_pfProjectionAngles[i]); - vectors[i].fDetUX = cosf(m_pfProjectionAngles[i]) * m_fDetectorWidth; - vectors[i].fDetUY = sinf(m_pfProjectionAngles[i]) * m_fDetectorWidth; - vectors[i].fDetSX = -0.5f * m_iDetectorCount * vectors[i].fDetUX; - vectors[i].fDetSY = -0.5f * m_iDetectorCount * vectors[i].fDetUY; - } - + SParProjection* vectors = genParProjections(m_iProjectionAngleCount, + m_iDetectorCount, + m_fDetectorWidth, + m_pfProjectionAngles, 0); + // TODO: ExtraOffsets? CParallelVecProjectionGeometry2D* vecGeom = new CParallelVecProjectionGeometry2D(); vecGeom->initialize(m_iProjectionAngleCount, m_iDetectorCount, vectors); delete[] vectors; diff --git a/src/ProjectionGeometry2D.cpp b/src/ProjectionGeometry2D.cpp index 03fbd22..a306b48 100644 --- a/src/ProjectionGeometry2D.cpp +++ b/src/ProjectionGeometry2D.cpp @@ -45,11 +45,10 @@ CProjectionGeometry2D::CProjectionGeometry2D() : configCheckData(0) CProjectionGeometry2D::CProjectionGeometry2D(int _iAngleCount, int _iDetectorCount, float32 _fDetectorWidth, - const float32* _pfProjectionAngles, - const float32* _pfExtraDetectorOffsets) : configCheckData(0) + const float32* _pfProjectionAngles) : configCheckData(0) { _clear(); - _initialize(_iAngleCount, _iDetectorCount, _fDetectorWidth, _pfProjectionAngles,_pfExtraDetectorOffsets); + _initialize(_iAngleCount, _iDetectorCount, _fDetectorWidth, _pfProjectionAngles); } //---------------------------------------------------------------------------------------- @@ -156,8 +155,7 @@ bool CProjectionGeometry2D::initialize(const Config& _cfg) bool CProjectionGeometry2D::_initialize(int _iProjectionAngleCount, int _iDetectorCount, float32 _fDetectorWidth, - const float32* _pfProjectionAngles, - const float32* _pfExtraDetectorOffsets) + const float32* _pfProjectionAngles) { if (m_bInitialized) { clear(); |