Operator syntax

The operator API provides a more compositional interface to gratopy’s projection operators. It is currently experimental and primarily focused on parallel-beam Radon transforms.

The legacy gratopy.ProjectionSettings API remains the main documented interface for the full feature set of gratopy. The operator API complements it with a syntax that is often more convenient when working with operator algebra, adjoint operators, and experimental kernels.

Warning

The operator API is experimental. Backward-incompatible changes may be introduced without a full deprecation cycle while the interface and internal abstractions are still settling.

Extent placeholders such as gratopy.utilities.ExtentPlaceholder are supported experimentally for Radon operators when exactly one of the image or detector extents is a placeholder. Passing placeholders for both extents at once remains unsupported and raises NotImplementedError.

Current scope

The current operator API supports in particular:

the gratopy.operator.projection.Radon operator,
adjoints via T,
operator composition and arithmetic,
custom OpenCL kernels via gratopy.operator.opencl.OpenCLKernelSpec.

At the moment, the operator API should be understood as Radon-first. Fanbeam support is not yet implemented at the same level.

Quick example

A Radon transform and its adjoint can be used as follows:

import numpy as np
import pyopencl as cl
import gratopy

ctx = cl.create_some_context(interactive=False)
queue = cl.CommandQueue(ctx)

Nx = 128
img = np.zeros((Nx, Nx), dtype=np.float32)

R = gratopy.operator.Radon(image_domain=Nx, angles=180)

sino = R.apply_to(img, queue=queue)
backprojection = R.T.apply_to(sino)

The same operations can be written with operator syntax:

sino = R * img
backprojection = R.T * sino

When the input is a NumPy array, a queue is still required for the first application so that the data can be transferred to the device.

Detailed geometry example

The operator API also supports explicit geometry helper objects from gratopy.utilities. This is often the clearest way to specify image extent, detector geometry, shifts, and angular sampling explicitly.

import numpy as np
import pyopencl as cl
import gratopy
from gratopy.utilities import Angles, Detectors, ImageDomain

ctx = cl.create_some_context(interactive=False)
queue = cl.CommandQueue(ctx)

img = np.zeros((192, 128), dtype=np.float32)

image_domain = ImageDomain(
    size=(192, 128),
    extent=3.0,
    center=(0.1, -0.2),
)

angles = Angles.uniform_interval(
    start=0.0,
    end=np.pi / 2,
    number=120,
)

detectors = Detectors(
    number=220,
    extent=3.0,
    center=0.15,
    reversed=False,
)

R = gratopy.operator.Radon(
    image_domain=image_domain,
    angles=angles,
    detectors=detectors,
)

sino = R.apply_to(img, queue=queue)
backproj = R.T.apply_to(sino)

The same setup can also be written directly inline when constructing the operator:

R = gratopy.operator.Radon(
    image_domain=ImageDomain(size=(192, 128), extent=3.0, center=(0.1, -0.2)),
    angles=Angles.uniform_interval(0.0, np.pi / 2, 120),
    detectors=Detectors(number=220, extent=3.0, center=0.15),
)

This explicit style is particularly useful when experimenting with geometry in Python code, because image domain, angles, and detector settings become first-class objects that can be reused and modified independently.

Extent placeholders

For Radon operators, one physical extent can be inferred from the other by using gratopy.utilities.ExtentPlaceholder. This is useful when one wants either the smallest detector covering a fixed image domain, or the largest image domain covered by a fixed detector.

For example, the detector extent can be inferred from a fixed image extent:

from gratopy.utilities import Detectors, ExtentPlaceholder, ImageDomain

R = gratopy.operator.Radon(
    image_domain=ImageDomain(size=128, extent=2.0),
    angles=180,
    detectors=Detectors(number=200, extent=ExtentPlaceholder.FULL),
)

Conversely, the image extent can be inferred from a fixed detector extent:

R = gratopy.operator.Radon(
    image_domain=ImageDomain(size=128, extent=ExtentPlaceholder.FULL),
    angles=180,
    detectors=Detectors(number=200, extent=2.0),
)

Only one side may use an extent placeholder at a time. Passing placeholders for both the image and detector extents is unsupported and raises NotImplementedError. If the requested placeholder semantics are geometrically impossible for the supplied centers and fixed extent, construction raises ValueError.

Operator algebra

Operators inherit from gratopy.operator.base.Operator, which supports basic arithmetic and composition. For example, one can form a Gram operator

G = R.T * R

and apply it to an image:

gram_img = G.apply_to(img, queue=queue)

This is one of the main motivations for the operator interface: projection operators can be combined with a syntax that mirrors the underlying linear algebra.

Class structure

The operator implementation is intentionally layered in a small number of classes.

As an experimental interface, the operator API may still change in backward-incompatible ways without a full deprecation cycle while the design is settling.

gratopy.operator.base.Operator: Provides generic operator arithmetic such as addition, scalar multiplication, composition, and application.
gratopy.operator.opencl._OpenCLOperator: Internal helper base for OpenCL-backed operators. It implements shared execution plumbing such as queue inference, array coercion, output allocation, program cacheing, and kernel lookup.
gratopy.operator.projection.Radon: Concrete Radon transform operator. It owns its operator state, geometry-specific preparation, and binds these to the OpenCL kernels.

This separation keeps the generic algebra in gratopy.operator.base backend-agnostic while concentrating OpenCL-specific behavior in a gratopy- specific internal layer.

Custom kernels

One goal of the new operator API is to make kernel experimentation easier. Kernel sources can be specified via gratopy.operator.opencl.OpenCLKernelSpec.

A custom kernel spec can be passed directly to the operator:

from gratopy.operator import OpenCLKernelSpec

spec = OpenCLKernelSpec.from_path("scratch/my_radon.cl", base_name="radon")
R = gratopy.operator.Radon(image_domain=128, angles=180, kernel_spec=spec)

This allows experimenting with alternative kernels while keeping the Python-side operator interface unchanged.

Subclassing _OpenCLOperator

For experimental custom operators, the internal class gratopy.operator.opencl._OpenCLOperator provides a default apply_to() implementation. Although _OpenCLOperator is internal, its documented hook methods form the intended customization surface for OpenCL-backed operators.

The default execution pipeline performs, in order:

queue inference,
coercion of array-like inputs to pyopencl.array.Array,
input validation,
output allocation (if needed),
output validation,
kernel lookup,
kernel invocation,
multiplication by the operator scalar.

Subclasses can adapt this behavior mostly via hooks instead of overriding the entire method.

Important hooks are:

gratopy.operator.opencl._OpenCLOperator._default_kernel_spec() for the default kernel bundle,
gratopy.operator.opencl._OpenCLOperator._expected_output_shape() for shape inference,
gratopy.operator.opencl._OpenCLOperator._validate_argument() and gratopy.operator.opencl._OpenCLOperator._validate_output() for validation,
gratopy.operator.opencl._OpenCLOperator._get_kernel() for choosing the compiled kernel,
gratopy.operator.opencl._OpenCLOperator._kernel_arguments() for supplying additional kernel arguments beyond output and input buffers,
gratopy.operator.opencl._OpenCLOperator._global_shape() for customizing the OpenCL launch shape.

In simple cases, a custom operator only needs to provide a kernel spec, adjoint / T behavior, and static input/output shapes.

Limitations and status

The operator API is still evolving. In particular:

the focus is currently on gratopy.operator.projection.Radon,
extent placeholders are currently supported experimentally for Radon operators when exactly one of the image or detector extents is a placeholder,
higher-level solver interfaces are still centered around the legacy API.

For the full and mature feature set of gratopy, the legacy API documented in Getting started and Reference manual remains the main reference.