Python support

Added interface to support Python.

Bindings written by Chaithya G R and Carole Lazarus.

Improved deapodization

• Apply deapodization by computing inverse FFT of gridding kernel (allows for correct scaling independent of kernel width)
• Apply sqrt of density compensation in both forward and adjoint operation (better usability in iterative algorithms)

2D performance upgrade

Processing of multiple coils at once depending on available resources

GPU Arrays

GPU Arrays

The gpuNUFFT operator now supports data vectors which are already located on the GPU device, which is exemplarily necessary when performing iterative reconstructions in standalone CPP applications. The Matlab interface is not affected by this extension.

See corresponding test case:

TEST(TestForwardBackward,Test_GpuArray)
{
  //Test the same as above but use GpuArray data structure  
  int kernel_width = 3;
  float osf = 1.25;//oversampling ratio
  int sector_width = 8;
  //Data
  int data_entries = 2;
  DType2* data = (DType2*) calloc(data_entries,sizeof(DType2)); //2* re + im
  data[0].x = 5;//Re
  data[0].y = 0;//Im
  data[1].x = 1;//Re
  data[1].y = 0;//Im

  //Coords
  //Scaled between -0.5 and 0.5
  //in triplets (x,y,z) as structure of array
  //p0 = (0,0,0)
  //p1 0 (0.25,0.25,0.25)
  DType* coords = (DType*) calloc(3*data_entries,sizeof(DType));//3* x,y,z
  coords[0] = 0.00; //x0
  coords[1] = 0.25; //x1

  coords[2] = 0.00; //y0
  coords[3] = 0.25; //y1

  coords[4] = 0.00; //z0
  coords[5] = 0.25; //z1

  //Input data array, complex values
  //and copy to GPU 
  gpuNUFFT::GpuArray<DType2> dataArray_gpu;
  dataArray_gpu.dim.length = data_entries;
  allocateAndCopyToDeviceMem<DType2>(&dataArray_gpu.data,data,data_entries);

  //Input array containing trajectory in k-space
  gpuNUFFT::Array<DType> kSpaceData;
  kSpaceData.data = coords;
  kSpaceData.dim.length = data_entries;

  gpuNUFFT::Dimensions imgDims;
  imgDims.width = 64;
  imgDims.height = 64;
  imgDims.depth = 64;

  //precomputation performed by factory
  gpuNUFFT::GpuNUFFTOperatorFactory factory; 
  gpuNUFFT::GpuNUFFTOperator *gpuNUFFTOp =    
  factory.createGpuNUFFTOperator(kSpaceData,kernel_width,sector_width,osf,imgDims);

  //Output Array
  gpuNUFFT::GpuArray<CufftType> imgArray_gpu;
  imgArray_gpu.dim = imgDims;
  allocateDeviceMem<CufftType>(&imgArray_gpu.data,imgArray_gpu.count());

  //Perform FT^H Operation
  gpuNUFFTOp->performGpuNUFFTAdj(dataArray_gpu, imgArray_gpu);

  //Perform FT Operation
  gpuNUFFTOp->performForwardGpuNUFFT(imgArray_gpu,dataArray_gpu);
  copyFromDevice(dataArray_gpu.data,data,data_entries);

  printf("contrast %f \n",data[0].x/data[1].x);
  EXPECT_NEAR(data[0].x/data[1].x,5.0,epsilon);

  free(data);
  free(coords);
  freeTotalDeviceMemory(dataArray_gpu.data,imgArray_gpu.data,NULL);

  delete gpuNUFFTOp;
}

Memory estimator

Additionally to the initialization step of the gpuNUFFT operator a rough memory estimation is performed, in order to avoid crashes when the problem data exceeds the total device memory. To handle this, an exception is thrown and a Matlab error is returned.

Automatic coil summation

Arbitrary sector widths

The sector width parameter can now be arbitrarily chosen in such manner that the total grid dimension (image dim * oversampling factor) does not have to be an integer multiple of the sector width, e.g.:

imgDim = [256,256,240]; 
sw = 8; 
osf = 1.125; 
% imgDim * osf ./ sw = [36, 36, 33.75] 
...
FT = gpuNUFFT(k',w,osf,kw,sw,imgDim); % still works
...

Automatic coil summation

If the gpuNUFFT operator is initialized with coil sensitivities following steps are performed automatically on the GPU:

%without sens data
FT = gpuNUFFT(k',w,osf,wg,sw,[N,N,nSl]);

%with sens data
FT = gpuNUFFT(k',w,osf,wg,sw,[N,N,nSl],sens);

Adjoint operation: sum over coils as final step

% the block
for ii = 1:nc
    y = y + conj(sens(:,:,:,ii) .* FH(data(:,ii)); 
end

% can be replaced by 
y = FH(data);

Forward operation: loop over coils

% the block
img = repmat(img,[1,1,1,nCh]);
for ii = 1:nc
    data(:,ii) = F( sens(:,:,:,ii) .* img(:,:,:,ii) ); 
end

% can be replaced by 
y = F(img);

Varying dimensions in x, y and z

Varying dimensions in all directions, e.g. 3d grids with 256x128x64 or 2d grids with 256x128 are supported now.

Texture lookup, balanced workload, sensitivity data

• Use textures for kernel lookup
• Use textures for data readout
• Balanced and textured gridding to work more efficiently
• Apply sens data on GPU
• Automatically export mex-files into gpuNUFFT/@gpuNUFFT/private folder
• Refactored initialization params of MATLAB Operator

Anisotropic grid support

Now anisotropic grid sizes, e.g. 128x128x64 are supported. Generated test data and MATLAB scripts are added in the demo folder.

Parallel precomputation

This release contains the precomputation step on the GPU. This will lead to better performance if it is neccessary to continuously create the Gridding Operator, e.g. when the k-Space trajectory is shifted per iteration step.