[Rtk-users] unable to use CudaParkerShortScanImageFilter in python

Thu Jul 13 09:53:16 CEST 2023

Hi,
As I said, the output type (set by ttype) can only be float or double even
if the projections are unsigned short.
But you need to change the component type for the input to
IOComponent_USHORT (2 bytes), IOComponent_INT is 4 bytes.
Simon

On Thu, Jul 13, 2023 at 7:27 AM Akshara P K <akshara at advitech.in> wrote:

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