data_processor only applied in forward pass for ParallelProj Projector

[samdporter]

Observed behaviour:

• ParallelProj: data_processor applied in forward projection only.
• RayTracing: data_processor applied in both forward and backprojection.
• Matrix (e.g. SPECTUB): data_processor applied in both forward and backprojection.
• NiftyPet: I haven't tested this as it's not included in my build

This leads to inconsistent forward/backward operators depending on the projector backend, and the current test (test_data_processor_projectors.cxx) only exercises RayTracing.

A minimal example in SIRF python:

#!/usr/bin/env python3

import os
import numpy as np

from sirf.STIR import (
    AcquisitionData,
    AcquisitionModelUsingParallelproj,
    AcquisitionModelUsingRayTracingMatrix,
    ImageData,
    SeparableGaussianImageFilter,
    examples_data_path,
    MessageRedirector,
)

# silence STIR output
msg = MessageRedirector()


def load_pet_data():
    data_root = examples_data_path("PET")
    acq_path = os.path.join(data_root, "simulated_data.hs")
    img_path = os.path.join(data_root, "test_image_PM_QP_6.hv")

    acq = AcquisitionData(acq_path)
    init = ImageData(img_path)
    return acq, init


def main():
    acq, init = load_pet_data()

    # same random projection data for all tests
    rng = np.random.default_rng(0)
    random_proj = acq.clone()
    random_proj.fill(
        rng.standard_normal(random_proj.as_array().shape).astype(np.float32)
    )

    projector_classes = [
        AcquisitionModelUsingParallelproj,
        AcquisitionModelUsingRayTracingMatrix,
    ]

    for cls in projector_classes:
        print(f"\n=== {cls.__name__} ===")

        # model without PSF
        am_no_psf = cls()
        am_no_psf.set_up(acq, init)

        # model with PSF as image_data_processor
        am_psf = cls()
        psf = SeparableGaussianImageFilter()
        psf.set_fwhms([21, 21, 21])
        am_psf.set_image_data_processor(psf)
        am_psf.set_up(acq, init)

        # check that PSF changes the forward projection
        x = init.clone()
        x.fill(1.0)

        f_no_psf = am_no_psf.forward(x)
        f_psf = am_psf.forward(x)

        f_diff = (f_psf - f_no_psf).norm()
        f_rel = f_diff / max(f_no_psf.norm(), 1e-6)
        print(f"forward: diff={f_diff:.3e}, rel={f_rel:.3e}")

        # check whether PSF changes the backward projection
        b_no_psf = am_no_psf.backward(random_proj)
        b_psf = am_psf.backward(random_proj)

        b_diff = (b_psf - b_no_psf).norm()
        b_rel = b_diff / max(b_no_psf.norm(), 1e-6)
        print(f"backward: diff={b_diff:.3e}, rel={b_rel:.3e}")


if __name__ == "__main__":
    main()

The backward difference is negligible for Parallelproj

[KrisThielemans]

Looks like we need to call this around here. Same probably for NiftyPET.