doi:10.1155/2010/575817 Research Article Point Spread Function Estimation for a Terahertz Imaging System

We present a method for estimating the point spread function of a terahertz imaging system designed to operate in reflection mode. The method is based on imaging phantoms with known geometry, which have patterns with sharp edges at all orientations. The point spread functions are obtained by a deconvolution technique in the Fourier domain. We validate our results by using the estimated point spread functions to deblur several images of natural scenes and by direct comparison with a point source response. The estimations turn out to be robust and produce consistent deblurring quality over the entire depth of the focal region of the imaging system.


Introduction
Imaging systems operating in the terahertz (THz) region of the spectrum have the potential to enable new applications due to the unique combination of properties that occur in this region, such as penetration through clothes, packaging and plastics, and also the fact that THz waves are nonionising and hence do not pose a health hazard for humans. Application domains such as security [1], medical imaging [2] and nondestructive testing [3] are likely to benefit from developments in this area. Despite these advantages, commercial systems in this spectral region have been slow to emerge, due to a lack of mature THz components and technology.
Imaging at THz frequencies poses a challenge to the resolution of the images that can be achieved [4], both because of the technology's immaturity and the long wavelengths employed (relative to wavelengths at optical frequencies), which are typically around or over the millimetre range. Due to the expensive nature of terahertz imaging systems and the likelihood of long acquisition times, there is ample scope for employing image processing techniques, without increasing the system cost nor image acquisition time. Knowledge of the point spread function (PSF) of the imaging system is very important for improving image quality.
The point spread function is the imaging system's response to an ideal, point-like source. In practical situations it may not be easy to find such ideal sources, and methods relying on direct measurement of the point spread function from the response of a point-source approximation will face the challenge of balancing resolution against sensitivity. Examples of approximations for point sources include standard stars or quasars when calibrating astronomical instruments [5,6], recording beads in microscopy [7,8], and pinholes into opaque materials for various optical systems [9]. However, most practical computational methods used for the estimation of the point spread function are indirect and rely on some measured output of the system and possibly some additional knowledge of the scene being imaged and imaging system parameters. In general, these methods are application dependent and fall into two categories: parametric and nonparametric methods. Parametric methods assume that the PSF belongs to a given shape class, modelled by a small number of parameters, such as a confusion disk or a Gaussian, and then focus on finding a robust method for estimating the parameters [10][11][12][13]. Nonparametric methods [14,15] allow for the point spread function to be of any shape although they may still impose some mild restrictions on it, such as not having a too large a support.
Here we propose an approach for calculating the PSF based on the imaging of phantom objects designed to take advantage of the imaging system characteristics. Before explaining this approach, the system design will be discussed along with properties of the PSF. The paper is organised as follows. In Section 2 we present the architecture of our experimental terahertz imaging system. In Section 3 we describe the phantoms used in our experiments and the alignment procedure. We present our PSF estimation procedure and experimental results in Section 4 and summarise our conclusions in Section 5.

System Design
In the CSIRO Wireless Technologies Laboratory, we have designed a 180 GHz coherent imaging system. The system operates in reflection mode. A reflection mode approach is required where the object being imaged does not allow THz waves to penetrate through the object whereas THz waves being scattered off the object may be detected. Practical scenarios requiring this approach include detection of skin lesions and cancers, explosive detection in packaging, and corrosion detection under paint.
The configuration of the reflection mode imaging system is shown in Figure 1(a), and the photo of the system is shown in Figure 1(b). A target to be imaged is oriented such that its surface is orthogonal to the direction of the THz beam incident on the surface. The imaging system focuses the THz beam onto the surface. Imaging proceeds by translating the target such that the target's surface remains at the focal point of the system.
The focused THz beam on the target surface is created by a quasioptical system that directs a THz beam generated from a THz source. The THz source employs a smooth-walled spline-profile horn to create a diverging Gaussian beam. This beam then strikes mirror M1 which collimates the beam. The collimated beam exhibits a Gaussian amplitude distribution and a constant phase distribution in the plane orthogonal to the direction of beam propagation. The beam's amplitude cross-section does not vary between mirrors M1 and M2. The collimated beam strikes mirror M2 and is transformed into a Gaussian beam converging towards the focal point on the target surface. A portion of the THz beam penetrates the target, and the remainder is reflected across a range of angles. The energy reflected off the target and captured by mirror M2 is coupled through the optical system back to mirror M1 and is focused towards the THz source. A silicon wafer partially reflects the energy into a THz receiver.
The amplitude of the signal captured at the receiver is used to determine the amount of energy reflected by the target. The THz source is a continuous sine wave oscillating at 180 GHz. Reflection at the target generates a reflected wave which differs from the incident sine wave in both amplitude and phase. Measurement of the phase proceeds by comparing the phase of the signal at the THz receiver with the phase of the THz source. The schematic diagram in Figure 1(c) shows the electronics that achieves this comparison. A 10 dB coupler is used to capture 10% of the transmitted signal. This signal is mixed down to an IF frequency of approximately 1.9 GHz. The signal at the THz receiver is also mixed down to 1.9 GHz. The two 1.9 GHz signals are then filtered and cross-correlated to determine the phase difference. The described system is physically large, occupying a region of approximately 1 m × 1 m. However, the system is based on electronic components which have the potential to be reconfigured in the future into a compact configuration.
The described system has a PSF that ideally should be Gaussian, have flat phase, and be invariant to the image coordinates. Invariance to image coordinates follows as the PSF does not vary as the target is translated through the fixed beam. The depth of focus and PSF size can be calculated from the properties of the source and receiver horns, which are similar to those described in [16], and the focal lengths of the mirrors M1 and M2. The resulting depth of focus is about 3.4 mm with a spot size of 2.5 mm.

Phantom Design
The core idea of our method is to evaluate the PSF by imaging "phantom" objects with precisely defined geometry. The phantom shape and alignment procedure are an extension of a design we have proposed in [17], based on a 2-value phase image, corresponding to π phase shifts. Here the phantoms were designed to produce true complex images under imaging with an ideal, delta PSF. To this end, we manufactured two phantoms of aluminium, which has very good reflective properties, representing the same geometric pattern, which consists of a series of elevated concentric disks. The radii of the disks follow a quadratic growth law and are 5.0, 5.3, 5.9, 6.8, 8.0, 9.5, 11.3, 13.4, 15.8, 18.5, and 21.5 mm, respectively. There were several reasons for choosing this particular phantom design. Firstly, we wanted a pattern that is easy to generate automatically and to be able to produce the ideal phantom images aligned with the real phantom data. We wanted it to have strong edges at all orientations, which would lead to strong components in all Fourier domain low frequencies. This would make it suitable for using a Wiener filter deconvolution technique, which we detail in Section 4. We also wanted a fair degree of variability with respect to the radial steps and, in particular, to have at least one annular ring thinner than the expected extent of the support of the point spread function. The elevation step between consecutive disks was kept constant, which means that, for an ideal, delta-like point spread function, the image of the phantom would be constant on every annular ring. The elevation step was 0.4 mm for the first phantom (which we will henceforth refer to as "Phantom 1") and 0.2 mm for the second phantom (which we will henceforth refer to as "Phantom 2"). The frequency of our system was set to 180 GHz. The elevation steps on both phantoms are not integral fractions of the wavelength, and therefore the phase variation on consecutive disks is nonperiodic. The phantoms were both placed within the focal zone of our imaging system but were not perfectly aligned, in order to test the variability of the estimated point spread function over the focal range. Pictures of the two aluminium phantoms are shown in Figure 2.

Phantom Image Registration.
To evaluate the point spread function, the procedure described in the next section requires both the actual phantom image acquired with our imaging system and a registered ideal phantom image of the phantom, under a delta point spread function. In order to generate the ideal phantom data, we firstly acquire the complex images corresponding to the two phantoms by scanning the aluminium phantoms pictured in Figure 2. From these, we extract the phase data, on which we manually identify several points lying on the circumference of the outermost circle, which corresponds to the edge of the largest disk. The coordinates of those points are fed into a procedure of least squares circle fitting, as described in [18].
and then from (1) one finds the coordinates of the circle center as (x c , y c ) = (u 1 /2, u 2 /2) and the circle radius as r = u 3 + (u 2 1 + u 2 2 )/4. Our experiments show that identifying around 15 pairs of points on the outermost circle on the phantom was sufficient to get both the circle center and radius with subpixel accuracy (i.e., feeding more points coordinates into the algorithm did not result in any significant variation.) In practice, about 20 pairs of point coordinates have been used to generate the ideal phase image data shown in Figure 3. The other inner circles are then easy to generate automatically form the known dimensions of the phantom. The phase value on the outmost flat area, outside the largest ring, is set to be equal to the dominant value of corresponding area in the acquired phase image, and the other phase values, constant on each ring, are easily computed from the elevation step and the value of the wavelength. Specifically, for every ring elevated by an additional d mm, the phase value decreases by 4πd/λ (because the distance d is travelled twice in reflection mode). In our case, at 180 GHz, the wavelength λ = 1.667 mm. The implicit assumption is that the phantom plane is perfectly perpendicular to the incident THz wave. The amplitude value on the ideal phantom data is set constant and equal to the dominant value on the area outside first ring, in the amplitude image of the measured data. Figure 3 displays the phantom phase images corresponding to the two phantoms, in the context of the imaging setup described in detail in Section 4, the ideal phantom images obtained using the data fitting procedure outlined in this section, and the difference images. Both sets of data are well aligned, with the main discrepancies occurring around the thin inner annular ring having a radius below the wavelength. However, the phase image of Phantom 2 exhibits slightly poorer alignment to its generated ideal phase data, compared to the same data from Phantom 1. This is due to the fact that, in the experimental setup, the alignment of the plane of Phantom 2 has had a tiny deviation from a 90 • angle to the incident THz beam.

Point Spread Function Estimation
If an image i(x, y) is captured with an imaging system having point spread function p(x, y) in the presence of independent additive noise n(x, y), then the resulting observed image c(x, y) satisfies the equation where • denotes convolution. Because of the commutativity of the convolution operation, the roles of i(x, y) and p(x, y) are dual to each other, which means that c and p can be used to estimate i (deblurring) or c and i can be used to estimate p (point spread function estimation). In the absence of any noise (n = 0), one could find either i in terms of c and p or p in terms of c and i from (2). The straightforward way is to take the Fourier transform on both sides, which maps convolution into multiplication, and then find the Fourier transform of the unknown (either i or p) by a simple pointwise division. However, in most practical situations, the assumption of negligible noise is unrealistic, and the direct approach suggested above would lead to strong amplification of high frequency noise. A Wiener filter approach can be used to counter the effects of noise [19]. The point spread function can be estimated from the equation where C, P, and I denote the Fourier transforms of c, p, and i and S i and S n denote the power spectra of the i and n. In most practical situations, the inverse of the signal to noise ratio is difficult to measure or estimate accurately and is EURASIP Journal on Advances in Signal Processing 5 Figure 4: Imaging scene at 180 GHz.
often approximated by a constant s, leading to the simplified version of the Wiener deconvolution: This approximation is particularly appropriate in the case of white Gaussian noise and a flat signal spectrum. In the case of reflections from a flat metal phantom, the spectrum is flat by design. Measurements of signal reflectance with no target have indeed confirmed that the noise in our system closely resembles white Gaussian noise [17]. Therefore, the simplified Wiener filter model of (4) is very suitable for our PSF estimation.

Experimental Setup.
The scene pictured in Figure 4 was imaged with our system at 180 GHz. The scene consists of the two phantoms previously described, a 50 cents dodecagonal coin, and a metal keyring in the shape of a kangaroo fixed on a Styroflex substrate. The base planes for the two phantoms are closely aligned; however, the base plane of Phantom 2 is about 0.5 mm in front of that of Phantom 1, while the coin and the keyring are about 1 to 2 mm above this plane. The entire scene fits well within the focal range of the imaging system, which is estimated to be about 3.4 mm. A tiny metal ball with diameter of about 2 mm has also been added to the scene and is located slightly above and to the right of the coin. The idea was to let this ball approximate a "point source" and use its response for a direct estimation of the PSF shape. The imaging data corresponding to the two phantoms were cropped out of the image, and the algorithm described in Section 3.1 has been applied to produce the ideal phantom data. Subsequently, the simplified version of the Wiener filter deconvolution of (4) has been applied to the measured and ideal phantom images, to produce a PSF estimation. The value of the parameter s was estimated from the measured value of our signal and the estimated power of the noise data to be in the range of 10 6 to 10 7 . The amplitude and phase of the estimated point spread functions for the two phantoms are shown in Figure 5. We remark that the phase of the point spread function is almost flat over the high intensity region of the PSF signal, which is in accordance with the phase variation of a Gaussian beam in its focal region. This is an expected result, since our source and receiver horns were designed to produce Gaussian beams.

Validation.
As a validation test, we use both estimated point spread functions, obtained on the basis of two phantom data, and the direct estimation of the point spread function from the pinball response to deblur the natural scene pictured in Figure 4, consisting of the keyring and the coin. We apply a Wiener deconvolution with s = 0.03. The amplitude data of the deblurred images, using the three estimations of the point spread function, are shown in Figure 6, and the corresponding phase images are shown in Figure 7.
The effects of the deblurring are most noticeable on the amplitude images of Figure 6. We remark the poor quality of the deblurring obtained using the PSF estimated directly from the pinball response. Admittedly, a ball with diameter around 2 mm is not a close match to an ideal point source, but it is, in practical terms, as close as we could get to it; attemps to use metal balls of smaller sizes in our experiments have resulted in response signals too weak for any reliable estimation of their shape. By contrast, both deblurred images using the PSFs estimated from phantom data show similar and remarkable detail improvement. On the coin image, the edges are sharper. Details in the area of the mouth and nose are enhanced and some texture areas of the hair become more prominent. The letters on the coin remain indistinguishable, which is to be expected, because their fine features have sizes smaller than the wavelength. However, the overall letters' blocks are still more visible on the two deblurred images. A tiny horizontal image misalignment, due to a mechanical lag of the translation table, also becomes visible in the deblurred images, especially across the upper half of the coin. The features on the kangaroo keyring are enhanced, particularly around the head and the paws areas, where the edges become better separated. The two PSFs estimated from the two phantom images are not identical, which is to be expected, given the fact that they are estimated at slightly different depths. In spite of this, the deblurring results on Figure 6(c) are remarkably similar. The patterns in the phase images, shown in Figure 7(b), are more consistent with depth geometry suggested by the optical image in Figure 7(a). By contrast, the phase image at the right of Figure 7(a) (corresponding to deblurring with the pinball PSF) is again the most inconsistent with the same depth geometry pattern.

Conclusions
We have presented a procedure for estimating the complexvalued point spread function of a terahertz imaging system which operates in reflection mode. A metal phantom with known geometry is placed at the focal region of the system and imaged. From this acquired phantom image and a computed version of the ideal phantom image, registered to the acquired image, the point spread function is estimated 6 EURASIP Journal on Advances in Signal Processing  using a Wiener filter technique. The quality of the estimated point spread function is tested by using it to deconvolve images of a scene containing manufactured objects with good reflective properties. The improvement in detail areas of the scene validates the point spread function estimation obtained using our technique. Our methodology was applied in a context where there was no practical alternative for estimating the PSF using direct measurement of a point source response. The only physical approximations of point sources for which we could obtain reasonably strong response signals were metal balls too large in size, and the PSFs estimated from such direct measurements turned out to be of poor quality. Estimations of the PSF at slightly different depths around the focal plane have produced deblurring of similar quality at all depths around the focal plane region. While in our experiments we have only tested a THz imaging system in reflection mode, the same technique could be applied for a THz imaging system operating in transmission mode; the only challenge would be the manufacturing of a phantom from materials with precisely controlled absorption coefficients. A study of such materials is presented in [20].