avior-cuda

A CUDA image denoiser. No third-party libraries: the kernels, the PNG codec, the DEFLATE implementation, the quality metrics and the benchmark harness are all in this repo. Not OpenCV, not NPP, not Thrust, not libpng, not zlib. Run it on one image or on a whole directory.

Every filter exists twice: once as a CUDA implementation and once as a multithreaded CPU reference. The reference is not a fallback, it is the definition — the test suite runs both on the same input and checks that they agree, and the benchmark uses the CPU as its baseline so the speedups below are measured against a real parallel implementation rather than a strawman.

The main filter is non-local means, and it is written twice again on each backend. One version is the direct formulation straight from the paper; the other replaces the patch loop with a box sum, which on the CPU removes the dependence on patch size completely.

Getting that second one onto a GPU is the interesting part of this project, and it does not survive the trip intact. The CPU does it with a summed-area table, which is the one shape in this problem that a GPU is bad at, and the thing you can build instead is a step worse asymptotically. It is also, at the default settings, worth far less than the CPU numbers would lead you to expect — 1.4x rather than 4.9x — and at small patch sizes it is a net loss. That is measured below, along with why.

Gallery

Kodak kodim07 from the Lossless True Color Image Suite with Gaussian noise at sigma 0.08, every filter run on the same frame and the same GPU, each one given its own best parameters rather than one shared setting. The source image, the noisy version and all five results are in data/kodak.

Reference	Noisy input	Non-local means
ground truth	22.04 dB · 0.4699	31.12 dB · 0.8916

Those are full frames at half size. The differences live at a smaller scale, so the crops below are 1:1 pixels from the hibiscus at (520, 330), 200 by 156.

Reference	Noisy input	Gaussian	Median
ground truth	22.04 dB · 0.4699	29.26 dB · 0.8301	28.01 dB · 0.7635

Bilateral	Non-local means	Non-local means, no summed-area table
30.66 dB · 0.8775	31.12 dB · 0.8916	31.12 dB · 0.8916

filter	settings	PSNR	SSIM	GPU	CPU
noisy input	sigma 0.08	22.04 dB	0.4699
Gaussian	--filter gaussian --sigma 0.9	29.26 dB	0.8301	0.27 ms	25.97 ms
median	--filter median --radius 1	28.01 dB	0.7635	0.79 ms	204.52 ms
bilateral	--filter bilateral --sigma 2.0 --sigma-range 0.25	30.66 dB	0.8775	0.93 ms	690.56 ms
non-local means	--filter nlm --h 0.085 --patch 2 --search 7	31.12 dB	0.8916	14.25 ms	3656.08 ms
non-local means, no SAT	--filter nlm-fast --h 0.085 --patch 2 --search 7	31.12 dB	0.8916	9.88 ms	1012.75 ms

Times are one run of the CLI on this 768x512 frame: the GPU column is kernel time on a Tesla T4, the CPU column is the same filter on four threads of the Xeon. The last two rows are the point of the project: the same estimator, agreeing to 5e-07 per sample, 1.44x apart on the GPU and 3.6x apart on the CPU.

Everything above was produced by the binary itself; docs/gallery.sh is the script that does it, and re-running it on a clean checkout reproduces all ten PNGs byte for byte. That is the determinism claim below, end to end and including the PNG encoder.

build/avior noise   --input data/kodak/kodim07.png \
                    --output data/kodak/kodim07_noisy.png \
                    --type gaussian --strength 0.08 --seed 7
build/avior denoise --input data/kodak/kodim07_noisy.png \
                    --output data/kodak/output/nlm-fast.png \
                    --reference data/kodak/kodim07.png \
                    --filter nlm-fast --h 0.085 --patch 2 --search 7
build/avior crop    --input data/kodak/output/nlm-fast.png \
                    --output docs/gallery/nlm-fast_crop.png --rect 520,330,200,156

Build

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build -j
ctest --test-dir build --output-on-failure

Needs a CUDA toolkit and CMake 3.18 or newer. With CMake 3.24 or newer the build targets whatever card is in the machine (CMAKE_CUDA_ARCHITECTURES native); on older CMake it builds a fat binary for sm_75 through sm_89, sm_70 being absent because CUDA 13 dropped it. Pass -DCMAKE_CUDA_ARCHITECTURES=89 to pin one.

-DAVIOR_CUDA=OFF builds the CPU reference alone, with no toolkit installed and no device present. That build still runs the whole test suite — it reports the device tests as skipped rather than passed, and --device cuda fails with a message instead of quietly running somewhere else. -DAVIOR_NATIVE_ARCH=ON adds -march=native to the host code, -DAVIOR_SANITIZE=ON adds ASan and UBSan to it.

The build does not turn on --use_fast_math. Two kernels want the fast transcendental path and they ask for it by calling __expf directly, which keeps the reciprocal and the divisions in the rest of the code IEEE.

Running it

Point --input at a directory and every image in it is denoised into <directory>/output:

build/avior denoise --input photos/ --filter nlm-fast --h 0.08

filter    nlm-fast
device    cuda:0 Tesla T4 (40 SMs, sm_75)
images    4 -> /tmp/photos/output

kodim01.png                768x512      15.86 ms     24.8 Mpix/s  (1.16 up / 9.91 kernel / 1.07 down)
kodim07.png                768x512      14.67 ms     26.8 Mpix/s  (1.10 up / 9.54 kernel / 0.99 down)
kodim13.png                768x512      12.74 ms     30.9 Mpix/s  (1.08 up / 9.54 kernel / 0.98 down)
kodim21.png                768x512      12.91 ms     30.5 Mpix/s  (1.08 up / 9.53 kernel / 0.99 down)

total     56.17 ms for 4 images
wrote     all outputs

Point it at one image and the result goes next to it with _output on the name, so shot.png becomes shot_output.png:

build/avior denoise --input shot.png --h 0.08

--output overrides either behaviour, taking a file or a directory, and any missing parent directories get created. A rerun on a directory ignores its own output folder, so it is safe to repeat.

--device cpu runs the same command on the CPU reference, and --device cuda is the default whenever a device is present. --gpu <n> picks between several.

build/avior devices

cuda:0  Tesla T4
        sm_75, 40 SMs, 1024 threads/block, 48 KB shared/block
        14.6 GB memory

The other commands work the same way:

build/avior noise   --input photos/ --type gaussian --strength 0.08
build/avior crop    --input shot.png --output detail.png --rect 520,330,200,156
build/avior metrics --reference clean.png --input shot_output.png
build/avior bench   --filter nlm-fast --width 1024 --height 1024 --runs 5
build/avior synth   --output target.png --width 512 --height 512

avior --help lists every option.

Formats

PNG in and out, 8- or 16-bit, grayscale, RGB, RGBA or paletted, plus binary PNM. The output format follows the extension you ask for. PNG support is first-party: src/deflate.cpp is an inflate/deflate implementation and src/png.cpp is the chunk, filter and colour handling, so the no-dependencies claim still holds and there is nothing to install besides the toolkit. Interlaced PNGs are rejected with a clear message rather than decoded wrongly.

Numbers

Tesla T4 (sm_75, 40 SMs, 16 GB, 48 KB shared per block), CUDA 13.0, against an Intel Xeon Platinum 8259CL at 2.50 GHz, g++ 13.3, -O3.

Read the speedup column with the CPU in mind: this box has two physical cores, so cpu x4 is four hyperthreads on two cores and the CPU baseline is a small one. The honest way to read the table is the split, not the ratio — for the cheap filters the kernel is already a rounding error and what is left is the bus.

Every filter on a 1024x1024 RGB frame, best of five. cpu x1 is the baseline for the speedup column, and the split shows where the GPU time actually went.

filter	cpu x1	cpu x4	cuda	up / kernel / down	vs cpu x1	vs cpu x4
Gaussian	145.9 ms	76.3 ms	7.25 ms	2.81 / 0.41 / 2.33	20.1x	10.5x
median	166.0 ms	124.3 ms	7.46 ms	2.77 / 0.28 / 2.63	22.2x	16.7x
bilateral	2770.3 ms	1326.7 ms	8.49 ms	3.05 / 0.93 / 2.98	326.2x	156.3x
non-local means	23746.5 ms	9859.8 ms	45.45 ms	2.83 / 39.36 / 2.48	522.5x	216.9x
nlm, no SAT	4881.1 ms	2492.9 ms	34.55 ms	3.05 / 28.09 / 2.72	141.3x	72.2x

The split is the interesting column. For the three cheap filters the kernel is 0.28 to 0.93 ms against 5.1 to 6.0 ms of PCIe, so 87 to 95 per cent of the timed GPU work is the bus and the "speedup" is really a statement about how long a 1024x1024 frame takes to cross it. Only the two NLM rows are actually compute-bound, and they are the only two where the number means what it looks like it means.

(The three legs are CUDA-event times averaged over the five runs; best(ms) is best-of-five wall clock and also carries the host-side allocate-and-decode work outside them, which is why the two do not add up.)

Frame size, fast NLM, patch radius 2, search radius 7, against four CPU threads:

frame	cpu x4	cuda	up / kernel / down	Mpix/s	transfers	speedup
256x256	135.8 ms	2.18 ms	0.30 / 1.48 / 0.28	30.1	27%	62.4x
512x512	569.4 ms	8.02 ms	0.84 / 6.12 / 0.73	32.7	20%	71.0x
1024x1024	2571.0 ms	33.49 ms	2.84 / 28.26 / 2.54	31.3	16%	76.8x
2048x2048	10837.0 ms	170.34 ms	11.78 / 117.53 / 10.98	24.6	13%	63.6x
4096x4096	43853.6 ms	717.52 ms	42.99 / 504.92 / 43.39	23.4	12%	61.1x

This filter is heavy enough that it never becomes transfer-bound: the copies fall from 27 per cent of the call to 12 as the frame grows, and throughput stays between 23.4 and 32.7 Mpix/s across a 256-fold range in pixel count. The decline at the top end is a 70 W passively cooled card under a sustained half-second kernel, plus a working set that stopped fitting the 4 MB L2 some time ago.

Patch size, 512x512 RGB, search radius 7. CUDA columns are kernel time, CPU columns are one thread. Not caring how big the patch is was the whole point of the second NLM, and it is the one property that did not survive the move to the GPU — the four columns below are four different growth laws:

patch radius	window	cpu x1 nlm	cpu x1 nlm-fast	cuda nlm	cuda nlm-fast
1	3x3	2622 ms	1216 ms	3.45 ms	5.14 ms
2	5x5	5908 ms	1214 ms	9.28 ms	5.99 ms
4	9x9	18901 ms	1202 ms	31.30 ms	8.69 ms
6	13x13	39278 ms	1214 ms	65.46 ms	12.14 ms

• cuda nlm divided by window area is 0.383, 0.371, 0.386, 0.387 — quadratic, to four per cent over a 19-fold range in area;
• cpu x1 nlm-fast varies by 1.2 per cent across the whole sweep — the summed-area table really is O(1) in the patch;
• cuda nlm-fast is neither: it fits 2.2 ms + 0.77 ms per unit of patch side, which is the separable box sum being linear where the table was constant.

What that costs in practice is in When the trick stops paying.

Search radius, 512x512 RGB, patch radius 2, against four CPU threads. Every displacement has to be visited by both formulations, so all four columns are linear in the number of displacements and no backend or algorithm changes that — what changes is the constant:

search radius	displacements	cpu x4 nlm	cpu x4 nlm-fast	cuda nlm	cuda nlm-fast
3	49	518 ms	183 ms	2.02 ms	1.54 ms
5	121	1465 ms	306 ms	4.97 ms	3.38 ms
7	225	2655 ms	590 ms	9.15 ms	6.06 ms
10	441	4774 ms	1256 ms	18.40 ms	12.05 ms

Dividing the two CUDA columns by the displacement count gives 0.041 ms each time for the direct version, to within 2.5 per cent over a 9-fold range, and 0.027 to 0.031 for the separable one. This is the axis where the GPU has nothing clever to offer: it is the same work either way.

How it works

The context

CudaContext owns the device, one stream, the event pair used for timing, and a small set of scratch slabs. The slabs grow to the high-water mark of the run and are then handed back out, so filtering the second image in a directory costs a kernel launch rather than a cudaMalloc. This is worth having but it is not dramatic, and it is worth being exact about which: on this card a cudaMalloc of the row-sum buffer measures 0.06 to 0.07 ms whether it is 15 MB or 240 MB, so reusing the four slabs a filter needs saves a couple of tenths of a millisecond per image. Against a 9.9 ms NLM kernel that is noise; against the 0.27 ms Gaussian kernel it is not.

The header the rest of the program sees is plain C++. cuda_runtime.h is included only by the .cu files, the implementation type is a pimpl, and the CLI and the test binary are ordinary host-compiler translation units. That is what lets -DAVIOR_CUDA=OFF swap in a 50-line stub and still build everything.

Padding

Filters that scan a neighbourhood work on a border-replicated copy, so the inner loop has no bounds tests and, more to the point on a GPU, no branch that only the edge warps take. The pad is built by a kernel rather than on the host, so only the original image crosses PCIe.

Reference NLM

One thread per output pixel, the paper's formulation transcribed. Every thread in a block walks the same search window over heavily overlapping data, so the block stages its whole footprint — its tile plus a halo of patch + search on each side — into shared memory once, and the inner loops never touch global memory again. At the default radii that is a 50x34 tile of RGB for a 32x16 block, about 20 KB.

The staged footprint grows as (block + 2(patch + search))^2, so a large enough search radius stops fitting. The host works out the shared-memory footprint before it launches, picks the largest block shape that fits, and falls back to a variant that reads global memory directly if none does. Both variants are the same arithmetic in the same order, so which one ran is not observable in the output.

An accident of the layout: with three interleaved channels, neighbouring threads read shared memory three floats apart, and 3 is coprime with 32, so every lane lands in a different bank. The natural indexing is already conflict-free.

Fast NLM, and the table that had to go

The CPU version fixes a displacement, builds a summed-area table over the squared-difference image for that displacement, and then reads every patch SSD out of it as a four-corner difference. Cost per pixel stops depending on the patch size entirely, which is the whole point, and the measurements below show it holding exactly: 1216, 1214 and 1202 ms at patch radius 1, 2 and 4.

A SAT is a prefix sum, and it is the one shape in this problem that a GPU is bad at. It is sequential along both axes; a parallel scan wants the fast-moving index to be the thread index, which fights the coalescing you want for the rest of the kernel; and the accumulation has to be done in double, because a float table over a megapixel image loses enough low bits by the far corner that the four-corner difference stops being accurate, which is why cpu::nlm_filter_fast accumulates its table in double. Doing that on a device is a scan kernel per axis per displacement, in double, 225 times.

So this version drops the table and uses the other property of a box sum: it is separable.

• one kernel writes, for a fixed dy and every dx at once, the horizontal patch_side-wide sum of squared differences at every position;
• one kernel walks those vertically, which finishes the box sum, turns each SSD into a weight, and accumulates.

This is not the same asymptotics, and the numbers say so. A summed-area table is O(1) in the patch area. A separable box sum computed this way is O(patch_side) — it removes one of the two patch dimensions, not both. The direct version is O(patch_side^2). All three show up cleanly in the patch-radius table under Numbers: cuda nlm is the quadratic, cpu x1 nlm-fast is the flat line, and cuda nlm-fast — the linear one in between — is mine.

Getting back to genuine O(1) means putting a running sum back into each pass, which is the sequential shape the GPU dislikes and which would cost the coalesced access pattern the current kernels enjoy. On this hardware that trade did not look worth making: at the patch radii anyone actually uses, the separable version is already several times faster than the direct one, and the remaining factor of patch_side is small next to the patch_side^2 it removed.

Two things did survive the move, and they are the ones worth having. Nothing is ever summed over more than patch_side terms, so the whole thing stays in float and is in fact more accurate than the double-precision table it replaces: the table's four-corner difference subtracts two large numbers to get a small one, and this never forms the large numbers at all. And the second kernel keeps its accumulators in registers across the whole row of displacements, so the sum and weight planes in global memory are read and written once per dy instead of once per displacement — fifteen times less traffic at the default search radius, on the buffer that is touched most.

The row-sum buffer is the one allocation that scales with the search radius. If it would not fit a fixed budget the image is cut into horizontal bands and the whole thing runs once per band; one band is the common case, and the 4096x4096 frame at the bottom of the size table takes four.

When the trick stops paying

On a CPU the integral-image formulation is simply better, at every patch size, by a lot. On a GPU that is not true, and the crossover is not where you would guess. The patch-radius table under Numbers has the measurements; this is what they mean.

At patch radius 1 the direct version wins on the GPU, by 1.5x: the separable version is doing more global memory traffic than the arithmetic it saves is worth. They cross between radius 1 and 2, and radius 2 happens to be the default. By radius 4 the separable version is ahead by 3.6x and by radius 6 by 5.4x — while on one CPU thread the same two radii give it 15.7x and 32.4x.

The reason is that the two machines are short of different things. The direct kernel is almost pure ALU work on data that a block stages once into shared memory, and a T4 has 8 TFLOP/s of fp32 to spend and only 320 GB/s of bandwidth. The separable version spends bandwidth to save arithmetic. That is exactly the right trade on a CPU, where the arithmetic is scarce and the caches are close, and close to a wash on a GPU at small patch sizes.

The practical consequence, at the default settings this README quotes everywhere: the whole integral-image apparatus buys 1.40x on the GPU (39.4 ms of kernel down to 28.1 at 1024x1024) against 4.87x on one CPU thread. It is still the right default, but it is not the landslide the CPU numbers suggest, and a port that assumed the CPU's answer carried over would have been wrong.

The lookup table that isn't

The bilateral filter's range term is one exp per neighbour per pixel. On a CPU the standard fix is a lookup table, and it is a real win. On a GPU it is the wrong move twice over: __expf is a couple of instructions on the special function unit, and a table index that depends on the pixel data makes the threads in a warp read different addresses, which is the access pattern shared and constant memory are worst at. Both weights are computed here, neither is looked up.

The CPU reference does not use a table either — it is the reference, and it calls std::exp so that it is obviously right. So the bilateral row of the benchmark flatters the GPU by whatever a well-implemented CPU table would have saved. That is not nothing, and it is measured rather than assumed: building the CPU-only predecessor of this code and running the same 1024x1024 frame through both, the table is worth 1.76x on one thread and 1.36x on fourteen — on a 14-thread Apple laptop rather than the Xeon above, so treat it as the right order of magnitude and not a fourth significant figure. Read the bilateral speedup with that in mind. The NLM rows, which are what the project is about, have no table on either side.

The Gaussian coefficients are the opposite case and do live in constant memory: every thread in a warp reads the same tap at the same time, which is exactly what constant memory broadcasts for free.

Median

Two strategies, picked by width. At radius 1 and 2 the window is 9 or 25 values and it stays in registers, sorted by an odd-even transposition network: N passes over compile-time index pairs, so there is no data-dependent addressing to force the array into local memory, and no divergence. Wider than that and the register file gives out, so each thread takes a slice of shared memory strided by the block size — lane k always lands in bank k — and runs a selection sort stopped at the middle element. The host picks the block size so the slices fit, and says so plainly if they cannot.

The median is the one filter where the two backends have to agree exactly: it selects an actual sample out of the window rather than computing a new value, so there is no rounding to differ over. The test asserts equality with tolerance zero.

Determinism

Every filter here is bit-identical run to run, and independent of block size and of how the blocks happen to get scheduled. There are no atomics anywhere.

For the four per-pixel filters that is easy: one thread owns one output pixel and computes it start to finish. Fast NLM is the one that could have gone the other way — the obvious parallelisation is over displacements, which is what the CPU version does, and on a GPU that would mean atomicAdd into the accumulator planes and a different summation order every run. Instead the displacement loop is on the host, in a fixed order, and a pixel's accumulators are only ever touched by the thread that owns it. The test runs each filter three times and asserts the results are byte-for-byte equal.

Agreement with the CPU

Measured on a T4, over both channel counts and a frame whose dimensions are multiples of nothing:

filter	worst absolute difference
median, any radius	0 (exact)
Gaussian, bilateral, both NLMs	4.8e-07

For scale, one step of 8-bit quantisation is 3.9e-03, so the two backends disagree by about one part in eight thousand of the smallest value the output format can represent. The tests assert 1e-5, which is 20x headroom for another architecture without letting a real regression through, and exactly 0 for the median.

The two backends differ only by __expf against std::exp, and by fused multiply-adds the two compilers contract differently. The summation orders were kept the same on purpose, including the order of the innermost patch loop, which is why the differences stay at that level instead of at the level of float associativity.

Layout

include/avior/   image, cuda context, filters, io, png, deflate, metrics, noise
src/             codecs, metrics, noise, thread pool, benchmark harness
src/cpu/         the CPU reference, split into basic and nlm
src/cuda/        context and scratch, basic kernels, nlm kernels
src/cuda_stub.cpp  built instead of src/cuda when -DAVIOR_CUDA=OFF
app/main.cpp     CLI: denoise, noise, synth, crop, metrics, bench, devices
tests/           test binary, no framework
data/kodak/      sample inputs and their outputs
docs/gallery/    the crops shown above; gallery.sh and sanitize.sh make them

Tests

tests/test_avior.cpp covers the host side — the thread pool (full index coverage under uneven chunking, slot exclusivity, exception propagation, reuse after a throw), DEFLATE round-tripping over empty, repetitive, random and mixed input, PNG round-tripping at 8 and 16 bits in gray and RGB, PNG and PNM decoding to identical pixels, output path planning for both the single-file and directory cases, PSNR and SSIM edge cases, noise reproducibility and argument validation — and then the part that matters here:

• every filter, CUDA against CPU, at one and three channels, on a frame whose dimensions are deliberately not multiples of any block size, so the partial blocks and the halo clamping get exercised;
• every filter run three times on the GPU, asserting byte-for-byte equality;
• CUDA fast NLM against CUDA reference NLM at four patch radii, including radius 4, where the two have completely different cost;
• context reuse across a small frame, then a large one, then the small one again, asserting the scratch slabs grew, were reused, and did not change the answer;
• fast NLM run at four progressively smaller band budgets, down to one output row per band, asserting every one is byte-for-byte equal to the single-band result — banding changes which pixels share a kernel launch, never the order in which one pixel's contributions are summed, so anything short of exact equality would be a bug;
• the timing split adds up and every leg is non-zero;
• invalid arguments and out-of-range device indices throw.

Without a device the CUDA tests report skipped, not ok, so a green run on a machine with no GPU cannot be mistaken for a verified backend.

Clean under compute-sanitizer:

sh docs/sanitize.sh

Every distinct kernel path — the register median network, the shared-memory median, the tiled reference NLM, the two fast-NLM kernels, the separable Gaussian and the bilateral — run under memcheck, racecheck, initcheck and synccheck, at one and three channels. 64 runs, and every one of them reports 0 errors and 0 hazards.

and the host code is clean under ASan, UBSan and TSan:

cmake -S . -B build-tsan -DAVIOR_CUDA=OFF -DCMAKE_CXX_FLAGS="-fsanitize=thread -g" \
      -DCMAKE_EXE_LINKER_FLAGS="-fsanitize=thread"
cmake --build build-tsan -j && ./build-tsan/avior_tests