lanl / scico

Scientific Computational Imaging COde
BSD 3-Clause "New" or "Revised" License
105 stars 17 forks source link

Apparent regression in `linop.Parallel2dProjector` #534

Closed bwohlberg closed 4 months ago

bwohlberg commented 5 months ago

The projection accuracy of linop.Parallel2dProjector seems to have significantly degraded since the notebook ct_multi_cs_tv_admm.ipynb currently in scico-data was built. Compare the performance of the scico projector reported in that notebook with a recent build, noting the roughly 0.75 dB SNR reduction in reconstruction quality.

bwohlberg commented 4 months ago

For the most part the output of the scico projector matches those of the other two projectors very closely, but there seems to be a significant difference at a boundary, as illustrated by this script:

from scico import functional, linop, loss, metric, plot
from scico.linop.xray import Parallel2dProjector, XRayTransform, astra, svmbir
from scico.optimize.admm import ADMM, LinearSubproblemSolver
from scico.util import device_info

N = 512
np.random.seed(1234)
x_gt = snp.array(discrete_phantom(Foam(size_range=[0.075, 0.0025], gap=1e-3, porosity=1), size=N))

n_projection = 45
angles = np.linspace(0, np.pi, n_projection)
projectors = {
    "astra": astra.XRayTransform(x_gt.shape, 1, N, angles - np.pi / 2.0),  # astra
    "svmbir": svmbir.XRayTransform(x_gt.shape, 2 * np.pi - angles, N),  # svmbir
    "scico": XRayTransform(Parallel2dProjector((N, N), angles, det_count=N)),  # scico
}

y = {}
algs = ("astra", "svmbir", "scico")
for p in algs:
    A = projectors[p]
    y[p] = A @ x_gt

y_gt = x_gt.sum(axis=1)
fig, ax = plot.plot(np.vstack([y_gt] + [y[p][0] for p in algs]).T,
                    lgnd=("row sum",) + algs)
ax.set_xlim((-0.5, 2))
ax.set_ylim((-5, 60))

input()