## Multidimensional Spatial Frequencies and Filtering

At a conceptual level, there is a great deal of similarity between 1D signal processing and signal processing in higher dimensions. For example, the intuition and knowledge gained from working with the 1D Fourier transform extends fairly straightforwardly to higher dimensions. For overviews of signal processing techniques in 2D see Lim [3], or Granlund and Knutsson for higher dimensional signal processing [4].

### 2.1 Spatial Frequency

The only difference between the classical notions of time frequencies and spatial frequencies is the function variable used: Instead of time, the variable is spatial position in the

All rights of reproduction in any form reserved.

latter case. A multidimensional sinusoidal function can be written as f = cos(xT e),

where x is the spatial position vector and e is a normalized vector defining the orientation of the wave (x, e e The signal is constant on all hyperplanes normal to e. For any ID function g, the multidimensional function f = g (*T e)

where 5 denotes the Dirac function, the Fourier transform will be concentrated on a line with orientation normal to the plane, and the function along this line will be constant.

Examples of 2D sinusoidal functions and their Fourier transforms are shown in Fig. 1. The top figures show a transform pair of a sinusoid of fairly low frequency. The Fourier transform, which contains two Dirac impulse functions, is shown to the left and may be expressed as will have a Fourier transform in which all energy is concentrated on a line through the origin with direction e. For example, for a plane in three dimensions given by

where u denotes the 2D frequency variable, m1 the spatial frequency of the signal. The bottom figures show the transform pair of a signal consisting of the same signal in F1 plus another sinusoidal signal with frequency ffl2.

FIGURE 1 (Top) Sinusoidal signal with low spatial frequency. (Bottom) Sum of the top signal and a sinusoidal with a higher spatial frequency.

FIGURE 1 (Top) Sinusoidal signal with low spatial frequency. (Bottom) Sum of the top signal and a sinusoidal with a higher spatial frequency.

### 2.2 Filtering

Linear filtering of a signal can be seen as a controlled scaling of the signal components in the frequency domain. Reducing the components in the center of the frequency domain (low frequencies), gives the high-frequency components an increased relative importance, and thus highpass filtering is performed. Filters can be made very selective. Any of the Fourier coefficients can be changed independently of the others. For example, let H be a constant function minus a pair of Dirac functions symmetrically centered in the Fourier domain with a distance | from the center,

Unsharp masking, an old technique known to photographers, is used to change the relative highpass content in an image by subtracting a blurred (lowpass filtered) version of the image [5]. This can be done optically by first developing an unsharp picture on a negative film and then using this film as a mask in a second development step. Mathematically, unsharp masking can be expressed as f = «f — ßfp

This filter, known as a notch filter, will leave all frequency components untouched, except the component that corresponds to the sinusoid in Fig. 1, which will be completely removed. A weighting of the Dirac functions will control how much of the component is removed. For example, the filter where a and 6 are positive constants, a > ¿6. When processing digital image data, it is desirable to keep the local mean of the image unchanged. If the coefficients in the lowpass filter flp are normalized, i.e., their sum equals 1, the following formulation of unsharp masking ensures unchanged local mean in the image:

By expanding the expression in the parentheses (af — 6fp) afip + a(f — fp) — 6ftp, we can write Eq. (9) as will reduce the signal component to 10% of its original value. The result of the application of this filter to the signal F1 + F2. (Fig. 1, bottom) is shown in Fig. 2. The lower-frequency component is almost invisible.

Filters for practical applications have to be more general than "remove sinusoidal component cos(raTx)." In image enhancement, filters are designed to remove noise that is spread out all over the frequency domain. It is a difficult task to design filters that remove as much noise as possible without removing important parts of the signal.

which provides a more intuitive way of writing unsharp masking. Further rewriting yields f = flp + y(f—flp) (11)

where y can be seen as a gain factor of the high frequencies. For

FIGURE 2 Notch filtering of the signal fl + f2, a sum of the sinusoids. The application of the filter h in Eq. (7) reduces the low-frequency component to one-tenth of its original value.

FIGURE 2 Notch filtering of the signal fl + f2, a sum of the sinusoids. The application of the filter h in Eq. (7) reduces the low-frequency component to one-tenth of its original value.

y = 1, the filter is the identity map, fp + fhp = fp +(f — fip) = f, and the output image equals the input image. For y > 1 the relative highpass content in the image is increased, resulting in higher contrast; for y < 1 the relative highpass content is decreased and the image gets blurred. This process can be visualized by looking at the corresponding filters involved. In the Fourier domain, the lowpass image fp can be written as the product of a lowpass filter Hp and the Fourier transform of the original image,

Figure 3 (top left) shows Hp and the highpass filter that can

Lowpass filter, Hp

be constructed thereof, 1 — Hp (top right). Figure 3 (bottom left) shows the identity filter from adding the two top filters, and (bottom right) a filter where the highpass component has been amplified by a factor of 2. Figure 4 shows a slice of a CT data set through a skull (left), and the result of unsharp masking with y = 2, i.e., a doubling of the highpass part of the image (right). The lowpass filter used (fp in Eq. (11)) was a Gaussian filter with a standard deviation of 4, implemented on a 21 x 21 grid.

A natural question arises. How should the parameters for the lowpass filter be chosen? It would be advantageous if the filter could adapt to the image automatically either through an a priori model of the data or by some estimation process. There

Highpass filter, 1 — Htp

Hip+2Hhp