Method and Apparatus for Improved Digital Image Processing

The present invention relates generally to digital image processing and more particularly to digital filtering techniques involving convolution with two-dimensional filters having finite impulse response.

Digital image processing refers generally to methods of performing various tasks upon images stored in digital form. A familiar example of digital image processing is the restoration improvement of photographs or other images taken by satellites or astronauts in outer space and then transmitted to Earth. In simple terms, digital images may be thought of as an array of dots (commonly referred to as "pixels"), with each pixel being assigned a value in accordance with its brightness or luminance. For color images, one pixel may be assigned a value for each primary color of the system. The pixel array may be viewed mathematically as a matrix of data for purposes of digital image processing.

As noted above, digital image processing encompasses restoration, or enhancement either by "smoothing" or "sharpening," and numerous other techniques well known in the art. Each of these techniques is accomplished by some kind of linear convolutional processing or filtering, whether it be convolution, recursive filtering, or Fourier transform filtering. These techniques generally employ two-dimensional filters having finite impulse response (FIR).

Unfortunately, while applications of two-dimensional filtering are becoming of increasing importance, as, for example, in geophysical exploration, earthquake/nuclear test detection, sonar, radar and radio astronomy, in addition to picture encoding, filtering of large size data sets requires a great amount of computation using conventional filtering methods. An example of the computational problems encountered with conventional convolution or filtering techniques can be seen from Digital Image Processing, by William K. Pratt, Wiley-Interscience, 1978, Library of Congress catalog no. TA 1632.P7 1978, particularly chapters 1 and 9. In discrete convolution a filter is defined by an impulse response, or operator, which generally comprises a plurality of values, or parameters, arranged in matrix form. The operator matrix may be viewed as a mask which is then scanned over, or convolved with a generally larger pixel matrix to achieve the desired processing. Typical dimensions of the impulse response array are on the order of 15.times.15 or 21.times.21, but may be as large as 61.times.61. The pixel matrix typically has dimensions on the order of 512.times.512. Thus a single discrete convolution involves approximately 15.multidot.15.multidot.512.multidot.512=58,982,400 computations. It can thus be seen that a more manageable method of digitally processing an image has been required.

There have been some background investigations into more efficient convolution processing methods. In an article by J. H. MacClellan entitled "The Design of Two-Dimensional Digital Filters by Transformations", Proc. 7th Annual Princeton Conf. Information Sciences and Systems, pp. 247-251 (1973), it was shown that a one-dimensional zero-phase FIR filter could be mapped into a two-dimensional zero-phase FIR filter by substitution of variables. In a subsequent article by Mecklenbrauker and Mersereau entitled "MacClellan Transformations for Two-Dimensional Digital Filtering: I-Design, II-Implementation," IEEE Transactions on Circuits and Systems--Vol. CAS-23, No. 7, pp. 405-422 (July 1976), the MacClellan model was generalized and methods were developed to implement such filters by ordinary sequential filtering techniques. A third and fourth section of the article by Mecklenbrauker and Mersereau provided algorithms for determining parameters, and also examples.

The present invention provides a processing method which substantially reduces the number of computations required for discrete convolution while at the same time providing a good approximation of the results reached by conventional convolution. The method provides an approximation of the transfer function of a two-dimensional FIR filter implemented by direct or conventional convolution. In the method of the present invention, a new operator matrix having dimensions smaller than the conventional operator is defined, and may be termed a "small generating kernel." The small generating kernel (SGK) may typically have dimensions on the order of 3.times.3 elements.

The SGK is then scanned over the data matrix in the manner of conventional convolution, except that the SGK is re-scanned over the previous output data matrix a plurality of times, and different element values may be employed for the operator mask during each scan. The repeated scans are then weighted and/or summed in accordance with the desired processing effect, yielding a close approximation to the results obtained with direct convolution. In this manner the data are sequentially convolved with a plurality of masks defined by the small generating kernels until the dimensions of the sequential convolutions equal those which would result from a direct convolution. The selection of element values of the SGK operator may generally be determined by conventional error minimization techniques, such as the mean square error or Chebyshev error methods.

The present invention also includes apparatus employing the above-described techniques to provide a digitally processed display. In one embodiment, the input data is provided to a data storage device, and is sequentially scanned or convolved with a mask. The output of the mask is then supplied to a point operator circuit, which may for example be configured to perform a linear operation on the mask output. The output of the point operator circuit may then be supplied to a switching circuit.

By actuating the switching circuit, the information which forms the output of the point operator circuit is supplied to the data storage device and replaces the information originally contained therein. The new data is now ready for a second scan through the mask. However, the elements of the mask may now be varied in accordance with the desired filtering effects. The information which forms the output of the mask is then operated upon in the manner described above, and the iterations are continued until the sequential convolution attains the dimensions of the conventional convolution.

It is one object of the present invention to provide an improved method of digitally processing data.

It is another object of the present invention to provide a more efficient method of processing digital images.

It is another object of the present invention to provide an improved method of processing digital images through two-dimensional linear filtering.

It is a further object of the present invention to provide a technique for processing digital images involving synthesis of a filter having a finite impulse response which employs a small generating kernel technique to approximate conventional convolution.

Other and further objects of the present invention will be apparent from the following detailed description of the invention, taken in conjunction with the attached drawings which are described briefly below.

FIG. 1 shows in schematic block diagram form a system for processing monochromatic digital images according to the present invention.

FIG. 2 shows in schematic block diagram form a system for processing polychromatic digital images according to the present invention.

Referring now to FIG. 1, there is illustrated therein, in schematic block diagram form, a system for processing monochromatic digital images in accordance with the present invention. Data is provided to the system via an input terminal to a memory device 10, where it is stored for subsequent operation. The input data is typically an array representative of the brightness of the pixels which comprise a given image. Conventionally each pixel is scaled by brightness on a 256-Bit (2.sup.8) grey scale.

The memory 10 may be any of a number of conventional storage devices suitable for storing large data arrays. For conventional digital images, the dimensions of the array will typically be on the order of 512.times.512 pixels. Thus, for example, a commercially available solid state memory of suitable size may be used.

Once the input data has been entered into the memory 10, the data array is scanned by a mask circuit 20. The mask circuit 20 is a two-dimensional digital operator having a finite impulse response, and may mathematically be viewed as a matrix array, the numerical values of which are supplied by a sequencing device 30. The scanning process can then be thought of as a convolution of the mask array with the data array. Thus it follows that the mask circuit 20 comprises a network of multipliers and adders, with the multipliers producing the appropriate products of elements from the mask array and data array, which products are then added. For example, if the desire filter is represented by a 3.times.3 mask array, nine multipliers and nine adders are required for the mask circuit 20.

The sequencer 30 may be any of a number of conventional devices for supplying data to a particular unit at a particular time. For example, the sequencer 30 may comprise a series of registers and counters in combination with memory, or may comprise a microprocessor. The sequencer 30 thus stores the values of the mask array and provides them at the appropriate times to the mask circuit 20 as needed, in addition to other functions described hereinafter. The values supplied by the sequencer 30 are determined by any one of several conventional methods described hereinafter.

After the scan of the data array by the mask circuit 20 has been completed by convolving the mask array with the data array, the output is passed through a point operator circuit 40. The point operator circuit 40 typically modifies an input in accordance with a predetermined criteria which may be either linear or nonlinear. For example, the point operator circuit 40 may comprise a look-up table which provides a predetermined output for each input. This may be conventionally implemented by the user of an input addressable ROM (read only memory) or other suitable device. The point operator circuit 40 is likewise typically controlled by the sequencer 30.

The output of the point operator circuit 40 is then provided to a switching circuit 50, which is controlled by the sequencer 30. The output of the switching circuit 50 is connected to the input of the memory 10. By actuating the switch 50, the data which has been modified by the mask circuit 20 and point operator circuit 40 is entered into the memory 10 and replaces the originally entered data.

It is at this point that the distinctions between the present invention and the conventional convolution filtering techniques known in the art emerge. Conventionally, the output data array from the point operator circuit 40 is displayed as a final product; i.e., a completed convolution and a completely processed image. However, in the present invention the output of the point operator circuit represents only a first iteration of processing. For the present invention, the above process is repeated a plurality of times, the exact number of iterations being determined by the dimensions of the matrix which defines the conventional convolution operator and therefore also defines an ideal, or prototype, filter response.

As noted previously, a conventional convolution operator typically has dimensions on the order of 15.times.15 or 21.times.21 and may have dimensions on the order of 61.times.61 . Because of the substantially smaller dimensions of the mask array of the present invention, and also because the dimensions of convolved arrays build at a rate of (L-1) for each iteration after the first for an L.times.L array, the iterations required to approximate a 15.times.15 conventional operator (seven iterations) can be completed far more efficiently than a single convolution of the conventional convolution operator. This can be seen from a comparison of the number of computations required for a conventional convolution of a 15.times.15 operator and an M.times.N data array, with the number of computations required in the method of the present invention wherein a 3.times.3 mask array, which may also be referred to as a small generating kernel (or SGK) is used. For a conventional convolution, the number of computations required is 15.15.M.N, or 225.M.N. In contrast, the number of computations required with the present invention is (3.3.M.N)7 or 63.M.N. The magnitude of advantage available with the present invention becomes more apparent when it is seen that the number of computations required for a conventional convolution increases quadratically with increases in the size of the operator matrix, whereas the compuations required under the system of the present invention increases only linearly, as a function of the required number of iterations.

From the foregoing it can be appreciated that the present method provides a substantially more efficient method of digitally filtering images than conventional convolution, once the values of the elements in the mask array are selected. As will be discussed below, the element values of the mask array may be assigned through any of a number of conventional error minimization techniques. Further, the assigned values may be altered with each iteration or, depending upon which of the below-described implementations is used, the results of each iteration may be weighted according to a predetermined weighting input through the use of a weighting multiplier 60, and the results of the iterations may be accumulated in an accumulator memory 70. The processed image is then displayed by a conventional display means 80.

A detailed description of two embodiments of FIG. 1 is given below.

In the first embodiment, input F(j,k) is convolved with K.sub.1 (j,k), and its output is convolved with K.sub.2 (j,k), and so on until the output of the next to last convolution is convolved with K.sub.Q (j,k). At each stage of the convolution, the output is multiplied by a weighting factor .lambda..sub.q and added to the weighted outputs of previous stages. The final output, after Q stages of convolution, is the accumulated weighted output.

Applying the above process to the description of FIG. 1, it can be seen that the input F(j,k) is provided to the memory 10, and is then scanned or convolved with a first mask comprising a SGK operator K.sub.1 (j,k). The output is then provided back to the memory 10, and also to the weighting circuit (.lambda..sub.1). The output of the weighting circuit is then stored in the accumulator memory 70 and displayed (if desired) in the display 80. On the second iteration, the revised data is scanned by the SGK operator K.sub.2 (j,k), and again stored in memory, weighted (by .lambda..sub.2), accumulated, weighted, and displayed. The iterations continue in the manner described in connection with FIG. 1 until the dimensions of the sequential convolution equal those of the direct or conventional implementation.

From the mathematical description given above, it can be seen that the values of the SGK operators may be varied for each iteration, and likewise the .lambda. weighting inputs may be varied for each iteration. It will be apparent that modifications can readily be made to the optimization process which will permit either the weighting inputs or the SGK operators to remain constant.

Alternatively the weighting inputs .lambda. could be disposed of entirely as described below in connection with a second implementation. In the second implementation, input F(j,k) is convolved with K.sub.1 (j,k), and its output is convolved with K.sub.2 (j,k), and so on until the output of the next to last convolution is convolved with K.sub.Q (j,k). The final output after Q stages of convolution is the direct output from convolution array K.sub.Q (j,k). Applying the above to FIG. 1, each convolution step represents one iteration of the input, mask, point operator, and switch as described with respect to FIG. 1. However, the weighting multiplier 60, the weighting inputs, and the accumulator memory 70 are eliminated, with the output of the point operator circuit 40 being connected directly to the display 80. Thus the output of the point operator circuit may be displayed during each sequential convolution. Also, since substantial circuitry is eliminated, economic advantages attach to this method of implementation.

From the above, it can be seen that the method of the present invention depends, at least in part, upon the accuracy with which the set of mask arrays approximates the ideal filter impulse response, as represented by the conventional convolution operator or its transform function. One of the conventional techniques for minimizing the error between an approximation and an ideal is the technique of mean square error. Applying this technique to the choices of values for the SGK operator, let:

F(j,k)=input array (image) (J.times.K);

H(j,k)=ideal impulse response array for a prototype two-dimensional filter (2M+1).times.(2M+1);

K.sub.q (j,k)=SGK impulse response array of q-th stage to the method of the present invention (2L+1).times.(2L+1);

A.sub.q (j,k)=output of the q-th convolution stage of the present invention where A.sub.0 (j,k)=F(j,k)

G(j,k)=ideal output response; and

G(j,k)=approximate output response according to the method of the present invention.

Then the conventional output response is

where the convolution operation is defined as ##EQU1## The output response for the first embodiment of this invention is ##EQU2## or equivalently ##EQU3## The output response for the second embodiment of this invention is

or equivalently

By using the technique of mean square error, the error is: ##EQU4## However, because F(j,k) appears in the space domain expansions for both G(j,k) and G(j,k) above, the error for the first embodiment is ##EQU5## and for the second embodiment the error is ##EQU6## where C is a proportionality constant.

The optimization method will now be described for the first embodiment of the invention in which each of the small generating kernels K.sub.q (j,k).ident.K(j,k) are identical. The optimization method generalizes directly for the case of different small generating kernels and also for the second embodiment.

For purposes of derivation, the desired impulse response H(j,k) can be expressed in the Z-transform domain as ##EQU7## The Z-transform of the equivalent impulse response array is ##EQU8## where ##EQU9## The mean square error may then be expressed in the Z-transform domain as ##EQU10## where .GAMMA..sub.1 and .GAMMA..sub.2 are unit circles of the z.sub.1 and z.sub.2 complex planes. The problem may then be split into two parts, first a nonlinear problem in K(i,j) where -L.ltoreq.i.ltoreq.+L,-L.ltoreq.j.ltoreq.+L; and second, a linear problem in .lambda..sub.q where q=1, . . . , Q.

Choosing N.sup.2 samples of H(z.sub.1,z.sub.2), K.sup.q (z.sub.1,z.sub.2), .beta.(z.sub.1,z.sub.2) on .GAMMA..sub.1 .times..GAMMA..sub.2 yields the N.times.N arrays H(m,n), K.sub.q (m,n), .beta.(m,n), respectively where the tilde overbar indicates the discretization process. The optimization problem then reduces to finding .lambda..sub.q and K(i,j) such that ##EQU11## is a minimum. By stacking the columns of H(m,n) and K.sup.q (m,n) into the vectors h and k.sup.q, respectively, and defining P as a N.sup.2 .times.Q matrix ##EQU12## the criterion becomes

where the underbar denotes vector and matrix symbols, the asterisk denotes the complex conjugate relation, and .lambda. is the vector constructed with {.lambda..sub.q, q=1, . . . , Q}.

One can easily see that for each P matrix, the optimum value of .lambda. is given by

where the superscript (-1) denotes a matrix inverse and the superscript (T) denotes a matrix transpose. The derivatives of the criterion with respect to K(m,n) are given by ##EQU13## where R.sub.e {.multidot.} denote the real part of the function and ##EQU14## The derivative d.epsilon./dK(j,k) is then obtained by the inverse discrete Fourier transform of E(m,n). Thus we have derived an analytic expression for the gradient which can be easily computed and makes the use of the descent gradient type algorithm extremely appealing.

Although the computations required for calculating the gradient necessary for optimization of the element values of the mask array may be calculated by hand, the more expedient alternative is to perform these calculations by computer. Thus, a computer program suitable for performing such calculations is attached hereto as Appendix A. Alternately, the results of the approach shown above may be achieved using the Chebyshev method. It is believed that those skilled in the art will have no difficulty in making the modifications necessary to implement the Chebyshev method, given the teachings herein.

As an example of the optimization, the following are the optimum values of a 3.times.3 small generating kernel and the 12 weighting values for approximation of a 27.times.27 conventional convolution operator. The mean square error is approximately 4.7 percent.

Referring now to FIG. 2, there is illustrated therein in schematic block diagram form a system for digitally processing multiple band images, such as color pictures. As will be apparent from a comparison with FIG. 1, FIG. 2 generally comprises three separate channels each of which is similar to the monochromatic single channel illustrated in FIG. 1. For conventional color images, the channels will represent blue, red, and green bands. However, other multiple band systems may employ different constituent bands. Further, more than three bands may be employed. It is believed apparent to those skilled in the art that the system of the present invention may be expanded to include as many bands as are required for a particular technique.

In the system of FIG. 2, input data is supplied to each of the three memory devices 100a, 110b, 110c. The input data is generally different for each band. The input data arrays are then sequentially scanned by their respective mask circuits 120a, 120b, 120c. The elements of the mask array are determined in the manner described above for FIG. 1. It should be noted that the SGK operator may vary for each band, as well as varying for each sequential convolution, since the desired filtering may require emphasis or deemphasis of a particular band. As with FIG. 1, the elements of the mask arrays are supplied to the respective mask circuits 120a-c by a sequencer 150 through conrols A, C and E, respectively.

The outputs of the mask circuits 120a, 120b, 120c are then passed through respective point operator circuits 130a, 130b, 130c as described in connection with FIG. 1. Again, the point operator circuits 130a-c are controlled by the sequencer 150 through buses B, D and F, respectively. The outputs of the point operator circuits 130a-c are then supplied to a display means 160, and are also fed back to their respective memories 110a-c through a switching circuit 140. As with the system of FIG. 1, the switching circuit 140 is controlled by the sequencer 150. As with the system of FIG. 1, the input data is sequentially scanned by the mask circuits until the dimensions of the output have been built from those of the SGK operators to those of the direct implementation operator, as discussed in detail with FIG. 1.

Having fully described one embodiment of the invention, it is to be understood that numerous variations and equivalents which do not depart from the present invention will be apparent to those skilled in the art based on the teachings herein. Therefore, it is not intended that the scope of the invention be limited to the details disclosed herein, but rather that the invention shall include the variation and equivalents which flow therefrom. While square small generating kernels have been disclosed herein, rectangular or other dimensioned small generating kernels may equally well be used. ##SPC1## ##SPC2## ##SPC3## ##SPC4## ##SPC5## ##SPC6## ##SPC7## ##SPC8## ##SPC9## ##SPC10## ##SPC11## ##SPC12## ##SPC13## ##SPC14## ##SPC15## ##SPC16## ##SPC17## ##SPC18## ##SPC19## ##SPC20## ##SPC21## ##SPC22## ##SPC23## ##SPC24## ##SPC25## ##SPC26## ##SPC27## ##SPC28## ##SPC29## ##SPC30## ##SPC31## ##SPC32## ##SPC33## ##SPC34##