diff options
46 files changed, 1328 insertions, 1370 deletions
diff --git a/astra_vc14.vcxproj b/astra_vc14.vcxproj index a28bbfa..36d37ec 100644 --- a/astra_vc14.vcxproj +++ b/astra_vc14.vcxproj @@ -512,12 +512,15 @@ <ClCompile Include="src\Float32Data.cpp" /> <ClCompile Include="src\Float32Data2D.cpp" /> <ClCompile Include="src\Float32Data3D.cpp" /> + <ClCompile Include="src\Float32Data3DGPU.cpp" /> <ClCompile Include="src\Float32Data3DMemory.cpp" /> <ClCompile Include="src\Float32ProjectionData2D.cpp" /> <ClCompile Include="src\Float32ProjectionData3D.cpp" /> + <ClCompile Include="src\Float32ProjectionData3DGPU.cpp" /> <ClCompile Include="src\Float32ProjectionData3DMemory.cpp" /> <ClCompile Include="src\Float32VolumeData2D.cpp" /> <ClCompile Include="src\Float32VolumeData3D.cpp" /> + <ClCompile Include="src\Float32VolumeData3DGPU.cpp" /> <ClCompile Include="src\Float32VolumeData3DMemory.cpp" /> <ClCompile Include="src\ForwardProjectionAlgorithm.cpp" /> <ClCompile Include="src\Fourier.cpp" /> @@ -624,12 +627,15 @@ <ClInclude Include="include\astra\Float32Data.h" /> <ClInclude Include="include\astra\Float32Data2D.h" /> <ClInclude Include="include\astra\Float32Data3D.h" /> + <ClInclude Include="include\astra\Float32Data3DGPU.h" /> <ClInclude Include="include\astra\Float32Data3DMemory.h" /> <ClInclude Include="include\astra\Float32ProjectionData2D.h" /> <ClInclude Include="include\astra\Float32ProjectionData3D.h" /> + <ClInclude Include="include\astra\Float32ProjectionData3DGPU.h" /> <ClInclude Include="include\astra\Float32ProjectionData3DMemory.h" /> <ClInclude Include="include\astra\Float32VolumeData2D.h" /> <ClInclude Include="include\astra\Float32VolumeData3D.h" /> + <ClInclude Include="include\astra\Float32VolumeData3DGPU.h" /> <ClInclude Include="include\astra\Float32VolumeData3DMemory.h" /> <ClInclude Include="include\astra\ForwardProjectionAlgorithm.h" /> <ClInclude Include="include\astra\Fourier.h" /> diff --git a/astra_vc14.vcxproj.filters b/astra_vc14.vcxproj.filters index dd7f574..591a4c7 100644 --- a/astra_vc14.vcxproj.filters +++ b/astra_vc14.vcxproj.filters @@ -321,6 +321,15 @@ <ClCompile Include="src\CudaSirtAlgorithm3D.cpp"> <Filter>CUDA\astra source</Filter> </ClCompile> + <ClCompile Include="src\Float32Data3DGPU.cpp"> + <Filter>CUDA\astra source</Filter> + </ClCompile> + <ClCompile Include="src\Float32ProjectionData3DGPU.cpp"> + <Filter>CUDA\astra source</Filter> + </ClCompile> + <ClCompile Include="src\Float32VolumeData3DGPU.cpp"> + <Filter>CUDA\astra source</Filter> + </ClCompile> </ItemGroup> <ItemGroup> <ClInclude Include="include\astra\Algorithm.h"> @@ -581,6 +590,15 @@ <ClInclude Include="include\astra\CudaSirtAlgorithm3D.h"> <Filter>CUDA\astra headers</Filter> </ClInclude> + <ClInclude Include="include\astra\Float32Data3DGPU.h"> + <Filter>CUDA\astra headers</Filter> + </ClInclude> + <ClInclude Include="include\astra\Float32ProjectionData3DGPU.h"> + <Filter>CUDA\astra headers</Filter> + </ClInclude> + <ClInclude Include="include\astra\Float32VolumeData3DGPU.h"> + <Filter>CUDA\astra headers</Filter> + </ClInclude> <ClInclude Include="cuda\2d\algo.h"> <Filter>CUDA\cuda headers</Filter> </ClInclude> diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index 371b656..9066f0a 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -199,7 +199,10 @@ CUDA_CXX_OBJECTS=\ src/CudaFDKAlgorithm3D.lo \ src/CudaSirtAlgorithm3D.lo \ src/CudaBackProjectionAlgorithm3D.lo \ - src/CudaForwardProjectionAlgorithm3D.lo + src/CudaForwardProjectionAlgorithm3D.lo \ + src/Float32Data3DGPU.lo \ + src/Float32ProjectionData3DGPU.lo \ + src/Float32VolumeData3DGPU.lo CUDA_OBJECTS=\ cuda/2d/algo.lo \ diff --git a/build/msvc/gen.py b/build/msvc/gen.py index 8a40c45..9c14ffe 100644 --- a/build/msvc/gen.py +++ b/build/msvc/gen.py @@ -274,6 +274,9 @@ P_astra["filters"]["CUDA\\astra source"] = [ "src\\CudaSartAlgorithm.cpp", "src\\CudaSirtAlgorithm.cpp", "src\\CudaSirtAlgorithm3D.cpp", +"src\\Float32Data3DGPU.cpp", +"src\\Float32ProjectionData3DGPU.cpp", +"src\\Float32VolumeData3DGPU.cpp", ] P_astra["filters"]["CUDA\\cuda headers"] = [ "4e17872e-db7d-41bc-9760-fad1c253b583", @@ -411,7 +414,9 @@ P_astra["filters"]["CUDA\\astra headers"] = [ "include\\astra\\CudaSartAlgorithm.h", "include\\astra\\CudaSirtAlgorithm.h", "include\\astra\\CudaSirtAlgorithm3D.h", - +"include\\astra\\Float32Data3DGPU.h", +"include\\astra\\Float32ProjectionData3DGPU.h", +"include\\astra\\Float32VolumeData3DGPU.h", ] P_astra["filters"]["Projectors\\inline"] = [ "0daffd63-ba49-4a5f-8d7a-5322e0e74f22", diff --git a/cuda/3d/mem3d.cu b/cuda/3d/mem3d.cu index 2b26fe1..97be8a4 100644 --- a/cuda/3d/mem3d.cu +++ b/cuda/3d/mem3d.cu @@ -118,6 +118,13 @@ MemHandle3D allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, Me return ret; } +bool zeroGPUMemory(MemHandle3D handle, unsigned int x, unsigned int y, unsigned int z) +{ + SMemHandle3D_internal& hnd = *handle.d.get(); + cudaError_t err = cudaMemset3D(hnd.ptr, 0, make_cudaExtent(sizeof(float)*x, y, z)); + return err == cudaSuccess; +} + bool freeGPUMemory(MemHandle3D handle) { size_t free = availableGPUMemory(); @@ -307,6 +314,23 @@ bool FDK(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, co } +MemHandle3D wrapHandle(float *D_ptr, unsigned int x, unsigned int y, unsigned int z, unsigned int pitch) +{ + cudaPitchedPtr ptr; + ptr.ptr = D_ptr; + ptr.xsize = sizeof(float) * x; + ptr.pitch = sizeof(float) * pitch; + ptr.ysize = y; + + SMemHandle3D_internal h; + h.ptr = ptr; + + MemHandle3D hnd; + hnd.d = boost::shared_ptr<SMemHandle3D_internal>(new SMemHandle3D_internal); + *hnd.d = h; + + return hnd; +} diff --git a/cuda/3d/mem3d.h b/cuda/3d/mem3d.h index a0829e2..7a87ae6 100644 --- a/cuda/3d/mem3d.h +++ b/cuda/3d/mem3d.h @@ -80,6 +80,8 @@ enum Mem3DZeroMode { size_t availableGPUMemory(); int maxBlockDimension(); +MemHandle3D wrapHandle(float *D_ptr, unsigned int x, unsigned int y, unsigned int z, unsigned int pitch); + MemHandle3D allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, Mem3DZeroMode zero); bool copyToGPUMemory(const float *src, MemHandle3D dst, const SSubDimensions3D &pos); @@ -88,6 +90,8 @@ bool copyFromGPUMemory(float *dst, MemHandle3D src, const SSubDimensions3D &pos) bool freeGPUMemory(MemHandle3D handle); +bool zeroGPUMemory(MemHandle3D handle, unsigned int x, unsigned int y, unsigned int z); + bool setGPUIndex(int index); @@ -97,7 +101,6 @@ bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, con bool FDK(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, bool bShortScan, const float *pfFilter = 0); - } #endif diff --git a/include/astra/CompositeGeometryManager.h b/include/astra/CompositeGeometryManager.h index c0acf4f..08eb7af 100644 --- a/include/astra/CompositeGeometryManager.h +++ b/include/astra/CompositeGeometryManager.h @@ -42,9 +42,9 @@ namespace astra { class CCompositeVolume; class CCompositeProjections; -class CFloat32Data3DMemory; -class CFloat32ProjectionData3DMemory; -class CFloat32VolumeData3DMemory; +class CFloat32Data3D; +class CFloat32ProjectionData3D; +class CFloat32VolumeData3D; class CVolumeGeometry3D; class CProjectionGeometry3D; class CProjector3D; @@ -77,7 +77,7 @@ public: PART_VOL, PART_PROJ } eType; - CFloat32Data3DMemory* pData; + CFloat32Data3D* pData; unsigned int subX; unsigned int subY; unsigned int subZ; @@ -88,8 +88,11 @@ public: virtual void splitY(TPartList& out, size_t maxSize, size_t maxDim, int div) = 0; virtual void splitZ(TPartList& out, size_t maxSize, size_t maxDim, int div) = 0; virtual CPart* reduce(const CPart *other) = 0; - virtual void getDims(size_t &x, size_t &y, size_t &z) = 0; - size_t getSize(); + virtual void getDims(size_t &x, size_t &y, size_t &z) const = 0; + size_t getSize() const; + + bool canSplitAndReduce() const; + bool isFull() const; }; class CVolumePart : public CPart { @@ -104,7 +107,7 @@ public: virtual void splitY(TPartList& out, size_t maxSize, size_t maxDim, int div); virtual void splitZ(TPartList& out, size_t maxSize, size_t maxDim, int div); virtual CPart* reduce(const CPart *other); - virtual void getDims(size_t &x, size_t &y, size_t &z); + virtual void getDims(size_t &x, size_t &y, size_t &z) const; CVolumePart* clone() const; }; @@ -120,7 +123,7 @@ public: virtual void splitY(TPartList& out, size_t maxSize, size_t maxDim, int div); virtual void splitZ(TPartList& out, size_t maxSize, size_t maxDim, int div); virtual CPart* reduce(const CPart *other); - virtual void getDims(size_t &x, size_t &y, size_t &z); + virtual void getDims(size_t &x, size_t &y, size_t &z) const; CProjectionPart* clone() const; }; @@ -150,23 +153,23 @@ public: bool doJobs(TJobList &jobs); SJob createJobFP(CProjector3D *pProjector, - CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData); + CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData); SJob createJobBP(CProjector3D *pProjector, - CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData); + CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData); // Convenience functions for creating and running a single FP or BP job - bool doFP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData); - bool doBP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData); - bool doFDK(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData, bool bShortScan, + bool doFP(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData); + bool doBP(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData); + bool doFDK(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData, bool bShortScan, const float *pfFilter = 0); - bool doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData); - bool doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData); + bool doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3D *>& volData, const std::vector<CFloat32ProjectionData3D *>& projData); + bool doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3D *>& volData, const std::vector<CFloat32ProjectionData3D *>& projData); void setGPUIndices(const std::vector<int>& GPUIndices); diff --git a/include/astra/CudaBackProjectionAlgorithm3D.h b/include/astra/CudaBackProjectionAlgorithm3D.h index 6738988..114d6f3 100644 --- a/include/astra/CudaBackProjectionAlgorithm3D.h +++ b/include/astra/CudaBackProjectionAlgorithm3D.h @@ -69,8 +69,8 @@ public: * @param _pReconstruction VolumeData3D object for storing the reconstructed volume. */ CCudaBackProjectionAlgorithm3D(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pProjectionData, - CFloat32VolumeData3DMemory* _pReconstruction); + CFloat32ProjectionData3D* _pProjectionData, + CFloat32VolumeData3D* _pReconstruction); /** Copy constructor. */ @@ -99,8 +99,8 @@ public: * @return initialization successful? */ bool initialize(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pSinogram, - CFloat32VolumeData3DMemory* _pReconstruction); + CFloat32ProjectionData3D* _pSinogram, + CFloat32VolumeData3D* _pReconstruction); /** Get all information parameters * diff --git a/include/astra/CudaFDKAlgorithm3D.h b/include/astra/CudaFDKAlgorithm3D.h index 386129e..1c4c622 100644 --- a/include/astra/CudaFDKAlgorithm3D.h +++ b/include/astra/CudaFDKAlgorithm3D.h @@ -81,8 +81,8 @@ public: * @param _pReconstruction VolumeData3D object for storing the reconstructed volume. */ CCudaFDKAlgorithm3D(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pProjectionData, - CFloat32VolumeData3DMemory* _pReconstruction); + CFloat32ProjectionData3D* _pProjectionData, + CFloat32VolumeData3D* _pReconstruction); /** Copy constructor. */ @@ -111,8 +111,8 @@ public: * @return initialization successful? */ bool initialize(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pSinogram, - CFloat32VolumeData3DMemory* _pReconstruction); + CFloat32ProjectionData3D* _pSinogram, + CFloat32VolumeData3D* _pReconstruction); /** Get all information parameters * diff --git a/include/astra/CudaForwardProjectionAlgorithm3D.h b/include/astra/CudaForwardProjectionAlgorithm3D.h index 95af73a..9dc889e 100644 --- a/include/astra/CudaForwardProjectionAlgorithm3D.h +++ b/include/astra/CudaForwardProjectionAlgorithm3D.h @@ -71,8 +71,8 @@ public: * @return initialization successful? */ bool initialize(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pSinogram, - CFloat32VolumeData3DMemory* _pReconstruction, + CFloat32ProjectionData3D* _pSinogram, + CFloat32VolumeData3D* _pReconstruction, int _iGPUindex = -1, int _iDetectorSuperSampling = 1); @@ -116,8 +116,8 @@ public: protected: CProjector3D* m_pProjector; - CFloat32ProjectionData3DMemory* m_pProjections; - CFloat32VolumeData3DMemory* m_pVolume; + CFloat32ProjectionData3D* m_pProjections; + CFloat32VolumeData3D* m_pVolume; int m_iGPUIndex; int m_iDetectorSuperSampling; diff --git a/include/astra/Float32Data2D.h b/include/astra/Float32Data2D.h index b33bcd0..1524683 100644 --- a/include/astra/Float32Data2D.h +++ b/include/astra/Float32Data2D.h @@ -41,7 +41,7 @@ public: }; /** - * This class represents a two-dimensional block of float32ing point data. + * This class represents a 2-dimensional block of 32-bit floating point data. * It contains member functions for accessing this data and for performing * elementary computations on the data. * The data block is "owned" by the class, meaning that the class is diff --git a/include/astra/Float32Data3D.h b/include/astra/Float32Data3D.h index 1cd9c20..226689e 100644 --- a/include/astra/Float32Data3D.h +++ b/include/astra/Float32Data3D.h @@ -32,10 +32,12 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. #include "Float32Data.h" #include "Float32Data2D.h" +#include "../../cuda/3d/mem3d.h" + namespace astra { /** - * This class represents a three-dimensional block of float32ing point data. + * This class represents a 3-dimensional block of 32-bit floating point data. */ class _AstraExport CFloat32Data3D : public CFloat32Data { @@ -107,22 +109,6 @@ public: */ int getDimensionCount() const; - /** - * Clamp data to minimum value - * - * @param _fMin minimum value - * @return l-value - */ - virtual CFloat32Data3D& clampMin(float32& _fMin) = 0; - - /** - * Clamp data to maximum value - * - * @param _fMax maximum value - * @return l-value - */ - virtual CFloat32Data3D& clampMax(float32& _fMax) = 0; - /** get a description of the class * * @return description string diff --git a/include/astra/Float32Data3DGPU.h b/include/astra/Float32Data3DGPU.h new file mode 100644 index 0000000..1247e65 --- /dev/null +++ b/include/astra/Float32Data3DGPU.h @@ -0,0 +1,107 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp + 2014-2016, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://www.astra-toolbox.com/ + +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/>. + +----------------------------------------------------------------------- +*/ + +#ifndef _INC_ASTRA_FLOAT32DATA3DGPU +#define _INC_ASTRA_FLOAT32DATA3DGPU + +#ifdef ASTRA_CUDA + +#include "Globals.h" +#include "Float32Data3D.h" + +#include "../../cuda/3d/mem3d.h" + +namespace astra { + + +astraCUDA3d::MemHandle3D wrapHandle(float *D_ptr, unsigned int x, unsigned int y, unsigned int z, unsigned int pitch); + + +/** + * This class represents a 3-dimensional block of 32-bit floating point data. + * The data block is stored on a GPU, and owned by external code. + * + * TODO: Store/remember which GPU the data is stored on + */ +class _AstraExport CFloat32Data3DGPU : public virtual CFloat32Data3D { + +protected: + /** Handle for the memory block */ + astraCUDA3d::MemHandle3D m_hnd; + + /** Clear all member variables, setting all numeric variables to 0 and all pointers to NULL. + */ + void _clear(); + + /** Un-initialize the object, bringing it back in the unitialized state. + */ + void _unInit(); + + /** Initialization. Initializes an instance of the CFloat32Data3DGPU class. + * Can only be called by derived classes. + * + * This function does not set m_bInitialized to true if everything is ok. + * + * @param _iWidth width of the 3D data (x-axis), must be > 0 + * @param _iHeight height of the 3D data (y-axis), must be > 0 + * @param _iDepth depth of the 3D data (z-axis), must be > 0 + * @param _hnd the CUDA memory handle + */ + + bool _initialize(int _iWidth, int _iHeight, int _iDepth, astraCUDA3d::MemHandle3D _hnd); + +public: + + /** Default constructor. Sets all numeric member variables to 0 and all pointer member variables to NULL. + * + * If an object is constructed using this default constructor, it must always be followed by a call + * to one of the initialize() methods before the object can be used. Any use before calling init() is not allowed, + * except calling the member function isInitialized(). + * + */ + CFloat32Data3DGPU(); + + /** Destructor. + */ + virtual ~CFloat32Data3DGPU(); + + /** which type is this class? + * + * @return DataType: ASTRA_DATATYPE_FLOAT32_PROJECTION or + * ASTRA_DATATYPE_FLOAT32_VOLUME + */ + virtual EDataType getType() const { return BASE; } + + astraCUDA3d::MemHandle3D getHandle() const { return m_hnd; } + +}; + +} // end namespace astra + +#endif + +#endif // _INC_ASTRA_FLOAT32DATA3DGPU diff --git a/include/astra/Float32Data3DMemory.h b/include/astra/Float32Data3DMemory.h index d885101..876aa37 100644 --- a/include/astra/Float32Data3DMemory.h +++ b/include/astra/Float32Data3DMemory.h @@ -34,7 +34,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. namespace astra { /** - * This class represents a three-dimensional block of float32ing point data. + * This class represents a 3-dimensional block of 32-bit floating point data. * It contains member functions for accessing this data and for performing * elementary computations on the data. * The data block is "owned" by the class, meaning that the class is @@ -106,8 +106,8 @@ protected: * This function does not set m_bInitialized to true if everything is ok. * * @param _iWidth width of the 3D data (x-axis), must be > 0 - * @param _iHeight height of the 3D data (y-axis), must be > 0 - * @param _iDepth depth of the 3D data (z-axis), must be > 0 + * @param _iHeight height of the 3D data (y-axis), must be > 0 + * @param _iDepth depth of the 3D data (z-axis), must be > 0 * @return initialization of the base class successfull */ bool _initialize(int _iWidth, int _iHeight, int _iDepth); @@ -121,9 +121,9 @@ protected: * object is reinitialized and memory is freed and reallocated if necessary. * This function does not set m_bInitialized to true if everything is ok. * - * @param _iWidth width of the 2D data (x-axis), must be > 0 - * @param _iHeight height of the 2D data (y-axis), must be > 0 - * @param _iDepth depth of the 2D data (z-axis), must be > 0 + * @param _iWidth width of the 3D data (x-axis), must be > 0 + * @param _iHeight height of the 3D data (y-axis), must be > 0 + * @param _iDepth depth of the 3D data (z-axis), must be > 0 * @param _pfData pointer to a one-dimensional float32 data block * @return initialization of the base class successfull */ @@ -138,9 +138,9 @@ protected: * object is reinitialized and memory is freed and reallocated if necessary. * This function does not set m_bInitialized to true if everything is ok. * - * @param _iWidth width of the 2D data (x-axis), must be > 0 - * @param _iHeight height of the 2D data (y-axis), must be > 0 - * @param _iDepth depth of the 2D data (z-axis), must be > 0 + * @param _iWidth width of the 3D data (x-axis), must be > 0 + * @param _iHeight height of the 3D data (y-axis), must be > 0 + * @param _iDepth depth of the 3D data (z-axis), must be > 0 * @param _fScalar scalar value to fill the data * @return initialization of the base class successfull */ @@ -157,9 +157,9 @@ protected: * object is reinitialized and memory is freed and reallocated if necessary. * This function does not set m_bInitialized to true if everything is ok. * - * @param _iWidth width of the 2D data (x-axis), must be > 0 - * @param _iHeight height of the 2D data (y-axis), must be > 0 - * @param _iDepth depth of the 2D data (z-axis), must be > 0 + * @param _iWidth width of the 3D data (x-axis), must be > 0 + * @param _iHeight height of the 3D data (y-axis), must be > 0 + * @param _iDepth depth of the 3D data (z-axis), must be > 0 * @param _pCustomMemory the custom memory handle */ @@ -237,35 +237,11 @@ public: * After the call p = getData3D(), use p[iy][ix] to access element (ix, iy, iz). * The data memory and pointer array are still "owned" by the CFloat32Data3DMemory * instance; this memory may NEVER be freed by the caller of this function. - * If changes are made to this data, the function updateStatistics() - * should be called after completion of all changes. * * @return pointer to the 3-dimensional 32-bit floating point data block */ const float32*** getData3DConst() const; - /** Update data statistics, such as minimum and maximum value, after the data has been modified. - */ - virtual void updateStatistics(); - - /** Get the minimum value in the data block. - * If the data has been changed after construction, the function - * updateStatistics() must be called at least once before - * a query can be made on this value. - * - * @return minimum value in the data block - */ - virtual float32 getGlobalMin() const; - - /** Get the maximum value in the data block - * If the data has been changed after construction, the function - * updateStatistics() must be called at least once before - * a query can be made on this value. - * - * @return maximum value in the data block - */ - virtual float32 getGlobalMax() const; - /** which type is this class? * * @return DataType: ASTRA_DATATYPE_FLOAT32_PROJECTION or @@ -306,26 +282,11 @@ inline CFloat32Data3DMemory::EDataType CFloat32Data3DMemory::getType() const } //---------------------------------------------------------------------------------------- -// Get the minimum value in the data block. -inline float32 CFloat32Data3DMemory::getGlobalMin() const -{ - ASTRA_ASSERT(m_bInitialized); - return m_fGlobalMin; -} - -//---------------------------------------------------------------------------------------- -// Get the maximum value in the data block -inline float32 CFloat32Data3DMemory::getGlobalMax() const -{ - ASTRA_ASSERT(m_bInitialized); - return m_fGlobalMax; -} - -//---------------------------------------------------------------------------------------- // Get a pointer to the data block, represented as a 1-dimensional array of float32 values. inline float32* CFloat32Data3DMemory::getData() { ASTRA_ASSERT(m_bInitialized); + return m_pfData; } @@ -334,6 +295,7 @@ inline float32* CFloat32Data3DMemory::getData() inline const float32* CFloat32Data3DMemory::getDataConst() const { ASTRA_ASSERT(m_bInitialized); + return (const float32*)m_pfData; } diff --git a/include/astra/Float32ProjectionData3D.h b/include/astra/Float32ProjectionData3D.h index ae0664b..1634eeb 100644 --- a/include/astra/Float32ProjectionData3D.h +++ b/include/astra/Float32ProjectionData3D.h @@ -85,55 +85,6 @@ public: */ virtual CFloat32Data3D::EDataType getType() const; - /** Fetch a COPY of a projection of the data. Note that if you update the 2D data slice, the data in the - * 3d data object will remain unaltered. To copy the data back in the 3D-volume you must return the data by calling 'returnProjection'. - * - * @param _iProjectionNr projection number - * @return Volume data object - */ - virtual CFloat32VolumeData2D* fetchProjection(int _iProjectionNr) const = 0; - - /** Return a projection slice to the 3d data. The data will be deleted. If the slice was fetched with - * 'fetchProjection', the data will be stored first. - * - * @param _iProjectionNr projection number - * @param _pProjection 2D Projection Data - */ - virtual void returnProjection(int _iProjectionNr, CFloat32VolumeData2D* _pProjection) = 0; - - /** Fetch a COPY of a sinogram slice of the data. Note that if you update the 2D data slice, the data in the - * 3d data object will remain unaltered. To copy the data back in the 3D-volume you must return the data by calling 'returnSlice'. - * - * @param _iSliceNr slice number - * @return Sinogram data object - */ - virtual CFloat32ProjectionData2D* fetchSinogram(int _iSliceNr) const = 0; - - /** Return a sinogram slice to the 3d data. The data will be stored in the 3D Data object. - * - * @param _iSliceNr slice number - * @param _pSinogram2D 2D Sinogram Object. - */ - virtual void returnSinogram(int _iSliceNr, CFloat32ProjectionData2D* _pSinogram2D) = 0; - - /** This SLOW function returns a detector value stored a specific index in the array. - * Reading values in this way might cause a lot of unnecessar__y memory operations, don't - * use it in time-critical code. - * - * @param _iIndex Index in the array if the data were stored completely in main memory - * @return The value the location specified by _iIndex - */ - virtual float32 getDetectorValue(int _iIndex) = 0; - - /** This SLOW function stores a detector value at a specific index in the array. - * Writing values in this way might cause a lot of unnecessary memory operations, don't - * use it in time-critical code. - * - * @param _iIndex Index in the array if the data were stored completely in main memory - * @param _fValue The value to be stored at the location specified by _iIndex - */ - virtual void setDetectorValue(int _iIndex, float32 _fValue) = 0; - /** * Overloaded Operator: data += data (pointwise) * diff --git a/include/astra/Float32ProjectionData3DGPU.h b/include/astra/Float32ProjectionData3DGPU.h new file mode 100644 index 0000000..cbf9b6f --- /dev/null +++ b/include/astra/Float32ProjectionData3DGPU.h @@ -0,0 +1,91 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp + 2014-2016, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://www.astra-toolbox.com/ + +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/>. + +----------------------------------------------------------------------- +*/ + +#ifndef _INC_ASTRA_FLOAT32PROJECTIONDATA3DGPU +#define _INC_ASTRA_FLOAT32PROJECTIONDATA3DGPU + +#include "Float32Data3DGPU.h" +#include "ProjectionGeometry3D.h" +#include "Float32ProjectionData3D.h" + +#ifdef ASTRA_CUDA + +namespace astra { + +/** + * This class represents three-dimensional Projection Data where the entire data block is stored in GPU memory. + */ +class _AstraExport CFloat32ProjectionData3DGPU : public CFloat32Data3DGPU, public CFloat32ProjectionData3D +{ +public: + + /** Default constructor. Sets all numeric member variables to 0 and all pointer member variables to NULL. + * + * If an object is constructed using this default constructor, it must always be followed by a call + * to one of the init() methods before the object can be used. Any use before calling init() is not allowed, + * except calling the member function isInitialized(). + * + */ + CFloat32ProjectionData3DGPU(); + + /** Construction. + * + * @param _pGeometry 3D volume geometry + * @param _hnd the CUDA memory handle + */ + + CFloat32ProjectionData3DGPU(CProjectionGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd); + + virtual ~CFloat32ProjectionData3DGPU(); + + /** Initialization. + * + * @param _pGeometry 3D volume geometry + * @param _hnd the CUDA memory handle + */ + + bool initialize(CProjectionGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd); + + /** Which type is this class? + * + * @return DataType: PROJECTION + */ + virtual CFloat32Data3D::EDataType getType() const { return PROJECTION; } + + /** Get the volume geometry. + * + * @return pointer to volume geometry. + */ + CProjectionGeometry3D* getGeometry() const { ASTRA_ASSERT(m_bInitialized); return m_pGeometry; } + +}; + +} // end namesProjection astra + +#endif + +#endif // _INC_ASTRA_FLOAT32PROJECTIONDATA3DGPU diff --git a/include/astra/Float32ProjectionData3DMemory.h b/include/astra/Float32ProjectionData3DMemory.h index 14bae36..732c31a 100644 --- a/include/astra/Float32ProjectionData3DMemory.h +++ b/include/astra/Float32ProjectionData3DMemory.h @@ -166,55 +166,6 @@ public: - /** Fetch a COPY of a projection of the data. Note that if you update the 2D data slice, the data in the - * 3D data object will remain unaltered. To copy the data back in the 3D-volume you must return the data by calling 'returnProjection'. - * - * @param _iProjectionNr projection number - * @return Volume data object - */ - virtual CFloat32VolumeData2D* fetchProjection(int _iProjectionNr) const; - - /** Return a projection slice to the 3D data. The data will be deleted. If the slice was fetched with - * 'fetchProjection', the data will be stored first. - * - * @param _iProjectionNr projection number - * @param _pProjection 2D Projection image - */ - virtual void returnProjection(int _iProjectionNr, CFloat32VolumeData2D* _pProjection); - - /** Fetch a COPY of a sinogram slice of the data. Note that if you update the 2D data slice, the data in the - * 3D data object will remain unaltered. To copy the data back in the 3D-volume you must return the data by calling 'returnSlice'. - * - * @param _iSliceNr slice number - * @return Sinogram data object - */ - virtual CFloat32ProjectionData2D* fetchSinogram(int _iSliceNr) const; - - /** This SLOW function returns a detector value stored a specific index in the array. - * Reading values in this way might cause a lot of unnecessary memory operations, don't - * use it in time-critical code. - * - * @param _iIndex Index in the array if the data were stored completely in main memory - * @return The value the location specified by _iIndex - */ - virtual float32 getDetectorValue(int _iIndex); - - /** This SLOW function stores a detector value at a specific index in the array. - * Writing values in this way might cause a lot of unnecessary memory operations, don't - * use it in time-critical code. - * - * @param _iIndex Index in the array if the data were stored completely in main memory - * @param _fValue The value to be stored at the location specified by _iIndex - */ - virtual void setDetectorValue(int _iIndex, float32 _fValue); - - /** Return a sinogram slice to the 3d data. The data will be stored in the 3D Data object. - * - * @param _iSliceNr slice number - * @param _pSinogram2D 2D Sinogram Object. - */ - virtual void returnSinogram(int _iSliceNr, CFloat32ProjectionData2D* _pSinogram2D); - /** Which type is this class? * * @return DataType: PROJECTION diff --git a/include/astra/Float32VolumeData3D.h b/include/astra/Float32VolumeData3D.h index 7a32efe..535e960 100644 --- a/include/astra/Float32VolumeData3D.h +++ b/include/astra/Float32VolumeData3D.h @@ -85,72 +85,6 @@ public: */ virtual CFloat32Data3D::EDataType getType() const; - /** Fetch a slice from the data in the x direction. Note that if you update the 2D data slice, the data in the - * 3d data object will remain unaltered. To copy the data you must return the data by calling 'returnSliceX'. - * You should not delete data fetched with this function yourself, instead call the 'returnSliceX' function. - * - * @param _iColumnIndex slice number - * @return Volume data object - */ - virtual CFloat32VolumeData2D* fetchSliceX(int _iColumnIndex) const = 0; - - /** Fetch a slice from the data in the y direction. Note that if you update the 2D data slice, the data in the - * 3d data object will remain unaltered. To copy the data you must return the data by calling 'returnSliceY'. - * You should not delete data fetched with this function yourself, instead call the 'returnSliceY' function. - * - * @param _iRowIndex slice number - * @return Volume data object - */ - virtual CFloat32VolumeData2D* fetchSliceY(int _iRowIndex) const = 0; - - /** Fetch a slice from the data in the z direction. Note that if you update the 2D data slice, the data in the - * 3d data object will remain unaltered. To copy the data you must return the data by calling 'returnSliceZ'. - * You should not delete data fetched with this function yourself, instead call the 'returnSliceZ' function. - * - * @param _iSliceIndex slice number - * @return Volume data object - */ - virtual CFloat32VolumeData2D* fetchSliceZ(int _iSliceIndex) const = 0; - - /** Return a slice from the data in the x direction to the 3d data. The data will be deleted. If the slice was - * fetched with 'fetchSliceX', the data will be stored first. - * - * @param _iColumnIndex slice number - */ - virtual void returnSliceX(int _iColumnIndex, CFloat32VolumeData2D * _pSlice) = 0; - - /** Return a slice from the data in the y direction to the 3d data. The data will be deleted. If the slice was - * fetched with 'fetchSliceY', the data will be stored first. - * - * @param _iRowIndex slice number - */ - virtual void returnSliceY(int _iRowIndex, CFloat32VolumeData2D * _pSlice) = 0; - - /** Return a slice from the data in the z direction to the 3d data. The data will be deleted. If the slice was - * fetched with 'fetchSliceZ', the data will be stored first. - * - * @param _iSliceIndex slice number - */ - virtual void returnSliceZ(int _iSliceIndex, CFloat32VolumeData2D * _pSlice) = 0; - - /** This SLOW function returns a voxel value stored at a specific index in the array. - * Reading values in this way might cause a lot of unnecessary memory operations, don't - * use it in time-critical code. - * - * @param _iIndex Index in the array if the data were stored completely in main memory - * @return The value stored at the location specified by _iIndex - */ - virtual float32 getVoxelValue(int _iIndex) = 0; - - /** This SLOW function stores a voxel value at a specific index in the array. - * Writing values in this way might cause a lot of unnecessary memory operations, don't - * use it in time-critical code. - * - * @param _iIndex Index in the array if the data were stored completely in main memory - * @param _fValue The value to be stored at the location specified by _iIndex - */ - virtual void setVoxelValue(int _iIndex, float32 _fValue) = 0; - /** * Overloaded Operator: data += data (pointwise) * diff --git a/include/astra/Float32VolumeData3DGPU.h b/include/astra/Float32VolumeData3DGPU.h new file mode 100644 index 0000000..ddeb469 --- /dev/null +++ b/include/astra/Float32VolumeData3DGPU.h @@ -0,0 +1,91 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp + 2014-2016, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://www.astra-toolbox.com/ + +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/>. + +----------------------------------------------------------------------- +*/ + +#ifndef _INC_ASTRA_FLOAT32VOLUMEDATA3DGPU +#define _INC_ASTRA_FLOAT32VOLUMEDATA3DGPU + +#include "Float32Data3DGPU.h" +#include "VolumeGeometry3D.h" +#include "Float32VolumeData3D.h" + +#ifdef ASTRA_CUDA + +namespace astra { + +/** + * This class represents three-dimensional Volume Data where the entire data block is stored in GPU memory. + */ +class _AstraExport CFloat32VolumeData3DGPU : public CFloat32Data3DGPU, public CFloat32VolumeData3D +{ +public: + + /** Default constructor. Sets all numeric member variables to 0 and all pointer member variables to NULL. + * + * If an object is constructed using this default constructor, it must always be followed by a call + * to one of the init() methods before the object can be used. Any use before calling init() is not allowed, + * except calling the member function isInitialized(). + * + */ + CFloat32VolumeData3DGPU(); + + /** Construction. + * + * @param _pGeometry 3D volume geometry + * @param _hnd the CUDA memory handle + */ + + CFloat32VolumeData3DGPU(CVolumeGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd); + + virtual ~CFloat32VolumeData3DGPU(); + + /** Initialization. + * + * @param _pGeometry 3D volume geometry + * @param _hnd the CUDA memory handle + */ + + bool initialize(CVolumeGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd); + + /** Which type is this class? + * + * @return DataType: VOLUME + */ + virtual CFloat32Data3D::EDataType getType() const { return VOLUME; } + + /** Get the volume geometry. + * + * @return pointer to volume geometry. + */ + CVolumeGeometry3D* getGeometry() const { ASTRA_ASSERT(m_bInitialized); return m_pGeometry; } + +}; + +} // end namespace astra + +#endif + +#endif // _INC_ASTRA_FLOAT32VOLUMEDATA3DGPU diff --git a/include/astra/Float32VolumeData3DMemory.h b/include/astra/Float32VolumeData3DMemory.h index 70f0939..4a70f65 100644 --- a/include/astra/Float32VolumeData3DMemory.h +++ b/include/astra/Float32VolumeData3DMemory.h @@ -159,55 +159,6 @@ public: CVolumeGeometry3D* getGeometry() const; /** - * Gets a slice, containing all voxels with a given x (= column) index. - */ - CFloat32VolumeData2D * fetchSliceX(int _iColumnIndex) const; - - /** - * Gets a slice, containing all voxels with a given y (= row) index. - */ - CFloat32VolumeData2D * fetchSliceY(int _iRowIndex) const; - - /** - * Gets a slice, containing all voxels with a given z (= slice) index. - */ - CFloat32VolumeData2D * fetchSliceZ(int _iSliceIndex) const; - - /** - * Gets a slice, containing all voxels with a given x (= column) index. - */ - void returnSliceX(int _iColumnIndex, CFloat32VolumeData2D * _pSliceData); - - /** - * Gets a slice, containing all voxels with a given y (= row) index. - */ - void returnSliceY(int _iRowIndex, CFloat32VolumeData2D * _pSliceData); - - /** - * Copies data from a 2D slice containing all voxels with a given z (= slice) index to the - * 3D memory stored in this class. - */ - void returnSliceZ(int _iSliceIndex, CFloat32VolumeData2D * _pSliceData); - - /** This SLOW function returns a volume value stored a specific index in the array. - * Reading values in this way might cause a lot of unnecessary memory operations, don't - * use it in time-critical code. - * - * @param _iIndex Index in the array if the data were stored completely in main memory - * @return The value the location specified by _iIndex - */ - virtual float32 getVoxelValue(int _iIndex); - - /** This SLOW function stores a voxel value at a specific index in the array. - * Writing values in this way might cause a lot of unnecessary memory operations, don't - * use it in time-critical code. - * - * @param _iIndex Index in the array if the data were stored completely in main memory - * @param _fValue The value to be stored at the location specified by _iIndex - */ - virtual void setVoxelValue(int _iIndex, float32 _fValue); - - /** * Overloaded Operator: data = data (pointwise) * * @param _dataIn r-value diff --git a/include/astra/Globals.h b/include/astra/Globals.h index dec978d..44a77b0 100644 --- a/include/astra/Globals.h +++ b/include/astra/Globals.h @@ -53,8 +53,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. #include <iostream> #include <fstream> #include <math.h> -#include <boost/static_assert.hpp> -#include <boost/throw_exception.hpp> //---------------------------------------------------------------------------------------- // macro's @@ -227,16 +225,6 @@ namespace astra { int m_iSliceIndex; }; } -//---------------------------------------------------------------------------------------- -// some toys - -// safe reinterpret cast -template <class To, class From> -To safe_reinterpret_cast(From from) -{ - BOOST_STATIC_ASSERT(sizeof(From) <= sizeof(To)); - return reinterpret_cast<To>(from); -} //---------------------------------------------------------------------------------------- // functions for testing diff --git a/matlab/mex/astra_mex_data3d_c.cpp b/matlab/mex/astra_mex_data3d_c.cpp index 7909a79..f0aa3c8 100644 --- a/matlab/mex/astra_mex_data3d_c.cpp +++ b/matlab/mex/astra_mex_data3d_c.cpp @@ -60,7 +60,7 @@ using namespace astra; //----------------------------------------------------------------------------------------- /** * id = astra_mex_io_data('create', datatype, geometry, data); - * datatype: ['-vol','-sino','-sinocone'] + * datatype: ['-vol','-sino'] */ void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArray* prhs[]) { @@ -103,7 +103,6 @@ void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra copyMexToCFloat32Array(data, pDataObject3D->getData(), pDataObject3D->getSize()); } - pDataObject3D->updateStatistics(); // step4: store data object int iIndex = CData3DManager::getSingleton().store(pDataObject3D); @@ -162,8 +161,6 @@ void astra_mex_data3d_link(int& nlhs, mxArray* plhs[], int& nrhs, const mxArray* return; } - //pDataObject3D->updateStatistics(); - // step4: store data object int iIndex = CData3DManager::getSingleton().store(pDataObject3D); @@ -234,7 +231,6 @@ void astra_mex_data3d_store(int nlhs, mxArray* plhs[], int nrhs, const mxArray* } copyMexToCFloat32Array(prhs[2], pDataObject->getData(), pDataObject->getSize()); - pDataObject->updateStatistics(); } void astra_mex_data3d_dimensions(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) diff --git a/matlab/mex/mexDataManagerHelpFunctions.cpp b/matlab/mex/mexDataManagerHelpFunctions.cpp index ebf6f8f..c47ee27 100644 --- a/matlab/mex/mexDataManagerHelpFunctions.cpp +++ b/matlab/mex/mexDataManagerHelpFunctions.cpp @@ -166,17 +166,6 @@ checkDataSize(const mxArray * const mArray, //----------------------------------------------------------------------------------------- void -updateStatistics(const std::vector<astra::CFloat32Data3DMemory *> & vecIn) -{ - const size_t tot_size = vecIn.size(); - for (size_t count = 0; count < tot_size; count++) - { - vecIn[count]->updateStatistics(); - } -} - -//----------------------------------------------------------------------------------------- -void getDataPointers(const std::vector<astra::CFloat32Data3DMemory *> & vecIn, std::vector<astra::float32 *> & vecOut) { diff --git a/matlab/mex/mexDataManagerHelpFunctions.h b/matlab/mex/mexDataManagerHelpFunctions.h index fff83a0..b7edb0f 100644 --- a/matlab/mex/mexDataManagerHelpFunctions.h +++ b/matlab/mex/mexDataManagerHelpFunctions.h @@ -50,8 +50,6 @@ bool checkDataSize(const mxArray * const, const astra::CProjectionGeometry3D * c bool checkDataSize(const mxArray * const, const astra::CVolumeGeometry3D * const, const mwIndex & zOffset); -void updateStatistics(const std::vector<astra::CFloat32Data3DMemory *> &); - void getDataPointers(const std::vector<astra::CFloat32Data3DMemory *> &, std::vector<astra::float32 *> &); void getDataSizes(const std::vector<astra::CFloat32Data3DMemory *> &, diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd index 39f9039..b40d787 100644 --- a/python/astra/PyIncludes.pxd +++ b/python/astra/PyIncludes.pxd @@ -30,242 +30,262 @@ from .PyXMLDocument cimport XMLNode include "config.pxi" cdef extern from "astra/Globals.h" namespace "astra": - ctypedef float float32 - ctypedef double float64 - ctypedef unsigned short int uint16 - ctypedef signed short int sint16 - ctypedef unsigned char uchar8 - ctypedef signed char schar8 - ctypedef int int32 - ctypedef short int int16 + ctypedef float float32 + ctypedef double float64 + ctypedef unsigned short int uint16 + ctypedef signed short int sint16 + ctypedef unsigned char uchar8 + ctypedef signed char schar8 + ctypedef int int32 + ctypedef short int int16 cdef extern from "astra/Config.h" namespace "astra": - cdef cppclass Config: - Config() - void initialize(string rootname) - XMLNode self + cdef cppclass Config: + Config() + void initialize(string rootname) + XMLNode self cdef extern from "astra/VolumeGeometry2D.h" namespace "astra": - cdef cppclass CVolumeGeometry2D: - bool initialize(Config) - int getGridColCount() - int getGridRowCount() - int getGridTotCount() - float32 getWindowLengthX() - float32 getWindowLengthY() - float32 getWindowArea() - float32 getPixelLengthX() - float32 getPixelLengthY() - float32 getPixelArea() - float32 getWindowMinX() - float32 getWindowMinY() - float32 getWindowMaxX() - float32 getWindowMaxY() - Config* getConfiguration() - bool isEqual(CVolumeGeometry2D*) + cdef cppclass CVolumeGeometry2D: + bool initialize(Config) + int getGridColCount() + int getGridRowCount() + int getGridTotCount() + float32 getWindowLengthX() + float32 getWindowLengthY() + float32 getWindowArea() + float32 getPixelLengthX() + float32 getPixelLengthY() + float32 getPixelArea() + float32 getWindowMinX() + float32 getWindowMinY() + float32 getWindowMaxX() + float32 getWindowMaxY() + Config* getConfiguration() + bool isEqual(CVolumeGeometry2D*) cdef extern from "astra/Float32Data2D.h" namespace "astra": - cdef cppclass CFloat32CustomMemory: - pass + cdef cppclass CFloat32CustomMemory: + pass cdef extern from "astra/Float32VolumeData2D.h" namespace "astra": - cdef cppclass CFloat32VolumeData2D: - CFloat32VolumeData2D(CVolumeGeometry2D*) - CFloat32VolumeData2D(CVolumeGeometry2D*, CFloat32CustomMemory*) - CVolumeGeometry2D * getGeometry() - int getWidth() - int getHeight() - void changeGeometry(CVolumeGeometry2D*) - Config* getConfiguration() + cdef cppclass CFloat32VolumeData2D: + CFloat32VolumeData2D(CVolumeGeometry2D*) + CFloat32VolumeData2D(CVolumeGeometry2D*, CFloat32CustomMemory*) + CVolumeGeometry2D * getGeometry() + int getWidth() + int getHeight() + void changeGeometry(CVolumeGeometry2D*) + Config* getConfiguration() cdef extern from "astra/ProjectionGeometry2D.h" namespace "astra": - cdef cppclass CProjectionGeometry2D: - CProjectionGeometry2D() - bool initialize(Config) - int getDetectorCount() - int getProjectionAngleCount() - bool isOfType(string) - float32 getProjectionAngle(int) - float32 getDetectorWidth() - Config* getConfiguration() - bool isEqual(CProjectionGeometry2D*) + cdef cppclass CProjectionGeometry2D: + CProjectionGeometry2D() + bool initialize(Config) + int getDetectorCount() + int getProjectionAngleCount() + bool isOfType(string) + float32 getProjectionAngle(int) + float32 getDetectorWidth() + Config* getConfiguration() + bool isEqual(CProjectionGeometry2D*) cdef extern from "astra/Float32Data2D.h" namespace "astra::CFloat32Data2D": - cdef enum TWOEDataType "astra::CFloat32Data2D::EDataType": - TWOPROJECTION "astra::CFloat32Data2D::PROJECTION" - TWOVOLUME "astra::CFloat32Data2D::VOLUME" + cdef enum TWOEDataType "astra::CFloat32Data2D::EDataType": + TWOPROJECTION "astra::CFloat32Data2D::PROJECTION" + TWOVOLUME "astra::CFloat32Data2D::VOLUME" cdef extern from "astra/Float32Data3D.h" namespace "astra::CFloat32Data3D": - cdef enum THREEEDataType "astra::CFloat32Data3D::EDataType": - THREEPROJECTION "astra::CFloat32Data3D::PROJECTION" - THREEVOLUME "astra::CFloat32Data3D::VOLUME" + cdef enum THREEEDataType "astra::CFloat32Data3D::EDataType": + THREEPROJECTION "astra::CFloat32Data3D::PROJECTION" + THREEVOLUME "astra::CFloat32Data3D::VOLUME" cdef extern from "astra/Float32Data2D.h" namespace "astra": - cdef cppclass CFloat32Data2D: - bool isInitialized() - int getSize() - float32 *getData() - float32 **getData2D() - int getWidth() - int getHeight() - TWOEDataType getType() + cdef cppclass CFloat32Data2D: + bool isInitialized() + int getSize() + float32 *getData() + float32 **getData2D() + int getWidth() + int getHeight() + TWOEDataType getType() cdef extern from "astra/SparseMatrixProjectionGeometry2D.h" namespace "astra": - cdef cppclass CSparseMatrixProjectionGeometry2D: - CSparseMatrixProjectionGeometry2D() + cdef cppclass CSparseMatrixProjectionGeometry2D: + CSparseMatrixProjectionGeometry2D() cdef extern from "astra/FanFlatProjectionGeometry2D.h" namespace "astra": - cdef cppclass CFanFlatProjectionGeometry2D: - CFanFlatProjectionGeometry2D() + cdef cppclass CFanFlatProjectionGeometry2D: + CFanFlatProjectionGeometry2D() cdef extern from "astra/FanFlatVecProjectionGeometry2D.h" namespace "astra": - cdef cppclass CFanFlatVecProjectionGeometry2D: - CFanFlatVecProjectionGeometry2D() + cdef cppclass CFanFlatVecProjectionGeometry2D: + CFanFlatVecProjectionGeometry2D() cdef extern from "astra/ParallelProjectionGeometry2D.h" namespace "astra": - cdef cppclass CParallelProjectionGeometry2D: - CParallelProjectionGeometry2D() + cdef cppclass CParallelProjectionGeometry2D: + CParallelProjectionGeometry2D() cdef extern from "astra/Float32ProjectionData2D.h" namespace "astra": - cdef cppclass CFloat32ProjectionData2D: - CFloat32ProjectionData2D(CProjectionGeometry2D*) - CFloat32ProjectionData2D(CProjectionGeometry2D*, CFloat32CustomMemory*) - CProjectionGeometry2D * getGeometry() - void changeGeometry(CProjectionGeometry2D*) - int getDetectorCount() - int getAngleCount() + cdef cppclass CFloat32ProjectionData2D: + CFloat32ProjectionData2D(CProjectionGeometry2D*) + CFloat32ProjectionData2D(CProjectionGeometry2D*, CFloat32CustomMemory*) + CProjectionGeometry2D * getGeometry() + void changeGeometry(CProjectionGeometry2D*) + int getDetectorCount() + int getAngleCount() cdef extern from "astra/Algorithm.h" namespace "astra": - cdef cppclass CAlgorithm: - bool initialize(Config) - void run(int) nogil - bool isInitialized() + cdef cppclass CAlgorithm: + bool initialize(Config) + void run(int) nogil + bool isInitialized() cdef extern from "astra/ReconstructionAlgorithm2D.h" namespace "astra": - cdef cppclass CReconstructionAlgorithm2D: - bool getResidualNorm(float32&) + cdef cppclass CReconstructionAlgorithm2D: + bool getResidualNorm(float32&) cdef extern from "astra/Projector2D.h" namespace "astra": - cdef cppclass CProjector2D: - bool isInitialized() - CProjectionGeometry2D* getProjectionGeometry() - CVolumeGeometry2D* getVolumeGeometry() - CSparseMatrix* getMatrix() + cdef cppclass CProjector2D: + bool isInitialized() + CProjectionGeometry2D* getProjectionGeometry() + CVolumeGeometry2D* getVolumeGeometry() + CSparseMatrix* getMatrix() cdef extern from "astra/Projector3D.h" namespace "astra": - cdef cppclass CProjector3D: - bool isInitialized() - CProjectionGeometry3D* getProjectionGeometry() - CVolumeGeometry3D* getVolumeGeometry() + cdef cppclass CProjector3D: + bool isInitialized() + CProjectionGeometry3D* getProjectionGeometry() + CVolumeGeometry3D* getVolumeGeometry() IF HAVE_CUDA==True: - cdef extern from "astra/CudaProjector3D.h" namespace "astra": - cdef cppclass CCudaProjector3D + cdef extern from "astra/CudaProjector3D.h" namespace "astra": + cdef cppclass CCudaProjector3D - cdef extern from "astra/CudaProjector2D.h" namespace "astra": - cdef cppclass CCudaProjector2D + cdef extern from "astra/CudaProjector2D.h" namespace "astra": + cdef cppclass CCudaProjector2D + + cdef extern from "astra/Float32Data3DGPU.h" namespace "astraCUDA3d": + cdef cppclass MemHandle3D: + pass + + cdef extern from "astra/Float32Data3DGPU.h" namespace "astraCUDA3d": + cdef MemHandle3D wrapHandle(float *D_ptr, unsigned int x, unsigned int y, unsigned int z, unsigned int pitch) cdef extern from "astra/SparseMatrix.h" namespace "astra": - cdef cppclass CSparseMatrix: - CSparseMatrix(unsigned int,unsigned int,unsigned long) - unsigned int m_iWidth - unsigned int m_iHeight - unsigned long m_lSize - bool isInitialized() - float32* m_pfValues - unsigned int* m_piColIndices - unsigned long* m_plRowStarts + cdef cppclass CSparseMatrix: + CSparseMatrix(unsigned int,unsigned int,unsigned long) + unsigned int m_iWidth + unsigned int m_iHeight + unsigned long m_lSize + bool isInitialized() + float32* m_pfValues + unsigned int* m_piColIndices + unsigned long* m_plRowStarts + +cdef extern from "astra/Float32Data3D.h" namespace "astra": + cdef cppclass CFloat32Data3D: + bool isInitialized() + int getSize() + int getWidth() + int getHeight() + int getDepth() + cdef extern from "astra/Float32Data3DMemory.h" namespace "astra": - cdef cppclass CFloat32Data3DMemory: - CFloat32Data3DMemory() - bool isInitialized() - int getSize() - int getWidth() - int getHeight() - int getDepth() - void updateStatistics() - float32 *getData() - float32 ***getData3D() - THREEEDataType getType() + cdef cppclass CFloat32Data3DMemory(CFloat32Data3D): + CFloat32Data3DMemory() + void updateStatistics() + float32 *getData() + float32 ***getData3D() + THREEEDataType getType() cdef extern from "astra/VolumeGeometry3D.h" namespace "astra": - cdef cppclass CVolumeGeometry3D: - CVolumeGeometry3D() - bool initialize(Config) - Config * getConfiguration() - int getGridColCount() - int getGridRowCount() - int getGridSliceCount() + cdef cppclass CVolumeGeometry3D: + CVolumeGeometry3D() + bool initialize(Config) + Config * getConfiguration() + int getGridColCount() + int getGridRowCount() + int getGridSliceCount() cdef extern from "astra/ProjectionGeometry3D.h" namespace "astra": - cdef cppclass CProjectionGeometry3D: - CProjectionGeometry3D() - bool initialize(Config) - Config * getConfiguration() - int getProjectionCount() - int getDetectorColCount() - int getDetectorRowCount() + cdef cppclass CProjectionGeometry3D: + CProjectionGeometry3D() + bool initialize(Config) + Config * getConfiguration() + int getProjectionCount() + int getDetectorColCount() + int getDetectorRowCount() cdef extern from "astra/Float32VolumeData3DMemory.h" namespace "astra": - cdef cppclass CFloat32VolumeData3DMemory: - CFloat32VolumeData3DMemory(CVolumeGeometry3D*) - CFloat32VolumeData3DMemory(CVolumeGeometry3D*, CFloat32CustomMemory*) - CVolumeGeometry3D* getGeometry() - void changeGeometry(CVolumeGeometry3D*) - int getRowCount() - int getColCount() - int getSliceCount() - bool isInitialized() - - + cdef cppclass CFloat32VolumeData3DMemory: + CFloat32VolumeData3DMemory(CVolumeGeometry3D*) + CFloat32VolumeData3DMemory(CVolumeGeometry3D*, CFloat32CustomMemory*) + CVolumeGeometry3D* getGeometry() + void changeGeometry(CVolumeGeometry3D*) + int getRowCount() + int getColCount() + int getSliceCount() + bool isInitialized() cdef extern from "astra/ParallelProjectionGeometry3D.h" namespace "astra": - cdef cppclass CParallelProjectionGeometry3D: - CParallelProjectionGeometry3D() + cdef cppclass CParallelProjectionGeometry3D: + CParallelProjectionGeometry3D() cdef extern from "astra/ParallelVecProjectionGeometry3D.h" namespace "astra": - cdef cppclass CParallelVecProjectionGeometry3D: - CParallelVecProjectionGeometry3D() + cdef cppclass CParallelVecProjectionGeometry3D: + CParallelVecProjectionGeometry3D() cdef extern from "astra/ConeProjectionGeometry3D.h" namespace "astra": - cdef cppclass CConeProjectionGeometry3D: - CConeProjectionGeometry3D() - bool initialize(Config) + cdef cppclass CConeProjectionGeometry3D: + CConeProjectionGeometry3D() + bool initialize(Config) cdef extern from "astra/ConeVecProjectionGeometry3D.h" namespace "astra": - cdef cppclass CConeVecProjectionGeometry3D: - CConeVecProjectionGeometry3D() + cdef cppclass CConeVecProjectionGeometry3D: + CConeVecProjectionGeometry3D() cdef extern from "astra/Float32ProjectionData3DMemory.h" namespace "astra": - cdef cppclass CFloat32ProjectionData3DMemory: - CFloat32ProjectionData3DMemory(CProjectionGeometry3D*) - CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*) - CFloat32ProjectionData3DMemory(CProjectionGeometry3D*, CFloat32CustomMemory*) - CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*, CFloat32CustomMemory*) - CProjectionGeometry3D* getGeometry() - void changeGeometry(CProjectionGeometry3D*) - int getDetectorColCount() - int getDetectorRowCount() - int getAngleCount() - bool isInitialized() + cdef cppclass CFloat32ProjectionData3DMemory: + CFloat32ProjectionData3DMemory(CProjectionGeometry3D*) + CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*) + CFloat32ProjectionData3DMemory(CProjectionGeometry3D*, CFloat32CustomMemory*) + CFloat32ProjectionData3DMemory(CConeProjectionGeometry3D*, CFloat32CustomMemory*) + CProjectionGeometry3D* getGeometry() + void changeGeometry(CProjectionGeometry3D*) + int getDetectorColCount() + int getDetectorRowCount() + int getAngleCount() + bool isInitialized() -cdef extern from "astra/Float32Data3D.h" namespace "astra": - cdef cppclass CFloat32Data3D: - CFloat32Data3D() - bool isInitialized() - int getSize() - int getWidth() - int getHeight() - int getDepth() - void updateStatistics() +IF HAVE_CUDA==True: + cdef extern from "astra/Float32VolumeData3DGPU.h" namespace "astra": + cdef cppclass CFloat32VolumeData3DGPU: + CFloat32VolumeData3DGPU(CVolumeGeometry3D*, MemHandle3D) + CVolumeGeometry3D* getGeometry() + void changeGeometry(CVolumeGeometry3D*) + int getRowCount() + int getColCount() + int getSliceCount() + bool isInitialized() + + cdef extern from "astra/Float32ProjectionData3DGPU.h" namespace "astra": + cdef cppclass CFloat32ProjectionData3DGPU: + CFloat32ProjectionData3DGPU(CProjectionGeometry3D*, MemHandle3D) + CProjectionGeometry3D* getGeometry() + void changeGeometry(CProjectionGeometry3D*) + int getRowCount() + int getColCount() + int getSliceCount() + bool isInitialized() diff --git a/python/astra/creators.py b/python/astra/creators.py index f45fd52..aa16221 100644 --- a/python/astra/creators.py +++ b/python/astra/creators.py @@ -258,7 +258,7 @@ This method can be called in a number of ways: elif intype == 'cone': if len(args) < 7: raise Exception('not enough variables: astra_create_proj_geom(cone, detector_spacing_x, detector_spacing_y, det_row_count, det_col_count, angles, source_origin, origin_det)') - return {'type': 'cone','DetectorSpacingX':args[0], 'DetectorSpacingY':args[1], 'DetectorRowCount':args[2],'DetectorColCount':args[3],'ProjectionAngles':args[4],'DistanceOriginSource': args[5],'DistanceOriginDetector':args[6]} + return {'type': 'cone','DetectorSpacingX':args[0], 'DetectorSpacingY':args[1], 'DetectorRowCount':args[2],'DetectorColCount':args[3],'ProjectionAngles':args[4],'DistanceOriginSource': args[5],'DistanceOriginDetector':args[6]} elif intype == 'cone_vec': if len(args) < 3: raise Exception('not enough variables: astra_create_proj_geom(cone_vec, det_row_count, det_col_count, V)') diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx index 203fde2..9c88073 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -83,8 +83,10 @@ def create(datatype, geometry, data=None, link=False): cdef CFloat32Data2D * pDataObject2D cdef CFloat32CustomMemory * pCustom - if link and data.shape!=geom_size(geometry): - raise Exception("The dimensions of the data do not match those specified in the geometry.") + if link: + geom_shape = geom_size(geometry) + if data.shape != geom_shape: + raise ValueError("The dimensions of the data do not match those specified in the geometry: {} != {}".format(data.shape, geom_shape)) if datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) @@ -92,7 +94,7 @@ def create(datatype, geometry, data=None, link=False): if not pGeometry.initialize(cfg[0]): del cfg del pGeometry - raise Exception('Geometry class not initialized.') + raise RuntimeError('Geometry class not initialized.') if link: pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) pDataObject2D = <CFloat32Data2D * > new CFloat32VolumeData2D(pGeometry, pCustom) @@ -114,7 +116,7 @@ def create(datatype, geometry, data=None, link=False): if not ppGeometry.initialize(cfg[0]): del cfg del ppGeometry - raise Exception('Geometry class not initialized.') + raise RuntimeError('Geometry class not initialized.') if link: pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) pDataObject2D = <CFloat32Data2D * > new CFloat32ProjectionData2D(ppGeometry, pCustom) @@ -123,11 +125,11 @@ def create(datatype, geometry, data=None, link=False): del ppGeometry del cfg else: - raise Exception("Invalid datatype. Please specify '-vol' or '-sino'.") + raise ValueError("Invalid datatype. Please specify '-vol' or '-sino'.") if not pDataObject2D.isInitialized(): del pDataObject2D - raise Exception("Couldn't initialize data object.") + raise RuntimeError("Couldn't initialize data object.") if not link: fillDataObject(pDataObject2D, data) @@ -138,6 +140,10 @@ cdef fillDataObject(CFloat32Data2D * obj, data): fillDataObjectScalar(obj, 0) else: if isinstance(data, np.ndarray): + obj_shape = (obj.getHeight(), obj.getWidth()) + if data.shape != obj_shape: + raise ValueError( + "The dimensions of the data do not match those specified in the geometry: {} != {}".format(data.shape, obj_shape)) fillDataObjectArray(obj, np.ascontiguousarray(data,dtype=np.float32)) else: fillDataObjectScalar(obj, np.float32(data)) @@ -150,18 +156,15 @@ cdef fillDataObjectScalar(CFloat32Data2D * obj, float s): @cython.boundscheck(False) @cython.wraparound(False) cdef fillDataObjectArray(CFloat32Data2D * obj, float [:,::1] data): - if (not data.shape[0] == obj.getHeight()) or (not data.shape[1] == obj.getWidth()): - raise Exception( - "The dimensions of the data do not match those specified in the geometry.") cdef float [:,::1] cView = <float[:data.shape[0],:data.shape[1]]> obj.getData2D()[0] cView[:] = data cdef CFloat32Data2D * getObject(i) except NULL: cdef CFloat32Data2D * pDataObject = man2d.get(i) if pDataObject == NULL: - raise Exception("Data object not found") + raise ValueError("Data object not found") if not pDataObject.isInitialized(): - raise Exception("Data object not initialized properly.") + raise RuntimeError("Data object not initialized properly.") return pDataObject @@ -180,15 +183,15 @@ def get_geometry(i): pDataObject3 = <CFloat32VolumeData2D * >pDataObject geom = utils.configToDict(pDataObject3.getGeometry().getConfiguration()) else: - raise Exception("Not a known data object") + raise RuntimeError("Not a known data object") return geom cdef CProjector2D * getProjector(i) except NULL: cdef CProjector2D * proj = manProj.get(i) if proj == NULL: - raise Exception("Projector not initialized.") + raise RuntimeError("Projector not initialized.") if not proj.isInitialized(): - raise Exception("Projector not initialized.") + raise RuntimeError("Projector not initialized.") return proj def check_compatible(i, proj_id): @@ -203,7 +206,7 @@ def check_compatible(i, proj_id): pDataObject3 = <CFloat32VolumeData2D * >pDataObject return pDataObject3.getGeometry().isEqual(proj.getVolumeGeometry()) else: - raise Exception("Not a known data object") + raise RuntimeError("Not a known data object") def change_geometry(i, geom): cdef Config *cfg @@ -227,12 +230,14 @@ def change_geometry(i, geom): if not ppGeometry.initialize(cfg[0]): del cfg del ppGeometry - raise Exception('Geometry class not initialized.') - if (ppGeometry.getDetectorCount() != pDataObject2.getDetectorCount() or ppGeometry.getProjectionAngleCount() != pDataObject2.getAngleCount()): + raise RuntimeError('Geometry class not initialized.') + geom_shape = (ppGeometry.getProjectionAngleCount(), ppGeometry.getDetectorCount()) + obj_shape = (pDataObject2.getAngleCount(), pDataObject2.getDetectorCount()) + if geom_shape != obj_shape: del ppGeometry del cfg - raise Exception( - "The dimensions of the data do not match those specified in the geometry.") + raise ValueError( + "The dimensions of the data do not match those specified in the geometry: {} != {}", obj_shape, geom_shape) pDataObject2.changeGeometry(ppGeometry) del ppGeometry del cfg @@ -243,17 +248,19 @@ def change_geometry(i, geom): if not pGeometry.initialize(cfg[0]): del cfg del pGeometry - raise Exception('Geometry class not initialized.') - if (pGeometry.getGridColCount() != pDataObject3.getWidth() or pGeometry.getGridRowCount() != pDataObject3.getHeight()): + raise RuntimeError('Geometry class not initialized.') + geom_shape = (pGeometry.getGridRowCount(), pGeometry.getGridColCount()) + obj_shape = (pDataObject3.getHeight(), pDataObject3.getWidth()) + if geom_shape != obj_shape: del cfg del pGeometry - raise Exception( - 'The dimensions of the data do not match those specified in the geometry.') + raise ValueError( + "The dimensions of the data do not match those specified in the geometry: {} != {}", obj_shape, geom_shape) pDataObject3.changeGeometry(pGeometry) del cfg del pGeometry else: - raise Exception("Not a known data object") + raise RuntimeError("Not a known data object") @cython.boundscheck(False) @cython.wraparound(False) @@ -274,7 +281,7 @@ def get_shared(i): def get_single(i): - raise Exception("Not yet implemented") + raise NotImplementedError("Not yet implemented") def info(): diff --git a/python/astra/data3d.py b/python/astra/data3d.py index 9c2bea4..ecb99d0 100644 --- a/python/astra/data3d.py +++ b/python/astra/data3d.py @@ -26,6 +26,8 @@ from . import data3d_c as d import numpy as np +from .pythonutils import GPULink + def create(datatype,geometry,data=None): """Create a 3D object. @@ -52,12 +54,13 @@ def link(datatype, geometry, data): :returns: :class:`int` -- the ID of the constructed object. """ - if not isinstance(data,np.ndarray): - raise ValueError("Input should be a numpy array") - if not data.dtype==np.float32: - raise ValueError("Numpy array should be float32") - if not (data.flags['C_CONTIGUOUS'] and data.flags['ALIGNED']): - raise ValueError("Numpy array should be C_CONTIGUOUS and ALIGNED") + if not isinstance(data,np.ndarray) and not isinstance(data,GPULink): + raise TypeError("Input should be a numpy ndarray or GPULink object") + if isinstance(data, np.ndarray): + if data.dtype != np.float32: + raise ValueError("Numpy array should be float32") + if not (data.flags['C_CONTIGUOUS'] and data.flags['ALIGNED']): + raise ValueError("Numpy array should be C_CONTIGUOUS and ALIGNED") return d.create(datatype,geometry,data,True) diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 915c60d..78ed620 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -45,12 +45,17 @@ from .PyXMLDocument cimport XMLDocument cimport utils from .utils import wrap_from_bytes -from .pythonutils import geom_size +from .pythonutils import geom_size, GPULink import operator from six.moves import reduce +include "config.pxi" + +cdef extern from "Python.h": + void* PyLong_AsVoidPtr(object) + cdef CData3DManager * man3d = <CData3DManager * >PyData3DManager.getSingletonPtr() @@ -65,12 +70,22 @@ def create(datatype,geometry,data=None, link=False): cdef Config *cfg cdef CVolumeGeometry3D * pGeometry cdef CProjectionGeometry3D * ppGeometry - cdef CFloat32Data3DMemory * pDataObject3D + cdef CFloat32Data3D * pDataObject3D cdef CConeProjectionGeometry3D* pppGeometry - cdef CFloat32CustomMemory * pCustom + cdef CFloat32CustomMemory * pCustom = NULL + IF HAVE_CUDA==True: + cdef MemHandle3D hnd - if link and data.shape!=geom_size(geometry): - raise Exception("The dimensions of the data do not match those specified in the geometry.") + if link: + geom_shape = geom_size(geometry) + if isinstance(data, np.ndarray): + data_shape = data.shape + elif isinstance(data, GPULink): + data_shape = ( data.z, data.y, data.x ) + else: + raise TypeError("data should be a numpy.ndarray or a GPULink object") + if geom_shape != data_shape: + raise ValueError("The dimensions of the data do not match those specified in the geometry: {} != {}".format(data_shape, geom_shape)) if datatype == '-vol': cfg = utils.dictToConfig(six.b('VolumeGeometry'), geometry) @@ -78,15 +93,25 @@ def create(datatype,geometry,data=None, link=False): if not pGeometry.initialize(cfg[0]): del cfg del pGeometry - raise Exception('Geometry class not initialized.') + raise RuntimeError('Geometry class not initialized.') if link: - pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) - pDataObject3D = <CFloat32Data3DMemory * > new CFloat32VolumeData3DMemory(pGeometry, pCustom) + if isinstance(data, np.ndarray): + pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) + pDataObject3D = <CFloat32Data3D * > new CFloat32VolumeData3DMemory(pGeometry, pCustom) + elif isinstance(data, GPULink): + IF HAVE_CUDA==True: + s = geom_size(geometry) + hnd = wrapHandle(<float*>PyLong_AsVoidPtr(data.ptr), data.x, data.y, data.z, data.pitch/4) + pDataObject3D = <CFloat32Data3D * > new CFloat32VolumeData3DGPU(pGeometry, hnd) + ELSE: + raise NotImplementedError("CUDA support is not enabled in ASTRA") + else: + raise TypeError("data should be a numpy.ndarray or a GPULink object") else: - pDataObject3D = <CFloat32Data3DMemory * > new CFloat32VolumeData3DMemory(pGeometry) + pDataObject3D = <CFloat32Data3D * > new CFloat32VolumeData3DMemory(pGeometry) del cfg del pGeometry - elif datatype == '-sino' or datatype == '-proj3d': + elif datatype == '-sino' or datatype == '-proj3d' or datatype == '-sinocone': cfg = utils.dictToConfig(six.b('ProjectionGeometry'), geometry) tpe = wrap_from_bytes(cfg.self.getAttribute(six.b('type'))) if (tpe == "parallel3d"): @@ -98,41 +123,38 @@ def create(datatype,geometry,data=None, link=False): elif (tpe == "cone_vec"): ppGeometry = <CProjectionGeometry3D*> new CConeVecProjectionGeometry3D(); else: - raise Exception("Invalid geometry type.") + raise ValueError("Invalid geometry type.") if not ppGeometry.initialize(cfg[0]): del cfg del ppGeometry - raise Exception('Geometry class not initialized.') + raise RuntimeError('Geometry class not initialized.') if link: - pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) - pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(ppGeometry, pCustom) + if isinstance(data, np.ndarray): + pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) + pDataObject3D = <CFloat32Data3D * > new CFloat32ProjectionData3DMemory(ppGeometry, pCustom) + elif isinstance(data, GPULink): + IF HAVE_CUDA==True: + s = geom_size(geometry) + hnd = wrapHandle(<float*>PyLong_AsVoidPtr(data.ptr), data.x, data.y, data.z, data.pitch/4) + pDataObject3D = <CFloat32Data3D * > new CFloat32ProjectionData3DGPU(ppGeometry, hnd) + ELSE: + raise NotImplementedError("CUDA support is not enabled in ASTRA") + else: + raise TypeError("data should be a numpy.ndarray or a GPULink object") else: pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(ppGeometry) del ppGeometry del cfg - elif datatype == "-sinocone": - cfg = utils.dictToConfig(six.b('ProjectionGeometry'), geometry) - pppGeometry = new CConeProjectionGeometry3D() - if not pppGeometry.initialize(cfg[0]): - del cfg - del pppGeometry - raise Exception('Geometry class not initialized.') - if link: - pCustom = <CFloat32CustomMemory*> new CFloat32CustomPython(data) - pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(pppGeometry, pCustom) - else: - pDataObject3D = <CFloat32Data3DMemory * > new CFloat32ProjectionData3DMemory(pppGeometry) else: - raise Exception("Invalid datatype. Please specify '-vol' or '-proj3d'.") + raise ValueError("Invalid datatype. Please specify '-vol' or '-proj3d'.") if not pDataObject3D.isInitialized(): del pDataObject3D - raise Exception("Couldn't initialize data object.") - - if not link: fillDataObject(pDataObject3D, data) + raise RuntimeError("Couldn't initialize data object.") - pDataObject3D.updateStatistics() + if not link: + fillDataObject(dynamic_cast_mem(pDataObject3D), data) return man3d.store(<CFloat32Data3D*>pDataObject3D) @@ -147,7 +169,7 @@ def get_geometry(i): pDataObject3 = <CFloat32VolumeData3DMemory * >pDataObject geom = utils.configToDict(pDataObject3.getGeometry().getConfiguration()) else: - raise Exception("Not a known data object") + raise RuntimeError("Not a known data object") return geom def change_geometry(i, geom): @@ -168,18 +190,18 @@ def change_geometry(i, geom): elif (tpe == "cone_vec"): ppGeometry = <CProjectionGeometry3D*> new CConeVecProjectionGeometry3D(); else: - raise Exception("Invalid geometry type.") + raise ValueError("Invalid geometry type.") if not ppGeometry.initialize(cfg[0]): del cfg del ppGeometry - raise Exception('Geometry class not initialized.') + raise RuntimeError('Geometry class not initialized.') del cfg - if (ppGeometry.getDetectorColCount() != pDataObject2.getDetectorColCount() or \ - ppGeometry.getProjectionCount() != pDataObject2.getAngleCount() or \ - ppGeometry.getDetectorRowCount() != pDataObject2.getDetectorRowCount()): + geom_shape = (ppGeometry.getDetectorRowCount(), ppGeometry.getProjectionCount(), ppGeometry.getDetectorColCount()) + obj_shape = (pDataObject2.getDetectorRowCount(), pDataObject2.getAngleCount(), pDataObject2.getDetectorColCount()) + if geom_shape != obj_shape: del ppGeometry - raise Exception( - "The dimensions of the data do not match those specified in the geometry.") + raise ValueError( + "The dimensions of the data do not match those specified in the geometry: {} != {}".format(obj_shape, geom_shape)) pDataObject2.changeGeometry(ppGeometry) del ppGeometry @@ -190,19 +212,19 @@ def change_geometry(i, geom): if not pGeometry.initialize(cfg[0]): del cfg del pGeometry - raise Exception('Geometry class not initialized.') + raise RuntimeError('Geometry class not initialized.') del cfg - if (pGeometry.getGridColCount() != pDataObject3.getColCount() or \ - pGeometry.getGridRowCount() != pDataObject3.getRowCount() or \ - pGeometry.getGridSliceCount() != pDataObject3.getSliceCount()): + geom_shape = (pGeometry.getGridSliceCount(), pGeometry.getGridRowCount(), pGeometry.getGridColCount()) + obj_shape = (pDataObject3.getSliceCount(), pDataObject3.getRowCount(), pDataObject3.getColCount()) + if geom_shape != obj_shape: del pGeometry - raise Exception( - "The dimensions of the data do not match those specified in the geometry.") + raise ValueError( + "The dimensions of the data do not match those specified in the geometry.".format(obj_shape, geom_shape)) pDataObject3.changeGeometry(pGeometry) del pGeometry else: - raise Exception("Not a known data object") + raise RuntimeError("Not a known data object") cdef fillDataObject(CFloat32Data3DMemory * obj, data): @@ -210,6 +232,10 @@ cdef fillDataObject(CFloat32Data3DMemory * obj, data): fillDataObjectScalar(obj, 0) else: if isinstance(data, np.ndarray): + obj_shape = (obj.getDepth(), obj.getHeight(), obj.getWidth()) + if data.shape != obj_shape: + raise ValueError( + "The dimensions of the data do not match those specified in the geometry: {} != {}".format(data.shape, obj_shape)) fillDataObjectArray(obj, np.ascontiguousarray(data,dtype=np.float32)) else: fillDataObjectScalar(obj, np.float32(data)) @@ -222,18 +248,15 @@ cdef fillDataObjectScalar(CFloat32Data3DMemory * obj, float s): @cython.boundscheck(False) @cython.wraparound(False) cdef fillDataObjectArray(CFloat32Data3DMemory * obj, float [:,:,::1] data): - if (not data.shape[0] == obj.getDepth()) or (not data.shape[1] == obj.getHeight()) or (not data.shape[2] == obj.getWidth()): - raise Exception( - "The dimensions of the data do not match those specified in the geometry.") cdef float [:,:,::1] cView = <float[:data.shape[0],:data.shape[1],:data.shape[2]]> obj.getData3D()[0][0] cView[:] = data cdef CFloat32Data3D * getObject(i) except NULL: cdef CFloat32Data3D * pDataObject = man3d.get(i) if pDataObject == NULL: - raise Exception("Data object not found") + raise ValueError("Data object not found") if not pDataObject.isInitialized(): - raise Exception("Data object not initialized properly.") + raise RuntimeError("Data object not initialized properly.") return pDataObject @cython.boundscheck(False) @@ -256,7 +279,7 @@ def get_shared(i): return np.PyArray_SimpleNewFromData(3,shape,np.NPY_FLOAT32,<void *>pDataObject.getData3D()[0][0]) def get_single(i): - raise Exception("Not yet implemented") + raise NotImplementedError("Not yet implemented") def store(i,data): cdef CFloat32Data3D * pDataObject = getObject(i) diff --git a/python/astra/experimental.pyx b/python/astra/experimental.pyx index b6c7881..0af3118 100644 --- a/python/astra/experimental.pyx +++ b/python/astra/experimental.pyx @@ -36,12 +36,20 @@ IF HAVE_CUDA==True: cdef extern from "astra/CompositeGeometryManager.h" namespace "astra": cdef cppclass CCompositeGeometryManager: - bool doFP(CProjector3D *, vector[CFloat32VolumeData3DMemory *], vector[CFloat32ProjectionData3DMemory *]) - bool doBP(CProjector3D *, vector[CFloat32VolumeData3DMemory *], vector[CFloat32ProjectionData3DMemory *]) + bool doFP(CProjector3D *, vector[CFloat32VolumeData3D *], vector[CFloat32ProjectionData3D *]) + bool doBP(CProjector3D *, vector[CFloat32VolumeData3D *], vector[CFloat32ProjectionData3D *]) cdef extern from *: - CFloat32VolumeData3DMemory * dynamic_cast_vol_mem "dynamic_cast<astra::CFloat32VolumeData3DMemory*>" (CFloat32Data3D * ) except NULL - CFloat32ProjectionData3DMemory * dynamic_cast_proj_mem "dynamic_cast<astra::CFloat32ProjectionData3DMemory*>" (CFloat32Data3D * ) except NULL + CFloat32VolumeData3D * dynamic_cast_vol_mem "dynamic_cast<astra::CFloat32VolumeData3D*>" (CFloat32Data3D * ) except NULL + CFloat32ProjectionData3D * dynamic_cast_proj_mem "dynamic_cast<astra::CFloat32ProjectionData3D*>" (CFloat32Data3D * ) except NULL + + cdef extern from "astra/Float32ProjectionData3D.h" namespace "astra": + cdef cppclass CFloat32ProjectionData3D: + bool isInitialized() + cdef extern from "astra/Float32VolumeData3D.h" namespace "astra": + cdef cppclass CFloat32VolumeData3D: + bool isInitialized() + cimport PyProjector3DManager from .PyProjector3DManager cimport CProjector3DManager @@ -52,9 +60,9 @@ IF HAVE_CUDA==True: cdef CData3DManager * man3d = <CData3DManager * >PyData3DManager.getSingletonPtr() def do_composite(projector_id, vol_ids, proj_ids, t): - cdef vector[CFloat32VolumeData3DMemory *] vol - cdef CFloat32VolumeData3DMemory * pVolObject - cdef CFloat32ProjectionData3DMemory * pProjObject + cdef vector[CFloat32VolumeData3D *] vol + cdef CFloat32VolumeData3D * pVolObject + cdef CFloat32ProjectionData3D * pProjObject for v in vol_ids: pVolObject = dynamic_cast_vol_mem(man3d.get(v)) if pVolObject == NULL: @@ -62,7 +70,7 @@ IF HAVE_CUDA==True: if not pVolObject.isInitialized(): raise Exception("Data object not initialized properly") vol.push_back(pVolObject) - cdef vector[CFloat32ProjectionData3DMemory *] proj + cdef vector[CFloat32ProjectionData3D *] proj for v in proj_ids: pProjObject = dynamic_cast_proj_mem(man3d.get(v)) if pProjObject == NULL: diff --git a/python/astra/extrautils.pyx b/python/astra/extrautils.pyx index 502cc78..248c6ee 100644 --- a/python/astra/extrautils.pyx +++ b/python/astra/extrautils.pyx @@ -26,19 +26,19 @@ def clipCircle(img): - cdef int i,j - cdef double x2,y2,mid,bnd - cdef long sz,sz2 - sz = img.shape[0] - sz2 = sz*sz - bnd = sz2/4. - mid = (sz-1.)/2. - nDel=0 - for i in range(sz): - for j in range(sz): - x2 = (i-mid)*(i-mid) - y2 = (j-mid)*(j-mid) - if x2+y2>bnd: - img[i,j]=0 - nDel=nDel+1 - return nDel + cdef int i,j + cdef double x2,y2,mid,bnd + cdef long sz,sz2 + sz = img.shape[0] + sz2 = sz*sz + bnd = sz2/4. + mid = (sz-1.)/2. + nDel=0 + for i in range(sz): + for j in range(sz): + x2 = (i-mid)*(i-mid) + y2 = (j-mid)*(j-mid) + if x2+y2>bnd: + img[i,j]=0 + nDel=nDel+1 + return nDel diff --git a/python/astra/pythonutils.py b/python/astra/pythonutils.py index 3bd3321..27fa8fd 100644 --- a/python/astra/pythonutils.py +++ b/python/astra/pythonutils.py @@ -61,3 +61,21 @@ def geom_size(geom, dim=None): s = s[dim] return s + +class GPULink(object): + """Utility class for astra.data3d.link with a CUDA pointer + + The CUDA pointer ptr must point to an array of floats. + + x is the fastest-changing coordinate, z the slowest-changing. + + pitch is the width in bytes of the memory block. For a contiguous + memory block, pitch is equal to sizeof(float) * x. For a memory block + allocated by cudaMalloc3D, pitch is the pitch as returned by cudaMalloc3D. + """ + def __init__(self, ptr, x, y, z, pitch): + self.ptr = ptr + self.x = x + self.y = y + self.z = z + self.pitch = pitch diff --git a/samples/python/s021_pygpu.py b/samples/python/s021_pygpu.py new file mode 100644 index 0000000..db6584b --- /dev/null +++ b/samples/python/s021_pygpu.py @@ -0,0 +1,72 @@ +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam +# +# Contact: astra@uantwerpen.be +# Website: http://www.astra-toolbox.com/ +# +# 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/>. +# +# ----------------------------------------------------------------------- + +import astra +import numpy as np +import pygpu +import pylab + +# Initialize pygpu +ctx = pygpu.init('cuda') +pygpu.set_default_context(ctx) + +vol_geom = astra.create_vol_geom(128, 128, 128) +angles = np.linspace(0, 2 * np.pi, 180, False) +proj_geom = astra.create_proj_geom('cone', 1.0, 1.0, 128, 192, angles, 1000, 0) + +# Create a simple hollow cube phantom, as a pygpu gpuarray +vol_gpuarr = pygpu.gpuarray.zeros(astra.functions.geom_size(vol_geom), dtype='float32') +vol_gpuarr[17:113, 17:113, 17:113] = 1 +vol_gpuarr[33:97, 33:97, 33:97] = 0 + +# Create a pygpu gpuarray for the output projection data +proj_gpuarr = pygpu.gpuarray.zeros(astra.functions.geom_size(proj_geom), dtype='float32') + +# Create the astra GPULink objects and create astra data3d objects from them +z, y, x = proj_gpuarr.shape +proj_data_link = astra.data3d.GPULink(proj_gpuarr.gpudata, x, y, z, + proj_gpuarr.strides[-2]) +z, y, x = vol_gpuarr.shape +vol_link = astra.data3d.GPULink(vol_gpuarr.gpudata, x, y, z, + vol_gpuarr.strides[-2]) + +proj_id = astra.data3d.link('-sino', proj_geom, proj_data_link) +vol_id = astra.data3d.link('-vol', vol_geom, vol_link) + +# Run a 3D FP +cfg = astra.astra_dict('FP3D_CUDA') +cfg['VolumeDataId'] = vol_id +cfg['ProjectionDataId'] = proj_id +alg_id = astra.algorithm.create(cfg) +astra.algorithm.run(alg_id) + +pylab.figure(1) +pylab.gray() +pylab.imshow(proj_gpuarr[:, 20, :]) +pylab.show() + +astra.algorithm.delete(alg_id) +astra.data3d.delete(vol_id) +astra.data3d.delete(proj_id) diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index c3af228..a4dace2 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -39,6 +39,8 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. #include "astra/CudaProjector3D.h" #include "astra/Float32ProjectionData3DMemory.h" #include "astra/Float32VolumeData3DMemory.h" +#include "astra/Float32ProjectionData3DGPU.h" +#include "astra/Float32VolumeData3DGPU.h" #include "astra/Logging.h" #include "../cuda/3d/mem3d.h" @@ -97,6 +99,127 @@ CCompositeGeometryManager::CCompositeGeometryManager() // (First approach: 0.5/0.5) + + + +class _AstraExport CFloat32CustomGPUMemory { +public: + astraCUDA3d::MemHandle3D hnd; // Only required to be valid between allocate/free + virtual bool allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, astraCUDA3d::Mem3DZeroMode zero)=0; + virtual bool copyToGPUMemory(const astraCUDA3d::SSubDimensions3D &pos)=0; + virtual bool copyFromGPUMemory(const astraCUDA3d::SSubDimensions3D &pos)=0; + virtual bool freeGPUMemory()=0; + virtual ~CFloat32CustomGPUMemory() { } +}; + +class CFloat32ExistingGPUMemory : public astra::CFloat32CustomGPUMemory { +public: + CFloat32ExistingGPUMemory(CFloat32Data3DGPU *d); + virtual bool allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, astraCUDA3d::Mem3DZeroMode zero); + virtual bool copyToGPUMemory(const astraCUDA3d::SSubDimensions3D &pos); + virtual bool copyFromGPUMemory(const astraCUDA3d::SSubDimensions3D &pos); + virtual bool freeGPUMemory(); + +protected: + unsigned int x, y, z; +}; + +class CFloat32DefaultGPUMemory : public astra::CFloat32CustomGPUMemory { +public: + CFloat32DefaultGPUMemory(CFloat32Data3DMemory* d) { + ptr = d->getData(); + } + virtual bool allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, astraCUDA3d::Mem3DZeroMode zero) { + hnd = astraCUDA3d::allocateGPUMemory(x, y, z, zero); + return (bool)hnd; + } + virtual bool copyToGPUMemory(const astraCUDA3d::SSubDimensions3D &pos) { + return astraCUDA3d::copyToGPUMemory(ptr, hnd, pos); + } + virtual bool copyFromGPUMemory(const astraCUDA3d::SSubDimensions3D &pos) { + return astraCUDA3d::copyFromGPUMemory(ptr, hnd, pos); + } + virtual bool freeGPUMemory() { + return astraCUDA3d::freeGPUMemory(hnd); + } + +protected: + float *ptr; +}; + + + +CFloat32ExistingGPUMemory::CFloat32ExistingGPUMemory(CFloat32Data3DGPU *d) +{ + hnd = d->getHandle(); + x = d->getWidth(); + y = d->getHeight(); + z = d->getDepth(); +} + +bool CFloat32ExistingGPUMemory::allocateGPUMemory(unsigned int x_, unsigned int y_, unsigned int z_, astraCUDA3d::Mem3DZeroMode zero) { + assert(x_ == x); + assert(y_ == y); + assert(z_ == z); + + if (zero == astraCUDA3d::INIT_ZERO) + return astraCUDA3d::zeroGPUMemory(hnd, x, y, z); + else + return true; +} +bool CFloat32ExistingGPUMemory::copyToGPUMemory(const astraCUDA3d::SSubDimensions3D &pos) { + assert(pos.nx == x); + assert(pos.ny == y); + assert(pos.nz == z); + assert(pos.pitch == x); + assert(pos.subx == 0); + assert(pos.suby == 0); + assert(pos.subnx == x); + assert(pos.subny == y); + + // These are less necessary than x/y, but allowing access to + // subvolumes needs an interface change + assert(pos.subz == 0); + assert(pos.subnz == z); + + return true; +} +bool CFloat32ExistingGPUMemory::copyFromGPUMemory(const astraCUDA3d::SSubDimensions3D &pos) { + assert(pos.nx == x); + assert(pos.ny == y); + assert(pos.nz == z); + assert(pos.pitch == x); + assert(pos.subx == 0); + assert(pos.suby == 0); + assert(pos.subnx == x); + assert(pos.subny == y); + + // These are less necessary than x/y, but allowing access to + // subvolumes needs an interface change + assert(pos.subz == 0); + assert(pos.subnz == z); + + return true; +} +bool CFloat32ExistingGPUMemory::freeGPUMemory() { + return true; +} + + +CFloat32CustomGPUMemory * createGPUMemoryHandler(CFloat32Data3D *d) { + CFloat32Data3DMemory *dMem = dynamic_cast<CFloat32Data3DMemory*>(d); + CFloat32Data3DGPU *dGPU = dynamic_cast<CFloat32Data3DGPU*>(d); + + if (dMem) + return new CFloat32DefaultGPUMemory(dMem); + else + return new CFloat32ExistingGPUMemory(dGPU); +} + + + + + bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div, TJobSet &split) { int maxBlockDim = astraCUDA3d::maxBlockDimension(); @@ -205,6 +328,7 @@ static std::pair<double, double> reduceProjectionVertical(const CVolumeGeometry3 // reduce self to only cover intersection with projection of VolumePart // (Project corners of volume, take bounding box) + assert(pProjGeom->getProjectionCount() > 0); for (int i = 0; i < pProjGeom->getProjectionCount(); ++i) { double vol_u[8]; @@ -280,7 +404,7 @@ CCompositeGeometryManager::CVolumePart::~CVolumePart() delete pGeom; } -void CCompositeGeometryManager::CVolumePart::getDims(size_t &x, size_t &y, size_t &z) +void CCompositeGeometryManager::CVolumePart::getDims(size_t &x, size_t &y, size_t &z) const { if (!pGeom) { x = y = z = 0; @@ -292,13 +416,28 @@ void CCompositeGeometryManager::CVolumePart::getDims(size_t &x, size_t &y, size_ z = pGeom->getGridSliceCount(); } -size_t CCompositeGeometryManager::CPart::getSize() +size_t CCompositeGeometryManager::CPart::getSize() const { size_t x, y, z; getDims(x, y, z); return x * y * z; } +bool CCompositeGeometryManager::CPart::isFull() const +{ + size_t x, y, z; + getDims(x, y, z); + return x == (size_t)pData->getWidth() && + y == (size_t)pData->getHeight() && + z == (size_t)pData->getDepth(); +} + +bool CCompositeGeometryManager::CPart::canSplitAndReduce() const +{ + return dynamic_cast<CFloat32Data3DMemory *>(pData) != 0; +} + + static bool testVolumeRange(const std::pair<double, double>& fullRange, const CVolumeGeometry3D *pVolGeom, @@ -334,6 +473,9 @@ static bool testVolumeRange(const std::pair<double, double>& fullRange, CCompositeGeometryManager::CPart* CCompositeGeometryManager::CVolumePart::reduce(const CPart *_other) { + if (!canSplitAndReduce()) + return clone(); + const CProjectionPart *other = dynamic_cast<const CProjectionPart *>(_other); assert(other); @@ -654,7 +796,7 @@ static CProjectionGeometry3D* getSubProjectionGeometryV(const CProjectionGeometr // - maybe all approximately the same size? void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div) { - if (true) { + if (canSplitAndReduce()) { // Split in vertical direction only at first, until we figure out // a model for splitting in other directions @@ -664,7 +806,7 @@ void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::T size_t blockSize = computeLinearSplit(m, div, sliceCount); int rem = blockSize - (sliceCount % blockSize); - if (rem == blockSize) + if ((size_t)rem == blockSize) rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -698,12 +840,14 @@ void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::T out.push_back(boost::shared_ptr<CPart>(sub)); } + } else { + out.push_back(boost::shared_ptr<CPart>(clone())); } } void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div) { - if (true) { + if (canSplitAndReduce()) { // Split in vertical direction only at first, until we figure out // a model for splitting in other directions @@ -713,7 +857,7 @@ void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::T size_t blockSize = computeLinearSplit(m, div, sliceCount); int rem = blockSize - (sliceCount % blockSize); - if (rem == blockSize) + if ((size_t)rem == blockSize) rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -747,12 +891,14 @@ void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::T out.push_back(boost::shared_ptr<CPart>(sub)); } + } else { + out.push_back(boost::shared_ptr<CPart>(clone())); } } void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div) { - if (true) { + if (canSplitAndReduce()) { // Split in vertical direction only at first, until we figure out // a model for splitting in other directions @@ -762,7 +908,7 @@ void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::T size_t blockSize = computeLinearSplit(m, div, sliceCount); int rem = blockSize - (sliceCount % blockSize); - if (rem == blockSize) + if ((size_t)rem == blockSize) rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -796,6 +942,8 @@ void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::T out.push_back(boost::shared_ptr<CPart>(sub)); } + } else { + out.push_back(boost::shared_ptr<CPart>(clone())); } } @@ -815,7 +963,7 @@ CCompositeGeometryManager::CProjectionPart::~CProjectionPart() delete pGeom; } -void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, size_t &z) +void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, size_t &z) const { if (!pGeom) { x = y = z = 0; @@ -831,6 +979,9 @@ void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, s CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::reduce(const CPart *_other) { + if (!canSplitAndReduce()) + return clone(); + const CVolumePart *other = dynamic_cast<const CVolumePart *>(_other); assert(other); @@ -868,7 +1019,7 @@ CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::re void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div) { - if (true) { + if (canSplitAndReduce()) { // Split in vertical direction only at first, until we figure out // a model for splitting in other directions @@ -878,7 +1029,7 @@ void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManage size_t blockSize = computeLinearSplit(m, div, sliceCount); int rem = blockSize - (sliceCount % blockSize); - if (rem == blockSize) + if ((size_t)rem == blockSize) rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -903,6 +1054,8 @@ void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManage out.push_back(boost::shared_ptr<CPart>(sub)); } + } else { + out.push_back(boost::shared_ptr<CPart>(clone())); } } @@ -914,7 +1067,7 @@ void CCompositeGeometryManager::CProjectionPart::splitY(CCompositeGeometryManage void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div) { - if (true) { + if (canSplitAndReduce()) { // Split in vertical direction only at first, until we figure out // a model for splitting in other directions @@ -924,7 +1077,7 @@ void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManage size_t blockSize = computeLinearSplit(m, div, sliceCount); int rem = blockSize - (sliceCount % blockSize); - if (rem == blockSize) + if ((size_t)rem == blockSize) rem = 0; ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); @@ -949,6 +1102,8 @@ void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManage out.push_back(boost::shared_ptr<CPart>(sub)); } + } else { + out.push_back(boost::shared_ptr<CPart>(clone())); } } @@ -959,8 +1114,8 @@ CCompositeGeometryManager::CProjectionPart* CCompositeGeometryManager::CProjecti } CCompositeGeometryManager::SJob CCompositeGeometryManager::createJobFP(CProjector3D *pProjector, - CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData) + CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData) { ASTRA_DEBUG("CCompositeGeometryManager::createJobFP"); // Create single job for FP @@ -992,8 +1147,8 @@ CCompositeGeometryManager::SJob CCompositeGeometryManager::createJobFP(CProjecto } CCompositeGeometryManager::SJob CCompositeGeometryManager::createJobBP(CProjector3D *pProjector, - CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData) + CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData) { ASTRA_DEBUG("CCompositeGeometryManager::createJobBP"); // Create single job for BP @@ -1022,8 +1177,8 @@ CCompositeGeometryManager::SJob CCompositeGeometryManager::createJobBP(CProjecto return BP; } -bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData) +bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData) { TJobList L; L.push_back(createJobFP(pProjector, pVolData, pProjData)); @@ -1031,8 +1186,8 @@ bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, CFloat32VolumeDat return doJobs(L); } -bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData) +bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData) { TJobList L; L.push_back(createJobBP(pProjector, pVolData, pProjData)); @@ -1041,8 +1196,8 @@ bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeDat } -bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, - CFloat32ProjectionData3DMemory *pProjData, bool bShortScan, +bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeData3D *pVolData, + CFloat32ProjectionData3D *pProjData, bool bShortScan, const float *pfFilter) { if (!dynamic_cast<CConeProjectionGeometry3D*>(pProjData->getGeometry())) { @@ -1061,11 +1216,11 @@ bool CCompositeGeometryManager::doFDK(CProjector3D *pProjector, CFloat32VolumeDa return doJobs(L); } -bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData) +bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3D *>& volData, const std::vector<CFloat32ProjectionData3D *>& projData) { ASTRA_DEBUG("CCompositeGeometryManager::doFP, multi-volume"); - std::vector<CFloat32VolumeData3DMemory *>::const_iterator i; + std::vector<CFloat32VolumeData3D *>::const_iterator i; std::vector<boost::shared_ptr<CPart> > inputs; for (i = volData.begin(); i != volData.end(); ++i) { @@ -1079,7 +1234,7 @@ bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector inputs.push_back(boost::shared_ptr<CPart>(input)); } - std::vector<CFloat32ProjectionData3DMemory *>::const_iterator j; + std::vector<CFloat32ProjectionData3D *>::const_iterator j; std::vector<boost::shared_ptr<CPart> > outputs; for (j = projData.begin(); j != projData.end(); ++j) { @@ -1115,12 +1270,12 @@ bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector return doJobs(L); } -bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData) +bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3D *>& volData, const std::vector<CFloat32ProjectionData3D *>& projData) { ASTRA_DEBUG("CCompositeGeometryManager::doBP, multi-volume"); - std::vector<CFloat32VolumeData3DMemory *>::const_iterator i; + std::vector<CFloat32VolumeData3D *>::const_iterator i; std::vector<boost::shared_ptr<CPart> > outputs; for (i = volData.begin(); i != volData.end(); ++i) { @@ -1134,7 +1289,7 @@ bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, const std::vector outputs.push_back(boost::shared_ptr<CPart>(output)); } - std::vector<CFloat32ProjectionData3DMemory *>::const_iterator j; + std::vector<CFloat32ProjectionData3D *>::const_iterator j; std::vector<boost::shared_ptr<CPart> > inputs; for (j = projData.begin(); j != projData.end(); ++j) { @@ -1188,14 +1343,25 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter if (L.begin()->eType == CCompositeGeometryManager::SJob::JOB_NOP) { // just zero output? if (zero) { - for (size_t z = 0; z < outz; ++z) { - for (size_t y = 0; y < outy; ++y) { - float* ptr = output->pData->getData(); - ptr += (z + output->subX) * (size_t)output->pData->getHeight() * (size_t)output->pData->getWidth(); - ptr += (y + output->subY) * (size_t)output->pData->getWidth(); - ptr += output->subX; - memset(ptr, 0, sizeof(float) * outx); + // TODO: This function shouldn't have to know about this difference + // between Memory/GPU + CFloat32Data3DMemory *hostMem = dynamic_cast<CFloat32Data3DMemory *>(output->pData); + if (hostMem) { + for (size_t z = 0; z < outz; ++z) { + for (size_t y = 0; y < outy; ++y) { + float* ptr = hostMem->getData(); + ptr += (z + output->subX) * (size_t)output->pData->getHeight() * (size_t)output->pData->getWidth(); + ptr += (y + output->subY) * (size_t)output->pData->getWidth(); + ptr += output->subX; + memset(ptr, 0, sizeof(float) * outx); + } } + } else { + CFloat32Data3DGPU *gpuMem = dynamic_cast<CFloat32Data3DGPU *>(output->pData); + assert(gpuMem); + assert(output->isFull()); // TODO: zero subset? + + zeroGPUMemory(gpuMem->getHandle(), outx, outy, outz); } } return true; @@ -1214,10 +1380,11 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter dstdims.subx = output->subX; dstdims.suby = output->subY; dstdims.subz = output->subZ; - float *dst = output->pData->getData(); - astraCUDA3d::MemHandle3D outputMem = astraCUDA3d::allocateGPUMemory(outx, outy, outz, zero ? astraCUDA3d::INIT_ZERO : astraCUDA3d::INIT_NO); - bool ok = outputMem; + CFloat32CustomGPUMemory *dstMem = createGPUMemoryHandler(output->pData); + + bool ok = dstMem->allocateGPUMemory(outx, outy, outz, zero ? astraCUDA3d::INIT_ZERO : astraCUDA3d::INIT_NO); + if (!ok) ASTRA_ERROR("Error allocating GPU memory"); for (CCompositeGeometryManager::TJobList::const_iterator i = L.begin(); i != L.end(); ++i) { const CCompositeGeometryManager::SJob &j = *i; @@ -1238,7 +1405,8 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter size_t inx, iny, inz; j.pInput->getDims(inx, iny, inz); - astraCUDA3d::MemHandle3D inputMem = astraCUDA3d::allocateGPUMemory(inx, iny, inz, astraCUDA3d::INIT_NO); + + CFloat32CustomGPUMemory *srcMem = createGPUMemoryHandler(j.pInput->pData); astraCUDA3d::SSubDimensions3D srcdims; srcdims.nx = j.pInput->pData->getWidth(); @@ -1251,9 +1419,11 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter srcdims.subx = j.pInput->subX; srcdims.suby = j.pInput->subY; srcdims.subz = j.pInput->subZ; - const float *src = j.pInput->pData->getDataConst(); - ok = astraCUDA3d::copyToGPUMemory(src, inputMem, srcdims); + ok = srcMem->allocateGPUMemory(inx, iny, inz, astraCUDA3d::INIT_NO); + if (!ok) ASTRA_ERROR("Error allocating GPU memory"); + + ok = srcMem->copyToGPUMemory(srcdims); if (!ok) ASTRA_ERROR("Error copying input data to GPU"); switch (j.eType) { @@ -1264,7 +1434,7 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing FP"); - ok = astraCUDA3d::FP(((CCompositeGeometryManager::CProjectionPart*)j.pOutput.get())->pGeom, outputMem, ((CCompositeGeometryManager::CVolumePart*)j.pInput.get())->pGeom, inputMem, detectorSuperSampling, projKernel); + ok = astraCUDA3d::FP(((CCompositeGeometryManager::CProjectionPart*)j.pOutput.get())->pGeom, dstMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pInput.get())->pGeom, srcMem->hnd, detectorSuperSampling, projKernel); if (!ok) ASTRA_ERROR("Error performing sub-FP"); ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FP done"); } @@ -1276,7 +1446,7 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing BP"); - ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, outputMem, voxelSuperSampling, densityWeighting); + ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, srcMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, dstMem->hnd, voxelSuperSampling, densityWeighting); if (!ok) ASTRA_ERROR("Error performing sub-BP"); ASTRA_DEBUG("CCompositeGeometryManager::doJobs: BP done"); } @@ -1292,7 +1462,7 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter } else { ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing FDK"); - ok = astraCUDA3d::FDK(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, outputMem, j.FDKSettings.bShortScan, j.FDKSettings.pfFilter); + ok = astraCUDA3d::FDK(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, srcMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, dstMem->hnd, j.FDKSettings.bShortScan, j.FDKSettings.pfFilter); if (!ok) ASTRA_ERROR("Error performing sub-FDK"); ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FDK done"); } @@ -1302,17 +1472,20 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter assert(false); } - ok = astraCUDA3d::freeGPUMemory(inputMem); + ok = srcMem->freeGPUMemory(); if (!ok) ASTRA_ERROR("Error freeing GPU memory"); + delete srcMem; } - ok = astraCUDA3d::copyFromGPUMemory(dst, outputMem, dstdims); + ok = dstMem->copyFromGPUMemory(dstdims); if (!ok) ASTRA_ERROR("Error copying output data from GPU"); - ok = astraCUDA3d::freeGPUMemory(outputMem); + ok = dstMem->freeGPUMemory(); if (!ok) ASTRA_ERROR("Error freeing GPU memory"); + delete dstMem; + return true; } @@ -1455,6 +1628,8 @@ void CCompositeGeometryManager::setGPUIndices(const std::vector<int>& GPUIndices bool CCompositeGeometryManager::doJobs(TJobList &jobs) { + // TODO: Proper clean up if substeps fail (Or as proper as possible) + ASTRA_DEBUG("CCompositeGeometryManager::doJobs"); // Sort job list into job set by output part diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp index 223a9a4..27bb968 100644 --- a/src/CudaBackProjectionAlgorithm3D.cpp +++ b/src/CudaBackProjectionAlgorithm3D.cpp @@ -60,8 +60,8 @@ CCudaBackProjectionAlgorithm3D::CCudaBackProjectionAlgorithm3D() //---------------------------------------------------------------------------------------- // Constructor with initialization CCudaBackProjectionAlgorithm3D::CCudaBackProjectionAlgorithm3D(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pProjectionData, - CFloat32VolumeData3DMemory* _pReconstruction) + CFloat32ProjectionData3D* _pProjectionData, + CFloat32VolumeData3D* _pReconstruction) { _clear(); initialize(_pProjector, _pProjectionData, _pReconstruction); @@ -145,8 +145,8 @@ bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg) //---------------------------------------------------------------------------------------- // Initialize - C++ bool CCudaBackProjectionAlgorithm3D::initialize(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pSinogram, - CFloat32VolumeData3DMemory* _pReconstruction) + CFloat32ProjectionData3D* _pSinogram, + CFloat32VolumeData3D* _pReconstruction) { // if already initialized, clear first if (m_bIsInitialized) { @@ -187,17 +187,21 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations) // check initialized ASTRA_ASSERT(m_bIsInitialized); - CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram); + CFloat32ProjectionData3D* pSinoMem = dynamic_cast<CFloat32ProjectionData3D*>(m_pSinogram); ASTRA_ASSERT(pSinoMem); - CFloat32VolumeData3DMemory* pReconMem = dynamic_cast<CFloat32VolumeData3DMemory*>(m_pReconstruction); + CFloat32VolumeData3D* pReconMem = dynamic_cast<CFloat32VolumeData3D*>(m_pReconstruction); ASTRA_ASSERT(pReconMem); const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry(); const CVolumeGeometry3D& volgeom = *pReconMem->getGeometry(); if (m_bSIRTWeighting) { - astraCudaBP_SIRTWeighted(pReconMem->getData(), - pSinoMem->getDataConst(), + CFloat32ProjectionData3DMemory* pSinoMemory = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram); + ASTRA_ASSERT(pSinoMemory); + CFloat32VolumeData3DMemory* pReconMemory = dynamic_cast<CFloat32VolumeData3DMemory*>(m_pReconstruction); + ASTRA_ASSERT(pReconMemory); + astraCudaBP_SIRTWeighted(pReconMemory->getData(), + pSinoMemory->getDataConst(), &volgeom, projgeom, m_iGPUIndex, m_iVoxelSuperSampling); } else { diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp index d02db6d..d503351 100644 --- a/src/CudaFDKAlgorithm3D.cpp +++ b/src/CudaFDKAlgorithm3D.cpp @@ -59,8 +59,8 @@ CCudaFDKAlgorithm3D::CCudaFDKAlgorithm3D() //---------------------------------------------------------------------------------------- // Constructor with initialization CCudaFDKAlgorithm3D::CCudaFDKAlgorithm3D(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pProjectionData, - CFloat32VolumeData3DMemory* _pReconstruction) + CFloat32ProjectionData3D* _pProjectionData, + CFloat32VolumeData3D* _pReconstruction) { _clear(); initialize(_pProjector, _pProjectionData, _pReconstruction); @@ -179,8 +179,8 @@ bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg) //---------------------------------------------------------------------------------------- // Initialize - C++ bool CCudaFDKAlgorithm3D::initialize(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pSinogram, - CFloat32VolumeData3DMemory* _pReconstruction) + CFloat32ProjectionData3D* _pSinogram, + CFloat32VolumeData3D* _pReconstruction) { // if already initialized, clear first if (m_bIsInitialized) { @@ -225,9 +225,9 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(conegeom); - CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram); + CFloat32ProjectionData3D* pSinoMem = dynamic_cast<CFloat32ProjectionData3D*>(m_pSinogram); ASTRA_ASSERT(pSinoMem); - CFloat32VolumeData3DMemory* pReconMem = dynamic_cast<CFloat32VolumeData3DMemory*>(m_pReconstruction); + CFloat32VolumeData3D* pReconMem = dynamic_cast<CFloat32VolumeData3D*>(m_pReconstruction); ASTRA_ASSERT(pReconMem); const float *filter = NULL; diff --git a/src/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp index 6783093..ce808eb 100644 --- a/src/CudaForwardProjectionAlgorithm3D.cpp +++ b/src/CudaForwardProjectionAlgorithm3D.cpp @@ -101,14 +101,14 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg) node = _cfg.self.getSingleNode("ProjectionDataId"); ASTRA_CONFIG_CHECK(node, "CudaForwardProjection3D", "No ProjectionDataId tag specified."); id = node.getContentInt(); - m_pProjections = dynamic_cast<CFloat32ProjectionData3DMemory*>(CData3DManager::getSingleton().get(id)); + m_pProjections = dynamic_cast<CFloat32ProjectionData3D*>(CData3DManager::getSingleton().get(id)); CC.markNodeParsed("ProjectionDataId"); // reconstruction data node = _cfg.self.getSingleNode("VolumeDataId"); ASTRA_CONFIG_CHECK(node, "CudaForwardProjection3D", "No VolumeDataId tag specified."); id = node.getContentInt(); - m_pVolume = dynamic_cast<CFloat32VolumeData3DMemory*>(CData3DManager::getSingleton().get(id)); + m_pVolume = dynamic_cast<CFloat32VolumeData3D*>(CData3DManager::getSingleton().get(id)); CC.markNodeParsed("VolumeDataId"); // optional: projector @@ -140,8 +140,8 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg) bool CCudaForwardProjectionAlgorithm3D::initialize(CProjector3D* _pProjector, - CFloat32ProjectionData3DMemory* _pProjections, - CFloat32VolumeData3DMemory* _pVolume, + CFloat32ProjectionData3D* _pProjections, + CFloat32VolumeData3D* _pVolume, int _iGPUindex, int _iDetectorSuperSampling) { m_pProjector = _pProjector; diff --git a/src/Float32Data3D.cpp b/src/Float32Data3D.cpp index cad1f18..cc824bd 100644 --- a/src/Float32Data3D.cpp +++ b/src/Float32Data3D.cpp @@ -28,6 +28,10 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. #include "astra/Float32Data3D.h" #include <sstream> +#ifdef ASTRA_CUDA +#include "../../cuda/3d/mem3d.h" +#endif + using namespace std; namespace astra { @@ -60,7 +64,5 @@ std::string CFloat32Data3D::description() const if (getType() == CFloat32Data3D::VOLUME) res << " volume data \t"; return res.str(); } -//---------------------------------------------------------------------------------------- - } // end namespace astra diff --git a/src/Float32Data3DGPU.cpp b/src/Float32Data3DGPU.cpp new file mode 100644 index 0000000..cd9c4ad --- /dev/null +++ b/src/Float32Data3DGPU.cpp @@ -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 "astra/Float32Data3DGPU.h" + +namespace astra { + +//---------------------------------------------------------------------------------------- +// Default constructor. +CFloat32Data3DGPU::CFloat32Data3DGPU() +{ + _clear(); + m_bInitialized = false; +} + +//---------------------------------------------------------------------------------------- +// Destructor. +CFloat32Data3DGPU::~CFloat32Data3DGPU() +{ + if (m_bInitialized) + { + _unInit(); + } +} + +//---------------------------------------------------------------------------------------- +// Initializes an instance of the CFloat32Data3DGPU class with pre-allocated memory +bool CFloat32Data3DGPU::_initialize(int _iWidth, int _iHeight, int _iDepth, astraCUDA3d::MemHandle3D _hnd) +{ + // basic checks + ASTRA_ASSERT(_iWidth > 0); + ASTRA_ASSERT(_iHeight > 0); + ASTRA_ASSERT(_iDepth > 0); + //ASTRA_ASSERT(_pCustomMemory != NULL); + + if (m_bInitialized) { + _unInit(); + } + + // calculate size + m_iWidth = _iWidth; + m_iHeight = _iHeight; + m_iDepth = _iDepth; + m_iSize = (size_t)m_iWidth * m_iHeight * m_iDepth; + + m_hnd = _hnd; + + // initialization complete + return true; +} +//---------------------------------------------------------------------------------------- +// Clear all member variables, setting all numeric variables to 0 and all pointers to NULL. +void CFloat32Data3DGPU::_clear() +{ + m_iWidth = 0; + m_iHeight = 0; + m_iDepth = 0; + m_iSize = 0; + + m_hnd.d.reset(); +} + +//---------------------------------------------------------------------------------------- +// Un-initialize the object, bringing it back in the unitialized state. +void CFloat32Data3DGPU::_unInit() +{ + ASTRA_ASSERT(m_bInitialized); + + _clear(); + m_bInitialized = false; +} + +} // end namespace astra diff --git a/src/Float32Data3DMemory.cpp b/src/Float32Data3DMemory.cpp index 95b1859..5c5c310 100644 --- a/src/Float32Data3DMemory.cpp +++ b/src/Float32Data3DMemory.cpp @@ -27,6 +27,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. #include "astra/Float32Data3DMemory.h" #include <iostream> +#include <cstdlib> namespace astra { @@ -75,10 +76,6 @@ bool CFloat32Data3DMemory::_initialize(int _iWidth, int _iHeight, int _iDepth) m_pCustomMemory = 0; _allocateData(); - // set minmax to default values - m_fGlobalMin = 0.0; - m_fGlobalMax = 0.0; - // initialization complete return true; @@ -166,7 +163,7 @@ bool CFloat32Data3DMemory::_initialize(int _iWidth, int _iHeight, int _iDepth, C ASTRA_ASSERT(_iWidth > 0); ASTRA_ASSERT(_iHeight > 0); ASTRA_ASSERT(_iDepth > 0); - ASTRA_ASSERT(_pCustomMemory != NULL); + //ASTRA_ASSERT(_pCustomMemory != NULL); if (m_bInitialized) { _unInit(); @@ -272,9 +269,6 @@ void CFloat32Data3DMemory::_clear() m_ppfDataRowInd = NULL; m_pppfDataSliceInd = NULL; m_pCustomMemory = NULL; - - //m_fGlobalMin = 0.0f; - //m_fGlobalMax = 0.0f; } //---------------------------------------------------------------------------------------- @@ -289,37 +283,6 @@ void CFloat32Data3DMemory::_unInit() } //---------------------------------------------------------------------------------------- -// Update data statistics, such as minimum and maximum value, after the data has been modified. -void CFloat32Data3DMemory::updateStatistics() -{ - _computeGlobalMinMax(); -} - -//---------------------------------------------------------------------------------------- -// Find the minimum and maximum data value. -void CFloat32Data3DMemory::_computeGlobalMinMax() -{ - // basic checks - ASTRA_ASSERT(m_bInitialized); - ASTRA_ASSERT(m_pfData != NULL); - ASTRA_ASSERT(m_iSize > 0); - - // initial values - m_fGlobalMin = m_pfData[0]; - m_fGlobalMax = m_pfData[0]; - - // loop - size_t i; - float32 v; - for (i = 0; i < m_iSize; ++i) - { - v = m_pfData[i]; - if (v < m_fGlobalMin) m_fGlobalMin = v; - if (v > m_fGlobalMax) m_fGlobalMax = v; - } -} - -//---------------------------------------------------------------------------------------- // Copy the data block pointed to by _pfData to the data block pointed to by m_pfData. void CFloat32Data3DMemory::copyData(const float32* _pfData, size_t _iSize) { diff --git a/src/Float32ProjectionData3D.cpp b/src/Float32ProjectionData3D.cpp index c751c58..ba41e10 100644 --- a/src/Float32ProjectionData3D.cpp +++ b/src/Float32ProjectionData3D.cpp @@ -47,232 +47,6 @@ CFloat32ProjectionData3D::~CFloat32ProjectionData3D() { m_pGeometry = 0; } -CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const CFloat32ProjectionData3D& _data) -{ - CProjectionGeometry3D * pThisGeometry = getGeometry(); - - int iProjectionCount = pThisGeometry->getProjectionCount(); - int iDetectorCount = pThisGeometry->getDetectorTotCount(); -#ifdef _DEBUG - CProjectionGeometry3D * pDataGeometry = _data.getGeometry(); - int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount(); - - ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount()); - ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount); -#endif - - for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; - - fThisValue += fDataValue; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnProjection(iProjectionIndex, pThisProjection); - - delete pThisProjection; - delete pDataProjection; - } - - return *this; -} - -CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const CFloat32ProjectionData3D& _data) -{ - CProjectionGeometry3D * pThisGeometry = getGeometry(); - - int iProjectionCount = pThisGeometry->getProjectionCount(); - int iDetectorCount = pThisGeometry->getDetectorTotCount(); -#ifdef _DEBUG - CProjectionGeometry3D * pDataGeometry = _data.getGeometry(); - int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount(); - - ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount()); - ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount); -#endif - - for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; - - fThisValue -= fDataValue; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnProjection(iProjectionIndex, pThisProjection); - - delete pThisProjection; - delete pDataProjection; - } - - return *this; -} - -CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const CFloat32ProjectionData3D& _data) -{ - CProjectionGeometry3D * pThisGeometry = getGeometry(); - - int iProjectionCount = pThisGeometry->getProjectionCount(); - int iDetectorCount = pThisGeometry->getDetectorTotCount(); -#ifdef _DEBUG - CProjectionGeometry3D * pDataGeometry = _data.getGeometry(); - int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount(); - - ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount()); - ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount); -#endif - - for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; - - fThisValue *= fDataValue; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnProjection(iProjectionIndex, pThisProjection); - - delete pThisProjection; - delete pDataProjection; - } - - return *this; -} - -CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const float32& _fScalar) -{ - CProjectionGeometry3D * pThisGeometry = getGeometry(); - - int iProjectionCount = pThisGeometry->getProjectionCount(); - int iDetectorCount = pThisGeometry->getDetectorTotCount(); - - for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - - fThisValue *= _fScalar; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnProjection(iProjectionIndex, pThisProjection); - - delete pThisProjection; - } - - return *this; -} - -CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator/=(const float32& _fScalar) -{ - CProjectionGeometry3D * pThisGeometry = getGeometry(); - - int iProjectionCount = pThisGeometry->getProjectionCount(); - int iDetectorCount = pThisGeometry->getDetectorTotCount(); - - for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - - fThisValue /= _fScalar; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnProjection(iProjectionIndex, pThisProjection); - - delete pThisProjection; - } - - return *this; -} - -CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const float32& _fScalar) -{ - CProjectionGeometry3D * pThisGeometry = getGeometry(); - - int iProjectionCount = pThisGeometry->getProjectionCount(); - int iDetectorCount = pThisGeometry->getDetectorTotCount(); - - for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - - fThisValue += _fScalar; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnProjection(iProjectionIndex, pThisProjection); - - delete pThisProjection; - } - - return *this; -} - -CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const float32& _fScalar) -{ - CProjectionGeometry3D * pThisGeometry = getGeometry(); - - int iProjectionCount = pThisGeometry->getProjectionCount(); - int iDetectorCount = pThisGeometry->getDetectorTotCount(); - - for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - - fThisValue -= _fScalar; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnProjection(iProjectionIndex, pThisProjection); - - delete pThisProjection; - } - - return *this; -} - void CFloat32ProjectionData3D::changeGeometry(CProjectionGeometry3D* _pGeometry) { if (!m_bInitialized) return; diff --git a/src/Float32ProjectionData3DGPU.cpp b/src/Float32ProjectionData3DGPU.cpp new file mode 100644 index 0000000..0e063d6 --- /dev/null +++ b/src/Float32ProjectionData3DGPU.cpp @@ -0,0 +1,71 @@ +/* +----------------------------------------------------------------------- +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/Float32ProjectionData3DGPU.h" + +using namespace std; + +namespace astra +{ + +//---------------------------------------------------------------------------------------- +// Default constructor +CFloat32ProjectionData3DGPU::CFloat32ProjectionData3DGPU() : + CFloat32Data3DGPU() +{ + m_pGeometry = NULL; + m_bInitialized = false; +} + +//---------------------------------------------------------------------------------------- +// Create an instance of the CFloat32ProjectionData2D class with pre-allocated data +CFloat32ProjectionData3DGPU::CFloat32ProjectionData3DGPU(CProjectionGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd) +{ + m_bInitialized = false; + m_bInitialized = initialize(_pGeometry, _hnd); +} + + +//---------------------------------------------------------------------------------------- +// Destructor +CFloat32ProjectionData3DGPU::~CFloat32ProjectionData3DGPU() +{ + delete m_pGeometry; + m_pGeometry = 0; +} + +//---------------------------------------------------------------------------------------- +// Initialization +bool CFloat32ProjectionData3DGPU::initialize(CProjectionGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd) +{ + m_pGeometry = _pGeometry->clone(); + m_bInitialized = _initialize(m_pGeometry->getDetectorColCount(), m_pGeometry->getProjectionCount(), m_pGeometry->getDetectorRowCount(), _hnd); + return m_bInitialized; +} + +} // end namespace astra diff --git a/src/Float32ProjectionData3DMemory.cpp b/src/Float32ProjectionData3DMemory.cpp index f42e216..69033d5 100644 --- a/src/Float32ProjectionData3DMemory.cpp +++ b/src/Float32ProjectionData3DMemory.cpp @@ -114,7 +114,7 @@ bool CFloat32ProjectionData3DMemory::initialize(CProjectionGeometry3D* _pGeometr //---------------------------------------------------------------------------------------- // Initialization -bool CFloat32ProjectionData3DMemory::initialize(CProjectionGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) +bool CFloat32ProjectionData3DMemory::initialize(CProjectionGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) { m_pGeometry = _pGeometry->clone(); m_bInitialized = _initialize(m_pGeometry->getDetectorColCount(), m_pGeometry->getProjectionCount(), m_pGeometry->getDetectorRowCount(), _pCustomMemory); @@ -131,102 +131,6 @@ CFloat32ProjectionData3DMemory::~CFloat32ProjectionData3DMemory() } //---------------------------------------------------------------------------------------- -// Fetch a projection -CFloat32VolumeData2D* CFloat32ProjectionData3DMemory::fetchProjection(int _iProjectionNr) const -{ - // fetch slice of the geometry - CVolumeGeometry2D volGeom(m_pGeometry->getDetectorColCount(), m_pGeometry->getDetectorRowCount()); - // create new volume data - CFloat32VolumeData2D* res = new CFloat32VolumeData2D(&volGeom); - // copy data - int row, col; - for (row = 0; row < m_pGeometry->getDetectorRowCount(); ++row) { - for (col = 0; col < m_pGeometry->getDetectorColCount(); ++col) { - res->getData()[row*m_pGeometry->getDetectorColCount() + col] = - m_pfData[_iProjectionNr * m_pGeometry->getDetectorColCount() + m_pGeometry->getDetectorColCount()* m_pGeometry->getProjectionCount() * row + col]; - } - } - // return - return res; -} - -//---------------------------------------------------------------------------------------- -// Return a projection -void CFloat32ProjectionData3DMemory::returnProjection(int _iProjectionNr, CFloat32VolumeData2D* _pProjection) -{ - /// TODO: check geometry - // copy data - int row, col; - for (row = 0; row < m_pGeometry->getDetectorRowCount(); ++row) { - for (col = 0; col < m_pGeometry->getDetectorColCount(); ++col) { - m_pfData[_iProjectionNr * m_pGeometry->getDetectorColCount() + m_pGeometry->getDetectorColCount()* m_pGeometry->getProjectionCount() * row + col] = - _pProjection->getData()[row*m_pGeometry->getDetectorColCount() + col]; - } - } -} - -//---------------------------------------------------------------------------------------- -// Fetch a sinogram -CFloat32ProjectionData2D* CFloat32ProjectionData3DMemory::fetchSinogram(int _iSliceNr) const -{ - CParallelProjectionGeometry3D * pParallelProjGeo = (CParallelProjectionGeometry3D *)m_pGeometry; - CParallelProjectionGeometry2D * pProjGeo2D = pParallelProjGeo->createProjectionGeometry2D(); - - // create new projection data - CFloat32ProjectionData2D* res = new CFloat32ProjectionData2D(pProjGeo2D); - // copy data - int row, col; - - int iDetectorColumnCount = m_pGeometry->getDetectorColCount(); - int iProjectionAngleCount = m_pGeometry->getProjectionCount(); - - for (row = 0; row < m_pGeometry->getProjectionCount(); ++row) { - for (col = 0; col < m_pGeometry->getDetectorColCount(); ++col) - { - int iTargetIndex = row * iDetectorColumnCount + col; - int iSourceIndex = _iSliceNr * iDetectorColumnCount * iProjectionAngleCount + row * iDetectorColumnCount + col; - - float32 fStoredValue = m_pfData[iSourceIndex]; - - res->getData()[iTargetIndex] = fStoredValue; - } - } - - delete pProjGeo2D; - - // return - return res; -} - -//---------------------------------------------------------------------------------------- -// Return a sinogram -void CFloat32ProjectionData3DMemory::returnSinogram(int _iSliceNr, CFloat32ProjectionData2D* _pSinogram2D) -{ - /// TODO: check geometry - // copy data - int row, col; - for (row = 0; row < m_pGeometry->getProjectionCount(); ++row) { - for (col = 0; col < m_pGeometry->getDetectorColCount(); ++col) { - m_pfData[_iSliceNr*m_pGeometry->getDetectorColCount()*m_pGeometry->getProjectionCount() + row*m_pGeometry->getDetectorColCount() + col] = - _pSinogram2D->getData()[row*m_pGeometry->getDetectorColCount() + col]; - } - } -} - -//---------------------------------------------------------------------------------------- -// Returns a specific value -float32 CFloat32ProjectionData3DMemory::getDetectorValue(int _iIndex) -{ - return m_pfData[_iIndex]; -} - -//---------------------------------------------------------------------------------------- -// Sets a specific value -void CFloat32ProjectionData3DMemory::setDetectorValue(int _iIndex, float32 _fValue) -{ - m_pfData[_iIndex] = _fValue; -} -//---------------------------------------------------------------------------------------- CFloat32ProjectionData3DMemory& CFloat32ProjectionData3DMemory::operator=(const CFloat32ProjectionData3DMemory& _dataIn) { diff --git a/src/Float32VolumeData3D.cpp b/src/Float32VolumeData3D.cpp index bf7eced..96e6940 100644 --- a/src/Float32VolumeData3D.cpp +++ b/src/Float32VolumeData3D.cpp @@ -43,228 +43,6 @@ CFloat32VolumeData3D::~CFloat32VolumeData3D() { } -CFloat32VolumeData3D& CFloat32VolumeData3D::operator+=(const CFloat32VolumeData3D& _data) -{ - CVolumeGeometry3D * pThisGeometry = getGeometry(); - - int iSliceCount = pThisGeometry->getGridSliceCount(); -#ifdef _DEBUG - CVolumeGeometry3D * pDataGeometry = _data.getGeometry(); - int iThisSlicePixelCount = pThisGeometry->getGridRowCount() * pThisGeometry->getGridColCount(); - int iDataSlicePixelCount = pDataGeometry->getGridRowCount() * pDataGeometry->getGridColCount(); - - ASTRA_ASSERT(iSliceCount == pDataGeometry->getGridSliceCount()); - ASTRA_ASSERT(iThisSlicePixelCount == iDataSlicePixelCount); -#endif - - for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex); - CFloat32VolumeData2D * pDataProjection = _data.fetchSliceZ(iSliceIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; - - fThisValue += fDataValue; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnSliceZ(iSliceIndex, pThisProjection); - - delete pThisProjection; - delete pDataProjection; - } - - return *this; -} - -CFloat32VolumeData3D& CFloat32VolumeData3D::operator-=(const CFloat32VolumeData3D& _data) -{ - CVolumeGeometry3D * pThisGeometry = getGeometry(); - - int iSliceCount = pThisGeometry->getGridSliceCount(); -#ifdef _DEBUG - CVolumeGeometry3D * pDataGeometry = _data.getGeometry(); - int iThisSlicePixelCount = pThisGeometry->getGridRowCount() * pThisGeometry->getGridColCount(); - int iDataSlicePixelCount = pDataGeometry->getGridRowCount() * pDataGeometry->getGridColCount(); - - ASTRA_ASSERT(iSliceCount == pDataGeometry->getGridSliceCount()); - ASTRA_ASSERT(iThisSlicePixelCount == iDataSlicePixelCount); -#endif - - for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex); - CFloat32VolumeData2D * pDataProjection = _data.fetchSliceZ(iSliceIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; - - fThisValue -= fDataValue; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnSliceZ(iSliceIndex, pThisProjection); - - delete pThisProjection; - delete pDataProjection; - } - - return *this; -} - -CFloat32VolumeData3D& CFloat32VolumeData3D::operator*=(const CFloat32VolumeData3D& _data) -{ - CVolumeGeometry3D * pThisGeometry = getGeometry(); - - int iSliceCount = pThisGeometry->getGridSliceCount(); -#ifdef _DEBUG - CVolumeGeometry3D * pDataGeometry = _data.getGeometry(); - int iThisSlicePixelCount = pThisGeometry->getGridRowCount() * pThisGeometry->getGridColCount(); - int iDataSlicePixelCount = pDataGeometry->getGridRowCount() * pDataGeometry->getGridColCount(); - - ASTRA_ASSERT(iSliceCount == pDataGeometry->getGridSliceCount()); - ASTRA_ASSERT(iThisSlicePixelCount == iDataSlicePixelCount); -#endif - - for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex); - CFloat32VolumeData2D * pDataProjection = _data.fetchSliceZ(iSliceIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; - - fThisValue *= fDataValue; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnSliceZ(iSliceIndex, pThisProjection); - - delete pThisProjection; - delete pDataProjection; - } - - return *this; -} - -CFloat32VolumeData3D& CFloat32VolumeData3D::operator*=(const float32& _fScalar) -{ - CVolumeGeometry3D * pThisGeometry = getGeometry(); - - int iSliceCount = pThisGeometry->getGridSliceCount(); - - for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - - fThisValue *= _fScalar; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnSliceZ(iSliceIndex, pThisProjection); - - delete pThisProjection; - } - - return *this; -} - -CFloat32VolumeData3D& CFloat32VolumeData3D::operator/=(const float32& _fScalar) -{ - CVolumeGeometry3D * pThisGeometry = getGeometry(); - - int iSliceCount = pThisGeometry->getGridSliceCount(); - - for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - - fThisValue /= _fScalar; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnSliceZ(iSliceIndex, pThisProjection); - - delete pThisProjection; - } - - return *this; -} - -CFloat32VolumeData3D& CFloat32VolumeData3D::operator+=(const float32& _fScalar) -{ - CVolumeGeometry3D * pThisGeometry = getGeometry(); - - int iSliceCount = pThisGeometry->getGridSliceCount(); - - for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - - fThisValue += _fScalar; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnSliceZ(iSliceIndex, pThisProjection); - - delete pThisProjection; - } - - return *this; -} - -CFloat32VolumeData3D& CFloat32VolumeData3D::operator-=(const float32& _fScalar) -{ - CVolumeGeometry3D * pThisGeometry = getGeometry(); - - int iSliceCount = pThisGeometry->getGridSliceCount(); - - for(int iSliceIndex = 0; iSliceIndex < iSliceCount; iSliceIndex++) - { - CFloat32VolumeData2D * pThisProjection = fetchSliceZ(iSliceIndex); - - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) - { - float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; - - fThisValue -= _fScalar; - - pThisProjection->getData()[iDetectorIndex] = fThisValue; - } - - returnSliceZ(iSliceIndex, pThisProjection); - - delete pThisProjection; - } - - return *this; -} - void CFloat32VolumeData3D::changeGeometry(CVolumeGeometry3D* _pGeometry) { if (!m_bInitialized) return; diff --git a/src/Float32VolumeData3DGPU.cpp b/src/Float32VolumeData3DGPU.cpp new file mode 100644 index 0000000..82a222f --- /dev/null +++ b/src/Float32VolumeData3DGPU.cpp @@ -0,0 +1,71 @@ +/* +----------------------------------------------------------------------- +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/Float32VolumeData3DGPU.h" + +using namespace std; + +namespace astra +{ + +//---------------------------------------------------------------------------------------- +// Default constructor +CFloat32VolumeData3DGPU::CFloat32VolumeData3DGPU() : + CFloat32Data3DGPU() +{ + m_pGeometry = NULL; + m_bInitialized = false; +} + +//---------------------------------------------------------------------------------------- +// Create an instance of the CFloat32VolumeData2D class with pre-allocated data +CFloat32VolumeData3DGPU::CFloat32VolumeData3DGPU(CVolumeGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd) +{ + m_bInitialized = false; + m_bInitialized = initialize(_pGeometry, _hnd); +} + + +//---------------------------------------------------------------------------------------- +// Destructor +CFloat32VolumeData3DGPU::~CFloat32VolumeData3DGPU() +{ + delete m_pGeometry; + m_pGeometry = 0; +} + +//---------------------------------------------------------------------------------------- +// Initialization +bool CFloat32VolumeData3DGPU::initialize(CVolumeGeometry3D* _pGeometry, astraCUDA3d::MemHandle3D _hnd) +{ + m_pGeometry = _pGeometry->clone(); + m_bInitialized = _initialize(m_pGeometry->getGridColCount(), m_pGeometry->getGridRowCount(), m_pGeometry->getGridSliceCount(), _hnd); + return m_bInitialized; +} + +} // end namespace astra diff --git a/src/Float32VolumeData3DMemory.cpp b/src/Float32VolumeData3DMemory.cpp index ba39eaa..27ae87b 100644 --- a/src/Float32VolumeData3DMemory.cpp +++ b/src/Float32VolumeData3DMemory.cpp @@ -114,103 +114,13 @@ bool CFloat32VolumeData3DMemory::initialize(CVolumeGeometry3D* _pGeometry, float } //---------------------------------------------------------------------------------------- // Initialization -bool CFloat32VolumeData3DMemory::initialize(CVolumeGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) +bool CFloat32VolumeData3DMemory::initialize(CVolumeGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) { m_pGeometry = _pGeometry->clone(); m_bInitialized = _initialize(m_pGeometry->getGridColCount(), m_pGeometry->getGridRowCount(), m_pGeometry->getGridSliceCount(), _pCustomMemory); return m_bInitialized; } - -//---------------------------------------------------------------------------------------- -// Fetch a slice -CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceZ(int _iSliceIndex) const -{ - // fetch slice of the geometry - int iRowCount = m_pGeometry->getGridRowCount(); - int iColumnCount = m_pGeometry->getGridColCount(); - CVolumeGeometry2D volGeom(iColumnCount, iRowCount); - - // create new volume data - CFloat32VolumeData2D* res = new CFloat32VolumeData2D(&volGeom); - - // copy data - float * pfTargetData = res->getData(); - for(int iRowIndex = 0; iRowIndex < iRowCount; iRowIndex++) - { - for(int iColumnIndex = 0; iColumnIndex < iColumnCount; iColumnIndex++) - { - int iSourceIndex = _iSliceIndex * iColumnCount * iRowCount + iRowIndex * iColumnCount + iColumnIndex; - int iTargetIndex = iRowIndex * iColumnCount + iColumnIndex; - float fStoredValue = m_pfData[iSourceIndex]; - pfTargetData[iTargetIndex] = fStoredValue; - } - } - // return - return res; -} - -//---------------------------------------------------------------------------------------- -// Return a slice -void CFloat32VolumeData3DMemory::returnSliceZ(int _iSliceIndex, CFloat32VolumeData2D * _pSlice) -{ - int iRowCount = _pSlice->getGeometry()->getGridRowCount(); - int iColumnCount = _pSlice->getGeometry()->getGridColCount(); - - assert(iRowCount == m_pGeometry->getGridRowCount()); - assert(iColumnCount == m_pGeometry->getGridColCount()); - - for(int iRowIndex = 0; iRowIndex < iRowCount; iRowIndex++) - { - for(int iColumnIndex = 0; iColumnIndex < iColumnCount; iColumnIndex++) - { - int iSourceIndex = iRowIndex * iColumnCount + iColumnIndex; - int iTargetIndex = _iSliceIndex * iColumnCount * iRowCount + iRowIndex * iColumnCount + iColumnIndex; - float fStoredValue = _pSlice->getDataConst()[iSourceIndex]; - m_pfData[iTargetIndex] = fStoredValue; - } - } -} - -CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceX(int _iColumnIndex) const -{ - // TODO: - assert(false); - return NULL; -} - -CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceY(int _iRowIndex) const -{ - // TODO: - assert(false); - return NULL; -} - -void CFloat32VolumeData3DMemory::returnSliceX(int _iColumnIndex, CFloat32VolumeData2D * _pSliceData) -{ - // TODO: - assert(false); -} - -void CFloat32VolumeData3DMemory::returnSliceY(int _iRowIndex, CFloat32VolumeData2D * _pSliceData) -{ - // TODO: - assert(false); -} - -//---------------------------------------------------------------------------------------- -// Returns a specific value -float32 CFloat32VolumeData3DMemory::getVoxelValue(int _iIndex) -{ - return m_pfData[_iIndex]; -} - -//---------------------------------------------------------------------------------------- -// Sets a specific value -void CFloat32VolumeData3DMemory::setVoxelValue(int _iIndex, float32 _fValue) -{ - m_pfData[_iIndex] = _fValue; -} //---------------------------------------------------------------------------------------- CFloat32VolumeData3DMemory& CFloat32VolumeData3DMemory::operator=(const CFloat32VolumeData3DMemory& _dataIn) |