[Rtk-users] [EXTERNAL] Help with CudaForwardProjectionImageFilter
Rahman, Obaid
rahmano at ornl.gov
Thu Apr 25 11:54:19 CEST 2024
Thanks, Nils for the swift response.
I have actually 80 GB memory per GPU which should be plenty.
The image is 8.65 GB (read as float32), and the projection size is 1.5 GB (read as float32).
Also, CudaBackProjectionImageFilter works fine for the given projection and image sizes.
I have tried other forward projectors (Astra) and they work fine with the given memory.
I was wondering if my code was suboptimal memory-wise.
Thanks!
Best,
Obaid
> On Apr 25, 2024, at 3:54 AM, krah <nils.krah at creatis.insa-lyon.fr> wrote:
>
> 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)
>>
>>
>> <Screenshot 2024-04-24 at 8.34.20 PM.png>
>>
>>
>> _______________________________________________
>> Rtk-users mailing list
>> rtk-users at openrtk.org
>> https://www.creatis.insa-lyon.fr/mailman/listinfo/rtk-users <https://urldefense.us/v2/url?u=https-3A__www.creatis.insa-2Dlyon.fr_mailman_listinfo_rtk-2Dusers&d=DwQFaQ&c=v4IIwRuZAmwupIjowmMWUmLasxPEgYsgNI-O7C4ViYc&r=J7uT21mkGp7aMwIrHQkTLGwy72wKx_bOB0IkoGp__bQ&m=mi_1gmtry3HzMfO0caEhCP-0wcNe3su2gocFndAoByaorId18bE57XL9ozhaj01k&s=K4sXSgvZVGEDHgEKoU168ighF9MBDAXWRLef6iJDkrs&e=><Screenshot 2024-04-24 at 8.34.20 PM.png>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20240425/88c8ee50/attachment-0001.htm>
-------------- 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/88c8ee50/attachment-0001.p7s>
More information about the Rtk-users
mailing list