Skip to content

SPECTGPU:starting point; created SPECTGPU_projector folder and schel… - #1687

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

danieldeidda wants to merge 16 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 on lines +123 to +124
float* dev_image;
cudaMalloc(&dev_image, stir_image.size_all() * sizeof(float));

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.

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.

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

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

ignoring Cuvec for now

@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

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

5 participants