CUDA for Engineers: 2D Grids and Interactive Graphics

This chapter is from the book

CUDA for Engineers: An Introduction to High-Performance Parallel Computing

Learn More Buy

This chapter is from the book

This chapter is from the book 

CUDA for Engineers: An Introduction to High-Performance Parallel Computing

Learn More Buy

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 ¹ and OpenGL ². 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 ³ 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 }

void kernel2D(uchar4 *d_out, int w, int h, ... )

#define W 500
#define H 500
kernel2D<<<gridSize, blockSize>>>(d_out, W, H, ... )

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.)