Depth-resolving THz imaging with tomosynthesis

We demonstrated a depth-resolved 3D imaging technique based on absorption contrast using tomosynthesis. Tomosynthesis is similar to computed tomography except that the number of projections is much smaller. We constructed a tomosynthesis imaging system, which detects a transmitted continuous THz wave. We applied a backprojection method that was suitable for the constructed detection configuration, to reconstruct an image. Using this system, we imaged a test sample made from paper and reproduced characters written by pencil. ©2009 Optical Society of America OCIS codes: (110.6795) Terahertz imaging; (110.6955) Tomographic imaging; (110.3010) Image reconstruction techniques.


Introduction
The terahertz (THz) wave is an electromagnetic wave with a frequency that lies between radio and infrared frequencies.The frequency of the THz wave ranges from 0.1 to 10 THz, corresponding to a wavelength of 3 mm -30 µm.The recent development of THz optical devices has accelerated research into THz imaging [1][2][3][4][5][6][7][8].A variety of imaging methods have been proposed to take advantage of the characteristics of THz waves.The high transmissivity of THz waves for soft materials enables tomographic imaging analogous to X-ray computed tomography (CT), because projected images can be obtained [9].
Early THz tomographic imaging was based on measuring the time-of-flight of reflected pulses [10].This technique reproduces the 3D refractive index profiles of objects consisting of well-separated layers of different refractive indices.Although this technique provides extremely high depth resolution in the order of 1 µm, it is restricted to samples without substantial refraction or multiple reflections.Subsequently, a multi-color THz-CT system using a time-domain spectroscopy was developed [11,12].In the experiments leading to this development, tomographic information was reproduced based on the same CT reconstruction algorithm used in X-ray CT.This imaging system simultaneously forms cross-sectional images by absorption and phase contrasts in a sample with a complicated structure and a spatial resolution in the order of sub millimeters.However, unlike X-ray CT, in the range of the THz wave frequency the effects of reflection and refraction are substantial when the sample is thin and wide, and the incident angle of the samples is large.Therefore, the beam transmission sharply decreases near the edge or in regions of the sample with a complex distribution of refractive indices.Thus, it is difficult to obtain a complete set of projections, which is indispensable for making CT measurements.In conventional CT imaging, it is necessary to obtain each line integral of the complex refractive index distribution for each incident beam over 180°.Moreover, it typically takes a long time to acquire all of the data necessary using this technique.
In the current study, we demonstrate THz tomosynthesis (TS) as an alternative to conventional CT.The technique was first proposed for X-ray imaging in the 1930s [9], and is a simple but useful method for reconstruction under limited conditions [13].TS does not require measurements in all directions, and can reconstruct tomographic images using only projections within small angles.In medical applications, X-ray TS is more suitable than CT for the diagnosis of some diseases.This is because TS involves less radiation exposure and physical demand for patients, due to its requirement of a smaller selected number of projections for image reconstruction.Furthermore, TS can be used in various imaging configurations.The ability of THz-TS to reconstruct tomographic images using only projections with small incident angles is especially important when there is a need to obtain CT images of thin and wide samples.Unlike conventional CT, TS does not require "low transmission" data at great incident angles, and the artifacts due to this problem are lessened compared with conventional CT.This property of TS also provides advantages compared with other high-speed THz-CT techniques that additionally require 180° projection angles, such as wavelet-based algorithm CT [14].Similar functionality can also be achieved with Fresnel-lens THz tomography [15].With our technique, however, it is possible to perform these functions using a monochromatic THz source and ordinary optical components.Moreover, this technique does not require a multi-color or broadband source or Fresnel lens.The tomosynthesis imaging system we constructed uses a continuous wave (CW) source.The main advantages of a CW source compared with the THz-TS system include the higher signal-to-noise ratio (SNR) [16][17][18][19][20][21], compactness, and stability of the system [22].Moreover, the measurement time to obtain raw absorption-based transmission images in this system can be dramatically reduced compared with other systems.We applied the backprojection method In the TS reconstruction algorithm, which is suitable for our configuration.
In this paper, we outline the principle of tomosynthesis and our novel experimental setup.We describe the results of a TS reconstruction of the image of a sheet of paper, and discuss the implications of these results.

Principle of Tomosynthesis
TS is a depth-resolving imaging technique that emphasizes a layer of interest against the out-of-focus layers.Our method uses absorption-contrast images of an object, which are projected from various views, similar to CT [9].However, TS uses a much lower number of views at selected angles, in contrast with CT, which commonly scans rotationally 180° around an axis.TS only requires projections from several to a few tens in number.This is one or two orders less than required for CT.Although the quality of TS images is not as high as CT images due to the limited number of projections, tomographic information can be obtained under various conditions.
TS typically adopts a "shift-and-add" method to reconstruct tomographic information (Fig. 1).In Fig. 1, three locations of a light source are depicted as '1', '2' and '3′.'T', 'H' and 'Z' are denoted in the upper, middle and lower layers, respectively.As the source moves, the relative projected locations of 'T', 'H' and 'Z' are changed.The three resultant projected images can be shifted and added so as to emphasize 'T', 'H' or 'Z', and to smear out the other layers.
In the actual measurements, we rotated the sample around the y-axis to obtain projection images with different incident angles as shown in Fig. 2. The rotation angle θ is defined as the angle between the beam direction and the normal to the surface (the z-axis); the x-axis was chosen parallel to the surface.The sample was rotated from the starting position θ = −K∆θ around the y-axis, where ∆θ is the rotational step and K is a predefined integer so as that the number of projected images is 2K + 1.The sample was raster-scanned in a plane perpendicular to the beam direction.After completing a raster scan, the sample was again rotated with a step of ∆θ around the y-axis, and the raster scan was performed again.The procedure was repeated until θ = K∆θ.

Experimental setup
In X-ray TS, a projection at a single exposure is acquired using a sheet-, fan-, or cone-beam from an X-ray tube as the incident beam.In THz-TS, however, the use of a spread beam from the source provides insufficient intensity for transmission imaging.Therefore, we focused the THz beam around the center of the object and obtained each projection image by scanning the object in horizontal and vertical directions (Fig. 3).A convex lens with a focal length of 50 mm and a diameter of 30 mm was used to focus the THz beam.The transmitted wave was directed to a Schottky-barrier diode detector.In order to obtain projection images with different incident angles, we rotated the sample around the y-axis.We used a frequency-multiplier CW source with a frequency of 540 GHz from Virginia Diodes, Inc.The amplitude of the source was modulated at 30 kHz with a rectangular waveform generated by a function generator.The detected signal was fed into a preamplifier (with a gain of 5000) followed by a lock-in amplifier (time constant: 30 ms), then fed to a personal computer through a data acquisition card.The dynamic range of this system was 73 dB.All data were processed after converting the absorbance.
In TS, the spot size of the beam at the sample should be held constant because the quality of the projection images depends on the uniformity of the beam.Therefore, we evaluated the beam profile using the knife-edge method.Figures 4(a   Post-it notes (50 sheets, about 6.0 mm in thickness) were used as a sample.The letters 'T', 'H', and 'Z' with a width of 2 mm were written with an ordinary pencil on the 2nd, 25th, and 50th sheet from the top, respectively.Figure 5 is a top view of the notes.The notes were first placed on the sample stage so that the beam was perpendicularly incident to the top surface.The initial rotation angle around the y-axis was θ = −50°.At this angle, we performed horizontal and vertical scans with steps of 0.1 and 0.5 mm respectively.The horizontal scan was parallel to the sample surface and was continuous, with a speed of 0.2 mm/s, while the sample was discretely scanned in the vertical direction.Because the speed of the horizontal continuous scan was sufficiently low compared with the data acquisition time at a single data-point (30 ms), we could neglect the blurring of the projected image.The scanning area was 10 × 22 mm 2 .After completing the raster scan, the sample was again rotated at a rotational step of ∆θ = 25° around the y-axis.This procedure was repeated until θ = + 50°, and we obtained five projected images (K = 2).The overall data acquisition time was about 20 min.

Image Reconstruction: Method and Results
From the collected projection images, we reconstructed a tomographic image.Figure 6 shows a schematic procedure of the reconstruction.As a preprocessing step before the reconstruction, we introduced the deconvolution of the raw projection images by the measured beam profile as shown in Fig. 4(a).This step was undertaken because the cross-sectional spatial resolution of a TS image strongly depends on the beam diameter of the THz wave.In the deconvolution, we adopted the Wiener filter to suppress artifacts in the resultant images [23].We confirmed that the filter improved the spatial resolution in the TS images.We also found that it reduced artifacts, which were likely to have been caused by the restricted number of projections.We will discuss these artifacts in Section 5.
As in CT, we define a "dataset" from the projected images as shown in Fig. 6(a).The x 0 -axis, y 0 -axis, and z 0 -axis show the direction of the θ = 0°.One data set consists of a combination of the data from a specific row of each projected image (i.e.x-θ plane).The data set is equivalent to a so-called "sinogram".In conventional TS used the medical field, the detector and the source are shifted in one direction to obtain the projected images as shown in Fig. 1.However, in this experiment, we rotated the sample instead of shifting the source and detector, and obtained projected images that differed from typical TS images.Therefore, we did not adopt the "shift-and-add" method.Instead, we applied the backprojection method to each data set [9].In this method, the reconstructed 3D image is made by superposition from the projected image to the space corresponding to the location of the object (Fig. 6 (b)).This method of TS reconstruction is mathematically equivalent to the "shift-and-add" method.Next, a three-dimensional image was formed by piling up all of the cross-sectional images generated by backprojection (Fig. 6 (c)).Each TS image was produced by clipping an x 0 -y 0 plane from a three-dimensional image.Note that the quality of cross-sectional images in the x 0 -y 0 plane is much better than that in the x 0 -z 0 plane in TS.
The obtained raw images at the rotation angles θ = −50°, −25°, 0°, + 25°, and + 50° are shown in Fig. 7, where each image was 100 × 45 pixels in size.It can be seen that each character was shifted and narrowed gradually as the angle was increased.At θ = −50° and + 50°, the image quality was lower because of the refraction at the large incident angle.Figures 8(a)-8(c) show the images of the characters 'T', 'H' and 'Z' corresponding to locations in the 2nd, 25th, and 50th layers, respectively.Note that these images were obtained from the reconstructed three-dimensional image without the Wiener filter.The reconstructed characters 'T', 'H' and 'Z' were located in planes at depths of 0.1, 3.0 and 6.0 mm from the top surface, respectively.The characters were successfully reconstructed, although some artifacts were present.Figures 9(a)-9(c) show the images produced with the Wiener filter.It can be seen that the image of each character became sharper, and that artifacts were more suppressed than those without the filter, although a small part of the image of the character 'Z' was still missing.We will evaluate the advantage of applying the Wiener filter in the next section.

Discussion
Using THz-TS, we were able to reconstruct tomographic images using a small number of projected images within limited angles.This technique allows the acquisition of higher quality tomographic images for thin and wide samples than THz-CT, which requires 180° projections.This is because TS is able to avoid the severe decrease of the beam transmission near the incident angle of ± 90°.In this study, we adopted a backprojection method instead of a "shift-and-add" method.Although these methods are equivalent, a "shift-and-add" method is suitable for configurations in which the detector and the source are shifted in one direction.The backprojection method, however, can be applied when the sample can be rotated.In the latter method, data conversion is not required and the obtained data can be used directly to reconstruct the tomographic image.
On the other hand, TS reconstruction causes diamond-shaped artifacts in the depth dimension (z 0 -axis) due to an insufficient number of projections (Fig. 6(b)).The size of these artifacts depends on the dimension of the sample details that are being imaged.When the horizontal length of the structure perpendicular to the y-axis is large, the diamond-shaped artifacts will be emphasized.This explains why the artifacts at the position of the other characters in the different sheets were observed as shown in Figs. 8 and 9. To evaluate the effect of these artifacts, we estimated the depth spatial resolution in the z 0 -direction.
Figure 10 shows the relationship between artifacts and the horizontal length of each character.We plotted the pixel values in the z 0 -direction for the designated positions of each character (positions 1-6).We evaluated spatial resolution in the beam direction using the full width at half maximum (FWHM) of the pixel value transition for character 'H', and we used twice the half width at half maximum (HWHM) instead of the FWHM for 'T' and 'Z'.We obtained spatial resolutions of 6.8, 3.6, 3.8, 8.4, 8.0 and 4.6 mm for positions 1-6, respectively.Figure 11 shows these values plotted as a function of the horizontal length for each position.It can be seen that there was a linear correlation between the horizontal length of the character and the spatial resolution in the z 0 -direction.This correlation arose because we rotated our sample only around the y-axis.The plot also shows that the spatial resolution does not strongly depend on the depth of the character.These findings suggest that we can expect to improve spatial resolution by adding the data obtained by rotating the sample around another rotation axis (e.g.around the z-axis).In our experiment, we used a sample with three characters written on different levels such that they were not overlapping.In TS, as in conventional CT, when features of interest are overlapped, the blurring of images at different levels is determined by the spatial resolution.
On the other hand, the spatial resolution for x 0 and y 0 dimensions would not be expected to be improved substantially with the addition of data obtained by rotating the sample around another rotation axis.In reconstructed images in the current study, the typical width in these dimensions was approximately 3 mm, corresponding ~2 mm resolution across both dimensions because of the width of the written characters being 2 mm.The main reason for this blurring is that the widths of the characters in the inline configuration (e.g.θ = 50°) were elongated in the projected images.The resolution would be improved by adding the number of projected images with small incident angles rather than those with large angles.This would, however, make the z 0 resolution poorer.Thus, the measurement configuration of the projected images should be optimized by the requested spatial resolution in each direction.
Next, we quantitatively evaluated the advantage of the Wiener filter.Figure 12(a) shows an example of a clipped image in the x 0 -z 0 plane from the reconstructed TS image without the Wiener filter.Figure 12(b) shows an example with the filter included.In these images, it can be seen that applying the filter reduced the width of the backprojected beam.Thus, the features at the center are sharper in Fig. 12(b) than in Fig. 12(a).Figures 12(c) and 12(d) show cross-sectional profiles along the red lines in Figs.12(a) and 12(b), respectively.The width of the profile with the filter was 1.5 times sharper than without it.This means that the spatial resolution was improved, and that the artifacts were effectively reduced by applying the Wiener filter, because blur in the z 0 -direction leads to artifacts.This effect can be seen in the experimental results as shown in Figs. 8 and 9.
Finally, we evaluated the depth resolution for a whole character.Denote the three-dimensional reconstructed image of the kth sheet (1≤k≤50) as f k (i, j), where i and j (1≤i≤45, 1≤j≤100) are the horizontal and vertical coordinates in the reconstructed images, respectively.We obtained the standard deviation of the difference between two images (i.e., the Kth and kth images from the top) using the formula where In these equations, K is the index of the reconstructed image of interest, R K is a square area including the character of the Kth image, and N is the number of pixels.The size of R K was set to 70 × 16 pixels.We then defined a normalized blur factor of the Kth image as , , , , , max max min where "max σ K,k " and "min σ K,k " are the maximum and minimum values of σ K,k respectively, when k is changed,.Figure 13 shows B K,k plotted against k for K = 2, 25 and 50, corresponding to 'T', 'H', and 'Z'.The recognizable range of an interesting character in the z-direction was evaluated using FWHM of B K,k for K = 25.For K = 2 and 50, we evaluated the ranges using twice the HWHM of B K,k instead of the FWHM, because the layers 2 and 50 were close to the top and bottom layers, respectively, and it was impossible to evaluate the exact FWHM.The obtained values were 2.8, 2.6, and 2.0 mm for K = 2 (twice the HWHM), 25 (FWHM), and 50 (twice of HWHM), respectively.These values are consistent with the width shown in Fig. 12(d).This confirms that the depth of the character can be determined if several characters are superimposed at intervals larger than 1.0-1.4mm.

Conclusion
We applied the TS method as a depth-resolved 3D imaging technique.We constructed an experimental setup with parallel beam geometry.We adopted the backprojection method for reconstruction of the 3D image.To reduce artifacts and sharpen the reconstructed image, we used the Wiener filter for the deconvolution of the projected images by the observed beam profile before the reconstruction.Using the system, we imaged a paper test sample, and reproduced characters written in pencil.The technique can be applied to various kinds of non-destructive detection, such as the analysis of fragile ancient texts.Future research should test whether TS can be useful in reconstructing a cross-sectional image in a real 3D material with more complex layer structures.For example, it would be valuable to test a depth resolution for a sample with a complex structure in terms of refractive index and absorption.Moreover, it is important for future studies to explore the optimum combination of incident projection angles to achieve better spatial resolution.This research will be important to optimize the measurement configuration and to aid the use of TS in various practical applications.

Fig. 1 .
Fig. 1.Schematic diagrams of the principle of TS.The diagram on the left shows the geometrical relationship between the light source and the detector plane, and the diagram on the right shows how the TS images are obtained from the measured projected images.
) and 4(b) show the beam profile across and along the beam, as a function of the knife position.A beam diameter less than 1.0 mm was realized in the 10 mm depth around the focus.The sample was measured within this part of the beam.

Fig. 4 .
Fig. 4. (a) Beam profile and (b) relationship between the beam diameter and the position in the direction of beam propagation.

Fig. 5 .
Fig. 5. Photograph of the paper test sample.The sheets with 'T', 'H', and 'Z' are the 2nd, 25th, and 50th from the top, respectively.

Fig. 7 .
Fig. 7. Projection images measured with the constructed TS system, without deconvolution.

Fig. 8 .
Fig. 8. Images reconstructed using the backprojection method without the Wiener filter, at levels corresponding to (a) the 2nd sheet, (b) the 25th sheet, and (c) the 50th sheet.

Fig. 9 .
Fig. 9. Images reconstructed using the backprojection method with the Wiener filter, at levels corresponding to (a) the 2nd sheet, (b) the 25th sheet, and (c) the 50th sheet (Media 1).

Fig. 10 .Fig. 11 .
Fig. 10.The spatial resolution in the depth direction for two positions within each character.

Fig. 12 .Fig. 13 .
Fig. 12.(a) Image obtained using the backprojection without the Wiener filter, (b) image with the Wiener filter, (c) and (d) profiles along the red line in (a) and (b), respectively.