Skip to content

SPECTGPU - #1687

Open
danieldeidda wants to merge 39 commits into
UCL:masterfrom
danieldeidda:GPUSPECT
Open

danieldeidda wants to merge 39 commits into
UCL:masterfrom
danieldeidda:GPUSPECT

Conversation

@danieldeidda

Copy link
Copy Markdown
Collaborator

…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).

danieldeidda and others added 9 commits March 4, 2026 17:01
…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)
Comment thread src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx Outdated
@KrisThielemans
KrisThielemans marked this pull request as draft March 23, 2026 15:46
@KrisThielemans KrisThielemans self-assigned this Mar 23, 2026
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.
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

codacy-production Bot commented Sep 9, 2026 •

Copy link
Copy Markdown

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
danieldeidda marked this pull request as ready for review September 9, 2026 17:28
@danieldeidda

Copy link
Copy Markdown
Collaborator Author
image 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?
image

adapted PairProjectorByBinUsingSPECTGPU to be able to run reconstruction
@danieldeidda

Copy link
Copy Markdown
Collaborator Author

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

  1. PSF+AC (PSF is just random to see the effect on recon and speed)
  2. AC
  3. no corrections
image image

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

@KrisThielemans

Copy link
Copy Markdown
Collaborator

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

Copy link
Copy Markdown
Collaborator Author

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

danieldeidda commented Sep 14, 2026 •

Copy link
Copy Markdown
Collaborator Author

comparison with SPECTUB
image

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.
image

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.
Jowett: Revert authorship added to the wrong file
@danieldeidda

Copy link
Copy Markdown
Collaborator Author

when OMP is on: the following debugging track is followed:
image
`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

Copy link
Copy Markdown
Collaborator

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

Copy link
Copy Markdown
Collaborator Author

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

Copy link
Copy Markdown
Collaborator

See also #1768 (comment)

@casperdcl

casperdcl commented Oct 7, 2026 •

Copy link
Copy Markdown
Collaborator

The .cu side is responsible for alloc (.resize(N)) so should also do dealloc (.resize(0)) on completion. This would make the .cxx destructors effectively do nothing.

@KrisThielemans

Copy link
Copy Markdown
Collaborator

yeah, but..

  • how do they know when to do this? (They could do it at the end of every run, but that defeats the purpose of having them in the class anyway, as it will re/deallocate at every run)
  • have a "allocate" and "deallocate` functions called in the .cxx, implemented in the .cu, which I guess is what @danieldeidda originally had anyway.

Either way, there seems to be little advantage now of having cuvec in the .h. Maybe it's still more exception safe (although it would depend on where the exception is thrown, and I suppose that'd call the .cxx destructor anyway. On the other hand, we've now introduced an easy way to do things that will crash, which isn't good.

As it caused problems for @vcosmi as well (he needed to guard with STIR_WITH_CUDA anyway), I feel like reverting cudautilities.h back again to only import cuvec if STIR_WITH_CUDA. (Of course, a user could still do what we attempted).

Unless @casperdcl has a safe way of doing this after all.

… a AllocatedStack structure was created:

- all the allocations happen at setup so that in the loop we only initialise and assigne values. this has significantly speeded up performonce.
- all arrays are now CuVec but at the end of the reconstruction we still have a free problem. The structure has a sestructor that calls resize(0) for all the cuvecs but the error is still there
… a AllocatedStack structure was created:

- all the allocations happen at setup so that in the loop we only initialise and assigne values. this has significantly speeded up performonce.
- all arrays are now CuVec but at the end of the reconstruction we still have a free problem. The structure has a sestructor that calls resize(0) for all the cuvecs but the error is still there
- the new stack structure create some race problem when using OMP because many threads where accessing the same object, so a thread_local AllocatedStack was used for everythread
in the BP part given that each thread contribute to the the same image it was also necessary to add an accumulation step as critical section, so a kernel for arrays accumulation was created;
updated timing given that spectUB becomes really fast once the matrix is created:
SPECTUB 2 iter 231 s (223 matrix + 4.2 s per full iter)
SPECTGPU 2 iter: 8.7s (6.4 matrix + 2.3 per full iter)
@codacy-production

codacy-production Bot commented Oct 8, 2026 •

Copy link
Copy Markdown

Not up to standards ⛔

🔴 Issues 2 high · 1 medium · 3 minor

Alerts:
⚠ 6 issues (≤ 0 issues of at least minor severity)

Results:
6 new issues

Category Results
BestPractice 1 medium
ErrorProne 2 high
CodeStyle 3 minor

View in Codacy

🟢 Metrics 98 complexity · 25 duplication

Metric Results
Complexity 98
Duplication 25

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

Copy link
Copy Markdown
Collaborator Author

I guess I will put a comment because from the commit it's a bit weird to read:
To make things faster and given that all arrays are of the same size a AllocatedStack structure was created:

  • all the allocations happen at setup so that in the loop we only initialise and assigne values. this has significantly speeded up performonce.
  • all arrays are now CuVec but at the end of the reconstruction we still have a free problem. The structure has a sestructor that calls resize(0) for all the cuvecs but the error is still there
  • the new stack structure create some race problem when using OMP because many threads where accessing the same object, so a thread_local AllocatedStack was used for everythread
    in the BP part given that each thread contribute to the the same image it was also necessary to add an accumulation step as critical section, so a kernel for arrays accumulation was created;
    updated timing given that spectUB becomes really fast once the matrix is created:
    SPECTUB 2 iter 231 s (223 matrix + 4.2 s per full iter)
    SPECTGPU 2 iter: 8.7s (6.4 matrix + 2.3 per full iter)

@casperdcl maybe have a look at my resize(0) destructor to see if I misinterpreted your comment.

creating OSMAPOSL and PpoissonLogLikelihood tests (these are running with the caveat that i had to make epsilon for the hessian product with input 0.001)
Comment thread src/recon_buildblock/SPECTGPU_projector/SPECTGPUForwardProjectorCUDA.cu Outdated

@KrisThielemans KrisThielemans left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

some changes in cuda_utilites for cuvec, which you're no longer using of course, but best to add them in anyway

Comment thread src/buildblock/cuda_utilities.cu
Comment thread src/buildblock/cuda_utilities.cu
Comment thread src/buildblock/cuda_utilities.cu
Comment thread src/buildblock/cuda_utilities.cu
Comment thread src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx Outdated

@KrisThielemans KrisThielemans left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

some comments. mostly house-keeping.

I'm still a bit puzzled by the threading and GPU interaction, but need to do other things...

Comment thread src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx Outdated
Comment thread src/recon_test/test_PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx Outdated
Comment thread src/recon_buildblock/SPECTGPU_projector/SPECTGPURotateAndGaussianInterpolate.cu Outdated
Comment thread src/recon_buildblock/SPECTGPU_projector/SPECTGPUPSFGaussian.cu Outdated
Comment thread src/recon_buildblock/SPECTGPU_projector/SPECTGPUProjection.cu Outdated
Comment thread src/recon_buildblock/SPECTGPU_projector/SPECTGPUForwardProjectorCUDA.cu Outdated
Comment thread src/recon_buildblock/SPECTGPU_projector/SPECTGPUForwardProjectorCUDA.cu Outdated
Comment thread src/recon_buildblock/SPECTGPU_projector/SPECTGPUForwardProjectorCUDA.cu Outdated
Comment thread src/recon_buildblock/SPECTGPU_projector/ProjectorByBinPairUsingSPECTGPU.cxx Outdated
general clean up;
moving cu_allocate and cu_free into  cuda utilities;
adressing all comments
@danieldeidda

Copy link
Copy Markdown
Collaborator Author

@KrisThielemans I think I am done

@KrisThielemans

Copy link
Copy Markdown
Collaborator

I think I am done

all tests are failing, so I guess not yet :-/

did you add the STIR_WITH_CUDA to the test?

@KrisThielemans KrisThielemans left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we have src/buildblock/cuda_utilities.cu and src/buildblock/CUDA/cuda_utilities.cu

Comment thread examples/samples/OSMAPOSL_osem_SPECTUB.par Outdated
Comment thread recon_test_pack/SPECT/SPECTGPU/org/FBP.hv
Comment thread recon_test_pack/SPECT/SPECTGPU/FBP2D.par Outdated
Comment thread recon_test_pack/SPECT/SPECTGPU/run_SPECTGPU_tests.sh
Comment thread src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h Outdated
Comment thread src/include/stir/recon_buildblock/SPECTGPU_projector/BackProjectorByBinSPECTGPU.h Outdated
Comment thread src/recon_buildblock/SPECTGPU_projector/ForwardProjectorByBinSPECTGPU.cxx Outdated
@danieldeidda

Copy link
Copy Markdown
Collaborator Author

I think I am done

all tests are failing, so I guess not yet :-/

did you add the STIR_WITH_CUDA to the test?

I thought i did bu maybe not everywher. will need to check

@KrisThielemans

Copy link
Copy Markdown
Collaborator

Compiling fine, but SPECTUB test segfaulting. Do you have the same problem when running it locally?

@danieldeidda

danieldeidda commented Oct 19, 2026 •

Copy link
Copy Markdown
Collaborator Author

Compiling fine, but SPECTUB test segfaulting. Do you have the same problem when running it locally?

I actually never run it till the end, because is quite slow but running now and getting some difference in the gradient test

@danieldeidda

Copy link
Copy Markdown
Collaborator Author

Compiling fine, but SPECTUB test segfaulting. Do you have the same problem when running it locally?

it is giving seg fault locally too

@danieldeidda

Copy link
Copy Markdown
Collaborator Author

i found why the segfault happens but the issue will still be that the gradient test does not pass

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

6 participants