cv::dft(): The Discrete Fourier Transform

The OpenCV function cv::dft() implements one such FFT algorithm. The cv::dft() function can compute FFTs for one- and two-dimensional arrays of inputs. In the latter case, the two-dimensional transform can be computed or, if desired, only the one-dimensional transforms of each individual row can be computed (this operation is much faster than calling cv::dft() several times):

void cv::dft(
  cv::InputArray    src,                        // Input array (real or complex)
  cv::OutputArray   dst,                        // Output array
  int               flags       = 0,            // for inverse, or other options
  int               nonzeroRows = 0             // number of rows to not ignore
);

The input array must be of floating-point type and may be single- or double-channel. In the single-channel case, the entries are assumed to be real numbers, and the output will be packed in a special space-saving format called complex conjugate symmetrical (CCS).² If the source and channel are two-channel matrices or images, then the two channels will be interpreted as the real and imaginary components of the input data. In this case, there will be no special packing of the results, and some space will be wasted with a lot of 0s in both the input and output arrays.³

The special packing of result values that is used with single-channel CCS output is as follows.

For a one-dimensional array:

For a two-dimensional array:

It is worth taking a moment to look closely at the indices of these arrays. Certain values in the array are guaranteed to be 0 (more accurately, certain values of f_k are guaranteed to be real). Also note that the last row listed in the table will be present only if N_y is even and that the last column will be present only if N_x is even. In the case of the two-dimensional array being treated as N_y separate one-dimensional arrays rather than a full two-dimensional transform (we’ll take a look at how to do this), all of the result rows will be analogous to the single row listed for the output of the one-dimensional array.

The third argument, called flags, indicates exactly what operation is to be done. As usual, flags is treated as a bit array, so you can combine any flags you need with Boolean OR. The transformation we started with is known as a forward transform and is selected by default. The inverse transform⁴ is defined in exactly the same way except for a change of sign in the exponential and a scale factor. To perform the inverse transform without the scale factor, use the flag cv::DFT_INVERSE. The flag for the scale factor is cv::DFT_SCALE, which results in all of the output being scaled by a factor of N^–1 (or (N_x N_y)^–1 for a two-dimensional transform). This scaling is necessary if the sequential application of the forward transform and the inverse transform is to bring us back to where we started. Because one often wants to combine cv::DFT_INVERSE with cv::DFT_SCALE, there are several shorthand notations for this kind of operation. In addition to just combining the two operations, you can use cv::DFT_INV_SCALE (or cv::DFT_INVERSE_SCALE if you’re not into brevity). The last flag you may want to have handy is cv::DFT_ROWS, which allows you to tell cv::dft() to treat a two-dimensional array as a collection of one-dimensional arrays that should each be transformed separately as if they were N_y distinct vectors of length N_x. This can significantly reduce overhead when you’re doing many transformations at a time. By using cv::DFT_ROWS, you can also implement three-dimensional (and higher) DFT.

Though the default behavior of the forward transform is to produce results in CCS format (which results in an output array exactly the same size as the input array), you can explicitly ask OpenCV to not do this with the flag cv::DFT_COMPLEX_OUTPUT. The result will be the full complex array (with all of the zeros in it). Conversely, when you’re performing an inverse transformation on a complex array, the result is normally also a complex array. If the source array had complex conjugate symmetry,⁵ you can ask OpenCV to produce a purely real array (which will be smaller than the input array) by passing the cv::DFT_REAL_OUTPUT flag.

In order to understand the last argument, nonzero_rows, we must digress for a moment to explain that, in general, DFT algorithms strongly prefer input vectors of some lengths over input vectors of other lengths; similarly, for arrays of some sizes over arrays of other sizes. In most DFT algorithms, the preferred sizes are powers of 2 (i.e., 2ⁿ for some integer n). In the case of the algorithm used by OpenCV, the preference is that the vector lengths, or array dimensions, be 2^p3^q5^r for some integers p, q, and r. Hence the usual procedure is to create a somewhat larger array and then to copy your array into that somewhat roomier zero-padded array. For convenience, there is a handy utility function, cv::getOptimalDFTSize(), which takes the (integer) length of your vector and returns the first equal or larger size that can be expressed in the form given (i.e., 2^p3^q5^r). Despite the need for this padding, it is possible to indicate to cv::dft() that you really do not care about the transform of those rows that you had to add down below your actual data (or, if you are doing an inverse transform, which rows in the result you do not care about). In either case, you can use nonzero_rows to indicate how many rows contain meaningful data. This will provide some savings in computation time.

cv::idft(): The Inverse Discrete Fourier Transform

As we saw earlier, the function cv::dft() can be made to implement not only the discrete Fourier transform, but also the inverse operation (with the provision of the correct flags argument). It is often preferable, if only for code readability, to have a separate function that does this inverse operation by default.

void cv::idft(
  cv::InputArray  src,                          // Input array (real or complex)
  cv::OutputArray dst,                          // Output array
  int             flags       = 0,              // for variations
  int             nonzeroRows = 0               // number of rows to not ignore
);

Calling cv::idft() is exactly equivalent to calling cv::dft() with the cv::DFT_INVERSE flag (in addition to whatever flags you supply to cv::idft(), of course).

cv::mulSpectrums(): Spectrum Multiplication

In many applications that involve computing DFTs, one must also compute the per-element multiplication of the two resulting spectra. Because such results are complex numbers, typically packed in their special high-density CCS format, it would be tedious to unpack them and handle the multiplication via the “usual” matrix operations. Fortunately, OpenCV provides the handy cv::mulSpectrums() routine, which performs exactly this function for us:

void cv::mulSpectrums(
  cv::InputArray  src1,                         // Input array (ccs or complex)
  cv::InputArray  src2,                         // Input array (ccs or complex)
  cv::OutputArray dst,                          // Result array
  int             flags,                        // for row-by-row computation
  bool            conj = false                  // true to conjugate src2
);

Note that the first two arguments are arrays, which must be either CCS packed single-channel spectra or two-channel complex spectra—as you would get from calls to cv::dft(). The third argument is the destination array, which will be of the same size and type as the source arrays. The final argument, conj, tells cv::mulSpectrums() exactly what you want done. In particular, it may be set to false for implementing the preceding pair multiplication or set to true if the element from the first array is to be multiplied by the complex conjugate of the corresponding element of the second array.⁶

Convolution Using Discrete Fourier Transforms

It is possible to greatly increase the speed of a convolution by using DFT via the convolution theorem [Titchmarsh26] that relates convolution in the spatial domain to multiplication in the Fourier domain [Morse53; Bracewell65; Arfken85].⁷ To accomplish this, we first compute the Fourier transform of the image and then the Fourier transform of the convolution filter. Once this is done, we can perform the convolution in the transform space in linear time with respect to the number of pixels in the image. It is worthwhile to look at the source code for computing such a convolution, as it will also provide us with many good examples of using cv::dft(). The code is shown in Example 12-1, which is taken directly from the OpenCV reference.

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;

int main(int argc, char** argv) {

  if(argc != 2) {
    cout << "Fourier Transform
Usage: " <<argv[0] <<" <imagename>" << endl;
    return -1;
  }

  cv::Mat A = cv::imread(argv[1],0);

  if( A.empty() ) { cout << "Cannot load " << argv[1] << endl; return -1; }

  cv::Size patchSize( 100, 100 );
  cv::Point topleft( A.cols/2, A.rows/2 );
  cv::Rect roi( topleft.x, topleft.y, patchSize.width, patchSize.height );
  cv::Mat B = A( roi );

  int dft_M = cv::getOptimalDFTSize( A.rows+B.rows-1 );
  int dft_N = cv::getOptimalDFTSize( A.cols+B.cols-1 );

  cv::Mat dft_A = cv::Mat::zeros( dft_M, dft_N, CV::F32 );
  cv::Mat dft_B = cv::Mat::zeros( dft_M, dft_N, CV::F32 );

  cv::Mat dft_A_part = dft_A( Rect(0, 0, A.cols,A.rows) );
  cv::Mat dft_B_part = dft_B( Rect(0, 0, B.cols,B.rows) );

  A.convertTo( dft_A_part, dft_A_part.type(), 1, -mean(A)[0] );
  B.convertTo( dft_B_part, dft_B_part.type(), 1, -mean(B)[0] );

  cv::dft( dft_A, dft_A, 0, A.rows );
  cv::dft( dft_B, dft_B, 0, B.rows );

  // set the last parameter to false to compute convolution instead of correlation
  //
  cv::mulSpectrums( dft_A, dft_B, dft_A, 0, true );
  cv::idft( dft_A, dft_A, DFT_SCALE, A.rows + B.rows - 1 );

  cv::Mat corr = dft_A( Rect(0, 0, A.cols + B.cols - 1, A.rows + B.rows - 1) );
  cv::normalize( corr, corr, 0, 1, NORM_MINMAX, corr.type() );
  cv::pow( corr, 3., corr );

  cv::B ^= cv::Scalar::all( 255 );

  cv::imshow( "Image", A );
  cv::imshow( "Correlation", corr );
  cv::waitKey();

  return 0;

}

In Example 12-1, we can see that the input arrays are first created and then initialized. Next, two new arrays are created whose dimensions are optimal for the DFT algorithm. The original arrays are copied into these new arrays and the transforms are then computed. Finally, the spectra are multiplied together and the inverse transform is applied to the product. The transforms are the slowest⁸ part of this operation; an N × N image takes O(N² log N) time and so the entire process is also completed in that time (assuming that N > M for an M × M convolution kernel). This time is much faster than O(N²M²), the non-DFT convolution time required by the more naïve method.

cv::dct(): The Discrete Cosine Transform

For real-valued data, it is often sufficient to compute what is, in effect, only half of the discrete Fourier transform. The discrete cosine transform (DCT) [Ahmed74; Jain77] is defined analogously to the full DFT by the following formula:

Of course, there is a similar transform for higher dimensions. Note that, by convention, the normalization factor is applied to both the cosine transform and its inverse (which is not the convention for the discrete Fourier transform).

The basic ideas of the DFT apply also to the DCT, but now all the coefficients are real-valued.⁹ The actual OpenCV call is:

void cv::dct(
  cv::InputArray  src,                          // Input array (even size)
  cv::OutputArray dst,                          // Output array
  int             flags = 0                     // for row-by-row or inverse
);

The cv::dct() function expects arguments like those for cv::dft() except that, because the results are real-valued, there is no need for any special packing of the result array (or of the input array in the case of an inverse transform). Unlike cv::dft(), however, the input array must have an even number of elements (you can pad the last element with a zero if necessary to achieve this). The flags argument can be set to cv::DCT_INVERSE to generate the inverse transformation, and may be combined with cv::DCT_ROWS with the same effect as with cv::dft(). Because of the different normalization convention, both the forward and inverse cosine transforms always contain their respective contribution to the overall normalization of the transform; hence, there is no analog to cv::DFT_SCALE for cv::dct().

As with cv::dft(), performance strongly depends on array size. In fact, deep down, the implementation of cv::dct() actually calls cv::dft() on an array exactly half the size of your input array. For this reason, the optimal size of an array to pass to cv::dct() is exactly double the size of the optimal array you would pass to cv::dft(). Putting everything together, the best way to get an optimal size for cv::dct() is to compute:

size_t optimal_dct_size = 2 * cv::getOptimalDFTSize( (N+1)/2 );

where N is the actual size of your data that you want to transform.

cv::idct(): The Inverse Discrete Cosine Transform

Just as with cv::idft() and cv::dft(), cv::dct() can be asked to compute the inverse cosine transform using the flags argument. As before, code readability is often improved with the use of a separate function that does this inverse operation by default:

void cv::idct(
  cv::InputArray  src,                          // Input array
  cv::OutputArray dst,                          // Output array
  int             flags = 0,                    // for row-by-row computation
);

Calling cv::idct() is exactly equivalent to calling cv::dct() with the cv::DCT_INVERSE flag (in addition to any other flags you supply to cv::idct()). 

cv::integral() for Squared Summation Integral

The squared sum is computed with the same function as the regular sum, except for the provision of an additional output argument for the squared sum.

void cv::integral(
  cv::InputArray  image,                        // Input array
  cv::OutputArray sum,                          // Output sum results
  cv::OutputArray sqsum,                        // Output sum of squares results
  int             sdepth = -1                   // Results depth (e.g., cv::F32)
);

The cv::OutputArray argument sqsum indicates to cv::integral() that the square sum should be computed in addition to the regular sum. As before, sdepth specifies the desired depth of the resulting images. sdepth can be any of CV::S32, CV::F32, or CV::F64.

cv::integral() for Tilted Summation Integral

Similar to the squared sum, the tilted sum integral is essentially the same function, with an additional argument for the additional result.

void cv::integral(
  cv::InputArray  image,                        // Input array
  cv::OutputArray sum,                          // Output sum results
  cv::OutputArray sqsum,                        // Output sum of squares results
  cv::OutputArray tilted,                       // Output tilted sum results
  int             sdepth = -1                   // Results depth (e.g., cv::F32)
);

The additional cv::OutputArray argument tilted is computed by this form of cv::integral(), in addition to the other sums; thus, all of the other arguments are the same.

cv::Canny()

The OpenCV implementation of the Canny edge detection algorithm converts an input image into an “edge image.”

void cv::Canny(
  cv::InputArray  image,                        // Input single channel image
  cv::OutputArray edges,                        // Output edge image
  double          threshold1,                   // "lower" threshold
  double          threshold2,                   // "upper" threshold
  int             apertureSize = 3,             // Sobel aperture
  bool            L2gradient   = false          // true=L2-norm (more accurate)
);

The cv::Canny() function expects an input image, which must be single-channel, and an output image, which will also be grayscale (but which will actually be a Boolean image). The next two arguments are the low and high thresholds. The next-to-last argument, apertureSize, is the aperture used by the Sobel derivative operators that are called inside of the implementation of cv::Canny(). The final argument L2gradient is used to select between computing the directional gradient “correctly” using the proper L₂-norm, or using a faster, less accurate L₁-norm-based method. If the argument L2gradient is set to true, the more accurate form is used:

If L2gradient is set to false, the faster form is used:

Hough Line Transform

The basic theory of the Hough line transform is that any point in a binary image could be part of some set of possible lines. If we parameterize each line by, for example, a slope a and an intercept b, then a point in the original image is transformed to a locus of points in the (a, b) plane corresponding to all of the lines passing through that point, of which it could potentially be a part (see Figure 12-5). If we convert every nonzero pixel in the input image into such a set of points in the output image and sum over all such contributions, then lines that appear in the input (i.e., (x, y) plane) image will appear as local maxima in the output (i.e., (a, b) plane) image. Because we are summing the contributions from each point, the (a, b) plane is commonly called the accumulator plane.

Figure 12-5. The Hough line transform finds many lines in each image; some of the lines found are expected, but others may not be

It might occur to you that the slope-intercept form is not really the best way to represent all the lines passing through a point (i.e. because of the considerably different density of lines as a function of the slope, and the related fact that the interval of possible slopes goes from –∞ to +∞). This is why the actual parameterization of the transform image used in numerical computation is somewhat different. The preferred parameterization represents each line as a point in polar coordinates (ρ, θ), with the implied line being the line passing through the indicated point but perpendicular to the radial from the origin to that point (see Figure 12-6). The equation for such a line is:

Figure 12-6. A point (x₀, y₀) in the image plane (a) implies many lines, each parameterized by a different ρ and θ (b); these lines in the (ρ, θ) plane, taken together, form a curve of characteristic shape (c)

The OpenCV Hough transform algorithm does not make this computation explicit to the user. Instead, it simply returns the local maxima in the (ρ, θ) plane. However, you will need to understand this process in order to understand the arguments to the OpenCV Hough line transform function.

OpenCV supports three different kinds of Hough line transform: the standard Hough transform (SHT) [Duda72], the multiscale Hough transform (MHT), and the progressive probabilistic Hough transform (PPHT).¹⁵ The SHT is the algorithm we just covered. The MHT algorithm is a slight refinement that gives more accurate values for the matched lines. The PPHT is a variation of this algorithm that, among other things, computes an extent for individual lines in addition to the orientation (as shown in Figure 12-7). It is called “probabilistic” because, rather than accumulating every possible point in the accumulator plane, it accumulates only a fraction of them. The idea is that if the peak is going to be high enough anyhow, then hitting it only a fraction of the time will be enough to find it; the result of this conjecture can be a substantial reduction in computation time.

Figure 12-7. The Canny edge detector (param1=50, param2=150) is run first, with the results shown in gray, and the progressive probabilistic Hough transform (param1=50, param2=10) is run next, with the results overlaid in white; you can see that the strong lines are generally picked up by the Hough transform

void cv::HoughLines(
  cv::InputArray  image,                  // Input single channel image
  cv::OutputArray lines,                  // N-by-1 two-channel array
  double          rho,                    // rho resolution (pixels)
  double          theta,                  // theta resolution (radians)
  int             threshold,              // Unnormalized accumulator threshold
  double          srn      = 0,           // rho refinement (for MHT)
  double          stn      = 0            // theta refinement (for MHT)
);

void cv::HoughLinesP(
  cv::InputArray  image,                  // Input single channel image
  cv::OutputArray lines,                  // N-by-1 4-channel array
  double          rho,                    // rho resolution (pixels)
  double          theta,                  // theta resolution (radians)
  int             threshold,              // Unnormalized accumulator threshold
  double          minLineLength = 0,      // required line length
  double          maxLineGap    = 0       // required line separation
);

Hough Circle Transform

The Hough circle transform [Kimme75] (see Figure 12-8) works in a manner roughly analogous to the Hough line transforms just described. The reason it is only “roughly” is that—if you were to try doing the exactly analogous thing—the accumulator plane would have to be replaced with an accumulator volume with three dimensions: one for x and one for y (the location of the circle center), and another for the circle radius r. This would mean far greater memory requirements and much slower speed. The implementation of the circle transform in OpenCV avoids this problem by using a somewhat trickier method called the Hough gradient method.

Figure 12-8. The Hough circle transform finds some of the circles in the test pattern and (correctly) finds none in the photograph

The Hough gradient method works as follows. First, the image is passed through an edge-detection phase (in this case, cv::Canny()). Next, for every nonzero point in the edge image, the local gradient is considered (we compute the gradient by first computing the first-order Sobel x- and y-derivatives via cv::Sobel()). Using this gradient, we increment every point along the line indicated by this slope—from a specified minimum to a specified maximum distance—in the accumulator. At the same time, the location of every one of these nonzero pixels in the edge image is noted. The candidate centers are then selected from those points in this (two-dimensional) accumulator that are both above some given threshold and larger than all of their immediate neighbors. These candidate centers are sorted in descending order of their accumulator values, so that the centers with the most supporting pixels appear first. Next, for each center, all of the nonzero pixels (recall that this list was built earlier) are considered. These pixels are sorted according to their distance from the center. Working out from the smallest distances to the maximum radius, we select a single radius that is best supported by the nonzero pixels. A center is kept if it has sufficient support from the nonzero pixels in the edge image and if it is a sufficient distance from any previously selected center.

This implementation enables the algorithm to run much faster and, perhaps more importantly, helps overcome the problem of the otherwise sparse population of a three-dimensional accumulator, which would lead to a lot of noise and render the results unstable. On the other hand, this algorithm has several shortcomings that you should be aware of.

First, the use of the Sobel derivatives to compute the local gradient—and the attendant assumption that this can be considered equivalent to a local tangent—is not a numerically stable proposition. It might be true “most of the time,” but you should expect this to generate some noise in the output.

Second, the entire set of nonzero pixels in the edge image is considered for every candidate center; hence, if you make the accumulator threshold too low, the algorithm will take a long time to run. Third, because only one circle is selected for every center, if there are concentric circles then you will get only one of them.

Finally, because centers are considered in ascending order of their associated accumulator value and because new centers are not kept if they are too close to previously accepted centers, there is a bias toward keeping the larger circles when multiple circles are concentric or approximately concentric. (It is only a “bias” because of the noise arising from the Sobel derivatives; in a smooth image at infinite resolution, it would be a certainty.)

void cv::HoughCircles(
  cv::InputArray  image,                  // Input single channel image
  cv::OutputArray circles,                // N-by-1 3-channel or vector of Vec3f
  int             method,                 // Always cv::HOUGH_GRADIENT
  double          dp,                     // Accumulator resolution (ratio)
  double          minDist,                // Required separation (between lines)
  double          param1    = 100,        // Upper Canny threshold
  double          param2    = 100,        // Unnormalized accumulator threshold
  int             minRadius = 0,          // Smallest radius to consider
  int             maxRadius = 0           // Largest radius to consider
);

#include <opencv2/opencv.hpp>
#include <iostream>
#include <math.h>

using namespace cv;
using namespace std;

int main(int argc, char** argv) {

  if(argc != 2) {
    cout << "Hough Circle detect
Usage: " <<argv[0] <<" <imagename>
" << endl;
    return -1;
  }

  cv::Mat src, image;

  src  = cv::imread( argv[1], 1 );
  if( src.empty() ) { cout << "Cannot load " << argv[1] << endl; return -1; }

  cv::cvtColor(src, image, cv::BGR2GRAY);
  cv::GaussianBlur(image, image, Size(5,5), 0, 0);

  vector<cv::Vec3f> circles;
  cv::HoughCircles(image, circles, cv::HOUGH_GRADIENT, 2, image.cols/10);

  for( size_t i = 0; i < circles.size(); ++i ) {
    cv::circle(
      src,
      cv::Point(cvRound(circles[i][0]), cvRound(circles[i][1])),
      cvRound(circles[i][2]),
      cv::Scalar(0,0,255),
      2,
      cv::AA
    );
  }

  cv::imshow( "Hough Circles", src);
  cv::waitKey(0);

  return 0;

}

cv::distanceTransform() for Unlabeled Distance Transform

When you call the OpenCV distance transform function, the output image will be a 32-bit floating-point image (i.e., CV::F32).

void cv::distanceTransform(
  cv::InputArray  src,                       // Input image
  cv::OutputArray dst,                       // Result image
  int             distanceType,              // Distance metric to use
  int             maskSize                   // Mask to use (3, 5, or see below)
);

cv::distanceTransform() has two parameters. The first is distanceType, which indicates the distance metric to be used. Your choices here are cv::DIST_C, cv::DIST_L1, and cv::DIST_L2. These methods compute the distance to the nearest zero based on integer steps along a grid. The difference between the methods is that cv::DIST_C is the distance when the steps are counted on a four-connected grid (i.e., diagonal moves are not allowed), and cv::DIST_L1 gives the number of steps on an eight-connected grid (i.e., diagonal moves are allowed). When distanceType is set to cv::DIST_L2, cv::distanceTransform() attempts to compute the exact Euclidean distances.

After the distance type is the maskSize, which may be 3, 5, or cv::DIST_MASK_PRECISE. In the case of 3 or 5, this argument indicates that a 3 × 3 or 5 × 5 mask should be used with the Borgefors method. If you are using cv::DIST_L1 or cv::DIST_C, you can always use a 3 × 3 mask and you will get exact results. If you are using cv::DIST_L2, the Borgefors method is always approximate, and using the larger 5 × 5 mask will result in a better approximation to the L2 distance, at the cost of a slightly slower computation. Alternatively, cv::DIST_MASK_PRECISE can be used to indicate the Felzenszwalb algorithm (when used with cv::DIST_L2).

cv::distanceTransform() for Labeled Distance Transform

It is also possible to ask the distance transform algorithm to not only calculate the distances, but to also report which object that minimum distance is to. These “objects” are called connected components. We will have a lot more to say about connected components in Chapter 14 but, for now, you can think of them as exactly what they sound like: structures made of continuously connected groups of zeros in the source image.

void cv::distanceTransform(
  cv::InputArray  src,                             // Input image
  cv::OutputArray dst,                             // Result image
  cv::OutputArray labels,                          // Connected component ids
  int             distanceType,                    // Distance metric to use
  int             maskSize,                        // (3, 5, or see below)
  int             labelType = cv::DIST_LABEL_CCOMP // How to label
);

If a labels array is provided, then as a result of running cv::distanceTransform() it will be of the same size as dst. In this case, connected components will be computed automatically, and the label associated with the nearest such component will be placed in each pixel of labels. The output “labels” array will basically be the discrete Voronoi diagram.

The argument labelType can be set either to cv::DIST_LABEL_CCOMP or cv::DIST_LABEL_PIXEL. In the former case, the function automatically finds connected components of zero pixels in the input image and gives each one a unique label. In the latter case, all zero pixels are given distinct labels.

Flood Fill

Flood fill [Heckbert90; Shaw04; Vandevenne04] is an extremely useful function that is often used to mark or isolate portions of an image for further processing or analysis. Flood fill can also be used to derive, from an input image, masks that can be used by subsequent routines to speed or restrict processing to only those pixels indicated by the mask. The function cv::floodFill() itself takes an optional mask that can be further used to control where filling is done (e.g., for multiple fills of the same image).

In OpenCV, flood fill is a more general version of the sort of fill functionality that you probably already associate with typical computer painting programs. For both, a seed point is selected from an image and then all similar neighboring points are colored with a uniform color. The difference is that the neighboring pixels need not all be identical in color.¹⁹ The result of a flood fill operation will always be a single contiguous region. The cv::floodFill() function will color a neighboring pixel if it is within a specified range (loDiff to upDiff) of either the current pixel or if (depending on the settings of flags) the neighboring pixel is within a specified range of the original seed value. Flood filling can also be constrained by an optional mask argument. There are two different prototypes for the cv::floodFill() routine, one that accepts an explicit mask parameter, and one that does not.

int cv::floodFill(
  cv::InputOutputArray image,                  // Input image, 1 or 3 channels
  cv::Point            seed,                   // Start point for flood
  cv::Scalar           newVal,                 // Value for painted pixels
  cv::Rect*            rect,                   // Output bounds painted domain
  cv::Scalar           lowDiff  = cv::Scalar(),// Maximum down color distance
  cv::Scalar           highDiff = cv::Scalar(),// Maximum up color distance
  int                  flags                   // Local/global, and mask-only
);

int cv::floodFill(
  cv::InputOutputArray image,                   // Input w-by-h, 1 or 3 channels
  cv::InputOutputArray mask,                    // 8-bit, w+2-by-h+2 (Nc=1)
  cv::Point            seed,                    // Start point for flood
  cv::Scalar           newVal,                  // Value for painted pixels
  cv::Rect*            rect,                    // Output bounds painted domain
  cv::Scalar           lowDiff  = cv::Scalar(), // Maximum down color distance
  cv::Scalar           highDiff = cv::Scalar(), // Maximum up color distance
  int                  flags                    // Local/global, and mask-only
);

The parameter image is the input image, which can be 8-bit or a floating-point type, and must either have one or three channels. In general, this image array will be modified by cv::floodFill(). The flood-fill process begins at the location seed. The seed will be set to value newVal, as will all subsequent pixels colored by the algorithm. A pixel will be colorized if its intensity is not less than a colorized neighbor’s intensity minus loDiff and not greater than the colorized neighbor’s intensity plus upDiff. If the flags argument includes cv::FLOODFILL_FIXED_RANGE, then a pixel will be compared to the original seed point rather than to its neighbors. Generally, the flags argument controls the connectivity of the fill, what the fill is relative to, whether we are filling only a mask, and what values are used to fill the mask. Our first example of flood fill is shown in Figure 12-10.

Figure 12-10. Results of flood fill (top image is filled with gray, bottom image with white) from the dark circle located just off center in both images; in this case, the upDiff and loDiff parameters were each set to 7.0

The mask argument indicates a mask that can function both as an input to cv::floodFill() (in which case, it constrains the regions that can be filled) and as an output from cv::floodFill() (in which case, it will indicate the regions that actually were filled). mask must be a single-channel, 8-bit image whose size is exactly two pixels larger in width and height than the source image.²⁰

In the sense that mask is an input to cv::floodFill(), the algorithm will not flood across nonzero pixels in the mask. As a result, you should zero it before use if you don’t want masking to block the flooding operation.

When the mask is present, it will also be used as an output. When the algorithm runs, every “filled” pixel will be set to a nonzero value in the mask. You have the option of adding the value cv::FLOODFILL_MASK_ONLY to flags (using the usual Boolean OR operator). In this case, the input image will not be modified at all. Instead, only mask will be modified.

Two possible values for flags have already been mentioned: cv::FLOODFILL_FIXED_RANGE and cv::FLOODFILL_MASK_ONLY. In addition to these, you can also add the numerical values 4 or 8.²¹ In this case, you are specifying whether the flood-fill algorithm should consider the pixel array to be four-connected or eight-connected. In the former case, a four-connected array is one in which pixels are only connected to their four nearest neighbors (left, right, above, and below). In the eight-connected case, pixels are considered to be connected to diagonally neighboring pixels as well.

The flags argument is slightly tricky because it has three parts that are possibly less intuitive than they could be. The low 8 bits (0–7) can be set to 4 or 8, as we saw before, to control the connectivity considered by the filling algorithm. The high 8 bits (16–23) are the ones that contain the flags such as cv::FLOODFILL_FIXED_RANGE and cv::FLOODFILL_MASK_ONLY. The middle bits (8–15), however, are used a little bit differently to actually represent a numerical value: the value with which you want the mask to be filled. If the middle bits of flags are 0s, the mask will be filled with 1s (the default), but any other value will be interpreted as an 8-bit unsigned integer. All these flags may be linked together via OR. For example, if you want an eight-way connectivity fill, filling only a fixed range, filling the mask (not the image), and filling using a value of 47, then the parameter to pass in would be:

flags = 8
      | cv::FLOODFILL_MASK_ONLY
      | cv::FLOODFILL_FIXED_RANGE
      | (47<<8);

Figure 12-11 shows flood fill in action on a sample image. Using cv::FLOODFILL_FIXED_RANGE with a wide range resulted in most of the image being filled (starting at the center). We should note that newVal, loDiff, and upDiff are prototyped as type cv::Scalar so they can be set for three channels at once. For example, loDiff = cv::Scalar(20,30,40) will set loDiff thresholds of 20 for red, 30 for green, and 40 for blue.

Figure 12-11. Results of flood fill (top image is filled with gray, bottom image with white) from the dark circle located just off center in both images; in this case, flood fill was done with a fixed range and with a high and low difference of 25.0

Watershed Algorithm

In many practical contexts, we would like to segment an image but do not have the benefit of any kind of separate background mask. One technique that is often effective in this context is the watershed algorithm [Meyer92], which converts lines in an image into “mountains” and uniform regions into “valleys” that can be used to help segment objects. The watershed algorithm first takes the gradient of the intensity image; this has the effect of forming valleys or basins (the low points) where there is no texture, and of forming mountains or ranges (high ridges corresponding to edges) where there are dominant lines in the image. It then successively floods basins starting from caller-specified points until these regions meet. Regions that merge across the marks so generated are segmented as belonging together as the image “fills up.” In this way, the basins connected to the marker point become “owned” by that marker. We then segment the image into the corresponding marked regions.

More specifically, the watershed algorithm allows a user (or another algorithm) to mark parts of an object or background that are known to be part of the object or background. Alternatively, the caller can draw a simple line or collection of lines that effectively tells the watershed algorithm to “group points like these together.” The watershed algorithm then segments the image by allowing marked regions to “own” the edge-defined valleys in the gradient image that are connected with the segments. Figure 12-12 clarifies this process.

Figure 12-12. Watershed algorithm: after a user has marked objects that belong together (left), the algorithm merges the marked area into segments (right)

The function specification of the watershed segmentation algorithm is:

void cv::watershed(
  cv::InputArray       image,                 // Input 8-bit, three channels
  cv::InputOutputArray markers                // 32-bit float, single channel
);

Here, image must be an 8-bit, three-channel (color) image and markers is a single-channel integer (CV::S32) image of the same (x, y) dimensions. On input, the value of markers is 0 except where the caller has indicated by using positive numbers that some regions belong together. For example, in the left panel of Figure 12-12, the orange might have been marked with a “1,” the lemon with a “2,” the lime with “3,” the upper background with “4,” and so on.

After the algorithm has run, all of the former zero-value pixels in markers will be set to one of the given markers (i.e., all of the pixels of the orange are hoped to come out with a “1” on them, the pixels of the lemon with a “2,” etc.), except the boundary pixels between regions, which will be set to –1. Figure 12-12 (right) shows an example of such a segmentation.

Grabcuts

The Grabcuts algorithm, introduced by Rother, Kolmogorov, and Blake [Rother04], extends the Graphcuts algorithm [Boykov01] for use in user-directed image segmentation. The Grabcuts algorithm is capable of obtaining excellent segmentations, often with no more than a bounding rectangle around the foreground object to be segmented.

The original Graphcuts algorithm used user-labeled foreground and user-labeled background regions to establish distribution histograms for those two classes of image regions. It then combined the assertion that unlabeled foreground or background should conform to similar distributions with the idea that these regions tend to be smooth and connected (i.e., a bunch of blobs). These assertions were then combined into an energy functional that gave a low energy (i.e., cost) to solutions that conformed to these assertions and a high energy to solutions that violated them. The algorithm obtained the final result by minimizing this energy function.²²

The Grabcuts algorithm extends Graphcuts in several important ways. The first is that it replaces the histogram models with a different (Gaussian mixture) model, enabling the algorithm to work on color images. In addition, it solves the energy functional minimization problem in an iterative manner, which provides better results overall, and allows much greater flexibility in the labeling provided by the user. Notably, this latter point makes possible even one-sided labelings, which identify either only background or only foreground pixels (where Graphcuts required both).

The implementation in OpenCV allows the caller to either just provide a rectangle around the object to be segmented, in which case the pixels “under” the rectangle’s boundary (i.e., outside of it) are taken to be background and no foreground is specified. Alternatively, the caller can specify an overall mask in which pixels are categorized as being either definitely foreground, definitely background, probably foreground, or probably background.²³ In this case, the definite regions will be used to classify the other regions, with the latter being classified into the definite categories by the algorithm.

The OpenCV implementation of Grabcuts is implemented by the cv::Grabcuts() function:

void cv::grabCut(
  cv::InputArray       img,
  cv::InputOutputArray mask,
  cv::Rect             rect,
  cv::InputOutputArray bgdModel,
  cv::InputOutputArray fgdModel,
  int                  iterCount,
  int                  mode     = cv::GC_EVAL
);

Given an input image img, the resulting labeling will be computed by cv::grabCut() and placed in the output array mask. This mask array can also be used as an input. This is determined by the mode variable. If mode contains the flag cv::GC_INIT_WITH_MASK,²⁴ then the values currently in mask when it is called will be used to initialize the labeling of the image. The mask is expected to be a single-channel image of cv::U8 type in which every value is one of the following enumerations.

The argument rect is used only when you are not using mask initialization. When the mode contains the flag cv::GC_INIT_WITH_RECT, the entire region outside of the provided rectangle is taken to be “definitely background,” while the rest is automatically set to “probably foreground.”

The next two arrays are essentially temporary buffers. When you first call cv::grabCut(), they can be empty arrays. However, if for some reason you should run the Grabcuts algorithm for some number of iterations and then want to restart the algorithm for more iterations (possibly after allowing a user to provide additional “definite” pixels to guide the algorithm), you will need to pass in the same (unmodified) buffers that were filled by the previous run (in addition to using the mask you got back from the previous run as input for the next run).

Internally, the Grabcuts algorithm essentially runs the Graphcuts algorithm some number of times (with the minor extensions mentioned previously). In between each such run, the mixture models are recomputed. The itercount argument determines how many such iterations will be applied. Typical values for itercount are 10 or 12, though the number required may depend on the size and nature of the image being processed.

Mean-Shift Segmentation

Mean-shift segmentation finds the peaks of color distributions over space [Comaniciu99]. It is related to the mean-shift algorithm, which we will discuss later when we talk about tracking and motion in Chapter 17. The main difference between the two is that the former looks at spatial distributions of color (and is thus related to our current topic of segmentation), while the latter tracks those distributions through time in successive frames. The function that does this segmentation based on the color distribution peaks is cv::pyrMeanShiftFiltering(). 

Given a set of multidimensional data points whose dimensions are (x, y, blue, green, red), mean-shift can find the highest density “clumps” of data in this space by scanning a window over the space. Notice, however, that the spatial variables (x, y) can have very different ranges from the color-magnitude ranges (blue, green, red). Therefore, mean-shift needs to allow for different window radii in different dimensions. In this case we should have, at a minimum, one radius for the spatial variables (spatialRadius) and one radius for the color magnitudes (colorRadius). As mean-shift windows move, all the points traversed by the windows that converge at a peak in the data become connected to or “owned” by that peak. This ownership, radiating out from the densest peaks, forms the segmentation of the image. The segmentation is actually done over a scale pyramid—cv::pyrUp(), cv::pyrDown()—as described in Chapter 11, so that color clusters at a high level in the pyramid (shrunken image) have their boundaries refined at lower levels in the pyramid.

The output of the mean-shift segmentation algorithm is a new image that has been “posterized,” meaning that the fine texture is removed and the gradients in color are mostly flattened. You can then further segment such images using whatever algorithm is appropriate for your needs (e.g., cv::Canny() combined with cv::findContours(), if in fact a contour segmentation is what you ultimately want).

The function prototype for cv::pyrMeanShiftFiltering() looks like this:

void cv::pyrMeanShiftFiltering(
  cv::InputArray   src,                       // 8-bit, Nc=3 image
  cv::OutputArray  dst,                       // 8-bit, Nc=3, same size as src
  cv::double       sp,                        // Spatial window radius
  cv::double       sr,                        // Color window radius
  int              maxLevel = 1,              // Max pyramid level
  cv::TermCriteria termcrit = TermCriteria(
    cv::TermCriteria::MAX_ITER | cv::TermCriteria::EPS,
    5,
    1
  )
);

In cv::pyrMeanShiftFiltering() we have an input image src and an output image dst. Both must be 8-bit, three-channel color images of the same width and height. The spatialRadius and colorRadius define how the mean-shift algorithm averages color and space together to form a segmentation. For a 640 × 480 color image, it works well to set spatialRadius equal to 2 and colorRadius equal to 40. The next parameter of this algorithm is max_level, which describes how many levels of scale pyramid you want used for segmentation. A max_level of 2 or 3 works well for a 640 × 480 color image.

The final parameter is cv::TermCriteria, which we have seen in some previous algorithms. cv::TermCriteria is used for all iterative algorithms in OpenCV. The mean-shift segmentation function comes with good defaults if you just want to leave this parameter blank.

Figure 12-13 shows an example of mean-shift segmentation using the following values:

cv::pyrMeanShiftFiltering( src, dst, 20, 40, 2);

Figure 12-13. Mean-shift segmentation over scale using cv::pyrMeanShiftFiltering() with parameters max_level=2, spatialRadius=20, and colorRadius=40; similar areas now have similar values and so can be treated as super pixels (larger statistically similar areas), which can speed up subsequent processing significantly