Faster Scientific Image Analysis with NVIDIA cuPhoton

Observatories and telescopes, lasers and X-ray light sources, and other high-throughput instruments generate image data faster than CPU-bound pipelines can process it to support timely decisions. The computational bottleneck is rarely one slow kernel. It’s the whole path from data in the sensor to a decision: read the raw data, match point-spread functions or detector responses, subtract or reduce, fit, classify, and alert.

Modern facilities routinely accumulate petabytes of multidimensional data across a single survey or experiment campaign. With CPU-first implementations, it can take hours to months—sometimes a year or more—to obtain scientifically useful results from data the instrument produced in seconds. These delays make it harder for scientists and engineers to extract insights and advance their research. What is missing is a GPU-native path that covers every stage from the sensor to the decision, not a faster version of any one step.

NVIDIA cuPhoton is an open-source NVIDIA CUDA-X toolkit of GPU-accelerated building blocks for that path, spanning use cases from spectral and optical astronomy to time-domain laser/X-ray analysis. By keeping image data on the GPU from sensor read through classification, cuPhoton collapses the wait to seconds, resulting in a tight loop that researchers can run interactively, shortening time to insight and letting facilities keep pace with the data their instruments are producing.

Accelerating scientific insights from months to minutes with cuPhoton

Across astronomy, X-ray science, and other high-throughput domains, NVIDIA cuPhoton provides tools for loading, processing, and analyzing massive multidimensional datasets. The following example combines cuPhoton modules into an end-to-end workflow.

The NSF-DOE Vera C. Rubin Observatory is a prominent example of modern scientific facilities with high demand for fast extraction of scientific insights. Every 39 seconds after dark, the Vera C. Rubin Observatory’s LSSTCam records a new 3.2-gigapixel exposure of the southern sky. The Prompt Processing pipeline compares each frame with a reference template of the same sky region. Within about 60 to 120 seconds, it must classify roughly 10,000 detections as astrophysical transients or artifacts such as cosmic rays, satellite trails, and processing errors. Across an entire night, that amounts to up to 20 terabytes of images and 10 million candidate objects.

NVIDIA cuPhoton scales across multi-GPU, multi-node NVIDIA Grace Blackwell and NVIDIA Vera Rubin systems. On representative workloads involving hundreds of terabytes of data using multiple GPUs, cuPhoton has accelerated image loading and reading by up to 14,900x and signal processing by up to 14,550x, turning months of computational time into minutes. Data analytics that previously took nine months have been demonstrated in four hours on GPU-accelerated Python. For smaller datasets on the order of kilobytes, the entire cuPhoton pipeline executes on millisecond or even microsecond timescales. Figure 3 compares speedups for individual cuPhoton operations against an x86 CPU baseline; these results are not an end-to-end pipeline speedup.

Walkthrough of an example scientific workflow

This post walks through a realistic scientific workflow, using a physically motivated synthetic pair of astronomical images. The following examples are kept small so they can run on NVIDIA DGX Spark or a single-GPU workstation. The same workflow scales across multi-GPU workstations and clusters; that launch is shown later in this post. Figure 4 shows the workflow stages covered in the following sections.

1. Load the current exposure and reference data into GPU memory.
2. Align: Place both visits on a shared sky grid so the same star falls on the same pixel.
3. PSF-match and subtract: Match the point-spread functions (PSFs)—the patterns that describe how light is distributed across pixels—in the current exposure and an earlier reference of the same region. Then subtract the reference. A moving object leaves a paired pattern of positive and negative pixels called a dipole.
4. Fit a dual-PSF model to the dipole to estimate the moving object’s positions.
5. Visualize the results: display the PSF-matching kernel, subtraction residual, and recovered dipole.

Execution of workflow

The following example uses small synthetic datasets so you can test it on your machine. Peak performance numbers quoted earlier vary by workload, dataset, implementation, and hardware. The largest speedups are individual operations and should not be read as an end-to-end pipeline factor.

This walkthrough covers five cuPhoton components: xDataReader, xRep, xPois, xFit, and xScan. The toolkit also includes xRay for time-domain X-ray detector analysis.

First, clone the repository and follow the README setup workflow. The setup command (uv sync --locked --extra dev --extra gpu --extra viz) provisions a locked Linux environment backing the NVIDIA CUDA 13 ecosystem: CuPy, PyTorch, Numba-CUDA, KvikIO, and NVIDIA nvCOMP, along with visualization libraries (Bokeh and Pillow). The same checkout includes the cuPhoton CLI, the modules used in this walkthrough, and examples/run_quickstarts.py, which generates synthetic inputs and writes run artifacts to quickstart-output/.

Set up a work directory

Follow the README to create the locked CUDA 13 GPU environment. xDataReader’s FITS path additionally needs its native extension; from the checkout run bash src/cuphoton/xdr/src/build.sh (xDataReader guide). cuPhoton 0.1.3 supports Python 3.12–3.14 and CUDA 13 on Linux.

Building xDataReader from source requires a C++17 compiler, CUDA and cuFile development headers, and a reentrant CFITSIO development installation. See the xDataReader guide for setup instructions.

git clone --branch v0.1.3 --depth 1 https://github.com/NVIDIA/cuPhoton.git
cd cuPhoton
export WORK_DIR="$PWD/blog-run"
mkdir -p "$WORK_DIR"/{fits,npz,figs}

Execute each Python snippet sequentially within a single interactive session (such as Jupyter, IPython, or an integrated script) starting from the root directory. Later sections reuse variables retained in memory: images, aligned, reference, target, fit, stamp, model, and result. Ensure WORK_DIR remains defined throughout execution. If the native xDataReader component is unbuilt, pause and execute bash src/cuphoton/xdr/src/build.sh.

The complete workflow script is available in the cuPhoton repository. For automated runs from a fresh terminal, invoke run_imaging_pipeline.py via uv run python. For an interactive notebook, open run_imaging_pipeline.ipynb in Jupyter.

Load FITS directly onto the GPU

The Flexible Image Transport System (FITS) is a standard data format in astronomy and astrophysics. Typically, you read the data with Astropy on the CPU and then call cupy.asarray to transfer the resulting array onto GPU memory. Parsing and decompression occur on the host, followed by a PCIe transfer before any science kernel runs. (At data rates of terabytes per hour, this sequence creates substantial data-movement overhead and is the first place the alert latency budget is consumed.)

To address this bottleneck, xDataReader accelerates decompression, reading, and loading FITS directly to the GPU. Specifically, xDataReader splits the workload as follows: CFITSIO plans byte ranges on the CPU, KvikIO reads them (NVIDIA GPUDirect Storage when available), nvCOMP decompresses GZIP tiles on the GPU, and CuPy holds the result. It supports uncompressed 2–16-axis image HDUs and two-dimensional images with GZIP_1 or GZIP_2 tile compression. Rice compression and dithered floating-point quantization are not supported.

GPUDirect Storage is available with a compatible driver and filesystem. On other CUDA 13 machines, KvikIO may take its PCIe compatibility path. The Python API does not change. batch_to_device still returns CuPy arrays; cuphoton.xdr.is_gds_active() reports the current path.

The public API is batch_to_device. It defaults to HDU index 1, the first extension, which must contain a supported image. For an image in the primary HDU, pass hdu_indices=(0,). The function returns one stacked CuPy array per requested HDU. Files grouped into one call must match in shape and dtype at each selected HDU.

Example:

The following code writes two synthetic 256×256 frames that mimic a reference template and a later science exposure: three static stars plus one fainter source that moved by a few pixels. The single moving source is the dipole you should recover after subtraction. The science frame is also a little blurrier, as if the seeing were worse. The example then shifts the science image to mimic a pointing offset and updates the header so both frames still describe the same sky. The two files do not share a detector pixel grid, but the World Coordinate System (WCS) in their headers maps both to the same sky coordinates.

The first Python block is scaffolding. Paste it once and keep the session open. It imports everything later sections need and sets WORK plus the scene constants those two FITS files were written with. The downloadable script (run_imaging_pipeline.py) writes the files with sky coordinates in the headers and defines the PNG helpers used at the end.

import os
from pathlib import Path

import numpy as np
from astropy.io import fits
from PIL import Image
from scipy.ndimage import gaussian_filter
from scipy.ndimage import shift as ndshift

from cuphoton.xdr import batch_to_device
from cuphoton.xfit import GaussianDipoleModel, fit_dipoles
from cuphoton.xpois import GaussianBasisComponent, solve_constant_kernel
from cuphoton.xrep import (
    BBox,
    Grid,
    build_stack_spec_from_fits,
    make_north_up_wcs,
    reproject_stack,
)

SHAPE = (256, 256)
STARS = (
    (70.0, 80.0, 45.0, 2.10),
    (175.0, 190.0, 32.0, 2.40),
    (40.0, 200.0, 22.0, 2.00),
)
MOVER_OLD = (120.0, 108.0, 14.0, 2.10)  # y, x, amp, sigma
MOVER_NEW = (126.5, 116.0, 14.0, 2.10)
CRVAL = (150.0, 2.0)
PIXEL_SCALE = 0.2  # arcsec / pixel
SCIENCE_SHIFT = (20.0, -32.0)  # dy, dx; large enough to see by eye
SEEING_SIGMA = 0.65  # extra Gaussian blur applied to the science frame

That scaffolding stays in memory. The product call is batch_to_device: it reads both image HDUs onto the GPU and returns a stacked CuPy array.

(images,) = batch_to_device(
    [WORK / "fits" / "template.fits", WORK / "fits" / "science.fits"],
    hdu_indices=(1,),
)

A successful run prints the array shape (2, 256, 256), its dtype, and the GPU device it landed on—confirming both frames made it onto the device. Those arrays are still in each file’s native pixel grid. The static stars should not line up; that is the planted mismatch, not a failure of the loader.

Interpolation and alignment: Align both visits on the same sky grid

Two visits to the same patch typically do not share a pixel grid. Dither, rotation, and plate scale put the same star on different (x, y) in the template and the science frame. The WCS in the FITS header is the map from those pixels to the sky. xPois matches point-spread functions on arrays that are already registered. Subtract first, and every static star becomes a residual; the mover is no longer a clean dipole.

xRep resamples two-dimensional images onto a shared north-up celestial grid. build_stack_spec_from_fits reads both WCS solutions and builds one destination footprint. reproject_stack warps every member onto that grid. backend="auto" prefers CuPy, then CUDA PyTorch, then CPU if no GPU is available. Interpolation is Lanczos-3 by default.

The next block keeps the FITS written above—native detector pixels plus the planted CRPIX offset—and puts both visits on one grid. After the warp, static stars should sit on the template. The mover should remain the only source that changed position in the sky.

native, spec = build_stack_spec_from_fits(
    [WORK / "fits" / "template.fits", WORK / "fits" / "science.fits"],
    hdu=1,
    grid=Grid.from_wcs(wcs_t),
    output_bbox=BBox(0, 0, SHAPE[1], SHAPE[0]),
    interpolation="lanczos3",
    mapping_grid_step=16,
)
aligned = reproject_stack(native, spec, backend="auto")

Print aligned.images.shape and aligned.backend. Expect a stack of two frames and a GPU backend on CUDA 13. If static stars still sit a few pixels apart, the WCS did not describe the shift you planted, or the destination grid did not cover both footprints.

Matching PSF and subtraction

A raw science-minus-template difference is dominated by residuals around every bright star caused by differences in atmospheric blurring (seeing). Optimal image subtraction, following Alard and Lupton, instead solves for a compact convolution kernel K and a differential background B so the reference template R acquires the science PSF and matches the science exposure, T:

T ≈ R ⊗ K + B,

where ⊗ denotes convolution. The resulting model is then subtracted from the exposure T to produce the difference D:

D = T − (R ⊗ K + B).

Using a Gaussian-polynomial kernel basis, xPois computes the convolution kernel required for optimal image subtraction. For the constant-kernel solve used here, backend="auto" selects the first available backend in this order: CuPy, Numba-CUDA, CPU. The small coefficient solve stays on the host; the pixels do not. Flux-conserving basis rewriting keeps the kernel sum near one. Check this sum when assessing flux conservation. A well-behaved kernel is compact and roughly symmetric. A bright ring or a large negative bowl means the basis cannot represent the seeing mismatch, or a saturated star entered the fit unmasked.

The next block keeps reference and target from the aligned stack and fits a 15×15 kernel. Because xRep puts the visits on one WCS, a compact kernel should absorb the seeing mismatch. The three static stars should cancel. The mover should remain as a positive lobe at t1 and a negative lobe at t0.

fit_mask = np.ones(SHAPE, dtype=bool)
fit_mask[105:142, 93:132] = False
fit = solve_constant_kernel(
    reference,
    target,
    [GaussianBasisComponent(sigma=SEEING_SIGMA, degree=0)],
    kernel_shape=(15, 15),
    background_degree=0,
    flux_conserve=True,
    fit_mask=fit_mask,
    backend="auto",
)

Print fit.backend, kernel.sum(), chi2, and fit_pixel_count. A kernel sum near 1.0 and a residual whose bright-star rings have collapsed means the match worked. If the stars still dominate D, either the images were not registered—use cuphoton.xrep to put them on a shared WCS—or the kernel basis is too narrow for the seeing difference. On a CUDA 13 GPU, fit.backend should be cupy and residual RMS should sit far below the peak star amplitude of about 45 in this scene.

Fit the dipoles

Moving objects, such as asteroids, and image-registration errors can produce dipole patterns. These are paired positive and negative flux regions in a difference image. Processing each candidate separately can become a bottleneck. A pipeline may need to solve thousands of least-squares optimization problems across small image stamps for each exposure.

In this example, xFit’s GaussianDipoleModel models each stamp as a difference between two rotated elliptical Gaussian components. Both components are evaluated at the same pixel coordinates (x,y) and share amplitude, widths, and orientation. The two centers (x⁺, y⁺) and (x-, y-) vary independently:

m(x, y) = G(x, y; x⁺, y⁺) − G(x, y; x⁻, y⁻)

xFit also supports models based on supplied sampled PSFs through StampDipoleModel.

xFit’s Levenberg–Marquardt solver factors many regularized Gauss–Newton systems at once on CuPy. Converged fits stop participating in subsequent solver iterations, reducing unnecessary work. The synchronous call returns the batch results after the remaining fits finish. Widths are optimized in log space and returned as positive pixel standard deviations; the analytic Jacobian is the default for the Gaussian model. Returned arrays are portable NumPy—parameters, covariance, and standard errors—so a classifier never imports the CuPy solver.

The following code cuts a 21×21 stamp around the planted mover in the xPois residual and starts the fit from a slightly wrong initial guess. Print the recovered centers alongside the known input positions. Compare the recovered centers with the known input positions to check fit accuracy.

STAMP = 21
half = STAMP // 2
cy, cx = 123, 112  # midpoint of the planted mover
stamp = fit.residual[cy - half : cy + half + 1, cx - half : cx + half + 1]

# Seed the fit from the observed stamp extrema, in centered x/y coordinates.
y_pos, x_pos = np.unravel_index(np.argmax(stamp), stamp.shape)
y_neg, x_neg = np.unravel_index(np.argmin(stamp), stamp.shape)
initial = np.array(
    [
        [
            stamp.max(),
            2.0,
            2.0,
            0.0,
            x_pos - half,
            y_pos - half,
            x_neg - half,
            y_neg - half,
        ]
    ]
)

model = GaussianDipoleModel((STAMP, STAMP), dtype=np.float64)
result = fit_dipoles(
    stamp[None], model=model, initial=initial, backend="auto"
)

Expect converged: True on this synthetic data. Recovered (x_pos, y_pos) and (x_neg, y_neg) should sit within a fraction of a pixel of the planted offset. A structured leftover in the fit residual usually means the stamp still contains a static star or the initial centers were on the wrong lobe.

Review the stamps

The examples above validate the workflow on a controlled synthetic scene with one planted moving object and known ground truth. Real observations provide neither that simplicity nor that certainty. A single telescope pointing (visit) can generate up to 10,000 candidates, and a full night yields millions, making individual inspection impractical. The next step is to prioritize candidates for human review.

cuphoton xscan review-queue ranks the candidates—by model uncertainty, by known errors, or as a straight dataset audit—and limits the number of candidates for human review. Use cuphoton xscan review-bokeh to inspect the search exposure (current observation), the reference template, the raw difference, and the xPois residual, with Real / Bogus / Unsure labels for each candidate.

Figure 9 shows the review interface with an example noise artifact, so you can inspect the layout before using it with an observation. Observational FITS data are not in the git repository. Copy an example under examples/xscan/, set the local paths, and follow the README plus the xScan guide to build a review queue and launch the app.

Scale the same workflow across multiple GPUs

The walkthrough above is one image pair on one GPU, small enough for an NVIDIA DGX Spark or a workstation. NVIDIA cuPhoton runs the same loading and processing path across multi-GPU, multi-node systems. After registration, cuPhoton distributes complete image pairs and their candidate coordinates across GPU workers, keeping image arrays on the device through subtraction, dipole fitting, and classification. This example uses synthetic inputs and an untrained model to verify a multi-GPU launch. Scaling that path is how the toolkit handles survey-sized streams such as Vera C. Rubin Observatory exposures. To run many pairs, use the following commands. Environment setup, Slurm and SSH launches, and the Dragon alternative are in the distributed execution guide.

python examples/distributed-pipeline/prepare_example.py \
  --output "$CUPHOTON_RUN_ROOT/input" --images 8

export CUDA_VISIBLE_DEVICES="${CUDA_VISIBLE_DEVICES-0,1,2,3}"
mpiexec -n 8 --map-by ppr:4:node --bind-to none \
  -x CUDA_VISIBLE_DEVICES \
  cuphoton-openmpi-rank-exec -- \
  cuphoton xscan run-pipeline --executor mpi \
  --manifest "$CUPHOTON_RUN_ROOT/input/pipeline.json" \
  --output-dir "$CUPHOTON_RUN_ROOT/runs" --name blog-mpi \
  --warmup-rounds 1 --measure-rounds 1

--images 8 creates eight image pairs, and -n 8 starts eight MPI workers, mapped one per GPU. The pairs must already share a pixel grid, which is the xRep step above. On a machine with one GPU per node, set both numbers to the node count and use --map-by ppr:1:node. Those commands come from the distributed-pipeline example, which also includes the host-file templates. The generated model only checks that the workers start and finish.

Get started

The examples were developed on NVIDIA DGX Spark and target Linux systems with compatible NVIDIA GPUs and drivers. Follow the cuPhoton README to clone the repository, check the supported Python and CUDA versions, and set up the environment.

• Build the native xDataReader extension for FITS input: xDataReader guide.
• Run data loading and processing across multiple GPUs and nodes: distributed execution guide.
• Create a review queue and launch the candidate review interface: xScan guide.
• Explore component APIs and output artifact specifications: cuPhoton documentation.

For related examples, read Accelerated X-Ray Analysis for Nanoscale Imaging (XANI) and Using Accelerated Computing to Live-Steer Scientific Experiments at Massive Research Facilities.

Watch the NVIDIA cuPhoton product overview for an introduction to the toolkit.