Wigner distribution representation of digital images

Pattern Recognition Letters 5 (1987) 215-221 North-Holland

March 1987

Wigner distribution representation of digital images G . CRISTOBAL*, J . BESCOS, J . SANTAMARIA and J . MONTES* Instituto de Optica (C.S.I .C.), Serrano 121,28006 Madrid, Spain

Received 17 March 1986 Abstract : The Wigner Distribution (WD) for discrete images is computed for different test images with gray level and spatial frequency information contents . The most relevant characteristics of the Wigner Distribution are analyzed from 2-D displays of the original 4-D distribution of the test images . This representation through the WD is shown to be specialty adequate for processing of textured information and of spatially variant degraded images .

Key words: Image processing, Wigner distribution, local spectrum representation .

1 . Introduction The Wigner Distribution (WD) introduced by Wigner [I] in 1932 as a phase space representation in quantum mechanics, gives a simultaneous representation of a signal in space and spatial-frequency variables . It belongs to a class of space-frequency functions that includes Woodward's [2] ambiguity function and Rihaczek's [31 complex energy density function, related to each other by linear transformations . The Wigner distribution might be interpreted as a local or regional spatial-frequency representation of an image . However, it presents two main advantages with respect to the local spectrum representation . First, the WD is a real valued function and encodes directly the Fourier phase information . Second, the election of the appropriate window size which depends upon the kind of the analyzed information, is not required for the computation of the WD . On the other hand, the WD of a 2-D image is a 4-D function that involves Fourier transformation for every pixel of the original image . As a con* Facultad de Informatica, Cra . de Valencia Km .7, 28031 Madrid, Spain .

sequence the calculation of the WD of discrete images requires an important amount of computations that can limit, in principle, the range of applications . The problem can be mitigated by using optical processors to obtain the WD . In this way, several optical methods have been proposed recently to produce the WD of 1-D [4,5] and 2-D [6-8] signals, and different optical visualizations have been shown to display such a four-dimensional function . Bamler and Glunder [6] display this 4-D function as two-dimensional frequency slides ; Easton et al . [7] use the Radon transformation to produce the Fourier transform, and display the WD as one-dimensional slices of the 4-D WD . On the other hand, Jacobson and Wechler [8] have proposed conformal mapping techniques to obtain invariant image representations . In this paper we compute the WD for discrete images . First, the distribution for different test images with spatial frequency and gray level information contents is computed, and the most important features of the WD are outlined . Afterwards the WD is calculated for 2-D discrete images that include textured information and spatial degradations . In general, the WD is displayed using a joint

0167-8655/87/53 .50 's? 1987, Elsevier Science Publishers B .V . (North-Holland)

215



Volume 5, Number 3

PATTERN RECOGNITION LETTERS

diagram, spatial frequency (y-axis) versus the spatial variable (x-axis), in order to obtain the variation of the transform along the image . Also, the autoconvolution of the input image in one direction is directly displayed on the x-axis . In the first section the most significant properties of the WD for image processing are included . In the second section the computation and representation characteristics for discrete images are considered, and finally the results obtained in the computation are shown and analyzed in the third section .

2 . The Wigner distribution

March 1987

tion . The integration of Wf over the spatial variable at a fixed frequency gives the energy spectral density at that frequency and the integration of Wf over the frequency variable at a fixed spatial point gives the local power at that point : 1

+m

2n

(3)

I _m Wf(x,co)d(a=[f(x)12,

+m

Wf(x, co)dx=~F(w)I 2.

(iii) Space and frequency shift . Shifts in space and frequency domains give corresponding shifts in the distribution : W8 (x,w)=Wj(x-xo,co),

g(X) =f(x-xo) In this section we consider the auto-Wigner distribution function of 1-D signals for notational simplicity, being the extension to 2-D signals straightforward . The WD Wf (x, w) of the signal f(x) is given by :

(4)

(5) . Wg (x,w)=Wj(x,w-w0) . G(w)=F(w-wo) (6) (iv) Evenness. The WD of a real signal is even in frequency :

Wf (x, w) m

=f

f x+

2

f* x

2 e

i

fo dxo .

(1)

The auto-Wigner distribution gives a generalized autoconvolution of non-zero frequency [9] . There exists an alternative definition of the frequency domain : w WAX,w)= 2n J F(w+ 2 1 . .odwo, .F* w_ w 2)e

(2)

A complete set of properties of the WD was formulated by Claseen and Mecklenbrauker [10] . The following are the most relevants for image processing applications : (i) Realness. The WD of any real or complex function is real as Fourier transforms of the Hermitian function f(x+xo /2)f *(x-x0 /2) are always real . However, it is not possible to interpret this distribution as a density energy distribution, because in general the WD is not always positive . This agrees with Heisenberg's uncertainty principle which prohibits sharp space discrimination with a sharp frequency discrimination . (ii) Space integration and frequency integra216

Wf (x, w) = Wf (x, - (0) •

(7)

(v) Moyal's formula . Moyal's formula [11] is related with the product of two WD's .

+m

1 " 27rI ~ J

Wj(xw)dxdw=l(x) j ' •

(8)

(vi) Interference property [121 . If we have a function f(x) for which Wf is large and positive around two points (x l ,w l ) and (x2 ,w 2 ), then we can expect I Wf I to be large around the point (2i (x1 +x2), L(w, + cu,)) : m

4I

+m

I

Wf (2a-x,2b-w)Wf (x,w)dxdw

m =(W1(a, b))2 •

( 9)

(vii) Autoconvolution . Wf(x,0)=I wf(x+ Z)

f*

dxo, \x- 2/ (10)

WF (0,w)= I

F

w+ 2 )F*(w)dw o . 2 ( / \ (I1)



Volume 5 . Number 3

PATTERN RECOGNITION LETTERS

(viii) Inversion [10] . The signal f(x) can be recovered from the WD by inverse Fourier transformation, up to a constant factor : fxf"(0)=Wf (

2 AA

x ,w ei- dw .

(12)

(ix) Product and convolution [10] . The WD of a convolution of two signals f and g is equal to the convolution W1 sW. i n the spatial variable : W,g (x,w)=I m Wj C o ,w)Wg (x-xo ,co)dxo . (13) The WD of a product of two signals f and g is equal to the convolution Wj *We in the frequency variable : IV) ;' (X, =2n

w )

lam

W (x,co o ) WW (x,w-wo)d coo .

3 . Computation and representation Suppose a function f(x) is space-limited in the interval [0,N-1] . The discrete space frequency WD can be obtained from : f(n+k)f*(n-k)e-'' x' iln-u W1 (n, ,n)=2 k

=2FFT[f(n+k)f 4 (n-k)] .

(15)

If we have N samples of the function f(x), the WD in expression (15) is completed only from the knowledge of even or odd samples separately . Therefore aliasing will occur unless it is overruled by any of the two following solutions : oversampling or interpolating . We have used the first solution for its simplicity in digital implementations . Therefore, considering n samples in the spatial domain, in in the frequency domain, an alternative definition of the discrete WD is given by : W'1 (2n,2w) = C f(n+k)f'(n-k)e

a

March 1987

Two problems could arise when we compute this function in practice : aliasing and leakage . Aliasing is produced in the presence of high resolution objects when they are sampled according to (16) . It can be reduced by lowpass filtering that can be performed by means of convolution in the spatial domain [13] . To reduce leakage it is necessary to employ a spatial apodization or truncation function in order to minimize the side lobes of the sine function [14] . An appropriate truncation function is the Harming window, given by 2nx

1-cos N-1 2 , f(X) =

0<_ .xSN-1 .

(17)

We have used both solutions to avoid the mentioned problems for the computed cases in which they were present . The procedure for the computation of expression (16) is described in the following sequence : (a) Shift of the function f(k)n samples to obtain f(n-k) . (b) Rotation of the function f(n+k) by 180° to obtain fit-k) ; (c) Product of the two functions ; and (d) FFT of the product . Steps (a), (b) and (c) have been implemented in a Vicom Digital Image Processor through the own Vicom commands, and (d) through a Pascal program written for the Motorola 68000 microprocessor of Vicom . Also, as it has been already mentioned, a reduced version of the W D is displayed, in which we only consider shifts in a fixed direction (x-direction) and display the local spectra corresponding to that direction in the v-axis . The simultaneous representation in space and spatial frequencies is obtained by means of successive shifts in the selected direction . The WD magnitude in each point has been normalized to the maximum of the whole WD, and a logarithmic representation has been used in the display .

4 . Experimental results (16)

We ha'e computed the WD of different test 2t7



Volume

Number 3

PATTERN RECOGNITION LETTERS

March 1987

b Figure L (a) Sinusoidal grating object of 64 pixels/period and 0-I gray level range ; (b) W D of Figure Ia for the x-direction, and autoconvolution profile (m=0), with zoom 2x .

images in order to analyze the distribution charac-

range, and its corresponding WD with zoom 2x .

teristics . First, sinusoidal gratings of diverse spatial

In this case, ech main spatial frequency associated

frequencies have been considered . Figure I a shows a 256 x 256 sinusoidal grating with 64 pixels/period

with the different regions of the object appears in

and 0-1 gray level range . Figure lb shows, with

zones of the two regions, while the two main fre-

zoom 2 x, the corresponding WD for the x-direc-

quencies are present in the W D of the centre of the

tion and with the local spectra displayed along the y-axis . The autoconvolution of the original grating

object, as expected . The profile for the autoconvolution also reflects the periodicity, in accordance

is displayed along the %7-axis (w=0), with the varia-

with that of the original object . In this case, this

tion of intensity represented on the profile curve .

multicomponent signal also shows the presence of

The WW' D is also periodic with period one half of the

armonics following property (vi) (Section 2) .

the WD . They are dominant mainly in the central

original grating object . For w=0 an important

Figures 3a and 3b show a rectangular window

variation for the autoconvolution function can be

test and its corresponding WD . In this case, the

observed . Along they-axis not only the fundamen-

WD can be seen as a sine's function composition

tal spatial frequency of the grating appears, but also harmonics corresponding to a (sin) 2 signal .

centre to the border (image product FFT for dif-

Figures 2a and 2b show a composite sinusoidal grating with 32 and 16 pixels/period, 0-1 gray level

where the zeros are shifted as we move from the ferent widths of the rectangle test) . Figure 4 shows a gray level test image, and its

a Figure 2 . (a) Composite sinusoidal crating of 32 and 16 pixels/period and 0-I gray level range ; lb) W D of Figure 2a for the .x-direction, and autoconvolution profile (w=0), with zoom 2x . 218

Volume

Figure

5,

3.

Number 3

(a) Rectangular window test ;

PATTERN

(b) WD

RECOGNITION

of Figure 3a for the x-direction . and autoconvolution profile (w=0), with zoom

corresponding W D with 2 x zoom . The presence of several steps of gray levels also produced distributions of sine's along the frequency axis, in accordance with the sizes of the image products for each object point . The asymmetry and coarseness of this test is also reflected on the autoconvolution, that shows an intensity variation smoother than those of the previous sinusoidal objects . Afterwards, the WD for conventional discrete images has been calculated . Figure 5a shows an aerial textured image (1 .2 x zoom) with two main tree species, 'Eucaliptus globulus' and 'Pinus Pinca', and Figure 5b shows its corresponding Fourier transform . Figure 5c shows on the top the image products at two representative points : -72 pixels from the centre (pinus region), and +72 pixels (eucalyptus region), and on the bottom the WD of the two points . The textured structure of both species is observed in the product images, and

(a)

Discrete gray level test image ; (b)

March

LETTERS

WD of

193?

2 x.

therefore it will also be kept in the WD of the two points . The display of the joint diagram for the WD in .v-direction is shown in Figure 5d . In the WD a more important contribution of high frequencies in the pinus region can be observed in contrast with that of the eucalyptus region . Finally, the WD for spatially variant degraded objects has been computed . Figure 6a shows a grating in focus on the right side, and progressively defocused in the left . Figure 6b shows the image products and WD at points x=-32, x=32 from the center, and their corresponding WD . The WD reflects the different frequency range for the object regions due to the spatially variant defocusing . The same behaviour is displayed in the representation of Figure 6c corresponding to the WD along the xdirection, in which the contribution of high frequencies is lost on the left region .

Figure 4a for the x-direction, and autoconvolution profile

(w=0),

with zoom

2s 219



Volume 5, Number 3

PATTERN RECOGNITION LETTERS

March 1987

.1 AAA _,

b

C Figure 5 . (a) Aerial test image with two main tree species : pinus (left) and eucalyptus (right) with zoom 1.2 x ; (b) Magnitude of the Fourier Transform (logarithmic representation) ; (e) Image products at the point x=-72 pixels of the pinus region (above left), and at the point c=-722 eucalyptus region (above right), and WD's in the two considered points (bottom) ; (d) WD's along the x-direction with zoom 2x .

5.

Conclusions

Advisory Commision for Scientific and Technical Research under Grant

The

WD

for discrete images with diverse infor-

mation contents of gray levels and spatial frequencies has been computed for analysis purposes . The

WD

We also thank Ana

encodes the main characteristics of the images

including the spectral local variation . From the computation of the

WD

for

2D discrete

References

images two

applications seems more appropriate : texture analysis and segmentation, and spatially variant filterin2 . Further research is been pursued in this direction, associated to a hybrid optical-digital system to reduce the great computational effort that requires this kind of representation .

Acknowledgements This work was partially supported by the Spanish 2 20

No . 358 .

Plaza for her collaboration in the general development of the work .

[I1 Wigner, E . (1932) . On the quantum correction for thermodinamic equilibrium . Phvs . Rev . 40, 749-759 . [2] Woodward, .M P . (1953) . Probability and Information Theory, with Applications to Radar . Pergamon, London . [3] Rihaczek, A . W . (1968) . Signal energy distributions in time and frequency . IEEE Trans . Inform . Theory 14, 369-374 . [41 Bartelt, H ., R .H . Brenner and A .W, Lehmann (1980) . The Wigner distribution function and its optical production . Opt. Comm . 32, 32-38 . [5] Baastians, M .J . (1980) . The Wigner distribution and its application to optics . [CO-Conference 'Optics in four dimensions', Ensenada, Mexico .



Volume 5, Number 3

PATTERN RECOGNITION LETTERS

X ~

f

March 198 -

*Y J11 # rY v

t

t

'

i iY R e 4 Y e d!1r s I d f

etEi*L

C I -igurc 6 .

(dl

Graimg with progressive defocus towards the left, as input for W'D ; (b) Image products at points .c= -32, x=32 from the centre and WD at those points ; (c) WD of Figure 6a along the .v-direction-

161 Hamlet . R . and H . Glunder (1983) . The Wigner distribution function of two-dimensional signals . Coherent-optical generation and display . Optfca Acta 30(12), 1789-1803 . (71 Jacobson, L . and H . Wechler (1982) . A paradigm for invariant object recognition of brightness, optical flow and binocular disparirv images . Pan- Recog. Lett . t, 61-68 . (8I Ea,ton, R .1 . ., A .J . Ticknor and H .H . Barret (1984) . Application of the Radon transform to optical production of the W iener distribution function . Opt. Eng . 23(6) . 3g--11 [91 Sru . H .H . (1982) . Two dimensional optical processing of one-dimensional acoustic data . Opt. Eng . 21(5), 804-813 . [101 Classen, T . A .C .M . and S1ecklenbrauker (1980) . The

Wigner distribution . A tool for time-frequency signal analysis, Part I, 11 and III, Philips J. Res. 35, 217-250 . 276-350 and 372-389 . [I1] lloyal, J .E . (1919) . Quantum mechanics as a statistical theory . Proc . Cambridge Philps . Soc . 15, 99-132 . [121 A .J .E .\1 . Janssen (1982) . On the locus and spread of timefrequency pseudo-density functions . Philips J. Res . 37, 79-110 . (131 Castelman, K .R . (1979) Datital Image Processing . Preniice-Hall . Englewood Cliffs, NJ . ([4] Brigham . EQ 119711 . The East Fourier Transform . Prentice-HALL Englew-ood Cliffs, NJ .

221