CUDA Sobel edge detection: from naive kernel to profiled pipeline

Prerequisites

This capstone case study pulls together the whole series: convolution in CUDA for the shared-memory halo-tiling pattern, advanced memory patterns for pinned memory, advanced stream patterns for overlapping transfers with compute, and performance case studies for the benchmarking and roofline methodology. If any of those feel unfamiliar, this is a good moment to go back; this article assumes you can read a shared-memory tiled kernel and a multi-stream pipeline without a line-by-line refresher.

Why Sobel is a good capstone

Sobel edge detection is small enough to implement from scratch in an afternoon and rich enough to exercise almost every optimization technique this series has covered: it is a stencil (each output depends on a small neighborhood of inputs, like the convolution case study), it has a border-handling problem (unlike a dense matrix multiply, edge pixels are a special case), its filter weights are tiny and read identically by every thread (a textbook use for constant memory), and it is naturally batchable across many independent images (a textbook use for streams and pinned double buffering). Building it end to end, from a naive kernel through a profiled, streamed pipeline, is a condensed tour of the entire series.

The operator

Sobel edge detection estimates the image gradient at each pixel using two small, fixed 3x3 convolution kernels: one that approximates the horizontal derivative (Gx) and one that approximates the vertical derivative (Gy).

Gx = | -1  0  +1 |        Gy = | -1  -2  -1 |
     | -2  0  +2 |             |  0   0   0 |
     | -1  0  +1 |             | +1  +2  +1 |

For a grayscale input image I, the gradient magnitude at pixel (x, y) is:

gx = sum over the 3x3 neighborhood of I * Gx
gy = sum over the 3x3 neighborhood of I * Gy
magnitude(x, y) = sqrt(gx^2 + gy^2)

sqrt(gx^2 + gy^2) is the mathematically correct magnitude; |gx| + |gy| is a common cheaper approximation that avoids a square root per pixel at the cost of slightly overestimating the true magnitude on diagonal edges. This article uses the exact sqrtf form throughout so the CPU oracle and every GPU variant compute the identical formula; swapping in the approximation is a one-line change you can make once correctness is established.

Grayscale conversion

Sobel operates on a single-channel image. If your input is RGB, convert first using a standard luminance-weighted formula (these particular weights come from the ITU-R BT.601 standard and approximate human luminance perception):

__device__ __forceinline__ float rgb_to_gray(unsigned char r, unsigned char g, unsigned char b) {
    return 0.299f * r + 0.587f * g + 0.114f * b;
}

This case study treats grayscale conversion as a separate, trivially parallel first pass (one thread per pixel, no neighborhood dependency) and focuses its optimization effort on the Sobel stencil itself, which is where the interesting memory-access patterns live.

Border policy

The 3x3 neighborhood needed by a pixel on the image’s edge extends outside the image. This article uses clamp-to-edge: any neighborhood coordinate outside [0, width) x [0, height) is clamped to the nearest valid coordinate, which is equivalent to treating the border as replicated. This is a deliberate choice, not the only correct one; zero-padding (treating out-of-bounds reads as 0) and mirrored padding are both used elsewhere in image processing, and would change the exact numeric output near the border. Whichever policy you pick, the CPU oracle and every GPU kernel variant must implement the same one, or your correctness check will fail at the borders even when the interior of the image matches perfectly.

One thread, one output pixel

Every kernel variant in this article maps one CUDA thread to one output pixel, using a 2D grid of 2D blocks, exactly the mapping introduced in CUDA thread hierarchy:

graph TD
  IMG["Input image
width x height"] --> BLK["2D grid of 2D blocks
e.g. 16x16 threads per block"]
  BLK --> PX["One thread per output pixel
reads its own 3x3 neighborhood"]
  PX --> OUT["Output edge magnitude image"]
  style PX fill:#e8744f,color:#fff

The CPU oracle

Before writing any GPU code, write a plain, obviously-correct CPU implementation. Its only job is to be a ground truth for validating every GPU variant against, so it favors clarity over speed:

void sobel_cpu(const float* gray, float* out, int width, int height) {
    static const int gx_kernel[3][3] = {{-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1}};
    static const int gy_kernel[3][3] = {{-1, -2, -1}, {0, 0, 0}, {1, 2, 1}};

    for (int y = 0; y < height; y++) {
        for (int x = 0; x < width; x++) {
            float gx = 0.0f, gy = 0.0f;
            for (int dy = -1; dy <= 1; dy++) {
                for (int dx = -1; dx <= 1; dx++) {
                    int sx = std::min(std::max(x + dx, 0), width - 1);   // clamp-to-edge
                    int sy = std::min(std::max(y + dy, 0), height - 1);
                    float pixel = gray[sy * width + sx];
                    gx += pixel * gx_kernel[dy + 1][dx + 1];
                    gy += pixel * gy_kernel[dy + 1][dx + 1];
                }
            }
            out[y * width + x] = sqrtf(gx * gx + gy * gy);
        }
    }
}

Every GPU kernel below is checked against this function with a tolerance-based comparison (never exact equality; see floating-point performance for why), because the summation order and rounding of gx/gy on the GPU need not match the CPU’s loop order bit-for-bit even when both are correct.

bool compare_to_oracle(const float* cpu, const float* gpu, int n, float atol = 1e-3f) {
    float max_abs_diff = 0.0f;
    for (int i = 0; i < n; i++) {
        max_abs_diff = std::max(max_abs_diff, std::fabs(cpu[i] - gpu[i]));
    }
    printf("Max abs diff vs CPU oracle: %e (tolerance %e)\n", max_abs_diff, atol);
    return max_abs_diff <= atol;
}

Variant 1: the naive kernel

The most direct translation of the operator into CUDA: each thread reads its own 3x3 neighborhood straight from global memory, with the filter weights hardcoded in the kernel body.

__global__ void sobel_naive(const float* gray, float* out, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x >= width || y >= height) return;   // bounds check for non-multiple grid sizes

    static const int gx_kernel[3][3] = {{-1, 0, 1}, {-2, 0, 2}, {-1, 0, 1}};
    static const int gy_kernel[3][3] = {{-1, -2, -1}, {0, 0, 0}, {1, 2, 1}};

    float gx = 0.0f, gy = 0.0f;
    #pragma unroll
    for (int dy = -1; dy <= 1; dy++) {
        #pragma unroll
        for (int dx = -1; dx <= 1; dx++) {
            int sx = min(max(x + dx, 0), width - 1);
            int sy = min(max(y + dy, 0), height - 1);
            float pixel = gray[sy * width + sx];
            gx += pixel * gx_kernel[dy + 1][dx + 1];
            gy += pixel * gy_kernel[dy + 1][dx + 1];
        }
    }
    out[y * width + x] = sqrtf(gx * gx + gy * gy);
}

Every thread issues 9 independent global memory loads. Neighboring threads (adjacent in x) read heavily overlapping neighborhoods: for a horizontal run of threads, each pixel is re-read from global memory by up to 3 different threads (once for each row of the 3x3 window it participates in). This redundant traffic is exactly the same problem the convolution case study identified, and it is the reason this naive version is not the end of the story.

Variant 2: constant memory for the filter weights

The gx_kernel/gy_kernel arrays above are tiny (9 ints each) and read identically by every single thread in the grid. This is the textbook use case for __constant__ memory: the constant cache broadcasts one read to an entire warp instead of each thread pulling the same value through the normal memory path.

__constant__ int c_gx[9];   // flattened 3x3, row-major
__constant__ int c_gy[9];

void upload_sobel_filters() {
    int gx_kernel[9] = {-1, 0, 1, -2, 0, 2, -1, 0, 1};
    int gy_kernel[9] = {-1, -2, -1, 0, 0, 0, 1, 2, 1};
    cudaMemcpyToSymbol(c_gx, gx_kernel, sizeof(gx_kernel));
    cudaMemcpyToSymbol(c_gy, gy_kernel, sizeof(gy_kernel));
}

__global__ void sobel_constant(const float* gray, float* out, int width, int height) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;
    int y = blockIdx.y * blockDim.y + threadIdx.y;
    if (x >= width || y >= height) return;

    float gx = 0.0f, gy = 0.0f;
    int k = 0;
    #pragma unroll
    for (int dy = -1; dy <= 1; dy++) {
        #pragma unroll
        for (int dx = -1; dx <= 1; dx++, k++) {
            int sx = min(max(x + dx, 0), width - 1);
            int sy = min(max(y + dy, 0), height - 1);
            float pixel = gray[sy * width + sx];
            gx += pixel * c_gx[k];
            gy += pixel * c_gy[k];
        }
    }
    out[y * width + x] = sqrtf(gx * gx + gy * gy);
}

This removes the filter weights from each thread’s register/immediate footprint and lets the compiler and hardware serve them from the constant cache, but it does not address the redundant global memory reads of the image data identified above; that requires shared memory.

Variant 3: shared memory tile with a one-pixel halo

Each block cooperatively loads a tile of the input image, padded by a one-pixel border (the “halo”) on every side, into shared memory once. Every thread in the block then computes its 3x3 neighborhood entirely from shared memory, with zero further global memory traffic per output pixel.

graph TD
  subgraph SharedTile["Shared memory tile (TILE+2) x (TILE+2)"]
  H["Halo row/col
loaded by border threads"]
  C["Core TILE x TILE
loaded by all threads"]
  end
  GMEM["Global memory
image data"] -->|"cooperative load"| SharedTile
  SharedTile -->|"__syncthreads()"| COMPUTE["Every thread computes
its own 3x3 neighborhood
from shared memory only"]
  style H fill:#9b59b6,color:#fff
  style C fill:#4a90d9,color:#fff

#define TILE 16
#define RADIUS 1   // 3x3 filter -> radius 1 halo on each side

__global__ void sobel_shared(const float* gray, float* out, int width, int height) {
    __shared__ float tile[TILE + 2 * RADIUS][TILE + 2 * RADIUS];

    int tx = threadIdx.x, ty = threadIdx.y;
    int x = blockIdx.x * TILE + tx;   // this thread's output pixel
    int y = blockIdx.y * TILE + ty;

    // Cooperative load: every thread loads one core element, and threads
    // near the block's edges load one additional halo element each.
    // Clamp-to-edge is applied at load time so the tile already contains
    // border-replicated values; the compute step below needs no further
    // bounds handling for the neighborhood itself.
    for (int ly = ty; ly < TILE + 2 * RADIUS; ly += TILE) {
        for (int lx = tx; lx < TILE + 2 * RADIUS; lx += TILE) {
            int gx_coord = blockIdx.x * TILE + lx - RADIUS;
            int gy_coord = blockIdx.y * TILE + ly - RADIUS;
            int sx = min(max(gx_coord, 0), width - 1);
            int sy = min(max(gy_coord, 0), height - 1);
            tile[ly][lx] = gray[sy * width + sx];
        }
    }
    __syncthreads();   // the whole tile, including halo, must be loaded before anyone reads it

    if (x >= width || y >= height) return;   // guard for images not a multiple of TILE

    float gx = 0.0f, gy = 0.0f;
    int k = 0;
    #pragma unroll
    for (int dy = 0; dy < 3; dy++) {
        #pragma unroll
        for (int dx = 0; dx < 3; dx++, k++) {
            float pixel = tile[ty + dy][tx + dx];   // shared memory only, no global read here
            gx += pixel * c_gx[k];
            gy += pixel * c_gy[k];
        }
    }
    out[y * width + x] = sqrtf(gx * gx + gy * gy);
}

Two details are easy to get wrong here and worth calling out explicitly:

• The cooperative load loop, not a single if per halo side. Writing separate if (tx < RADIUS) / if (tx >= TILE - RADIUS) branches for each halo edge (as a first attempt often does) works but is harder to get right for corner cells, which need both a horizontal and vertical halo offset simultaneously. The strided for loop above handles core, edge, and corner halo cells uniformly, at the cost of some threads doing more than one load iteration.
• The bounds check comes after the load, not before. Every thread, including ones whose output pixel (x, y) is outside the image (which happens whenever width/height are not exact multiples of TILE), still participates in the cooperative shared-memory load and the following __syncthreads(). Returning early before the load or the sync would leave __syncthreads() reached by only some threads in the block, which is undefined behavior: every thread in a block must reach the same __syncthreads() call.

Launch this kernel with a block size that matches TILE, and compute grid dimensions with ceiling division exactly as in every prior kernel in this series:

dim3 block(TILE, TILE);
dim3 grid((width + TILE - 1) / TILE, (height + TILE - 1) / TILE);
sobel_shared<<<grid, block>>>(d_gray, d_out, width, height);

Unrolling only after measurement

Both kernels above already carry #pragma unroll on their fixed 3x3 loops, which is a reasonable default for a compile-time-known trip count this small. But do not treat unrolling, or any other micro-optimization, as free: more aggressively unrolled code can increase register usage per thread, which (as the resource-allocation worked example in Inside a modern NVIDIA GPU showed) can reduce how many blocks fit resident on an SM at once, trading instruction-level parallelism within a thread for occupancy across threads. The correct process is: implement the clear version first, measure it with Nsight Compute, and only then try unrolling variations (#pragma unroll 1 to force no unrolling, #pragma unroll for full unrolling, or an explicit unroll factor) and re-measure. A pragma that reads well in a tutorial is not a substitute for a profiler telling you it actually helped on your kernel, on your data size, on your GPU.

Batching across images: streams and pinned double buffering

A single image’s pipeline is inherently sequential: you cannot start the Sobel kernel before its input has finished transferring to the device, and you cannot start the device-to-host copy of the result before the kernel has finished writing it. Declaring a stream for one image does not create overlap by itself, because there is nothing independent to overlap; H2D copy, kernel, and D2H copy for a single image form a strict dependency chain regardless of which stream they run on.

Overlap becomes possible only when you have multiple independent images in flight simultaneously, so that one image’s kernel execution can run concurrently with a different image’s transfer. This is the classic double-buffered producer/consumer pipeline, and it requires pinned host memory (from advanced memory patterns) so the transfers are actually asynchronous in the first place:

gantt
  dateFormat X
  axisFormat %L

  section Stream A
  H2D img0      :a1, 0, 10
  Kernel img0   :a2, after a1, 15
  D2H img0      :a3, after a2, 10
  H2D img2      :a4, after a3, 10

  section Stream B
  H2D img1      :b1, 5, 10
  Kernel img1   :b2, after a2, 15
  D2H img1      :b3, after b2, 10

The diagram is schematic, not a measured timeline; the important shape is that stream B’s H2D img1 overlaps with stream A’s Kernel img0, and stream B’s Kernel img1 overlaps with stream A’s D2H img0, because those pairs use different hardware engines (the copy engine and the SMs respectively) and belong to independent streams with no artificial dependency between them.

const int NUM_STREAMS = 2;
cudaStream_t streams[NUM_STREAMS];
for (int i = 0; i < NUM_STREAMS; i++) cudaStreamCreate(&streams[i]);

// Pinned host buffers, one pair per stream slot, reused across the batch.
float* h_pinned_in[NUM_STREAMS];
float* h_pinned_out[NUM_STREAMS];
float* d_in[NUM_STREAMS];
float* d_out[NUM_STREAMS];
for (int i = 0; i < NUM_STREAMS; i++) {
    cudaMallocHost(&h_pinned_in[i], image_bytes);
    cudaMallocHost(&h_pinned_out[i], image_bytes);
    cudaMalloc(&d_in[i], image_bytes);
    cudaMalloc(&d_out[i], image_bytes);
}

for (int img = 0; img < num_images; img++) {
    int slot = img % NUM_STREAMS;
    cudaStream_t s = streams[slot];

    // A host write is not ordered by a CUDA stream. Wait before reusing this
    // slot so DMA cannot still be reading its pinned input buffer.
    if (img >= NUM_STREAMS) {
        cudaStreamSynchronize(s);
        consume_result(h_pinned_out[slot], img - NUM_STREAMS);
    }

    load_image_into(h_pinned_in[slot], img);   // host-side decode/read, not shown

    cudaMemcpyAsync(d_in[slot], h_pinned_in[slot], image_bytes,
                     cudaMemcpyHostToDevice, s);
    sobel_shared<<<grid, block, 0, s>>>(d_in[slot], d_out[slot], width, height);
    cudaMemcpyAsync(h_pinned_out[slot], d_out[slot], image_bytes,
                     cudaMemcpyDeviceToHost, s);
}

// Drain and consume the final image queued in each active slot.
int first_pending = num_images > NUM_STREAMS ? num_images - NUM_STREAMS : 0;
for (int img = first_pending; img < num_images; img++) {
    int slot = img % NUM_STREAMS;
    cudaStreamSynchronize(streams[slot]);
    consume_result(h_pinned_out[slot], img);
}

Using img % NUM_STREAMS means image 2 eventually reuses slot 0. Operations enqueued on the GPU are ordered within that stream, but load_image_into is a regular host write and is not ordered by the stream. The host must therefore wait for the slot before overwriting its pinned input, and it must consume the prior output before the next D2H copy reuses that output buffer. The example uses cudaStreamSynchronize at reuse because it is easy to audit. A production pipeline can record a completion event per slot and wait or poll that event instead, allowing the CPU to decode into a deeper pool of host buffers while two GPU streams remain in flight.

Benchmark methodology

Reuse the warmup, CUDA-event, and repeated-sample methodology from performance case studies directly; this section only adds what is specific to a stencil pipeline like this one.

• Warm-up. Run the full pipeline (transfer, kernel, transfer back) several times before timing, exactly as before; the first iteration also pays for context setup and any JIT compilation, on top of the usual driver warm-up cost.
• Kernel-only timing with CUDA events. Bracket only the kernel launch with cudaEventRecord/cudaEventElapsedTime to isolate compute time from transfer time. This is what you compare across the naive, constant-memory, and shared-memory variants, since they differ only in kernel implementation.
• Wall-clock end-to-end timing. Separately, time the entire per-image pipeline (host-side image load, H2D, kernel, D2H, any host-side postprocessing) with a host timer (std::chrono::steady_clock) around the whole batch loop. This is the number that answers “how long does processing N images actually take,” including transfer and any streaming overlap, and it is not the same question the kernel-only CUDA event timing answers.
• Repeated samples. Run each configuration many times (20-50, per the earlier article’s guidance) and report median and interquartile range, not a single run or a mean, for exactly the reasons given there: run-to-run variance from clock boosting, OS scheduling, and thermal state is real and a single sample can mislead you.
• Size sweep. Repeat every measurement across a range of image resolutions (for example, a small thumbnail size, a common “HD” size, and a large size). The relative ranking of the naive, constant-memory, and shared-memory kernels can change with image size, because the constant-memory cache and the ratio of halo overhead to tile size behave differently at different scales.
• Correctness at every step. Run the CPU oracle comparison from the beginning of this article against every kernel variant, at every size in the sweep, every time you change the kernel. A kernel that is faster but no longer correct is not an optimization; it is a bug. Never trust a timing number from a variant you have not also validated for that same run’s input.

Profiling with Nsight

This article deliberately does not re-derive how to use Nsight Systems or Nsight Compute; that workflow is covered in full in debugging and profiling CUDA programs and performance case studies. What is specific to this pipeline is what to look for:

• Nsight Systems, run across the whole batched pipeline (nsys profile ./sobel_pipeline), should show the H2D/kernel/D2H bars for adjacent images overlapping on the timeline once double buffering is working, in roughly the pattern sketched in the Gantt diagram above. If the bars are fully sequential with no overlap, check first that you allocated pinned (not pageable) host memory, and second that consecutive images are genuinely alternating streams rather than accidentally sharing one.
• Nsight Compute, run on a single kernel invocation (ncu --set full ./sobel_pipeline), places each kernel variant on the roofline chart introduced in floating-point performance. Sobel’s arithmetic intensity is low (9 multiply-adds and one sqrtf per pixel, against at least one 4-byte read and one 4-byte write per pixel, before accounting for the shared-memory version’s reuse), so expect these kernels to sit well to the memory-bound side of the roofline boundary; the shared-memory variant’s entire purpose is to move the achieved memory throughput closer to the device’s peak memory bandwidth by eliminating redundant reads, not to change which side of the roofline the kernel sits on.
• Nsight Compute’s warp-state statistics are the right tool to confirm the naive kernel’s expected bottleneck (long-scoreboard/memory-dependency stalls from repeated global loads of overlapping neighborhoods) and to confirm the shared-memory kernel actually reduced that specific stall reason, rather than assuming it did because the code looks like it should.

Result table template

This article does not publish specific timing numbers, because they depend entirely on your GPU, image sizes, and driver/toolkit version, and a number copied from a tutorial into your own performance report would be a fabricated benchmark claim, not a measurement. Use the harness and methodology above to fill in a table shaped like this for your own hardware:

Fill every cell from your own harness output before drawing any conclusion about which variant “wins,” and expect the ranking to depend on image size, exactly as the size-sweep guidance above anticipates.

An analogy: stencils and thermal diffusion

The shared-memory halo-tiling pattern in this article is not specific to image processing. Any stencil computation, where each output grid point depends on a small, fixed neighborhood of its inputs, has the identical structure. The most common example outside image processing is a finite-difference solver for the heat (thermal diffusion) equation, where each grid cell’s next-timestep temperature is computed from its current temperature and its immediate neighbors:

T_new[i][j] = T[i][j] + alpha * (T[i-1][j] + T[i+1][j] + T[i][j-1] + T[i][j+1] - 4*T[i][j])

This is a 5-point stencil instead of Sobel’s 3x3, 9-point stencil, but the CUDA implementation strategy is identical: tile the domain across thread blocks, cooperatively load each tile plus a halo of neighboring cells into shared memory, synchronize, then compute every interior point from shared memory alone. The same non-multiple-of-tile-size bounds handling, the same “load before you check bounds, or __syncthreads() diverges” hazard, and the same benchmarking questions (kernel-only vs. end-to-end, size sweep, roofline placement) all carry over directly. If you understand why the Sobel shared-memory kernel in this article is built the way it is, you already understand the core of how a GPU-accelerated PDE solver is structured, which is one of the reasons stencil optimization is treated as its own well-studied subfield of high-performance computing rather than an image-processing-specific trick.

In practice

Build this pipeline incrementally and validate at every step: CPU oracle first, naive GPU kernel checked against it, constant memory checked against the naive kernel’s output (they must match exactly, since only the memory space of the filter weights changed), shared memory checked against both, and the streamed batch pipeline checked image-by-image against the single-image shared-memory kernel’s output. Skipping intermediate correctness checks and only validating the final, most-optimized version makes it much harder to tell whether a discrepancy came from the halo-loading logic, the border policy, or the streaming buffer reuse. Profile before assuming an optimization helped; a kernel that “should” be faster because it uses shared memory can still be measured to confirm it, and Nsight Compute’s stall-reason breakdown will tell you specifically whether the change moved the needle on the bottleneck it was meant to address.

What comes next

You have now built a complete, profiled CUDA pipeline from a mathematical operator down to a streamed, pinned-memory, multi-image batch, using the correctness discipline, memory hierarchy, streaming, and profiling tools developed across this entire series. The final article, Where to go from here, maps the broader CUDA ecosystem beyond hand-written kernels: the libraries, compilers, and frameworks you would reach for once a pattern like this one is proven correct and you want to decide whether to keep it as custom CUDA or replace it with a higher-level tool.