[Rtk-users] Backprojection in 2D/3D and adjointness concerns

[Simon Rit]

Dear Laura,
I don't remember what was NormalizedJoseph but it's been removed a while
ago:
https://github.com/RTKConsortium/RTK/pull/171/commits/892e5cc5c765a13a7205ad2cabeb6193a75631b5
To answer your questions:
1/ Possibly. The adjoint of CUDAForwardProjectionImageFilter is
CudaRayCastBackProjectionImageFilter. The adjoint of
JosephBackProjectionImageFilter is JosephFowardProjectionImageFilter.
2/ Good!
3/ This is not the expected behavior and I have opened an issue for this:
https://github.com/RTKConsortium/RTK/issues/670
It would be great if you have time to check what's going on. The forward
projection code is here:
https://github.com/RTKConsortium/RTK/blob/master/src/rtkCudaForwardProjectionImageFilter.cu#L61-L153
The backprojection code is here:
https://github.com/RTKConsortium/RTK/blob/master/src/rtkCudaRayCastBackProjectionImageFilter.cu#L63-L238
The trickiest part is that the wired interpolation in the forward
projection has to be transposed to "splat" the values.
Let us know if you spot something wrong. We'll try to investigate if you
don't manage to do it.
Thanks for the report in any case!
Simon


On Thu, Jan 9, 2025 at 1:19 PM Laura Hellwege <l.hellwege at uni-luebeck.de>
wrote:

> Hello!
>
>
>
> I have some questions related to this post from 2015:
>
> https://public.kitware.com/pipermail/rtk-users/2015-February/009750.html
>
>
>
> First, the referenced website “
> http://www.openrtk.org/Doxygen/classrtk_1_1NormalizedJosephBackProjectionImageFilter.html”
> seems to not exist anymore.
>
>
>
> Second, I have some issues concerning backprojection and adjointness:
>
> For the following I am using these parameters:
>
> sid = 500
>
> sdd = 1000
>
> n_projections = 4
>
> angular_increment = 90.0
>
> I am using the python version of ITK and tried forward projection
> (CUDAForwardProjectionImageFilter) followed by backprojection
> (JosephBackProjectionImageFilter/CudaBackProjectionImageFilter/CudaRayCastBackProjectionImageFilter)
> of an image with a single pixel value set to one:
>
> img = np.zeros(image_size)
>
> img[image_size[0] // 2, image_size[1] // 2, image_size[2] // 2] = 1.0
>
>
>
>    1. I do not get expected results for the 2D case:
>
> image_size = [3, 1, 3]
>
> pixel_size_mm = 1.0
>
> slice_thickness_mm = 1.0
>
>
>
> projection_size = [5, 1]
>
> projection_spacing = [1.0, 1.0]
>
>
>
> After forward projection, the backprojected image by
> JosephBackProjectionImageFilter is all *zeros*. Is this expected behavior?
>
>
>
>    2. According to the definition of adjoint operator, the following
>    should hold:
>
> <Ax, y> = <x, adjoint(A)y>            (with <.,.> indicating the scalar
> product)
>
>
>
> This works for 3D *CPU* case (JosephBackProjectionImageFilter):
>
> image_size = [3, 3, 3]
>
> pixel_size_mm = 1.0
>
> slice_thickness_mm = 1.0
>
>
>
> projection_size = [5, 5]
>
> projection_spacing = [1.0, 1.0]
>
>
>
> <Ax, y> = <fp, fp> = 9.0
>
> <x, adjoint(A)y> = <img, bp> = 9.000001907348633
>
> This discrepancy is within tolerance.
>
>
>
>    3. However, in the 3D *GPU* case:
>
> <x, adjoint(A)y> = <img, bp> = 4.0
>                                                         for
> CudaBackProjectionImageFilter
>
> and
>
> <x, adjoint(A)y> = <img, bp> = 0.33841705322265625                for
> CudaRayCastBackProjectionImageFilter
>
>
>
> What kind of filtering or normalization has been used on the
> projections/images?  I could not verify adjointment with the aforementioned
> method of backprojecting a sinogram of ones and applying division.
>
> Is there a GPU backprojection filter for which adjointness holds?
>
>
>
> Thank you in advance for your reply and suggestions.
>
> Best regards, Laura
>
>
> _______________________________________________
> 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/20250110/43584db9/attachment-0001.htm>