Refactor FBP and FDK implementations to utilize new weighting and filtering functions

sypsyp97 · sypsyp97 · commit 34ed5bff96c2 · 2026-02-14T09:11:31.000+01:00
- Updated `fbp_fan.py` to replace manual ramp filtering with `ramp_filter_1d` and added cosine weighting using `fan_cosine_weights`.
- Integrated Parker weighting option in fan-beam reconstruction.
- Refactored `fbp_parallel.py` to use `ramp_filter_1d` and added angular integration weights.
- Enhanced `fdk_cone.py` to implement cosine weighting with `cone_cosine_weights` and replaced 3D ramp filtering with `ramp_filter_1d`.
- Introduced new tests for CUDA functionality and weight calculations in `test_cuda_smoke.py` and `test_weights.py`.
- Added `pytest.ini` for test configuration.
diff --git a/diffct/__init__.py b/diffct/__init__.py
@@ -1,2 +1,17 @@
 # diffct/__init__.py
-from .differentiable import ParallelProjectorFunction, ParallelBackprojectorFunction, FanProjectorFunction, FanBackprojectorFunction, ConeProjectorFunction, ConeBackprojectorFunction
+from .differentiable import (
+    ParallelProjectorFunction,
+    ParallelBackprojectorFunction,
+    FanProjectorFunction,
+    FanBackprojectorFunction,
+    ConeProjectorFunction,
+    ConeBackprojectorFunction,
+    detector_coordinates_1d,
+    angular_integration_weights,
+    fan_cosine_weights,
+    cone_cosine_weights,
+    parker_weights,
+    ramp_filter_1d,
+    fan_weighted_backproject,
+    cone_weighted_backproject,
+)
diff --git a/diffct/differentiable.py b/diffct/differentiable.py
diff --git a/examples/fbp_fan.py b/examples/fbp_fan.py
@@ -3,7 +3,14 @@
 import torch
 import matplotlib.pyplot as plt
 import torch.nn.functional as F
-from diffct.differentiable import FanProjectorFunction, FanBackprojectorFunction
+from diffct.differentiable import (
+    FanProjectorFunction,
+    angular_integration_weights,
+    fan_cosine_weights,
+    fan_weighted_backproject,
+    parker_weights,
+    ramp_filter_1d,
+)
 
 
 def shepp_logan_2d(Nx, Ny):
@@ -34,19 +41,6 @@ def shepp_logan_2d(Nx, Ny):
     phantom = np.clip(phantom, 0.0, 1.0)
     return phantom
 
-def ramp_filter(sinogram_tensor):
-    device = sinogram_tensor.device
-    num_views, num_det = sinogram_tensor.shape
-    freqs = torch.fft.fftfreq(num_det, device=device)
-    omega = 2.0 * torch.pi * freqs
-    ramp = torch.abs(omega)
-    ramp_2d = ramp.reshape(1, num_det)
-    sino_fft = torch.fft.fft(sinogram_tensor, dim=1)
-    filtered_fft = sino_fft * ramp_2d
-    filtered = torch.real(torch.fft.ifft(filtered_fft, dim=1))
-    
-    return filtered
-
 def main():
     Nx, Ny = 256, 256
     phantom = shepp_logan_2d(Nx, Ny)
@@ -55,45 +49,61 @@ def main():
 
     num_detectors = 600
     detector_spacing = 1.0
+    detector_offset = 0.0
     voxel_spacing = 1.0
     sdd = 800.0
     sid = 500.0
+    apply_parker = False
 
     device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
-    image_torch = torch.tensor(phantom, device=device, dtype=torch.float32, requires_grad=True)
+    image_torch = torch.tensor(phantom, device=device, dtype=torch.float32)
     angles_torch = torch.tensor(angles_np, device=device, dtype=torch.float32)
 
     sinogram = FanProjectorFunction.apply(image_torch, angles_torch, num_detectors,
                                           detector_spacing, sdd, sid, voxel_spacing)
 
     # --- FBP weighting and filtering ---
-    # For fan-beam FBP, projections must be weighted before filtering.
-    # Weight = cos(gamma), where gamma is the fan angle for each detector.
-    u = (torch.arange(num_detectors, dtype=image_torch.dtype, device=device) - (num_detectors - 1) / 2) * detector_spacing
-    gamma = torch.atan(u / sdd)
-    weights = torch.cos(gamma).unsqueeze(0)  # Shape (1, num_detectors) for broadcasting
-
-    # Apply weights before filtering
+    # 1) Optional Parker redundancy weighting for short-scan trajectories
+    if apply_parker:
+        parker = parker_weights(angles_torch, num_detectors, detector_spacing, sdd, detector_offset)
+        sinogram = sinogram * parker
+
+    # 2) Fan-beam cosine pre-weighting
+    weights = fan_cosine_weights(
+        num_detectors,
+        detector_spacing,
+        sdd,
+        detector_offset=detector_offset,
+        device=device,
+        dtype=image_torch.dtype,
+    ).unsqueeze(0)
     sino_weighted = sinogram * weights
-    sinogram_filt = ramp_filter(sino_weighted)
-
-    reconstruction = F.relu(FanBackprojectorFunction.apply(sinogram_filt, angles_torch,
-                                                    detector_spacing, Ny, Nx,
-                                                    sdd, sid, voxel_spacing)) # ReLU to ensure non-negativity
-    
-    # --- FBP normalization ---
-    # The backprojection is a sum over all angles. To approximate the integral,
-    # we need to multiply by the angular step d_beta.
-    # The fan-beam FBP formula also includes a factor of 1/2 when integrating over [0, 2*pi].
-    # d_beta = 2 * pi / num_angles
-    # Normalization factor = (1/2) * d_beta = pi / num_angles
-    reconstruction = reconstruction * (math.pi / num_angles)
+
+    # 3) Ramp filter along detector axis
+    sinogram_filt = ramp_filter_1d(sino_weighted, dim=1)
+
+    # 4) Angle-integration weights
+    d_beta = angular_integration_weights(angles_torch, redundant_full_scan=(not apply_parker)).view(-1, 1)
+    sinogram_filt = sinogram_filt * d_beta
+
+    # 5) Weighted fan-beam backprojection
+    reconstruction = F.relu(
+        fan_weighted_backproject(
+            sinogram_filt,
+            angles_torch,
+            detector_spacing,
+            Ny,
+            Nx,
+            sdd,
+            sid,
+            voxel_spacing=voxel_spacing,
+            detector_offset=detector_offset,
+        )
+    )
 
     loss = torch.mean((reconstruction - image_torch)**2)
-    loss.backward()
 
     print("Loss:", loss.item())
-    print("Center pixel gradient:", image_torch.grad[Ny//2, Nx//2].item())
 
     sinogram_cpu = sinogram.detach().cpu().numpy()
     reco_cpu = reconstruction.detach().cpu().numpy()
@@ -119,4 +129,4 @@ def main():
     print("Reco range:", reco_cpu.min(), reco_cpu.max())
 
 if __name__ == "__main__":
-    main()
+    main()
diff --git a/examples/fbp_parallel.py b/examples/fbp_parallel.py
@@ -2,7 +2,12 @@
 import torch
 import matplotlib.pyplot as plt
 import torch.nn.functional as F
-from diffct.differentiable import ParallelProjectorFunction, ParallelBackprojectorFunction
+from diffct.differentiable import (
+    ParallelProjectorFunction,
+    ParallelBackprojectorFunction,
+    angular_integration_weights,
+    ramp_filter_1d,
+)
 
 
 def shepp_logan_2d(Nx, Ny):
@@ -33,19 +38,6 @@ def shepp_logan_2d(Nx, Ny):
     phantom = np.clip(phantom, 0.0, 1.0)
     return phantom
 
-def ramp_filter(sinogram_tensor):
-    device = sinogram_tensor.device
-    num_views, num_det = sinogram_tensor.shape
-    freqs = torch.fft.fftfreq(num_det, device=device)
-    omega = 2.0 * torch.pi * freqs
-    ramp = torch.abs(omega)
-    ramp_2d = ramp.reshape(1, num_det)
-    sino_fft = torch.fft.fft(sinogram_tensor, dim=1)
-    filtered_fft = sino_fft * ramp_2d
-    filtered = torch.real(torch.fft.ifft(filtered_fft, dim=1))
-    
-    return filtered
-
 def main():
     Nx, Ny = 256, 256
     phantom = shepp_logan_2d(Nx, Ny)
@@ -63,18 +55,12 @@ def main():
     sinogram = ParallelProjectorFunction.apply(image_torch, angles_torch,
                                                num_detectors, detector_spacing, voxel_spacing)
     
-    sinogram_filt = ramp_filter(sinogram)
+    sinogram_filt = ramp_filter_1d(sinogram, dim=1)
+    d_theta = angular_integration_weights(angles_torch, redundant_full_scan=True).view(-1, 1)
+    sinogram_filt = sinogram_filt * d_theta
 
     reconstruction = F.relu(ParallelBackprojectorFunction.apply(sinogram_filt, angles_torch,
                                                          detector_spacing, Ny, Nx, voxel_spacing)) # ReLU to ensure non-negativity
-    
-    # --- FBP normalization ---
-    # The backprojection is a sum over all angles. To approximate the integral,
-    # we need to multiply by the angular step d_theta.
-    # The FBP formula also includes a factor of 1/2 when integrating over [0, 2*pi].
-    # d_theta = 2 * pi / num_angles
-    # Normalization factor = (1/2) * d_theta = pi / num_angles
-    reconstruction = reconstruction * (np.pi / num_angles)
 
     loss = torch.mean((reconstruction - image_torch)**2)
     loss.backward()
@@ -107,4 +93,4 @@ def main():
     print("Reco range:", reco_cpu.min(), reco_cpu.max())
 
 if __name__ == "__main__":
-    main()
+    main()
diff --git a/examples/fdk_cone.py b/examples/fdk_cone.py
@@ -3,7 +3,13 @@
 import torch
 import matplotlib.pyplot as plt
 import torch.nn.functional as F
-from diffct.differentiable import ConeProjectorFunction, ConeBackprojectorFunction
+from diffct.differentiable import (
+    ConeProjectorFunction,
+    angular_integration_weights,
+    cone_cosine_weights,
+    cone_weighted_backproject,
+    ramp_filter_1d,
+)
 
 
 def shepp_logan_3d(shape):
@@ -59,19 +65,6 @@ def shepp_logan_3d(shape):
     shepp_logan = np.clip(shepp_logan, 0, 1)
     return shepp_logan
 
-def ramp_filter_3d(sinogram_tensor):
-    device = sinogram_tensor.device
-    num_views, num_det_u, num_det_v = sinogram_tensor.shape
-    freqs = torch.fft.fftfreq(num_det_u, device=device)
-    omega = 2.0 * torch.pi * freqs
-    ramp = torch.abs(omega)
-    ramp_3d = ramp.reshape(1, num_det_u, 1)
-    sino_fft = torch.fft.fft(sinogram_tensor, dim=1)
-    filtered_fft = sino_fft * ramp_3d
-    filtered = torch.real(torch.fft.ifft(filtered_fft, dim=1))
-    
-    return filtered
-
 def main():
     Nx, Ny, Nz = 128, 128, 128
     phantom_cpu = shepp_logan_3d((Nz, Ny, Nx))
@@ -81,53 +74,65 @@ def main():
 
     det_u, det_v = 256, 256
     du, dv = 1.0, 1.0
+    detector_offset_u = 0.0
+    detector_offset_v = 0.0
     sdd = 900.0
     sid = 600.0
 
     voxel_spacing = 1.0
 
     device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
-    phantom_torch = torch.tensor(phantom_cpu, device=device, dtype=torch.float32, requires_grad=True).contiguous()
+    phantom_torch = torch.tensor(phantom_cpu, device=device, dtype=torch.float32).contiguous()
     angles_torch = torch.tensor(angles_np, device=device, dtype=torch.float32)
 
     sinogram = ConeProjectorFunction.apply(phantom_torch, angles_torch,
                                            det_u, det_v, du, dv,
                                            sdd, sid, voxel_spacing)
 
     # --- FDK weighting and filtering ---
-    # For FDK, projections must be weighted before filtering.
-    # Weight = D / sqrt(D^2 + u^2 + v^2), where D is source_distance
-    # and (u,v) are detector coordinates.
-    u_coords = (torch.arange(det_u, dtype=phantom_torch.dtype, device=device) - (det_u - 1) / 2) * du
-    v_coords = (torch.arange(det_v, dtype=phantom_torch.dtype, device=device) - (det_v - 1) / 2) * dv
-
-    # Reshape for broadcasting over sinogram of shape (views, u, v)
-    u_coords = u_coords.view(1, det_u, 1)
-    v_coords = v_coords.view(1, 1, det_v)
-    
-    weights = sdd / torch.sqrt(sdd**2 + u_coords**2 + v_coords**2)
-    
-    # Apply weights and then filter
+    # 1) FDK cosine pre-weighting
+    weights = cone_cosine_weights(
+        det_u,
+        det_v,
+        du,
+        dv,
+        sdd,
+        detector_offset_u=detector_offset_u,
+        detector_offset_v=detector_offset_v,
+        device=device,
+        dtype=phantom_torch.dtype,
+    ).unsqueeze(0)
     sino_weighted = sinogram * weights
-    sinogram_filt = ramp_filter_3d(sino_weighted).contiguous()
-
-    reconstruction = F.relu(ConeBackprojectorFunction.apply(sinogram_filt, angles_torch, Nz, Ny, Nx,
-                                                    du, dv, sdd, sid, voxel_spacing)) # ReLU to ensure non-negativity
-    
-    # --- FDK normalization ---
-    # The backprojection is a sum over all angles. To approximate the integral,
-    # we need to multiply by the angular step d_beta.
-    # The FDK formula also includes a factor of 1/2 when integrating over [0, 2*pi].
-    # d_beta = 2 * pi / num_views
-    # Normalization factor = (1/2) * d_beta = pi / num_views
-    reconstruction = reconstruction * (math.pi / num_views)
+
+    # 2) Ramp filter along detector-u rows
+    sinogram_filt = ramp_filter_1d(sino_weighted, dim=1).contiguous()
+
+    # 3) Angle-integration weights
+    d_beta = angular_integration_weights(angles_torch, redundant_full_scan=True).view(-1, 1, 1)
+    sinogram_filt = sinogram_filt * d_beta
+
+    # 4) Weighted cone-beam backprojection
+    reconstruction = F.relu(
+        cone_weighted_backproject(
+            sinogram_filt,
+            angles_torch,
+            Nz,
+            Ny,
+            Nx,
+            du,
+            dv,
+            sdd,
+            sid,
+            voxel_spacing=voxel_spacing,
+            detector_offset_u=detector_offset_u,
+            detector_offset_v=detector_offset_v,
+        )
+    )
 
     loss = torch.mean((reconstruction - phantom_torch)**2)
-    loss.backward()
 
     print("Cone Beam Example with user-defined geometry:")
     print("Loss:", loss.item())
-    print("Volume center voxel gradient:", phantom_torch.grad[Nz//2, Ny//2, Nx//2].item())
     print("Reconstruction shape:", reconstruction.shape)
 
     reconstruction_cpu = reconstruction.detach().cpu().numpy()
@@ -155,4 +160,4 @@ def main():
     print("Reco data range:", reconstruction_cpu.min(), reconstruction_cpu.max())
 
 if __name__ == "__main__":
-    main()
+    main()
diff --git a/pyproject.toml b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "hatchling.build"
 
 [project]
 name = "diffct"
-version = "1.2.7"
+version = "1.2.8"
 description = "A CUDA-based library for computed tomography (CT) projection and reconstruction with differentiable operators"
 readme = "README.md"
 authors = [
diff --git a/pytest.ini b/pytest.ini
@@ -0,0 +1,3 @@
+[pytest]
+markers =
+    cuda: CUDA-dependent smoke tests
diff --git a/tests/test_cuda_smoke.py b/tests/test_cuda_smoke.py
diff --git a/tests/test_weights.py b/tests/test_weights.py

Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,3 @@`
	`1`	`+[pytest]`
	`2`	`+markers =`
	`3`	`+ cuda: CUDA-dependent smoke tests`