diff --git a/MATLAB/Algorithms/FDK.m b/MATLAB/Algorithms/FDK.m index a238a54d..bbf44510 100644 --- a/MATLAB/Algorithms/FDK.m +++ b/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? diff --git a/MATLAB/Demos/d26_ShortScanWangWeights.m b/MATLAB/Demos/d26_ShortScanWangWeights.m new file mode 100644 index 00000000..1a74bf96 --- /dev/null +++ b/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)'); diff --git a/Python/demos/d26_ShortScanWangWeights.py b/Python/demos/d26_ShortScanWangWeights.py new file mode 100644 index 00000000..67b24109 --- /dev/null +++ b/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() diff --git a/Python/tests/test_fdk_short_scan_wang.py b/Python/tests/test_fdk_short_scan_wang.py new file mode 100644 index 00000000..2aefa017 --- /dev/null +++ b/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) diff --git a/Python/tigre/algorithms/single_pass_algorithms.py b/Python/tigre/algorithms/single_pass_algorithms.py index c2c3c1ed..2d8a7c2a 100644 --- a/Python/tigre/algorithms/single_pass_algorithms.py +++ b/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")