[Rtk-users] Help with CudaForwardProjectionImageFilter
krah
nils.krah at creatis.insa-lyon.fr
Thu Apr 25 09:54:36 CEST 2024
Hey,
without knowing further details nor having verified the code, I guess your image is just too large to fit into the memory of your GPU.
Naïvely, I would calculate the size in GBytes of your image and compare with the memory of your GPU. Increasing the spacing means using fewer pixels = less memory. On top of that, I guess there are the projections which require space on the GPU.
If you need a small spacing, maybe you can reconstruct in chunks along the rotation axis and merge the chunks later?
Or the RTK experts know other tricks ...
just my non-expert two-cents ...
Nils
On Apr 25 2024, at 3:21 am, Rahman, Obaid <rahmano at ornl.gov> wrote:
>
>
> 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)
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> _______________________________________________
>
> Rtk-users mailing list
>
> rtk-users at openrtk.org
>
> https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20240425/d4625478/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/d4625478/attachment-0001.png>
More information about the Rtk-users
mailing list