CUDA image convolution and edge detection

This capstone combines a 5 by 5 Gaussian blur, a 3 by 3 Sobel edge detector, and a threshold. The output is a PGM image with white coastlines on a black background.

Capability

Given the shipped program and one CUDA GPU, you will account for each stage's global traffic, add a two-kernel separable blur, and validate it against the CPU reference. Your work passes when both required image sizes match, the timing row reports the correct byte count, and you explain the measured ratio against a prediction written before the run.

Submit blurRows, blurCols, one timing-table row, and the predicted ratio beside the measured ratio. Allow 40 to 60 minutes for the independent work.

Retrieve, then predict

Use these facts from the previous lessons:

1. Day 11: one warp reading 32 consecutive floats coalesces across 128 bytes and requires four 32-byte sectors.
2. Day 14: every thread that has not exited must reach a block barrier, which orders earlier shared-memory writes before later reads.
3. Day 15: shared-memory word w maps to bank w % 32; a warp reading consecutive words has no bank conflict.
4. Day 18: constant memory serves a warp efficiently when every lane reads the same address.

Write four predictions before you run the program:

• How many global values must a 32 by 32 output tile load for a 5 by 5 filter?
• What traffic ratio follows from 13 bytes per pixel for staged Sobel plus threshold and 5 bytes per pixel for the combined kernel?
• Should changing a shared row pitch from 36 to 37 help when each warp reads along a row?
• A one-kernel separable blur moves 8 bytes per pixel. What ratio do you predict for a two-kernel version that moves 16?

Do not use the results table yet. The run should test these predictions, not supply them after the fact.

Minimal model

A halo, also called an apron, is the border a block loads around its output tile. A 5 by 5 filter has radius 2, so a block producing 32 by 32 outputs needs a 36 by 36 shared tile.

A separable filter is a 2D filter that equals the outer product of two 1D filters. The Gaussian here can run as a 5-tap row pass followed by a 5-tap column pass, reducing 25 multiply-adds per pixel to 10.

The Sobel operator estimates horizontal and vertical image change with two 3 by 3 filters. This program computes hypotf(gx, gy) and writes 255 when that magnitude exceeds the threshold, otherwise it writes 0.

The filter coefficients live in constant memory. Every lane executes the same loop index, so a warp requests one coefficient address at a time and uses the broadcast path measured on day 18.

Worked example: account for the halo

The direct blur reads a 5 by 5 window from global memory for every output, so it issues 25 global loads. Adjacent outputs request many of the same pixels.

A 32 by 32 block instead loads one 36 by 36 apron into shared memory. It performs 36 * 36 = 1296 global loads for 32 * 32 = 1024 outputs, or 1296 / 1024 = 1.265625 loads per output.

The capstone budget is 1.6 loads per output. A static_assert enforces the geometry, so a tile shape that exceeds the budget fails the build. The 1.27 design leaves room for another tile size while rejecting 25 direct loads by a factor of more than fifteen.

The apron load uses signed coordinates rather than size_t because its top-left corner can be two pixels above or left of the image. It guards each load and leaves the __syncthreads() outside the guard, so all live threads reach the barrier.

    // 1024 threads fill 1296 apron slots, so 272 of the loads are second
    // helpings, and they fall on the 240 threads whose tx or ty is under 4.
    // The loops step by kTileDim rather than dividing by the pitch, because a
    // runtime integer division inside the load would cost more than the load,
    // and would land differently on the padded row than on the unpadded one.
    for (int sy = ty; sy < kBlurShared; sy += kTileDim) {
        for (int sx = tx; sx < kBlurShared; sx += kTileDim) {
            const long long gy = y0 + sy;
            const long long gx = x0 + sx;
            const bool inside = gy >= 0 && gy < static_cast<long long>(rows) &&
                                gx >= 0 && gx < static_cast<long long>(cols);
            tile[sy * pitch + sx] = inside ? in[static_cast<size_t>(gy) * cols +
                                                static_cast<size_t>(gx)]
                                           : 0.0f;
        }
    }
    __syncthreads();

The apron starts two floats before an aligned image row, so its load can cross an extra 32-byte sector. L1 and L2 caches also recover some overlap in the direct kernel, so 25 divided by 1.27 is not a runtime prediction.

Worked example: keep the intermediate on chip

The one-kernel separable blur writes the row-pass result to a second shared array. A second barrier makes those values visible before the column pass.

    // The row pass runs over all 36 apron rows, because the column pass below
    // needs two rows above and two below its own. It produces only 32 columns,
    // which is the width of the output tile.
    for (int sy = ty; sy < kBlurShared; sy += kTileDim) {
        float sum = 0.0f;
        for (int q = 0; q < kBlurDim; ++q) {
            sum += cGauss1d[q] * apron[sy][tx + q];
        }
        rowPass[sy][tx] = sum;
    }
    __syncthreads();

    const size_t row =
        blockIdx.y * static_cast<size_t>(kTileDim) + static_cast<size_t>(ty);
    const size_t col =
        blockIdx.x * static_cast<size_t>(kTileDim) + static_cast<size_t>(tx);
    if (row < rows && col < cols) {
        float sum = 0.0f;
        for (int p = 0; p < kBlurDim; ++p) {
            sum += cGauss1d[p] * rowPass[ty + p][tx];
        }
        out[row * cols + col] = sum;
    }

This kernel reads one input float and writes one output float, so it moves 8 bytes per output pixel. Splitting the passes across two launches adds a scratch-image write and read, raising the total to 16 bytes.

The arithmetic remains 10 multiply-adds in both forms. If memory traffic limits the blur, the byte counts predict that the two-kernel version will take about twice as long.

Worked example: remove an intermediate image

The staged edge path moves 13 bytes per pixel. Sobel reads and writes a float, then threshold reads that float and writes one byte.

The combined kernel moves 5 bytes per pixel because the magnitude stays in a register. The model predicts 13 / 5 = 2.6 for staged time over combined time if both paths are memory-bound.

__global__ void sobelThreshold(const float* __restrict__ in,
                               unsigned char* __restrict__ out, size_t rows,
                               size_t cols, float t) {
    __shared__ float tile[kSobelShared][kSobelShared];
    loadSobelApron(in, tile, rows, cols);
    __syncthreads();

    const size_t row = blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.y;
    const size_t col = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.x;
    if (row < rows && col < cols) {
        const float mag = sobelMagnitudeAt(tile, static_cast<int>(threadIdx.y),
                                           static_cast<int>(threadIdx.x));
        out[row * cols + col] = (mag > t) ? 255 : 0;
    }
}

Day 48 develops this combination into a general method. Here it supplies one prediction that the timing table can test.

Check before you run

Answer these without looking at the results:

1. Does pitch 37 remove a bank conflict from the row-major convolution access?
2. Can the program return early for an out-of-range output before the apron barrier?
3. If the combined edge kernel measures only 1.9 times faster instead of 2.6, does that make the byte count wrong?

Check answers

1. No. Within a warp, threadIdx.x changes by one, so the addresses already land in 32 distinct banks for either pitch.
2. No. The out-of-range thread may still own apron cells that an in-range neighbour reads, so guard the load and store but not the barrier.
3. No. The 2.6 ratio assumes traffic is the only limit; arithmetic, instruction issue, cache effects, and launch cost can reduce the measured ratio.

How the program checks the work

The full program is in code/day20-convolution/convolution.cu. Build it from that directory with nvcc -std=c++17 -O3 -arch=sm_75 -o convolution convolution.cu, then run ./convolution.

The program needs about 270 MiB of device memory and 400 MiB of host memory. Most of its runtime is the single-threaded CPU reference over 16,777,216 pixels.

Every timing row prints milliseconds, effective memory bandwidth, bytes moved, and a ratio to a copy over the same buffers. It uses CUDA events, three warm-ups, and the mean of ten runs, with copies outside the timed region.

One dynamic-shared-memory kernel tests pitch 36 and pitch 37 with the same instruction mix. The timed case is 4096 by 4096, while the correctness case is 613 rows by 611 columns so both grid axes end in partial tiles.

The width 611 is 13 * 47, and 613 is prime, so neither dimension divides by 32. The grid therefore rounds up on both axes and forces every bounds guard to run.

The run also checks a 256 by 256 impulse. Its peak must stay in place, its values must sum to 1.0, and all 25 taps must match the filter.

The blur uses relative tolerance 1e-5 and absolute tolerance 1e-6. Sobel uses absolute tolerance 1e-4 because gx and gy subtract nearly equal sums, leaving an answer near zero while input rounding remains.

Threshold byte differences are allowed only when the reference magnitude is within 1e-4 of the threshold. The program prints how many pixels used that allowance, so the tolerance cannot hide a broad error.

The output files are day20_input.pgm, day20_blur.pgm, and day20_edges.pgm. Open day20_edges.pgm in an image viewer; it is the visual deliverable.

Hardware note

The 32 by 32 output tile uses 1024 threads, the maximum block size on every compute capability this course targets. On a T4 this gives 32 resident warps and full occupancy, but only one resident block per SM.

A 32 by 8 block with four output rows per thread can keep the same occupancy while allowing more blocks to reside on an SM. That experiment is described in the code README and is not part of this gate.

The source enforces the 1.6 halo budget with tile geometry, not a runtime traffic audit. It also checks that each shared layout fits the T4's 48 KiB per-block limit and that the test dimensions exercise the intended grid edges.

Record the evidence

Now compare your predictions with the recorded T4 run. These measurements came from the exact shipped program with CUDA 12.6.

GPU: Tesla T4 (compute capability 7.5)
timed case 4096 x 4096, 64.0 MiB per float buffer
case odd      613 x 611, grid (20, 20), block (32, 32)
  every kernel matches the CPU reference at 374543 pixels, 0 pixels forgiven inside the threshold band
case impulse  256 x 256, peak in place, sums to 1.000000, all 25 taps equal the filter
wrote day20_blur.pgm and day20_edges.pgm, 2.43% of the edge image is white
case perf     4096 x 4096, grid (128, 128), block (32, 32)
  every kernel matches the CPU reference at 16777216 pixels, 0 pixels forgiven inside the threshold band

timed on 4096 x 4096, mean of 10 runs after 3 warm-ups, copies not included
stage                                     ms      GB/s         bytes   x copy
copy (the baseline)                    0.566     237.0     134217728     1.00
blur, naive                            3.242      41.4     134217728     0.17
blur, tiled, pitch 36                  3.033      44.3     134217728     0.19
blur, tiled, pitch 37 (padded)         3.028      44.3     134217728     0.19
blur, separable, one kernel            2.619      51.2     134217728     0.22
sobel magnitude                        2.221      60.4     134217728     0.25
threshold                              0.344     244.0      83886080     1.03
sobel + threshold, two kernels         2.565      85.0     218103808     0.36
sobel + threshold, fused               1.323      63.4      83886080     0.27
pipeline, blur + fused                 2.667      81.8     218103808     0.34
CPU, one thread, whole pipeline      978.101       n/a           n/a      n/a

CPU / GPU pipeline speedup: 366.7x
pipeline at 34.5% of this card's copy bandwidth
padding the shared tile changed the blur by 0.17%

Pitch 37 changed the blur by 0.17 percent, consistent with the prediction that the row access had no bank conflict. This differs from day 15's transpose, where changing pitch 32 to 33 improved copy-relative bandwidth from 57.1 to 82.5 percent, a 45 percent gain, because each warp read down a column.

The common declaration __shared__ float tile[36][36 + 1]; applies padding without first proving a conflict. Here the bank is (byteAddress / 4) % 32, consecutive lanes already reach distinct banks, and the + 1 costs 144 bytes of shared memory per block.

For a column read, the conflict degree is gcd(stride, 32). A pitch of 32 has gcd(32, 32) = 32, which is the 32-way transpose case that pitch 33 fixes.

The staged edge path took 2.565 ms and moved about 218 MB, while the combined path took 1.323 ms and moved about 84 MB. That is a 1.9x speedup rather than the traffic ceiling of 2.6x, which shows that traffic was not the only limit.

The one-kernel separable blur took 2.619 ms against 3.033 ms for the direct 2D tiled blur. Both moved 8 bytes per pixel, so reducing 25 multiply-adds to 10 cut the arithmetic by 2.5x and helped even though neither blur exceeded 22 percent of copy bandwidth.

The blur's arithmetic intensity is low, but the copy comparison alone cannot identify the remaining limit. Module 5 uses profiler counters to separate arithmetic, halo loads, and instruction limits.

The whole GPU pipeline took 2.667 ms against 978.101 ms for one CPU thread, a 366.7x ratio. This ratio is reported rather than graded because CPU and GPU models change the denominator.

The original run used driver 595.84 and CUDA 12.6. CUDA 13.0 V13.0.88 with driver 580.173.02 reproduced the behavior on 2026-09-02, reported a 0.00 percent pitch difference, and returned exit code 0; both transcripts remain in evidence/.

Faded practice

First, recompute the complete apron example for a 16 by 16 output tile with radius 2. The shared tile is 20 by 20, so the load count is 400 for 256 outputs, or 1.5625 loads per output, which still fits the 1.6 budget.

Next, complete this traffic account before opening the answer:

Variant	Input read	Intermediate write	Intermediate read	Output write	Total
one-kernel separable	4	0	0	4	8 bytes/pixel
two-kernel separable	4	fill in	fill in	4	calculate

Traffic answer

The two-kernel version writes and reads one float intermediate, so the missing entries are 4, 4, and 16 bytes per pixel. With equal arithmetic, a purely memory-bound run would take about twice as long.

Independent gate

Implement blurRows and blurCols with a full-image scratch buffer. Do not change the Gaussian coefficients, border rule, CPU reference, tolerances, warm-up count, timed-run count, or existing kernels.

Your work passes only when all of these statements are true:

• the build command above succeeds and ./convolution returns exit code 0;
• your two-kernel result matches the CPU blur at all 374,543 pixels in the 613-row by 611-column case and all 16,777,216 pixels in the 4096 by 4096 case;
• a failure prints the first mismatching row, column, value found, and value expected;
• your timing row reports 268435456 bytes for the 4096 by 4096 case, which is 16 bytes per pixel;
• you write the predicted ratio before timing, then report the measured ratio two-kernel ms / one-kernel ms and explain any gap from 2.0 using the byte model, launch cost, or cache behavior.

An exit code alone is not proof of correctness unless you add the new variant to both CPU-reference checks. A timing row without the byte count also fails because milliseconds from variants with different traffic do not compare the same work.

Hint 1

The row pass reads one float image and writes one float image. The column pass does the same, and the second pass must filter the other axis.

Hint 2

Size the scratch allocation for the full intermediate image. Reuse the program's grid, bounds checks, error checks, CPU comparison, and event-timing helper.

Common failures

If pitch 37 appears much faster, check that both rows launch the same kernel with only the pitch argument changed. For this row-major access, lanes read consecutive words and already use 32 different memory banks.

Do not return before the apron barrier described on day 14. A thread whose output is out of range can still own a shared cell needed by a neighbouring output, so guard the apron load and output store separately.

In particular, do not place if (row >= rows) return; above the barrier. compute-sanitizer --tool synccheck checks barrier use, but the odd-size CPU comparison is the gate that catches a missing apron value.

A kernel that passes 1024 by 1024 can still fail at 613 by 611 because the first size never executes a partial tile. Day 7 uses 611 by 397 for the same reason.

Do not compare thresholded bytes without the threshold band, and do not grade the Sobel magnitude with only a relative tolerance. Near zero, cancellation makes the absolute input error larger than a relative limit based on the small output.

Constant memory helps here because every lane reads the same filter tap for a given loop iteration. Indexing the coefficient by threadIdx.x would make the warp request different addresses and serialize the constant-memory accesses.

Run on a smaller machine

There is no Compiler Explorer embed because the CPU reference exceeds its 20-second cap before the first kernel launches. If that is your only option, use 512 by 512, skip the timing gate, and keep the correctness checks.

Colab's free T4 or any local CUDA GPU can run the full case. Without a GPU, use the smaller Compiler Explorer path above and skip only the timing gate.

Sources

• CUDA C++ Best Practices Guide: coalesced global access and shared-memory banks, checked 2026-08-30.
• CUDA C++ Programming Guide: writing kernels with shared memory and barriers, checked 2026-08-30.
• NVIDIA cuda-samples: cpp/2_Concepts_and_Techniques/convolutionSeparable, checked 2026-08-30.
• Programming Massively Parallel Processors, fourth edition, chapter 7.
• CUDA 12.6 T4 transcript and CUDA 13.0 T4 transcript.

Transfer

For another square filter, use radius r and output width T: the apron is (T + 2r)², and its global-load cost is (T + 2r)² / T² per output. Then ask whether the filter is separable and whether an intermediate image can stay on chip.

Day 21 begins module 3 by measuring warp behavior directly. Day 60 reuses this blur, Sobel, and threshold chain across 900 frames with overlapped copies.