Multiscale Feature Detector Using Filter Banks

The field of the invention is the detection of features in images, and particularly, the detection of features such as edges in two-dimensional images.

There are numerous applications in which images are analyzed to locate specific features. For example, an object may be scanned with a camera and the edges of the object may be detected and used to locate the object and determine its orientation. Edge detection is also used in establishing the boundaries of an object so that it may be identified using pattern recognition methods, or measured. For example, an ultrasound image of the human heart may be analyzed by first locating the edges of the various heart chambers and then calculating the sizes of those chambers and the thicknesses of the chamber walls.

There are many approaches used to detect features in images. One approach is to fit a model of the feature to the image. The model is moved around the image and where there is a good match the feature location is indicated. Other methods use probabilistic measures, and most methods use filtering techniques which attempt to remove noise without obscuring the feature being detected. Filter operators such as the Roberts operator, Haar operator, Sobel operator, Canny's detector and wavelet transform are used to enhance images as a first step in feature detection, but the effectiveness of such filtering is determined by the choice of the filter scale. The choice of the filter scale may be arbitrary, or it may be based on results from prior applications of a similar nature. If the filter scale is large, filtering is intense and noise is smoothed out at the expense of lost localization of a detected feature or even failure to detect a feature. If the scale is small, on the other hand, less filtering is used and the results may be corrupted by noise which produces false alarms.

The present invention is a method and apparatus for detecting features in an acquired image, and particularly, it includes: a multiscale filter bank which applies a filter operator to the image data to produce a plurality of filtered outputs at a corresponding plurality of scales; a detector for each filtered output which indicates the locations of a selected feature therein, and a MAP correlator which correlates the filtered output values at each indicated feature location with a model response and produces an output signal indicative of the presence of the feature at the location if the correlation exceeds a preset threshold.

A general object of the invention is to provide a feature detection system which is robust. Nearly all state-of-the art feature detectors work satisfactorily when the image is not corrupted with noise.

However, using the multiscale filter bank and the map correlator of the present invention enables the same features to be accurately located even when there is substantial noise in the image. Noise is a very significant factor, for example, with medical images and the present invention is particularly useful in such applications.

Another object of the present invention is to match the MAP correlator with the multiscale filter so that optimal performance is achieved. The model response and the preset threshold used by the MAP correlator to determine the presence of a feature are precisely determined by the operating parameters of the chosen multiscale filter bank.

The foregoing and other objects and advantages of the invention will appear from the following description. In the description, reference is made to the accompanying drawings which form a part hereof, and in which there is shown by way of illustration a preferred embodiment of the invention. Such embodiment does not necessarily represent the full scope of the invention, however, and reference is made therefore to the claims herein for interpreting the scope of the invention.

FIG. 1 is diagram of the elements of a feature detection system;

FIG. 2 is a graphic representation of discrete multiscale filters at different sampling rates and their corresponding frequency transforms derived from the Haar wavelet;

FIG. 3 is a diagram of a multiscale filter;

FIG. 4 is a schematic representation of the relationship between a multiscale filter and a corresponding multiscale filter bank;

FIG. 5 is a graphic representation of a sampling lattice produced by a decimator D used in a 2D filter bank of the present invention;

FIG. 6 is an electrical block diagram of the preferred embodiment of the invented feature detection system;

FIG. 7 is a schematic representation of an indicator response table which is a data structure employed in the system of FIG. 6;

FIG. 8 is a schematic representation of a scale space tree table which is a data structure employed in the system of FIG. 6;

FIG. 9 is a flow chart of the functions performed by a feature generator which forms part of the system of FIG. 6; and

FIG. 10 is a flow chart of the functions performed by the MAP correlation which forms part of the system of FIG. 6.

The present invention is a system for detecting and locating features in an input signal. Referring to FIG. 1 it includes a multiscale filter 10 which receives the input signal 12 and filters the signal to produce a set of filtered output signals 14a-c. The filtered output signals 14a-c represent the same input signal 12, but they are filtered at different "scales". The lower scale signal 14a is the highest resolution signal and it has undergone the least filtering. The highest scale signal 14c has been filtered the most and it is the lowest resolution signal. While 1D signals are shown in FIG. 1, 2D signals such as images and 3D signals such as image sequences (video) may also be processed.

The multiscale filtered signals 14a-c are applied to a multiscale detector 16. The detector 16 provides an indication of the location in the filtered signals 14a-c of the desired feature. "Features" to be detected may include, for example, step edges, roof edges, and pulse edges, and these may be "indicated" by maxima, zero-crossing, or a constant interval in the filtered signals. For example, maxima and zero-crossings are two common indicators for detecting edges and the choice of indicator will depend on the feature to be located and the type of filter used. For example, when the feature is a step edge, a maxima detector may be paired with a filter such as that described by Canny, "first difference operator, first derivative of Gaussian" or a quadratic spline filter. On the other hand, a zero-crossing detector may be paired with a Laplacian of Gaussian (LOG) filter. Regardless of the choice, the multiscale detector 16 produces samples of the filtered signals 14a-c at each indicated location at corresponding outputs 18a-c. The multiscale detector 16 in essence identifies locations in the filtered signals 14a-c where the feature may be present and masks out the remainder of the signal. This step significantly reduces the processing burden on a multiscale MAP correlator 20 which makes the final determination if the feature is present.

The multiscale MAP correlator 20 makes the final decision of where the feature occurs in the signal 12 and produces an output 22. As will be described in more detail below, it makes this determination by correlating the detector output signals 18a-c with a reference vector that represents the expected response when the feature is present. This reference vector as well as the threshold value of the correlation required to indicate the presence of the feature is precisely defined. In a 1D system the output 22 indicates the times at which the feature occurs in the signal 12, whereas in a 2D system the output 22 is a bit map showing where in the 2D input signal 12 the feature is located.

The construction of the multiscale filter for this system must meet a number of design criteria. The highest resolution filter is the smallest filter in the multiscale set, and it is the most critical for localization of the feature. To determine the highest resolution scale, one needs to examine the Fourier transforms of the desired set of filters. The filters in the multiscale set have to preserve their shape and symmetry, and this should be apparent in their responses. This is a necessary condition for the extension of one-scale filters to multiscale. Starting from the least number of samples, two for even-length filters and three for odd-length filters, the z transforms of the filters in the set are computed. The lowest scale is then chosen where:

1. The frequency transform starts showing consistency at the higher scales. This is verified by checking the sampling relation between successive scales. The "pulse" waveform is an example of a one-scale filter where consistency is not found, and accordingly the samples are unsuitable for multiscale feature detection.

2. The filter's response must be bandlimited. This is necessary since the multiscale filters will be implemented using analysis filter bank. The cascade combination of filters becomes equivalent to reconstructing a larger support filter from a smaller filter. This is a special case of reconstruction, and aliasing must be avoided.

3. The highest resolution filter corresponds to the filter with the smallest support. Therefore, this filter's response satisfies the above two constraints, and has its bandwidth closest to .pi.. This results in a maximally decimated filter, and higher scales are oversamples.

The details of this procedure are now described for the Haar filter, and similar methods are applied to other filter operators without detailing the steps.

The sampling rules described above are applied to the Haar wavelet as follows: ##EQU1##

Samples of this wavelet give a type 1 filter that is symmetric with even length. Assume that the largest support filter has length 16. The Haar wavelet is uniformly sampled to get those 16 samples as shown in FIG. 2 at rate 4 (rate and scale are used interchangeably here). The smaller support filters are derived by further decimation as seen at lower rates in FIG. 2. To find the highest resolution filter G.sub.0 (z) which has the least number of samples, one needs to examine the frequency responses. We observe the following results of the different sampling rates: The first filter, of length two at rate 1, passes components at frequency .pi.. However, the second filter of length four at rate 2 does not. This is a first sign of inconsistency where the first filter is not bandlimited. In addition, the z-transforms of both filters do not seem to correspond to filters related by sampling. This is another critical sign of inconsistency. On the other hand, comparing the second filter to the third, and so on, consistency is observed in the responses. The conclusion is that the smallest filter, where consistency starts to show, must be of length four (at rate 2). Detection of features using a multiscale filter is reinforced from the consistency of response at all scales, and if the filter does not provide this property, then it is not suitable for multiscale feature detection according to the present invention. Other filters that may be employed in the present invention include:

Laplacian of Gaussian (LOG), where the smallest filter has length 17;

first derivative of Gaussian (FOG), where the smallest filter has a length of 9;

triangular waveform, where the smallest filter has a length of 5;

sawtooth waveform, where the smallest filter has a length of 9.

The filters in this set of acceptable filters are denoted herein by G.sub.k (z). These filters can be implemented in parallel as shown in FIG. 3, but a more efficient structure is to implement these filters in a filter bank as shown in FIG. 4. For some choices of G.sub.k (z) an exact filter bank implementation is possible, and in other cases an approximation is used, as will be discussed below.

The concepts discussed above may be extended to cover the 2D domain. It is assumed that a one-scale filter has been derived for optimal detection of the desired feature. The filter can be either separable or non-separable. For example, separable and non-separable Roberts and Sobel operators may be used to detect vertical and diagonal edges respectively. The size of such filter operators is chosen to cover the largest scale of interest.

This 2D filter is then used to generate a set of multiscale 2D filters by applying different sampling rates. The sampling lattice can be characterized by a 2.times.2 nonsingular integer matrix. In 1D, the relation between a signal x (n) and its decimated version x.sub.s (n) by a factor M is given by: x.sub.s (n)=x(Mn). In 2D, the decimation factor is given by the determinant of the decimation matrix D. The relation between the 2D signal and its decimated version is given by: x.sub.s (n)=x(Dn).

Two conditions are imposed on the decimation matrix D:

In the 1D approach described above, the decimation rate is 2. This choice is common in multiresolution analysis for dyadic scales, and can be easily implemented using two channel filter banks. By analogy to the 1D case, the 2D rate should achieve a sampling density as power of four. This results in the requirement that the determinant of D is 4.

However, the choice is not unique. This leads to the second condition which limits the alternatives.

In the development of the detector, it is assumed that an isolated feature indicator persists at all the scales. This implies that the decimator can not alter the shape of the filter.

Combining these two constraints leads to the unique decimation matrix given by: ##EQU2## The generated sampling lattice resulting from this decimation matrix is shown in FIG. 5. Note that this lattice does not necessarily lead to separable filters. For example, 2D filters with symmetry along diagonals can be designed to detect features at 45 degree angles. This is not possible with the separable processing approach of 1D filters which restricts the detection to horizontal and vertical directions only. Similarly, other non-separable filters can also be derived using the 2D method.

With this choice of sampling lattice, a one-scale 2D filter generates a set of 2D discrete filters, denoted by G.sub.i for each scale i. These filters are related by the decimation matrix D: G.sub.i (z)=G.sub.i+1 (z.sup.D). This decimation relation lends itself to a natural implementation as a filter bank. To derive this set, the one-scale 2D filter provides support large enough to cover the largest scale of observation and the rest of the filters are successively derived by decimation on the lattice of FIG. 5. As with the 1D filter, the multiscale 2D filter G.sub.i may be implemented in parallel, but a better solution is a tree structure filter bank as will be described below.

In designing the 1D filter bank, the relation between the parallel and serial multiscale filter implementation in FIG. 4 is discussed first. Given the desired set of filters G.sub.i (z), one needs to find the FIR filters H.sub.0 (z) and H.sub.1 (z) such that: ##EQU3## H.sub.1 (z) must match G.sub.0 (z), and since application restricts the choice of G.sub.0 (z), H.sub.1 (z) must follow. This gives the first part of the filter bank, and to provide symmetry in the response, it needs to be a linear phase FIR. For the other part H.sub.0 (z), several restrictions can be drawn by examining the structure. It should be noted that a cascade of filters, as shown in FIG. 4, is equivalent to one of the G.sub.k (z) filters. The properties of the equivalent filter at any scale should not be altered by the cascade leading to the next scale. For example, if H.sub.1 (z) is even-length and H.sub.0 (z.sup.2) is even-length, then the combination would not have an even-length and this would contradict the desired G.sub.1 (z). Therefore, H.sub.0 (z) needs to be symmetric of odd length to preserve characteristics, such as symmetry and anti-symmetry, across scales.

An exact solution for H.sub.0 (z) exists and is given by: ##EQU4## if the roots of ##EQU5## are also roots of G.sub.0 (z).

This is demonstrated as follows. Since G.sub.0 can be obtained by subsampling G.sub.1, then ##EQU6## In addition, G.sub.1 (z)=H.sub.0 (z.sup.2)H.sub.1 (z). Combining the two equations, yields equation (2).

The above result may not always be satisfied. In addition, H.sub.0 (z) is not guaranteed to be linear phase. The Haar filter is an example where this requirement is satisfied, and a linear phase FIR solution can be found. In this case, ##EQU7## and the exact solution is:

The Haar filter is one example of an exact solution, but others are possible. In general, it is difficult to extend one-scale filters to multiscale filters satisfying the filter bank constraints. However, the reverse process is straightforward. A one-scale filter can be easily drawn from a filter bank by a simple collapsing of the cascaded filters. This is achieved by convolution of frequency multiplication.

In cases, where the requirement in equation (2) is not satisfied, an approximate solution based on least-squares approximation can be made. The desired filter H.sub.0 (z) should be chosen to minimize the error function: ##EQU8## where h.sub.0i is the upsampled version of h.sub.0, r is the number of scales used in the analysis, and L.sub.i is the length of g.sub.i.

The following can also be proven. Let N+1 be the size of filter h.sub.0, E be an extension matrix of size (2N+1).times.(N+1), U.sub.i be an upsampling matrix of size (L.sub.i+1 -L.sub.i -1).times.(L.sub.i+1 -L.sub.i)/2, and G.sub.i the convolution matrix of size L.sub.i+1 .times.L.sub.i. They are ##EQU9## Furthermore, ##EQU10## V.sub.r (h.sub.0) is minimum if h.sub.0 =B.sub.0 C.sup.-1 a.sup.t.

The above procedure may be used to design a filter bank for 1D feature detection based on the set of multiscale filters that are derived from a one-scale filter. This filter bank can also be used for 2D applications by applying separable processing, however it is not suitable for 2D applications where non-separable filters are required. In addition, 2D filters have different support regions than 1D filters, and therefore the responses may be different.

Conceptually, the formulation of a 2D filter bank is similar to the one dimensional case discussed above. Decimation matrices replace integer decimation, and matrices replace vectors in the signal representation. The set of filters G.sub.i (z) is generated from a given one-scale 2D filter as explained in the previous section. The index i corresponds to the scale of observation, and it reflects the support size of the filter G.sub.i. One needs to design H.sub.0 (z) and H.sub.1 (z) such that: ##EQU11## Moreover, the phase and length of the desired filters need to satisfy certain constraints as follows. The length of H.sub.1 is the same as that of G.sub.0, and the length of H.sub.0 can be derived from the lengths of G.sub.0 and G.sub.1. Since symmetry and shape are preserved across the filters G.sub.i, the same is expected when upsampled versions of H.sub.0 are cascaded with H.sub.1. This accordingly constrains H.sub.0 to be symmetric.

An exact solution H.sub.0 (.omega.) exists and is given by ##EQU12## if the roots of the denominator are also roots of the numerator. J(D) represents the sampling density of the decimeter D, and N(D.sup.T) is the set of integer vectors that lie inside the fundamental parallelepiped of D.sup.T. This can be demonstrated as follows. G.sub.0 is the resulting filter when G.sub.1 is cascaded with the decimator D. In other words, it is the decimated filter version of G, and ##EQU13## On the other hand, from the equivalence of the filterbank structure to the parallel implementation G.sub.1 (.omega.)=H.sub.0 (D.sup.T .omega.)H.sub.1 (.omega.). Combining these two equations yields equation (5).

The above hypothesis may not always be satisfied and hence the exact solution of equation (5) may not work. In addition, H.sub.0 (.omega.) is not guaranteed to satisfy the desired symmetry constraints. The 2D Haar filter is an example where the hypothesis is satisfied, and an exact FIR filter solution can be found.

In cases where the hypothesis fails, an alternative solution is proposed based on least-square approximation. The desired filter H.sub.0 (.omega.) needs to minimize the error function: ##EQU14## where H.sub.0i is the upsampled version of H.sub.0, r is the number of scales used in the analysis, L.sub.i =M.sub.i N.sub.i is the size of the filter G.sub.i, and ** denotes 2D convolution.

It can be shown that V.sub.r (H.sub.0) is minimum if H.sub.0 of size M.times.N is given by: ##EQU15## with ##EQU16##

.sub.i is the 2D convolution matrix corresponding to the 2D filter G.sub.i, U.sub.i is the upsampling matrix corresponding to the rectangular lattice decimator, ##EQU17## E is the matrix used to impose the symmetry constraints on the filter, where h.sub.0 =Eh.sub.00 with h.sub.00 =C.sup.-1 a.

The dimensions of these matrices are as follows: v.sub.n is N.times.1, N.sub.n is MN.times.1 h.sub.00 is M.sub.00 .times.N.sub.00, E is MN.times.M.sub.00 N.sub.00, g.sub.i is M.sub.i N.sub.i .times.1, a.sub.i and a.sub.i a are M.sub.00 N.sub.00 .times.1, C.sub.i and C are M.sub.00 N.sub.00 .times.M.sub.00 N.sub.00, U.sub.i is (M.sub.i+1 +1-M.sub.i)(N.sub.i+1 +1-N.sub.i).times.(M.sub.i +1-M.sub.i-1)(N.sub.i +1-N.sub.i-1), and B.sub.i is (M.sub.i+1 +1-M.sub.i)(N.sub.i+1 +1-N.sub.i).times.(M.sub.00 N.sub.00).

In summary, a 2D filter bank is designed by the following steps.

Choose or design a 2D filter with a desired indicator response for the feature of interest.

Generate a set of 2D filters sampled on a rectangular lattice as described above.

Apply equations (5) and (6a) or (6b) described above to determine the equivalent filter bank.

These steps were followed to implement a 2D filter using a 2D vertical Haar filter bank, using a horizontal Haar filter bank, a 2D vertical Sobel filter bank, and a 2D diagonal Sobel filter bank. The results of this implementation are summarized in Table 1.

Once the filter bank is derived, the multiscale MAP correlator 20 that works with it can be designed not only to detect the feature in the presence of noise, but also to distinguish several characteristics of the feature such as edges and end-points. Whereas in the past decisions based on the outputs of the filter were often heuristic, the multiscale MAP correlator 20 according to the present invention is well structured and optimized to minimize the detection error.

To reduce the computation burden on the correlator 20 the chosen "indicators" in each filtered signal 14a-c are first located by the multiscale detector 16. Starting from the lowest scale, the detected indicator locations are then traced across the different scales. The rule for tracing is to connect an indicator location at the lower scale to the spatially closest indicator at the next higher scale. This is repeated for all the observable scales. The result is a "scale-space tree" giving us the locations of the indicators at each scale. This scale space tree is the input to the MAP correlator 20.

An optimal and structured MAP correlator 20 is derived with little or no information about the signal or the noise. The decision of the MAP correlator 20 is based on the actual filter outputs corresponding to the locations of the indicators in the scale-space tree. The correlator 20 is designed to minimize the probability of error. To simplify the design of the MAP correlator, it is assumed that the noise samples are independent, and identically distributed with a Gaussian distribution. The hypothesis testing problem becomes one of deciding whether a located feature corresponds to noise or an actual signal. Depending on the type of desired indicators, the correlator 20 may combine more than one response from each scale. For example, a zero-crossing is marked by checking a sign change in two consecutive locations. The two hypotheses are shown to have multivariate Gaussian distributions. The difference between the two densities lies in the mean. The test does not depend on signal to noise ratio, but only the desired signal. Details on the design of a 1D multiscale MAP correlator and a 2D multiscale MAP correlator will now be given.

Assume that the input signal to the multiscale filter bank 10 has the desired feature embedded in noise. Let us denote the input signal samples by x(k), where k represents spatial or temporal position. Denote the noise samples an n(k), and the feature samples by s(k),

The noise samples are assumed to be "independent and identically distributed" ("i.i.d.") with marginal distribution N(0,.sigma..sup.2). Therefore, the input samples x(k) are i.i.d. whose distribution P can be represented as: ##EQU18## It is an i.i.d. sequence with Gaussian distribution.

The outputs at different scales are linear combinations of the Gaussian input. Denote the output of a filter h.sub.i at scale i by the samples y.sub.i (k). Collect these samples from different scales into a vector y.sub.k =[y.sub.1 (k) y.sub.2 (k) . . . ].sup.T. Let H represent the convolution matrix associated with the filters at different scales h.sub.i, and x.sub.k the set of samples convolved with H at time (position) k,

H represents a linear transformation applied to the input.

The number of input samples needed has to match the largest filter. Let n be the size of the largest h.sub.i filter, and p the number of scales used in the filter. It can be shown that the distribution of the linear combination of the n-tuples x.sub.k is a Gaussian random variable. In addition, samples of the vector y.sub.k are jointly Gaussian. Their joint density function f is characterized by the first two moments. The mean of the observation vector y.sub.k can be derived by:

where m.sub.y.sbsb.x y and m.sub.x.sbsb.k, are the means of the input and the output respectively. The output correlation matrix R.sub.y.sbsb.k can be written in terms of the input correlation matrix R.sub.x.sbsb.k and H:

Since the noise samples at the input are i.i.d. with zero mean, the covariance matrix is the same as the correlation matrix. The mean and the correlation give the density at the output: ##EQU19## Denote f(y.sub.k) by f.sub.0 (y.sub.k) when only noise is present, and denote it by f.sub.1 (y.sub.k) in the presence of the desired signal.

If we denote by H.sub.0 the hypothesis when only noise is present, and by H.sub.1 in the presence of the desired feature, the two hypotheses are associated with the densities f.sub.0 and f.sub.1 respectively. To determine that a feature is present at time position k, the likelihood ratio test is applied to the two densities. L(y.sub.k)=f.sub.1 (y.sub.k)/f.sub.0 (y.sub.k). Assuming uniform costs in erroneous detections, the Bayes detector can be reformulated as: ##EQU20## Substituting by the density given above, the MAP correlator 20 becomes: ##EQU21## Note that the test performed by the MAP correlator 20 does not depend on the signal to noise ratio, but only on the signal/feature sought. This test is useful, when the desired feature is indicated by a single position, such as a signal maxima.

The test performed by the MAP correlator 20 is different when the feature indication is a zero crossing. In this case two consecutive positions are indicated for each feature, and these must be examined. The two consecutive positions are at k and k+1 to give y.sub.k and y.sub.k+1. The two vectors are combined by interleaving to give us Y.sub.k =[y.sub.1 (k)y.sub.1 (k+1)y.sub.2 (k)y.sub.2 (k)y.sub.2 (k+1) . . . ].sup.T. The indexing represents the different scales. Again here, Y.sub.k can be written as a linear combination of the input samples x.sub.k. Recall that y.sub.k =Hx.sub.k. Assume that the largest filter h.sub.i has size n. Take n+.sub.1 samples of the input centered at k, i.e. one extra sample added to x.sub.k, denote it by X.sub.k. Let ##EQU22## Then

The MAP correlator 20 for zero crossings can then be derived using the above equation (8), where H is replaced by H', y.sub.k by Y.sub.k, and x.sub.k by X.sub.k.

Detection of features in 2D signals, such as an image, is based on the indications collected across all scales of the multiscale filter 10 and detector 16. In this case, locating the indicators also depends on the direction of the feature. For example, maxima occur along the gradient direction of edges. Once these locations are detected at all the scales, the observations at those locations are linked to form the scale space tree. This tree is then fed to the MAP correlator 20 to decide whether the feature is present or not.

Assume that the input samples to the filter bank in FIG. 1 has the desired 2D feature (signal) embedded in noise. Denote those input samples by x(k,l), where k and l represent the spatial coordinates, the noise samples by n(k,l), and the feature samples corresponding to hypothesis j by s.sub.j (k,l), we have

The noise samples are assumed to be independent and identically distributed ("i.i.d.") with marginal distribution N(0,.sigma..sup.2). Consequently, the input samples x(k,l) are also i.i.d. and the distribution P.sub.j corresponding to hypothesis j can be represented as: ##EQU23## It is an i.i.d. sequence with Gaussian distribution, the difference between the hypotheses lies in the mean.

The filter outputs at different scales are linear combinations of the Gaussian input. Denote the output of the filter H.sub.i at scale I by the samples y.sup.(i) (k,l). Collecting these samples from different scales into a vector y(k,l)=[y.sup.(1) (k,l)y.sup.(2) (k,l) . . . ].sup.T. Let represent the convolution matrix associated with the filters H.sub.i from different scales, then ##EQU24## since represents a linear transformation applied to the input.

The input samples size matches the size of the largest filter. Let N be the size of the largest H.sub.i filter, and r be the number of scales used in collecting the observations. It can be shown that the distribution of the linear combination of the n-tuples x(k,l) is a Gaussian random variable, and the samples of the vector y(k,l) are jointly Gaussian. Their joint density function f is characterized by the first two moments. The mean of the observation vector y(k,l) is ##EQU25## where m.sub.y(k,l), and m.sub.x(k,l) are the means of the input and the output respectively. The output correlation matrix R.sub.yy can be written in terms of the input correlation matrix R.sub.xx. and : ##EQU26## Since the noise samples at the input are i.i.d. with zero mean, the covariance matrix is the same as the correlation matrix. The mean and the correlation give the density function of y(k,l). ##EQU27## Denote f(y(k,l)) by f.sub.0 (y(k,l)) when only noise is present, and denote it by f.sub.j (y(kl)) in the presence of the signal from hypothesis j.

Assuming uniform costs in erroneous detections, the Bayes detector can be reformulated as: ##EQU28## With the density given above, the 2D MAP correlator 20 is ##EQU29## The hypothesis testing problem becomes one of deciding whether a located feature corresponds to noise or actual signal. The two hypothesis have multivariate Gaussian distributions with different mean. The test does not depend on signal to noise ratio, but only on the desired signal. The 2D MAP correlator 20 decides that a feature is present at position (k,l) if ##EQU30## where: y(k,l).sup.T is an indicator location vector at k,l;

R.sub.yy.sup.-1 Hs.sub.1 (k,l) is a correlation vector at k,l; and ##EQU31## is a threshold value. Note that the test does not depend on the signal to noise ratio, but only on the feature signal sought.

The preferred embodiment of the invention is employed to detect edges in images produced by an ultrasonic imaging system. Referring particularly to FIG. 6 the image to be processed is received at input 100 as a 400 by 400 pixel array of image intensity data. This image data is applied to one input of a convolver 102 that forms the first stage of a multiscale filter bank indicated by dashed line 104. The multiscale filter 104 includes three additional convolvers 106, 108 and 110 which are connected in series to form the remaining three stages of the multiscale filter 104. Decimators 112, 114 and 116 connect the successive filter stages and sample every other pixel value to change the scale of successive filter stages. The outputs of the convolvers 102, 106, 108 and 110 form the respective multiscale outputs 120-123, where:

Scale 1 is at 120;

Scale 2 is at 121;

Scale 3 is at 122;

Scale 4 is at 123;

The convolvers 102, 106, 108 and 110 are identical and are constructed using digital signal processing ("DSP") components. The convolvers perform a convolution operation on the image array received at one of its inputs with a convolution kernel applied to its other input. The convolver 102 receives a convolution kernel H1 from an EPROM 126 and the remaining convolvers 106, 108 and 110 receive a convolution kernel H.sub.0. In this preferred embodiment a Haar filter function kernel is employed and H.sub.0 and H.sub.1 are derived as described above to yield the following: ##EQU32##

It can be appreciated that other filter functions can be implemented and their kernels stored in the EPROM 126 for application to the multiscale filter 104 when the need arises.

Each convolver 102, 106, 108, 110 can be implemented by Fourier transforms of the image array and the convolution kernel applied to its inputs, matrix multiplies them, and inverse Fourier transforms the result to produce the filtered output image array.

The decimators 112, 114 and 116 are also constructed using DSP components. They implement the following decimation matrix: ##EQU33## to sample every other pixel as illustrated by the sampling lattice in FIG. 5.

The filtered image arrays produced at the filter outputs 120-123 are applied to the inputs of a multiscale detector indicated generally at 130. Each of the filtered images is processed by a separate detector channel comprised of respective detectors 131-134 and masks 136-139. Each detector 131-134 scans the filtered image array applied to its input and indicates the location therein of the "indicator response" of interest. In the preferred embodiment a peak detector is employed, and the detectors 131-134 indicate the location of peaks, or maxima in the filtered image intensity data. These locations are output to the corresponding mask 136-139, which forms an "indicator response table" by reading the filtered image intensity value at each indicated "peak" location. Each channel of the multiscale detector 130 thus produces an indicator response table such as that shown at 140 in FIG. 7 which lists the location of each peak and the value of the peak. In the preferred embodiment, the functions performed by the multiscale detector 130 are carried out by a general purpose processor under the direction of a stored program.

The indicator response tables 140 are input to a feature signature generator 150. In the preferred embodiment the function performed by this generator 150 is carried out by a stored program which is executed to perform the steps indicated by the flow chart in FIG. 9. This program examines the indicator response table 140 for each scale, and produces a scale space tree ("SST") table which is illustrated at 142 in FIG. 8.

Referring to FIG. 9, the first step in the procedure is to initialize the relevant data structures as indicated by process block 160. This includes, for example, clearing the SST table 142. The indicator responses 140 for the first scale are then stored in the first row of the SST table 142 as indicated at process block 162. This is the highest resolution scale and the indicator response table 140 for scale 1 will contain the most entries. A scale pointer is then set to read the indicator response table 140 for scale 2, as indicated at process block 164, and a loop is entered in which the indicated peaks in scales 2, 3 and 4 are examined.

Referring still to FIG. 9 each indicator response table 140 is examined by setting a response pointer to "1" at process block 166 to read the location of the first indicated peak. If this location corresponds (plus or minus one pixel location) with a peak location indicated in the previous scale, as determined at decision block 168, the value for that location and scale is loaded into the SST table 142 as indicated at process block 170. Otherwise, the system loops back to examine the next indicator response after incrementing the response pointer at 172. After the last indicator response in the table 140 has been examined, as determined at decision block 174, the scale pointer is incremented at 176 and the indicator response table 140 for the next scale is examined. The process terminates when the last scale is finished as determined at decision block 178.

Referring particularly to FIG. 8, the resulting SST table 142 stores a set of n pixel locations at which peaks were found in the scale I indicator response table 140. Where corresponding peaks were found at higher scales, values are entered for those scales. It should be apparent that locations where peaks were found at all scales have a higher probability of being the location of an edge, whereas locations with indicator values at only scale I are likely noise or weaker edges. The values at each SST table location (i.e. each column in FIG. 8) is referred to herein as an "indicator location vector", and it is this vector which is examined to determine if an edge is present at that location. The values in each indicator location vector can be viewed as a "signature" of some feature in the image produced by the multiscale filter 104.

Referring again to FIG. 6 the SST table 142 produced by the feature signature generator 150 is examined by a MAP correlator 180 to determine which indicator location vectors therein represent true edges. The functions performed by the MAP correlator 180 are carried out by a general purpose computer which executes a program illustrated by the flow chart in FIG. 10.

Referring particularly to FIG. 10, the first step performed by the MAP correlator 180 is to initialize the data structures as indicated by process block 182. Such structures include, for example, a feature image array which is set to zero, or "black". As indicated by the above equation 11, a correlation vector is then calculated at block 184. This correlation vector may be viewed as the ideal indicator location vector produced by the multiscale filter 104 in response to an edge in the input image. Next, a threshold value is calculated at process block 186. This calculation is also indicated in the above equation 11. A loop is then entered in which each indicator location vector in the SST table 142 is read at process block 188, cross correlated with the correlation vector at process block 190, and the calculated correlation value compared with the threshold at decision block 192. If the indicator location vector is similar enough to the correlation vector that the calculated correlation value exceeds the threshold, an edge is present and the pixel at the vector location is set to a bright value at process block 194. Otherwise, process block 194 is skipped. This process continues until each indicator location vector in the SST table 142 is examined. When the last vector has been examined as determined at decision block 196, the feature image is complete and it is output at process block 198 to a display 200 (FIG. 6). The feature image indicates each detected edge as a bright line. This may be the final objective, or further steps may be carried out using this feature image to locate and identify shapes, and calculate their areas or volumes.