[Rtk-users] Issue spr/Boellaard? + First recon from raw C++ vs python

[Simon Rit]

Hi,
1/ This implementation is known to have some issues and I think it should
be revisited. I have opened an issue to track this
https://github.com/SimonRit/RTK/issues/454
If I remember correctly, someone identified at some point a problem with
the patient mask, which was selecting air instead. In any case, the SPR
might have no influence because one hits the non negativity constraints
<http://www.openrtk.org/Doxygen/classrtk_1_1BoellaardScatterCorrectionImageFilter.html#ae8752f068ecb8b3208b360c1375ae0a6>.
If you're willing to contribute, the code
<https://github.com/SimonRit/RTK/blob/master/include/rtkBoellaardScatterCorrectionImageFilter.hxx#L64-L107>
is fairly simple so feel free to investigate.
2/ I think that the problem is that FDK expect a full scan and you have a
short scan. So the projections next to the gap get a huge weight. You need
to apply Parker weighting similarly to rtkfdk, see here
<https://github.com/SimonRit/RTK/blob/master/applications/rtkfdk/rtkfdk.cxx#L102-L117>
.
>From what I see in the code, the volume is not centered on the isocenter. I
would also correct it with

originVol = [-(sizeOutputVol[0] / 2 - 0.5) * spacingVol[0],
-(sizeoutputVol[1] / 2 - 0.5) * spacingVol[1], -(sizeOutputVol[2] / 2
- 0.5) * spacingVol[2]]

I hope it helps. Good luck with all this!
Simon

On Fri, Nov 5, 2021 at 4:07 PM Maspero, M. (Matteo) <M.Maspero at umcutrecht.nl>
wrote:

> Dear rtk-users,
>
>
> First of all a huge thanks to the developers of RTK, what a very powerful
> tool. Here I am am with a couple of points:
>
>
> 1) Boellaard scatter: after having performed multiple rtkfkd recons with
> different "spr" settings (from the complied C++ version), I could not find
> any differences between the reconstructed images. It seems like the
> scatter-to-primary ratio does not impact the recon. Is this to be
> expected/does anyone experience the same?
>
>
> 2) Finally I had some time to dig in the python wrapped side of RTK and
> run some very basic recons on Elekta raw (XVI release > 5) data - based on
> the FirstRecon example - for benchmarking purposes. However, I do not get
> to the desired result. I am pretty sure that is something quite superficial
> on my side, but any help is very welcome.
>
> So the reconstruction starts from the following projections:
>
>
> I get with the compiled version of C++
>
>
> while when trying my own recon, I get
>
>
>
> To me, it seems that the geometry/axis orientation is somewhat incorrect.
> Trying to facilitate coming up with some suggestions, I put at the end of
> the mail an excerpt of the code, maybe you spot the issue right away. If
> you may appreciate a dataset to reproduce the problem, just let me know.
> Again, any suggestion is welcome, especially if a quick look from someone
> else may speed up things on this side.
>
>
> Thank you very much in advance.
>
>
> Best Regards and have a great weekend,
>
>
> Matteo
>
>
>   Matteo Maspero
> | Clinical scientist/Postdoc
> | Radiotherapy, Imaging and Cancer
> | Computational Imaging Group www.compimag.org
> | Center for Image Sciences
> | University Medical Center Utrecht
> | Free on Wednesdays
> | Room Q.02.2.315
> | P.O. 85500
> | Heidelberglaan 100, 3508 GA Utrecht, Netherlands
> | E: m.maspero at umcutrecht.nl
> | Tel        +31-8875 67492 / +31-614956534
> | Fax       +31-8875 55850
>
> --------------------- code excerpt ----
>
> import sys
> import os
> import matplotlib.pyplot as plt
> import glob
> import itk
> from itk import RTK as rtk
> # image definition
> CPUImageType = rtk.Image[itk.F, 3]
>
> # TMP, pt_dir are the folders with the data
> # Flags
> Flag_saveProj = "on"
>
> # geometry -> read xml
> geometryReader = rtk.ThreeDCircularProjectionGeometryXMLFileReader.New()
> geometryReader.SetFilename(TMP + "elektaGeometry")
> geometryReader.GenerateOutputInformation()
> geometry = geometryReader.GetGeometry()
>
> # List of filenames
> fileNames = list()
> # Get list of all files in a given directory sorted by name
> fileNames = sorted( filter( os.path.isfile,
>                             glob.glob(pt_dir + '*his')))
>
> projReader = rtk.ProjectionsReader[CPUImageType].New()
> projReader.SetFileNames(fileNames)
> projections = projReader.GetOutput()
> # plot the projections
> plt.subplot(1, 3, 1); plt.imshow(projections[55, :, :], cmap='gray'); plt.title('Proj X')
> plt.subplot(1, 3, 2); plt.imshow(projections[:, 55, :], cmap='gray'); plt.title('Proj Y')
> plt.subplot(1, 3, 3); plt.imshow(projections[:, :, 50], cmap='gray'); plt.title('Proj Z')
> plt.show()
>
> # Writer
> if Flag_saveProj == "on":
>     print('Writing projections image...')
>     itk.imwrite(projections, TMP + 'projections.mha')
>     itk.imwrite(projections, TMP + 'projections.gipl')
>
> # Reconstruction parameter
> spacingVol = [1, 1, 1]           # resolution of the recon
> sizeOutputVol = [256, 256, 256]  # FOV = S ->  [270, 270, 264]
> originVol = [-(sizeOutputVol[0] / 2 + 0.5) * spacingVol[0], 0., -(sizeOutputVol[2] / 2 + 0.5) * spacingVol[2]]
>
> # Create (an empty) reconstructed image
> constantImageSource = rtk.ConstantImageSource[CPUImageType].New()
> constantImageSource.SetOrigin(originVol)
> constantImageSource.SetSpacing(spacingVol)
> constantImageSource.SetSize(sizeOutputVol)
> constantImageSource.SetConstant(0.)
> source = constantImageSource.GetOutput()
>
> # FDK reconstruction
> print('Reconstructing...')
> FDKCPUType = rtk.FDKConeBeamReconstructionFilter[CPUImageType]
> feldkamp = FDKCPUType.New()
> feldkamp.SetInput(0, source)
> feldkamp.SetInput(1, projections)
> feldkamp.SetGeometry(geometry)
> feldkamp.GetRampFilter().SetTruncationCorrection(0.0)
> feldkamp.GetRampFilter().SetHannCutFrequency(0.0)
> imageFDK = feldkamp.GetOutput()
>
> # plot the reconstructed image
> im = imageFDK
> plt.subplot(1, 3, 1); plt.imshow(im[55, :, :], cmap='gray'); plt.title('X')
> plt.subplot(1, 3, 2); plt.imshow(im[:, 55, :], cmap='gray'); plt.title('Y')
> plt.subplot(1, 3, 3); plt.imshow(im[:, :, 50], cmap='gray'); plt.title('Z'); plt.show()
>
>
> # Reconstruction parameter
> spacingVol = [1, 1, 1]           # resolution of the recon
> sizeOutputVol = [256, 256, 256]  # FOV = S ->  [270, 270, 264]
> originVol = [-(sizeOutputVol[0] / 2 + 0.5) * spacingVol[0], 0., -(sizeOutputVol[2] / 2 + 0.5) * spacingVol[2]]
>
> # Create (an empty) reconstructed image
> constantImageSource = rtk.ConstantImageSource[CPUImageType].New()
> constantImageSource.SetOrigin(originVol)
> constantImageSource.SetSpacing(spacingVol)
> constantImageSource.SetSize(sizeOutputVol)
> constantImageSource.SetConstant(0.)
> source = constantImageSource.GetOutput()
>
> # FDK reconstruction
> print('Reconstructing...')
> FDKCPUType = rtk.FDKConeBeamReconstructionFilter[CPUImageType]
> feldkamp = FDKCPUType.New()
> feldkamp.SetInput(0, source)
> feldkamp.SetInput(1, projections)
> feldkamp.SetGeometry(geometry)
> feldkamp.GetRampFilter().SetTruncationCorrection(0.0)
> feldkamp.GetRampFilter().SetHannCutFrequency(0.0)
> imageFDK = feldkamp.GetOutput()
>
> # plot the reconstructed image
> im = imageFDK
> plt.subplot(1, 3, 1); plt.imshow(im[55, :, :], cmap='gray'); plt.title('X')
> plt.subplot(1, 3, 2); plt.imshow(im[:, 55, :], cmap='gray'); plt.title('Y')
> plt.subplot(1, 3, 3); plt.imshow(im[:, :, 50], cmap='gray'); plt.title('Z'); plt.show()
>
>
> ------------------------------
>
> * De informatie opgenomen in dit bericht kan vertrouwelijk zijn en is
> uitsluitend bestemd voor de geadresseerde. Indien u dit bericht onterecht
> ontvangt, wordt u verzocht de inhoud niet te gebruiken en de afzender
> direct te informeren door het bericht te retourneren. Het Universitair
> Medisch Centrum Utrecht is een publiekrechtelijke rechtspersoon in de zin
> van de W.H.W. (Wet Hoger Onderwijs en Wetenschappelijk Onderzoek) en staat
> geregistreerd bij de Kamer van Koophandel voor Midden-Nederland onder nr.
> 30244197. *
>
> * Denk s.v.p aan het milieu voor u deze e-mail afdrukt. *
> ------------------------------
>
> * This message may contain confidential information and is intended
> exclusively for the addressee. If you receive this message unintentionally,
> please do not use the contents but notify the sender immediately by return
> e-mail. University Medical Center Utrecht is a legal person by public law
> and is registered at the Chamber of Commerce for Midden-Nederland under no.
> 30244197. *
>
> * Please consider the environment before printing this e-mail. *
> _______________________________________________
> Rtk-users mailing list
> Rtk-users at public.kitware.com
> https://public.kitware.com/mailman/listinfo/rtk-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20211110/2b436c40/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Projection.png
Type: image/png
Size: 54921 bytes
Desc: not available
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20211110/2b436c40/attachment-0003.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ReconC++.png
Type: image/png
Size: 86799 bytes
Desc: not available
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20211110/2b436c40/attachment-0004.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ReconPy.png
Type: image/png
Size: 63784 bytes
Desc: not available
URL: <https://public.kitware.com/pipermail/rtk-users/attachments/20211110/2b436c40/attachment-0005.png>