[Rtk-users] unable to use CudaParkerShortScanImageFilter in python

[Simon Rit]

Check this archive on the mailing list:
https://www.creatis.insa-lyon.fr/pipermail/rtk-users/2023-January/001921.html
ProjectionsReader supports any input raw but converts it to the required
output, i.e., float or double.

On Wed, Jul 12, 2023 at 2:06 PM Akshara P K <akshara at advitech.in> wrote:

> 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/c9a7f4f0/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/c9a7f4f0/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/c9a7f4f0/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/c9a7f4f0/attachment-0005.png>