Analyzing Oriented Patterns

I Introduction

A central focus in recent computational vision has been the de­ composition of the original intensity image into intrinsic images (Horn 1977; Barrow & Tenenbaum, 1978; Marr, 1982), repre­ senting such properties as depth, reflectance, and illuminance. These intrinsic properties are believed to be more meaningful than image intensity because they describe basic independent constituents of the image formation process. Thus, for exam­ ple, in separating shape from illumination, we can recognize an invariance of shape regardless of changing illumination.

The advantages of decomposing what we see into its more-or­ less independent parts extends beyond the image formation pro­ cess to the shapes and patterns on which that process operates. For instance, decomposing a bent rod into a straight rod and a bending transformation reveals the similarity between a bent rod and one that hasn't been bent, or some other solid that's been bent the same way (Barr, 1984).

Just as we need to understand the image-forming process to decompose an image into intrinsic images, we need to under­ stand the processes that generate patterns to decompose them into their intrinsic parts. But, while there is only one image­ forming process, a staggering variety of processes shape and color the world around us. Our only hope of dealing with this com­ plexity is to begin with some basic pattern classes that recur in nature, and understand how to decompose and describe them.

One such dass are oriented patterns, notably those produced by propagation, accretion, or deformation. To understand an oriented pattern we must be able to say (1) what is propagat­ ing, accreting, or deforming, and (2) which way and how much. More precisely, we must estimate everywhere the direction and magnitude of anistropy (which we will call the flow field,) and describe the residual pattern, independent of that field. Why this decomposition Ieads to simpler, more regular descriptions is best illustrated by example:

• A typical oriented pattern created by propagation is the streaked trail left by a paint brush dipped in variegated paint. The flow field describes the trajectory of the brush, the residual pattern depending only on the distribution of paint on the brush.
• Accretion typically results in laminar structures, such as wood grain. Here, the flow field gives isochrones (the mov­ ing accretion boundary,) and the residual pattern describes the change in color or brightness of the accreting material over time.
• If an isotropic body is deformed, the flow field principally describes the bending and stretching it has undergone, while the residual pattern describes the undeformed body.

In all these cases, separate descriptions of the flow field and the residual pattern are appropriate because they describe different processes. The path of propagation for many physical processes is controlled by very different mechanisms than control the col­ oration of the trailleft behind. Similarly, the mechanisms which control the shape of an accretion boundary are frequently unre­ lated to the processes controlling the color of the accreted ma­ terial. Finally, the forces which deform a piece of material are often completely unrelated to the process which created the piece of material in the first place. By separately describing these pro­ cesses, we can create descriptions of the whole which are often simpler than is possible without the separation because each of the pieces may have different regularities.

Orientation selective mechanisms have been extensively stud­ ied by physiologists since Hubel and Wiesel's (1962) discovery of orientation selective cells in mammalian visual cortex (see Schiller et. al. (1976) for a comprehensive example). There has also been considerable interest among psychologists in the perception of oriented patterns, particularly dot patterns (Glass, 1969). Only recently have the computational issues involved re­ ceived attention. Stevens (1978) examined the grouping of to­ kens in Glass pattems based on orientation. While successful with Glass patterns, his methods were never extended to natural imagery. Zucker (1983) investigated the estimation of orientation by combining the outputs of linear operators. Zucker's estima­ tion method for what he calls "Type II" patterns, while differing in many respects, is quite close in spirit to our own.

Little progress has been made in using local orientation esti­ mates to interpret patterns, perhaps because reliable estimates have proved difficult to obtain. The key difference between our work and earlier efforts lies in our use of the flow field to build a natural coordinate system for analyzing the pattern.

The remainder of the paper covers the computation of the flow field by local estimation of orientation, the construction of a coordinate system using the flow field, and some examples of analysis and description using flow coordinates.

In Section 2 we develop an estimator for the local flow direc­ tion, that direction in which intensity tends to vary most slowly due to an underlying anisotropic process. The estimator, based on the direction of least spatial variance in the output of an ori­ ented filter, is compu ted as follows: After initial filtering, the intensity gradient is measured at each point in the image. The gradient angle, 8, is then doubled (by treating the gradient vec­ tors as complex numbers and squaring them) to map directions differing by w into a single direction. The transformed vectors are then summed over a weighted neighborhood araund the point of interest. The angle of the summed vector is halved, undoing the previous transformation. This gives an estimate for the direction of greatest variance, which is then rotated by 1r / 2 to yield the ßow direction.

In section 3, we describe the construction and use of coordi­ nate systems based on the result of local estimation. Integral curves in the ßow field are computed numerically, by following the estimated vectors from point to point. A coordinate system is constructed in which the integral curves are parameter lines. Transforming the image into these "ßow coordinates" straight­ ens the pattern, removing the effects of changing orientation. We present several examples of analysis and description of the ßow field and the straightened pattern.

II Flow Computation

For intensity patterns created by anisotropic processes such as propagation, accretion, or deformation, variation in the ßow direction is much slower than variation in the perpendicular di­ rection. Anisotropy in such patternswill be evident in the local power spectrum. The high frequency energy will tend to duster along the line in the Fourier domain perpendicular to the ßow orientation.

A simple way to detect this dustering is to sum the energy in an appropriate region of the power spectrum and examine how the sum is affected by rotations. This can be done by examining the energy in the output of an appropriate orientation-selective linear filter. Theorientation at which the energy is maximal can be expected to be perpendicular to the ßow orientation.

Selection of the filter involves a nurober of tradeoffs. Very low spatial frequencies are affecte.d more strongly by illumina­ tion effects than surface coloration, so they are inappropriate for measuring textural anisotropy. Very high spatial frequencies are sensitive to noise and aliasing effects so they too are inappropri­ ate . Hence some type of roughly bandpass filtering is required. The orientation specificity of the filter is also quite important. H the filter is too orientation-specific then a !arge spatial neighbor­ hood will be required in order to make a reliable measurement of the energy. Conversely, if the filter responds over a wide range of orientations then it will be difficult to localize the orientation very accurately. Thus there is a trade-off between angular- and spatial- resolution.

One reasonable choice for the frequency response of the filter is

The filter is bandpass with passband determined by u1 and o-2· In our experience, ratios of the sigmas in the range of 2.0 to 10.0 work weil. The orientation specificity or tuning curve is provided by the cosine dependence of the filter on 8. This appears to strike a reasonable balance between angular- and spatial- resolutions for the range of patterns we have examined. The filter's power spectrum is shown in figure 1

The cosine orientation tuning-curve of the filter has some un­ usually good properties for computing the filter output at dif­ ferent orientations. The impulse response S(x, y) of the filter is where

is an isotropic filter. Let C = H · I and Iet Ro[SJ denote a counter-dockwise rotation of S by an angle 8. Then the convo­ lution R6 [SJ · I is just the directional derivative of H · I in the 8 direction. The directional derivative can easily be written in terms of the gradient so we have

Thus a single convolution suffices for all orientations.

Since the filter S severely attenuates very low frequencies Ro[SJ · I can be safely regarded as zero-mean. Thus the vari­ ance in its output can be estimated by the expression

where W(x, y) is a local weighting function with unit integral. We use Gaussian weighting functions W ( x, y) because approximate Gaussian convolutions can be computed efficiently (Burt 1979) .

Using the gradient formulation of the filter output in equation 2, we can write the variance V(O) as

A. Interpretation of Filter Output

There remains the issue of interpreting V (8). Assurne that there is only one axis of anisotropy. Then V (8) will have two extrema 1r apart corresponding to that axis. Let V2 ( 8) = V (8 / 2). Then V2(8) will have a single extremum in the interval 0 < {I < 21r. A computationally inexpensive way of estimating the po­ sition of this extremum is to consider V2 as a distribution and compute its mean. Since {I is periodic, V2 should be considered as a distribution on the unit circle. Hence its mean is the vector integral (a,ß) = J~"V2(8)(cos8,sin8)d8. The angle tan - 1 (ß/a) is an estimate of the angle of the peak in V2 and hence twice the angle of the peak in V . Thus the angle 4> of greatest variance can be written

These integrals a:re evaluated in Appendix A to show that the angle of anisotropy 4> can be written

which directly yields a simple algorithm for computing ,P.

B. Combining Gradient Orientations

Notice that the right hand side of equation 5 can be rega:rded as the orientation of a locally weighted sum of the vectors of the form J(x, y) = (c;- c;, 2C,Cvl · These vectors are related in a simple way to the gradient vectors G(x,y) = (C"Cv) · The mag­ nitude of J(x, y) is just the squa:re of the magnitude of G(x, y) and the angle between J ( x, y) and the x-axis is twice the an­ gle between G(x, y) and the x-axis. This follows easily from the observation that (C, + Cvi) 2 = c;- c; + 2C,Cvi.

One might be tempted to believe that smoothing the gradient vectors G(x, y) would be nea:rly as good a measure of anisotropy as smoothing the rotated squa:red gradient vectors J(x, y). This is emphatically not the case. Consider an intensity ridge such as l(x , y) = exp( -x 2 ). The gradient vectors on the left half­ plane all point to the right and the gradient vectors on the right half-plane all point to the left. Adding them tagether results in cancellation. By contrast, if they a:re first rotated to form the J vectors, they reinforce. The types of oriented patterns we a:re concerned with often have nea:rly symmetric distributions of gra­ dient directions a:round the axis of anisotropy. In such patterns, if the gradients a:re added tagether directly, the cancellation is so severe that the result often has little relation to the direction of anisotropy. Thus the difference between rotating the gradient vectors or leaving them be is often the difference between being able or unable to detect the anisotropy. Note also that smooth­ ing the image first and then computing the gradients is exactly the same as computing the gradients and then smoothing. lt will not avoid the difficulties of cancellation.

C. Coherence

In addition to finding the direction of anisotropy, it is im­ portant to determine how strong an anisotropy there is. If the orientation of the local J vectors a:re nea:rly uniformly distributed between 0 and 271", then the orientation 4> of slight anisotropy is not very meaningful. Conversely, if all the J vectors a:re pointing the same way then the indication of anisotropy is quite strong and 4> is very meaningful. A simple way of measuring the strength of the peak in the distribution of J vectors is to look at the ra­ tio x(x, y) = lW · JI /W · IJI which we will call the coherence of the ßow pattern. If the J vectors are close to uniformly dis­ tributed, then the ratio will be nearly zero. If the J vectors all point the same way, the ratio will be one. In between, the ratio will increase as the peak gets narrower.

D. Summary

The computation of the ßow direction and local coherence can be summa:rized as follows . First the image I(x, y) is

.....

------------

Figure 2: An image of wood grain with its ßow field . Estimated ßow directions a:re given by the black needles. Thelength of the needle encodes coherence. Notice that coherence is low within the knot at the center.

convolved with the isotropic portion H(x, y) of the filter re­ sponse. The result C(x, y) is then differentiated (by finite dif­ ferences) to form C,(:z:,y) and Cv(:z:,y) . The resulting vec­ tors (C,(x,y),Cv(x,y)) a:re rotated by computing JI(x,y) = 2C,(x,y)Cv(x,y) and J2(x,y) = c;(x,y)- c;(x,y). The gra­ dient magnitude J3(x,y) = [c;(x,y) + c;(x,y)] 1 1 2 also has to be computed in order to measure the coherence. The next step is to convolve J1(:z:,y),J2(:z:,y), and J3(x,y) with the weighting function W(:z:,y) to obtain J{(x,y),J2(x,y), and J;(x,y) . The angle ,P(x, y) of anisotropy and the coherence x(x, y) can then be computed from the formulas

and

An example of this computation applied to a picture of a piece ofwood is shown in figure 2. The ßow direction ,P(x,y) +1r / 2 is displayed by the orientation of small needles superimposed on the image. The lengths of the needles is proportional to the coherence x(x, y) . Note that the pattern is strongly oriented except near the knot in the middle.

E. Relation To Prior Work

The flow computation just described bears an interesting rela­ tion to an early proposal of David Marr that information about local distributions of oriented edge elements be included in the primal sketch (Marr 76). If this proposal is combined with bis later work with Hildreth on edge detection (Marr & Hildreth, 1980) it results in a special case of the above computation. Marr and Hildreth define edges as zero-crossings in the Laplacian of the Gaussian smoothed image. The natural combination of Marr's proposal with this definition of edge elements calls for examining the local density of zero crossings as a function of orientation. For stationary zero-mean Gaussian processes the square of the ori­ ented zero-erossing density is approximately V(8) (see appendix B). Thus in the special case where the point spread function of the filter is S = (a ; az)V 2 exp(- (z 2 + y 2 )/ 2u 2 ) our computation can be viewed as compu ting the direction of minimal edge density in the Marr-Hildreth theory.

Zucker's work on flow (Zucker 1983) is also related to a spe­ cial case of the above computation. For biological reasons, he prefers to use oriented second derivatives of Gaussians as the ini­ tial filters . These have F(r,8) = r 2 exp(-r 2 / 2)cos 2 (8). lnstead of looking at the variance of the filter outpu ts as the o~ientation is changed, he combines the outputs in a biologically motivated relaxation process. Although quite different in detail, the com­ putation described here has much in common with bis technique.

111 Flow coordinates.

The orientation field is an abstraction from the anisotropic pattern that defines it. We can, for example, get the same spiral field from a pattern composed of bands, irregular streaks, dot pairs, etc. In addition to measuring the orientation field, it is useful to be able to produce a description of the underlying pat­ tern independent of the changing direction of anisotropy. Such a description would make it possible to recognize, for example, that two very different orientation fields are defined by the same kind of bands or streaks.

A powerful way to remove the effects of changing orientation is to literally "straighten" the image, .subjecting it to a deformation that maps the flow lines into straight, parallel lines in a canon­ ical ( e.g. horizontal) orientation. Performing this deformation is equivalent to viewing the image in a coordinate system ( u, t1), with u = u(z,y) and t1 = v(z,y) that everywhere satisfies

Equation 6 does not determine a unique coordinate system. An additional constraint may be imposed by choosing lines of constant t1 orthogonal to those of constant u, i.e,

which has the desirable effect of avoiding the introduction of spurious shear in the deformation.

Even with equation 7, an additional constraint is needed, be­ cause we are free to specify arbitrary scaling functions for the u and t1 axes. In the spirit of equation 7, we want to choose these functions to avoid the introduction of spurious stretch or dila­ tion. Although difficult to do globally ( one might minize total stretch,) we will usually want to construct a fairly local coordi­ nate frame araund some point of interest . For this purpose, it suffices to take that point as the origin, scaling the axes u = 0 and t1 = 0 to preserve arc-length along them.

Intuitively, the flow field describes the way the pattern is bent, and viewing the image in these fiow coordinates straightens the pattern out. Figure 3 shows the flow coordinate grid for the wood-grain image from figure 1. The grid lines were computed by taking steps of fixed length in the direction ( cost/>, sint/>) or ( -sint/>, cost/>) for lines across and along the direction of flow re­ spectively, using bilinear interpolation on the orientation field. Since t/> is always computed between 0 and 1r, we must assume that there are no spurious discontinuities in direction to track smoothly.

For many purposes it is unnecessary to compute the deformed image explicitly, but doing so vividly illustrates the ßow Coor­ dinates' ability to simplify the pattern. Figure 4 shows the de­ formation from image Coordinates to ßow coordinates in several stages. As the grain lines straighten the knot shrinks and finally vanishes. The deformed images were anti-aEased using texture­ map techniques (Williams, 1983) The deformed image shows, to a reasonable approximation, what the grain would have looked like bad it not been subjected to the deforming inßuence of the knot.

Thus far, we have separated the image into a flow field, and a pattern derived by viewing the image in ßow coordinates. We argued earlier that the advantage of this decomposition, like the decomposition of an image into intrinsic images, is that the com­ ponents are liable to be simpler and more closely tied to indepen­ dent parts of the pattern-generating process than is the original image. To exploit the decomposition, we need ways of analyzing, describing, and comparing both the ßow field and the straight· ened pattern. These are difficult problems. In the remainder of this section, we present several examples illustrating the utility of the decomposition.

A. A coordinate frame for edge detection.

Oriented measurements have been widely used in edge detec­ tion. For example Marr & Poggio (1979) employed directional second derivative operators, whose zero-crossings were taken to denote rapid intensity changes. Due to the difficulty in select­ ing an orientation, Marr & Hildreth (1980) later abandoned this scheme, in favor of zero-crossings of the Laplacian, a non­ directional operator.

The flow field provides two meaningful dirP.ctions- along and across the direction of flow- in which to Iook for edges within an oriented pattern. Zero·crossings in the second directional derivative in the direction oft/> (again·t the grain) should high­ light edges that contribute to defining the flow field, while zerocrossings in the second derivative perpendicular to 4> (with the grain) should highlightanomalaus elements or terminations. The sum of these two derivatives is the Laplacian.

The two directional derivatives of the wood grain image are shown, with the Laplacian, in figure 5. Indeed, the derivative against the grain captures a.ll the elements comprising the grain pattern, while the derivative with the grain does not appear meaningful. The Laplacian confuses these very different signa.ls by adding them together.

The derivative along the grain can also be meaningful, where anomalaus elements are present. In addition to being perceptu­ ally salient, such anomalies are often physically significant, with origins such as cracks, intrusions, or occlusions, that are dis­ tinct from those of the main pattern. In ma.n-made structures, anamolies are often important because they indicate some variety of ßaw.

Figure 6 shows a pattern of aligned elements (straw) with some anomalaus elements. The directional second derivatives along and across the ßow direction are shown, tagether with their sum (which is just the Laplacian.) Differentiating along the grain highlights anomalaus elements, attenuating the rest ( thus finding the "needles" in the haystack.) Differentiating across the grain supresses the anomalies. The Laplacian shows both.

A related demonstration is shown in figure 7, in which the anomalaus elements have actua.lly been removed by directional median filtering in the ßow direction.

B. Singularities.

We have shown several ways in which viewing an oriented pat­ tern in ßow coordinates facilitates analysis and description of the pattern. Describing and anlyzing the ßow field itself is the other side of the coin. The topology of a ßow field, as of any vector field, is determined by the structure of its singularities, those points at which the field vanishes. Identifying and describing singularities is therefore basic to describing the ßow field. The singularities provide the framework araund which metric properties, such as curvature, may be described. Singularities are also perceptually salient (see figure 8.)

A robust basis for identifying singularities is the index or wind­ ing number (Spivak, 1979.) Suppose we follow a closed curve on a vector field . As we traverse the circuit, the vector rotates con­ tinuously, returning to its original orientation when the circuit is completed. The index or winding number of the curve is the number ofrevolutions made by the vector in traversing the curve. The index of a point is the index of a small circle as we shrink it araund the point:

To compute the winding number numerica.lly, we divide the ßow field into suitably sma.ll rectangles, summing the rotation of </> araund each rectangle. As in computing the ßow lines, we assume that </> has no spurious discontinuities. Where the result is non­ zero, the reetangle surrounds a singularity. Figure 9 shows an example of the detection of singularities using winding number, for a fingerprint . We are currently working on classifying the singularities, and using them to describe the topology of the ßow field.

IV Conclusion

We addressed the problern of analyzing oriented patterns by decomposing them into a ßow field, describing the direction of anisotropy, and describing the pattern independent of changing ßow direction.

A specific computation for estimating the ßow direction was proposed. The computation can be viewed as a) finding the di­ rection of maximal variance in the output of a linear filter, b) combining gradient directions locally, or c) finding the direction of maximal edge density. The computation has been applied to a number of natural and man-made patterns with consistent suc­ cess.

The ßow field was then used to form a coordinate system in which to view the pattern. Two orthogonal familes of curves­ along and across the direction of ßow-form the coordinate sys­ tem's parameter lines. Viewing the pattern in these ßow coor­ dinates amounts to deforming the pattern so that the ßow lines become parallel straight lines. This deformation produces a pat­ tern that is simpler, more regular, and therefore more amenable to analysis and description than the original one.

Several examples of the use of this decomposition were pre­ sented. These included the use of the ßow Coordinates to provide preferred directions for edge detection, detection of anomalies, fitting simple models to the straightened pattern, and detecting singularities in the ßow field.

Our ongoing work focuses on the analysis of patterns with multiple axes of anisotropy, statistical modeling and resynthesis of straightened patterns, and richer description of the structure of the ßow field.

Appendix A: Derivation of Equation 5

The integrals in equation 4 can be evaluated fairly easily by expansion using equation 3 for V . The numerator of tan(2</>) can be written

The first and third integrals are zero and the middle integral is 1r / 8. Hence

·~

Figure 7: Left: the straw picture from figure 6. Right: the anomalous elements have been removed by directional median filtering in the flow direction. (Following a suggestion by Richard Szeliski.)

.. · ..

The denomenator of tan(24>) similarly is

Here the first integral is 1r I 4, the second integral is zero and the third is - w: I 4. Hence the above expression can be simplified to

Substituting for both the numerator and denomenator of tan(24>) in equation 4 yields

There is a close connection between the fiow computation described in section 2 and the density of edges in the Marr­ Hildreth theory of edge detection. Suppose the image /(x,y) is a stationary zero-mean Gaussian process and Iet D(x, y) be the Laplacian of a Gaussian. The density of zero-crossings in Bo(x, y) = D * Ro[I] along the x-axis can be approximated by (Papoulis, 1965; Rice, 1944-45)

where K ( w) is the power spectrum of a slice through B8 along the x-axis. Using Papoulis' Fourier transform conventions, the numerator can be converted to a spatial integral through the identity

Similarly, the denomenator can be converted with the identity

H the integrals are computed locally with the windowing function W, then we have the following estimate for zero crossing density:

H W is radially symmetric, W * B 2 will not depend on 8 so the maximum zero-crossing density will occur at the maximum of W * B 2 . By assumption, the mean of the process is zero, so W * B! is the variance of Dx * Ro[IJ. Thus, for stationary zero­ mean Gaussian processes, if S = D,. A 2 (8) ~ V(O).