Open Access
Add to BibSonomy
Abstract
A 2D tomographic terahertz imaging set-up using a single pixel imaging approach is realized, where a liquid helium cooled bolometer is utilized as a bucket detector and a mercury-arc lamp serves as a broadband terahertz source. The different patterns of the terahertz radiation, which are necessary for the single pixel imaging approach, are realized by spatially addressed photodoping of a high resistivity float zone silicon window, employing a near-infrared laser diode, which is spatially modulated by a digital micromirror device. The two investigated sample objects have cylindrical and cuboid shapes and consist of polypropylene. Both sample shapes cause strong influences of refraction, reflection and diffraction, which distort the measured projections and thus have to be considered in the tomographic reconstruction. In order to consider these effects, a model is developed which combines refraction and diffraction effects by a hybrid approach using ray tracing and scalar diffraction theory yielding finally projections of the sample objects. These simulated projections are compared to the measured projections and show a good agreement between the experimental results and the developed model. In accordance with this result, an optimization problem is formulated, which offers an approach for tomographic reconstruction using the developed model.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Terahertz imaging techniques are evolving towards excellent tools for a plurality of industrial as well as societal applications. For instance, terahertz imaging techniques can be applied to security tasks [1], medical diagnostics [2, 3] or quality control [4]. In these fields of application, the spectrum of available implementations include optical [5], electro-optical [2] and electronic [6] schemes for the generation and detection of the terahertz radiation. Also the exploited imaging techniques range from mechanically raster scanning systems [6] to multi-pixel detector arrays with a large number of pixels [7]. Each of these imaging techniques provides its individual advantages and disadvantages, as for example multi-pixel detector arrays commonly increase the cost of the imaging system, but offer high frame rates of up to 450 frames per second [8]. Another promising imaging method is the so-called ’single pixel camera’, which uses a non-spatially resolving ’bucket’ detector and does not require any mechanically raster scanning of the object. This method reconstructs an image of the object by a set of measurements of the total intensity transmitted, reflected or scattered from the object under illumination with differently spatial patterned radiation. Real time videos, simultaneously taken in the visible and short-wave infrared have been realized [9] and also 3D imaging [10] has been demonstrated using a ’single pixel camera’ in the visible spectral range. Both random [11] and Hadamard [12] based patterns have been successfully used to perform imaging in the terahertz frequency range. Approaches for the generation of the patterns in the terahertz frequency domain range from masks printed on printed-circuit boards [11], all optical techniques [12] and metamaterial based spatial light modulators [13]. A promising progress of the single-pixel imaging method has been achieved in the terahertz frequency domain, leading for example to sub wavelength spatial resolution [14].
In order to realize an insight into a hidden object, rather than imaging a projection, tomographic reconstruction techniques can be applied. This is useful in security applications and a plurality of industrial processes [15]. Since the photon energy of terahertz radiation is non-ionizing, there are significant lower safety requirements compared to X-rays, easing the implementation of terahertz systems.
However, tomographic imaging in this spectral region entails characteristic challenges. Optical effects like refraction, reflection and diffraction influence the measured projections and can lead to distortions in the reconstructed image. For this reason the effects of refraction on the measured projections is discussed using ray tracing simulations in several articles [16,17]. Therefore, special image reconstruction algorithms, which for example compensate the influence of refraction for a cylindrical object, have been developed [18]. In cases where amplitude and phase of the terahertz radiation are available, image reconstruction methods known from the field of seismic imaging can be applied [19].
In this publication a ’single pixel camera’ 2D tomographic terahertz imaging concept is developed and used to reconstruct 2D images of two different sample geometries, demonstrating that raster scanning of the object and multi-pixel detector arrays are not mandatory for tomographic terahertz imaging. The image reconstruction is performed using a standard X-ray algorithm, which does not include the discussed optical effects and hence the reconstructed images show some characteristic distortions. In order to improve the image quality and to incorporate optical effects a model based on ray tracing and scalar diffraction theory is developed. With this model an optimization problem is formulated and thereby an object recognition is realized.
2. Experimental set-up
The experiment is integrated in a Bruker Vertex 80v Fourier transform spectrometer as shown in the photograph in Fig. 1(a), which is in line with the schematic drawing of the experimental set-up depicted in Fig. 1(b). A built-in mercury-arc lamp serves as a broadband terahertz source emitting radiation between 0.45 and 3 THz. This frequency range is caused by the combination of a beamsplitter installed in the spectrometer, a far-infrared cut-off filter incorporated into the detecting bolometer (Infrared Laboratories, 4.2 K Silicon Bolometer) and the outer vacuum window of the bolometer, which is a wedged white polyethylene window. The filter is a 1 mm thick diamond wedged crystal quartz coated on one face with garnet powder and ensures that only radiation below 3 THz is measured by the bolometer. The terahertz radiation is collimated and passes the Michelson interferometer, whose mirrors are fixed during the measurements, but offer the opportunity to realize spectral imaging in the future. Behind the Michelson interferometer the radiation is focused into the sample compartment as shown in the schematic of the experimental set-up in Fig. 1(b). Here it is transmitted through a 1 mm thick high resistivity float zone silicon (HRFZ-Si) window (Tydex W-HRFZ-SI-D25.4-T1) and the object under investigation. Finally a parabolic mirror focuses the radiation onto a liquid helium cooled bolometer, which measures the total optical intensity collected by this mirror. The HRFZ-Si window is irradiated by a 1 W 808 nm laser diode, whose collimated beam is spatially modulated by a digital micromirror device (DMD). To realize a lock-in detection of the bolometer signal, additional to the spatial modulation, the laser diode is temporally modulated by a chopper wheel. Because of the relatively low frequency bandwidth of the bolometer, the modulation frequency amounts to 200 Hz. The so-called photodoping of the HRFZ-Si window leads to a change of its transmission characteristics in the terahertz spectral region [20], transferring the temporal and spatial modulation of the near-infrared radiation into the terahertz radiation, yielding the spatial resolution of the imaging scheme. By the photodoping technique it is possible to achieve a modulation of the terahertz radiation in a broad spectral range from below 0.1 to over 3 THz [21]. Finally, the signal of the bolometer is fed to a lock-in amplifier, which extracts the signal amplitude and is readout by a computer.
Fig. 1 (a) Photograph of the experimental set-up. (b) Schematic depiction of the 2D tomographic imaging set-up. The inset shows the 32nd pattern as it is radiated on the HRFZ-Si window by the 808 nm laser diode.
The used DMD (Texas Instruments, DLP300) possesses 608 mirrors in x- and 684 mirrors in y-direction, which are combined to 32 stripes consisting of 19 mirrors in x-direction and 684 mirrors in y-direction. Thus a spatial resolution of 32 × 1 pixels is realized in the terahertz frequency range. The maximum number of used pixels is determined by the measurement time (more pixels correspond to longer measurement times) and is limited by a minimum reasonable width of each pixel at the HRFZ-Si window, which on the other hand is determined by the investigated wavelength of the radiation. In the realized set-up, the area concurrently covered by the terahertz radiation and the near-infrared radiation on the HRFZ-Si window amounts to approximately 19 mm. This results in a pixel width of 0.6 mm, which is in the order of the wavelength of the used radiation. In order to gain spatial informations about the transmission properties of the sample object, the HRFZ-Si window and thus the terahertz radiation is spatially modulated by a multitude of patterns, consisting of 32 pixels in x-direction, defined by a 32 × 32 Hadamard matrix H. Since the Hadamard matrix only consists of the values +1 and −1, each row of H represents two patterns. As an example for the pattern generation, the instruction for the creation of the first two and last two patterns is shown in Fig. 2. The H+ patterns are created by transparent regions represented by +1 entries and opaque regions for the −1 entries, while this is reversed for the H− patterns. As an example, the resulting 32nd H− pattern is shown as inset in Fig. 1(b), as it is radiated on the HRFZ-Si window by the 808 nm laser diode. For each of the 64 pattern (32 × H+ and 32 × H−) the total transmitted terahertz radiation is measured by the bolometer, resulting in two measurement vectors M⃗+ and M⃗−, which are related to H by
where n is the order of the Hadamard matrix, 𝟙nn is the matrix of ones (aij = 1) and O⃗ represents the spatially resolved transmission of the sample object. The 64 measured intensities are then used to calculate the spatially resolved transmission of the object. The reconstruction procedure gets apparent by rearranging Eq. (1) Hence, the spatially resolved transmission of the sample object can directly be calculated from the two measurement vectors and the Hadamard matrix. The measured intensities in M⃗+ and M⃗− are always a sum of intensities arising from multiple transparent regions defined by the corresponding pattern and are therefore higher as in case of an intensity measurement for each pixel on its own. These higher intensities lead to a better signal to noise ratio at the detector, which is one of the advantages of Hadamard imaging [14]. In addition, the origin of a ray does not get lost due to refraction since the patterns are generated in front of the sample object, which does not hold true for imaging using a multi-pixel detector array.Fig. 2 Illustration of the pattern generation instructions of a 32 × 32 Hadamard matrix used for the single pixel imaging. Each row of the Hadamard matrix represents two patterns. The H+ patterns are created by transparent regions represented by +1 entries and opaque regions for the −1 entries, while this is reversed for the H− patterns.
Before every tomographic measurement an image of the terahertz beam profile is taken using the formalism described above, while no object placed in the sample compartment. This beam profile is used to normalize the transmission through the sample object, which is realized by dividing the object image by the measured beam profile as shown in Fig. 3. In this way a normalized object transmission is received, which is then converted into absorption values and used for the tomographic image reconstruction. Figure 3 shows the beam profile of the utilized mercury-arc lamp at the sample compartment of the spectrometer. On the left (pixel position 1 – 6) and right side (pixel position 22–29) of the beam profile the intensity rises, reaching an intensity maximum of 12.5 and 13 (in arb. u.), respectively. In between these two maxima a drop of intensity is observed, which is caused by a window for an adjustment laser in the center of the beam splitter installed inside the spectrometer. This drop of the intensity is removed from the object transmission by the normalization procedure as can be seen in Fig. 3.
Fig. 3 Typical example of the measured beam profile and the raw object transmission for a cuboid sample object under a rotation angle of 120°. The object is divided by the beam profile to achieve the normalized transmission of the object as shown on the right. The normalized transmissions are then converted into absorption values, which are used for the tomographic image reconstruction. Due to the low intensity on the left and right side of the beam profile the outer pixels are discarded.
3. Experimental results
The concept discussed in the previous section is now applied to two dimensional tomographic imaging, for the investigation of cylindrical and cuboid shapes of sample objects consisting of polypropylene (PP). Since tomographic imaging requires spatially resolved transmission measurements of the sample objects under different rotational angles θy, the sample objects are placed on a rotational table and x-direction resolved projections of the sample object are measured using single pixel imaging at each rotation angle from 0° to 180° in 10° iterations. A projection (Fig. 4(a)) and the dependency of the projections on the rotation angle of the sample object (Fig. 4(b)), the so-called sinogram, of a PP cylinder with a diameter of 14 mm is shown in Fig. 4. In order to avoid aliasing in the rasterisation of the projections, occurring in the reconstruction algorithm for the 2D reconstruction, the resolution of the received transmission projections for every rotation angle θy is increased by a factor of 10 using a linear interpolation in between the measured points.
Fig. 4 (a) Measured projection under a rotation angle of 0° and (b) measured sinogram of a cylindrical PP sample with a diameter of 14 mm. From the projection (a), the diameter of the cylinder can be estimated between 13 and 14.5 mm.
At first, the projection shows a low absorption through the atmosphere in the periphery at the positions 0 to 3 and 27 to 28. The experiments are performed under ambient conditions and the contrast between the sample and the surrounding atmosphere is lower due to atmospheric absorption compared to the contrast between sample and vacuum. In the area in between, the absorption rises to a maximum of approximately 1 forming plateaus of no terahertz transmission. In the center (position 15) of the projection a reduced absorption peak can be observed. In the terahertz spectral region, this kind of projection shape is typical for a cylindric object, and has been observed several times [16–18]. Here, the plateaus do not represent regions of the sample object with a high absorption coefficient, they are rather the result of combined refraction and reflection effects. These optical effects have a particularly large impact on the measured projections in the terahertz frequency spectrum, leading to a lens like focusing of the radiation after a circular shape. A rotation of the cylindrical sample objects shows no change in the shape of the measured projections, as can be seen from the sinogram in Fig. 4(b). Only a slight shift of the central transmission peak can be observed, indicating a small difference of the rotational table axis and the cylinder center.
While the measured projections of the cylindrical sample object show strong influences of refraction and reflection effects caused by the bended surface, a cuboid sample object with one of its flat faces orthogonal to the beam propagation direction should only show some additional losses due to reflection, since the rays incident parallel to the surface normal. Two different projections measured at an rotation angle of 180° (Fig. 5(a)) and 120° (Fig. 5(b)) and the measured sinogram (Fig. 5(d)) of a cuboid PP sample with side lengths of 14 × 7 mm are shown in Fig. 5. It is noticeable that in the case of a 180° rotational angle of the sample (Fig. 5(a)), which corresponds to a geometric arrangement, where the rays incident parallel to the surface normal, a multitude of features can be observed. These features do not agree with the previous assumption that only some additional reflection losses occur, but rather can be attributed to diffraction effects. The outer dimensions of the cuboid sample object can be recognized due to the two absorption peaks located at position 4 and 25, which indicate the side edges. In contrast to the high absorption plateaus of the cylindrical sample object, there is still some transmission of terahertz radiation at these two absorption peaks. Next to each of the two peaks, there are two regions of slowly decreasing absorption (positions 7–11 and 17–22) representing the inner of the sample. In the center of the measured projection a small dip of low absorption is recognized. The described shape of the projection changes quite a lot, if the sample object is rotated as can be seen in Fig. 5(b). Here, the first positions (1–3) show a low absorption through the atmosphere, followed by a broad region (positions 7 – 11) of no transmission. Subsequently, a region with a constant lower absorption coefficient follows, which is slightly higher than the absorption of the atmosphere. The last characteristic is a broad peak at the positions 20 – 26, in whose maximum no terahertz radiation is transmitted. These three differently sized regions of transmittance can be understood as various refraction zones of the tilted cuboid, as can be seen in Fig. 5(c). Figure 5(c) shows the results for a ray-tracing simulation of a 14 × 7 mm large rectangle with an refractive index of 1.5 and under a rotation of 120°. First, there is radiation which incidents on the long front side, but is refracted to the short side (I). Secondly, there is a region, where the radiation incidents on the long side and leaves the object on the long rear side (II) and lastly the radiation incidents on the short side of the sample (III). Total reflection (red rays) has a large impact on the absorption coefficients in the first and third case, since radiation inside the object hits the surface under large angles. The size of each of the three regions vary with the rotation angle of the object, which can be seen in the measured sinogram depicted in Fig. 5(d).
Fig. 5 Two measured projections under a rotation angle of (a) 180° and (b) 120° and (d) measured sinogram of a cuboid PP sample with edge lengths of 14 × 7 mm. (c) Ray tracing simulation for a rectangular object with edge lengths of 14 × 7 mm and a refractive index of 1.5 under a rotation of 120°. Red colored rays indicate rays which get totally reflected on the rear side, while green rays can pass the object.
In order to reconstruct a 2D image of the two investigated PP samples the ASTRA toolbox [22] with its implementation of the simultaneous algebraic reconstruction technique (SART) is used. This iterative algorithm is chosen because it achieves a better image quality with a small number of measured projections compared to the inverse radon transformation [6]. The reconstructed images are shown in Fig. 6(a) for the cylindrical and cuboid (Fig. 6(b)) PP sample. The red circular shape of the cylindrical PP sample (Fig. 6(a)) with absorption values between 0.75 and 1 is clearly visible, whereas the surrounding atmosphere is represented by green and yellow coloring, representing absorption values between 0.5 and 0.65. As already indicated in the sinogram (Fig. 4(b)), the cylinder is not centered and possesses an additional low absorbing region in its center. This fully transparent blue area in the center is an artifact caused by the lens like behavior of the PP cylinder and leads to a distortion of the homogeneous sample object as discussed before.
Fig. 6 Reconstructed images of (a) a cylindrical and (b) a cuboid PP sample using the measured sinograms and performing the reconstruction by applying the SART algorithm.
The reconstructed 2D image of the cuboid PP sample shows a lot of distortions and its shape is difficult to recognize. The shorter 7 mm sides of the cuboid can be seen as red areas with an absorption between 0.75 and 1 in Fig. 6(b), while the longer 14 mm sides are yellow areas with absorption values around 0.65, which can be hardly differentiated via the color from the surrounding atmosphere. The horizontal and vertical axes of the rectangle are indicated by the dotted white lines in Fig. 6(b). By taking into account a pixel width of 0.6 mm the dimensions of the reconstructed rectangle amounts to 15.1 × 6.3 mm and is in the range of dimension measured using a caliper (14 × 7 mm). Actually, the object consists of a homogeneous PP cuboid and therefore the reconstructed image should show a rectangle of constant absorption. However, the interior is represented as a large dark blue area enclosed by the red and yellow surfaces of the sample object. The main reason for the distortions of the reconstructed images originate in the used reconstruction algorithm, which actually has been developed for X-ray tomography, i.e. wavelengths between 5 and 60 pm, thus much smaller than the object dimensions. Moreover, the refractive index of matter differs very little from unity in the X-ray spectral region, hence the algorithm assumes straight rays, which penetrate the sample object and get attenuated depending on the spatially distributed absorption coefficient and the optical path length through the objects. Although, the algorithm does not directly depend on the wavelength of the used radiation, it ignores optical effects like refraction, reflection and diffraction. But this circumstance is not given in case of tomographic terahertz imaging, since refraction of the rays, losses due to reflection and diffraction effects occur and strongly influence the measured projections. For this purpose, both new concepts for the field of tomographic terahertz imaging have to be developed, and improvements of the experimental set-ups, like for example refractive index matching [23], or special reconstruction algorithms [18] are needed. Also a good understanding of the impact of refraction, reflection and diffraction effects on terahertz tomography is required in the latter case, in order to include these effects to the tomographic reconstruction.
4. Refraction, reflection and diffraction modeling for tomographic terahertz imaging
The measured projections from the previous section clearly demonstrate the impact of refraction, reflection and diffraction effects. The influence of refraction can be well observed in case of the cylindrical object in Fig. 4(a). The effect of reflection can be seen in Fig. 5(b) and Fig. 5(a) shows strong diffraction effects. In oder to get a better understanding and to verify these observations a model is developed, which takes all these effects into account. The model is designed to take the same two sample geometries (cuboid and circular), as in the case of the experimental results and calculate the resulting sinogram for the selected sample object using the single pixel imaging technique taking into account the previously mentioned optical effects. Therefore, the model will be capable to reproduce the measured sinograms, which can also serve as a benchmark for the developed model.
The basic concept of the model is taken from [24] and it is adapted to the single-pixel tomographic imaging set-up. The model is based on a hybrid ansatz combining ray tracing for reflection and refraction effects with diffractive light propagation using the Rayleigh-Sommerfeld diffraction solution. A sketch for the case of a cylindrical sample object is depicted in Fig. 7. First, the initial constant electric field E0 passes the Hadamard mask Hi. Therefore, the electric field after each of the 64 binary Hadamard transmission masks is calculated using the Rayleigh-Sommerfeld diffraction solution [25]:
Here d is the distance between the Hadamard pattern and the object plane, λ is the wavelength, k is the wavenumber and r12 is the distance between a position on the Hadamard mask and the object plane. This electric field Eobj(x) serves as source field in the object plane at z = 0 and is divided into one part Eobj,i(x) incident on the object and another part Eobj,ni(x) not interacting with the object. The strength of the electric field Eobj,i(x) on the object is corrected by the Fresnel formulas and the Lambert-Beer law for each x1 coordinate, taking into account x1 dependent incident angles and optical path lengths through the object. Furthermore, the phase of this electric field incident on the object is corrected to match the phase at the position (z2, x2) (which is the position on which the ray leaves the object), while the ray seems to come from the position (0, x′). The coordinates (0, x′) and (z2, x2) depend on the original x-displacement x1, the geometry of the investigated object and its rotation angle. In order to obtain these coordinates optical ray tracing is performed. Each ray gets refracted on the ambient-to-object and object-to-ambient interface before reaching the detector plane. For an observer remaining in this plane, the ray seems to originate from the coordinate (0, x′). Hence, the x-coordinates of the electric field Eobj,i(x) at the position (0, x1) are shifted to the coordinates (0, x′) for the calculation of the electric field Edet(x) at the detector plane. The total electric field at the detector surface is again calculated using the Rayleigh-Sommerfeld diffraction solution using the unchanged electric field Eobj,ni(x) plus the corrected and shifted electric field Eobj,i(x). Finally, the intensity at the detector surface is summed up for each Hadamard mask and Eq. (2) is used to calculate the projection.Fig. 7 Schematic sketch of the model combining Hadamard imaging and refraction, reflection and diffraction effects for projection calculations of a circular PP object.
The model is used to calculate the projections of circular and rectangular objects under arbitrary rotation angles. As in the measurements, the objects consist of PP, which possess a refractive index of approximately 1.5 and an absorption coefficient in the order of 1 cm−1 at the terahertz frequency under investigation [26]. The object under investigation is rotated from 0° to 180° in 1° iterations and for each rotation angle the projection is calculated at a frequency of 1 THz. As in the case of the measured projections, which are normalized by the beam profile shown in Fig. 3, first of all a calculation of the intensity distribution is performed without an object, which is used to normalize the calculated projections including an object. Thus, absorption values of 1 correspond to fully opaque regions for the terahertz radiation, whereas an absorption of 0 implies a fully transparent region. Since the cylindrical object possesses a rotational symmetry only one projection is calculated and used for all rotation angles. The calculated projection for a PP cylinder with a diameter of 14 mm is shown in Fig. 8(a) together with the corresponding sinogram (Fig. 8(b)). Both sides of the projection show no absorption through the atmosphere, which changes abruptly at the object borders. Here, the object appears as fully opaque plateaus, which are only interrupted in the center of the projection by a valley of reduced absorption. The calculated projection resembles the measured one shown in Fig. 4(a), only the width and depth of the central peak differ slightly. This may be caused by the broad frequency range of the used terahertz radiation, which influences the measured sinogram, since parameters like the refraction index and absorption coefficient depend on the wavelength.
Fig. 8 (a) Calculated projection under an arbitrary rotation angle and (b) sinogram of a cylindrical PP sample with diameter of 14 mm at a frequency of 1 THz.
The model is now applied to a rectangular shaped PP sample object, for which two calculated projections at a rotation angle of 180° (Fig. 9(a)) and 120° (Fig. 9(b)) as well as the sinogram (Fig. 9(c)) are shown in Fig. 9. The effects of diffraction are especially apparent in the case of the 180° projection, because no refraction effects occur, which is caused by the plane incidence of the radiation on the surface of the rectangular samples object. Nevertheless, the projection exhibits some characteristics which do not depend on the absorption of the sample object. This demonstrates the importance to incorporate diffraction effects to the modeling of the tomographic imaging set-up in order to accurately reproduce the measured projections. The same absorption peaks between positions 9 and 25 are visible like it has been shown in case of the measured projection shown in Fig. 5(a). As before, these two peaks occur at the transition from the surrounding atmosphere to the sample object. Furthermore, the interior of the object’s projection at 180° also resembles the measured one. In between the two absorption peaks, the absorption decreases slowly and leads to an absorption minimum in the center of the object. By rotating the object the contributions of refraction and reflection effects to the projections increase and are the dominant source of distortions, as shown in the sinogram by yellow regions of high absorption (Fig. 9(c)). The yellow x-shaped structure in the sinogram can be attributed to the three regions as described in section 3 and depicted in Fig. 5(c). The radiation incident on the object’s top and bottom sides (region I and III) leads to non-transmitting zones of the object, while the sizes of these zones are dependent on the rotational angle. In Fig. 9(b) the projection of the rectangular PP object at a rotational angle of 120° clearly shows these three different transmission zones.
Fig. 9 Calculated projections under a rotation angle of (a) 180° and (b) 120° and (c) sinogram of a rectangular PP sample with edge lengths of 14 × 7 mm at a frequency of 1 THz.
In order to reconstruct the 2D images of both the circular and rectangular PP objects, again the SART algorithm from the ASTRA toolbox [22] is used. The resulting images are shown in Fig. 10(a) for the circular and (Fig. 10(b)) rectangular samples at a frequency of 1 THz. Here the circular PP object is shown as a hollow circle, a fact which can be attributed to the lensing effect of the object, but clearly stands out from the surrounding atmosphere. The good contrast between the atmosphere and object reduces in case of the rectangular shaped sample as can be seen in Fig. 10(b). The left and right sides of the rectangular are each represented by a thick red vertical line and are clearly emphasized in relation to the surrounding atmosphere, whereas the front and rear sides are only visible as thin yellow horizontal lines. Finally, the interior of the rectangle does not even contrasts from the atmosphere, which makes the object appear to be hollow as also observed in the measured 2D image of the rectangular sample depicted in Fig. 6(b). In conclusion, the results of the developed modeling approach for 2D Hadamard terahertz tomography are in good agreement with the observations obtained from the experimentally realized tomographic imaging. Calculated projections of both circular and rectangular objects under different rotation angles reproduce the findings of the measurements and clarify the effects of refraction, reflection and diffraction in the field of terahertz imaging.
Fig. 10 Reconstructed images of modeled (a) cylindrical and (b) rectangular PP samples using the calculated sinograms for a frequency of 1 THz and performing the reconstruction by applying the SART algorithm.
5. Towards tomographic reconstruction using the developed model
After the development of a model, which successfully reproduces the measurements, this model will be now used to reconstruct the shape of the investigated sample objects. For this purpose a genetic algorithm, implemented by the Matlab® function “ga”, is used to minimize the mean squared error (MSE), which is defined by
Here n is the number of measurement points, θy and x are the rotational angle and the x-position of the object respectively, Data(θy, x) represents a measurement point of the sinogram and Model(θy, x, a, b, s) is the function of the developed model. The parameters a, b and s are the fitting parameters of the model, by which the genetic algorithm tries to minimize the MSE. For this purpose, the genetic algorithm creates a sequence of populations after an initial random population. A population contains many individuals and each individual represents a set of parameters for the optimization problem. From the individuals of the current population a new population is generated by three different mechanisms. First, elite children are selected, which are given by individuals in the current generation with the best MSE values. Second, crossover children are generated by combining two parents and third, children are created by mutations like random changes. By these mechanisms, the individuals approach the global solution of the optimization problem with increasing population and do not get stuck in local minima.The parameter s has an integer value of 0 in case the object is of rectangular shape or 1 in case the object is of circular shape. Depending on the value of s, the parameter a represents the diameter of the circle or the first edge length of the rectangle, while the length of the second edge is represented by the parameter b, which is unused in case of a circular shape.
The results for the fitting of the measured sinograms of both samples to the developed model are shown in table 1. At first, it is apparent from the results, that the optimization algorithm results in the correct shape for each sample object. In case of the rectangular sample object (Fig. 5(d)) the s parameter is 0, while the parameter is 1 in case of the circular object (Fig. 4(b)). However, the estimated size of the sample objects are too small in both cases, due to a multitude of different factors. Size affecting factors are the refractive index, the absorption coefficient of the sample object and the collection angle of the detector, which are all predetermined in the model. Also, the broad frequency range of the used terahertz radiation influences the measured sinogram, since parameters like the absorption coefficient and optical effects like diffraction depend on the wavelength. Here, an extension of the model from the currently used single frequency calculations towards calculations using a broad frequency spectrum is needed. So far, the presented method can calculate the sinogram of rectangular and circular objects of different sizes and thereby is a proof-of-concept demonstration of tomographic image reconstruction in the terahertz frequency range, which has to be extended in the future, to allow tomographic imaging of arbitrarily shaped sample objects. Also the required time to perform the projection measurements (about 1 hour) and the subsequent reconstruction (several hours) need to be decreased. One possibility to reduce the measurement time would be the compressive sensing concept [11, 13, 27]. Compressive sensing offers the possibility to perform a reconstruction of the image with a reduced number of measurements. The number of needed measurements can be reduced to 25 – 30% of the total number of pixels [11], which enables shorter measurement times. This of course is achieved by a trade-off between the number of measurements and the quality of the reconstructed image. By increasing the number of measurements, the image quality increases quite fast and with a higher amount of measurements slowly converges like it is shown in [11]. The time required for the object reconstruction could be reduced by performing the genetic algorithm on a graphics processing unit [28] and by a Fourier optics based wave propagation model [25].
Table 1. Results for the fitting of the measured sinograms to the developed model. The initial values of the fit parameters are a = 10 and b = 10 for both objects, while s has an initial value of 0, in case of the circular object and 1, in case of the rectangular object, in order to prevent a bias of the algorithm.
6. Conclusion
In conclusion, a 2D terahertz tomography concept based on single pixel Hadamard imaging has been successfully realized, demonstrating that raster scanning of the object and multi-pixel detector are not mandatory. The investigated sample objects consist of PP and are of circular and cuboid shape. For both samples a sinogram from 0° to 180° in 10° iterations has been recorded and the SART algorithm has been used to calculate the corresponding 2D image. These images clearly show the outer borders of the objects, but leave the interior partly (circular object) or fully (rectangular object) hollow. These distortions are not ascribable to the applied single pixel concept, rather than the optical effects of refraction, reflection and diffraction, which in general severely affect terahertz imaging. These effects have been considered in a developed model, reproducing the experimental observed sinograms, such that a fitting algorithm was capable to reconstruct the geometry of the sample objects using the model. Finally the presented concept can also be extended to 3D tomographic imaging by a 2D spatial modulation of the near-infrared laser diode, which is easily realized by the DMD, and further spectrally resolved imaging can be realized by using the full capability of the Fourier transform spectrometer.
Acknowledgment
We thank Daniel Laukhardt for his work on the initialization and programming of the digital micro mirror device.
References and links
1. H.-B. Liu, H. Zhong, N. Karpowicz, Y. Chen, and X.-C. Zhang, “Terahertz Spectroscopy and Imaging for Defense and Security Applications,” Proc. IEEE 95(8), 1514–1527 (2007).
2. K. J. Siebert, T. Löffler, H. Quast, M. Thomson, T. Bauer, R. Leonhardt, S. Czasch, and H. G. Roskos, “All-optoelectronic continuous wave THz imaging for biomedical applications,” Phys. Med. Biol. 47(21), 3743–3748 (2002).
3. P. H. Siegel, W. Dai, R. A. Kloner, M. Csete, and V. Pikov, “First millimeter-wave animal in vivo measurements of L-Glucose and D-Glucose: Further steps towards a non-invasive glucometer,” in Proceedings of 41st International Conference on Infrared, Millimeter, and Terahertz waves (IEEE, 2016).
4. R. Gente, S. F. Busch, E.-M. Stübling, L. M. Schneider, C. B. Hirschmann, J. C. Balzer, and M. Koch, “Quality Control of Sugar Beet Seeds With THz Time-Domain Spectroscopy,” IEEE Trans. THz Sci. Techn. 6(5), 754–756 (2016).
5. P. Dean, Y. L. Lim, A. Valavanis, R. Kliese, M. Nikolić, S. P. Khanna, M. Lachab, D. Indjin, Z. Ikonić, P. Harrison, A. D. Rakić, E. H. Linfield, and A. G. Davies, “Terahertz imaging through self-mixing in a quantum cascade laser,” Opt. Lett. 36(13), 2587–2589 (2011).
6. B. Recur, A. Younus, S. Salort, P. Mounaix, B. Chassagne, P. Desbarats, J.-P. Caumes, and E. Abraham, “Investigation on reconstruction methods applied to 3D terahertz computed tomography,” Opt. Express 19(6), 5105–5117 (2011).
7. R. Al Hadi, H. Sherry, J. Grzyb, Y. Zhao, W. Forster, H. M. Keller, A. Cathelin, A. Kaiser, and U. R. Pfeiffer, “A 1 k-Pixel Video Camera for 0.7–1.1 Terahertz Imaging Applications in 65-nm CMOS,” IEEE J. Solid-State Circuits 47(12), 2999–3012 (2012).
8. J. Zdanevičius, M. Bauer, S. Boppel, V. Palenskis, A. Lisauskas, V. Krozer, and H. G. Roskos, “Camera for high-speed THz imaging,” J. Infrared Millim. Te. 36(10), 986–997 (2015).
9. M. P. Edgar, G. M. Gibson, R. W. Bowman, B. Sun, N. Radwell, K. J. Mitchell, S. S. Welsh, and M. J. Padgett, “Simultaneous real-time visible and infrared video with single-pixel detectors,” Sci. Rep. 5, 10669 (2015).
10. B. Sun, M. P. Edgar, R. Bowman, L. E. Vittert, S. Welsh, A. Bowman, and M. J. Padgett, “3D Computational Imaging with Single-Pixel Detectors,” Science 340(6134), 844–847 (2013).
11. W. L. Chan, K. Charan, D. Takhar, K. F. Kelly, R. G. Baraniuk, and D. M. Mittleman, “A single-pixel terahertz imaging system based on compressed sensing,” Appl. Phys. Lett. 93(12), 121105 (2008).
12. D. Shrekenhamer, C. M. Watts, and W. J. Padilla, “Terahertz single pixel imaging with an optically controlled dynamic spatial light modulator,” Opt. Express 21(10), 12507–12518 (2013).
13. C. M. Watts, D. Shrekenhamer, J. Montoya, G. Lipworth, J. Hunt, T. Sleasman, S. Krishna, D. R. Smith, and W. J. Padilla, “Terahertz compressive imaging with metamaterial spatial light modulators,” Nat. Photonics 8, 605–609 (2014).
14. R. I. Stantchev, B. Sun, S. M. Hornett, P. A. Hobson, G. M. Gibson, M. J. Padgett, and E. Hendry, “Noninvasive, near-field terahertz imaging of hidden objects using a single-pixel detector,” Sci. Adv. 2(6), e1600190 (2016).
15. J. P. Guillet, B. Recur, L. Frederique, B. Bousquet, L. Caniono, I. Manek-Hönninger, P. Desbarats, and P. Mounaix, “Review of Terahertz Tomography Techniques,” J. Infrared Millim. Te. 35(4), 382–411 (2014).
16. A. Brahm, A. Wilms, M. Tymoshchuk, C. Grossmann, G. Notni, and A. Tünnermann, “Optical Effects at projection measurements for Terahertz tomography,” Opt. Laser Technol. 62, 49–57 (2014).
17. T. Mohr, S. Breuer, G. Giuliani, and W. Elsäßer, “Two-dimensional tomographic terahertz imaging by homodyne self-mixing,” Opt. Express 23(21), 27221–27229 (2015).
18. S. Mukherjee, J. Federici, P. Lopes, and M. Cabral, “Elimination of Fresnel Reflection Boundary Effects and Beam Steering in Pulsed Terahertz Computed Tomography,” J. Infrared Millim. Te. 34(9), 539–555 (2013).
19. W. L. Chan, J. Deibel, and D. M. Mittleman, “Imaging with terahertz radiation,” Rep. Prog. Phys. 70, 1325 (2007).
20. P. Weis, J. L. Garcia-Pomar, M. Höh, B. Reinhard, A. Brodyanski, and M. Rahm, “Spectrally Wide-Band Terahertz Wave Modulator Based on Optically Tuned Graphene,” ACS Nano 6(10), 9118–9124 (2012).
21. S. Busch, B. Scherger, M. Scheller, and M. Koch, “Optically controlled terahertz beam steering and imaging,” Opt. Lett. 37(8), 1391–1393 (2012).
22. W. van Aarle, W. J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. De Beenhouwer, K. J. Batenburg, and J. Sijbers, “Fast and flexible X-ray tomography using the ASTRA toolbox,” Opt. Express 24(22), 25129–25147 (2016).
23. J. B. Perraud, J. Bou Sleiman, B. Recur, H. Balacey, F. Simoens, J. P. Guillet, and P. Mounaix, “Liquid index matching for 2D and 3D terahertz imaging,” Appl. Optics 55(32), 9185–9192 (2016).
24. D. J. Butler and G. W. Forbes, “Fiber-diameter measurement by occlusion of a Gaussian beam,” Appl. Optics 37(13), 2598–2607 (1998).
25. D. G. Voelz, Computational Fourier Optics: A MATLAB® Tutorial, SPIE (2011).
26. B. Ung, A. Mazhorova, A. Dupuis, M. Rozé, and M. Skorobogatiy, “Polymer microstructured optical fibers for terahertz wave guiding,” Opt. Express 19(26), B848–B861 (2011).
27. S. Augustin, J. Hieronymus, P. Jung, and H.-W. Hübers, “Compressed Sensing in a Fully Non-Mechanical 350 GHz Imaging Setting,” J. Infrared Millim. Te. 36(5), 496–512 (2015).
28. S. Cavuoti, M. Garofalo, M. Brescia, A. Pescape, G. Longo, and G. Ventre, “Genetic Algorithm Modeling with GPU Parallel Computing Technology,” in Proceedings of Neural Nets and Surroundings: 22nd Italian Workshop on Neural Nets19, 29–39 (2013).
