[Rtk-users] unable to use CudaParkerShortScanImageFilter in python

Wed Jul 12 14:06:25 CEST 2023

Hi Simon,
Thank you for pointing that out. We made some corrections in the gantry
angle.
There was a bug in assuming the axis setting of ProjectionsReader origin
too. All those are corrected and reconstruction works perfectly fine now.

Currently we are reading projections as .tif images but ultimately we need
to read .raw images and apply reconstruction on them.
We have 512 unsigned short projections in .raw format each of with
1024x1024 dimensions.
It seems unsigned short is not supported by ProjectionsReader.

Two possible methods that we can think of:
1. Read projections with ImageSeriesReader which supports unsigned short,
converts data to float and feeds to ProjectionsReader.
   If so, is there a method to set manipulated input to ProjectionsReader?
>From the literature, only SetFileNames is possible.
2. Is there any alternative method to convert projections to line integral
without using ProjectionsReader.

Looking forward to your suggestions.

Regards,
Akshara

On Tue, 11 Jul 2023 at 02:09, Simon Rit <simon.rit at creatis.insa-lyon.fr>
wrote:

> Hi,
> The folder does not contain the angle.bin file. I looked at the
> projections. Maybe it's because the teeth are far from the central slice,
> where data are missing? This is known to cause problems, in particular for
> short scans.
> Simon
>
> On Mon, Jul 10, 2023 at 9:14 PM Akshara P K <akshara at advitech.in> wrote:
>
>> Hi Simon,
>>
>> Cuda PSSF is working fine with the wheel file you shared. But it looks
>> like something I missed from the reconstruction process and the result is
>> not up to mark as I expected.
>> The result looks like as below:
>> [image: image.png]
>> [image: image.png]
>> and the expected is more or less something like below:
>>
>> [image: image.png]
>>
>>
>> The input projections and the parameter file (Which I used to fill the
>> input parameters) are available here:
>>
>> https://drive.google.com/drive/folders/13EC2Lusyrx3HvkKJL80aSBHfled0FduT?usp=sharing
>>
>> Can you please check and provide your valuable suggestions on how I can
>> improve the result.
>>
>> Here is the updated code:
>>
>> # Defines the image type
>> ImageType = itk.Image[itk.F,3]
>> GPUImageType = itk.CudaImage[itk.F,3]
>>
>> # Defines the RTK geometry object
>> geometry = rtk.ThreeDCircularProjectionGeometry.New()
>> numberOfProjections = 512
>> firstAngle = -0.54853
>> angularArc = 270.66079712
>> sid = 537.09972 # source to isocenter distance
>> sdd = 908.00674 # source to detector distance
>> size_x = 1024
>> size_y = 1024
>> projOffsetX = -2.06667995
>> projOffsetY = -0.19235000
>> outOfPlaneAngle = 0
>> inPlaneAngle = 0
>> sourceOffsetX = 0
>> sourceOffsetY = 0
>>
>> angle_path = os.path.join("/skull/","angle.bin")
>> angles2 = []
>> with open(angle_path, 'rb') as file:
>>     while True:
>>         data = file.read(4)  # Read 4 bytes (32 bits) at a time
>>         if not data:
>>             break  # Reached end of file
>>         value = struct.unpack('f', data)[0]  # Unpack binary data as a
>> float
>>         angles2.append(value)
>> print("Read angles from bin completed")
>>
>> for x in angles2:
>>     geometry.AddProjection(sid,sdd,x,projOffsetX, projOffsetY,
>> outOfPlaneAngle, inPlaneAngle, sourceOffsetX, sourceOffsetY)
>>
>> tiffio = itk.TIFFImageIO.New()
>>
>> extension = ".tif"
>> fileNameFormat = "/binary/raw" + "%04d" + extension
>> fnames = itk.NumericSeriesFileNames.New()
>> fnames.SetStartIndex(0)
>> fnames.SetEndIndex(numberOfProjections-1)
>> fnames.SetIncrementIndex(1)
>> fnames.SetSeriesFormat(fileNameFormat)
>>
>> ProjectionsReaderType = rtk.ProjectionsReader[ImageType]
>> projectionsSource = ProjectionsReaderType.New()
>> projectionsSource.SetImageIO(tiffio)
>> projectionsSource.SetFileNames(fnames.GetFileNames())
>> projOrigin = [ -0.42*(size_x-1)/2, 0, -0.42*(size_y-1)/2 ] #input images
>> are size_x x size_y pixels with a 0.42mm pixel size
>> projSpacing = [ 0.42, 1, 0.42 ]
>> projectionsSource.SetOrigin( projOrigin )
>> projectionsSource.SetSpacing( projSpacing )
>>
>> # Graft the projections to an itk::CudaImage
>> projections = GPUImageType.New()
>> projectionsSource.Update()
>>
>> projections.SetPixelContainer(projectionsSource.GetOutput().GetPixelContainer())
>> projections.CopyInformation(projectionsSource.GetOutput())
>>
>> projections.SetBufferedRegion(projectionsSource.GetOutput().GetBufferedRegion())
>>
>> projections.SetRequestedRegion(projectionsSource.GetOutput().GetRequestedRegion())
>>
>> ParkerShortScanImageFilterType = rtk.CudaParkerShortScanImageFilter
>> ParkerFilter = ParkerShortScanImageFilterType.New()
>> ParkerFilter.SetInput(projections)
>> ParkerFilter.SetGeometry(geometry)
>> ParkerFilter.InPlaceOff()
>>
>> ConstantImageSourceType = rtk.ConstantImageSource[GPUImageType]
>> # Create reconstructed image
>> constantImageSource2 = ConstantImageSourceType.New()
>> render_x = 512
>> render_y = 512
>> sizeOutput = [ render_x, render_y, numberOfProjections ]
>> origin = [ -(render_x/2.0), 0 , -(render_y/2.0) ]
>> spacing = [ 1.0, 0.5, 1.0 ]
>> constantImageSource2.SetOrigin( origin )
>> constantImageSource2.SetSpacing( spacing )
>> constantImageSource2.SetSize( sizeOutput )
>> constantImageSource2.SetConstant(0.)
>>
>> #Displaced detector weighting
>> ddf = rtk.CudaDisplacedDetectorImageFilter.New()
>> ddf.SetInput(ParkerFilter.GetOutput())
>> ddf.SetGeometry(geometry)
>>
>> # FDK reconstruction
>> print("Reconstructing...")
>> FDKGPUType = rtk.CudaFDKConeBeamReconstructionFilter
>> feldkamp = FDKGPUType.New()
>> feldkamp.SetInput(0, constantImageSource2.GetOutput())
>> feldkamp.SetInput(1, ddf.GetOutput())
>> feldkamp.SetGeometry(geometry)
>> feldkamp.GetRampFilter().SetTruncationCorrection(0.0)
>> feldkamp.GetRampFilter().SetHannCutFrequency(0.0)
>>
>> # Field-of-view masking
>> FOVFilterType = rtk.FieldOfViewImageFilter[ImageType, ImageType]
>> fieldofview = FOVFilterType.New()
>> fieldofview.SetInput(0, feldkamp.GetOutput())
>> fieldofview.SetProjectionsStack(feldkamp.GetOutput())
>> fieldofview.SetGeometry(geometry)
>>
>> # Writer
>> print("Writing output image...")
>> WriterType = rtk.ImageFileWriter[ImageType]
>> writer = WriterType.New()
>> writer.SetFileName("image_output_skull00_ver1.3.tiff")
>> writer.SetInput(fieldofview.GetOutput())
>> writer.Update()
>>
>> Regards,
>> Akshara
>>
>> On Thu, 6 Jul 2023 at 18:23, Akshara P K <akshara at advitech.in> wrote:
>>
>>> Hi Simon,
>>>
>>> Thank you for the quick response.
>>> It works perfectly fine except some warnings as you mentioned.
>>> Thank you again for the support.
>>>
>>> Regards,
>>> Akshara
>>>
>>> On Wed, 5 Jul 2023 at 18:36, Simon Rit <simon.rit at creatis.insa-lyon.fr>
>>> wrote:
>>>
>>>> Hi,
>>>> Following a recent similar request (see here
>>>> <https://www.creatis.insa-lyon.fr/pipermail/rtk-users/2023-June/001979.html>),
>>>> I have added CudaParkerShortScanImageFilter with this PR
>>>> <https://www.creatis.insa-lyon.fr/pipermail/rtk-users/2023-June/001979.html>.
>>>> New Python packages are available as artifacts of the CI here
>>>> <https://github.com/RTKConsortium/RTK/actions/runs/5448263998>. I
>>>> haven't tested them, let me know whether they work... In any case, I have
>>>> noticed new warnings when using these wheels which must be fixed (see issue
>>>> #553 <https://github.com/RTKConsortium/RTK/issues/553>).
>>>> Simon
>>>>
>>>> On Wed, Jul 5, 2023 at 2:20 PM Akshara P K <akshara at advitech.in> wrote:
>>>>
>>>>> Hi all,
>>>>>
>>>>> I was trying to use FirstCudaReconstruction.py to reconstruct a set of
>>>>> tiff projections images using the python binary of itk-rtk. Since the angle
>>>>> coverage is around 268 degrees, I had to implement a
>>>>> ParkerShortScanImageFilter between the output of the ProjectionsReader and
>>>>> the FDK reconstruction and it is working fine with the CPU version of API.
>>>>> However, the use of CudaParkerShortScanImageFilter causes an error saying
>>>>> module 'RTK' has no attribute 'CudaParkerShortScanImageFilter' even though
>>>>> CudaParkerShortScanImageFilter is added and implemented in the cpp version.
>>>>> Can anyone provide an insight on how to solve this issue?
>>>>>
>>>>> Regards,
>>>>> Akshara
>>>>> _______________________________________________
>>>>> 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/20230712/86bd7d24/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 185685 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20230712/86bd7d24/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 10665 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20230712/86bd7d24/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 187374 bytes
Desc: not available
URL: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20230712/86bd7d24/attachment-0005.png>