[Rtk-users] [EXTERNAL] Help with CudaForwardProjectionImageFilter
Simon Rit
simon.rit at creatis.insa-lyon.fr
Fri Apr 26 09:03:29 CEST 2024
Hi Obaid,
If you remove the "ForwardProj.Update()" line and connect instead a
StreamingImageFilter
<https://itk.org/Doxygen/html/classitk_1_1StreamingImageFilter.html>, your
projections can be processed piece by piece.
I'm surprised by your 80 GB memory, can you display the output of
nvidia-smi?
Simon
On Thu, Apr 25, 2024 at 11:56 AM Rahman, Obaid <rahmano at ornl.gov> wrote:
> 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>
>
>
> _______________________________________________
> 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/20240426/6675953a/attachment-0001.htm>
More information about the Rtk-users
mailing list