[Rtk-users] unable to use CudaParkerShortScanImageFilter in python

Thu Jul 13 07:26:40 CEST 2023

Hi Simon,

Thank you for the quick response.

We already tried that conversion. Here is the updated snippet:
size_x = 1024
size_y = 1024
numberOfProjections = 512
rawio = itk.RawImageIO.New()
rawio.SetNumberOfDimensions(2)
rawio.SetPixelType(itk.CommonEnums.IOPixel_SCALAR)
rawio.SetComponentType(itk.CommonEnums.IOComponent_FLOAT)
rawio.SetDimensions(0, size_x)
rawio.SetDimensions(1, size_y)
rawio.SetSpacing(0, 1)
rawio.SetSpacing(1, 1)
rawio.SetByteOrderToLittleEndian()

# List of filenames
fileNames = []
for i in range(numberOfProjections-1):
    fileNames.append(f'/subvolume00/{i:04d}.raw')

pixelSpacing = 0.420002
projOrigin = [ -pixelSpacing*(size_x-1)/2, 0, -pixelSpacing*(size_y-1)/2 ]
#input images are size_x x size_y pixels with a 0.42mm pixel size
projSpacing = [ pixelSpacing, 1, pixelSpacing ]
# Read projections
CPUImageType = itk.Image[itk.F,3]
proj = rtk.projections_reader(ttype=CPUImageType,
file_names=fileNames,ImageIO=rawio, Origin=projOrigin, Spacing=projSpacing)

Then it is throwing below error:
*ITK ERROR: RawImageIO(0000027078EE8D00): File seek failed*

Then tried setting header as below
rawio.SetHeaderSize(rawio.GetImageSizeInPixels() * 0)

and got another error saying
*ITK ERROR: RawImageIO(000001A2C02FAB20): Read failed: Wanted 4194304
bytes, but read 2097152 bytes.*

To read appropriate number of bytes from each file, we need to use an
unsigned int (16 bit) type:
rawio.SetComponentType(itk.CommonEnums.IOComponent_INT)
# Read projections
CPUImageType = itk.Image[itk.US,3]
proj = rtk.projections_reader(ttype=CPUImageType,
file_names=fileNames,ImageIO=rawio, Origin=projOrigin, Spacing=projSpacing)

the above changes throws error saying template not found.

*No suitable template parameter can be found.Please specify an input via
the first argument, the 'ttype' keyword parameter,or via one of the
following keyword arguments: Input, InputImage, Input1*

Regards,
Akshara

On Wed, 12 Jul 2023 at 22:22, Simon Rit <simon.rit at creatis.insa-lyon.fr>
wrote:

> 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/20230713/593a4a8a/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/20230713/593a4a8a/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/20230713/593a4a8a/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/20230713/593a4a8a/attachment-0005.png>