IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-24, NO. 4, JULY 1977

[22] [23]

[241 [25]

tus measuring oxygen consumption and expired 02 and CO2 in small subjects," J. Appl. Phys., vol. 29, pp. 731-739, November 1970. F. G. Benedict, "A portable respiration apparatus for clinical use," Boston Medical & Surgical J., vol. 128, pp. 667-678, May 1918. M. Pettenkofer, "Uber die Respiration," Annal. ChemL Supple., vol. 2, pp. 1-5 2, 1862. R. H. Murray, A. Marko, A. T. Kissen, and D. W. McGuire, "A New Miniaturized Multichannel, Personal Radiotelemetry System," J. Appl. Physiol., vol. 24, pp. 588-592, April 1968. A. C. Guyton and C. A. Farish, "A rapidly responding continuous

Digital

337

consumption recorder," J. Appl. Phys., vol. 14, pp. 143-145, January 1959.

[26] C. T. Kappagoda, J. B. Stoker, and R. J. Linden, "A method for the continuous measurement of oxygen consumption," J. Appl.

Phys., vol. 37, pp. 604-607, October 1974. [27] D. Baum, A. C. Baum, and S. C. Church, "Effect of sedation on

oxygen consumption of children undergoing cardiac catheteriza-

tion," Pediatrics, vol. 39, pp. 891-895, June 1967.

[28] D. B. Dill and A. Folling, "Studies in muscular exercise," J. Physiol., London, Vol. 66, pp. 133-135, 1928. [29] D. K. Gehmlich and S. B. Hammond, Electromechanical Systems, New York: McGraw-Hill Book Co., 1967, pp. 91-100.

Processing Applied to Scintillation Images from Biomedical Systems

Image

PANOS P. VAROUTAS, MEMBER, IEEE, LOUIS R. NARDIZZI, ERNEST M. STOKELY, MEMBER, IEEE Abstract-Several nonlinear filtering methods have been developed for two-dimensional image signals obtained from an Anger Scintillation camera. The first method termed "pruning" is implemented in the frequency domain and filters out certain components of the twodimensional amplitude spectrum. The other methods called "Efiltering" and "P-filtering" utilize a nonlinear transformation of the spatial coordinates which is proportional to the two-dimensional signal amplitude and energy spectrum, respectively. These nonlinear filtering techniques have been applied to phantom scintigrams as well as brain and heart scans of normal and abnormal patients. In addition, the detection of boundaries using the hexagonal Golay transform on filtered and unfiltered images is included. Finally, these filtering and edge detection techniques were for the first time applied to a new technique for the detection of acute myocardial infarcts.

I. INTRODUCTION

SEVERAL WORKERS [1]-[6] have shown that digital

image processing can enhance the diagnostic value of the nuclear medicine scintigrams (i.e., images from a scintillation camera). It has been shown that image processing can imManuscript received October 31, 1975; revised March 26, 1976, and May 24, 1976. This work was supported in part by the National Science Foundation under Grant ENG 73-03789 A02, in part by grants from the Southwestern Medical Foundation, the Southwestern Medical School Institutional Grant Fund, and the Moss Heart Fund, and in part by the National Institutes of Health under Grant HL 13625. P. P. Varoutas was with the Department of Electrical Engineering, Southern Methodist University, Dallas, TX 75275. L. R. Nardizzi is with the Department of Electrical Engineering, Southern Methodist University, Dallas, TX 75275. E. M. Stokely is with the Department of Radiology, University of Texas Health Science Center, Dallas, TX 75235.

MEMBER, IEEE, AND

prove the detectability of lesions in the scintigram both with no change in the false positive probability [6] and in one case with an increase in the false positive rate [7], [8]. Figure 1 shows the anterior view of a myocardial study using a radioisotope which is sequestered by infarcted myocardium. The top panel shows undigitized, unprocessed scintigrams (two exposures of the same image), of the damaged heart muscle which is not detectable and was judged normal by a nuclear medicine specialist. In the bottom panel, following simple background subtraction and contract enhancement, the damaged area is clearly visible and denoted by the arrow. The improvement in lesion detectability following this simple manipulation demonstrates the usefulness of image processing of nuclear medicine scintigrams. This paper investigates two types of nonlinear filters which were applied to scintigraphic 2-dimensional images. The fact that nonlinear filtering techniques may be potentially more useful than linear techniques, even when the image degradation is linear, provides the motivation and rationale for the present investigation. Nonlinear filters also may be computationally faster, may enhance certain salient features of the particular signal being processed, or may for certain image classes be easily implemented. These techniques have been applied to images obtained by an Anger scintillation camera, such as brain or heart scintigrams. The first class of nonlinear filters under investigation is called "pruning" filters which are implemented in the frequency domain. Certain components of the two-dimensional amplitude spectrum, due to noise, are selectively set equal to zero, while keeping the phase spectrum constant. A similar

338

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, JULY 1977

The detection of boundaries has also been investigated using the hexagonal Golay transform. A comparison of boundary detection for filtered and unfiltered images is included. Bounded areas were compared visually and also objectively in order to quantitate the comparative results. Image enhancement and edge detection are, for the first time, applied to a new technique for detection of acute myocardial infarcts. Finally, numerous examples of abnormal heart and brain scans are presented, including the filtered images and bounded areas of infarcted and abnormal tissue.

el Fig. 1. Anterior view of Technetium-99m pyrophosphate myocardial scintigrarn of a patient with a myocardial infarction.

technique was first applied to scintigrams by Kirch and Brown [9]. An inverse image transform generates the filtered result. The second class, termed E-filters [10] [11] utilizes a nonlinear expansion of the spatial coordinates which is proportional to the two-dimensional signal amplitude. A third class of nonlinear filters, termed P-filters, utilize a nonlinear expansion of the spatial coordinates which is proportional to the two-dimensional signal energy. A computational scheme, starting from the image center and resulting in an optimal rectangular grid, has been implemented for both the E- and P-filters. In the examples chosen for this investigation, the filters were not matched to the statistics of the noise (nonlinear Poisson process) presented by the radionuclide disintegration. The only assumption made was that the noise is of "low" amplitude with respect to the original signal. This very broad assumption may be useful for many applications where the noise statistics are not a priori known, and makes the techniques investigated important from both practical and theoretical points of view. In the case of nuclear medicine images, digitized scintigrams usually contain 100-300 counts per image element. Because of the Poisson nature of radioisotope decay, the variance of the "noise" in each image element is equal to the mean of the counts in that area. Thus the standard deviation of the noise component of each image element is VK, where K is the number of counts in the element. Taking the observed counts as an estimate of the mean and variance of the true number of counts, the standard deviation of the noise as a percent of the signal, or true number of counts, is l/rK/. For 100-300 counts per element, the standard deviation of the noise is 6-10% of the true count value; thus, the low amplitude noise assumption holds for scintigrams [1]. ,

,

II. LINEAR AND NONLINEAR IMAGE PROCESSING A. Image Functions An image function for a black-and-white picture can be thought of as real valued nonnegative function of two variables f(x, y). This function is proportional to the light intensity emerging from a picture at the point (x,y). The intensity of a picture at a point is also called the "gray level" or the "brightness." Since the amount of light reaching the observer from a picture is finite and nonnegative, it can be assumed that any image function is bounded and nonnegative. That is, for any image function there is an upper bound K to the possible brightness of any picture and thus for any image function

f(x, y)

O< f(x,y) < K,

V (x,y).

It is reasonable to assume that any image function is integrable since physically the integration corresponds to measuring the total amount of light reaching the observer from the picture. In computer image processing, one represents images as discrete arrays of nonnegative numbers; that is, matrices instead of continuous image functions. Any such matrix can represent a piecewise constant image function and will be called a "digital image function." An example of a 64 X 64 point digital image function would be the brain scans shown in Figure 2(a). This image has been produced on a cathode-ray tube. The white portion of the image represents a gray level of one, while the darkest portion represents a gray level of zero. Any digital picture can have its gray levels scaled between 0 to 1 by subtracting the smallest intensity from all points and then dividing by the intensity range. B. Image Transforms Similar to the case of one-dimensional signals, an image or two-dimensional signal has an equivalent representation in the frequency domain. It has been shown by Andrews [12] that information is neither lost nor gained by going from the spatial to the frequency domain and one may consider the frequency domain as a video scene source without any loss of information compared with the spatial domain. When an image function f (m, n) is represented by a square matrix of N X N points (in our case 64 X 64) then the twodimensional discrete Fourier transform can be defined by equation (1)

F(u, v) =

f(m, n) e(i27/N)(mu+nv) E N-1i-2rN(un

N-i

m=O n=O

(1)

VAROUTAS et al.: DIGITAL IMAGE PROCESSING AND SCINTILLATION IMAGES

or

G(u,v)= F

u-

N,

-

)

(5)

Therefore, the origin moves to the center of the transform domain. The Fourier transform F(u, v) of the image function f(m, n) is in general complex, even though f(m, n) is a real positive function. One can write (6) F(u,v)= F(u,v)Ieio(u, ) where F(u, v) is the magnitude or "amplitude spectrum" of f(m, n) and 0(u, v) is its "phase spectrum." Another important property of F(u, v) is that of conjugate symmetry. That is because f(m, n) is real and the cosine and sine functions are even and odd functions, respectively, such that it can easily be deduced that

F(u, v) = F*(-u, -v)

Fig. 2. (a) Left lateral brainscan with lesion. (b) Log of- the magnitude of Fourier transform of left lateral brainscan. (c) Ph.ase of Fourier transform of left lateral brainscan. (d) Vertex brainsc;in with lesion. (e) Log of the magnitude of Fourier transform of ver rtex bramscan. (f) Phase of Fourier transform of vertex brainscan.

The terms u and v are called the spatial freque ncies of the image. The inverse Fourier transform which will reconstruct the original image is given by N-1 nv rN ( F(u, v)e(127r/N)(m+nv f(mn)= i2 E N-1ii u=0 u=O

(2)

The two-dimensional transform can be computerd as two sequential one-dimensional transforms, because tl ie transform kernels are separable and symmetric. If the Fourier transform of an image is calculateed by the use of equation (1), then the origin or zero spatial fre quency term appears in the corner of the transform plane. In c)rder to get a more natural display, it is convenient to shift the origin to the center of the transform domain. This procedur e can be accomplished by multiplying the image function f( m, n) by the term (- l)1 +' before transforming. Let G (u, v)=

N-i N-i

1:

E

m=O n=O

(- 1)r

+'f(m, n) e( 2rN( nu+nv) (3)

and

(-

l)

+n = e -j7r(m

+n).

Subsequently, equation (3) may be rewritten as N-i N-i

G(u,v)= £ £ f(m,n) m =O n=O

,

e(j2r/N)[n(u-N/2)+n(v-N/2)]

(4)

(7)

where the * symbol denotes the complex conjugate operation. The display medium for Figure 2 was a cathode-ray tube (CRT). The images were read off 9 track magnetic tape onto a PDP-8I minicomputer, which then drove a special gray level CRT. Next, polaroid pictures of the display were taken. The origin of the frequencies in the transform domain, of all image transforms, has been translated to the center of the images by multiplying the original image functions by (- 1 )' +n The logarithm of the magnitude is displayed rather than the magnitude itself in order to nonlinearly expand the gray level contrast of the image. Thus some of the low magnitude components are displayed; whereas, if the pure magnitude was displayed then only the large magnitude components around the center would be observed in the image. Figure 2(a) shows an original scintigram of a left lateral brain-scan. This image is a matrix of 64 X 64 nonnegative terms representing a two-dimensional projection of the threedimensional radioisotope distribution in a human brain. This radiation is detected and measured by a collimated detector in a plane external to the source. Figure 2(b) displays the logarithm of the magnitude of the two-dimensional Fourier transform of the image function in Figure 2(a). Note the highest magnitude is at the center of the picture and the starlike shape and symmetry of the amplitude spectrum. Figure 2(c) displays the phase of the two-dimensional Fourier transform of Figure 2(a). Again note the symmetry of the phase spectrum. Figure 2(d) is a vertex (top) view of a human brain. Again the logarithm of the magnitude and the phase of the twodimensional Fourier transform are shown, respectively, in Figure 2(e) and (f). It is interesting to observe the differences in shape of the amplitude spectra of Figure 2(b) and (c). Each spectrum is different from the other in the four pictures shown (Figure 2(b), (c), (e), and (f)). Note that one obtains image transforms which are unique and invertible representations of an image in the frequency domain. Also, it is interesting to note that there is considerable information at all frequencies along the x axis of Figure 2(e). Higher spatial frequencies in this example result from abrupt gray scale

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, JULY 1977

340

TABLE I GOLAY TRANSFORM PATTERNS U

PATTERN

INDEX

PATTERN

INDEX

1

)

1

1

1 1 1 1 o 01 0 1 0 X X X X X 1 0 0 0 0 0 00 00 0 0 0 0 0

ZERO

ONE

0

1

11 0 X X 0 0 0 1 0 1

1

SEVEN

EIGHT

TWO

1

1

0 X 0 0

1

NINE

THREE

1 0

1 X

1

TEN

1 0

FOUR 0

0

1

X 1

1

1

1

1

1

1 0

0 X

0

1

X 1

1 1

SIX

FIVE

0 0 X

0

1

0

1

1 1 X 1 0

0

1

ELEVEN TWELVE THIRTEEN

changes, such as occur at the edges of an image. Therefore, linear low-pass filtering of this spectrum would distort the edges of the image obtained from the inverse transform of the filtered spectrum. III. GOLAY HEXAGONAL PATTERN TRANSFORM

Golay [13] developed a hexagonal modular transformation to be used in two-dimensional parallel computers. These pattern transforms are position independent local operators useful in transforming a binary image. This two dimensional logic is based on the assumption that any digital image may be represented by a registered stack of planes, usually three, each containing an n X n array of binary image elements. The elements of each plane are generated by thresholding the image at a specified gray level of interest and letting the picture elements darker than this level take the binary value ONE and all other elements are binary ZEROS. Next, the resulting binary image data points are arranged in a hexagonal tessellation (pattern) where each point is surrounded by six nearest neighbors. This arrangement is preferred to the square or checkerboard tessellation in order to avoid the connectivity ambiguity related to square arrays. Instead of having eight neighboring points, i.e., four immediate neighbors at a distance equal to 1 and four secondary neighbors at a distance of If-, each matrix point now has six equidistant nearest neighbors. A point of the hexagonal tessellation with its six nearest binary neighbors forms a Golay pattern. There are 14 orientation-independent, binary, hexagonal, nearest neighbor configurations of a Golay pattern [13]. These patterns are shown in Table I. Each pattern may represent more than one end and up to six configurations by rotating in one sense (clockwise or counterclockwise, but not both). Thus the operations on the digital image elements are isotropic in a sense that they are invariant to rotations by 600, 1200, etc. When an image is transformed into one plane of N X N binary data points, then the following sets may be distinguished among the binary sets of elements. 1. The set of all points in the plane that have the value one. 2. The set of all points in the plane that have the value zero. 3. The set of all points in the plane which have a specified property; i.e., their six nearest neighbors coincide with one or more of the above Golay patterns. 4. The set of all points which do not have the specified property. The intersection of the above sets represents certain operations, such as erasure of the contents of the plane, duplica-

tion, inversion, marking of points in the plane that have the specified six nearest neighbors, contracting the number of points which have the value one and propagating "waves" out from any picture element in the plane that has value one. Many procedures have been written for other types of feature extraction related to images such as those obtained for identification of white blood cells. The procedure followed by the authors to detect the boundaries of brain tumors or heart infarcts was as follows. After collecting the initial scintigraphic data the image was transformed into a binary matrix by setting all the values above a threshold equal to one and the rest equal to zero. The points were then arranged into a hexagonal (rhombic) pattern and an operation was performed at each point that tested its six nearest neighbors and compared it to the fourteen basic Golay patterns. Certain combinations of Golay indices are used to detect, edges. Matrix points found to be on an edge are replaced by one, otherwise they are set equal to zero. Consequently, if the point had an original value of one and its six nearest neighbors had a Golay index not equal to zero, one, or six, then its value remained unchanged. Otherwise, it is erased (set equal to zero). The above technique has proved useful in detecting the ventricle and lesion edges for brain scans as well as areas of infarcted tissue in heart scans. It was chosen primarily for its simplicity and minimum computer storage capabilities. Nevertheless, it is not the only method that could be used, and future investigations concerning image processing and edge detection for radioisotope scanning images is required [14]. IV. ANALYSIS OF THE THREE NONLINEAR FILTERING TECHNIQUES FOR IMAGE PROCESSING A. Pruning Filter The "pruning" filter can be explained from the block diagram of Figure 3. Let f(x,y) be the received signal and noise, i.e.,

f(x,y) = s(x,y) + n(x,y)

(8) where s(x,y) is the original uncorrupted signal and n(x,y), the noise. It is also assumed that the amplitude of the noise n(x,y) is "low" compared to the amplitude of the signal s(x,y), which is very often the case in scintillation camera images [1] as described earlier. The Fourier transform of f(x,y) will consist of the superposition of terms belonging to s (x, y) and n (x, y). Therefore, if one subtracts a constant value or threshold in the range of the amplitude of the noise spectrum from the total discrete amplitude spectrum, then most of the noise corrupting the image could be filtered. This would be ideal for noise with a more or less constant amplitude spectrum. The above method has been termed spectrum "pruning" and consists of the following procedure. The discrete two-dimensional Fourier transform F(u, v) of the digitized image is obtained and then the amplitude F(u, v) and phase O(u, v) spectra are determined. The mean of the magnitude of the Fourier transform termed M is calculated and a multiple (PER) of the mean M is subtracted from each spectral term, while keeping the phase un-

VAROUTAS et al.: DIGITAL IMAGE PROCESSING AND SCINTILLATION IMAGES

341

PRUNING FILTER

IG(u,v)I

f(XY)

FOURIER

TRANS FORM

tFCumv)Ie

(UV)

=

IF(u,v)l

If IG(UV)I < then I G (uv)

-T I G CuV) I

0 -

0

(u

v)

INVERSE

FOURIER TRANSFORM

g(xy)

Fig. 3. Pruning filter.

changed. The "pruned" amplitude spectrum is equal to IG(u,v) I= F(u,v)I- T (9) where T is the pruning threshold and equal to (PER) X (M). Also, if G(u, v) < 0, then I G(u, v) = 0. Next the inverse discrete Fourier transform of G(u, v) is calculated and the filtered image g(x,y) is obtained, as shown in Figure 3. The pruning threshold T to be used in obtaining an enhanced filtered image has been found empirically. Experimental facts indicate that above a certain threshold of pruning the image is oversmoothed and/or distorted. The pruning technique has been applied initially to artificial (phantom) images. A computer generated matrix of 64 X 64 points representing a "square pattern" image was created (Fig. 4(a)). The gray level value in the overlapping portion of the two central squares was set equal to five hundred. The non-overlapping sections have a value of two hundred and fifty. The two smaller squares at the upper left- and lower right-hand corners have a gray level value of one hundred and twenty-five and the rest of the points were set equal to zero. The gray level values were weighted in the above manner in order for the image to have a high and intermediate spatial frequency content. A pseudorandom noise of variable amplitude was added algebraically point by point to the phantom image (Fig. 4(b)). The noise is uniformly distributed between 0-500 with a mean value of 253, almost equal to the value of the non-overlapping portion of the central squares. Figure 4(c) displays the logarithm of the amplitude spectrum of the original "square pattern," and Figure 4(d) displays the logarithm of the amplitude spectrum of the "square pattern" plus "noise" image. In order to quantify some of our conclusions concerning the filtered images, a figure of merit for the filtering methods employed will be defined using a sum square error criteria. The sum square error (SSE) of the brightness function s(x,y) of an image is defined as N-1 N-1

SSE =

E E [S

i=0

i=0

-

(ij)]2

N-1 N-1

E E [s(ij)]2 i=O

j=O

where i, j are the spatial coordinates, s(i, j) is the uncorrupted image brightness value, and f(i, j) is the filtered image brightness value. This performance index was applied to the known values of a square pattern plus pseudorandom noise displayed in Fig. 5 with the result shown in Table II.

It is noted from Table II that pruning with a threshold of 5 times the mean of the amplitude spectrum (PER = 5.0) gives the best sum square error for the image of Figure 5(b) (i.e., square pattern + noise). For the square pattern image plus noise/5, the best result is obtained with a pruning threshold equal to the mean (PER = 1.0). Comparing the results of the first column with the second, it is seen that the threshold values that give the smallest SSE are different. In the case of pruning with PER = 5.0, the sum square error is higher than the SSE of the original image because of the excessive filtering of the high frequencies and distortion of the image (i.e., 136% SSE). In order to assess the pruning filter technique quantitatively, the Golay transform operator was applied to the phantom square pattern images with or without additive pseudorandom noise. The resulting measured areas of unfiltered patterns were plotted versus different percentages of additive noise. Similarly, the measured area for filtered images at a specified level of noise was graphed versus the pruning threshold. The results are tabulated in Table III and plotted in Figure 6. The true or actual area is indicated by a dotted line. The measured area versus noise level curve lies above the true area for amounts of additive noise comparable in amplitude with the original signal; whereas, it shows a decrease in area for smaller amounts of additive noise. The effect of varying the pruning threshold at a specified noise level is indicated on the other three curves for the original amount of noise, 2/3 and 1/5 the amount of additive noise. The full noise curve lies well above the actual area line as expected and changing the pruning level does reduce the error between the measured and the true area. The 2/3 noise curve is above and the 1/5 noise curve is below the true area line. Thus adjustment of the pruning level yields a slight decrease in the area error. The results demonstrate the need of prefiltering before applying a boundary detection algorithm to scintigraphic data when the additive noise has low amplitude with respect to the original signal. This study for the particular phantom patterns examined shows that the difference between the measured area and the actual one is not significantly decreased by prefiltering for large magnitudes of additive noise. However, prefiltering does reduce the area error for additive noise in the range of the mean of the original signal. Adjustment in the percentage of pruning (changing threshold levels) does affect the measured area for amplitudes of additive noise comparable to the original signal, while the threshold change exhibits small sensitivity for moderate amounts of additive noise. The purpose of obtaining improved estimates of areas of damaged tissue, especially for myocardial infarcts, is that the area is a

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, JULY 1977

342

TABLE 11 SUM SQUARE ERROR (SSE)

Square Pattern + Noise

IMAGE

X

Noise/5

+

Original

2.595

100

0.1038

100

Pruned (PER = 1.0)

2.042

79

0.0889

86

Pruned (PER = 5.0)

1.969

76

0.1412

136

AREA

Fig. 4. (a) Original square pattern. (b) Square pattern with added pseudorandom noise. (c) Logarithm of the amplitude spectrum of the square pattern. (d) Logarithm of the amplitude spectrum of the pattern plus noise.

Square Pattern

.

OF

TABLE III DETECTED BOUNDARY BY GOLAY TRANSFORM ALGORITHM Unfiltered

TWO-DIfENSIONAL SIGNAL

In2

Original square pattern

1.3325

Square pattern + noise

1.3615 -

Square pattern + noise/2

PRUNING-IMAGE

FILTERED BY

Image Area

TYPE

PER=1 PER=1/2 --

1

AREA IN

PER=1/3 PER=1/4

_

_

_

1.6125

1.4050

1.4025

1.4020

1.3350

1.3425

1.3400

1.3400

1.3115

__

__

Square pattern + noise/3

1.2765

__

Square pattern + noise/4

1.2650

__

Square pattern + noise/5

1.2650

1.2980

1.3020

Square pattern + noise/1.5

_

_

__

__

__

__

__ 1.3030

__

1.3030

1.6000 .

1.5500

ok

E

Measured Area vs Noise Level

a

Measured Area vs Pruning Threshold level

._

I.,

81

Fig. 5. (a) Boundary detection after square pattern and noise. (b) Boundary detection of filtered by pruning (PER = 0.5) square pattern and noise. (c) Boundary detection of square pattern and noise/S. (d) Boundary detection after filtered by pruning (PER = 0.5) square pattern and noise/5.

4 SQRPAT+Nois. 1.4000 F

valuable index for estimating total infarct weight [6] and thus has important diagnostic and prognostic implications. B. Two-Dimensional E-Filter The two-dimensional E-filter may be represented in block diagram form as shown in Figure 7(a). A transformation Tis defined on the spatial variables x and y to the new variables d and e and the history of the input function f( ,) up to the spatial coordinates (xl, y 1)

Tj f(x,y)} g(d, e). _

Referring to Figure 7(a) again, the output of the transform block T, g(d, e), is filtered by a linear filter H to produce

SQRPAT+Noine/1.5

Aetual Area

SQRPAT+Noise/5 1.3000 t vSQRPAT+Nose/.

1

2

3

4

5

1

Fig. 6. Noise or pruning threshold level in multiples (PER) of the mean of the amplitude spectrum.

VAROUTAS et al.: DIGITAL IMAGE PROCESSING AND SCINTILLATION IMAGES

SPACE

343

IE-Doman

SPACE

(a)

x

(b)

Fig. 7. (a) Block diagram representation of the two-dimensional E-filter. (b) Two-dimensional E-transform implementation of the function f(x,y).

k(d, e) as its output.

The inverse transformation T-1 is now applied to g(d, e) to transfer back into the spatial domain giving u(x, y). The overall transformation is designed as (. For the particular case considered in this investigation, the transformation of T is defined by letting g(d1, el) be equal to the value of f(xi, y1) at the spatial coordinates (d1, el) such that the area of the rectangle d1 el is equal to the surface area of the two variable function f(x, y), which is between the x- and y-axis planes and the planes xl = 0 andy1 = 0 as shown in Figure 7(b). Observe that the area of the surface f(x,y) which is bounded by the x- and y-axis planes and the vertical planes to the x- and y-axis, passing through the coordinates (xI ,y') is expanded (stretched) on a plane such that it fills a dl el rectangle of equal area. The value of the transformed function g(d1, el), though, is set equal to the original value of the function at (xl,yl). That is, the transformation T is the identity transformation as far as the value of the function is concerned. In mathematical terms, this transformation is defined as

d1.el~

f 0

~

1+(iaL) +1 ay ~ ax

JLx/1+1II Vf 112 dxdy and de AY el

yl

where I Vf i is the norm of the gradient of the function f(x,y). The value of the transformed function g(d, e) is set to the original function value at (x, y), i.e.,

g(dl, el) =f(xl,yl).

(13)

The digital realization of the two-dimensional E-filter involves the summation of elemental surface areas whose projections are termed elemental picture cells. An elemental picture cell is bounded by four adjacent sampled picture points forming a square. In order to calculate the elemental surface areas, the value of the norm of the gradient 1l Vf 11 is needed at the center of each picture cell. After the Etransformed function g(d, e) is defined (Figure 7(a)), then the linear filter H may be implemented, for example, by a bicubic spline interpolation scheme [151]. If the two orthogonal directions, along which the directional derivatives are known, are selected to be the lines whose slopes are plus and minus one, then one may approximate the gradient as

Vf (k, 1)11 R(k, 1) = { [f(k, 1)- f(k + 1,1+ 1)]2 + [f(kl+ 1)- f(k+ 1,)12}1/2

dxdy (1 1)

(14)

This can be demonstrated diagramatically if one considers at cell (k, 1) a 2 X 2 window whose diagonal elements are subtracted in the manner shown in Figure 8(a). The directional (12) derivative in each direction is approximated by simply subtracting adjacent elements. The operator R(k, 1) is sometimes

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, JULY 1977

344 f

f

k+1,

k+, Z+1

0

1

kE

k, Z+

1

0

(a)

1

(b)

0

0

1 (C)

Fig. 8.

referred to in the literature as a "Roberts cross operator" [16]. If the point (k, 1) is in a region of uniform intensity (constant), then the value of R(k, 1) is zero as one expects it to be. If there is a discontinuity in intensity between column 1 and column (1 + 1), and similarly if there is a discontinuity between row k and row (k + 1), then R(k, 1) has a large value. Equation (14) can be used to reasonably approximate the elemental surface area since it represents a discrete approximation to the gradient at the center of the elemental square. The corresponding masks for the directional derivatives f, and fy are shown in Figure 8(b) and (c), respectively. The surface area Ai of an elemental cell of the digitized picture using equation (14) will then be equal to (15) Ai = [1 + 0.5 11 Vf (i,j) 112] 1/2. Note that the magnitude of the gradient has been normalized since the distance between the diagonal elements of the elemental cell is equal to I/2. The E-transformed coordinates may be calculated now, as indicated in Figure 9. Starting from the center of the digitized image, one calculates first the surface area bounded by the center most elemental cell. The value of this area defines the first central square of the E-transformed grid. Next, starting from the center and along the horizontal direction, one adds the areas defined by the cells to the right and to the left of the central cell. Each new area defines another four points, two on each side of the e axis. This process is continued until the first horizontal strip of elemental cells bisected by the x axis is covered. The summing of areas is repeated starting from the center with the three cells nearest to the x axis until the threestrips are covered and so on with the next 5, 7, * -, N - 1 horizontal strips. The original grid is shown in Figure 9(a) and the E-transformed grid in Figure 9(b). Observe that the points that fall on the 450 and 1350 lines are symmetric about the origin and the d and e axes. Also points that are in the sectors between the 450 and 1350 lines are symmetric about the d and e axes. Thus this method will present minimal distortion for image scenes that are located at the central portion of the picture square boundary and are more or less circular in shape. Since the origin of the x-y coordinate system has been translated to the center of the image in this method, the limits of integration in equation (1 1) will have to change such that -

IJ

dlel =4

X1/2 x I/

rYI/2

J

V1

/

1Vidx dy

IJf 12

(16)

a,/

the coordinates in the translated d and e (di, ei) axes, as shown in Figure 9(b). If one now calculates the mean of all coordinate points that lie on the same gridline curve of Figure 9(b), the nonlinear grid where

are

(a)

I

I

(b) Fig. 9. (a) Original square grid. (b) Nonlinear E-transformed grid starting from the image center.

is transformed into an optimum rectangular grid. That is, instead of using the points along the 450 line as bench marks for constructing the rectangular grid, an optimum grid is obtained if one fits a straight line parallel to the d or e axis through the points lying on the same gridline curve, such that the distances of the points from the straight line have a least squares error. For example, one would like to pass through the coordinates di n, d2n, d3ns * * * , dmn in a straight line which is described by the equation d = c, and which minimizes the quantity m

I=y, i=l

(c - din)2

(17)

with respect to c. Setting (alac) = 0, one obtains mc

m

=E di i=i

=

0

(18)

and C=

m i=1

din/m

(19)

Each straight line is parallel to the e axis and has an intercept with the d axis equal to the mean of the coordinates din, d2l, * * *, dmn (and vice versa for the d axis). In this way, one transforms the nonlinear grid of Figure 9(b) into an optimum rectangular grid.

VAROUTAS et al.: DIGITAL IMAGE PROCESSING AND SCINTILLATION IMAGES

C Two-Dimensional P-Filters A different nonlinear transformation in the class of the Efilters could be defined by "stretching" the spatial coordinates proportional to the energy of the signal. This class of nonlinear filters will be called P-filters and can be defined as follows for the two-dimensional case. Let

f(xy) = f(x + 'Jx,Y + Ty)

(20)

be a periodic function with spatial periods of rT and ry in the x and y axes, respectively. A nonlinear P-transformation of the spatial coordinates is defined as

rxiY1

d1 el

f 2(x,y) dx dy

(21)

345

becomes

fc >Max E IT

P

'PP]

(30)

The terms r, P and ryP are proportional to the root mean square energy of the signal. It is obvious that the P-filter will process signals in a way which depends on the energy content of the signal. Signals with "high" energy will pass unattenuated, while "low" energy signals will be attenuated by the Pfilter. This type of nonlinear filter would be useful when one wants to separate two signals that have been superimposed by addition, but ate of different energy levels.

V. EXAMPLES OF BIOMEDICAL SCINITIGRAMS A new method for visualizing acute myocardial infarcts in and humans has been tested at the Department of Radiology, Uni(22) versity of Texas Health Science Center in Dallas [18] . Patients with chest pain or suspected myocardial infarcts were given an el Yi Note that equation (21) represents the energy of the two- intravenous injection- of 15 mCi of Technetium tagged to 5 dimensional signal in the area defined by the x and y axes and mg of stannous pyrophosphate (99 mTc-PYP). One hour after the x1 = 0, Yi = 0 lines. The periods Td and Te of the P- injection the distribution of radioactivity in the patient is imaged with a Nuclear-Chicago HP scintillation camera. Imtransformed signal will be from equation (21) equal to ages are collected in the anterior, lateral, and one or more left TOX rTY anterior oblique projections using a 16 000 hole high resoluf f2(x,y) dx dy (23) tion Td . Te= f collimator. Images were obtained from 12 hours to 14 days after the onset of chest pain. Positive scintigrams in the and using equation (22) patients with myocardial infarcts are thought to be due to inof pyrophosphate into the crystalline structure of corporation Td = TX (24) the hydroxyapatite found within the organelles (mitochondria) Te TY of irreversibly damaged myocardial cells [ 16]. Parseval's relationship yields In order to make a qualitative comparison between the Efilters and P-filters as applied to biomedical image processing, 00 1 rtx Iry '2 some myocardial scintigrams are presented. Figures 10(a) and t (25) (X,y) dx dy m=-oo A, A, TxTy L 11(a) represent original heartscan images obtained by an Anger n=-oo where the Cmn represent the two-dimensional Fourier series scintillation camera. Figure 10(b) is a two-dimensional Ecoefficients corresponding to an expansion of the function filtered image of the scintigram shown in 10(a), using an optimum rectangular grid as described in Section III(B). f(x,y). The term Figure 10(c) and (d) demonstrate the hexagonal Golay transform edge detection scheme applied to Figure 10(a) and (b), X Z (26) pA m=-oo E ICmnI respectively. There are three important structures to note in ln=-oo Figure 10(c) and (d). The large vertical crescent shaped strucrepresents the average power of the two-dimensional signal ture in each figure represents the sternum. The other two f (x, y). It can be deduced using equations (23), (24), (25), structures indicated by a large and small arrow represent two and (26) that the periods 'Td and re of the P-transformed sig- myocarcial infarcts with different times of occurrence. The nal are equal to larger infarct occurred earlier in time. The filtered image of =d Tx P (27) 10(b) clearly shows separation of the large infarct and sternum, well as a more discernible small infarct and reduced inte = Ty P. (28) as tensity at the center of the large infarct when compared to Referring to the block diagram representation of the two- Figure 10(a). The reduced intensity represents pathologically dimensional E-filters (Figure 7(a)), the P-transformed periodic a necrotic center with decreased uptake of tracer agents. signal with spatial periods rd and 'Te will pass unattenuated Finally, the rib structure and background noise is filtered out through a linear circularly symmetric filter H with cutoff fre- in Figure 10(b). quency f, provided that it satisfies the condition Figure 1 1(a) represents an original left lateral view of an abnormal heartscan. E-filtering and P-filtering of Figure l (a) f > Max[ j (29) are presented in Figure 1 (b) and (c), respectively. The hexagonal Golay transforms are shown in Figure 1 (d), (e), and -] where the expression Max [ , means the greater of the two (f). A smaller and separable infarcted area is discernible in arguments. The above condition from equations (27) and (28) Figure 1 1(e) and (f) by the smaller arrow and the crescent

IICmn

00

00

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, JULY 1977

346

(a)

(c)

(b)

(d)

(b) Two-dimensional E-filtering of heartscan in (a). (c) Boundary detection of original heartscan in (a). (d) Boundary detection of E-filtered heartscan in (a).

Fig. 10. (a) Original heartscan.

shape representing the sternum. It is important to note here that the P-filter enlarged and distorted the original heartscan while at the same time improving the contrast between areas of uptake and no uptake. These preliminary results would 11. (a) Original left lateral view of a heartscan. (b) Twosuggest that the P-filter would not be a suitable filtering Fig. dimensional E-filtering of (a). (c) Two-dimensional P-filtering of (a). scheme to use for the heart scintigrams described herein be(d) Boundary detection of original heartscan in (a). (e) Boundary detection of E-filtered heartscan in (b). (f) Boundary detection of cause the energy levels of the signal and noise are not sigP-filtered heartscan in (c). nificantly different to be properly filtered. In contrast, the E-filtered images of Figures 10(b) and 11(b) demonstrate considerable improvement over their original counterparts and This very broad assumption may be useful for some applicasuggest that applying E-filters to original scintigrams can be tions where the noise statistics are not known a priori. The detection of boundaries has also been investigated using useful in interpreting scintigrams which might otherwise be the hexagonal Golay transform. A comparison of boundary read incorrectly. detection for filtered and unfiltered images was included. Bounded areas were compared visually and also objectively in VI. CONCLUSIONS order to quantitate the comparative results. Image enhanceSeveral nonlinear filtering methods have been developed for edge detection were, for the first time, applied to a ment and two-dimensional signals. The first of these methods was for detection of acute myocardial infarcts. new technique termed "pruning" and is implemented in the frequency doexamples of abnormal heart- and brainscans Finally, numerous main. Certain components of the two-dimensional amplitude the filtered images and bounded including were presented, spectrum, due to noise are selectively set equal to zero, while tissue. and abnormal infarcted areas of keeping the phase spectrum constant. An inverse image transform generates the filtered result. The second and third REFERENCES methods called E-filtering and P-filtering utilize a nonlinear transformation of the spatial coordinates which is proportional [1] Metz, C. E., "A Mathematical Investigation of Radioisotope Scan Image Processing," Doctoral Thesis, University of Pennsylvania to the two-dimensional signal amplitude, or energy spectrum. (Ann Arbor, Michigan, University Microfilms, Pub. #70-16, 186), These filters were realized digitally by starting the transforma1969. tion from the image center and using an optimum rectangular [2] Brown, D. W., Kirch, D. L., et al., "Computer Processing of Scans Using Fourier and Other Transforms,"J. Nucl. Med., 12:287-291, grid structure. The scheme presents reduced distortion for 1971. images that are approximately symmetric about their center. [3] Pizer, S. M., Correia, J. A., Chesler, D. A., and Metz, C. E., "ReThe nonlinear filtering techniques have been applied to sults of Nonlinear and Nonstationary Image Processing,"Information Processing in Scintigraphy, U. S. Energy Research and Dephantom images as well as brain- and heartscan images obvelopment Administration, Technical Information Center, Conf. tained by an Anger scintillation camera. The only assumption 730687,pp.93-101,1973. concerning the two-dimensional signals was that the noise is of [4] Helstrom, C. W., "Image Restoration by the Method of Least Squares,"J. Opt. Soc. Am., 57:297,1967. relatively "low" amplitude with respect to the original signal.

347

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. BME-24, NO. 4, JULY 1977

[51 Tanaka, E., linuma, T. A., "Approaches to Optimal Data Process-

ing in Radioisotope Imaging," Phys. Med. Biol., 15 :603-694, 1970. [6] Knowles, L. G., Kohlenstein, L. C., Schulz, A. G., "Observer Performance as a Measure of 'Goodness' in the Evaluation of Scintigraphic Image Processing," Information Processing in Scintigraphy, U. S. Energy Research and Development Administration, Technical Information Center, Conf.-730687, pp. 15 3-158, 1973. [7] Kuhl, D. E., Sanders, T. D., Edwards, R. Q., and Makler, P. T., "Failure to Improve Observer Performance with Scan Smoothing,"J. Nucl. Med., 13:752-757, 1972. [81 Metz, C. E., Goodenough, D. J., "Failure to Improve Observer Performance with Smoothing: A Rebuttal," J. Nucl. Med.,

14:873-876, 1973.

[9] Kirch, D. L., Brown, D. W., "Nonlinear Frequency Domain Techniques for Enhancement of Radionuclide Images," Information Processing in Scintigraphy, U. S. Energy Research and Development Administration, Technical Information Center, Conf.730687, 102-114, 1973. [10] Moore, D. J. H. and Parker, D. J., "On Nonlinear Filters Involving Transformation of the Time Variable," IEEE Trans. Information Theory (July 1973), Vol. IT-19, No. 4, pp. 415-422.

[11] Varoutas, P., Nardizzi, [12]

[131 [141 [15]

[161 [17]

[18]

L. and Stokely, E., "Nonlinear Filtering and Boundary Detection for Digital Image Processing," 28th ACEMB, New Orleans, Louisiana, September 1975. Andrews, H. C., "Entropy Considerations in the Frequency Domain," Proceedings of the IEEE, January 1968, pp. 113-114. Golay, M. J. E., "Hexagonal Parallel Pattern Transformations," IEEE Trans. on Computers (August 1969), Vol. C-18, No. 8, pp. 7 33-740. King, R. E., "Digital Image Processing in Radioisotope Scanning," IEEE Transactions on Biomedical Engineering (September 1974), Vol. BME-21, No. 5, pp. 414-416. DeBoor, C., "Bicubic Spline Interpolation," Journal of Mathematics and Physics, Vol. 41, pp. 212-218, 1962. Roberts, L. G., "Machine Perception of Three Dimensional Solids," Optical and Electrooptical Information Processing, by J. J. Tippett et al., M.I.T. Press,Cambridge, Mass., pp. 159-197. Stokely, E., Buja, L. M., Lewis, S. E., Parkey, R. W., et al., "Measurement of Acute Myocardial Infarcts in Dogs with Technetium-99 Stannous Pyrophosphate Scintigrams," J. Nucl. Med., Vol. 17, pp. 1-5, January 1976. Bonte, F. J., Parkey, R. W., Graham, K. D., et al., "A New Method for Radionuclide Imaging of Myocardial Infarcts," Radiology, Vol. 110, pp. 473-474, January 1974.

A Computerized Breast Thermographic Interpreter MICHAEL NEGIN, MEMBER, IEEE, MARVIN C. ZISKIN, SENIOR MEMBER, AND MARC S. LAPAYOWKER

IEEE,

CHARLES PINER,

from breast cancer has remained essentially constant at 22 per 100,000 female population for the past 40 years [3], [47]. The incidence of this disease increases with age. It is rarely seen before age 30, but its incidence rises rapidly following the menopause. The 5 year survival rate depends on the size of the cancer performance. and metastatic spread existing at the time of diagnosis. Breast The CBTI thus developed has been embedded in a minicomputer cancer found by self-examination results in a 5 year survival (Linc-8) and requires 5 min of minicomputer time to interpret one rate of 45-55% [15], [47]. However, it has been shown that case and issue its findings on a line printer. the use .of methods which aid in the early detection of breast cancer can yield five year survival rates of 80-85% [1, 15]. INTRODUCTION These findings have stimulated research and development of BREAST CANCER is the leading cause of death by cancer techniques which can increase the probability of early in women. In the United States alone, it kills more than detection. 27,000 women each year. One out of every 15 women will at some time develop breast cancer and in spite of new innova- Thermography tions in surgery, radiotherapy, and diagnosis, the death rate Thermography is the science of recording the infrared radiation emitted from the surface of a patient. A thermogram is Manuscript received January 12, 1976; revised March 1, 1976, and therefore a map of the skin temperature of an object or April 9, 1976. This work was supported in part by the National Insti- patient. The use of breast thermograms to detect breast tutes of Health under Grants GM 14548 and GM 54466. M. Negin is with the Biomedical Engineering and Science Program, cancer has been studied by many groups and has been found Department of Electrical Engineering, Drexel University, Philadelphia, to be an imaging modality which can play an important role in PA 19104. screening programs. [13-17], [19], [25], [45]-[46], [48]. M. C. Ziskin, C. Piner, and M. S. Lapayowker are with the Diagnostic Radiology Research Laboratory, School of Medicine, Temple Univer- Although breast thermography has a useful role in screening, sity, Philadelphia, PA 19140. it must be emphasized that breast thermography should never

Abstract-A computerized breast thermographic interpreter (CBTI) has been developed which outperforms human thermographers on a series of 81 frontal thermograms. The CBTI is composed of feature extraction algorithms followed by linear discriminant classifiers. Human interpretatiois and computer interpretations were compared to a thermographic panel decision and to biopsy findings. In each case, the CBTI performance was better than individual human