Merge pull request #52 from andyschwarzl/gpu-arrays

Adapt gridding operations to support gpu arrays. Fixes #50. Fixes #20.
andyschwarzl · Feb 17, 2015 · a71702b · a71702b
2 parents dc980e5 + 57fa3ba
commit a71702b
Show file tree

Hide file tree

Showing 34 changed files with 759 additions and 312 deletions.
diff --git a/.clang-format b/.clang-format
@@ -0,0 +1,48 @@
+---
+# BasedOnStyle:  LLVM
+AccessModifierOffset: -1
+ConstructorInitializerIndentWidth: 2
+AlignEscapedNewlinesLeft: false
+AlignTrailingComments: true
+AllowShortFunctionsOnASingleLine: false
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+AlwaysBreakTemplateDeclarations: false
+AlwaysBreakBeforeMultilineStrings: false
+BreakBeforeBinaryOperators: false
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BinPackParameters: true
+ColumnLimit:     80
+ConstructorInitializerAllOnOneLineOrOnePerLine: false
+DerivePointerBinding: false
+ExperimentalAutoDetectBinPacking: false
+IndentCaseLabels: false
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: None
+ObjCSpaceBeforeProtocolList: true
+PenaltyBreakBeforeFirstCallParameter: 19
+PenaltyBreakComment: 60
+PenaltyBreakString: 1000
+PenaltyBreakFirstLessLess: 120
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 60
+PointerBindsToType: false
+SpacesBeforeTrailingComments: 2
+Cpp11BracedListStyle: false
+Standard:        Cpp11
+IndentWidth:     2
+TabWidth:        2
+UseTab:          Never
+BreakBeforeBraces: Allman
+IndentFunctionDeclarationAfterType: false
+SpacesInParentheses: false
+SpacesInAngles:  false
+SpaceInEmptyParentheses: false
+SpacesInCStyleCastParentheses: false
+SpaceAfterControlStatementKeyword: true
+SpaceBeforeAssignmentOperators: true
+ContinuationIndentWidth: 4
+...
+
diff --git a/.gitignore b/.gitignore
@@ -43,5 +43,7 @@ inc/*config.hpp
 *.mat
 tmp/
 *.swp
+*.hdr
+*.cfl
 
 matlab/demo/results/
diff --git a/CUDA/.ycm_extra_conf.py → .ycm_extra_conf.py b/CUDA/.ycm_extra_conf.py → .ycm_extra_conf.py
diff --git a/CUDA/CMakeLists.txt b/CUDA/CMakeLists.txt
@@ -78,6 +78,9 @@ endif(GPU_DOUBLE_PREC)
 
 MESSAGE(STATUS "setting NVCC FLAGS to: ${CUDA_NVCC_FLAGS}")
 
+# TODO build type dependent
+SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -g")
+
 #LIB and MEX-file names
 SET(GRID_LIB_NAME "gpuNUFFT${PREC_SUFFIX}")
 SET(GRID_LIB_ATM_NAME "gpuNUFFT_ATM${PREC_SUFFIX}")

diff --git a/CUDA/cmake_modules/FindGPUNUFFT.cmake b/CUDA/cmake_modules/FindGPUNUFFT.cmake
@@ -0,0 +1,24 @@
+INCLUDE(FindPackageHandleStandardArgs)
+
+IF (NOT DEFINED GPUNUFFT_ROOT_DIR)
+  MESSAGE(FATAL_ERROR "GPUNUFFT_ROOT_DIR not defined!")
+ELSE()
+  MESSAGE("Searching for GPUNUFFT in ${GPUNUFFT_ROOT_DIR}")
+
+  SET(GPUNUFFT_IncludeSearchPaths ${GPUNUFFT_ROOT_DIR}/CUDA/inc ${GPUNUFFT_ROOT_DIR}/inc) #/usr/local/include
+  SET(GPUNUFFT_LibrarySearchPaths ${GPUNUFFT_ROOT_DIR}/CUDA/bin ${GPUNUFFT_ROOT_DIR}/bin) #/usr/local/lib
+
+  FIND_PATH(GPUNUFFT_INCLUDE_DIRS NAMES gpuNUFFT_operator_factory.hpp PATHS ${GPUNUFFT_IncludeSearchPaths})
+
+  FIND_LIBRARY(GPUNUFFT_LIBRARIES NAMES gpuNUFFT gpuNUFFT_ATM_f gpuNUFFT_f  PATHS ${GPUNUFFT_LibrarySearchPaths})
+
+  #Handle the REQUIRED argument and set the < UPPERCASED_NAME > _FOUND variable
+  FIND_PACKAGE_HANDLE_STANDARD_ARGS(GPUNUFFT "Could NOT find GPUNUFFT." GPUNUFFT_LIBRARIES GPUNUFFT_INCLUDE_DIRS)
+  IF(GPUNUFFT_FOUND)
+    FIND_PACKAGE_MESSAGE(
+        GPUNUFFT      
+        "Found GPUNUFFT ${GPUNUFFT_LIBRARIES}"
+        "[${GPUNUFFT_LIBRARIES}][${GPUNUFFT_INCLUDE_DIRS}]")
+    ENDIF(GPUNUFFT_FOUND)
+  MARK_AS_ADVANCED(GPUNUFFT_INCLUDE_DIRS GPUNUFFT_LIBRARIES)
+ENDIF()
diff --git a/CUDA/cmake_modules/FindGpuNUFFT.cmake b/CUDA/cmake_modules/FindGpuNUFFT.cmake
diff --git a/CUDA/inc/balanced_gpuNUFFT_operator.hpp b/CUDA/inc/balanced_gpuNUFFT_operator.hpp
@@ -21,47 +21,47 @@ namespace gpuNUFFT
   public:
 
     BalancedGpuNUFFTOperator(IndType kernelWidth, IndType sectorWidth, DType osf, Dimensions imgDims): 
-    GpuNUFFTOperator(kernelWidth,sectorWidth,osf,imgDims,true,BALANCED)
+      GpuNUFFTOperator(kernelWidth,sectorWidth,osf,imgDims,true,BALANCED)
     {
     }
 
-    ~BalancedGpuNUFFTOperator()
+    virtual ~BalancedGpuNUFFTOperator()
     {
     }
 
     Array<IndType2>  getSectorProcessingOrder(){return this->sectorProcessingOrder;}
-    void setSectorProcessingOrder(Array<IndType2> sectorProcessingOrder)	{this->sectorProcessingOrder = sectorProcessingOrder;}
+    void setSectorProcessingOrder(Array<IndType2> sectorProcessingOrder)  {this->sectorProcessingOrder = sectorProcessingOrder;}
 
     // OPERATIONS
     void performGpuNUFFTAdj(Array<DType2> kspaceData, Array<CufftType>& imgData, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
     void performForwardGpuNUFFT(Array<DType2> imgData,Array<CufftType>& kspaceData, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
-        
+
     OperatorType getType() {return gpuNUFFT::BALANCED;}
 
   protected:
-	
+
     // sectorProcessingOrder
     Array<IndType2> sectorProcessingOrder;
-    
+
     IndType2* sector_processing_order_d;
-    
+
     GpuNUFFTInfo* initAndCopyGpuNUFFTInfo();
 
     void adjConvolution(DType2* data_d, 
-      DType* crds_d, 
-      CufftType* gdata_d,
-      DType* kernel_d, 
-      IndType* sectors_d, 
-      IndType* sector_centers_d,
-      gpuNUFFT::GpuNUFFTInfo* gi_host);
-			
-    void forwardConvolution(CufftType*		data_d, 
-      DType*			crds_d, 
-      CufftType*		gdata_d,
-      DType*			kernel_d, 
-      IndType*		sectors_d, 
-      IndType*		sector_centers_d,
-      gpuNUFFT::GpuNUFFTInfo* gi_host);
+        DType* crds_d, 
+        CufftType* gdata_d,
+        DType* kernel_d, 
+        IndType* sectors_d, 
+        IndType* sector_centers_d,
+        gpuNUFFT::GpuNUFFTInfo* gi_host);
+
+    void forwardConvolution(CufftType*    data_d, 
+        DType*      crds_d, 
+        CufftType*    gdata_d,
+        DType*      kernel_d, 
+        IndType*    sectors_d, 
+        IndType*    sector_centers_d,
+        gpuNUFFT::GpuNUFFTInfo* gi_host);
   };
 }
 

diff --git a/CUDA/inc/balanced_texture_gpuNUFFT_operator.hpp b/CUDA/inc/balanced_texture_gpuNUFFT_operator.hpp
@@ -34,13 +34,16 @@ namespace gpuNUFFT
       initKernel();	
     }
 
-    ~BalancedTextureGpuNUFFTOperator()
+    virtual ~BalancedTextureGpuNUFFTOperator()
     {
     }
 
     // OPERATIONS
     void performGpuNUFFTAdj(Array<DType2> kspaceData, Array<CufftType>& imgData, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
+    void performGpuNUFFTAdj(GpuArray<DType2> kspaceData_gpu, GpuArray<CufftType>& imgData_gpu, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
+
     void performForwardGpuNUFFT(Array<DType2> imgData,Array<CufftType>& kspaceData, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
+    void performForwardGpuNUFFT(GpuArray<DType2> imgData_gpu, GpuArray<CufftType>& kspaceData, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
 
     //Getter and Setter for Processing Order
     Array<IndType2>  getSectorProcessingOrder(){return this->sectorProcessingOrder;}

diff --git a/CUDA/inc/cuda_utils.hpp b/CUDA/inc/cuda_utils.hpp
@@ -82,6 +82,18 @@ inline void allocateAndSetMem(TypeName** device_ptr, IndType num_elements,int va
   HANDLE_ERROR(cudaMemset(*device_ptr,value,num_elements*sizeof(TypeName)));
 }
 
+/** \brief CUDA memcpy call to copy data from device ptr to device ptr
+ * 
+ * @param device_ptr_src   source device pointer 
+ * @param device_ptr_dest  destination device pointer
+ * @param num_elements     amount of elements of size TypeName
+ */
+template<typename TypeName>
+inline void copyDeviceToDevice(TypeName* device_ptr_src, TypeName* device_ptr_dest, IndType num_elements)
+{
+  HANDLE_ERROR(cudaMemcpy(device_ptr_dest, device_ptr_src, num_elements*sizeof(TypeName),cudaMemcpyDeviceToDevice ));
+}
+
 /** \brief Copy CUDA memory from device to host
  * 
  * @param device_ptr    device pointer

diff --git a/CUDA/inc/gpuNUFFT_operator.hpp b/CUDA/inc/gpuNUFFT_operator.hpp
@@ -38,58 +38,61 @@ namespace gpuNUFFT
       * @param operatorType Overwrite defalt operator type
       */
     GpuNUFFTOperator(IndType kernelWidth, IndType sectorWidth, DType osf, Dimensions imgDims, bool loadKernel = true, OperatorType operatorType = DEFAULT): 
-        osf(osf), kernelWidth(kernelWidth), sectorWidth(sectorWidth),imgDims(imgDims),operatorType(operatorType), gpuMemAllocated(false), debugTiming(DEBUG)
+      operatorType(operatorType),osf(osf), kernelWidth(kernelWidth), sectorWidth(sectorWidth),imgDims(imgDims), gpuMemAllocated(false), debugTiming(DEBUG),
+      sens_d(NULL),crds_d(NULL),density_comp_d(NULL),deapo_d(NULL),gdata_d(NULL),sector_centers_d(NULL),sectors_d(NULL),
+      data_indices_d(NULL),data_sorted_d(NULL)
     {
       if (loadKernel)
-        initKernel();	
+        initKernel();
 
       sectorDims.width = sectorWidth;
       sectorDims.height = sectorWidth;
       if (imgDims.depth > 0)
         sectorDims.depth = sectorWidth;
     }
 
-    ~GpuNUFFTOperator()
+    virtual ~GpuNUFFTOperator()
     {
       free(this->kernel.data);
+      freeDeviceMemory();
     }
 
     friend class GpuNUFFTOperatorFactory;
 
     // SETTER 
-    void setOsf(DType osf)			{this->osf = osf;}
+    void setOsf(DType osf)      {this->osf = osf;}
 
-    void setKSpaceTraj(Array<DType> kSpaceTraj)				{this->kSpaceTraj = kSpaceTraj;}
-    void setSectorCenters(Array<IndType> sectorCenters)		{this->sectorCenters = sectorCenters;}
-    void setSectorDataCount(Array<IndType> sectorDataCount)	{this->sectorDataCount = sectorDataCount;}
-    void setDataIndices(Array<IndType> dataIndices)			{this->dataIndices = dataIndices;}
-    void setSens(Array<DType2> sens)						{this->sens = sens;}
-    void setDens(Array<DType> dens)							{this->dens = dens;}
+    void setKSpaceTraj(Array<DType> kSpaceTraj)               {this->kSpaceTraj = kSpaceTraj;}
+    void setSectorCenters(Array<IndType> sectorCenters)       {this->sectorCenters = sectorCenters;}
+    void setSectorDataCount(Array<IndType> sectorDataCount)   {this->sectorDataCount = sectorDataCount;}
+    void setDataIndices(Array<IndType> dataIndices)           {this->dataIndices = dataIndices;}
+    void setSens(Array<DType2> sens)                          {this->sens = sens;}
+    void setDens(Array<DType> dens)                           {this->dens = dens;}
 
-    void setImageDims(Dimensions dims)						{this->imgDims = dims;}
-    void setGridSectorDims(Dimensions dims)						{this->gridSectorDims = dims;}
+    void setImageDims(Dimensions dims)                        {this->imgDims = dims;}
+    void setGridSectorDims(Dimensions dims)                   {this->gridSectorDims = dims;}
 
     // GETTER
-    Array<DType>  getKSpaceTraj()	{return this->kSpaceTraj;}
+    Array<DType>  getKSpaceTraj()     {return this->kSpaceTraj;}
 
-    Array<DType2>	getSens()			{return this->sens;}
-    Array<DType>	getDens()			{return this->dens;}
-    Array<DType>    getKernel()			{return this->kernel;}
+    Array<DType2>  getSens()          {return this->sens;}
+    Array<DType>  getDens()           {return this->dens;}
+    Array<DType>    getKernel()       {return this->kernel;}
     Array<IndType>  getSectorDataCount(){return this->sectorDataCount;}
 
-    IndType getKernelWidth()		{return this->kernelWidth;}
-    IndType getSectorWidth()		{return this->sectorWidth;}
+    IndType getKernelWidth()    {return this->kernelWidth;}
+    IndType getSectorWidth()    {return this->sectorWidth;}
 
     Dimensions getImageDims() {return this->imgDims;}
-    Dimensions getGridDims() {return this->imgDims * osf;}
+    Dimensions getGridDims()  {return this->imgDims * osf;}
 
-    Dimensions getGridSectorDims() {return this->gridSectorDims;}
-    Dimensions getSectorDims() {return this->sectorDims;}
+    Dimensions getGridSectorDims()  {return this->gridSectorDims;}
+    Dimensions getSectorDims()      {return this->sectorDims;}
 
     Array<IndType> getSectorCenters() {return this->sectorCenters;}
-    IndType* getSectorCentersData() {return reinterpret_cast<IndType*>(this->sectorCenters.data);}
+    IndType* getSectorCentersData()   {return reinterpret_cast<IndType*>(this->sectorCenters.data);}
 
-    Array<IndType>  getDataIndices()		{return this->dataIndices;}
+    Array<IndType>  getDataIndices()  {return this->dataIndices;}
 
     bool is2DProcessing() {return this->imgDims.depth == 0;}
     bool is3DProcessing() {return !is2DProcessing();}
@@ -130,6 +133,19 @@ namespace gpuNUFFT
       */
     virtual void     performGpuNUFFTAdj(Array<DType2> kspaceData, Array<CufftType>& imgData, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
 
+    /** \brief Perform Adjoint gridding operation on kspaceData already residing in GPU memory 
+     *
+     * This may be the case in iterative reconstructions, when k-Space and image
+     * data is already residing on the GPU.
+     *
+     *
+     * @param k-space data 
+     * @param preallocated image data array
+     * @param Stop gridding operation after gpuNUFFT::GpuNUFFTOutput 
+     * @return Regridded Image
+     */
+    virtual void     performGpuNUFFTAdj(GpuArray<DType2> kspaceData_gpu, GpuArray<CufftType>& imgData_gpu, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
+
     /** \brief Perform Adjoint gridding operation on given kspaceData 
       *
       * @param k-space data 
@@ -150,7 +166,7 @@ namespace gpuNUFFT
       * Basic steps: - apodization correction
       *              - zero padding with osf
       *              - FFT
-      *							 - convolution and resampling
+      *               - convolution and resampling
       *
       * The memory for the output array is allocated automatically but has to be freed 
       * manually. 
@@ -171,6 +187,17 @@ namespace gpuNUFFT
       */
     virtual void     performForwardGpuNUFFT(Array<DType2> imgData,Array<CufftType>& kspaceData, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
 
+    /** \brief Perform forward gridding operation on kspaceData already residing in GPU memory
+     *
+     * This may be the case in iterative reconstructions, when k-Space and image
+     * data is already residing on the GPU.
+     * 
+     * @param image data 
+     * @param Stop gridding operation after gpuNUFFT::GpuNUFFTOutput 
+     * @return k-space data
+     */
+    virtual void     performForwardGpuNUFFT(GpuArray<DType2> imgData_gpu,GpuArray<CufftType>& kspaceData_gpu, GpuNUFFTOutput gpuNUFFTOut = DEAPODIZATION);
+
     /** \brief Perform forward gridding operation on given kspaceData 
       * 
       * The memory for the output array is allocated automatically but has to be freed 
@@ -301,12 +328,12 @@ namespace gpuNUFFT
     /** \brief Virtual forward convolution call, which can be used by sub-classes to add behaviour to the gridding steps
      *
      */
-    virtual void forwardConvolution(CufftType*		data_d, 
-      DType*			crds_d, 
-      CufftType*		gdata_d,
-      DType*			kernel_d, 
-      IndType*		sectors_d, 
-      IndType*		sector_centers_d,
+    virtual void forwardConvolution(CufftType*    data_d, 
+      DType*      crds_d, 
+      CufftType*    gdata_d,
+      DType*      kernel_d, 
+      IndType*    sectors_d, 
+      IndType*    sector_centers_d,
       gpuNUFFT::GpuNUFFTInfo* gi_host);
 
     /** \brief Virtual method to allow different methods of generation of the lookup table
@@ -367,7 +394,7 @@ namespace gpuNUFFT
     void initDeviceMemory(int n_coils);
 
     /** \brief Function to free the neccessary device memory used by the GriddingOperator. */
-    void freeDeviceMemory(int n_coils);
+    void freeDeviceMemory();
 
   };
 

diff --git a/CUDA/inc/gpuNUFFT_operator_factory.hpp b/CUDA/inc/gpuNUFFT_operator_factory.hpp
@@ -265,6 +265,13 @@ namespace gpuNUFFT
       */
     GpuNUFFTOperator* createNewGpuNUFFTOperator(IndType kernelWidth, IndType sectorWidth, DType osf, Dimensions imgDims);
 
+    /**
+     * \brief Function to check if the problem will fit into device memory
+     *
+     * @throws Exception in case of too much required memory 
+     */
+    void checkMemoryConsumption(Dimensions& kSpaceDims, const IndType& sectorWidth, const DType& osf, Dimensions& imgDims, Dimensions& densDims, Dimensions& sensDims);
+
   private:
     /** \brief Flag to indicate texture interpolation */
     bool useTextures;