[Rtk-users] Help with CudaForwardProjectionImageFilter

[Rahman, Obaid]

Hello all,

I am getting “out of memory” error when I run CudaForwardProjectionImageFilter.
Please refer to the attached screenshot.

These are the array sizes with which I get the error:
Projection size: (slice/row, view, column)=(1456,145,1840)
Recon size: (z, y, x)=(1264, 1356, 1356)

When I increase the voxel size to twice i.e. Recon size: (z, y, x)=(632, 678, 678), it works fine.
Maybe I’m not being very efficient with the memory?

If anyone has suggestion to make it more memory efficient, please let me know.
Thanks.

Best,
Obaid

Here’s what I’m doing:

def rtk_projection(recon_arrar, proj_params, miscalib, vol_params, return_itk=False, return_RVC=True):
    '''
    inputs:
      recon_arrar: image numpy array (ZXY format is expected)
      proj_params: projection parameter dictionary
      miscalib: miscalibation parameter dictionary
      vol_params: image volume parameter dictionary
      return_itk: True if itk image is expected, False if a numpy array is expected
    output:
      itk image (c,r,v format internally) or image numpy array (r,v,c format) of projection data
    '''
    ImageType = itk.Image[itk.F,3]; CudaImageType = itk.CudaImage[itk.F,3]

    # Convert the recon to itk format
    ###########################################################################
    # proj_data is in (z,x,y); rtk requires  numpy array in (x,z,y)
    recon_arrar = np.transpose(recon_arrar, [1,0,2]).astype('float32') # now in (x,z,y)
    recon_shape = recon_arrar.shape
    recon_arrar = itk.GetImageFromArray(recon_arrar) # internally img has format (y,z,x)

    # Center the image around 0 which is the default center of rotation
    recon_arrar.SetOrigin([-0.5*(recon_arrar.GetLargestPossibleRegion().GetSize()[0]-1)*recon_arrar.GetSpacing()[0],
                    -0.5*(recon_arrar.GetLargestPossibleRegion().GetSize()[1]-1)*recon_arrar.GetSpacing()[1],
                    -0.5*(recon_arrar.GetLargestPossibleRegion().GetSize()[2]-1)*recon_arrar.GetSpacing()[2]])

    # Graft the projections to an itk::CudaImage
    ###########################################################################
    vol_xy = float(vol_params['vox_xy']); vol_z = float(vol_params['vox_z'])
    rtk_recon = CudaImageType.New()# img will have format crv
    rtk_recon.SetPixelContainer(recon_arrar.GetPixelContainer()) # img has format crv
    rtk_recon.SetLargestPossibleRegion(recon_arrar.GetLargestPossibleRegion())
    rtk_recon.SetBufferedRegion(recon_arrar.GetBufferedRegion())
    rtk_recon.SetRequestedRegion(recon_arrar.GetRequestedRegion())
    rtk_recon.SetSpacing([vol_xy, vol_z, vol_xy]) # spacing for xzy
    rtk_recon.SetOrigin([-0.5*(recon_shape[2]-1)*vol_xy, # origin for x direction
                  -0.5*(recon_shape[1]-1)*vol_z, # origin for z direction
                  -0.5*(recon_shape[0]-1)*vol_xy])  # origin for y direction
    # print('GPU recon:', rtk_recon)
    del recon_arrar
    ###########################################################################

    # Defines the RTK geometry object
    geometry = rtk.ThreeDCircularProjectionGeometry.New()
    for i in range(len(proj_params['angles'])):
        angle_deg = proj_params['angles'][i]*180/np.pi # angle in degree
        geometry.AddProjection(proj_params['cone_params']['src_orig'],
                               proj_params['cone_params']['src_orig']+proj_params['cone_params']['orig_det'],
                               angle_deg, miscalib['delta_u'], miscalib['delta_v'])

    # define some parameters (to make code more readable)
    det_pix_x = proj_params['cone_params']['pix_x'] # row direction pixel size
    det_pix_y = proj_params['cone_params']['pix_y'] # channel direction pixel size
    proj_shape = [int(proj_params['dims'][0]), int(proj_params['dims'][1]), int(proj_params['dims'][2])] # r,v,c
    ###########################################################################

    # define a zero projection (GPU)
    constantImageSource = rtk.ConstantImageSource[CudaImageType].New()
    constantImageSource.SetSpacing( [det_pix_y, det_pix_x, 1.] ) # c,r,v
    constantImageSource.SetSize( [proj_shape[2], proj_shape[0], proj_shape[1]] ) # c,r,v
    constantImageSource.SetOrigin([ -0.5*(proj_shape[2]-1)*det_pix_y,   # origin for channel direction
                                    -0.5*(proj_shape[0]-1)*det_pix_x,   # origin for row direction
                                    -0.5*(proj_shape[1]-1)*1])          # origin for view direction (will be ignored anyway)
    constantImageSource.SetConstant(0.)
    ###########################################################################

    # forward project the recon to fill out the zero projection image
    ForwardProj = rtk.CudaForwardProjectionImageFilter[CudaImageType].New()
    ForwardProj.SetGeometry( geometry )
    ForwardProj.SetInput(0, constantImageSource.GetOutput()) # projection image
    ForwardProj.SetInput(1, cpu2gpu(rtk_recon)) # recon volume
    ForwardProj.SetStepSize(float(vol_params['vox_xy'])/4) # step size along the ray (default is 1mm)
    ForwardProj.Update()
    ###########################################################################

    # graft the projection to CPU / extract the array
    rtk_reprojection = gpu2cpu(ForwardProj.GetOutput()) # array inside the image has c,r,v format
    if return_itk:
        return rtk_reprojection # array inside the image has c,r,v format
    else:
        rtk_reprojection = itk.GetArrayViewFromImage(rtk_reprojection) # now v,r,c format
        if return_RVC:
            rtk_reprojection = np.transpose(rtk_reprojection, [1,0,2]) # r,v,c format to match Astra projection
        return rtk_reprojection # numpy array (r,v,c)





-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20240425/a18acb74/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screenshot 2024-04-24 at 8.34.20?PM.png
Type: image/png
Size: 175206 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20240425/a18acb74/attachment-0001.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 1611 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20240425/a18acb74/attachment-0001.p7s>