FDK: do not apply Wang detector-offset weights on short scans (Python…

MATLAB/Algorithms/FDK.m

@@ -43,6 +43,15 @@
 geo=checkGeo(geo,angles);
 geo.filter=filter;
 
+if dowang && (max(angles(1,:))-min(angles(1,:)) < 2*pi-1.5*max(abs(diff(angles(1,:)))))
+    % Wang's displaced-detector weights assume a FULL circle: they ramp one
+    % side of the detector down and rely on the opposing (beta+pi) views to
+    % bring the coverage back to uniform. On a short scan those views do not
+    % exist and the ramp survives in the image as one-sided shading. Same
+    % short-scan test as the Parker default in parse_inputs.
+    warning('FDK: short scan, Wang detector-offset weights not applied (they assume a full 360-degree circle)');
+    dowang=false;
+end
 if dowang
     % Zero-padding to avoid FFT-induced aliasing %TODO: should't this be
     % for all cases, not just wang?

MATLAB/Demos/d26_ShortScanWangWeights.m

@@ -0,0 +1,77 @@
+%% Demo 26: Wang detector-offset weights on a short scan
+%
+% Wang weighting compensates a displaced detector on a FULL circular scan: it
+% ramps one side of the detector down and relies on the opposing (beta + pi)
+% views to bring the coverage back to uniform. On a short scan those opposing
+% views do not exist, so the ramp survives into the reconstruction as a
+% one-sided shading (a "teardrop"). Any non-zero offDetector switches the
+% weights on - a sub-pixel calibrated offset is enough - so a short scan with a
+% calibrated geometry was silently shaded.
+%
+% FDK now skips the Wang weights when the angles do not cover a full circle
+% (the same test it already uses to decide Parker weighting). This demo
+% reconstructs a uniform cylinder from a 200-degree scan three ways and prints
+% the left/right interior ratio:
+%
+%   centred detector                    -> balanced (reference)
+%   0.1 mm offset, Wang applied by hand -> the teardrop (what happened before)
+%   0.1 mm offset, FDK default          -> Wang skipped, balanced again
+%
+% Wang weights are still applied on full scans, where they belong.
+%--------------------------------------------------------------------------
+% This file is part of the TIGRE Toolbox
+%
+% Copyright (c) 2015, University of Bath and
+%                     CERN-European Organization for Nuclear Research
+%                     All rights reserved.
+%
+% License:            Open Source under BSD.
+%                     See the full license at
+%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
+%
+% Contact:            tigre.toolbox@gmail.com
+% Codes:              https://github.com/CERN/TIGRE/
+%--------------------------------------------------------------------------
+clear; close all;
+
+%% Geometry: default, 64^3 image
+geo = defaultGeometry('nVoxel', [64; 64; 64]);
+geo.sVoxel = [64; 64; 64];
+geo.dVoxel = geo.sVoxel ./ geo.nVoxel;
+
+%% Uniform cylinder, mu = 0.02
+[xx, yy, zz] = meshgrid(1:64, 1:64, 1:64);
+mu = 0.02;
+phantom = single(((xx - 32.5).^2 + (yy - 32.5).^2 <= 20^2) * mu);
+core = (xx - 32.5).^2 + (yy - 32.5).^2 <= 14^2;
+
+%% Short scan: 200 degrees
+angles = linspace(0, deg2rad(200), 200);
+proj = Ax(phantom, geo, angles);
+
+lr = @(vol) [mean(vol(core & xx <= 32 & zz >= 17 & zz <= 48)), ...
+             mean(vol(core & xx >  32 & zz >= 17 & zz <= 48))];
+
+%% 1. centred detector (reference)
+vol_ref = FDK(proj, geo, angles);
+r = lr(vol_ref);
+fprintf('centred detector           : left %.5f right %.5f ratio %.3f\n', r(1), r(2), r(1) / r(2));
+
+%% 2. 0.1 mm offset, Wang weights applied regardless of the arc (the old behaviour)
+geo_off = geo;
+geo_off.offDetector = [0.1; 0];          % u offset, a fraction of a pixel
+proj_w = proj .* redundancy_weighting(geo_off);
+vol_wang = FDK(proj_w, geo_off, angles, 'wang', false);
+r = lr(vol_wang);
+fprintf('offset, Wang applied (old) : left %.5f right %.5f ratio %.3f\n', r(1), r(2), r(1) / r(2));
+
+%% 3. 0.1 mm offset, FDK default: Wang skipped on the short scan (with a warning)
+vol_gated = FDK(proj, geo_off, angles);
+r = lr(vol_gated);
+fprintf('offset, FDK default (new)  : left %.5f right %.5f ratio %.3f\n', r(1), r(2), r(1) / r(2));
+
+%% Show the central axial slice of each
+figure('Name', '200-degree scan: Wang detector-offset weights need a full circle');
+subplot(1, 3, 1); imshow(squeeze(vol_ref(:, :, 32)),   [0 1.5 * mu]); title('centred (reference)');
+subplot(1, 3, 2); imshow(squeeze(vol_wang(:, :, 32)),  [0 1.5 * mu]); title('offset, Wang applied (before)');
+subplot(1, 3, 3); imshow(squeeze(vol_gated(:, :, 32)), [0 1.5 * mu]); title('offset, Wang skipped (after)');

Python/demos/d26_ShortScanWangWeights.py

@@ -0,0 +1,77 @@
+"""Demo 26: Wang detector-offset weights on a short scan.
+
+Wang weighting compensates a displaced detector on a FULL circular scan: it
+ramps one side of the detector down and relies on the opposing (beta + pi)
+views to bring the coverage back to uniform. On a short scan those opposing
+views do not exist, so the ramp survives into the reconstruction as a
+one-sided shading (a "teardrop"). Any non-zero offDetector switches the
+weights on - a sub-pixel calibrated offset is enough - so a short scan with a
+calibrated geometry was silently shaded.
+
+FDK now skips the Wang weights when the angles do not cover a full circle
+(the same test MATLAB's FDK uses to decide Parker weighting). This demo
+reconstructs a uniform cylinder from a 200-degree scan three ways and reports
+the left/right interior ratio:
+
+  centred detector            -> balanced (reference)
+  0.1 mm offset, Wang forced  -> the teardrop (what happened before)
+  0.1 mm offset, default      -> Wang skipped, balanced again
+
+Wang weights are still applied on full scans, where they belong.
+"""
+import numpy as np
+import tigre
+import tigre.algorithms as algs
+import matplotlib.pyplot as plt
+
+geo = tigre.geometry_default(high_resolution=False)
+geo.nVoxel = np.array([64, 64, 64])
+geo.sVoxel = np.array([64.0, 64.0, 64.0])
+geo.dVoxel = geo.sVoxel / geo.nVoxel
+
+zz, yy, xx = np.mgrid[:64, :64, :64]
+mu = 0.02
+phantom = (((yy - 31.5) ** 2 + (xx - 31.5) ** 2 <= 20 ** 2) * mu).astype(np.float32)
+core = (yy - 31.5) ** 2 + (xx - 31.5) ** 2 <= 14 ** 2
+
+angles = np.linspace(0, np.deg2rad(200), 200)      # short scan
+proj = tigre.Ax(phantom, geo, angles)
+
+
+def lr(vol):
+    left = float(vol[16:48, :, :32][core[16:48, :, :32]].mean())
+    right = float(vol[16:48, :, 32:][core[16:48, :, 32:]].mean())
+    return left, right
+
+
+vol_ref = algs.fdk(proj, geo, angles)                 # centred detector
+
+geo_off = tigre.geometry_default(high_resolution=False)
+geo_off.nVoxel, geo_off.sVoxel, geo_off.dVoxel = geo.nVoxel, geo.sVoxel, geo.dVoxel
+geo_off.offDetector = np.array([0.0, 0.1])            # 0.1 mm: a fraction of a pixel
+
+# The old behaviour: Wang weights applied because offDetector != 0. Reproduce
+# it by calling the Wang path directly on a full-circle-shaped call.
+from tigre.algorithms import single_pass_algorithms as spa
+_orig = spa.is_short_scan
+spa.is_short_scan = lambda a: False                   # force the pre-fix path
+vol_wang = algs.fdk(proj, geo_off, angles)
+spa.is_short_scan = _orig
+vol_gated = algs.fdk(proj, geo_off, angles)           # new default: skipped
+
+fig, axes = plt.subplots(1, 3, figsize=(13, 4.4))
+for ax, vol, title in zip(
+        axes, (vol_ref, vol_wang, vol_gated),
+        ("centred detector (reference)",
+         "0.1 mm offset, Wang weights applied (before)",
+         "0.1 mm offset, Wang skipped on the short scan (after)")):
+    l, r = lr(vol)
+    ax.imshow(vol[32], cmap="gray", vmin=0, vmax=1.5 * mu)
+    ax.set_title(f"{title}\nleft/right interior ratio {l / r:.2f}", fontsize=9)
+    ax.axis("off")
+    print(f"{title}: left {l:.5f} right {r:.5f} ratio {l / r:.3f}")
+fig.suptitle("200-degree scan of a uniform cylinder - Wang detector-offset weights need a full circle",
+             fontsize=10)
+plt.tight_layout()
+plt.savefig("d26_short_scan_wang.png", dpi=120)
+plt.show()

Python/tests/test_fdk_short_scan_wang.py

@@ -0,0 +1,82 @@
+"""FDK must not apply Wang detector-offset weights on a short scan.
+
+Wang weighting (displaced-detector redundancy weights) assumes a full circle:
+it ramps one side of the detector down and relies on the opposing (beta + pi)
+views to restore uniform coverage. On a short scan those views do not exist,
+so the ramp survives as a one-sided "teardrop" shading. Any non-zero
+offDetector - a sub-pixel calibration value is enough - switches the weights
+on, so short scans with a calibrated detector offset were silently wrong.
+
+The classification tests are pure NumPy. The reconstruction test needs a GPU
+and is skipped without one.
+"""
+import numpy as np
+import pytest
+
+from tigre.algorithms.single_pass_algorithms import is_short_scan
+
+
+class TestIsShortScan:
+    def test_full_circle_endpoint_false_is_not_short(self):
+        assert not is_short_scan(np.linspace(0, 2 * np.pi, 360, endpoint=False))
+
+    def test_full_circle_endpoint_true_is_not_short(self):
+        assert not is_short_scan(np.linspace(0, 2 * np.pi, 361))
+
+    def test_200_degrees_is_short(self):
+        assert is_short_scan(np.linspace(0, np.deg2rad(200), 200))
+
+    def test_half_circle_is_short(self):
+        assert is_short_scan(np.linspace(0, np.pi, 180, endpoint=False))
+
+    def test_euler_triplets_use_first_column(self):
+        a = np.linspace(0, 2 * np.pi, 360, endpoint=False)
+        full = np.stack([a, np.zeros_like(a), np.zeros_like(a)], axis=1)
+        short = np.stack([a[:200], np.zeros(200), np.zeros(200)], axis=1)
+        assert not is_short_scan(full)
+        assert is_short_scan(short)
+
+    def test_descending_angles(self):
+        assert not is_short_scan(np.linspace(2 * np.pi, 0, 360, endpoint=False))
+        assert is_short_scan(np.linspace(np.deg2rad(200), 0, 200))
+
+    def test_single_angle_is_short(self):
+        assert is_short_scan(np.array([0.3]))
+
+
+def _gpu_available():
+    try:
+        from tigre.utilities.gpu import getGpuNames
+        return len(getGpuNames()) > 0
+    except Exception:
+        return False
+
+
+@pytest.mark.skipif(not _gpu_available(), reason="needs a CUDA GPU")
+def test_short_scan_with_subpixel_offset_has_no_teardrop():
+    """A 200-degree scan of a uniform cylinder with a 0.1 mm (sub-pixel)
+    detector offset: the interior must reconstruct left/right balanced, as
+    the centred short scan does. Before the gate the Wang ramp made one half
+    much brighter than the other."""
+    import tigre
+    import tigre.algorithms as algs
+
+    geo = tigre.geometry_default(high_resolution=False)
+    geo.nVoxel = np.array([64, 64, 64])
+    geo.sVoxel = np.array([64.0, 64.0, 64.0])
+    geo.dVoxel = geo.sVoxel / geo.nVoxel
+    zz, yy, xx = np.mgrid[:64, :64, :64]
+    phantom = (((yy - 31.5) ** 2 + (xx - 31.5) ** 2 <= 20 ** 2) * 0.02).astype(np.float32)
+    core = (yy - 31.5) ** 2 + (xx - 31.5) ** 2 <= 14 ** 2
+    angles = np.linspace(0, np.deg2rad(200), 200)
+    proj = tigre.Ax(phantom, geo, angles)
+
+    def lr_ratio(vol):
+        left = float(vol[16:48, :, :32][core[16:48, :, :32]].mean())
+        right = float(vol[16:48, :, 32:][core[16:48, :, 32:]].mean())
+        return left / right
+
+    ref = lr_ratio(algs.fdk(proj, geo, angles))                  # centred
+    geo.offDetector = np.array([0.0, 0.1])                       # sub-pixel
+    gated = lr_ratio(algs.fdk(proj, geo, angles))
+    assert gated == pytest.approx(ref, rel=0.05)

Python/tigre/algorithms/single_pass_algorithms.py

@@ -9,6 +9,23 @@
 from tigre.utilities.filtering import filtering
 
 
+def is_short_scan(angles):
+    """True when the projection angles do not cover a full circle.
+
+    Same test as MATLAB's FDK.m uses to decide Parker weighting: the arc
+    ``max - min`` of the in-plane angle (first Euler column for (n, 3) input)
+    is shorter than ``2*pi`` by more than 1.5 angular steps, so a full circle
+    sampled with ``endpoint=False`` (arc = 2*pi - one step) is NOT short.
+    """
+    a = np.asarray(angles, dtype=float)
+    if a.ndim > 1:
+        a = a[:, 0]
+    if a.size < 2:
+        return True
+    step = float(np.max(np.abs(np.diff(a))))
+    return float(np.max(a) - np.min(a)) < 2 * np.pi - 1.5 * step
+
+
 def FDK(proj, geo, angles, **kwargs):
     """
     solves CT image reconstruction.
@@ -165,6 +182,18 @@ def preweighting2(proj, geo, theta):
     if not np.any(geo.offDetector):
         dowang = False
 
+    # Wang's displaced-detector weights assume a FULL circle: they ramp one
+    # side of the detector down and rely on the opposing (beta + pi) views to
+    # bring the coverage back to uniform. On a short scan those views do not
+    # exist, so the ramp survives into the image as a one-sided "teardrop"
+    # shading - and ANY non-zero offDetector (a sub-pixel calibration value is
+    # enough) switches the weights on. Skip them on short scans.
+    if dowang and is_short_scan(angles):
+        if verbose:
+            print("FDK: short scan - Wang detector-offset weights not applied "
+                  "(they assume a full 360-degree circle)")
+        dowang = False
+
     if dowang:
         if verbose:
             print("FDK: applying detector offset weights")