The practical significance of two-dimensional deconvolution in echography

ULTRASONIC

THE

IMAGING

9,

PRACTICAL

106-116

(1987)

SIGNIFICANCE

T.J.M.

Jeurens,

OF J.C.

6200

TWO-DIMENSIONAL

Somer:

DECONVOLUTION

F.A.M.

Department University P.O. MD Maastricht,

of of

Smeets,

and

A.P.G.

IN

ECHOGRAPHY

Hoeks

Biophysics Limburg

BOX

616

The

Netherlands

This paper evaluates deconvolution (inverse filtering) as applied to ultrasonic imaging systems, and discusses the obstacles which are encountered employing the technique in practice. A minicomputer is used to generate artifical echo signals, simulating rf signals resulting from a set of point reflectors in a homogeneous medium, as recorded by an electronically focused group-steered linear array scanner. Two-dimensional deconvolution in combination with a Wiener noise reduction filter (i.e., a Wiener-Inverse filter) is applied to these simulated rf signals, which were The efficacy of the Wiener-Inverse filter is contaminated with white noise. defined in te,rms of its ability to resolve two point reflectors with a lateral spacing equal to the local -6 dB width of the ultrasonic beam. In favorable circumstances, the targets are resolved at signal-to-noise ratios (SNR) better than 20 dB, where SNR is defined as the maximum signal power divided by the average noise power level. Nonlinear effects due to quantization or signal clipping are investigated. In order to improve the resolution of an rf signal with a dynamic range of 40 dB, the input signal should be digitized at a minimum of 12 bits. The problem of signal clipping can be circumvented by oversampling. The two-dimensional Wiener-Inverse filter is defined in terms of both temporal and spatial properties of the insonification. Effects of wave diffraction give rise to a depth-dependent ultrasonic beam. As a result of a misfit of the Wiener-Inverse filter and the local properties of the ultrasonic beam, erroneous noisy texture arises in the image. Adaptation of the Wiener-Inverse filter with respect to the beam at the expense of a rather large comproperties gives acceptable results, o 1987 Academic Press. Inc. putational effort. Key

words:

Computer spatial

I.

INTRODUCTION

simulation; filtering.

deconvolution;

linear

array;

resolution;

In the past decade, many papers have been published about deconvolution (or inverse filtering) in the field of echography. In spite of some promising results on computer simulated data and on echographic signals the technique has not, to our knowledge, been recorded from phantoms, applied in commercial echo systems. Theoretically, the resolving power of an ultrasonic imaging system is determined by the frequency content of the emitted signal and the apertureangle of the transducer. Selection of these system parameters depends on and is always a tradeoff between the desired depth range the application, and the resolving power. In the field of echocardiography, for instance, small transducers operating over a large depth range are required. Conse' Author

to

whom

correspondence

0161-7346/87 $3.00 Copyright 0 1987 by Academic Press, Inc. All rights of reproduction in any form reserved.

should

be

106

addressed.

quently, transducer lable spatial

the unfavorable size demand information

TWO-DIMENSIONAL

DECONVOLUTION

conditions careful processing, optimally.

concerning in order

emission to bring

frequency out the

and avai-

As far as deconvolution is concerned, mostly one-dimensional signal processing is investigated in the literature, either in the axial [1,21 or lateral [3-i'] direction. However, an ultrasonic imaging system maps the three-dimensional object information into a two-dimensional data structure. The received data are not independent in the axial and lateral directions, so signal processing should be carried out two-dimensionally. The concept of deconvolution is based on the assumption that the information about the interrogated medium and the characteristics of the measuring configuration are mixed in the received echo signals in a convolved way [El. The frequency characteristic of an echographic system in both axial and lateral directions is comparable with that of a bandpass It limits the broad band frequency content of the observed medium, filter. resulting in images with poor resolution. Moreover, different tissues produce rather similar textures in the image, and can hardly be distinguished by eye. Two-dimensional inverse filtering of the received echo signals by the frequency characteristic of the measuring system attempts to remove its band limiting properties, in order to restore the image of the medium. The deconvolution process amplifies the frequency components with less power relative to frequency components at higher power levels within the passband of the system. However, the low power frequency components have a rather bad signal-to-noise ratio so the deconvolved signals exhibit an unacceptable noise level. Therefore, a noise filter should be applied, which, unfortunately, reduces the efficacy of the deconvolution filter.

a posterior-i two-dimensional This paper evaluates the features of deconvolution as applied to ultrasonic imaging equipment. A PDP 11/73 minicomputer is used to generate artificial echo signals, simulating rf signals resulting from a set of point reflectors in a homogeneous medium, as recorded by an electronically-focused group-steered linear array scanner . Section II describes the two-dimensional deconvolution process. In section III, the basic principles of the computer programs generating the artificial echo signals are explained. Section IV discusses the results of the two-dimensional deconvolution procedure considering several practical problems, such as noise, guantization effects, and spatial invariance of the ultrasonic beam.

II.

TWO-DIMENSIONAL The

infinitely

DECONVOLUTION

characteristic short,

of 3-dimensional

&(x,y,zf

= lim N3 rect(Nx).rect(Ny).rect(Nz) NW

The &function applied as the system impulse racteristic.

a

OF ECHO-SIGNALS system is symmetrical

contains all frequencies input signal to a system, response, which defines

specified pulse,

by its the Dirac

response &-function,

to

an

(1)

with the the

equal output system

amplitude signal frequency

and phase. represents transfer

If the cha-

For ultrasonic imaging systems, an echo signal, returned by a highly reflective point target in a nonabsorbing medium provides, after proper amplification, the system impulse response as a two-dimensional signal. In our model, we will assume that no reflections from outside the scanning plane are involved, so the medium can be described in two dimensions. The

107

JEURENS

lateral t

(z)

1

a. Single point target, signal s,(x,z), returned Multiple point targets. returned by the point

system returned

impulse response, by a point target s

where sional given cated signal

AL

(xl

axial

Fig.

ET

1

(x,2)

*

denotes plot of by Zz=ct, configuration depicted

located at (x , z ). b. Two-dimensional by a single Opoi& target at (x0, zo). d. Two-dimensional rf signal sE(x,z), target configuration of figure lc.

h(x,z), located

= h(x,z)

is

defined (x0,

at

* 6(x-x

o1

by zo)

z-zo)

the (Fig.

echo signal, la) according

rf c.

s,(x,z), to

(2)

I

two-dimensional convolution. Figure lb shows a two-dimens,(x,z). The relation between axial distance z and time t is where c is the sound propagation velocity. A more compliof point targets, given in figure lc results in the rf in figure Id. The latter signal can be described by N

SN(X,Z)

=

1

h(x,z)

*

Ai.5(x-x.,

1

i=l = h(x,z)

*

i? Ai. i=l

6(x-x.,

1

Z-zi)

z-zi)

=

.

(3)

Because of the blurring property of h(x,z), the single point targets often can no longer be distinguished. It is clear from Eq. (3) that a two-dimensional deconvolution of sB(x,z) by the function h(x,z) provides the desired information about the location and amplitude of the targets. The function h(x,z) can be determined from Eq. (21, using the a priori knowledge that s,(x,z) represents the echo signal, reflected by only one point target at a known position in a homogeneous, lossless medium. Deconvolution convolution in the complex spectra in formed by dividing is defined by

S(kx,kz)

is executed in the ,two-dimensional Fourier-domain. As spatial domain is equivalent to multiplication of the the Fourier domain, spatial deconvolution can be pertwo spectra. The two-dimensional Fourier transformation

=

+,I!

s(x,z).exp(-j2a(kxx+kZz)dxdz

108

,

(4)

TWO-DIMENSIONAL

DECONVOLUTION

where k and k represent the spatial frequencies corresponding to x and in the gpatialxdomain, respectively. Let the Fourier transforms of s,(x,z), h(x,z) and sR(x,z) be denoted by S,(kx,ke), H(kx,ks) and SR(kx,ks), respectively. Eqs. (2) and (3) can be written U-I the Fourier domain as S,(kx,ks)

= H(kx,kr).exp(-j2a(kxxo+kzso))

SN(kx,kz)

= H(kx,kz).

Knowing (xo,zo) as a deconvolution H(kx,kz), we get

in

Eq. filter,

SD(kx,kZ)

(5)

N 1 i=l H(k beinG'&e

(S),

z

Ai

k

exp(-j2a(kxxi+kzzi))

)

can be inverse

(61

.

solved. of the

If we system

choose R(k character%&

,k

)

= R(kx,kx).SN(kx,kZ) 1 = H(kx,kr) =LA F!

Inverse

Fourier

'SN(kx'kx) .

: exp(-j2Wkxxi+kxzi)

transformation

of

SD(kx,kZ)

provides

the

(7)

deconvolved

signal

N y,(x,z)

=

1 i=l

Ai.6(x-xi,

(8)

.

Z-Zi)

,kr) become very large where H"( kx,kx) tends to zero. Consenoise outside the system passband

However, the values of the inverse filter the system frequency transfer characteristic quently, the frequency components of the will be enhanced unacceptably.

R(k

The optimal noise filter, which combines a minimal signal distortion and a maximal reduction of the noise, is known as the Wiener filter [9). Its spectral characteristic is given by the ratio of the signal power density spectrum and the sum of the power density spectra of signal and s qyared noise. The signal power density spectrum is estimated by the amplitude spectrum of the system transfer characteristic, 1 H(kx,kx)I , so the noise filter becomes

tWkx,kZ)I2 Hw(kx,kz)

where IN(kx,k combination the restoration

=

IH(kx,kz)

)I 2 represents o H the inverse filter

R(kx,ks)

=

I*

+ IWkx,kz)

the filter

power and

density spectrum the Wiener noise

- -/H(kx,kz)12

H* (kx,kz) =

IH(kx,kz)

(9)

,

of the filter

noise. results

The in

jH(kx,kz) I2

1 H(kx,kz)

I2

1'

+

IN(kx,ka)l'

109

+

IN(kx,kz)12

(10)

JEURENS

which is called the conjugate of H(kx,ks). filter if the spectral to the signal power (k ,k ). On the other signa I to noise ratio

ET AL

Wiener-Inverse filter. H*(kx,ks) denotes the complex The restoration filter ten$s to a pure inverse noise power density IN(kx,kZ)[ is small as compared density for a given spatial frequency component hand, R(k x, ks) turns into a matched filter when the is Poor:

H* (kx,kz) R(kx,ks)

In this case, frequency band, the expense of III.

z

filter suppresses unwanted and induces an improvement of degrading the resolution. OF

s(t)

P r t PO C

pou Equation per

to

unit

the

TWO-DIMENSIONAL

=

=

s(t-r/c) r

strength

phenomena outside the signal-to-noise

the

at

of The by

,

-$&30u(t)~

states that the at the position of the monopole

system

ratio

simulation program uses the description in a homogeneous, lossless medium. in a homogeneous medium is determined

(13)

= pressure variation = distance from the point =ti.me = density of the medium = sound propagation velocity = mass flux of the monopole. (13) time

(11)

RF SIGNALS

algorithm of the field of a monopole of a point source

p(r,t)

where

.

I2

the

SIMULATION

The basic the pressure pressure field

IN(kx,kz)

function of the

source

in

s(t) monopole,

the

medium

represents and

in

the this

total way

is

mass flux related

[lo].

Using Eq. (12), the pressure field of any transmitting surface can be calculated by integrating the contributions of the single monopoles, composing the surface. In the same way, Eq. (12) can be used to describe the pressure field of a wave reflected by a hard point scatterer. The point reflector is treated as a point source, the strength of which is given by the local strength of the incident wave field. Accordingly, the pressure field of an arbitrary point on the transducer surface can be calculated if the distance between that point and the point reflector is known. The superposition of the received pressure fields of all the points covering the transducer surface provides the received RF signal. The advantage of constructing an echo response in the elementary way as described above, is the extreme flexibility in choosing transducer size as well as the shape and dimensions of the reflecting oband curvature, A drawback is that the program is rather inefficient and consumes a ject. lot of computer time when large objects are involved. However, it was not the intention of this study to optimize the computer simulation of ultrasonic scanners with respect to efficiency and computational speed, but to create a class of two-dimensional rf signals, which is representative for To save computation time, signals assessed by an ultrasonic imaging system. it was decided to simulate purely two-dimensional targets: no information from outside the scanning plane is involved. Moreover, the distance between

110

TWO-DIMENSIONAL

Fig.

2

individual number surface,

B-scan image of an artificial 2.5 MHz group-steered linear reflectors, axially spaced by a x 1 = 170 x 85 mm.

of

monopoles monopoles thus reducing

DECONVOLUTION

signal, array 5 mm.

is set to a quarter to be considered for computer processing

simulating scanner The total

of a wavelength. the description time.

on size

the a

response of a set of point of the image is

This limits the of a transducer

The features of Wiener-Inverse filtering will be evaluated using artificial two-dimensional rf signals simulating an electronically-focussed (number of elements per group is group-steered linear array scanner seven, Figure 2 presents the B-scan image of a element spacing is 1.9 mm). reflected by a set of point targets in a homogeneous simulated signal, medium. The envelopes of the signals are detected using a Hilbert transformation. The Hilbert pairs are squared, added and displayed logarithmically on a 30 dB scale. IV.

RESULTS

AND

DISCUSSION

As stated before, the concept of deconvolution is based on the assumption that the information about the interrogated medium and the characteristics of the recording system are convolved in the received echo signals. Convolution is a linear process, complying with the principle of superposition. deconvolution can only be applied to signals which are Therefore, processed in a linear way like the uncompressed rf signals rather than the video-detected envelopes. Here, the first practical problem is encountered: the uncompressed rf signals which exhibit a large dynamic range must be digitized, which is a nonlinear process causing quantization errors comparable with white noise superimposed on the signals. Our first experiment concerned the efficacy of the Wiener-Inverse filter in the presence of white noise. A returned tance sional

two-dimensional rf signal was computed simulating an echo signal by a point target in the far field of the transducer, at a disof 128 mm (Fig. 3a). This signal was defined to be the two-dimensystem impulse response. The local -6 dB width of the ultrasonic

111

JEURENS

Fig.

3

ET AL

The axial range of these B-scan images is expanded by a factor of 4 relative to the lateral dimensions. The display range is 30 dB. impulse response: a. System echo returned by a point target at a distance of 128 mm from the transducer. b. Image of two point targets at the same depth of 128 mm, laterally spaced by 8 mm. Gaussian noise with an average power level of 40 dB below the signal maximum is added to the rf signal. c. rf signal of figure 3b deconvolved by the rf signal of figure 3a. SNR better than 30 dB.

beam could be measured from the lateral pulse width, being 8 mm, equivalent to 13 wavelengths (emission frequency 2.5 MHz). Next, an rf signal of two point targets was created, located at the same depth, laterally spaced by 8 mm. Figure 3b shows the B-scan image of this signal after adding white Gaussian noise (zero mean) with an average power level of 40 dB below the signal maximum. In the B-scan image, the echo can not be recognized as originating from two point targets, because the signal amplitude of the central echo line between the targets (Al) exceeds the amplitude of the echo lines at the exact positions of the targets (A21 by 1 The ratio Al/A2 was taken as a measure of filter efficacy. Table I dB. shows the filter efficacy as a function of the spatial signal-to-noise ratio, SNR, which is defined as the ratio between the maximum signal power and the average noise power2 level. As white noise was applied, its power in Eq. (10) is represented by a constant. The density spectrum IN(kx,kz!l filter output signals exhibit a SNR better than 30 dB in all cases. As expected, a better SNR in the input signals allows for a better filter performance. We may conclude that under favorable circumstances (a homogeneous medium and an inverse filter that "fits" exactly), the two points can be resolved at signal-to-noise ratios better than 20 dB. Figure 3c shows the deconvolved result for SNR=IO dB. The of additive adjusted efficacy, The filter

using quantized signals instead previous experiment was repeated, The noise level of the Wiener-Inverse filter was white noise. according to half a discretization step. Table II shows the filter as a function of the quantization intensity. as defined before, results indicate, that the non-linear effect due to quantization

112

TWO-DIMENSIONAL

Table

I.

Filter

DECONVOLUTION

efficacy

20 unfiltered

vs.

SNR

0 +1

may be treated in a similar way as additive white noise. Suppose the rf Table II indicates that the filter input signals are digitized at 12 bits. efficacy of full range signals amounts to -17 dB. However, rf signals at a level of 40 dB below the signal maximum are quantized at 4 bits. According to table II, weak signals in the least significant 4 bits cannot be improved. Thus, a restoration of the rf input signals over a dynamic range of 40 Weak signals can be procesdB requires a digitization of at least 12 bits. sed satisfactorily if they are amplified before guantization. Then, however, strong signals are pushed into saturation after digitization, due to the limited dynamic range of the A/D-converter. This is another nonlinear artifact. Figure 4 shows an example of the deconvolution of a two-point target echo signal, which was clipped at a level of 80 percent of the signal maximum before signal processing: the nonlinear filtering errors which are clipped by the appear as noisy texture in the image. Signals A/D-converter can be restored if the signals are oversampled. The clipping of the rf signals induces higher harmonics , which can be removed by digital filtering if a sampling frequency is used high enough to describe these The oversampling should prevent higher order frequency bands unambiguously. a fold-over of the higher order distortion bands into the signal band, which makes a reconstruction of the clipped signals impossible. Thus, the application of a sampling frequency high enough to include the distortion frequencies due to signal clipping might alleviate the saturation problem. The removal of the higher order bands by digital filtering should be done before the Wiener-Inverse filter is applied. is the fact that the deconvolution The second fundamental problem process does not account for the diffraction effects. In Eq. (21, the system impulse response is a function of lateral position, x, and depth, z. In the case of a linear array transducer, the system characteristic is independent of a lateral shift and is valid over the entire lateral scanning width. However, the effects of wave diffraction give rise to a variation of the ultrasonic beam, which is a function of depth (Fig. 2). Consequently, the two-dimensional system characteristic is depth dependent. To evaluate this depth dependancy, two arrays of point reflectors were simu-

Table

II.

Filter

efficacy

vs. quantization

113

intensity

JEURENS ET AL

Fig.

4

Deconvolution of a cl .ipped two-point target echo signal. Display range 30 dB.

lated, containing targets positioned at various distances from the transducer. The lateral spacing of the arrays was 8 mm, equivalent to 13 wavelengths at 2.5 MHz (Fig. 5a). The detected envelopes of the generated rf signals are shown in figure 5b. If we define the system characteristic in terms of a response from a point target in the center of the image, WienerInverse filtering of this image provides the result depicted in figure SC, assuming a SNR of 40 dB. Only one small area in the center of the image, where the system characteristic was defined, shows an acceptable filter result. This means that it is necessary to divide the ultrasonic field into small segments, each having its own Wiener-Inverse filter. The size of a segment is related to the local inhomogeneity of the ultrasonic beam. As an example, we have defined five equally-sized overlapping segments with five Wiener-Inverse filters. Figure 5d shows the result of the depth-dependent deconvolution. Note that in practice the system impulse response as well as the local SNR depend on the distance from the transducer. Both the spatial variance of the system impulse response inside the segments as well as the local SNR ratio affect the attainable resolution of the image. V.

CONCLUSIONS

In practice, the two-dimensional deconvolution process is a complicated and time consuming technique, and is therefore, at least for the moment, not suitable for real-time echo systems. Under ideal circumstances, using a rather homogeneous medium and an inverse filter that fits exactly, the technique provides an improvement of the resolution for an SNR better than 20 dB. In order to avoid erroneous texture due to nonlinear effects during data acquisition, we recommend a sampling rate at at least 6 times the One of the central frequency, and a quantization of more than 12 bits. severest problems employing the technique is the spatial variance of the array and mechanical sector scanners, the ultrasonic beam. For linear Wiener-Inverse filter has to be adapted in depth, while phased array scanners require an angular dependence in addition. A misfit of the restoration filter to the local beam properties causes noisy texture in the image, which might hamper proper tissue differentiation. Improvements of the image with regard to resolution and tissue differentiation are rather disappointing considering the computational effort which is required for this technique.

114

TWO-DIMENSIONAL

DECONVOLUTION

lateral axial

t

Fig.

5

a. Double array of point targets. lateral spacing is 8 mm. b. B-scan signal, simulated on the target c.

Deconvolution

by

an

impulse

Axial spacing is image corresponding configuration of

response,

which

is

5 mm, the to the rf figure 5a.

defined

in

the

center of the image. A SNR of 40 dB in the input signals was assumed. Display range is 30 dB. d. Depth-dependent deconvolution, using five equally-sized overlapping segments, where a local system impulse response is defined. A SNR of 40 dB in the input signals was assumed. Display range is 30 dB.

ACKNOWLEDGEMENT The authors manuscript.

are

most

grateful

to Mrs.

R.M.

Hanssen

for

preparing

the

REFERENCES [II

Martens, "Optimal"

W.L.J., Inverse

Somer, Filtering

J.C., Hoeks, and Matched

115

A.P.G., Filtering

and Smeets, F.A.M., and their Relations

JEURENS ET AL

to In Vivo Echographic Tissue Differentiation, in Ultrasonic Tissue Characterization, J.M. Thijssen, ed., pp. 245-256 (Stafleu's Scientific Publishing Company, Alphen aan den Rijn/ Brussels, 1980). [2]

Liu, C.N., vement of 66-75

[3]

Fatemi, ultrasonic

M.,

and Waag, R.C., abdominal images,

Digital processing IEEE Trans. Med.

for improImag. MI-2,

(1983).

Hundt, E.E. and Trautenberg, E.A., data by deconvolution, IEEE Trans.

processing Ultrason.

Digital

Sonics

of ultrasonic SU-27, 249-252

(1980).

[41 Vollmann, W., deconvolution, [5]

enhancement of Sonics Ultrason.

ultrasonic SU-29,

B-scan images 78-83 11982).

t71 Vaknine, B-scans

and Wing, M., Imaging 2, 1-12

R. and Lorenz, W.J., before image generation,

Fatemi, M., and Kak, A.C., formation and a technique

Lateral

deconvolution

of

ultrasonic

(1984).

Lateral filtering of medical ultrasonic Ultrasonic Imaging 2, 152-158 (1984).

Ultrasonic B-scan imaging: for restoration, Ultrasonic

theory of image Imaging 2, 1-47

(1980).

tg]

[lOI

by

inverse filtering Schomberg, H., Vollmann, W., and Mahnke, G., Lateral of ultrasonic B-scan images, Ultrasonic Imaging 5-, 38-54 (1983).

[61 Robinson, D.E. beams, Ultrasonic

[81

Resolution IEEE Trans.

Wiener, N., Extrapolation, Interpolation and Smoothing Time Series (MIT Press, Cambridge, MA: Wiley, New York, Berkhout, Company,

A-J., Seismic Migration. Amsterdam/Oxford/New York,

116

(Elsevier 1982).

Scientific

of Stationary 1949).

Puhlishing