UCL / UCL/STIR

data_processor only applied in forward pass for ParallelProj Projector

Open
#1,669 1 comment 0 reactions 1 assignee View on GitHub

@KrisThielemans is already working on this.

Since Nov 26, 2025.

bug
Dominant language
C++
Stars
160
Forks
113
Avg merge
12d 15h
Merged PRs (30d)
1

Description

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

Contributor guide

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.