SPECTGPU:starting point; created SPECTGPU_projector folder and schel… by danieldeidda · Pull Request #1687 · UCL/STIR

danieldeidda · 2026-03-04T17:02:27Z

…eton classe, atm a lot of copy from niftyGPU

Changes in this pull request

Testing performed

Related issues

Checklist before requesting a review

[] I have performed a self-review of my code
[] I have added docstrings/doxygen in line with the guidance in the developer guide
[] I have implemented unit tests that cover any new or modified functionality (if applicable)
[] The code builds and runs on my machine
[] documentation/release_XXX.md has been updated with any functionality change (if applicable)

Contribution Notes

Please tick the following:

The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in STIR (the Work) under the terms and conditions of the Apache-2.0 License.
I (or my institution) have signed the STIR Contribution License Agreement (not required for small changes).

…eton classe, atm a lot of copy from niftyGPU

…calls the kernels; we need to double check that the kernel can write into the sino but probably we need to dodevice_to_host(stir_sino,cuda_array_created_by_forward)

… that calls the kernels; we need to double check that the kernel can write into the sino but probably we need to dodevice_to_host(stir_sino,cuda_array_created_by_forward)

KrisThielemans · 2026-03-23T15:36:00Z

+        float* dev_image;
+        cudaMalloc(&dev_image, stir_image.size_all() * sizeof(float));


could use

CuVec<float> dev_image(stir_image.size_all());

then you don't need to free at the end. If so, you can pass it to a kernel as dev_image.data().

Also, you probably want to have this out of the loop.

ignoring Cuvec for now

addition of wrappers for actual_FP/BP removal of most of helper currently FP and BP work pretty well and are very close to spectUB rough speed estimate on debug mode were: SPECTUB 3min36s SPECTGPU 20s 3x3x3 gaussian kernel was fine for FP but introduced holes artifact in BP it was necessary to change the kernel to a 5x5x7 with sigma 2 to make a BPimage very close to SPECTUB. Speed above is related to 5x5x7.

codacy-production · 2026-09-09T17:26:28Z

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 65 complexity · 15 duplication

Metric Results

Complexity 65

Duplication 15

View in Codacy

_{NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer}
_{TIP This summary will be updated as you push new changes.}

danieldeidda · 2026-09-10T10:18:22Z

the forward projction of a 150x150x200 image into a sino of 150, 80, 200 (tg,view,ax) values are pretty much identical.

the back projection on the other hand shows some scaling factor which I suppose we expect?

adapted PairProjectorByBinUsingSPECTGPU to be able to run reconstruction

danieldeidda · 2026-09-14T09:16:14Z

Very shallow tests on Attenuation and PSF but is seems that it is doing what I would expect it to be doing. But of course more tests are needed. Actual unit test are also needed like adjointness, etc.

reconstruction also seem to work fairly well:
the following is a reconstruction of simulated data (no noise) 2 iteration with

PSF+AC (PSF is just random to see the effect on recon and speed)
AC
no corrections

currently running with OMP off. when turning OMP on I get an assert aerror on ptr_access

KrisThielemans · 2026-09-14T09:31:28Z

Hi. looks like great progress! Obviously, I'm interested to know some timings :-)

the back projection on the other hand shows some scaling factor which I suppose we expect?

no, we don't. The SPECTUB projector backprojector is adjoint (as it uses the matrix), so has the correct scale.

danieldeidda · 2026-09-14T10:53:55Z

I thought there was some scaling in SPECTUB as I always need to multiply by the number of views to get images in the right scale. Anyway in debug mode with no OMP around 300s 2 iterations fully 3D PSF and AC I am reconstructing with SPECTUB now but it's been running for a while and not finished the first iter yet. GPU recon with no correction is 90 s and 160s only AC

danieldeidda · 2026-09-14T14:00:45Z

comparison with SPECTUB

I though SPECTUB would have some sort of scaling factor as we always need to multiply by the number of views to get "quantitative" values interestingly enough the reconstruction has the same scale.

Rough times in debug mode with OMP=off 2 iterations are:
GPU:
PSFAC 300s (CYAN)
AC 150s
no corr 90s
SPECTUB:
PSFAC 3400s (purple)

Added Hei Yin Jowett Chan as an author to the file.

pdate authorship

Jowett: Revert authorship added to the wrong file

danieldeidda · 2026-09-14T15:22:49Z

when OMP is on: the following debugging track is followed:

`array_to_device(elemT* dev_data, const Array<num_dimensions, elemT>& stir_array)
{
if (stir_array.is_contiguous())

{
  info("array_to_device contiguous", 100);

  cudaMemcpy(dev_data, stir_array.get_const_full_data_ptr(), stir_array.size_all() * sizeof(elemT), cudaMemcpyHostToDevice);

  **stir_array.release_const_full_data_ptr();**

}`

the assert happens when calling realease_const_full_data_ptr() which is called by SPECTGPUBackwardProjectorCUDA.cu: array_to_device(dev_image, stir_image);

KrisThielemans · 2026-09-14T16:11:41Z

interesting. only with OMP?

I think it's indeed the case that the get/release pair currently is not thread-safe, i.e. would have to be in a critical section to make it safe. I don't know about cudaMemCpy and thread-safety.

There's various issues here:

do we want to make it thread-safe?
does it matter for const access?
why are we copying the same image multiple times to the GPU? (That seems wrong/inefficient).

I guess the last one is the main one.

danieldeidda · 2026-09-15T11:20:47Z

Yes that's it! I think that is because the image needs to be updated as the function is called for every view. But then I should copy outside of this call maybe in start_accumulating and free in get_output? so that the BP just updates the output image

KrisThielemans · 2026-09-15T12:38:33Z

Definitely try to do as much of allocation/transfer as you can upfront, or at the back, e.g. set_input() for the forward projector (transfer image and atten image), and start_accumulating_in_new_target/get_output() for the backprojector. Unfortunately, we currently don't have an "end" function for the forward projector (we should!). It will also be a lot faster.

For the forward projector, you'd still run into trouble I think with the output viewgrams, not sure.

For this reason, the parallelprojand NiftyPET interface do all the work upfront/at the end. This however creates trouble as well, as it means we'll be forward projecting all of the data, even if we only do a subset. There's considerably overhead there.
This needs serious thought and re-design. Maybe not for this PR though...

fix problem with OMP by allocating images outside of the kernel

danieldeidda · 2026-09-17T14:58:50Z

Update on release comparison (8 subsets, 16 iterations):
SPECTUB with 3DPSF 1.3 sigma (views in cache out of memory with more) 250 s
SPECTGPU with 3DPSF 3sigmas 101 s
if views are not in cache with 3 sigmas SPECTUB take more than 1h and 30min

codacy-production · 2026-09-17T16:39:52Z

Up to standards ✅

🟢 Issues 0 issues

Results:
0 new issues

View in Codacy

🟢 Metrics 65 complexity · 15 duplication

Metric Results

Complexity 65

Duplication 15

View in Codacy

_{NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer}
_{TIP This summary will be updated as you push new changes.}

KrisThielemans

Looks good, but some comments (mostly house-keeping)

add SPECTGPU in doxygengroups.h
indeces->indices
run pre-commit
I strongly recommend to use CuVeC for all your pointers. No more need for cudaFree, no memory leaks.
It'd be nice to have more const int var = .... could even speed things up (probably not though)
If you copied code from NiftyRec, you need to adapt your copyright header.

KrisThielemans · 2026-09-17T16:51:10Z

+void copy_im_to_stir(
+    DiscretisedDensity<3,float>& image,
+    const float* dev_image);
+
+void copy_stir_im_to_dev(
+    float* dev_image,
+    const DiscretisedDensity<3,float>& image);


there's no reason for these. Just directly call array_to_host etc.

KrisThielemans · 2026-09-17T16:51:35Z

+
+void run_backward_projection_cuda(
+        const RelatedViewgrams<float>& stir_sino,
+        float* dev_image,


This is the output arg,so in STIR convention it needs to come first.

KrisThielemans · 2026-09-17T16:52:41Z

+  std::vector<float> _crs;
+  std::vector<short> _s2c;


I'd prefer less cryptic variable names

KrisThielemans · 2026-09-17T16:53:13Z

+    See STIR/LICENSE.txt for details
+*/
+
+#pragma once


surprising, but this is not portable. Use the usual #ifdef stuff

KrisThielemans · 2026-09-17T16:53:25Z

+
+#pragma once
+
+#include <cuda_runtime.h>


doesn't seem to be needed to include it here

KrisThielemans · 2026-09-17T16:55:19Z

+
+      if (o1!=o2)
+      {
+          std::cout << "image origin: "


can't use std::cout.I suppose put them all in a string, but better to use has_same_characteristics()

KrisThielemans · 2026-09-17T16:55:46Z

+      std::fabs(this->spacing_x - tg_spacing> 1e-5f) )
+  {
+      error(
+          "SPECTGPU: expected axial and tangential dimensions/spacings "


use stir/format

KrisThielemans · 2026-09-17T16:56:52Z

+      _slope(-1),
+      _sigma0(-1),


strange defaults. might be ok, but dangerous?

KrisThielemans · 2026-09-17T16:58:32Z

+//    array_to_device(dev_image, stir_image);
+
+
+    float3 spacing = make_float3(spacing_x,


I find it somewhat confusing to switch from STIR's z,y,x to x,y,z

KrisThielemans · 2026-09-17T16:59:26Z

+//            << "  stir view = " << vg.get_view_num()
+//            << std::endl;
+        //the following sign is introduced to match SPECTUB
+        float angle_rad = -vg.get_view_num() * 2.f * M_PI / num_views;


I would vastly prefer to use the STIR function to get the angle (vg.get_phi() I guess)

danieldeidda and others added 9 commits March 4, 2026 17:01

SPECTGPU:starting point; created SPECTGPU_projector folder and schel…

b06adca

…eton classe, atm a lot of copy from niftyGPU

forward line integrals

c948b16

$JCHAN\jchan$

implemented a backward projection kernel in cuda

5d300cf

GPUSPECT: forward rotation and interpolation in cuda

d4cd638

Adding kernel push

9ca0e11

fixup! Adding kernel push

d731b5c

SPECTGPU: run pre-commit hooks on SPECTGPURotateAndGaussianInterpolate

8dfbdd1

joining all contributions; created an actual_forward_project() that …

aec9adf

…calls the kernels; we need to double check that the kernel can write into the sino but probably we need to dodevice_to_host(stir_sino,cuda_array_created_by_forward)

fixup! joining all contributions; created an actual_forward_project()…

edc60a8

… that calls the kernels; we need to double check that the kernel can write into the sino but probably we need to dodevice_to_host(stir_sino,cuda_array_created_by_forward)

casperdcl moved this to In Progress in 2026-03 AIRBI Hackathon Mar 16, 2026

casperdcl added this to 2026-03 AIRBI Hackathon Mar 16, 2026

KrisThielemans reviewed Mar 23, 2026

View reviewed changes

KrisThielemans marked this pull request as draft March 23, 2026 15:46

KrisThielemans self-assigned this Mar 23, 2026

danieldeidda added 2 commits September 9, 2026 18:21

danieldeidda marked this pull request as ready for review September 9, 2026 17:28

SPECTGPU added Attenuation and depth dependent PSF;

74a840c

adapted PairProjectorByBinUsingSPECTGPU to be able to run reconstruction

jchan192 added 3 commits September 14, 2026 10:48

Update authorship in SPECTGPURotateAndGaussianInterpolate.cu

e7cbf28

Added Hei Yin Jowett Chan as an author to the file.

Update SPECTGPUProjection.cu

5d162c3

pdate authorship

Update SPECTGPURotateAndGaussianInterpolate.cu

864c871

Jowett: Revert authorship added to the wrong file

update contributors;

8f0ab94

fix problem with OMP by allocating images outside of the kernel

relax condition of spacing

f79f532

KrisThielemans requested changes Sep 17, 2026

View reviewed changes

		float* dev_image;
		cudaMalloc(&dev_image, stir_image.size_all() * sizeof(float));

		// array_to_device(dev_image, stir_image);


		float3 spacing = make_float3(spacing_x,


		#pragma once

		#include <cuda_runtime.h>

Conversation

danieldeidda commented Mar 4, 2026

Changes in this pull request

Testing performed

Related issues

Checklist before requesting a review

Contribution Notes

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

codacy-production Bot commented Sep 9, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Up to standards ✅

Uh oh!

danieldeidda commented Sep 10, 2026

Uh oh!

danieldeidda commented Sep 14, 2026

Uh oh!

KrisThielemans commented Sep 14, 2026

Uh oh!

danieldeidda commented Sep 14, 2026

Uh oh!

danieldeidda commented Sep 14, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

danieldeidda commented Sep 14, 2026

Uh oh!

KrisThielemans commented Sep 14, 2026

Uh oh!

danieldeidda commented Sep 15, 2026

Uh oh!

KrisThielemans commented Sep 15, 2026

Uh oh!

danieldeidda commented Sep 17, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

codacy-production Bot commented Sep 17, 2026

Up to standards ✅

Uh oh!

KrisThielemans left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

codacy-production Bot commented Sep 9, 2026 •

edited

Loading

danieldeidda commented Sep 14, 2026 •

edited

Loading

danieldeidda commented Sep 17, 2026 •

edited

Loading