1.

Introduction

Digital processing of medical images and computer-aided diagnosis (CAD) have been studied extensively in order to improve the sensitivity and specificity of radiography for clinical diagnosis.^{1 2 3} However, radiographic images contain large amounts of information. A state-of-the-art film digitizer may create digital images with dimensions of 4K×5K×12 bit (40 megabytes); a prototype system for full-field digital mammography acquires digital mammograms in the 4K×6K×12 bit format.^{4 5 6 7} Digital computers process image data in series, therefore requiring either fast computers or long processing times. Without future development that would dramatically increase the speed and reduce the cost, use of the current digital processing techniques can only be applied to certain select cases in academic institutions; it may only be suitable for images with a small number of pixels or for procedures that do not require real-time processing. On the other hand, optical processing manipulates image data in parallel, thus the processing speed is not a limiting factor. Of course, lens aberrations and vignetting affect the accuracy of optical processing. Therefore, it is important to understand the advantages and disadvantages of both techniques for two-dimensional image processing in order to select the right method for medical image analysis and clinical diagnosis.

Two-dimensional (2D) image processing has been shown to be a powerful technology in many fields, such as in pattern recognition, climate reconnaissance, ocean monitoring, and medical diagnoses.^{8 9 10} In general, there are two ways to process a 2D image. The first is to process the 2D images directly in the spatial domain. Image processing in the spatial domain can be performed using a variety of methods, such as the mathematical morphological image processing, statistical image processing, differential image processing, and fuzzy logic image processing.^{11 12 13} The second is to process the images in the frequency domain. Mathematical operations are used to transform an image into a frequency distribution and processing is performed in the frequency domain. Next, the processed data is transformed into the spatial domain again through the inverse mathematical transform. The mathematical transformation and frequency domain processing can also be performed using different methods, such as Fourier transform, Hadamard transform, Walsh transform, hoteling transform, and wavelet transform image processing.¹⁴

Image processing in the spatial domain is straightforward and simpler than in the frequency domain, because no mathematical transformations between the spatial domain and frequency domain are needed.¹⁵ Spatial-domain processing usually involves operations such as scale and intensity refinements of the input images. Certain applications such as pseudocolor density encoding and phase-only filtering for 2D images cannot be processed in the spatial domain; such applications require processing in the frequency domain. Image processing in the frequency domain corresponds to manipulation of the frequency information of images.

Among the mathematical transformations used in 2D image processing, the Fourier transform is one of the most important techniques. Optical Fourier transformation can be divided into two diffraction processes during image processing, according to Abbe theory.¹⁶ First, the light falls normally upon a diffraction grating (the object), and several orders of diffraction are produced and converged to a series of points in the frequency domain. Second, since the diffraction orders are continuous, the intensities of their inverse Fourier transforms will overlap in the spatial domain. Therefore, the relationship between the object and image can in general be formulated as a double-Fourier transform. The exact image can be obtained only if the Fourier transform lenses are aberration free and if all the diffraction orders of the input object can be recovered during Fourier transform. The veracity of the 2D Fourier transform will influence the output results, especially in the case of medical imaging.

The main advantage of the optical (analog) Fourier transform is its speed. The input image is represented by a 2D intensity function and the intensities of all pixels of the input images are processed in parallel. Therefore, the analog method is suitable for real-time processing. The veracity of Fourier transform mainly depends on the aberration and the size of the Fourier transform lens. If the aberrations are corrected and the aperture-limitation effect of the lens is minimized, the veracity of 2D Fourier transform could be ideal. However, it is impossible to achieve such perfection because of the existence of aberrations and limitations of the aperture size of all lenses.

Different from the analog method, Fourier transform using a computer needs to digitize the input image first.¹⁷ The input image can also be represented by a 2D density function, but in this case each complex number will be represented by a binary series at a certain point, and then complete set of points distribution can be used as representation of the image. At the same time, some warping errors exist during image distribution processing.¹⁸ Furthermore, the aliasing effect, which is produced during sampling processing, will strongly influence the veracity of digital Fourier transforms.¹⁹ However, benefits of the digital method include high accuracy, programmability, and flexibility. Furthermore, digital processes can be easily modified in order to implement different kinds of functions and algorithms. Both methods, optical (analog) and digital (computer) Fourier transforms, can be used for 2D image processing. Each method has its own advantages and disadvantages.

In this article, the errors in both optical and digital Fourier transforms are discussed. To compare the two methods experimentally, a set of x-ray images was processed using both optical and digital methods and the results are analyzed. The selection of the processing method in real applications, based on different limiting factors, is also discussed in this article.

2.

Properties and Errors of the Optical Fourier Transform

2.1.

Optical Fourier Transform Properties of Lenses

The Fourier transform of a two-dimensional image can be implemented using a positive lens in an optical system. Such a system is shown in Figure 1(a), where the object plane P₁ serves as a light source, P₂ is the back focal plane of the Fourier transform lens, and P₃ is the principal plane of the thin Fourier transform lens. For a given complex light field f(x,y) on P₁ , its conjugate complex light field G_L(α,β) on P₂ can be determined by the Huygen’s principle²⁰

Eq. (1)

$$\begin{array}{ccc}{G}_L(\alpha ,\beta )& =& C{\displaystyle \int }{\displaystyle {\int }_{S_3}}{\displaystyle \int }{\displaystyle {\int }_{S_1}}\left[f(x,y)\text{exp}(ik{r}^{\prime })dxdy\right]\\ & & \times T(p,q)\text{exp}(ik{r}^{″})dpdq,\end{array}$$

where S₁ and S₃ denote the P₁ and P₃ surface integrals, respectively, C is an arbitrary complex constant, and k(=2π/λ) is the wave number. The distance r^′ between the point (x,y) on P₁ and the point (p,q) on P₃ and the distance r^″ between (p,q) and (α,β) on P₂ are defined as

$$\begin{array}{ccc}{\mathit{r}}^{\prime }\hfill & =\hfill & {[{\mathit{L}}^2+{(p-x)}^2+{(q-y)}^2]}^{\mathrm{1/2}}\hfill \\ \hfill & \approx \hfill & \mathit{L}+\frac{1}{2L}[{(p-x)}^2+{(q-y)}^2]\hfill \end{array}$$

Eq. (2)

$$\begin{array}{ccc}{\mathit{r}}^{″}\hfill & =\hfill & {[{\mathit{f}}^2+{(\alpha -p)}^2+{(\beta -q)}^2]}^{\mathrm{1/2}}\hfill \\ \hfill & \approx \hfill & \mathit{f}+\frac{1}{2f}[{(\alpha -p)}^2+{(\beta -q)}^2]\text{,}\hfill \end{array}$$

where L is the distance between P₁ and P₃ , and the distance between P₃ and P₂ is the focal length f of the Fourier transform lens. The phase distribution function of the thin Fourier transform lens is defined as T(p,q)=exp[−i(π/λf)(p²+q²)] .²¹ Note that the approximation in Eq. (2) is taken under the assumption that the perpendicular distance L between P₁ and P₃ and the focal length f of the lens are both much larger than the dimensions of the input and output images. Ignoring the spatial quadratic phase variation, one can obtain²²

Eq. (3)

$$\begin{array}{ccc}{\mathit{G}}_{\mathit{L}}(\alpha ,\beta )\hfill & =\hfill & {\mathit{C}}_1\text{exp}[-i\frac{\pi }{\lambda \mathit{f}}{\displaystyle \left(\frac{1-\mathit{v}}{\mathit{v}}\right)}({\alpha }^2+{\beta }^2)]\hfill \\ \hfill & \hfill & \times {\displaystyle \int }{{\displaystyle \int }}_{{\mathit{S}}_1}f(\mathit{x},\mathit{y})\text{exp}[-i\frac{\mathit{k}}{\mathit{f}}(\alpha x+\beta y)]\mathit{dxdy}\text{,}\hfill \end{array}$$

where v=f/L and C₁ is an arbitrary complex constant.

Figure 1

(a) Schematic of an optical Fourier transform. The function f(x,y) is the complex light field on object plane P₁ and G_L (α,β) is the light field on the back focal plane of lens P₂ [see Eq. (1) for the definition of G_L (α,β) ]. (p,q) is a typical point on the principal plane of thin Fourier transform lens P₃ . The lens has a focal length of f and an aperture of D. The separation between the object plane and the lens is L and the separation between the back focal plane and the lens is f. (b) Schematic of the optical point spread function. The function f(x,y) is the complex light field on object plane P₁ and g(α,β) is the light field on image plane P₂ [see Eq. (5) for the definition of g(α,β) ]. The principal plane of thin Fourier transform lens P₃ is between the two conjugate planes. The lens has a focal length of f and an aperture of D. The separation between the object plane and the lens is L and the separation between the image plane and the lens is d_i.

Equation (3) shows that G_L(α,β) is the 2D Fourier transform of f(x,y). In fact, if L=f the quadratic phase factor in Eq. (3) equals to unity. Thus Eq. (3) can be written as

Eq. (4)

$${\mathit{G}}_{\mathit{L}}(\alpha ,\beta )={\mathit{C}}_1\int {\int }_{{\mathit{S}}_1}f(\mathit{x},\mathit{y})\text{exp}[-i\frac{\mathit{k}}{\mathit{f}}(\alpha x+\beta y)]\mathit{dxdy}\text{.}$$

It must be emphasized that the exact Fourier transform takes place under the condition L=f only. In other words, the input image must be placed in the front focal plane of the lens. However, errors are inevitable due to the well-known vignetting effect²³ and to aberrations of the Fourier lens.²⁴ Such errors are analyzed next.

3.

Properties and Errors of the Digital Fourier Transform

When using a computer for digital Fourier transforms, the integration in Eq. (4) must be rewritten as a dispersed function.²⁷ For a continuous input function f(x,y) , f_s(x,y) can represent its sampling function:

We can use the comb function to represent the sampling mathematically,

We assume that F(u,v) is band limited, with a bandwidth B_x=u _max on the x axis and B_y=v _max on the y axis.²⁸ One must choose the sampling distances according to the following rule:²⁹

The specific sampling distances D_x and D_y in Eq. (15), as the reciprocal of the spatial bandwidth, are called the Nyquist rate. The Nyquist rate is the largest sampling distance that can be chosen without introducing the aliasing effect.²⁹

Following Eqs. (11 12 13 14 15), the inverse Fourier transform of the input function can be rewritten as

The errors in the digital Fourier transform spectra are mainly caused by the discrete processing of data. If the sampling distances are larger than the reciprocal of the spatial bandwidth, the aliasing effect will appear.

4.

Comparison Between Analogy and Digital Fourier Spectra

4.3.

Comparison Between Optical and Digital Fourier Spectra

A set of medical images was processed using both optical and digital Fourier transform methods. The original medical images on x-ray films were optically reduced and recorded on 36 mm×24 mm photographic films, as shown in Figures 2(a), 3(a), and 4(a). Having completed the processing procedures described in Secs. 4.1 and 4.2, the optical Fourier spectra of the medical images are given in Figures 2(b), 3(b), and 4(b), and the digital Fourier spectra are given in Figures 2(c), 3(c), and 4(c).

Figure 2

(a) Image of the letter E with a black background used for (b) optical and (c) digital Fourier transforms. Because of the rectangular pattern of the letter E, the Fourier spectrum is distributed along the x and y axes only. Using the optical method, the spectrum contains more distinctive diffraction orders, shown in (b).

Figure 3

(a) X-ray medical image of breast tissue. The original image was optically reduced and recorded on a film with a dimension of 36 mm×24 mm and with resolution of 200 lines/mm. (b) Optical Fourier spectrum. (c) Digital Fourier spectrum. In the center region, both methods provide the same spectrum, as evidenced by the bright spots at the centers of both (b) and (c). In the high frequency region (the periphery of the image), the optical method clearly provides a wider distribution of the spectrum.

Figure 4

(a) X-ray medical image of bones. The original image was optically reduced and recorded on a film of 36 mm×24 mm and with resolution of 200 lines/mm. (b) Optical Fourier spectrum. (c) Digital Fourier spectrum. Because of the image pattern of the bone, the Fourier spectrum is distributed along the bisection of the first and the third quadrants. The two images show a significant difference in the peripheral regions though they are very similar in the center region. The optical method provides more information of higher frequency.

Figure 2(a) is an image of letter the E with a black background. As expected, due to the rectangular shape of the image, the spectra in both cases showed a distribution along the x axis and the y axis. The optical Fourier spectrum in Figure 2(b) shows more details of the structure than the digital spectrum in Figure 2(c) under the same display conditions.

The optical and digital Fourier spectra of an x-ray image of breast tissue complex in structure [Figure 3(a)] are given in Figures 3(b) and 3(c). Apparently, the center region (low frequency region) of both spectra has the same distribution pattern, evidenced by the close resemblance of the center bright spots. However, differences in the periphery of the spectra are seen clearly. The optical method was able to process signals of higher frequency, as shown by the wider spectrum distribution in Figure 3(b). Optical and digital Fourier spectra of an x-ray image of bones [Figure 4(a)] are shown in Figures 4(b) and 4(c). These bone images were used to show the result when the image with more sharp contrast was inputted. Both Figures 4(b) and 4(c) show similar patterns in the center region with the spectrum distribution mainly along the bisection of the first and third quadrants. Figure 4(b), however, shows more details of the structure in the high frequency region (peripheral region).

Figures 2 3 4 showed the following common features of the two transformation methods: (1) the results in the center regions (low-frequency region) of the spectra are nearly identical; (2) the optical method provides more details in the high frequency regions.

The shortcomings of the digital Fourier transform in this experiment are largely due to the fact that the number of sampling points in the processing is limited. If enough pixels were sampled to satisfy the same frequency cutoff like in optical method, the veracity of both transformations would be similar. Under the specific conditions of our experiments, the digital processing points required should be 9000×4000 based on the required sampling distances given by Eq. (15). However, the digital Fourier transform in our experiments required about 30 s to process an image with 3000×1333 points and it will require an estimated 500 s to process one of 9000×4000 points using a personal computer. The spectra in Figures 2(c), 3(c), and 4(c) obtained using fewer sampling points—one ninth of the points required for the same optical cutoff frequency—clearly demonstrated the aliasing effect. If enough pixels were sampled, the digital method can provide highly accurate results, provided that the processing time is not a limiting factor.

The disadvantages of optical method mainly lie in the complex structure of the optical system, which requires an assembly of various optical components. Thus, the level of difficulty of data processing increases. Furthermore, all the analyses here are based on the assumption of a thin lens. In reality, the thickness of the lens will always introduce noticeable errors. Based on our analysis of errors introduced by lens aberration [see Eq. (10)], optical processing can provide the correct power of the spectrum, but the accuracy of the phase distribution is limited.

On the positive side, the optical method has the advantages of real-time parallel processing for a Fourier transform of larger size images. The computer can provide greater flexibility and high accuracy.

5.

Discussions

The Fourier transform is one of the most important methods in 2D image processing. It has several properties that make it attractive for image processing applications. In this article, the characteristics of analog (optical processing) and discrete (digital processing) Fourier transforms are analyzed.

Specifically, we discussed the effect of lens aberration on the optical Fourier transform by introducing the general phase distribution function. As a consequence, the lens aberrations only contribute a phase correction factor, and do not affect the overall intensity and the cutoff frequency of the Fourier spectrum. The vignetting effect contributes to errors of the optical spectrum by introducing a cutoff frequency that is dependent on the size of the lens.

The main factors that contribute to the errors in the digital Fourier transform are the aliasing effect and the band limit. The former is due to sampling of a small number of pixels and it affects the veracity of the spectrum in the high frequency region, as shown in our experiments (see the spectra in Figures 2 3 4). The latter is due to the requirement of digital processing; the band limit automatically restricts the high-frequency components in the Fourier spectrum, as determined by Eq. (15). In principle, digital processing can provide a high degree of veracity if there is a large enough number of sampling points. However, a long processing time is required in such cases. Therefore, the difficult task is to keep a balance between the high accuracy desired (more pixels) and the processing speed desired (fewer pixels).

Overall, the error analysis and our experimental results indicate that the outcomes of the two Fourier transform methods for medical images are similar in the center region of the Fourier plane (the lower-frequency region), but different in the peripheral region (the higher-frequency region). In general, the analog method (optical processing) is suitable for medical images with a large number of pixels and for clinical procedures that require real-time operations, such as fluoroscopy because of its required high processing speed. Digital processing, on the other hand, can provide greater flexibility and programmability.

In summary, optical processors have good analog quality and high speed through parallel processing. In order to achieve high accuracy, optical components with minimal aberration are needed. In addition, current optical processors are less flexible since they are more limited in programmability than digital processors. One technical solution with which to improve their flexibility and programmability is to develop hybrid devices, namely, parallel optical processors equipped with optimal digital–optical interfaces.