From 4e9610b6b74b7b9ac807733f14492b25a284d447 Mon Sep 17 00:00:00 2001
From: zez666 <40331996+zezisme@users.noreply.github.com>
Date: Wed, 12 Aug 2026 07:16:25 +0800
Subject: [PATCH] Fix: deg/rad unit mismatch in removeProjections corrupts
 chunked helical FBP
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

### Bug Summary

In helical cone-beam FBP, when the projection + volume data exceed GPU
memory, LEAP triggers `backproject_FBP_multiGPU_splitViews`, which creates
per-chunk parameter copies via `removeProjections()`. Two bugs in this
function caused incorrect reconstruction values:

**Bug 1 (Major): phis_full stored in degrees, consumed as radians**

`removeProjections()` calls `get_angles(phis_full)` to save the full-scan
angle array. `get_angles()` converts from internal radians to degrees:
  phis_full[i] = (phis[i] + 0.5*PI) * 180.0 / PI

Later, `T_phi()` preferentially reads `phis_full` over `phis`:
  if (phis_full != NULL) phis_local = phis_full;  // ← uses degrees!

Since `T_phi()` returns angular step without deg→rad conversion, and
`FBPscalar ∝ T_phi()`, the normalization applied during Hilbert filtering
is scaled by 180/π ≈ 57.3×. The derivative kernel (`parallelRay_derivative`)
uses `phis[]` (always radians) for its local T_phi, so the 1/T_phi and ×T_phi
do NOT cancel — the error propagates directly to reconstructed CT values.

**Bug 2 (Minor): phi_start/phi_end restored to full-scan range**

After `set_angles()`, `phi_start`/`phi_end` correctly reflect the chunk's
angular range, but are immediately restored to full-scan values:
  phi_start = phi_start_save;  // full-scan
  phi_end = phi_end_save;

These are copied to GPU constant memory (`d_phi_start`, `d_phi_end`) and
used in the helical weighted backprojector kernel to bound the redundancy
`sumWeights` computation. With full-scan bounds, `sumWeights` includes
turns not present in the chunk's data, causing incorrect per-view
normalization.

### Repro Conditions

- Helical cone-beam scan (pitch ≠ 0)
- Projection data + volume data exceed `available_GPU_memory` (triggers
  `backproject_FBP_multiGPU_splitViews`)
- Example: 6000 views × 505 rows × 2063 cols = 23.3 GiB, vol 1527³ × 445
  = 3.9 GiB, required ≈ 50.7 GiB > typical 48 GB GPU → triggered

Small datasets that fit entirely on GPU are NOT affected (no chunking,
`phis_full` remains NULL, `T_phi()` uses `phis[]` in radians).

### Verification

FBPscalar can be printed before FBP (via `get_FBPscalar()`) to verify:
  FBPscalar = 1/(2π) × T_φ × pixW × (sod/sdd)² × pixH / (voxW² × voxH)

With the fix, `FBPscalar ≈ 0.03–0.04` for typical micro-CT geometries
(corresponding to T_φ ≈ 0.004 rad). Before the fix, chunked FBP would
report the same FBPscalar value (computed from the pre-chunk params),
but internally use a 57× larger value for each chunk.

### Files Changed

- `src/parameters.cpp`: `removeProjections()` — store `phis_full` in
  radians, manually set `phis` and `phi_start`/`phi_end` without
  degree-expecting `set_angles()`.
---
 src/parameters.cpp | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/src/parameters.cpp b/src/parameters.cpp
index aeaad95..8264c92 100644
--- a/src/parameters.cpp
+++ b/src/parameters.cpp
@@ -2252,7 +2252,8 @@ bool parameters::removeProjections(int firstProj, int lastProj)
 		if (phis_full != NULL)
 			delete[] phis_full;
 		phis_full = new float[numAngles];
-		get_angles(phis_full);
+		for (int i = 0; i < numAngles; i++)
+			phis_full[i] = phis[i];  // radians
 		phis_new = new float[numAngles_new];
 		for (int i = firstProj; i <= lastProj; i++)
 			phis_new[i - firstProj] = phis_full[i];
@@ -2287,11 +2288,11 @@ bool parameters::removeProjections(int firstProj, int lastProj)
 	}
 	if (phis_new != NULL)
 	{
-		float phi_start_save = phi_start;
-		float phi_end_save = phi_end;
-		set_angles(phis_new, numAngles_new);
-		phi_start = phi_start_save;
-		phi_end = phi_end_save;
+		delete[] phis;
+		phis = phis_new;
+		phis_new = NULL;
+		phi_start = min(phis[0], phis[numAngles_new - 1]);
+		phi_end   = max(phis[0], phis[numAngles_new - 1]);
 		angularRange = angularRange_save;
 	}
 	numAngles = numAngles_new;