NVIDIA cuPhoton Accelerates Scientific Image Analysis Pipelines on GPUs

NVIDIA cuPhoton Accelerates Scientific Image Analysis Pipelines on GPUs NVIDIA cuPhoton is an open-source CUDA-X toolkit built to eliminate CPU bottlenecks in high-throughput scientific imaging by executing pipelines—from raw disk I/O through dipole fitting and transient classification—entirely on GPUs. cuPhoton accelerates image reading by up to 14,900x and signal processing by up to 14,550x over an x86 CPU baseline on representative Grace Blackwell and Vera Rubin workloads. The modular pipeline includes xDataReader for GPUDirect and GPU-decompressed I/O, xRep for celestial grid reprojection, and xPois for optimal image subtraction. Candidate analysis is powered by xFit's concurrent Levenberg–Marquardt dipole solver and xScan's interactive review tools, scaling across multi-GPU and multi-node clusters via MPI. Throughput Challenges and Acceleration Benchmarks Modern high-throughput instruments, such as the Vera C. Rubin Observatory's LSSTCam, produce petabytes of multidimensional data that quickly overwhelm traditional CPU pipelines. LSSTCam records a 3.2-gigapixel exposure every 39 seconds, generating up to 20 terabytes and 10 million candidates per night that must be processed within 60 to 120 seconds. Testing on multi-GPU workloads across hundreds of terabytes on NVIDIA Grace Blackwell and Vera Rubin systems demonstrated that cuPhoton accelerated image reading by up to 14,900x and signal processing by up to 14,550x, shrinking nine months of CPU-bound analytics down to four hours. GPU-Resident Data Ingestion and Image Subtraction cuPhoton keeps data resident on the device through modular components. To bypass host decompression bottlenecks, xDataReader uses CFITSIO for CPU byte-range planning, KvikIO for GPUDirect Storage transfers, and nvCOMP for GPU-based GZIP tile decompression. Images are aligned to celestial grids using xRep with Lanczos-3 interpolation. Following alignment, xPois performs optimal image subtraction to compute convolution kernels and differential background offsets; pixel data remains on the GPU while the smaller coefficient linear system is solved on the host. Dipole Fitting, Transient Inspection, and Distributed Execution For transient candidate analysis, xFit models difference image stamps using Gaussian or PSF-sampled dipole models, executing concurrent Levenberg–Marquardt optimization on CuPy. xScan provides candidate review queues and an interactive Bokeh interface for human-in-the-loop validation. cuPhoton scales across multi-GPU and multi-node environments by distributing pre-aligned image pairs across MPI workers, maintaining GPU residency through subtraction, fitting, and candidate classification.

NVIDIA推出开源科学图像分析工具包 cuPhoton NVIDIA发布开源CUDA-X工具包cuPhoton，旨在消除天文与时域科学仪器图像分析中的CPU性能瓶颈，通过全GPU管线加速从数据加载、图像减差到偶极拟合与候选目标分类的完整处理流程。 针对现代天文与X射线等高通量仪器的海量数据，cuPhoton将端到端分析保留在GPU上，在代表性多GPU工作负载中将9个月的数据分析缩短至4小时。 包含xDataReader、xRep、xPois、xFit、xScan和xRay核心模块，支持FITS瓦片解压、重投影、最优化图像减差及GPU偶极拟合。 cuPhoton 0.1.3基于Linux，支持Python 3.12–3.14与CUDA 13，可结合MPI在多节点多GPU集群上实现横向扩展。 性能瓶颈与cuPhoton定位 高通量仪器产生的数据量往往超出CPU流水线处理能力，以智利薇拉·鲁宾天文台为例，LSSTCam每39秒拍摄32亿像素曝光，需在60至120秒内对约1万个候选体完成分类，单夜处理量可达20 TB。cuPhoton通过GPU全流程常驻设计，在NVIDIA Grace Blackwell与Vera Rubin系统上测试显示，图像读取速度提升高达14,900倍，信号处理提升达14,550倍，大幅缩短端到端科研分析耗时。 核心模块与技术实现 xDataReader通过CFITSIO规划CPU字节范围、KvikIO读取数据并借助nvCOMP在GPU端解压GZIP瓦片，绕过主机解压和PCIe瓶颈；xRep通过CuPy或PyTorch后端将图像重采样至共享天球网格；xPois利用高斯多项式基底求解卷积核与差异背景实现图像相减；xFit则利用CuPy并发执行Levenberg–Marquardt算法进行偶极拟合，并向主机返回NumPy数组；xScan提供了审查队列及交互式标注工具。 系统支持与多节点扩展 cuPhoton 0.1.3适用于Linux环境下的Python 3.12–3.14和CUDA 13，源码构建xDataReader需C++17、cuFile开发头文件和重入CFITSIO。工具包支持通过OpenMPI在多GPU及多节点集群上分发图像对和候选坐标，并在DGX Spark等系统上实现了多卡并行处理管线。

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. Figure 1. NVIDIA cuPhoton package components 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. Figure 2. (Top panel) Every night the Vera C. Rubin Observatory telescope captures 20 TB of data (Credit: Original images from NSF-DOE Vera C. Rubin Observatory/NOIRLab/SLAC/AURA, modified using AI). (Bottom panel) NVIDIA cuPhoton modules accelerate the Rubin Observatory processing and analysis workflow 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. Figure 3. NVIDIA cuPhoton speedup on NVIDIA GB200 NVL72 compared to an x86 CPU baseline across 1, 32, and 64 GPUs. Operations scale significantly with multi-GPU execution 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. Load the current exposure and reference data into GPU memory. Align: Place both visits on a shared sky grid so the same star falls on the same pixel. 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. Fit a dual-PSF model to the dipole to estimate the moving object’s positions. Visualize the results: display the PSF-matching kernel, subtraction residual, and recovered dipole . Figure 4. Scientific workflow from a new exposure to an AI-ready product 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. Figure 5. Reference template at t0 (left) and science exposure at t1 (right) before alignment. The science WCS offset prevents static stars from registering Figure 6. Reference template at t0 (left) and science exposure at t1 (right) after reproject_stack . Static stars register while the planted moving source remains displaced 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. Figure 7. xPois products: fitted kernel K, matched template R⊗K+B, and diverging residual D. Static stars should cancel; the planted mover should appear as a dipole 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. Figure 8. Difference stamp, xFit model, and fit residual. The recovered lobe centers should match the planted mover at t0 and t1 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 x Scan guide to build a review queue and launch the app. Figure 9. xScan review ( cuphoton xscan review-bokeh ) showing search, template, difference, and xPois panels and Real, Bogus, and Unsure labeling controls. Run the same command on your LSSTComCam or HSC products; the git tree does not ship those pixels 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.