# CUDA for Engineers: 2D Grids and Interactive Graphics

- Launching 2D Computational Grids
- Live Display via Graphics Interop
- Application: Stability
- Summary
- Suggested Projects

In this chapter, we see that the CUDA model of parallelism extends readily to two dimensions (2D). We go through the basics of launching a 2D computational grid and create a skeleton kernel you can use to compute a 2D grid of values for functions of interest to you. We then specialize the kernel to create `dist_2d`, an app that computes the distance from a reference point in the plane to each member of a uniform 2D grid of points. By identifying the grid of points with pixels in an image, we compute data for an image whose shading is based on distance values.

Once we are generating image data, it is only natural to take advantage of CUDA’s **graphics interoperability** (or **graphics interop** for short) capability, which supports cooperation with standard graphics **application programming interfaces** (**APIs**) including **Direct3D** ^{1} and **OpenGL** ^{2}. We’ll use OpenGL, and maintaining our need-to-know approach, we’ll very quickly provide just the necessities of OpenGL to get your results on the screen at interactive speeds.

By the end of this chapter you will have run a `flashlight` app that interactively displays an image with shading based on distance from a reference point that you can move using mouse or keyboard input and a `stability` app that interactively displays the results of several hundred thousand numerical simulations of the dynamics of an oscillator. This experience should get you to the point where you are ready to start creating your own CUDA-powered interactive apps.

## Launching 2D Computational Grids

Here we expand on our earlier examples that involved a 1D array (points distributed regularly along a line segment) and move on to consider applications involving points regularly distributed on a portion of a 2D plane. While we will encounter other applications (e.g., simulating heat conduction) that fit this scenario, the most common (and likely most intuitive) example involves digital image processing. To take advantage of the intuitive connection, we will use image-processing terminology in presenting the concepts—all of which will transfer directly to other applications.

A digital raster image consists of a collection of picture elements or **pixels** arranged in a uniform 2D rectangular grid with each pixel having a quantized intensity value. To be concrete, let’s associate the width and height directions with the `x` and `y` coordinates, respectively, and say that our image is `W` pixels wide by `H` pixels high. If the quantized value stored in each pixel is simply a number, the data for an image matches exactly with the data for a matrix of size `W x H`.

As we move on from 1D to 2D problems in CUDA, we hope you will be pleasantly surprised by how few adjustments need to be made. In 1D, we specified integer values for block and grid sizes and computed an index `i` based on `blockDim.x`, `blockIdx.x`, and `threadIdx.x` according to the formula

`int i = blockIdx.x*blockDim.x + threadIdx.x;`

Here we reinterpret the expression on the right-hand side of the assignment as the specification of a new index `c` that keeps track of what column each pixel belongs to. (As we traverse a row of pixels from left to right, `c` increases from its minimum value `0` to its maximum value `W-1`.) We also introduce a second index `r` to keep track of row numbers (ranging from `0` to `H-1`). The row index is computed just as the column index is, but using the `.y` components (instead of the `.x` components), so the column and row indices are computed as follows:

```
int c = blockIdx.x*blockDim.x + threadIdx.x;
int r = blockIdx.y*blockDim.y + threadIdx.y;
```

To keep data storage and transfer simple, we will continue to store and transfer data in a “flat” 1D array, so we will have one more integer variable to index into the 1D array. We will continue to call that variable `i`, noting that `i` played this role in the 1D case, but in other places (including the CUDA Samples) you will see variables named `idx`, `flatIdx`, and `offset` indexing the 1D array. We place values in the 1D array in row major order—that is, by storing the data from row 0, followed by the data from row 1, and so on—so the index `i` in the 1D array is now computed as follows:

`int i = r*w + c;`

To describe the 2D computational grid that intuitively matches up with an image (or matrix or other regular 2D discretization), we specify block and grid sizes using `dim3` variables with two nontrivial components. Recall that an integer within the triple chevrons of a kernel call is treated as the `.x` component of a `dim3` variable with a default value of 1 for the unspecified `.y` and `.z` components. In the current 2D context, we specify nontrivial `.x` and `.y` components. The `.z` component of the `dim3`, which here has the default value 1, will come into play when we get to 3D grids in Chapter 7, “Interacting with 3D Data.”

Without further ado, let’s lay out the necessary syntax and get directly to parallel computation of pixel values with a 2D grid.

### Syntax for 2D Kernel Launch

The 2D kernel launch differs from the 1D launch only in terms of the execution configuration. Computing data for an image involves `W` columns and `H` rows, and we can organize the computation into 2D blocks with `TX` threads in the `x`-direction and `TY` threads in the `y`-direction. (You can choose to organize your 2D grid into 1D blocks, but you will run into limits on both maximum block dimension and total number of threads in a block. See the CUDA C Programming Guide ^{3} for details.)

We specify the 2D block size with a single statement:

`dim3 blockSize(TX, TY); // Equivalent to dim3 blockSize(TX, TY, 1);`

and then we compute the number of blocks (`bx` and `by`) needed in each direction exactly as in the 1D case.

```
int bx = (W + blockSize.x - 1)/blockSize.x ;
int by = (H + blockSize.y – 1)/blockSize.y ;
```

The syntax for specifying the grid size (in blocks) is

`dim3 gridSize = dim3 (bx, by);`

With those few details in hand, we are ready to launch:

kernelName<<<gridSize, blockSize>>>(args)

### Defining 2D Kernels

The prototype or declaration of a kernel to be launched on a 2D grid will look exactly as before: it starts with the qualifier `__global__` followed by return type `void` and a legal name, such as `kernel2D`, and ends with a comma-separated list of typed arguments (which better include a pointer to a device array `d_out` where the computed image data will be stored, along with the width and height of the image and any other required inputs). The `kernel2D` function begins by computing the row, column, and flat indices and testing that the row and column indices have values corresponding to a pixel within the image. All that is left is computing the value for the pixel.

Putting the pieces together, the structure of a typical 2D kernel is given in Listing 4.1.

#### Listing 4.1 “Skeleton” listing for a kernel to be launched on a 2D grid. Replace `INSERT_CODE_HERE` with your code for computing the output value.

```
1 __global__
2 void kernel2D(float *d_out, int w, int h, ... )
3 {
4 // Compute column and row indices.
5 const int c = blockIdx.x * blockDim.x + threadIdx.x;
6 const int r = blockIdx.y * blockDim.y + threadIdx.y;
7 const int i = r * w + c; // 1D flat index
8
9 // Check if within image bounds.
10 if ((c >= w) || (r >= h))
11 return;
12
13 d_out[i] = INSERT_CODE_HERE; // Compute/store pixel in device array.
14 }
```

One detail worth dealing with at this point is a common data type for images. The quantized value stored for each pixel is of type `uchar4`, which is a vector type storing four unsigned character values (each of which occupies 1 byte of storage). For practical purposes, you can think of the four components of the `uchar4` (designated as usual by suffixes `.x`, `.y`, `.z`, and `.w`) as specifying integer values ranging from 0 to 255 for the red, green, blue, and alpha (opacity) display channels. This format for describing pixel values in an image is often abbreviated as **RGBA**.

Putting the pieces together, the structure of a typical 2D kernel for computing an image is given in Listing 4.2.

#### Listing 4.2 “Skeleton” listing for computing data for an image. `RED_FORMULA`, `GREEN_FORMULA`, and `BLUE_FORMULA` should be replaced with your code for computing desired values between 0 and 255 for each color channel.

```
1 __global__
2 void kernel2D(uchar4 *d_output, int w, int h, ... )
3 {
4 // Compute column and row indices.
5 int c = blockIdx.x*blockDim.x + threadIdx.x;
6 int r = blockIdx.y*blockDim.y + threadIdx.y;
7 int i = r * w + c; // 1D flat index
8
9 // Check if within image bounds.
10 if ((r >= h) || (c >= w)) {
11 return;
12 }
13
14 d_output[i].x = RED_FORMULA; //Compute red
15 d_output[i].y = GREEN_ FORMULA; //Compute green
16 d_output[i].z = BLUE_ FORMULA; //Compute blue
17 d_output[i].w = 255; // Fully opaque
18 }
```

`dist_2d`

Let’s tie the general discussion of 2D grids together with our earlier examples involving distance apps by coding up an app that produces a 2D array of distances from a reference point, and then we’ll adapt the app to produce an array of data for an RGBA image. Listing 4.3 provides all the code for computing distances on a 2D grid.

#### Listing 4.3 Computing distances on a 2D grid

```
1 #define W 500
2 #define H 500
3 #define TX 32 // number of threads per block along x-axis
4 #define TY 32 // number of threads per block along y-axis
5
6 __global__
7 void distanceKernel(float *d_out, int w, int h, float2 pos)
8 {
9 const int c = blockIdx.x*blockDim.x + threadIdx.x;
10 const int r = blockIdx.y*blockDim.y + threadIdx.y;
11 const int i = r*w + c;
12 if ((c >= w) || (r >= h)) return;
13
14 // Compute the distance and set d_out[i]
15 d_out[i] = sqrtf((c - pos.x)*(c - pos.x) +
16 (r - pos.y)*(r - pos.y));
17 }
18
19 int main()
20 {
21 float *out = (float*)calloc(W*H, sizeof(float));
22 float *d_out; // pointer for device array
23 cudaMalloc(&d_out, W*H*sizeof(float));
24
25 const float2 pos = {0.0f, 0.0f}; // set reference position
26 const dim3 blockSize(TX, TY);
27 const int bx = (W + TX - 1)/TX;
28 const int by = (W + TY - 1)/TY;
29 const dim3 gridSize = dim3(bx, by);
30
31 distanceKernel<<<gridSize, blockSize>>>(d_out, W, H, pos);
32
33 // Copy results to host.
34 cudaMemcpy(out, d_out, W*H*sizeof(float), cudaMemcpyDeviceToHost);
35
36 cudaFree(d_out);
37 free(out);
38 return 0;
39 }
```

The kernel, lines 6–17, is exactly as in Listing 4.1 but with a result computed using the Pythagorean formula to compute the distance between the location {`c`, `r`} and a reference location `pos`. (Note that we have defined `pos` to have type `float2` so it can store both coordinates of the reference location {`pos.x`, `pos.y`}.) The rest of the listing, lines 19–39, gives the details of `main()` starting with declaration of an output array of appropriate size initialized to zero. Lines 22–23 declare a pointer to the device array `d_out` and allocate the memory with `cudaMalloc()`. Line 25 sets the reference position, and lines 26–29 set the kernel launch parameters: a 2D grid of `bx` × `by` blocks each having `TX` × `TY` threads. Line 31 launches the kernel to compute the distance values, which are copied back to `out` on the host side on line 34. Lines 36–37 free the allocated device and host memory, then `main()` returns zero to indicate completion.

Next we make a few minor changes to produce an app that computes an array of RGBA values corresponding to a distance image. The full code is provided in Listing 4.4.

#### Listing 4.4 Parallel computation of image data based on distance from a reference point in 2D

```
1 #define W 500
2 #define H 500
3 #define TX 32 // number of threads per block along x-axis
4 #define TY 32 // number of threads per block along y-axis
5
6 __device__
7 unsigned char clip(int n) { return n > 255 ? 255 : (n < 0 ? 0 : n); }
8
9 __global__
10 void distanceKernel(uchar4 *d_out, int w, int h, int2 pos)
11 {
12 const int c = blockIdx.x*blockDim.x + threadIdx.x;
13 const int r = blockIdx.y*blockDim.y + threadIdx.y;
14 const int i = r*w + c;
15 if ((c >= w) || (r >= h)) return;
16
17 // Compute the distance (in pixel spacings)
18 const int d = sqrtf((c - pos.x) * (c - pos.x) +
19 (r - pos.y) * (r - pos.y));
20 // Convert distance to intensity value on interval [0, 255]
21 const unsigned char intensity = clip(255 - d);
22
23 d_out[i].x = intensity; // red channel
24 d_out[i].y = intensity; // green channel
25 d_out[i].z = 0; // blue channel
26 d_out[i].z = 255; // fully opaque
27 }
28
29 int main()
30 {
31 uchar4 *out = (uchar4*)calloc(W*H, sizeof(uchar4));
32 uchar4 *d_out; // pointer for device array
33 cudaMalloc(&d_out, W*H*sizeof(uchar4));
34
35 const int2 pos = {0, 0}; // set reference position
36 const dim3 blockSize(TX, TY);
37 const int bx = (W + TX - 1)/TX;
38 const int by = (W + TY - 1)/TY;
39 const dim3 gridSize = dim3(bx, by);
40
41 distanceKernel<<<gridSize, blockSize>>>(d_out, W, H, pos);
42
43 // Copy results to host.
44 cudaMemcpy(out, d_out, W*H*sizeof(uchar4), cudaMemcpyDeviceToHost);
45
46 cudaFree(d_out);
47 free(out);
48 return 0;
49 }
```

Here the distance is computed in pixel spacings, so the reference position, `pos`, now has type `int2`, and the distance `d` has type `int`. The distance value is then converted to `intensity` of type `unsigned char`, whose value is restricted to the allowed range of 0 to 255 using the function `clip()`. The output arrays, `out` and `d_out`, have the corresponding vector type `uchar4`. The assignments `d_out[i].x = intensity` and `d_out[i].y = intensity` store the intensity value in the red and green channels to produce a yellow distance image. (We set the blue component to zero and the alpha to 255, corresponding to full opacity, but you should experiment with other color specifications.)