Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Demonstration of a terahertz multi-spectral 3×3 Mueller matrix polarimetry system for 2D and 3D imaging

Open Access Open Access

Abstract

Mueller matrix polarimetry (MMP) has been demonstrated and recognized as an effective approach to attaining imaging enhancement as well as revealing polarization properties of an imaged sample. Generally, a minimum of 16 combinations of intensity-only measurements involving both linear and circular polarizations are required to completely and accurately determine the 4 × 4 Mueller matrix (MM) and comprehensively describe the polarization properties of the sample. However, broadband circular polarizations (CP) are rather difficult to obtain for design and fabrication limitations in the terahertz region, which poses a challenge to the acquisition of the 4 × 4 MM. In this circumstance, the 3 × 3 MM degradation using only linear polarizations (LP) is preferred and sufficient for characterization of non-depolarizing samples. In this paper, a multi-spectral 3 × 3 MMP system based on the THz time-domain spectroscopy (THz-TDS) is established from 0.1 to 1 THz. The system demonstrated is capable of fulfilling the accurate determination of the 3 × 3 MM. The Mueller matrix polar decomposition (MMPD), modified to be compatible with the MM degradation, is employed to explore the fine details and properties of the sample. By signal post-processing techniques, the MM elements in the time domain are retrieved, and the time dimension reflecting the depth information facilitates the 3D reconstruction of the sample. This work provides a prototype for 3D imaging of biological samples at higher frequencies in the future.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Studies on Mueller matrix polarimetry have come a long way and received considerable attention in recent years owing to its remarkable performance in characterizing the polarization properties of the sample compared with conventional imaging [13]. In addition, it was demonstrated a practical mechanism for achieving high-contrast imaging performance with fine details of the microstructures of the sample [46]. Substantial progress has been made on the study of MMP in a wide diversity of fields including remote sensing [7], botany [8], and especially biomedical applications [912]. Lu and Chipman proposed the Mueller matrix polar decomposition algorithm (MMPD), in which an MM is decomposed into three sub-matrices describing depolarization, retardance and diattenuation, respectively [13]. These fundamental polarization properties of a sample provide us with a new perspective to interpret the Mueller images. In general, a complete MM has 4 × 4 elements representing different meanings of polarization property of the sample, which could be further extracted using certain decomposition methods such as MMPD. A minimum of 16 combinations of intensity-only measurement involving both linear and circular polarizations are required to completely and accurately determine the 4 × 4 MM. However, CP is rather difficult to realize due to design and fabrication limitations in the terahertz regime, which makes broad spectral bandwidth for CP hard to obtain accordingly. This will bring about the inconsistency between LP and CP in terms of polarization purity, insertion loss and available bandwidth. This kind of inconsistency will produce error in elements concerning CP, which consequently propagate to cause chain errors in other elements on account that MM elements are correlated with each other for the case of 16 combinations [14]. These together could become major obstacles to the accurate acquisition of a complete MM. In this circumstance, the 3 × 3 MM degradation using only LP emerges as a preferred alternative for analysis of non-depolarizing samples. The 3 × 3 MM, exclusive of the last row and column of the conventional 4 × 4 MM, only involves LP and can be determined using only nine LP combinations of measurement. Without the spectral limitation of CP, the system can maintain excellent performance across an extended wide bandwidth. Several works have led the way and contributed to the 3 × 3 MMP [15,16]. Forward et al. and Qi et al. proposed narrow band endoscopic probes for 3 × 3 MM measurements and demonstrated their initial feasibility on ex vivo tissue [17,18]. Swami et al. presented a method based on polar decomposition for quantification of the different polarization parameters of a scattering medium, and elaborated the parameters comprehensively [19]. Wang et al. implemented the Monte Carlo simulations to compare the polarization parameters derived from the 3 × 3 and 4 × 4 decomposition methods, and revealed their similar qualitative relations in the microstructure of the sample, demonstrating the validity of 3 × 3 MMP [20]. It is noticed that these works are all proposed in optics and focused on the characterization of biological samples. Very few studies have been reported on the topic of terahertz MM experiments or measurements [21,22]. Extending MMP from optics to terahertz (THz) region should be expected, because THz radiation on biological samples may own some distinct advantages over optics for its sensitivity to water molecules and non-ionizing effect. Our previous work implemented the 4 × 4 MMP using 36 combinations of polarization involving both LP and CP at a single frequency of 0.3 THz [23].

In the present work, a multi-spectral 3 × 3 MMP system based on the THz-TDS is established from 0.1 to 1 THz for 2D and 3D reconstruction using only nine LP measurements, which significantly lessens the measurement time. It should be noticed that the MMP is an intensity-only formalism where only the intensity of the detection is recorded, which would reduce the complexity of measurement and post-processing in comparison with complex tensor analysis techniques [24]. The THz-TDS is configured in transmission mode and modified with our own polarization state generator (PSG) and polarization state analyzer (PSA) modules. The 3 × 3 MMs of the samples are calculated at multiple frequencies, and the modified MMPD is employed for extraction of polarization properties of the samples, which could map the fine details of the microstructure of the samples in a way. Meanwhile, the time-domain MMs are retrieved for 3D reconstruction of the samples.

This paper is organized as follows. Section 2 elaborates the 3 × 3 MM theory and its decomposition, followed by how the theory is specifically put into practice. Then, Section 3 explains in detail the experiments where the configuration, calibration and experimental procedure are elaborated in a row. Section 4 discusses about how signal post-processing is implemented to improve the imaging results and retrieve the 3D time-domain data, which is used to reconstruct the 3D structure of the sample. Afterwards, the results are discussed and interpreted. At last, the performance of the work is concluded in Section 5.

2. Theory and method

2.1 3×3 Mueller matrix

The Stokes vector describing the polarization state of electromagnetic radiation is generally composed of four components denoted as I, Q, U, V, which respectively represent the total intensity, the difference of intensities between horizontal and vertical states, 45° and 135° states, and the right-hand and left-hand CPs. Therefore, the fourth element of the Stokes vector could be ignored when circular polarizations are not taken into account, which would thus degrade the Stokes-Mueller formalism into its linear-polarization-only version

$${S_{out}} = \left[ {\begin{array}{c} {{I_{out}}}\\ {{Q_{out}}}\\ {{U_{out}}} \end{array}} \right] = M{S_{in}} = \left[ {\begin{array}{ccc} {{M_{11}}}&{{M_{12}}}&{{M_{13}}}\\ {{M_{21}}}&{{M_{22}}}&{{M_{23}}}\\ {{M_{31}}}&{{M_{32}}}&{{M_{33}}} \end{array}} \right]\left[ {\begin{array}{c} {{I_{in}}}\\ {{Q_{in}}}\\ {{U_{in}}} \end{array}} \right],$$
where Sin, Sout and M are the Stokes vectors of the incidence and reception and the Mueller matrix (MM) of the system including the PSG, PSA and measured sample, the relation between which could be written as
$$M = {M_{PSA}}{M_{sample}}{M_{PSG}}.$$
When different polarization combinations are measured, the MM of the sample could be calculated by solving the linear equation. In this convention, there are nine unknowns, i.e. the nine components of the 3 × 3 Mueller matrix to be determined, which means at least nine polarization combinations are required to solve the matrix accurately. In this work, 45° (P), 90° (V) and 135° (M) linear polarizations are adopted both for the PSG and PSA modules. The Stokes vectors for these polarizations and Mueller matrices for these polarizers in theory are enumerated in Fig. 1.

 figure: Fig. 1.

Fig. 1. Normalized Stokes vectors and Mueller matrices for 45° (P), 90° (V) and 135° (M) linear polarizations. (a) Stokes vectors for P, V, M polarizations; (b), (c), (d) Mueller matrices for P, V and M polarizations.

Download Full Size | PDF

For every incidence generated by the PSG, the Stokes vector is denoted as $S_{in}^{(m)}$, where m could be one of the three polarizations of P, V, and M. Therefore, the Stokes vector $S_{out}^{(n)}$ (n = P, V, M) of the detection could be expressed by the following Stokes-Mueller formalism

$$S_{out}^{(n)} = M_{PSA}^{(n)}{M_{sample}}M_{PSG}^{(m)}S_{in}^{(m)},$$
where Msample is the 3 × 3 MM of the sample, $M_{PSG}^{(m)}$ and $M_{PSA}^{(n)}$ are the MMs of the PSG which generates the polarization of m and that of the PSA which receives the polarization of n, respectively. In this paper, the PSG and PSA are simplified as linear polarizers illustrated in Figs. 1(b)–1(d). As previously described, only the first component I of the Stokes vector representing the total intensity of the radiation needs to be collected. For the incidence of m and reception of n, the intensity of detection is denoted as Imn. The linear algebra in Eq. (3) tells us that Imn is only relevant to the first row of $M_{PSA}^{(n)}$. If we build a matrix MPSA with each row being the first row of$M_{PSA}^{(n)}$written as$M_{PSA}^{(n,1)}$and a matrix MSout with each column representing the Stokes vector of the radiation before entering the PSA written as $S_{PSA}^{(p)}$ (p = P, V, M), the relation in Eq. (3) could be reshaped as
$${I_{meas}} = {M_{PSA}}{M_{Sout}},$$
where Imeas is the measured intensity matrix, and
$$\begin{array}{l} {M_{PSA}} = {\left[ {\begin{array}{ccc} {M_{PSA}^{(V,1)}}&{M_{PSA}^{(P,1)}}&{M_{PSA}^{(M,1)}} \end{array}} \right]^T},\\ {M_{Sout}} = \left[ {\begin{array}{ccc} {S_{PSA}^{(V)}}&{S_{PSA}^{(P)}}&{S_{PSA}^{(M)}} \end{array}} \right], \end{array}$$
Thus, MSout could be inverted as MSout = (MPSA)-1Imeas from Eq. (4). One might also notice that $S_{PSA}^{(p)} = {M_{sample}}M_{PSG}^{(p)}S_{in}^{(p)}$, therefore the MM of the sample could be retrieved as
$${M_{sample}} = {({{M_{PSA}}} )^{ - 1}}I{({M_{Sin}})^{ - 1}},$$
where ${M_{Sin}} = [M_{PSG}^{(V)}S_{in}^{(V)},M_{PSG}^{(P)}S_{in}^{(P)},M_{PSG}^{(M)}S_{in}^{(M)}]$. In our practice using polarizations of V, P and M, the composite matrices for PSG and PSA are calculated as illustrated in Fig. 2.

 figure: Fig. 2.

Fig. 2. The composite matrices for PSG and PSA using V, P, M polarizations. (a) Composite matrix for PSG; (b) Composite matrix for PSA.

Download Full Size | PDF

It should be aware of that each pixel of the scanned sample corresponds to an MM, and all the MM elements of the sample form a 2D MM image of the sample.

2.2 3×3 Mueller matrix polar decomposition

Generally, the MM could be further decomposed to acquire the polarization parameters of the sample, which could characterize the properties as well as reveal the fine structures of the sample. Several decomposition methods have been proposed in previous works, including the MMPD, reverse polar decomposition [25], symmetric decomposition [26], differential decomposition [27,28], and Mueller matrix transformation (MMT) [4,29]. Of all these methods, the MMPD gives the most comprehensive and widely applicable interpretation of the sample, extracting the three fundamental parameters of depolarization, diattenuation and retardance by decomposing the MM into three sub-matrices representing the three parameters, respectively. Depolarization describes the ratio of incident wave that becomes unpolarized upon interaction with the sample, while diattenuation characterizes the dependence of transmission (or reflection) upon the incident polarization state. Retardance refers to the phenomenon of phase difference between two orthogonally polarized components of radiation when propagating through the sample. In this paper, the MMPD is modified to adapt to the 3 × 3 MMP. For the 3 × 3 MMPD, a Mueller matrix is similarly decomposed into three sub-matrices of a depolarizer, a retarder and a diattenuator in sequence written as

$$M = {M_\varDelta }{M_R}{M_D},$$
where M is the Mueller matrix, MΔ, MR, and MD denote the depolarizer, retarder and diattenuator, respectively. The three parameters of depolarization, diattenuation and retardance denoted as Δ, D and R respectively could be further broken down to their linear and circular components. For the 3 × 3 MMPD, the circular components will be degraded except for the retardance. For brevity and clarity, the decomposition process will be omitted, and the expressions of the parameters are simply listed as
$$ \begin{aligned}{D_L} &= \sqrt {M_{12}^2 + M_{13}^2} /{M_{11}},\\ {{\Delta }_L} &= 1 - \left( {\left| {{m_{{\Delta }(11)}}} \right| + \left| {{m_{{\Delta }(22)}}} \right|} \right)/2,\\ {R_L} &= {\cos ^{ - 1}}\left( {\sqrt {\begin{array}{l} {{{\left[ {{M_R}(2,2) + {M_R}(3,3)} \right]}^2}}\\ { + {{\left[ {{M_R}(3,2) - {M_R}(2,3)} \right]}^2}} \end{array}} - 1} \right) \\{R_C} &= {\tan ^{ - 1}}\frac{{{M_R}(3,2) - {M_R}(2,3)}}{{{M_R}(2,2) + {M_R}(3,3)}}, \end{aligned} $$
where the subscripts “L” and “C” represent the linear and circular components of the parameters, respectively. And mΔ is the sub-matrix of MΔ with only its last two rows and columns.

3. Experiments

3.1 Experiment setup

The experimental configuration is based on the THz-TDS, which is a spectroscopic technique where properties of the sample are explored with short pulses of terahertz radiation. The THz-TDS directly measures the time domain pulses of transient electric fields passing through the sample. The measured THz electric field could then be Fourier transformed to yield the spectral response of the sample. Typically, the probed THz electric field at the detector is on the order of 10∼100 V/cm and has a time duration of a few picoseconds [30]. Our model is T-Spec series from EKSPLA, and the detected pulse signal is in the form of electric current, which could reach as high as 200∼300 nA when there is no sample placed on the stage. Figure 3(a) shows the typical time domain pulses comparison with and without sample inserted between the emitter and detector. It could be observed that the pulse of the sample would present time delay as well as attenuation in amplitude compared with that of air. This time delay reflects the depth information of the sample, which will be used to retrieve the 3D structure of the sample. The time domain pulses in Fig. 3(a) are then Fourier-transformed to obtain the spectrums illustrated in Fig. 3(b). The usable spectrum could extend from 0.1 to 1 THz. Figure 3(c) shows the mean pulse of 10 measurements and their standard deviation (STDEV), from which the signal-to-noise ratios (SNR) [31] are calculated to be as high as 100 plus, which is sufficiently strong for measurement of samples of ordinary thickness.

 figure: Fig. 3.

Fig. 3. Typical signal comparison between air and the sample in (a) time domain and (b) frequency domain; (c) The SNR of the system.

Download Full Size | PDF

The experimental setup in transmission configuration is illustrated in Fig. 4. The Pump Source provides pulses of 760∼840 nm wavelength, 50∼150 fs pulse duration and 50∼150 mW output power at approximate 80 MHz pulse repetition rate. The beam height adapter (BHA) is used to adjust the beam height. The single ultrashort pulse is then divided by the pellicle beam splitter (BS1) into two separate beams with the splitting ratio of 55:45 called pumping beam and probing beam, respectively. The pumping beam goes through the fast delay line (HLR1) to the emitter guided by the mirror (M2), and the lens (L1) is used to adjust the power of the pumping beam entering the emitter. The probing beam goes through the slow delay line (HLR2) via the prism (PR3) to the detector guided by the mirror (M3) instead. The THz emitter, which is a photo-conductive antenna, would then transform the pumping beam into THz radiation of sub-picosecond pulses with vertical linear polarization, which would be converted by our self-designed linear polarizer (LP1, PSG) to produce three polarizations (V, P, M) and then focused by the parabolic mirror (M4) to the sample. The beam passing through the sample will be guided by the parabolic mirror (M5) via LP2 (PSA) to the detector. For brevity, details of the LP lens will be given in Supplement 1. The focal length of the mirrors (M4, M5) is 50.8 mm, and the averaged spot size of them is about 2 mm. The THz signal is detected by the digital signal processing (DSP) card integrated into the control unit. The pulse waveform of the THz radiation is built by scanning the fast delay line in 10 Hz frequency. The sample, fixed on the X-Y stage, is raster-scanned in the x-y plane to acquire the 2D transmission data. For each scan point, the time domain pulse passing through the point is recorded, eventually forming the 3D data, which will be post-processed for imaging of the sample.

 figure: Fig. 4.

Fig. 4. The experimental setup in transmission geometry. “BHA” is the beam height adapter, “BS” means beam splitter, “HLR” means hollow retro-reflector, “PR” means prism, “M” means mirror, “ID” means iris diaphragm, “L” means lens, and “LP” means linear polarizer. LP1 and LP2 are the self-designed linear polarizers acting as the PSG and PSA respectively.

Download Full Size | PDF

3.2 Calibration

Due to the impurity, fabrication errors and insertion loss, etc., there is no way that the polarizers can be perfect. And the imperfections will give rise to errors in the process of the measurement, thus eventually contributing to the inaccuracy for extraction of the MM, as well as the related property parameters of the sample. As a consequence, a calibration procedure needs to be carried out to make up for the imperfections of the polarizers. The procedure will take advantage of the transmission properties in the air and is composed of three steps : the first step is to record the transmission data for each single polarization combination when there is no sample (air) placed at the stage; the second is to theoretically calculate the transmittances in the air for all polarization combinations, which could be realized by considering the MM of air as an identity matrix; finally, the measured transmittances are compared with the theoretical ones to obtain the compensations, which will be implemented in real measurements of the sample for better imaging performance. Figure 5 shows the transmittances of the nine polarization combinations with and without calibration evaluated with the theoretical ones from 0.2 to 1 THz. It could be observed that the performance of the polarizers without calibration is less satisfactory, but when the calibration algorithm is employed, the calibrated data are in excellent agreement with the theory. To validate the feasibility of the measurement system and the calibration approach, a 45° linear polarizer is used as the test sample to evaluate the accuracy for MM obtainment. The calculated MMs with and without calibration are illustrated in Fig. 6. The simulations are proceeded in COMSOL Multiphysics using the Ray Optics module, and the result in Fig. 6(e) presents perfect agreement with the theory in Fig. 6(a). Figures 6(b)–6(d) show the raw MMs for the 45° polarizer without calibration at three specific frequencies, where the errors for some elements like M21 and M23 are as high as 20%, and are exacerbated to 38% at higher frequencies. While after the calibration, the averaged error could be reduced to within 5% for the same frequencies, which is considerably negligible in the process of measurement. The measurements of the sample could be all set when the system is well calibrated. In our experiment, a 3D-printing sample is studied and elaborated in Subsection 3.3.

 figure: Fig. 5.

Fig. 5. Transmittances of the nine polarization combinations with and without calibration evaluated with the theoretical ones from 0.2 to 1 THz.

Download Full Size | PDF

 figure: Fig. 6.

Fig. 6. MMs for 45° polarizer. (a) theoretical MM; (e) simulated MM; (b)-(d) raw MMs at 0.57 THz, 0.75 THz, and 0.92 THz respectively; (f)-(h) calibrated MMs at 0.57 THz, 0.75 THz, and 0.92 THz respectively.

Download Full Size | PDF

3.3 Measurement of the sample

The sample we use is a 3D-printing sample illustrated in Fig. 7(a), the structure of which is inspired by the oriented fibrous collagen with retardance property. The sample is 3D-printed using high temperature resin of ɛr = 2.66 with the honeycomb-like hollow structure shown in Fig. 7(a). The thickness of the sample is 3 mm, and the hollow pattern is formed by connected circular tubings with different diameters centered at the central plane. There are three eigen-structures in the pattern: one is indicated by the red triangle with the diameter of 2 mm, the other two are indicated by the green and purple arrows of 1 mm and 0.7 mm diameters, respectively. The scan area is restricted within the blue square with the side length 23 mm. The scan steps along x- and y-axis are both 0.27 mm. The whole scan includes 85 × 85 pixels, which takes about 80 minutes for one combination of measurement. The slower raster scan provides more uniform imaging quality than typical beam deflection approach used in optical experiment. Figure 7(b) presents the temporal maximum value mapping at each pixel for the nine polarization combinations normalized by the VV measurement, which shows the general structure of the sample. Figure 7(c) shows the transmittances of the nine combinations at 0.44 THz normalized by the transmittance of VV, indicating more details of the sample compared with Fig. 7(b).

 figure: Fig. 7.

Fig. 7. (a) the profile of the sample; (b) the temporal maximum mapping for nine polarization combinations; (c) the transmittances of nine polarization combinations at 0.44 THz.

Download Full Size | PDF

4. Post-processing and results

4.1 Signal post-processing

The raw data collected are the results of convolution among the sample, the system and the ambient environment. Therefore, in order to eliminate the influence of the system errors and the ambient environment on the imaging performance, certain signal processing approaches are employed for data deconvolution. In this work, the detected time domain pulse transmitted through the sample y(t) are expressed as

$$y(t) = x(t) \ast h(t) = \int_{ - \infty }^\infty {x(t - \tau )} h(\tau )d\tau ,$$
where x(t) and h(t) are the detected signal when no sample (air) is inserted and the impulse response of the sample, respectively. Thus, the transfer function in the frequency domain could be written as
$$H(\omega ) = {{Y(\omega )} / {X(\omega )}}.$$
In addition, the Weiner filter W(ω) [32] and double Gaussian bandpass filter G(ω) [33] are used to avoid divergence and denoise the lower and higher frequency noise expressed as
$$\begin{array}{l} W(\omega ) = \frac{{\bar{X}(\omega )}}{{{{|{X(\omega )} |}^2} + 1/SNR}},\\ G(\omega ) = b{e^{ - {\omega ^2}/2\alpha _2^2}} - a{e^{ - {\omega ^2}/2\alpha _1^2}},\;\;\;\;b \ge a,\;{\alpha _2} > {\alpha _1}, \end{array}$$
where SNR is the signal-to-noise ratio, α1 and α2 are the lower and upper cutoff frequencies set as 0.1 and 1 THz respectively, a and b are constants set as a = b = 1. And the deconvolved pulse in time domain could be expressed as
$$h(t) = {F^{ - 1}}\{{H(\omega )W(\omega )G(\omega )} \},$$
where F-1{–} is an inverse Fourier transform operator. Figure 8 shows the typical pulses before and after deconvolution, where the insets correspond to their spectrums respectively. It could be seen that the signal after deconvolution removes the lower and higher frequency noise, and presents a smoother curve. The deconvolution process will be applied to every pixel of the sample.

 figure: Fig. 8.

Fig. 8. The time domain pulses before and after deconvolution along with their spectrums. (a) before deconvolution; (b) after deconvolution.

Download Full Size | PDF

4.2 Results and discussions

After deconvolution, the data are used to calculate the Mueller matrix across 0.2∼1 THz on the basis of the method presented in Subsection 2.1. For an isotropic material like our 3D-printing sample, only diagonal elements of the Mueller matrix are occupied with about the same intensities [34], while the off-diagonal ones are considerably smaller to be neglected. The averaged MM elements of all pixels of the sample are calculated and illustrated in Fig. 9. The results in Fig. 9 demonstrate that the 3D-printing sample is nearly isotropic. The normalized MM values are almost independent of frequency across 0.3∼0.9 THz. In spite of this, when the frequency goes higher, the imaging resolution of MMP could be enhanced and the fine details of the sample could be revealed at different frequencies. The diagonal elements of the raw MM are then integrated over several frequency intervals of 0.1THz to improve the resolution and displayed in Fig. 10. The contrasts of the images are much better enhanced compared with that of the images at single frequency points. With the frequency goes higher, the power drops with increasing noise, but the images present more and more clear and sharp edges of the structure. The narrowest channels (0.7 mm wide) denoted along the purple arrow in Fig. 7(a) are not much clearer as the others but can still be clearly observed, especially at higher frequencies. The hollow cavity indicated by the red triangle in Fig. 7(a) emerges over 0.7∼0.8 THz, which cannot be observed at lower frequencies. In general, the Mueller matrix is normalized by M11 pixelwise before the MMPD. The normalized MM elements are listed in Fig. 11, where the fine eigen-structures of the sample could be separated by frequency intervals. Over 0.2∼0.3 THz, only the hollow cavities could be seen in M33 and M21, while the green channels (1 mm in diameter) could be clearly exhibited in M32. Over 0.3∼0.4 THz, the hollow cavities are clearer, and M23 and M21 show the purple channels (0.7 mm in diameter) and green channels, respectively. When the frequency goes higher, the SNR decreases but the fine structures are sharper. Over 0.6∼0.7 THz, M33 and M23 present the hollow cavities of the sample. Generally, M33, M21, M12, and M23 (or M32) are respectively related to the depolarization, polarizance, diattenuation and anisotropy of the sample [4], which could be further extracted using the MMPD method.

 figure: Fig. 9.

Fig. 9. The averaged MM elements of all pixels of the sample versus frequency from 0.3 to 0.9 THz.

Download Full Size | PDF

 figure: Fig. 10.

Fig. 10. Integrated diagonal elements of the raw MM over frequency intervals of 0.1 THz from 0.3∼0.4 THz to 0.7∼0.8 THz.

Download Full Size | PDF

 figure: Fig. 11.

Fig. 11. The integrated normalized elements over frequency intervals of 0.1 THz from 0.2∼0.3 THz to 0.6∼0.7 THz.

Download Full Size | PDF

Since the sample is non-depolarizing, the depolarization parameter may not be a good indicator for imaging. The MMPD parameters except for depolarization over several frequency intervals are illustrated in Fig. 12. The reported literature [4,23] has demonstrated that the retardance is very sensitive to the oriented structure of the sample. The linear component of retardance generally resonates with the oriented but symmetric structure, while the circular component echoes with the oriented but asymmetric structure like our sample. Therefore, the circular retardance shows sharper contrast than the linear one. The linear diattenuation, characterizing the intensity attenuation of linear polarizations, shows clear outline of the hollow cavities at lower frequencies and finer channels at higher frequencies.

 figure: Fig. 12.

Fig. 12. The MMPD parameters over frequency intervals of 0.3∼0.4 THz and 0.6∼0.7 THz.

Download Full Size | PDF

4.3 3D reconstruction

The time delay information in the 3D deconvolved data reveals the optical length, which reflects the depth information of the sample in a way [32,35]. In our approach, the time delay data are used to map the third dimension (depth) of the 3D-printing sample. Figure 13(a) shows the main pulse distribution of all pixels of the 3D printing sample, and two characteristic peaks are distributed over 2∼4 ps and 6∼8 ps respectively, representing the two distinct structures of the sample. When time delay is confined within 2∼4 ps indicated by the blue frame Fig. 13(a), the peak corresponds to the hollow cavities of the sample, the data of which are used for 3D reconstruction illustrated in Figs. 13(b) and 13(c) from side and back views, respectively. It could be seen that only the cavities are retrieved and the hollow feature can be spotted in green and yellow. When time delay is restricted within 6∼8 ps framed by the green rectangle in Fig. 13(a), the main peak corresponds to the whole structure of the sample and the reconstruction is shown in Figs. 13(d) and 13(e), where all the channels as well as the hollow cavities could be observed and the orange hexagonal prisms surrounded by the channels could also be outlined. The time delay cross sections and B-scans are displayed in Figs. 13(f)–13(i) to show the fine inner-structures of the sample. Figure 13(f) shows the cross section at t = 3.04 ps, where the hollow cavities appear very clear. And the cross section at t = 7.04 ps exhibits the general structure of the sample similar to Fig. 10. The B scans of X = 12 mm and Y = 12 mm positioned by red and blue frames in Fig. 13(e) are presented in Figs. 13(h)–13(i), which show the distribution of hexagonal prisms along Y and X positions, respectively.

 figure: Fig. 13.

Fig. 13. The 3D reconstruction of the 3D printing sample. (a) the main pulses distribution of all pixels; (b) the side view of the 3D image during 2∼4 ps; (c) the back view of the 3D image during 2∼4 ps; (d) the back view of the 3D image during 6∼8 ps; (e) the side view of the 3D image during 6∼8 ps; (f) the cross-section of the 3D image when t = 3.04 ps; (g) the cross-section of the 3D image when t = 7.04 ps; (h) the B-scan of Y/t profile when X = 12 mm; (i) the B-scan of X/t profile when Y = 12 mm.

Download Full Size | PDF

5. Conclusion

In this paper, a terahertz multi-spectral 3 × 3 Mueller matrix polarimetry system in transmission configuration from 0.1 to 1 THz is established and demonstrated, capable of achieving 2D and 3D imaging. The high-performance LP polarizers are incorporated into the THz-TDS to generate different linear polarizations with focusing beams. The raw polarization data are measured to determine the 3 × 3 Mueller matrix of the sample, followed by the MMPD to extract the polarization properties of the sample. The diagonal MM elements show much sharper contrast of the sample compared with the raw polarization images. The normalized MM elements by M11 are demonstrated to show the separated eigen-structures respectively. The MMPD parameters can not only show the fine details of the sample at different frequencies, but also characterize the polarization properties of the sample. The circular retardance and linear diattenuation show the most prominent performance among all the parameters, which demonstrates the fact that the circular retardance is very sensitive to the oriented but asymmetric structure of our sample. Generally, the contrast of the image is better enhanced and more fine details emerge with the increasing frequency. The multi-spectral results demonstrate a good way to improve the resolution of imaging and exhibit the whole picture of the sample. At last, the multi-spectral data of the Mueller matrix are post-processed and inverted to produce the 3D time domain data, in which the time delay information could facilitate the retrieval of the depth profile of the sample, with which the 3D view of the sample could be reconstructed to provide an overall perspective of the sample. We are making efforts to design high performance wideband LP-to-CP converter to extract a more comprehensive Mueller matrix for imaging of more complex samples like biological tissues. This work provides a prototype and lays the foundation for what follows in our next-step in 2D and 3D biomedical imaging.

Funding

Research Grants Council, University Grants Committee (CityU11217720, General Research Fund CityU11217720, T42-103/16-N, Theme Based Research Scheme T42-103/16-N).

Disclosures

The authors declare no conflicts of interest.

Supplemental document

See Supplement 1 for supporting content.

References

1. R. A. Chipman, “Polarimetry,” in Handbook of Optics, M. Bass, ed. (McGraw-Hill, 1995).

2. E. Garcia-Caurel, R. Ossikovski, M. Foldyna, A. Pierangelo, B. Drevillon, and A. D. Martino, “Advanced Mueller Ellipsometry Instrumentation and Data Analysis,” in Ellipsometry at the Nanoscale, M. Losurdo and K. Hingerl, eds. (Springer, 2013).

3. J. J. Gil and R. Ossikovski, Polarized Light and the Mueller Matrix Approach (CRC, 2016).

4. M. H. Sun, H. H. He, N. Zeng, E. Du, Y. H. Guo, S. X. Liu, J. Wu, Y. H. He, and H. Ma, “Characterizing the microstructures of biological tissues using Mueller matrix and transformed polarization parameters,” Biomed. Opt. Express 5(12), 4223–4234 (2014). [CrossRef]  

5. D. S. Chen, N. Zeng, Q. L. Xie, H. H. He, V. V. Tuchin, and H. Ma, “Mueller matrix polarimetry for characterizing microstructural variation of nude mouse skin during tissue optical clearing,” Biomed. Opt. Express 8(8), 3559–3570 (2017). [CrossRef]  

6. B. Liu, Y. Yao, R. Q. Liu, H. Ma, and L. Ma, “Mueller polarimetric imaging for characterizing the collagen microstructures of breast cancer tissues in different genotypes,” Opt. Commun. 433, 60–67 (2019). [CrossRef]  

7. S. R. Cloude, Polarisation: Applications in Remote Sensing (Oxford University, 2009), pp. 91–98.

8. C. H. L. Patty, D. A. Luo, F. Snik, F. Ariese, W. J. Buma, I. L. T. Kate, R. J. M. V. Spanning, W. B. Sparks, T. A. Germer, G. Garab, and M. W. Kudenov, “Imaging linear and circular polarization features in leaves with complete Mueller matrix polarimetry,” Biochim. Biophys. Acta 1862(6), 1350–1363 (2018). [CrossRef]  

9. M. R. Antonelli, A. Pierangelo, T. Novikova, P. Validire, A. Benali, B. Gayet, and A. D. Martino, “Mueller matrix imaging of human colon tissue for cancer diagnostics: how Monte Carlo modeling can help in the interpretation of experimental data,” Opt. Express 18(10), 10200–10208 (2010). [CrossRef]  

10. M. Dubreuil, P. Babilotte, L. Martin, D. Sevrain, S. Rivet, Y. L. Grand, G. L. Brun, B. Turlin, and B. L. Jeune, “Mueller matrix polarimetry for improved liver fibrosis diagnosis,” Opt. Lett. 37(6), 1061–1063 (2012). [CrossRef]  

11. E. Du, H. H. He, N. Zeng, M. H. Sun, Y. H. Guo, J. Wu, S. X. Liu, and H. Ma, “Mueller matrix polarimetry for differentiating characteristic features of cancerous tissues,” J. Biomed. Opt. 19(7), 076013 (2014). [CrossRef]  

12. S. Alali and I. A. Vitkin, “Polarized light imaging in biomedicine: emerging Mueller matrix methodologies for bulk tissue assessment,” J. Biomed. Opt. 20(6), 061104 (2015). [CrossRef]  

13. S. Y. Lu and R. A. Chipman, “Interpretation of Mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13(5), 1106–1113 (1996). [CrossRef]  

14. N. Ghosh, A. Banerjee, and J. Soni, “Turbid medium polarimetry in biomedical imaging and diagnosis,” Eur. Phys. J.: Appl. Phys. 54(3), 30001 (2011). [CrossRef]  

15. Y. F. Fu, Z. W. Huang, H. H. He, H. Ma, and J. Wu, “Flexible 3 × 3 Mueller Matrix Endoscope Prototype for Cancer Detection,” IEEE Trans. Instrum. Meas. 67(7), 1700–1712 (2018). [CrossRef]  

16. M. K. Swami, H. S. Patel, and P. K. Gupta, “Conversion of 3 × 3 Mueller matrix to 4 × 4 Mueller matrix for non-depolarizing samples,” Opt. Commun. 286, 18–22 (2013). [CrossRef]  

17. S. Forward, A. Gribble, S. Alali, A. A. Lindenmaier, and I. A. Vitkin, “Flexible polarimetric probe for 3 × 3 Mueller matrix measurements of biological tissue,” Sci. Rep. 7(1), 11958 (2017). [CrossRef]  

18. J. Qi, M. L. Ye, M. Singh, N. T. Clancy, and D. S. Elson, “Narrow band 3 × 3 Mueller polarimetric endoscopy,” Biomed. Opt. Express 4(11), 2433–2449 (2013). [CrossRef]  

19. M. K. Swami, S. Manhas, P. Buddhiwant, N. Ghosh, A. Uppal, and P. K. Gupta, “Polar decomposition of 3 × 3 Mueller matrix: a tool for quantitative tissue polarimetry,” Opt. Express 14(20), 9324–9337 (2006). [CrossRef]  

20. Y. F. Wang, Y. H. Guo, N. Zeng, D. S. Chen, H. H. He, and H. Ma, “Study on the validity of 3 × 3 Mueller matrix decomposition,” J. Biomed. Opt. 20(6), 065003 (2015). [CrossRef]  

21. Z. Ahmed, T. Dmitri, Y. Colin, H. Martin, E. Henry, P. Matteo, and K. Junichiro, “Carbon nanotube fiber terahertz polarizer,” Appl. Phys. Lett. 108(24), 242907 (2016). [CrossRef]  

22. T. Hofmann, C. M. Herzinger, A. Boosalis, T. E. Tiwald, J. A. Woollam, and M. Schubert, “Variable-wavelength frequency-domain terahertz ellipsometry,” Rev. Sci. Instrum. 81(2), 023101 (2010). [CrossRef]  

23. S. X. Huang, Y. S. Zeng, G. B. Wu, K. F. Chan, B. J. Chen, M. Y. Xia, S. W. Qu, and C. H. Chan, “Terahertz Mueller Matrix Polarimetry and Polar Decomposition,” IEEE Trans. THz Sci. Technol. 10(1), 74–84 (2020). [CrossRef]  

24. R. Peng, O. Yurduseven, T. Fromenteze, and D. R. Smith, “Advanced Processing of 3D Computational Microwave Polarimetry Using a Near-Field Frequency-Diverse Antenna,” IEEE Access 8, 166261–166272 (2020). [CrossRef]  

25. R. Ossikovski, A. D. Martino, and S. Guyot, “Forward and reverse product decompositions of depolarizing Mueller matrices,” Opt. Lett. 32(6), 689–691 (2007). [CrossRef]  

26. R. Ossikovski, “Analysis of depolarizing Mueller matrices through a symmetric decomposition,” J. Opt. Soc. Am. A 26(5), 1109–1118 (2009). [CrossRef]  

27. N. Ortega-Quijano and J. L. Arce-Diego, “Mueller matrix differential decomposition,” Opt. Lett. 36(10), 1942–1944 (2011). [CrossRef]  

28. N. Ortega-Quijano, B. Haj-Ibrahim, E. García-Caurel, J. L. Arce-Diego, and R. Ossikovski, “Experimental validation of Mueller matrix differential decomposition,” Opt. Express 20(2), 1151–1163 (2012). [CrossRef]  

29. N. Zeng, X. Y. Jiang, Q. Gao, Y. H. He, and H. Ma, “Linear polarization difference imaging and its potential applications,” Appl. Opt. 48(35), 6734–6739 (2009). [CrossRef]  

30. J. Neu and C. A. Schmuttenmaer, “Tutorial: An introduction to terahertz time domain spectroscopy (THz-TDS),” J. Appl. Phys. 124(23), 231101 (2018). [CrossRef]  

31. W. Withayachumnankul and M. Naftaly, “Fundamentals of Measurement in Terahertz Time-Domain Spectroscopy,” J. Infrared, Millimeter, Terahertz Waves 35(8), 610–637 (2014). [CrossRef]  

32. J. Takayanagi, H. Jinno, S. Ichino, K. Suizu, M. Yamashita, T. Ouchi, S. Kasai, H. Ohtake, H. Uchida, N. Nishizawa, and K. Kawase, “High-resolution time-of-flight terahertz tomography using a femtosecond fiber laser,” Opt. Express 17(9), 7533–7539 (2009). [CrossRef]  

33. R. M. Woodward, B. E. Cole, V. P. Wallace, R. J. Pye, D. D. Arnone, E. H. Linfield, and M. Pepper, “Terahertz pulse imaging in reflection geometry of human skin cancer and skin tissue,” Phys. Med. Biol. 47(21), 3853–3863 (2002). [CrossRef]  

34. J. Freudenthal, “Intuitive interpretation of Mueller matrices of transmission,” Tech. Rep. (Hinds Instruments, Inc., 2018).

35. V. P. Wallace, E. MacPherson, J. A. Zeitler, and C. Reid, “Three-dimensional imaging of optically opaque materials using nonionizing terahertz radiation,” J. Opt. Soc. Am. A 25(12), 3120–3133 (2008). [CrossRef]  

Supplementary Material (1)

NameDescription
Supplement 1       Supplemental document for the self-designed linear polarizer

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (13)

Fig. 1.
Fig. 1. Normalized Stokes vectors and Mueller matrices for 45° (P), 90° (V) and 135° (M) linear polarizations. (a) Stokes vectors for P, V, M polarizations; (b), (c), (d) Mueller matrices for P, V and M polarizations.
Fig. 2.
Fig. 2. The composite matrices for PSG and PSA using V, P, M polarizations. (a) Composite matrix for PSG; (b) Composite matrix for PSA.
Fig. 3.
Fig. 3. Typical signal comparison between air and the sample in (a) time domain and (b) frequency domain; (c) The SNR of the system.
Fig. 4.
Fig. 4. The experimental setup in transmission geometry. “BHA” is the beam height adapter, “BS” means beam splitter, “HLR” means hollow retro-reflector, “PR” means prism, “M” means mirror, “ID” means iris diaphragm, “L” means lens, and “LP” means linear polarizer. LP1 and LP2 are the self-designed linear polarizers acting as the PSG and PSA respectively.
Fig. 5.
Fig. 5. Transmittances of the nine polarization combinations with and without calibration evaluated with the theoretical ones from 0.2 to 1 THz.
Fig. 6.
Fig. 6. MMs for 45° polarizer. (a) theoretical MM; (e) simulated MM; (b)-(d) raw MMs at 0.57 THz, 0.75 THz, and 0.92 THz respectively; (f)-(h) calibrated MMs at 0.57 THz, 0.75 THz, and 0.92 THz respectively.
Fig. 7.
Fig. 7. (a) the profile of the sample; (b) the temporal maximum mapping for nine polarization combinations; (c) the transmittances of nine polarization combinations at 0.44 THz.
Fig. 8.
Fig. 8. The time domain pulses before and after deconvolution along with their spectrums. (a) before deconvolution; (b) after deconvolution.
Fig. 9.
Fig. 9. The averaged MM elements of all pixels of the sample versus frequency from 0.3 to 0.9 THz.
Fig. 10.
Fig. 10. Integrated diagonal elements of the raw MM over frequency intervals of 0.1 THz from 0.3∼0.4 THz to 0.7∼0.8 THz.
Fig. 11.
Fig. 11. The integrated normalized elements over frequency intervals of 0.1 THz from 0.2∼0.3 THz to 0.6∼0.7 THz.
Fig. 12.
Fig. 12. The MMPD parameters over frequency intervals of 0.3∼0.4 THz and 0.6∼0.7 THz.
Fig. 13.
Fig. 13. The 3D reconstruction of the 3D printing sample. (a) the main pulses distribution of all pixels; (b) the side view of the 3D image during 2∼4 ps; (c) the back view of the 3D image during 2∼4 ps; (d) the back view of the 3D image during 6∼8 ps; (e) the side view of the 3D image during 6∼8 ps; (f) the cross-section of the 3D image when t = 3.04 ps; (g) the cross-section of the 3D image when t = 7.04 ps; (h) the B-scan of Y/t profile when X = 12 mm; (i) the B-scan of X/t profile when Y = 12 mm.

Equations (12)

Equations on this page are rendered with MathJax. Learn more.

$${S_{out}} = \left[ {\begin{array}{c} {{I_{out}}}\\ {{Q_{out}}}\\ {{U_{out}}} \end{array}} \right] = M{S_{in}} = \left[ {\begin{array}{ccc} {{M_{11}}}&{{M_{12}}}&{{M_{13}}}\\ {{M_{21}}}&{{M_{22}}}&{{M_{23}}}\\ {{M_{31}}}&{{M_{32}}}&{{M_{33}}} \end{array}} \right]\left[ {\begin{array}{c} {{I_{in}}}\\ {{Q_{in}}}\\ {{U_{in}}} \end{array}} \right],$$
$$M = {M_{PSA}}{M_{sample}}{M_{PSG}}.$$
$$S_{out}^{(n)} = M_{PSA}^{(n)}{M_{sample}}M_{PSG}^{(m)}S_{in}^{(m)},$$
$${I_{meas}} = {M_{PSA}}{M_{Sout}},$$
$$\begin{array}{l} {M_{PSA}} = {\left[ {\begin{array}{ccc} {M_{PSA}^{(V,1)}}&{M_{PSA}^{(P,1)}}&{M_{PSA}^{(M,1)}} \end{array}} \right]^T},\\ {M_{Sout}} = \left[ {\begin{array}{ccc} {S_{PSA}^{(V)}}&{S_{PSA}^{(P)}}&{S_{PSA}^{(M)}} \end{array}} \right], \end{array}$$
$${M_{sample}} = {({{M_{PSA}}} )^{ - 1}}I{({M_{Sin}})^{ - 1}},$$
$$M = {M_\varDelta }{M_R}{M_D},$$
$$ \begin{aligned}{D_L} &= \sqrt {M_{12}^2 + M_{13}^2} /{M_{11}},\\ {{\Delta }_L} &= 1 - \left( {\left| {{m_{{\Delta }(11)}}} \right| + \left| {{m_{{\Delta }(22)}}} \right|} \right)/2,\\ {R_L} &= {\cos ^{ - 1}}\left( {\sqrt {\begin{array}{l} {{{\left[ {{M_R}(2,2) + {M_R}(3,3)} \right]}^2}}\\ { + {{\left[ {{M_R}(3,2) - {M_R}(2,3)} \right]}^2}} \end{array}} - 1} \right) \\{R_C} &= {\tan ^{ - 1}}\frac{{{M_R}(3,2) - {M_R}(2,3)}}{{{M_R}(2,2) + {M_R}(3,3)}}, \end{aligned} $$
$$y(t) = x(t) \ast h(t) = \int_{ - \infty }^\infty {x(t - \tau )} h(\tau )d\tau ,$$
$$H(\omega ) = {{Y(\omega )} / {X(\omega )}}.$$
$$\begin{array}{l} W(\omega ) = \frac{{\bar{X}(\omega )}}{{{{|{X(\omega )} |}^2} + 1/SNR}},\\ G(\omega ) = b{e^{ - {\omega ^2}/2\alpha _2^2}} - a{e^{ - {\omega ^2}/2\alpha _1^2}},\;\;\;\;b \ge a,\;{\alpha _2} > {\alpha _1}, \end{array}$$
$$h(t) = {F^{ - 1}}\{{H(\omega )W(\omega )G(\omega )} \},$$
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.