Subspace-Based Line Detection

The invention relates generally to signal procssing for line detection, and more particularly the invention relates to fitting of multiple lines in a two-dimensional image.

The background of the invention is described with reference to prior art publications referenced in Appendix A attached hereto.

The problem of fitting lines to a given set of data points is a well-known scientific and engineering task. In computer vision, for example, a digitized image may contain a number of discrete, `1` pixels lying on a `0` background, and the goal is the detection and parameter extraction of a number of straight lines that fit groups of collinear `1` pixels. A number of standard methods exist for such problems. Ordinary least squares methods, for example, seek to minimize the sum of squared vertical (or horizontal) distances of all points to the desired line. However, the resulting solutions are coordinate-system dependent; moreover in ordinary least squares, all measurement error is attributed to only one coordinate of the data points. To avoid these problems, the total least squares method (TLS) [1, 2] uses a criterion in which the sum of squared normal distances of all data points to the line is minimized. However, there are two major limitations that are inherited in all least squares approaches. One is the sensitivity of these methods to outliers. In other words, although the least squares methods can efficiently handle measurement errors that appear as small deviations in the point set, they cannot tolerate outliers; the existence of even a relatively small number of outliers tends to produce large residuals, leading to inadequate and unacceptable line fitting. The second main limitation of least squares methods lies in their inability to handle cases where there is more than one underlying line, especially crossing lines, in an image. The least squares method will result in an obviously wrong answer in these cases; see FIG. 1 for a simple example.

Another now classical method of finding line parameters to fit a given data point set, which is capable of handling multi-line problems, is the Hough transform [3, 4, 5], which is a special case of the Radon transform [6]. In this approach, first a special transformation is applied to all the points and then a two dimensional search is accomplished to find the maxima in the transform plane. The line parameters are .rho., the normal distance of the desired line from the origin, and .theta., the angle that the normal to the line makes with positive x-axis; .rho. and .theta. are also the coordinates in the Hough transform output plane. Each data point (x,y) is mapped under the Hough transform to the curve

in the .rho.-.theta. plane. This equation also represents the line in the x-y plane that has a distance .rho. to the origin and the normal to it makes an angle .theta. with the x axis. Therefore, all the points in the x-y plane located on the line .rho..sub.0 =x cos .theta..sub.0 +y sin .theta..sub.0 are mapped to curves in the .rho.-.theta. plane that all pass through the point (.rho..sub.0,.theta..sub.0).

To fit a straight line to a set of data points, both .rho. and .theta. axes have to be quantized and hence a two-dimensional accumulator array be constructed in the .rho.-.theta. plane. The Hough transform equation is applied to each point in the data set and the contents of all the cells in the transform plane that the corresponding curve passes through are increased by one. This is done until all the data points are transformed and then, an exhaustive search will be needed to locate a number of maxima points in the .rho.-.theta. plane.

The Hough transform method is capable of handling a fairly high amount of noise contamination; however, it has certain drawbacks that may drastically limit its use. One problem, mentioned earlier, is the quantization of the .rho. and .theta. axes: to obtain an acceptable resolution, these coordinates need to be quantized finely, but reducing the quantization steps will broaden the peaks in the transform plane; in other words, since the actual measurement points for a line are not exactly collinear, the corresponding curves in the .rho.-.theta. plane do not pass through the same point. Certain procedures have been investigated to compensate for this difficulty in the quantization step of the Hough transform [7, 8, 9]. Secondly the burden of the two dimensional exhaustive search in the .rho.-.theta. plane for finding the maxima points is a serious drawback of the Hough transform, and this becomes more severe when the number of lines to be fitted to a given data set is unknown. In such cases, the existence of local maxima may cause the search to provide incorrect answers.

The present invention is directed to a new systematic multi-line fitting and straight edge detection procedure.

The Subspace-Based Line Detection (SLIDE) algorithm in accordance with the invention, reformulates the line parameter estimation problem into a subspace fitting framework by applying a propagation (phase-preserving projection) scheme on the image pixels. More specifically, SLIDE is based on the introduction of a perfect mathematical analogy between the problem of estimating angles of multiple lines in an image and the problem of determining the directions of arrival of planar electromagnetic waves impinging on an antenna array. Another analogy is made between the problem of estimating the offsets of the lines and the problem of time series harmonic retrieval. The direction-of-arrival estimation and the harmonic retrieval problems have been extensively studied, e.g. for radar and sonar applications and time series analysis, and in particular, in the last decade several so called super-resolution algorithms (see [10],[11],[12],[13]) have been developed. In the line fitting problem, it is shown that the conditions for using a certain computationally efficient (ESPRIT) algorithm are met. We will introduce these analogies by adopting two signal generation scenarios from the image pixels that we call the constant-.mu. propagation for line angle estimation, and the variable-.tau. propagation for line offset estimation. Another particularly appealing benefit of the analogy is that an efficient statistical technique from the field of sensor array processing can be used to provide an objective estimate of the number of lines in an image.

SLIDE yields closed-form and high resolution estimates of the line parameters, as opposed to the standard method of Hough transform in which the line parameters are subject to limited resolution due to the chosen bin sizes and are found by implementing a search procedure. Computational complexity of SLIDE is an order of magnitude less than that of the Hough transform method. Also SLIDE does not require a large memory space, whereas the Hough transform method needs a large space to store the output grid values.

In the new approach, estimation of line parameters is accomplished in two phases. First, the directions (or angles) of the lines are estimated and then estimates of the offsets are obtained. There exist important applications where only the line orientations are of concern. For example, when an autonomous vehicle is programed to follow a road, the scene that its camera pictures includes primarily (after preprocessing,) some `1` pixels representing the two sides of the road, as in FIG. 2(a), for example. By determining the orientations of these two lines, the vehicle can decide on its next movement. As another application, one can mention the task of precise rotational alignment of wafer and mask in microlithography for VLSI fabrication. In fact, it was the latter problem that inspired the present work. An example of a preprocessed picture from lithography is shown in FIG. 2(b).

We can now introduce a new formalism for the multiple line angle estimation problem that casts this problem into the general framework of array processing problems.

FIG. 1 illustrates an example of line fitting where the least squares method results in a wrong answer.

FIGS. 2(a) and 2(b) illustrate two examples of practical use in which only the orientation of the lines is of concern; FIG. 2(a) is a typical scene pictured by an autonomous vehicle commanded to move along a road where by estimating the direction of the road sides, the vehicle decides on its next movement; and FIG. 2(b) illustrates, in precision wafer-mask rotational alignment, the direction of the feature-lines on the wafer is to be estimated.

FIG. 3 illustrates passive sensor array geometry.

FIG. 4(a) illustrates the image matrix and the hypothetical sensors located along the y-axis, and FIG. 4(b) illustrates measured signal at the l.sup.th receiver due to a single line.

FIG. 5 illustrates an overlapped subarray configuration in a uniform linear array, the subarray separation .DELTA., equal to 1 times the element spacing.

FIGS. 6(a)-6(c) are illustrations of measurement vectors in the line angle and offset estimation problems: (a) image matrix, (b) projection vector obtained in the constant-.mu. propagation for angle estimation, (c) projection vector obtained in variable-r propagation for offset estimation, (top) before, and (bottom) after dechirping. Estimation results: .theta.=-39.9905.degree., x0=100.9645 pixels.

FIGS. 7(a)-7(d) illustrate (a) image matrix, (b) the minimum description length function, (c) estimated line angles via different array procssing techniques, and (d) regenerated image using the estimated line parameters by SLIDE.

FIGS. 8(a) and 8(b) illustrate another example of line fitting by the proposed method, (a) original image with parallel and crossing lines, and (b) reconstructed image using the estimated line parameters: .theta.=30.0636.degree., -60,0041.degree., x0=98.5226, 199.6088, 300.1067, 400,0804, 500.7668, and 100.7426 pixels.

FIGS. 9(a)-9(f) illustrate (a) an example of a noisy image, (b) Hough transform output plane, (c)-(f) comparison of the results obtained by SLIDE and the Hough transform method, with dotted lines the results of SLIDE and solid curves are cross sections of the Hough transform plane where maxima occur, (c), (d) angle estimates, (e), (f) offset estimates.

FIGS. 10(a)-10(d) illustrate (a) original image showing part of a semiconductor wafer, (b) edge-enhanced image, and estimates of the lines angle (c), and offsets (d) by SLIDE and the Hough transform method.

FIGS. 11(a)-11(c) illustrate (a) aerial image, (b) high-contrast image, (c) detected straight roads superimposed on the contour plot of the aerial image. The angle estimates are -21.2685.degree. and 24.7177.degree. relative to the vertical axis. The offset estimates are 48.5027 and 379.6986 pixels, respectively.

FIGS. 12(a)-12(c) illustrate an example of detecting linearly aligned objects: (a) aerial airport image, (b) high-contrast image, (c) estimated support line of the airplanes superimposed on the contour plot of the aerial airport image. The angle estimate is 12.6081.degree. relative to the vertical axis. The offset estimate is 123.8449 pixels.

FIG. 13 illustrates corners of the rectangle which are detected by fitting lines to its sides.

SLIDE reformulates the multiple line angle estimation into the framework of sensor array processing by the virtue of introducing a propagation scenario for the image under consideration. In the so-called constant-.mu. propagation scenario [14], each pixel of the image in a row contributes to a received signal at a hypothetical sensor located in front of the corresponding row of the image matrix. This contribution is in effect equivalent to a phase delay if the pixel is regarded to be a source of electromagnetic wave. More formally, let's assume that all pixels in the image start to propagate narrowband electromagnetic waves with zero initial phases. Furthermore, assume that the waves estimated from pixels in a given row of the image matrix are confined to travel only along that row towards the corresponding sensor. In this propagation environment, each straight line in the image will be in effect equivalent to a wavefront of a travelling planar wave, because the waves of pixels on a line retain constant relative phase along paths parallel to that line. As we will see later, this propagation scenario creates measurements at the sensors similar to the measurements obtained in real antenna array processing.

We assume that a uniform linear array of sensors is located along the vertical axis of the image matrix (4(a)), each sensor receiving signals only from its corresponding row in the matrix according to the following rule. If there are p nonzero pixels on the l-th row of the image matrix located on columns q.sub.1, . . . q.sub.p respectively, then the sensor in front of the l-th row will receive the signal ##EQU1## where .mu. is a constant parameter to be chosen at our convenience [15]. Now assume that there is one straight line in an N.times.N image with offset x.sub.0 and angle .theta., as in FIG. 4. Provided that the line width is such that the line gives rise to only one nonzero pixel per row, and that there are no noise pixels, the measured signal at the l-th receiver will be

In a more general case, where there are d straight lines and also noise pixels in the image, the signal will be ##EQU2## where n.sub.l includes the effect of noise in the l-th row, which may consist of individual noise pixels and also takes care of a possible displacement of a line pixel in that row.

Now let us define ##EQU3##

We will call this signal generation scenario the constant-.mu. propagation. This scenario is effectively encoding the line angles as frequencies of cisoidal components in the measurement vector z=[z.sub.0, . . . , z.sub.N-1 ].sup.T.

As Eq. (3) shows, the measurement vector z includes d cisoidal components contaminated with noise. Any spectral estimation technique (including FET) can be applied in this stage to extract the frequencies of the cisoids. However, since high resolution estimates of the angles are desired, we pursue with exploiting the subspace-based techniques of array processing framework.

Using a window of length M, we can arrange the measurements in a matrix form as ##EQU4## where P=N+1-M and the new snapshot vector z.sub.p .epsilon. C.sup.M.times.1 can be written as

In the above, A(.theta.)=[a(.theta..sub.1), . . . , a(.theta..sub.d)] is an M.times.d matrix with k-th column a(.theta..sub.k)=[1..lambda..sub.k, . . . , .lambda..sub.k.sup.M-1 ].sup.T with .lambda..sub.k =e.sup.j.mu. tan .theta..sbsp.k. This rearrangement of the data is called spatial smoothing.

Combining Eqs. (5) and (6), we can write ##EQU5## contains the offsets of the lines, and .PHI.=diag[.lambda..sub.1, . . . , .lambda..sub.d ].

Now we are ready to apply the well-studied array processing techniques to the multi-line angle estimation problem. In sensor array processing framework, the goal is to estimate the directions of arrival .theta. of d impinging planar wavefronts on the array of M sensors, given P snapshots z.sub.i, i=1, . . . , P. There are many algorithms for solving this problem (see e.g. [10]-[13]), but for the present application a particular one called TLS-ESPRIT appears to be the most appropriate.

In the ideal case of P.fwdarw..infin., the covariance matrix of the measurements has the following structure ##EQU6## and it is assumed that ##EQU7## Now let the eigendecomposition of R.sub.zz be ##EQU8## where .lambda..sub.1 .gtoreq..lambda..sub.2 .gtoreq. . . . .gtoreq..lambda..sub.M. The matrix E.sub.s contains the d eigenvectors of R.sub.zz corresponding to its d largest eigenvalues, i.e. E.sub.s =[e.sub.1, . . . , e.sub.d ]. The range space of E.sub.s is called the signal subspace. Its orthogonal complement is called the noise subspace and is spanned by the columns of E.sub.n =[e.sub.d+1, . . . , e.sub.M ]. It is easily seen from (11) that the smallest eigenvalue of R.sub.zz has multiplicity M-d and is equal to the noise variance. Thus, .LAMBDA..sub.n =.sigma..sup.2 I. Moreover, it can be seen that the signal subspace is equal to span {A(.theta.)}, i.e. span {E.sub.s }=span {A(.theta.)}, thus explaining the name. This is the basic observation in all subspace-based array processing techniques, including ESPRIT.

In practice, only the same covariance matrix R.sub.zz.sbsb.P is available, and its eigendecomposition is ##EQU9## where now .LAMBDA..sub.n .noteq..sigma..sup.2 I, except asymptotically as P.fwdarw..infin..

Determining the value of d is crucial here because by knowing d we can identify the eigenvectors corresponding to the d largest eigenvalues of R.sub.zz as a basis for the estimate of the signal subspace. Under the stochastic signal assumption, the minimum description length (MDL) criterion is defined by ##EQU10## where .lambda..sub.i are the eigenvalues of R.sub.zz. The use of this criterion was first proposed in [16]. The estimate of the number of sources d is given by the value of k for which the MDL function is minimized.

Now define E.sub.1 as the submatrix E.sub.s formed by rows 1 to M-1, and E.sub.2 as the submatrix of rows 2 to M. This is a special case of partitioning the receiver array into two maximally overlapping subarrays with a displacement .DELTA.; here we have chosen .DELTA.=1 (see FIG. 5). Next form the eigendecomposition of the 2d.times.2d matrix ##EQU11## Partition F into d.times.d submatrices, ##EQU12## and then let .lambda..sub.k, k=1, . . . , d be the eigenvalues of -F.sub.12 F.sub.22.sup.-1. Finally the line angles are obtained as ##EQU13##

In this section, we briefly discuss the method used in SLIDE algorithm for estimating the line offsets.

After the angles have been estimated by the method of Sec. 3.1, the offsets can be estimated by either of the three methods described in this section.

As was mentioned above, the vector ##EQU14## contains the offsets of the lines. From the analogy with sensor array processing framework, solving the so-called signal copy problem yields estimates of values of elements of s. This is done by solving the following problem

In a noise-free case, this criterion tries to find appropriate phases for cisoidal signals in order to fit their superposition to the given measurement vector y [17]. It is easy to show that the distribution of the measurement noise is gaussian if the outliers are uniformly distributed over the image [15]. Therefore, the above least squares problem provides the maximum likelihood estimate for the signal vector. Then, the offsets are computed via the relation ##EQU15##

However, the least squares technique is highly sensitive to outliers and its result rapidly degrades if the number of outliers increases. Also in many practical situations, the image may contain only a part of a straight line. In such cases, the above least squares estimate will be biased and unsatisfactory. There also exists a degenerate case in using this technique for offset estimation. Namely, since each straight line is assumed to represent a travelling planar wavefront, it can be shown that the phases of the signals originating from parallel lines cannot be resolved at the sensors. In other words, the received signal will be a coherent superposition of the two or more parallel wavefronts.

Another method for estimating the offsets is to project the image along each of the detected angles and search for a peak in the resulting projection. This method is robust to outliers but suffers from limited resolution problem as we have to choose bin sizes for the projection.

In order to obtain closed-form and high resolution estimates for the offsets, SLIDE proposes another propagation scenario that leads to encoding of the line offsets as frequencies of cisoidal components in the measurement vector. We will call this propagation scheme the variable-.tau. propagation for reasons that will become clear shortly. The propagation part of this procedure can be done in parallel with the constant-.mu. propagation that encodes the angles, but the offset estimation procedure that follows the new propagation scenario has to be done after the line angles have been estimated.

Assume that there are d.sub.1 parallel lines in the image with angle .theta..sub.1. Performing the same propagation scheme of above to the pixels on the l-th row and replacing the name of propagation parameter .mu. with .tau., the signal received at the l-th sensor is ##EQU16## Now we let .tau. vary along the rows of the matrix. For simplicity of the analysis, here we choose .tau. to be proportional to the row number l

where .alpha. is a constant. Thus, the measurement vector z will now have elements that are functions of .tau. ##EQU17## As we observe in Eq. (24), in the new measurements .tau. appears in quadratic form with the angle and in linear form with the offsets. If there is only one line in the image (d.sub.1 =1), then the sequence z.sub..tau. will be of the form

which has the form of a chirp signal. In general, z.sub..tau. in Eq. (24) represents the sum of d.sub.1 chirp signals with the same quadratic term but with different linear terms. The quadratic term can be factored out of the summation and as before be called a.sub..tau. (.theta..sub.1) as in Eq. (25). Since we know the angle .theta..sub.1 from the previous analysis (constant-.mu. propagation), we can divide z.sub..tau. by a.sub..tau. (.theta..sub.1) and get another sequence w.sub..tau. which is free from the quadratic term in .tau.. This is basically equivalent to dechirping the measurement sequence. The new measurement vector w has elements of the following form ##EQU18## where n'.sub..tau. =n.sub..tau. /a.sub..tau. (.theta..sub.1) is the new additive noise term.

The sequence w.sub..tau. in Eq. (28) is a superposition d.sub.1 cisoids plus noise. The frequency of each of the cisoids is proportional to one of the line offsets. Therefore, we have reformulated the offset estimation problem as a time series harmonic retrieval problem. By exploiting the eigenstructure of the covariance matrix of the measurements and subspace fitting techniques of sensor array processing, high resolution estimates of the offsets can be obtained. Also the number of parallel lines for each angle can be estimated using the minimum description length criterion.

In general, an image may contain parallel and non-parallel lines. In order to analyze the nature of the measurement vector w in that case, let's assume that there are d.sub.1 parallel lines in the image with angle .theta..sub.1 and d.sub.2 parallel lines with angle .theta..sub.2. We order the offsets of all the lines such that ##EQU19## correspond to the lines with angle .theta..sub.1 and ##EQU20## correspond to the lines with angle .theta..sub.2. The z vector will have elements ##EQU21## Now consider the estimation of the offsets of the new lines having angle .theta..sub.1. As before, by dividing the sequence z.sub..tau. by a.sub..tau. (.theta..sub.1), we dechirp the components of the measurement sequence that are contributions from angle .theta..sub.1. This results in ##EQU22## The terms corresponding to angle .theta..sub.1 have been dechirped and transformed to pure cisoids as before. The terms corresponding to angle .theta..sub.2 however, have been transformed to another chirp form. The basic assumption for applying the ESPRIT algorithm to the harmonic retrieval problem is that the components s(t) of the underlying analog signal be narrowband, i.e. for each component we can write s(t-t.sub.0).apprxeq.e.sup.-j.omega.t.sbsp.0 s(t). If only the pure cisoids corresponding to the liens with angle .theta..sub.1 were present, this condition would be directly satisfied. The sequence w.sub..tau..sup.(.theta..sbsp.1.sup.) in Eq. (33) contains both pure cisoids and chirp components. However, the energy of the chirp components is distributed over the frequency axis and they can be considered as noise in the analysis. Therefore, the harmonic retrieval procedure will be able to detect the pure cisoids corresponding to the angle under consideration. More clarification and examples for this point will be presented later. In a general case, the variable-.tau. propagation needs to be performed only once, and the dechirping and harmonic retrieval procedures have to be done for each of the detected angles.

The procedure of above is then used to rearrange the sequence w.sub..tau. into a Hankel matrix W.sub.P as in Eq. (6) and then to implement the subspace-based estimation techniques such as ESPRIT to solve for the frequency components in w.sub..tau. which correspond to the line offsets.

Generalizing to Grey-Scale Images

The formalism introduced above can also be generalized to handle the problem of straight edge detection in grey-scale images. Namely, we can start from a grey-scale image and detect parameters of straight lines in that image. This is accomplished by simply assigning an amplitude to the propagated signal from each pixel proportional to the grey-scale value of that pixel. A primitive edge enhancement procedure may be used first in this case to attenuate background contributions. This edge enhancement can in general be as simple as a shirt-and-subtract or a contrast enhancement procedure applied to the image.

Mathematically, if the edge strength in the first row (after edge enhancement) is represented by a function f(x)=g(x-x.sub.0), then the sensor in front of the first row will receive ##EQU23## where G(.mu.) is equal to the Fourier transform of g(x) at a single frequency point, and is a constant. In binary images we have G(.mu.)=1. The edge strength in row l (in a one-line and noise-free image) will be equal to g(x-x.sub.0 +l tan .theta.), and the receiver in front of that row will read ##EQU24## The only difference between the above equation and Eq. (2) is the appearance of a constant (G(.mu.)) as the amplitude. This treatment shows that identical procedures as those that are proposed for line fitting in binary images can be performed on edge-enhanced grey-scale images.

A case of interest in line fitting in grey-scale images is non-uniformity in edge strength over different image lines. Mathematically, the deviation in value from the nominal (or average) edge strength can be lumped into the noise term for each line. Simulation results show that as long as the signal to noise ratio remains in a reasonable level, the algorithm is still able to detect the edge in these cases.

Aliasing Issues

As was mentioned above, the solution is provided in two steps. First, the line directions are estimated and then an estimate of the offsets is obtained. We start with the estimation of line angles .theta..sub.k, k=1, . . . d. Using the TLS-ESPRIT algorithm, the angles are estimated via the relation ##EQU25## where .DELTA. is the displacement between the two subarrays used in the ESPRIT algorithm. In the above equation, normalization is done in order to compensate for possible perturbations in magnitude of .lambda..sub.k =e.sup.j.mu..DELTA. tan .theta..sbsp.k from unity. The range of non-ambiguous angles is determined by the following formula ##EQU26## Values of .mu. and .DELTA. can be chosen appropriately in order to avoid the ambiguity problem. It is also worth mentioning that since the value of .mu. is not dictated by the physics as it is the case in sensor array processing, improved resolution can be achieved by increasing its value. To be more specific, from Eq. (37) we must choose ##EQU27## in order to avoid ambiguity. In practice however, due to the fact that the image is in fact discretized (sampled), the value of .mu. cannot be chosen higher than .pi.. For angle estimation stage, .mu. is better to be chosen to have a value around unity to provide both non-ambiguity and good resolution. In estimating the offsets, the value of .mu. should be chosen such that the line offsets do not become ambiguous. In other words, since the received signal due to a pixel at a receiver has a period equal to 2.pi., and offsets can have values in the range [0, N-1], .mu. has to be chosen such that .vertline..mu..vertline..ltoreq.2.pi./N. Therefore, it is recommenced that the propagations performed in angle and offset estimation steps be performed with different values of .mu..

The aliasing issue arises also in the line offset estimation problem. Here, the factor that needs to be controlled is the step size chosen for the sequence .tau., which we call .delta..tau.. Larger values of .delta..tau. result in more periods of the cisoids in the data records which is desirable if aliasing does not occur. In general, the value of .delta..tau. should satisfy ##EQU28## in order that aliasing does not happen. However, if an estimate of the maximum possible value for all line offsets is available, the N in the above formula can be replaced by that value and .delta..tau. can be increased accordingly.

Combining Different Estimation Results

Another flexibility of the proposed algorithm that can be exploited to improve the signal to noise ratio of the angle estimation procedure is the possibility of doing estimation using hypothetical sensors along the x-axis as well and then combining the results of the two estimations.

In offset estimation procedure, the mentioned combination of estimates along the x-axis and the y-axis can help to improve the results. Moreover, we can observe that in this case the lines that have larger offsets are transformed to cisoids with higher frequencies. These frequencies can normally be detected more accurately due to the presence of many cycles of the signal in the data record. Therefore, the propagation can be performed along both positive and negative directions of the coordinate axes, and in each case the more accurate results be chosen according to the location of the estimated offsets. For example, with propagating along negative x-axis, offset estimates with values larger than N/2 are more accurate, while with propagating along positive x-axis, the estimates that are less than N/2 are preferred.

Partial Extent Lines

Another practical issue arises when there are line segments in the image rather than long lines. The spatial smoothing procedure provides a level of robustness for the process of angle estimation by sliding a small size window over the measurement vector. This way, a line with a support long enough to be included in many windows will be adequately represented in the covariance matrix. However, small support lines cannot be represented adequately and won't be resolved in the estimation. A remedy for such a case is to divide the image into several smaller sub-images and perform the analysis on each sub-image. The results of individual sub-image estimations can be put together to obtain parameters of both long and short lines in the image. In order to clarify why lines that have partial support in the image can be detected, the following observation is also worth mentioning. The basic assumption for applying the subspace techniques to a set of subsequent snapshots is that there exists a certain phase relation between subsequent snapshots. The matrix .PHI. in the ESPRIT algorithm is in fact a diagonal matrix containing these phase relations for different signal components. In the SLIDE algorithm, this phase relation is retained between subsequent subvectors made from measurements that result from lines with partial extent. More specifically, snapshots (subvectors) that are fully covered by the lateral extent of the line clearly satisfy the phase relationship. Snapshots that come from the outside of the line support have zero contribution to the estimate (although in the presence of noise they tend to lower the SNR). Finally, for snapshots that are partially covered by the line lateral extent, the phase relation is preserved between each two subsequent snapshot in all but two points that correspond to the ends of the line support. This is the case simply because the same phase relationship that exists for the nonzero part of the snapshots can be assumed to exist for the identically-zero parts of them. Another viewpoint about this observation follows shortly.

Choice of the Window Size

In the procedure of windowing the long N.times.1 measurement vector z (and w) to produce the matrix Z.sub.P (and W.sub.P), the values of P and M need to be chosen such that P plays the role of number of snapshots and M is the dimension of the snapshots. Obviously, for the sample covariance matrix to approximate the true covariance matrix, we need to choose a large P. On the other hand, choosing the dimension M of the new vectors large helps to smooth out the eigenvalues corresponding to the noise subspace of the covariance matrix. Therefore, in order to make the best use of the measurement vector, sub-vectors are chosen in such a way to maintain maximum overlap with their adjacent neighbors. With this choice, there will be a fixed relation between the number of snapshots P, and the dimension of the new measurement vector M, namely P+M=N+1. It should also be mentioned that choosing the dimension of new vectors to be very large increases the dimension of the sample covariance matrix. This results in dramatic increase in computation load if one tries to obtain the whole eigendecomposition of the sample covariance matrix; however, if the so called Fast Subspace Decomposition (FSD) technique of Xu and Kailath [18] which is based on the Lanczos algorithm [19, 1] is used to extract only the signal subspace of the sample covariance matrix, the computation load will increase linearly with M.

In order to compare the computational effort needed in SLIDE with that of the Hough transform method, two cases are considered here. In the first case, it is assumed that the image is binarized and that there is one pixel in each row of the N.times.N image matrix for every one of the d lines. Without the presence of any outliers, the computational effort needed for applying the Hough transformation to all pixels will be of O(N.sup.2 d) provided that we choose to discretize the .theta. axis to N bins; in some applications, we may have prior information about the possible range of the angles, and discretize only that range into N bins. In addition to the above, a search procedure is needed for estimating the location of peaks in the transform plane. On the other hand, SLIDE needs O(Nd) calculations for the propagations and O(N M.sup.2) for the complete eigendecomposition of the sample covariance matrix. However, if the structure of the data matrix is exploited and also only the signal subspace of the sample covariance matrix is extracted as in the new method of Xu and Kailath [18] (note that this subspace contains all needed information for estimating the parameters), then the total computation load of SLIDE will be further decreased to O(N(M+2d)). In practice, the value of M is chosen around .sqroot.N, so the computational complexity of SLIDE in this case is O(N.sup.3/2) compared to O(N.sup.2 d) for the Hough transform method. In the second case, it is assumed that all the N.sup.2 pixels in the image participate in the computations (e.g. we have a grey-scale image). In this case, the Hough transform method requires O(N.sup.3) calculations plus search. The computational effort for SLIDE in this case will be O(N.sup.2) for the propagations and again O(N(M+d)) for the rest of the calculations. So the computational complexity of SLIDE in this case is O(N.sup.2) as compared to O(N.sup.3) for the Hough transform method. This reduction in computational complexity could allow real-time implementation of SLIDE.

In this section, an outline of the SLIDE algorithm as applied for estimating angles and offsets of straight lines and edges in an image is presented. Again, the image size is assumed to be N.times.N.

If the image is binary, go to the next step. If the image is grey-scale, produce an edge-enhanced image from the digitized grey-scale image. This can be done by a simple shift an subtract (one-pixel gradient) operator or a contrast enhancement procedure.

Perform the constant-.mu. propagation as discussed above on the image and obtain the measurement vector z. In edge-enhanced grey-scale images, this can be done by multiplying each row of the image matrix by the vector [1,e.sup.-j.mu., e.sup.-j2.mu., . . . , e.sup.-j(N-1).mu. ].sup.T. However, in binary images savings in computations can be made by performing calculations only on nonzero pixels of the image.

In this stage, any spectral estimation method can be applied to estimate the cisoidal components of z. However, in order to obtain high resolution and closed-form estimates, use subspace-based techniques of array processing as follows.

Rearrange the measurement vector using Eq. (6) such that P snapshot vectors of size M are produced.

Compute the sample covariance matrix of the snapshots as in Eq. (15); and perform its eigendecomposition. By applying the MDL method, the number of present angles can be estimated.

Estimate the angle by any subspace-based method. An outline of how to exploit the ESPRIT algorithm is presented above. ESPRIT provides high resolution estimates without the need for a search procedure in the parameter space.

After knowing the angles of the lines, estimate the offsets by any of the three methods that follow. Using least-squares minimization as in Eq. (20), or projecting the image along each detected angle and searching for peaks, or applying the variable-.tau. propagation stage of SLIDE as described in Sec. 3.2 and summarized below.

For exploiting the offset estimation part of SLIDE, perform the variable-.tau. propagation (as described in Sec. 3.2), by multiplying each line of the image matrix by a vector as follows. Numbering the rows of the image matrix from 0 to N-1, multiply the l-th row by the vector [1, e.sup.-jl/.alpha., e.sup.-j2l/.alpha., . . . , e.sup.-j(N-1)l/.alpha. ].sup.T to obtain measurement vector z. (Again, in binary images computational savings can be made by only performing calculations only on nonzero pixels.)

If there are d different angles in the image, perform the following stages for dechirped measurements for each angle obtained from the measurement vector.

Divide each element of the measurement vector by the corresponding element response value as in Eq. (28) to get the dechirped measurement vector w.

Rearrange the elements of w the same way the measurement vector z was rearranged in angle estimation stage to get P snapshots each of length M.

Compute the sample covariance matrix of the new snapshots and perform its eigendecomposition. Use the a priori knowledge or MDL to decide on the number of present parallel lines.

Exploit the ESPRIT algorithm (or other spectral estimation methods) to get estimates of the offsets.

In this section, several experimental results of applying the SLIDE algorithm to estimating the parameters of lines in images are provided. In order to make comparisons, the method of Hough transform which is essentially a projection technique is also applied to images and the results of both approaches are presented. It is worth re-emphasizing here that the estimates provided by the SLIDE algorithm are high resolution while the Hough transform (and projection) results are limited in resolution to the bin sizes chosen for the parameters. Also, SLIDE is computationally much more efficient than the Hough transform and there is no need for a search procedure in the parameter space or curve reading in its implementation.

FIG. 6 illustrates the measurement vectors obtained in constant-.mu. propagation and variable-.tau. propagation. FIG. 6(a) shows an image containing a straight line and outliers. In FIG. 6(b) the measurement vector obtained in constant-.mu. propagation is shown. In FIG. 6(c), the measurement vector obtained in variable-.tau. propagation is shown (top) before, and (bottom) after dechirping. The line angle and offset estimation problems are transformed to the problem of estimating frequencies of cisoids contaminated in noise.

In FIG. 7(a), an image containing five nearly collinear pixel sets is shown. The line angles are {-40, -35, 25, 30, 50} degrees and the line offsets are [25, 45, 150, 175, 190} pixels, respectively. After performing the constant-.mu. propagation and computing the sample covariance matrix of the measurements, the minimum description length function is calculated. This function is shown in FIG. 7(b); the location of its minimum is an estimate of the number of lines in the image. In FIG. 7(c), the results of applying different sensor array processing techniques to this example are presented; the ESPRIT algorithm requires no search procedure and is by far the fastest one to use. The estimated angles by SLIDE are {-40.0424, -35.0399, 25.2412, 30.2075, 50.0380} degrees and the estimated offsets are {25.7130, 44.9882, 151.1676, 174.0602, 190.0925} pixels, respectively. Using these values, the image of FIG. 7(d), is regenerated.

Another example is presented in FIG. 8 in which parameters of the parallel and crossing lines in an image are estimated.

To compare the efficiency and speed of the proposed technique with the standard method of Hough transform, we consider the rather noisy image of FIG. 9(a). This image is 400.times.400 pixels and was generated to contain two lines with random angles .theta.=19.8418.degree., 27.2582.degree. and random offsets x.sub.0 =186.8439, 277.2624 pixels, respectively. Note that due to discretization effects, the actual values may be slightly different from the above. For applying the Hough transform, we need to choose a search range for the angle and a bin size for discretizing the normal distance (.rho.) of lines from the origin. For this example, a uniform range of angles was chosen between -35.degree. and 35.degree. with 0.1.degree. steps. The bin size for .rho. was chosen to be 0.5 pixels. Note that .rho. is the normal distance of lines from the origin, and the corresponding x-axis offset is computed by ##EQU29## FIG. 9(b) shows the Hough transform output plane, where the coordinates of the peaks yield estimates for .theta. and .rho.. Cross sections of the Hough transform plane at the locations of the peaks are depicted in FIG. 9(c)-(f) together with estimates obtained by SLIDE. The numerical values are as follows. The estimates by SLIDE are .theta.=-19.7485.degree., 27.2491.degree. and x.sub.0 =187.1408, 277.6463 pixels. Performing a search in the Hough transform output plane, yields the angle estimates as .theta.=-19.8.degree., 27.2.degree. and the offset estimates as x.sub.0 =187.6, 277.7 pixels. It is seen that both techniques yield rather accurate estimates, however, the effort of computations of SLIDE is about an order of magnitude (400 times in this example) less than that of the Hough transform. Also, the Hough transform method needs a large memory space to store and accumulate the output.

As is discussed above, a straight-forward generalization of the proposed algorithm for being applied to grey-scale images can be readily derived. Consider the grey-scale image shown in FIG. 10(a). First, a primitive edge enhancement was done on this image in order to attenuate the background contributions. This was done by shifting the image matrix by one pixel, subtracting the result from the original matrix, and taking the absolute value of the resulting matrix. The image of FIG. 10(b) shows the result of this edge enhancement procedure. Then, SLIDE and the Hough transform method were applied to this image. FIGS. 10(c), (d) show the angle and offset estimation results, respectively.

In FIG. 11, the result of applying SLIDE to a real aerial image is presented. FIG. 11(a) shows the original aerial image. The high-contrast image of FIG. 11(b) is used for line detection. In FIG. 11(c), the results of line detection are superimposed on the contour plot of the image.

FIG. 12 is an example of detecting linearly aligned objects by SLIDE. FIGS. 12(a), (b) show the original and high-contrast aerial images in which several airplanes are aligned along a certain direction which we might like to determine. This is an illustrative example of non-uniform edge strength over different image rows. Applying the SLIDE algorithm to the high-contrast image results in line parameters that are used to generate the dotted line in FIG. 12(c). As this figure suggests, the direction along which the airplanes lie seems to be well estimated by the algorithm.

In FIG. 13, SLIDE is applied to detection of the sides (and hence the corners) of a rectangle. Although the sides of the rectangle have a limited extent in the image (and in the measurement vector), the estimates of the corners are reasonably accurate.

While the invention has been described with reference to specific embodiments, the description is illustrative of the invention and is not to be construed as limiting the invention. Various modifications and applications may occur to those skilled in the art without departing from the true spirit and scope of the invention as defined by the appended claims.

[1] G. H. Golub and C. F. Van Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, Md., 1984.

[2] S. Van Huffel and J. Vandewalle, "The Total Least Squares Technique: Computation, Properties and Applications", In E. F. Deprettere, editor, SVD and Signal processing: Algorithms, Applications and Architectures, pages 189-207. Elsevier Science Publishers, 1988.

[3] P. Hough, "Method and Means for Recognizing Complex Patterns", U.S. Pat. No. 3,069,654, 1962.

[4] R. O. Duda and P. E. Hart, "Use of the Hough Transform to Detect Lines and Curves in Pictures", Communications of the ACM, 15:11-15, 1972.

[5] A. K. Jain, Fundamentals of Digital Image Processing, Prentice-Hall, Englewood Cliffs, N.J., 1989.

[6] S. R. Deans, The Radon Transform and Some of Its Applications, John Wiley & Sons, New York, N.Y., 1983.

[7] W. Niblack and D. Petkovic, "On Improving the Accuracy of the Hough Transform: Theory, Simulations and Experiments", In Proc. IEEE Conf. on Computer Vision and Pattern Recognition CVPR'88, pages 574-579, 1988.

[8] E. R. Davies, Machine Vision: Theory, Algorithms, Practicalities, Academic Press, London, 1990.

[9] N. Kiryati and A. M. Bruckstein, "On Navigation Between Friends and Foes", IEEE Trans. on Pattern Analysis and Machine Intelligence, 13(6):602-606, 1991.

[10] R. O. Schmidt, "Multiple Emitter Location and Signal Parameter Estimation", In RADC Spectrum Estimation Workshop, Griffiss AFB, NY, 1979.

[11] A. Paulraj, R. Roy, and T. Kailath, "Estimation of Signal Parameters by Rotational Invariance Techniques (ESPRIT)", In Proc. of 19th Asilomar Conference on Circuits, Systems and Comp., 1985.

[12] M. Viberg, B. Ottersten, and T. Kailath, "Detection and Estimation in Sensor Arrays Using Weighted Subspace Fitting", IEEE Trans. on SP, 39(11):2436-2449, Nov. 1991.

[13] R. Roy and T. Kailath, "ESPRIT: Estimation of Signal Parameters via Rotational Invariance Techniques", IEEE Trans. on ASSP, 37(7):984-995, July 1989.

[14] H. K. Aghajan and T. Kailath, "SLIDE: Subspace-based Line Detection", In Proc. of IEEE ICASSP, Minneapolis, Minn., to appear, April 93.

[15] H. K. Aghajan and T. Kailath, "Sensor Array Processing Techniques for Super Resolution Multi-Line Fitting and Straight Edge Detection", IEEE Trans. on Image Processing, to appear, 93.

[16] M. Wax and T. Kailath, "Detection of Signals by Information Theoretic Criteria", IEEE Trans. on ASSP, 33(2):387-392, April 1985.

[17] H. K. Aghajan and T. Kailath, "A Subspace Fitting Approach to Super Resolution Multi-Line Fitting and Straight Edge Detection", In Proc. of IEEE ICASSP, pages III:121-124, San Francisco, Calif., 1992.

[18] G. Xu and T. Kailath, "A Fast Algorithm for Signal Subspace Decomposition and Its Performance Analysis", In Proc. of IEEE ICASSP, pages 3069-3072, Toronto, Canada, 1991.

[19] B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice-Hall, Englewood Cliffs, N.J., 1980.