Repository navigation
Conversation
This adds an AMD GPU build of TIGRE's Python backend with ROCm/HIP, alongside the existing CUDA build. The CUDA path is selected by default and the AMD support is added behind USE_HIP / BUILD_WITH_HIP, so every effort has been made to leave the CUDA build untouched; the ROCm build is opt-in through environment variables. The CUDA sources were compile-checked with nvcc but were not run, as the development hosts have no NVIDIA GPU.
How to review: setup.py first (the build wiring), then Common/CUDA/cuda_to_hip.h (the CUDA->HIP symbol mapping and the texture path), then the per-source changes (texture filter mode, the Siddon ray-loop bound, and the cuRAND -> hipRAND include), then the install-doc update.
setup.py keeps the existing hand-rolled per-.cu compiler driver and adds a BUILD_WITH_HIP branch (both the Unix and the MSVC/Windows paths). With BUILD_WITH_HIP=1 it compiles the .cu sources with hipcc (-x hip --offload-arch=$HIP_ARCH, default gfx90a; comma-separated for multiple targets), links amdhip64 instead of cudart, and skips the nvcc-at-import probing that is absent on a ROCm-only host. On Windows it drives clang++ directly with --hip-path and -fms-runtime-lib=dll to match the MSVC /MD CRT.
The main porting concern is the texture path. The backprojection kernels and the interpolated forward projector sample a 3D float array with trilinear filtering (cudaFilterModeLinear + cudaReadModeElementType). Whether AMD hardware can do that fetch is architecture-dependent: gfx90a (CDNA2) rejects creation of such a texture, while RDNA (gfx1100, gfx1201) accepts it and filters in hardware. Rather than hard-code one path, cuda_to_hip.h decides at runtime: a one-time self-test creates the real texture configuration over a small ramp and confirms a known sample actually interpolates -- not merely that creation succeeded, since some hardware accepts the texture but silently point-samples, which would ship wrong results. When supported, the texture is created Linear and tex3D_TIGRE forwards to the hardware fetch; otherwise it is created Point and tex3D_TIGRE interpolates in software, point-sampling the eight neighbours and lerping with CUDA's unnormalized -0.5 texel-center convention (out-of-array neighbours read 0, matching cudaAddressModeBorder). The texture's creation mode and the sampling path are driven from the single cached verdict so they cannot disagree. On CUDA the runtime helper reports supported without testing and the AMD-only code is excluded by the preprocessor, so the original hardware-filtered textures are used.
A separate robustness fix bounds the Siddon ray-marching loop. The loop length Np is computed from index bounds chosen by exact float-equality tests on quantities formed with __fdividef, an approximate division whose result can differ slightly between back ends; a flipped equality test could select an out-of-range bound and make Np astronomical (a data-dependent hang). A straight ray crosses at most Nx+Ny+Nz voxel-boundary planes, so Np is now capped at that geometric maximum. This cap is the one change on a shared code path; it is a no-op for every valid ray and on CUDA.
RandomNumberGenerator.cu uses device-side cuRAND, mapped to hipRAND in the compat header; the directed-rounding intrinsics __fsqrt_rd / __frcp_rd (absent in HIP) map to the round-to-nearest variants, immaterial for an nRMSE/adjointness-graded reconstruction.
This work was authored with the assistance of Claude, an AI assistant by Anthropic.
Test Plan:
Linux, gfx90a (AMD Instinct MI250X), ROCm 7.2.1 -- software trilinear fallback (hardware linear fp32 textures rejected on CDNA2):
export BUILD_WITH_HIP=1 ROCM_PATH=/opt/rocm HIP_ARCH=gfx90a
pip install -e . --no-build-isolation
Linux, gfx1100 (RDNA3), ROCm 7.2.1, and Windows 11, gfx1201 (RDNA4), ROCm 7.14 -- hardware linear texture filtering (the self-test reports supported and the hardware fetch is used). Windows drives clang++ directly:
BUILD_WITH_HIP=1 HIPCC=<rocm>/lib/llvm/bin/clang++.exe ROCM_PATH=<rocm> HIP_ARCH=gfx1201
pip install -e . --no-build-isolation
On the 256^3 head phantom (cone geometry): forward project (Siddon and interpolated), back project (matched and FDK), reconstruct with FDK and OS-SART; all outputs finite. Interpolated (trilinear) projection matches Siddon to within 0.5% relative norm; adjointness <Ax(x),y> ~= <x,Atb(y)> relative residual ~1.2e-05; sart/ossart/sirt/cgls/fista and the TV-regularized asd_pocs/awasd_pocs/os_asd_pocs solvers all produce finite reconstructions. Verified with both the software-fallback path (gfx90a) and the hardware-filtering path (gfx1100, gfx1201).
The CUDA build was compile-checked with nvcc (CUDA 13.3, IS_FOR_PYTIGRE, BUILD_WITH_HIP off) -- all changed sources compile; it was not linked or run, as no NVIDIA GPU was available. Other AMD architectures are buildable via HIP_ARCH but were not exercised here.
Member
|
hi @jeffdaily, spectacular PR, thanks a ton! I don't have access to an AMD GPU at the moment, but I will try to get one ASAP and test this in my side to understand all the changes better. |
Closed
# Conflicts: # Common/CUDA/Siddon_projection_parallel.cu
The Linux build appends -std=c++11 to every compile that does not pass its own -std, including the HIP compile of the .cu sources. Current rocRAND device headers (reached through hiprand_kernel.h) use std::enable_if_t, a C++14 feature, so RandomNumberGenerator.cu and the other .cu sources fail on current ROCm with "no template named 'enable_if_t' in namespace 'std'". C++14 is the lowest standard those headers accept. The HIP branch of unix_wrap_compile now appends -std=c++14 when no -std was given, and the existing C++11 default then leaves it alone. The nvcc path and the host C++ sources keep -std=c++11, and the Windows HIP path passes no -std (clang's default is newer), so neither changes. Authored with Claude. Test Plan: Linux, gfx90a (AMD Instinct MI250X), ROCm 10.1: ``` export BUILD_WITH_HIP=1 ROCM_PATH=<rocm> HIP_ARCH=gfx90a python setup.py build_ext --inplace -j 16 cd Python && python -m pytest -q tests/test_compute_cor.py tests/test_demo_angles.py tests/test_denoise_2d_guard.py tests/test_filtering_workers.py tests/test_gpuids_default.py tests/test_krylov_restart_and_norms.py tests/test_krylov_smoke.py ``` All 8 extensions build (with -std=c++11 the .cu sources do not compile on this ROCm), and the 7 test files pass (70 tests). On the 256^3 head phantom (cone geometry, 32 angles) the interpolated projection matches Siddon to 0.40% relative norm, FDK and OS-SART reconstructions are finite with OS-SART below FDK error, and the adjointness relative residual is 1.0e-05; the sart/ossart/sirt/cgls/fista/asd_pocs/awasd_pocs/os_asd_pocs/ossart_tv solvers produce finite reconstructions.
CreateTextureParallel (Siddon_projection_parallel.cu) and CreateTextureParallelInterp (ray_interpolated_projection_parallel.cu) upload the image into the texture array with cudaMemcpy3DAsync on stream[1], and the caller then launches the first projection kernel on stream[0]. Nothing orders the two streams, so the first kernel batch can read the array before the upload has finished.
On CUDA the race is hidden: the image is a pageable host buffer, and an async copy from pageable memory is staged synchronously, so the upload is complete by the time cudaMemcpy3DAsync returns. HIP does not guarantee that, and on an MI250X (gfx90a) the parallel Siddon projector intermittently returned wrong values in projections 0 to 8 (the first PROJ_PER_BLOCK batch), with errors up to 22% of the maximum, in about 1 run in 7.
Both functions now call cudaStreamSynchronize(stream[1]) right after the upload, as the cone projector (cudaDeviceSynchronize at the end of CreateTexture) and the parallel backprojector (cudaStreamSynchronize after CreateTextureParallel) already do. The interpolated case was not observed to fail but has the same ordering. On CUDA the copy has already completed when the sync is reached, so results there are unchanged. The other host-to-device uploads in the built sources are already synchronized against the stream that consumes them.
Authored with Claude.
Test Plan:
Linux, gfx90a (AMD Instinct MI250X), ROCm 10.1:
```
export BUILD_WITH_HIP=1 ROCM_PATH=<rocm> HIP_ARCH=gfx90a
python setup.py build_ext --inplace -j 16
cd Python && python -m pytest -q tests/test_compute_cor.py tests/test_demo_angles.py tests/test_denoise_2d_guard.py tests/test_filtering_workers.py tests/test_gpuids_default.py tests/test_krylov_restart_and_norms.py tests/test_krylov_smoke.py
python -c "
import numpy as np, tigre
from tigre.utilities import sample_loader
geo = tigre.geometry(mode='parallel', nVoxel=np.array([128]*3), default=True)
x = sample_loader.load_head_phantom(geo.nVoxel).astype(np.float32)
ang = np.linspace(0, np.pi, 64, endpoint=False) + 0.3
for method in ('Siddon', 'interpolated'):
ref = tigre.Ax(x, geo, ang, method)
print(method, sum(not np.array_equal(tigre.Ax(x, geo, ang, method), ref) for _ in range(20)), 'of 20 differ')
"
```
The 7 test files pass (70 tests). Repeating Ax 20 times on identical input (head phantom and a random volume, parallel and cone geometry, run three times for Siddon and twice for interpolated) gives bitwise identical projections every time; before the change, parallel Siddon differed in up to 7 of 20 runs. Parallel-beam Siddon matches the interpolated projector to 0.43%, 0.41% and 0.56% at axis-aligned, offset and full Euler angle sets, and the cone-beam checks (interpolated vs Siddon 0.40%, FDK nRMSE 0.286, adjointness residual 1.0e-05) are unchanged.
The HIP port replaced the <cuda.h>, <curand_kernel.h> and <curand.h> includes in RandomNumberGenerator.cu with the compatibility header. That header's HIP branch includes hipRAND, but its CUDA branch includes only the CUDA runtime headers, so nvcc failed on this file (`identifier "curandState" is undefined`, 38 errors) and the _RandomNumberGenerator extension did not build for NVIDIA GPUs. Include <curand_kernel.h> and <curand.h> again when USE_HIP is not defined, in the same place as before, rather than adding cuRAND to the compatibility header: only this file uses cuRAND, so every other CUDA translation unit keeps exactly the headers it had. <cuda.h> still comes from the compatibility header's CUDA branch. The HIP preprocessed output of this file is unchanged apart from blank lines. Authored with Claude. Test Plan: CUDA (nvcc), from the repository root: ``` python setup.py build_ext --inplace ``` HIP, from the repository root: ``` BUILD_WITH_HIP=1 HIP_ARCH=gfx90a python setup.py build_ext --inplace cd Python && python -m pytest -q tests/test_compute_cor.py tests/test_demo_angles.py tests/test_denoise_2d_guard.py tests/test_filtering_workers.py tests/test_gpuids_default.py tests/test_krylov_restart_and_norms.py tests/test_krylov_smoke.py ```
Contributor
Author
|
Note: this reply was drafted by an AI assistant. The branch as I previously pushed it did not compile with nvcc: RandomNumberGenerator.cu had lost its cuRAND includes on the CUDA path. That is fixed; the CUDA side of that file includes the cuRAND headers again, in the same place as master. I also merged current master into this branch, keeping the new parallel Siddon projector from master where it conflicted with the port, and added a small setup.py fix so the HIP build works with current ROCm, whose rocRAND headers need C++14. The branch also fixes an upload/kernel ordering race in the parallel-beam projectors: on an MI250X, 28 of 120 repeated parallel Siddon projections of identical input differed, with errors up to 44% of the maximum, and with the fix every repeat is bitwise identical. I re-ran the cone and parallel-beam projection, reconstruction and adjointness checks on an AMD Instinct MI250X, an AMD Radeon PRO W7800 and an AMD Radeon 8060S (Windows), and checked that all 8 extensions build with nvcc 12.8. |
maxdistanceCuboid() is host code in a .cu file and calls unqualified max() on floats. CUDA declares host min/max overloads for float and double, but HIP's host math header only declares max(int, int), so under hipcc the float arguments are silently converted to int. In the 3D-angle branch (theta or psi nonzero) the largest voxel size becomes an integer: for voxel sizes above 1 the interpolated projector (cone and parallel) starts sampling too far along the ray and skips part or all of the object, and for voxel sizes below 1 it divides by zero. With dVoxel 1.5, theta 0.3 and psi 0.2 the interpolated projection of a 64^3 phantom is all zeros where Siddon gives a sum of 2.1e5. Call fmaxf, which is what CUDA's max(float, float) resolves to, at the four float max() calls in the function. CUDA results, NaN handling included, are unchanged. The 2D-angle branch gets the same change; its result is floored by the caller, so its output does not change. Device code is unchanged. Add a test comparing the interpolated and Siddon projectors with voxel sizes 0.8 and 1.5 and nonzero theta and psi, for cone and parallel beam. Before this change the two 1.5 cases fail on ROCm with a relative error of 1.0. Authored with Claude. Test Plan: ``` BUILD_WITH_HIP=1 HIP_ARCH=gfx90a python setup.py build_ext --inplace cd Python python -m pytest tests/test_interpolated_voxel_size.py tests/test_compute_cor.py tests/test_demo_angles.py tests/test_denoise_2d_guard.py tests/test_filtering_workers.py tests/test_gpuids_default.py tests/test_krylov_restart_and_norms.py tests/test_krylov_smoke.py ```
Contributor
Author
|
Note: this reply was drafted by an AI assistant. I added a follow-up commit that fixes the interpolated forward projector on AMD GPUs when theta or psi is nonzero and the voxel size is not a whole number: a host max() call cut the voxel size to an integer, so for sizes above 1 part or all of the object was skipped. It also adds a small test that checks the interpolated projector against Siddon for that case. On NVIDIA nothing changes (the call now names fmaxf, which CUDA's max(float, float) already called); the new test and the existing Python tests pass on a V100 before and after. |
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This adds an AMD GPU build of TIGRE's Python backend with ROCm/HIP, alongside the existing CUDA build. The CUDA path is selected by default and the AMD support is added behind USE_HIP / BUILD_WITH_HIP, so every effort has been made to leave the CUDA build untouched; the ROCm build is opt-in through environment variables. The CUDA sources were compile-checked with nvcc but were not run, as the development hosts have no NVIDIA GPU.
How to review: setup.py first (the build wiring), then Common/CUDA/cuda_to_hip.h (the CUDA->HIP symbol mapping and the texture path), then the per-source changes (texture filter mode, the Siddon ray-loop bound, and the cuRAND -> hipRAND include), then the install-doc update.
setup.py keeps the existing hand-rolled per-.cu compiler driver and adds a BUILD_WITH_HIP branch (both the Unix and the MSVC/Windows paths). With BUILD_WITH_HIP=1 it compiles the .cu sources with hipcc (-x hip --offload-arch=$HIP_ARCH, default gfx90a; comma-separated for multiple targets), links amdhip64 instead of cudart, and skips the nvcc-at-import probing that is absent on a ROCm-only host. On Windows it drives clang++ directly with --hip-path and -fms-runtime-lib=dll to match the MSVC /MD CRT.
The main porting concern is the texture path. The backprojection kernels and the interpolated forward projector sample a 3D float array with trilinear filtering (cudaFilterModeLinear + cudaReadModeElementType). Whether AMD hardware can do that fetch is architecture-dependent: gfx90a (CDNA2) rejects creation of such a texture, while RDNA (gfx1100, gfx1201) accepts it and filters in hardware. Rather than hard-code one path, cuda_to_hip.h decides at runtime: a one-time self-test creates the real texture configuration over a small ramp and confirms a known sample actually interpolates -- not merely that creation succeeded, since some hardware accepts the texture but silently point-samples, which would ship wrong results. When supported, the texture is created Linear and tex3D_TIGRE forwards to the hardware fetch; otherwise it is created Point and tex3D_TIGRE interpolates in software, point-sampling the eight neighbours and lerping with CUDA's unnormalized -0.5 texel-center convention (out-of-array neighbours read 0, matching cudaAddressModeBorder). The texture's creation mode and the sampling path are driven from the single cached verdict so they cannot disagree. On CUDA the runtime helper reports supported without testing and the AMD-only code is excluded by the preprocessor, so the original hardware-filtered textures are used.
A separate robustness fix bounds the Siddon ray-marching loop. The loop length Np is computed from index bounds chosen by exact float-equality tests on quantities formed with __fdividef, an approximate division whose result can differ slightly between back ends; a flipped equality test could select an out-of-range bound and make Np astronomical (a data-dependent hang). A straight ray crosses at most Nx+Ny+Nz voxel-boundary planes, so Np is now capped at that geometric maximum. This cap is the one change on a shared code path; it is a no-op for every valid ray and on CUDA.
RandomNumberGenerator.cu uses device-side cuRAND, mapped to hipRAND in the compat header; the directed-rounding intrinsics __fsqrt_rd / __frcp_rd (absent in HIP) map to the round-to-nearest variants, immaterial for an nRMSE/adjointness-graded reconstruction.
This work was authored with the assistance of Claude, an AI assistant by Anthropic.
Test Plan:
Linux, gfx90a (AMD Instinct MI250X), ROCm 7.2.1 -- software trilinear fallback (hardware linear fp32 textures rejected on CDNA2):
Linux, gfx1100 (RDNA3), ROCm 7.2.1, and Windows 11, gfx1201 (RDNA4), ROCm 7.14 -- hardware linear texture filtering (the self-test reports supported and the hardware fetch is used). Windows drives clang++ directly:
On the 256^3 head phantom (cone geometry): forward project (Siddon and interpolated), back project (matched and FDK), reconstruct with FDK and OS-SART; all outputs finite. Interpolated (trilinear) projection matches Siddon to within 0.5% relative norm; adjointness <Ax(x),y> ~= <x,Atb(y)> relative residual ~1.2e-05; sart/ossart/sirt/cgls/fista and the TV-regularized asd_pocs/awasd_pocs/os_asd_pocs solvers all produce finite reconstructions. Verified with both the software-fallback path (gfx90a) and the hardware-filtering path (gfx1100, gfx1201).
The CUDA build was compile-checked with nvcc (CUDA 13.3, IS_FOR_PYTIGRE, BUILD_WITH_HIP off) -- all changed sources compile; it was not linked or run, as no NVIDIA GPU was available. Other AMD architectures are buildable via HIP_ARCH but were not exercised here.