[Rtk-users] Creating DRRs with RTK in python

[Simon Rit]

Hi,
Your center and semiprincipalaxis variables are not used, you can delete
them. I think they are leftovers from the example you used, which was
probably the projection of an ellipsoid.
The center of rotation is the location of point with coordinates (0,0,0)
mm. You must modify the origin of ctpreop.mha to change it. For example
with ChangeInformationImageFilter.
Your results are indeed weird. I would check ctpreop.mha because it looks
like a thin CT, probably something is wrong in the antero-posterior
direction. Or is it just a slice?
Yes, you can set any source and detector positions so you can change the
rotation axis. You can use one of the AddProjection functions of the
geometry object to do so. This is also a way to change the rotation center.
Simon

On Wed, Aug 18, 2021 at 12:06 PM Milou de Boer <
m.m.deboer at student.utwente.nl> wrote:

> Dear Simon,
>
>
>
> After modifying and trying out the script you sent me, I have some
> questions since I don’t get the desired result. My desired result is to get
> projections from a patient CT from all 360 degrees. However, enclosed are
> some pictures from my actual result (projection_result.png). It seems that
> the isocenter around which the projections are made is wrong. I think I
> should change the center variable, but whatever I change it to, the results
> are the same. Wat is going wrong?
>
>
>
> Changing the variables origin and semiprincipalaxis seem to do nothing as
> well to my result. Just for clarification, what do these variables mean? To
> my understanding, the origin variable sets the world coordinates for where
> the detector (and thus projection images) initially is. The center defines
> the isocenter around which the source and detector rotate.  I don’t really
> understand the semiprincipalaxis variable, it probably defines the shape of
> the path which is made by the source/detector. Are these explanations
> correct or am I misinterpreting something?
>
>
>
> Lastly, I have another question: Is it also possible to change around
> which axis the projections are made? And is it also possible to first do a
> rotation around the z axis and then the y axis leading to 360x360
> projection images?
>
>
>
> Thank you very much.
>
>
>
> Best regards,
>
> Milou
>
>
>
> *Van: *Simon Rit <simon.rit at creatis.insa-lyon.fr>
> *Verzonden: *vrijdag 13 augustus 2021 08:25
> *Aan: *Milou de Boer <m.m.deboer at student.utwente.nl>
> *CC: *rtk-users at public.kitware.com
> *Onderwerp: *Re: [Rtk-users] Creating DRRs with RTK in python
>
>
>
> Hi,
>
> Both your input volume and your stack of projections should be 3D. If you
> simulate only one projection, you can always make it 2D afterwards by
> slicing the stack. The third dimension of spacing and origin of the stack
> are never used in RTK but must be provided since it's a 3D image.
>
> I might be wrong but I don't think SimpleITK and the RTK wrappings are
> compatible. I would suggest the use of the ITK wrappings described in their
> examples. See here
> <https://itkpythonpackage.readthedocs.io/en/master/Quick_start_guide.html>
> for an example.
>
> Enclosed is a slightly modified version of your script which produces a
> projection stack. The geometry is probably not correct but I let you play
> with the script to set them as you wish.
>
> Cheers,
>
> Simon
>
>
>
> On Mon, Aug 9, 2021 at 9:55 AM Milou de Boer <
> m.m.deboer at student.utwente.nl> wrote:
>
>
>
> Dear all,
>
>
>
> I’m quite new to Python (and the use of RTK and ITK), but I have a project
> for which I want to create digitally reconstructed radiographs (DRRs) of a
> CT. I therefore have successfully installed RTK for python and I have run
> the FirstReconstruction.py example. Now I want to edit this so that it is
> possible to create DRRs, thus forwardprojection.
>
>
>
> I have read in the archives of RTK-users and people are mostly directed
> towards the C++ code
> https://github.com/SimonRit/RTK/blob/b32cffdc6f9d7a432c50023c370ed996a7985b69/applications/rtkforwardprojections/rtkforwardprojections.cxx
> for forward projections. However, I find this very hard to read since I
> don’t know C++ so I hope you can help. I have also read
> https://discourse.itk.org/t/itk-to-rtk-migration-forward-projection-issues-questions/2107
> regarding this matter, but I cannot reproduce this.
>
>
>
> What I do in my code:
>
>    - Load 3D CT image and change it to float Image Type (
>    sitk.sitkFloat32)
>    - Define the image type I want for DRR
>    - Create a stack of empty projection images using ConstantImageSource,
>    using ImageType of DRR
>    - Define the RTK geometry object (ThreeDCircularProjectionGeometry)
>    - Then create JosephForwardProjectionImageFilter
>
> Here, the problems arise.
>
> The filter only seems to accept 3D imageTypes, so I have given it the
> [itk.F,3], [itk.F,3] input.
>
> Then, if I try to setInput(0) with constantImageSource.GetOutput(), it
> throws the error:
>
>
>
> Expecting argument of type itkImageF3 or itkImageSourceIF3.
>
> Additional information:
>
> Wrong number or type of arguments for overloaded function 'itkImageToImageFilterIF3IF3_SetInput'.
>
>
>
> Probably since I have used the ImageType of the DRR (2D) for the
> constantImageSource.
>
> Additionally, it seems that the setInput(1) = CT also doesn’t work, maybe
> because I have used sitk.sitkFloat32) to get it to a float image?
>
>
>
> Does anyone know what I have to do here to solve this? And/or is there
> some python script available which does forward projection so I can have a
> look how it is done in Python?
>
>
>
> I also don’t know what I would have to do after this step. Do I just do an
> FDK reconstruction etc. just as in the FirstReconstruction.py script?
>
>
>
> Thank you in advance for anyone who can answer (some) of my questions.
>
>
>
> M
>
>
>
>
>
>
>
>
>
>
>
> My code:
>
> import sys
>
> import itk
>
> from itk import RTK as rtk
>
> import SimpleITK as sitk
>
>
>
> # Utility method that either downloads data from the Girder repository or
>
> # if already downloaded returns the file name for reading from disk
> (cached data).
>
> %run update_path_to_download_script
>
> from downloaddata import fetch_data as fdata
>
>
>
> # Always write output to a separate directory, we don't want to pollute
> the source directory.
>
> import os
>
> OUTPUT_DIR = 'Output_test'
>
>
>
> # Loading 3D CT image
>
> CT = sitk.ReadImage(fdata("training_001_ct.mha"), sitk.sitkFloat32)
>
>
>
> # Defines the image type CT
>
> Dimension_CT = 3
>
> PixelType = itk.F
>
> ImageType_CT = itk.Image[PixelType, Dimension_CT]
>
>
>
> # Defines the image type DRR
>
> Dimension_DRR = 2
>
> ImageType_DRR = itk.Image[PixelType, Dimension_DRR]
>
>
>
> # Create a stack of empty projection images
>
> ConstantImageSourceType = rtk.ConstantImageSource[ImageType_DRR]
>
> constantImageSource = ConstantImageSourceType.New()
>
> # Define origin, sizeOutput and spacing (still need to change these)
>
> origin = [ -127, -127]
>
> sizeOutput = [ 512, 512]
>
> spacing = [ 0.653595, 2.0]
>
> constantImageSource.SetOrigin( origin )
>
> constantImageSource.SetSpacing( spacing )
>
> constantImageSource.SetSize( sizeOutput )
>
> constantImageSource.SetConstant(0.)
>
>
>
> # Defines the RTK geometry object
>
> geometry = rtk.ThreeDCircularProjectionGeometry.New()
>
> numberOfProjections = 360
>
> firstAngle = 0.
>
> angularArc = 360.
>
> sid = 600 # source to isocenter distance
>
> sdd = 1200 # source to detector distance
>
> for x in range(0,numberOfProjections):
>
>   angle = firstAngle + x * angularArc / numberOfProjections
>
>   geometry.AddProjection(sid,sdd,angle)
>
>
>
> # Writing the geometry to disk
>
> xmlWriter = rtk.ThreeDCircularProjectionGeometryXMLFileWriter.New()
>
> xmlWriter.SetFilename ( sys.argv[2] )
>
> xmlWriter.SetObject ( geometry )
>
> xmlWriter.WriteFile()
>
>
>
> REIType = rtk.JosephForwardProjectionImageFilter[ImageType_CT,
> ImageType_CT]
>
> rei = REIType.New()
>
> semiprincipalaxis = [ 50, 50, 50]
>
> center = [ 0, 0, 10]
>
> # Set GrayScale value, axes, center...
>
> #rei.SetDensity(2)
>
> #rei.SetAngle(0)
>
> #rei.SetCenter(center)
>
> #rei.SetAxis(semiprincipalaxis)
>
> rei.SetGeometry( geometry )
>
> rei.SetInput(0, constantImageSource.GetOutput())
>
> rei.SetInput(1, CT)
>
>
>
>
>
> _______________________________________________
> 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: <http://www.creatis.insa-lyon.fr/pipermail/rtk-users/attachments/20210824/28848874/attachment.htm>