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

Co-registered combined OCT and THz imaging to extract depth and refractive index of a tissue-equivalent test object

Open Access Open Access

Abstract

Terahertz (THz) imaging and optical coherence tomography (OCT) provide complementary information with similar length scales. In addition to OCT’s extensive use in ophthalmology, both methods have shown some promise for other medical applications and non-destructive testing. In this paper, we present an iterative algorithm that combines the information from OCT and THz imaging at two different measurement locations within an object to determine both the depth of the reflecting layers at the two locations and the unknown refractive index of the medium for both the OCT wavelengths and THz frequencies. We validate this algorithm using a silicone test object with embedded layers and show that the depths and refractive index values obtained from the algorithm agreed with the measured values to within 3.3%. We further demonstrate for the first time that OCT and THz images can be co-registered and aligned using unsupervised image registration. Hence we show that a combined OCT/THz system can provide unique information beyond the capability of the separate modalities alone, with possible applications in the medical, industrial and pharmaceutical sectors.

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

1. Introduction

Optical coherence tomography (OCT) and terahertz (THz) imaging both provide the potential for sub-surface characterization of materials and biological tissue that may advance medical, pharmaceutical and industrial applications [13].

OCT is sensitive to the gradient of the refractive index at near-infrared wavelengths which, in biological specimens, depends on scattering from structural arrangements of cells and organelles [4], whereas THz imaging, which uses much longer wavelengths, is much more dependent on absorption mechanisms than scattering [5,6]. THz can detect, for example, variations in hydration [79] while OCT offers clear indications of the structure and morphology of the specimen [3]. Furthermore, THz imaging and OCT both provide 3D information over a similar size, scale and depth; both provide diagnostic information correlated to histopathology [10,11] and both can use the same ultrafast laser source for signal generation. The complementary nature of the two imaging techniques and their common use of ultrafast laser systems poses the intriguing question: could they be combined into a single imaging system?

There are several common applications where the two techniques are currently being researched or applied separately. The penetration depth of THz and OCT, which are both on the scale of millimeters or less in most materials, including biological tissues, tends to guide the nature of these applications. Some of the areas showing potential are near-surface (epithelial) or intra-operative cancer imaging [1215], skin hydration changes [16,17], corneal hydration [18,19], pharmaceutical quality control [20,21] and cultural heritage [2227]. There are possible applications where combining the techniques and overlaying complementary properties of structure and morphology from OCT with functionality and hydration properties from THz imaging, may lead to improved diagnostic sensitivity and specificity. For example, in ophthalmology, this could lead to an improvement in the assessment of certain eye diseases such as Fuch’s dystrophy [28], a degenerative condition of the corneal endothelium which causes swelling of the cornea; or the combined technique could aid clinical decisions regarding corneal grafts [29].

Previous researchers have proposed that imaging a sample using THz and OCT may be of benefit for diagnosis [30]; however, as yet there have been no studies that combined the two. For example, Chernomyrdin et al. [14] used both OCT and THz to examine malignant brain tissue, but did so separately and did not use either OCT to inform the THz properties, or vice versa. They do propose ‘multiplexing the information’ without providing methods to do so. Taylor et al. [31] discuss the potential of using THz and OCT for imaging the cornea to provide structure and hydration, but again, do not discuss the way in which the OCT and THz data can be combined to inform each other, and enhance the information that can be obtained. This study is the first where the information from OCT and THz has been combined to provide enhanced information that is not available using either technique independently.

In this study, we present an algorithm for using the time-of-flight information obtained from OCT and THz images of the same areas as a way to determine the depth of a layer within an object and the refractive index values at THz frequencies and OCT wavelengths. The algorithm uses measurements of OCT and THz at two locations within the object with two different depths to give the depths of the layers and the refractive index values. The algorithm was tested on OCT and THz images of a silicone test object with titanium dioxide (TiO2) scatterers as contrast agents to create reflecting layers at different depths. Silicone has been used previously as a tissue-equivalent test object for OCT [32,33] and we verify its use with THz imaging in this study. The aim was to evaluate the algorithm to determine the accuracy with which it can provide both depth and refractive indices and examine factors that affect its performance.

Further, we demonstrate for the first time that OCT and THz images can be processed to give quantitatively accurate distance information using this algorithm. The images can then be co-registered to overlay the complementary OCT and THz information. Finally, we show that it is possible to determine the depths of layers at other locations in the images using the information obtained from the algorithm, assuming a uniform refractive index throughout the medium. The ability to overlay co-registered THz and OCT images and establish the depths of layers and properties of the materials may support clinical and non-destructive testing assessments for improved pathological decisions, and improved non-destructive characterization of materials.

2. Methods

2.1 Development of the tissue-equivalent test object

A test object was developed that provided tissue-equivalent contrast at optical wavelengths and was an appropriate material in terms of absorption and refractive index for THz imaging. Silicone is frequently used in OCT imaging as a test object due to its low scattering properties and its high transparency in the near-infrared. To test whether silicone was suitable as a THz test object material, we used THz spectroscopy to determine its refractive index (see section 2.1.2). The silicone test object was formed by mixing a silicone base liquid with a hardener according to the manufacturer's recipe (ELASTOSIL RT 601 A/B, Wacker Chemie AG, Germany) [34].

2.1.1 Embedded objects using a contrast agent

Titanium dioxide (TiO2) powder, which is frequently used as a scattering agent at near-infrared wavelengths [35], was tested to determine whether it could create a refractive index boundary that could be detected using THz imaging. After testing a variety of concentrations, it was found that 2% by weight gave absorption and refractive index changes that were detectable with both THz and OCT. The TiO2 powder was added to silicone at 2% concentration to form a layer approximately 1 mm thick. After setting, a square shape and an approximate rectangle shape were cut from the TiO2 layer. The surfaces of both the square and the rectangle were not flat but had shapes that changed across the profile to provide a variety of depths. These shapes were then embedded into the silicone test object before it set at a depth of just over 1 mm. The final dimensions were 26.6 mm × 27.4 mm × 3 mm (Fig. 1(a)).

 figure: Fig. 1.

Fig. 1. (a) Photograph of the front surface of the silicone main test object. The white areas are the embedded highly-scattering silicone-TiO2 contrast objects. Areas A, B and C were the locations where the depth of the silicone back interface was calculated and represent the full thickness of the test object. Areas D and E were the locations where the depth of the embedded square was calculated. The red dotted line represents the profile region along which the depth profile of the rectangle was calculated. (b) An en-face THz image, which takes the absolute value of the integral of the reflected THz pulse over a time range of 12 ps for each time-domain waveform and assigns that value to the pixel to form the 2D image and (c) several of the OCT en-face images manually stitched together. The image property used for each pixel in the OCT image was the maximum of the log intensity of the signal returned from the embedded objects.

Download Full Size | PDF

In a combined system, we would envisage THz imaging with normal incidence [36] or very near normal incidence [37] similar to the OCT set up. However, for this study, the THz imaging system used had a 30° angle of incidence and a larger scanning range of 20 mm, compared to 5 mm for the OCT system. For these reasons, we designed a test object with multiple embedded objects that were similar in size to the field of view of the THz system. For OCT imaging the test object was moved after each image acquisition so that relevant sections could be imaged with OCT that overlapped with the THz image. The regions of the test object that were imaged using both THz and OCT are shown in Fig. 1(a), while Fig. 1(b) and Fig. 1(c) show the THz and OCT en-face images of the test object, respectively.

2.1.2 Refractive index of the silicone test object at THz frequencies

The refractive index of the silicone test object was measured at THz frequencies using a commercial THz time domain system (TeraPulse4000, Teraview Ltd, Cambridge, UK). This system was used for both spectroscopy and imaging of the test object. The system provides a broadband pulse with a range of frequencies 0.06 - 4 THz. Since the refractive index of silicone varies as a function of THz frequency, the measured refractive index was determined using a weighted mean refractive index to represent the bulk properties of silicone across this frequency range. The weighting was based on the spectral amplitude of each frequency component of the pulse, as determined by spectral analysis. The bulk refractive index for the silicone across this frequency range determined, by transmission spectroscopy on a piece of silicone with known thickness, was 1.53 ± 0.08.

2.1.3 Refractive index of the silicone test object at OCT frequencies

OCT measurements were performed using a fiber-based, swept-source polarization-sensitive (PS) OCT system with a standard bench-top configuration. The light source used was a wavelength-swept laser (Axsun Technologies, Bellerica, MA, USA), centered at 1310 nm, with a full sweep range of 130 nm at a 50 kHz sweep rate [38,39]. For comparison of the refractive index obtained for the algorithm for the OCT wavelengths, a calibrated refractive index of the silicone was measured by using a known thickness of several millimeters and measuring the optical path length. The refractive index was found to be 1.43 ± 0.04, which agrees with published values of 1.4 [32]. Small differences may arise based on the formulation of the silicone.

2.2 Imaging the test object

2.2.1 THz imaging of the silicone test object

We used the imaging module of the Terapulse4000 to collect a 3D data set from the test object. The imaging module raster scans a THz pulse over a defined area (x, y) in which each pixel contains the reflected time-domain signal (amplitude and phase). The THz pulse was incident on the sample at 30° and focused into the sample using f2 optics. The THz images were scanned over an area of 20 mm by 20 mm with a sample spacing of 0.2 mm. The front surface of the test object (the surface at which the THz was incident) was placed at the focal point and then the z-position of the test object was focused by maximizing the reflected signal from the embedded square TiO2 object which is approximately 1 mm deep within the silicone. Thus a three-dimensional dataset (x, y, t) is obtained where the third dimension is time which, if the refractive index of the sample is known, can be converted into a spatial dimension. A reference measurement is acquired from a flat gold coated mirror, and is deconvolved from the sample pulse prior to further analysis of the signal. Table 1 compares THz imaging characteristics to OCT.

Tables Icon

Table 1. Comparison of the OCT and THz imaging systems

2.2.2 OCT imaging of the silicone test object

The OCT system had a field of view of 5 mm by 5 mm. The angle of incidence for the OCT system was perpendicular to the surface of the test object. The test object was mapped in the x-y direction in steps of 5 µm. In the z-direction, 1152 pixels were obtained corresponding to an optical path length of 8 mm in air. Due to the smaller field of view for OCT, the test object was imaged at several locations around the boundaries of the embedded square and rectangle (see Fig. 1(c)). The intensity value obtained from processing the OCT images is used for all images and data analysis. The resulting processed images were used to calculate the thickness of the silicone and the depth of the embedded object.

2.3 Image analysis and registration

2.3.1 Determining depth and refractive index values for normal incidence of OCT and THz

Given two locations, 1 and 2, with unknown depths d1 and d2, (as illustrated in Fig. 2) we can calculate d1 and d2, and the unknown refractive indices, nOCT, and nTHz of the material using OCT and THz measurements at these two locations. THz and OCT measurements at each of the locations 1 and 2, provide four known values, tOCT1, tOCT2, tTHz1, tTHz2, where t represents the ‘time-of-flight’ for the time travelled by the electromagnetic radiation for the OCT measurement and the THz measurement to the reflecting surface within the material, assuming normal incidence. In the measurements performed, the THz signal was not normally incident on the silicone, and the method to account for this is given in section 2.3.2. These four values can be used to solve four simultaneous equations that relate to the depth d1 and d2 and the refractive index in the optical regime (nOCT) and the THz range (nTHz). The simultaneous equations were obtained as follows:

$${t_{OCT1}} = {d_1} \ast {n_{OCT}}/c$$
$${t_{OCT2}} = {d_2} \ast {n_{OCT}}/c$$
$${t_{THz1}} = {d_1} \ast {n_{THz}}/c$$
$${t_{THz2}} = {d_2} \ast {n_{THz}}/c$$

 figure: Fig. 2.

Fig. 2. Illustration of the THz pulse incident on the front surface of the silicone at 30° and refracted due to the refractive index nTHz before being reflected from the embedded TiO2 layer at location 1 and silicone back interface at location 2. The angle ϑ2 is unknown in this case since it depends on nTHz, and nTHz is assumed to be unknown. The time-of-flight for a normally incident THz pulse is shown as tTHz, while the measured time-of-flight of the THz pulse, incident at 30° is tTHz30.

Download Full Size | PDF

Rearranging Eq. (1)–(4) gives four simultaneous equations with four unknowns:

$${n_{OCT}} = \; \frac{{c^{\ast} ({{t_{OCT1}} + {t_{OCT2}}} )}}{{{d_1} + {d_2}}}$$
$${n_{THz}} = \; \frac{{c^{\ast} ({{t_{THz1}} + {t_{THz2}}} )}}{{{d_1} + {d_2}}}$$
$${d_1} = \; \frac{{c^{\ast} ({{t_{OCT1}} + {t_{THz1}}} )}}{{{n_{OCT\; }} + {n_{THz}}}}$$
$${d_2} = \; \frac{{c^{\ast} ({{t_{OCT2}} + {t_{THz2}}} )}}{{{n_{OCT}} + {n_{THz}}}}$$
Solving these four simultaneous equations yields the depths d1 and d2 of the two locations, and the refractive index nOCT and nTHz for both OCT wavelengths and THz frequencies. For this study the equations were solved using Matlab’s (Matlab, version 14, Mathworks, Natick USA) fsolve function.

2.3.2 Correction of THz time-of-flight for non-normal THz incidence

Unlike the OCT signal, which is at normal incidence to the surface, the THz pulse is incident at an angle of 30° (ϑ1). The implication of this is that the measured time-of-flight for the THz pulse, tTHz30 must be corrected to obtain the time-of-flight for normally incident THz radiation, tTHz before being input into Eq. (6), (7) and (8). In order to use Snell’s law to determine tTHz at normal incidence, the angle ϑ2 must be known (Fig. 2). However, determining ϑ2 requires knowledge of the refractive index, which we assume is not known in this case. Instead an iterative approach can be used to estimate ϑ2 using an initial first guess for nTHz that then allows tTHz to be calculated from tTHz30 with knowledge of tOCT and an initial guess for nOCT. This updated tTHz is then used in Eq. (6)–(8) to determine the next iteration of depth values and the refractive index values nOCT and nTHz. The updated nTHz and nOCT are used to improve the determination of tTHz from tTHz30. The procedure is repeated iteratively until the values of tTHz, and subsequently d1, d2, nOCT and nTHz converge to a desirable limit. The mathematical relationship between tTHz and tTHz30 is shown in Eq. (9)–(12), and the full algorithm is shown schematically in Fig. 3.

 figure: Fig. 3.

Fig. 3. Schematic flowchart of the iterative algorithm to determine the depth and refractive index values from two locations using OCT and THz images, accounting for the angle of incidence of the THz pulse if it is non-normal. The algorithm assumes initial values for the refractive index at OCT wavelengths and THz frequencies and successively updates these and determines the depth of the layers simultaneously.

Download Full Size | PDF

From Snell’s law we have

$${\vartheta _2} = {\sin ^{ - 1}}{\bigg(}{{{{\sin ({30} )}} {{\bigg /} {\vphantom {{\sin ({30} )} {{n_{THz}}}}}}{{{n_{THz}}}}}} {\bigg)}$$
where nair = 1 and ϑ1 for the THz pulse is 30°.

Using the trigonometric relationships and substituting for ϑ2 we obtain

$${t_{THz}} = {t_{THz30}}^{\ast} \cos {\bigg(}{{{\sin }^{ - 1}}{\bigg(}{{{{\sin ({30} )}} {{\bigg /} {\vphantom {{\sin ({30} )} {{n_{THz}}}}}}{{{n_{THz}}}}}} {\bigg)}} {\bigg)}$$
Given that the depth of a location is the same for THz as it is for OCT, we have
$${n_{THz}} = {n_{THz\; }}^{\ast} {\bigg(}{\; {\raise0.7ex\hbox{${{t_{THz}}}$} \!\mathord{{\bigg /} {\vphantom {{{t_{THz}}} {{t_{OCT}}}}}}\!\lower0.7ex\hbox{${{t_{OCT}}}$}}} {\bigg)}$$
Substituting Eq. (11) into (10), we determine the series for tTHz
$${t_{THz\; }}({i + 1} )= {t_{THz30}}.\cos {\bigg(}{{{\sin }^{ - 1}}{\bigg(}{{\raise0.7ex\hbox{${({\sin ({30} ).{t_{OCT}}} )}$} \!\mathord{{\bigg /} {\vphantom {{({\sin ({30} ).{t_{OCT}}} )} {({{t_{THz}}(i ).{n_{OCT}}} )}}}}\!\lower0.7ex\hbox{${({{t_{THz}}(i ).{n_{OCT}}} )}$}}} {\bigg)}} {\bigg)}$$
Equation (12) can be solved iteratively by inputting an initial first guess for the nOCT and using the iterative previous value of tTHz(n) as well as the known value of tOCT measured from the two locations. The solution tTHz was successively updated until an acceptable convergence was reached, defined arbitrarily as the limit in which the change from the previous value was less than 0.0001 ps. In the study we performed, this was typically after 4-5 iterations. As mentioned previously, the adjusted value for tTHz was then used to solve Eq. (5)–(8) and the procedure repeated, using the updated refractive index for THz and OCT in Eq. (11) and (12) until the refractive index and depth from Eq. (5)–(8) converge. The flow of the full algorithm is given in Fig. 3.

2.3.3 Calculation of layer depths in the silicone using combined THz and OCT

The iterative method of obtaining the depth and refractive index outlined in sections 2.3.1 and 2.3.2 was applied to measurements on the silicone object using five different regions labelled A, B, C, D and E in Fig. 1. Regions A, B and C had no buried layers, so the reflected signal at depth represented the full depth of the silicone and could be compared with a measured depth. Regions D and E had the embedded TiO2 square shape that led to reflections with shorter times-of-flight so it could be used to test the efficacy of the algorithm with a range of depth locations. All regions were 2.5 mm by 2.5 mm in the x-y direction. The actual thickness of the silicone at locations A, B and C was measured using a micrometer. The micrometer thickness measurements were repeated ten times in each of the three locations A, B and C and then averaged to give a mean measured silicone thickness at each location with a standard deviation of the mean.

To measure the depth of the embedded square, the test object was sliced into sections and viewed under a microscope at eight times magnification. By calibration of the microscope image, the distance from the front surface of the silicone to the front surface of the embedded square could be determined at 10 locations within region D and region E. These were averaged and a standard deviation obtained. The mean calculated depths and the refractive index values determined from the algorithm for these ten unique pair combinations are given in Table 2 compared to the measured values for all regions A-E along with a percentage difference between the two.

Tables Icon

Table 2. Values of the depth and refractive index calculated from the iterative algorithm for each of the ten unique pair combinations for the 5 locations A, B, C, D and E. The mean of each value is calculated and presented with the standard deviation. These calculated values are compared to the measured values and the percentage difference between the measured and calculated is given in the bottom row of the table.

2.3.4 Factors affecting the performance of the algorithm

To understand the limitations and performance of the algorithm, we examined several factors that influenced the outcome of the final values. These were the depth of the locations, the refractive index values used for the initial guess, and the number of iterations for convergence. The effect of depth was examined by using a range of five depths, two shallower around 1 mm, and three deeper around 3 mm.

To understand the influence of refractive index initial starting values, we tested the algorithm with a range of starting points of nOCT ranging from 1.3-1.5, and nTHz from 1.4-1.7, which represented a range either side of the expected values. The results for this range of starting refractive index initial estimates are given in Table 3. To observe the speed with which the algorithm converges with each iteration, we present a case study for pair-combination [A, E] where the values at each stage of the iteration are shown in Fig. 4 for 9 iterations compared to a line showing the convergence value.

 figure: Fig. 4.

Fig. 4. Plots showing the convergence of the algorithm over 9 iterations for (a) the THz time-of-flight for region A (circles) and region E (triangles). A full line is plotted at tTHz = 15.665 ps and dotted line at tTHz = 6.647 ps showing the convergence value. (b) The depth determined from each iteration of the algorithm for region A (circles) and region E (triangles) with the full line at d1 = 3.116 mm, and dotted line at d2 = 1.296 mm. (c) The refractive index from the algorithm after each iteration for nTHz (circles) and nOCT (triangles), with a full line at nTHz = 1.517 and nOCT = 1.418 showing convergence.

Download Full Size | PDF

Tables Icon

Table 3. The impact of the initial starting values of the refractive index for OCT (nOCT) and THz (nTHz) on the outcome of the algorithm for the case example of pair combination [A,E]. A percentage difference was calculated for each value compared to the measured value. Note that the measured values for d1 and d2 are 3.05 mm and 1.27 mm respectively for regions A and E. The measured refractive index for nOCT is 1.43, and measured refractive index for THz, nTHz, is 1.53.

2.3.5 Calibration and co-registration of the OCT and THz image information

The outcome from the algorithm is a set of four pieces of information: d1 and d2, nOCT, and nTHz. By assuming the medium is uniform prior to the embedded layers, the OCT and THz images can then be scaled according to these refractive index values. These ‘calibrated’ images provide the possibility of overlaying the images by co-registration since they are on equivalent distance scales. To test this, 2D cross-sectional OCT and THz images, known as B-scan images, were co-registered using Matlab’s imregconfig and imregister commands. The ‘multimodal’ option was selected, which maximizes the mutual information between the two images, ignoring the incommensurate intensities of the two modalities, using a one-plus-one evolutionary optimization configuration. Rigid registration was selected since the scale was determined and calibrated by the refractive index. The co-registered OCT and THz B-scan images are shown in Fig. 5 and Fig. 6 together with the photograph of the arrangement and the individual THz and OCT B-scans.

 figure: Fig. 5.

Fig. 5. (a) Side-view photograph of the test object showing the embedded square silicone-TiO2 test object. THz and OCT were incident on the top surface. (b) B-scan from the OCT image of the test object, showing a scattered signal from the top surface of the test object and increased scatter from the embedded square. (c) B-scan THz image showing reflections from the top surface of the test object, the top surface of the embedded square, and the bottom surface of the test object. (d) The co-registered overlay of the THz and OCT images showing the agreement between the two imaging methods, where the OCT image is in ‘hot’ colourmap for contrast.

Download Full Size | PDF

 figure: Fig. 6.

Fig. 6. (a) Side-view photograph of the test object showing the embedded rectangular test object. THz and OCT were incident on the top surface. (b) OCT B-scan image, along the profile shown by the red dotted line in Fig. 1, showing a reflection from the top surface of the test object and increased scatter from the embedded rectangle. (c) B-scan THz image showing reflections from the top surface of the test object, the top surface of the embedded rectangle, and the bottom surface of the test object. (d) The co-registered overlay of the THz and OCT B-scan images showing the agreement between the two imaging methods, where the OCT image is in ‘hot’ colourmap for contrast.

Download Full Size | PDF

As a further step to test the co-registration of the images, the depth of the embedded rectangle was determined over a profile that was 1 mm wide along the direction shown in Fig. 1(a) by the red dotted line. The depths extracted from the OCT and THz images along this profile were combined to give a depth profile, shown in Fig. 7.

 figure: Fig. 7.

Fig. 7. (a) Graph showing the depth of the embedded rectangle along the profile marked in Fig. 1 with a red dotted line. The depth of the layer identified from the THz image is shown as blue triangles, while the depth obtained with OCT images is shown as red circles. ‘OCT 1’ represents the data from the left hand OCT image in Fig. 6(b), while ‘OCT 2’ represents the data from the right hand OCT image

Download Full Size | PDF

3. Results

The results for the algorithm are presented in Tables 2, 3 and Fig. 4. Table 2 gives the values output by the algorithm for the range of pair locations used in this study compared to the appropriate measured value. Table 3 displays the manner in which the initial starting values for the refractive index of OCT and THz affect the algorithm. Figure 4 demonstrates the convergence of the algorithm over 9 recorded iterations, for the example, pair combination of [A, E], showing how the values converge in under 6 steps.

Figures 5 and 6 show the B-scans for the THz and OCT images of regions in the embedded square and rectangle. The OCT images co-registered by automated image registration are overlaid on the THz images to show the agreement in the two imaging modalities in Fig. 5(d) and Fig. 6(d).

Given the way the phantom was made, the rectangular inclusion was not parallel with the surface. We took advantage of this to determine how accurately the depth of features that were not used to calibrate the algorithm could be determined. The depths were first calibrated using measurements taken from areas A-E as described above, and the results of the algorithm were then applied to the varying depth of the top surface of the rectangle. These depth profiles, taken from the THz and two OCT images are compared in Fig. 7 and show good agreement over depths ranging from 1.3–1.8 mm.

4. Discussion

Overall the algorithm to extract depths and refractive index values using OCT and THz imaging measurements at two unknown depth locations performed well. The mean depth of the silicone test object and depth of the embedded layers calculated from the algorithm agreed within 3.3% of the measured depths. The refractive index values agreed even more closely, within 1.5%. Interestingly the depths from the algorithm were slightly overestimated compared to the measured depth. This corresponded with the refractive index calculated from the algorithm being slightly lower than the measured values.

It might be expected that the wider the depth range, the better the results. However, this was not observed with the analysis. For example, the algorithm gave the same depth for region C when the pair combination of [A, C] and [C, D] were used. The difference in measured depth between regions A and C is 30 µm, while the difference in depth between regions C and D was measured as 1.82 mm. This robustness to the depth of the two locations is promising for the technique to be used with samples that only exhibit a small but detectable variation of depths, for example in skin layers.

The refractive index obtained from the algorithm for OCT was 1.422, which compared favorably to the measured value of 1.43. The difference was just 0.6%. Similarly, for THz frequencies, the mean refractive index obtained from all the pair combinations was 1.507 compared to the value of 1.53 determined from the THz spectroscopy measurements. The difference, in this case, was 1.5%.

The results were almost invariant to the initial nTHz value, as can be seen from Table 2, but did depend on the initial value of nOCT. This may be because nOCT is used in both the determination of the corrected tTHz as well as the solution of the four simultaneous equations. The outcome for depth and refractive index was constant if nTHz varied from 1.5 - 1.7. However, for nOCT ranging from 1.3 - 1.5, the calculated values differed from the measured ones by between 7.7% and 0.2% for the example pair [A, E].

The two examples of co-registration shown in Fig. 5(d) and Fig. 6(d) demonstrates very good visual agreement between the OCT and THz B-scans images. This indicates the scaling of images for the refractive index calculated from the algorithm works well, assuming the medium is uniform, to establish the images on equivalent distance scales that then enables the overlay of the images. This similarity is supported by the fact the images can be co-registered using an automated, unsupervised, rigid registration algorithm, after this scaling.

As a further demonstration of the ability of the algorithm to improve the characterization of the samples, Fig. 7 clearly shows that the profile of the rectangle agrees closely for both the OCT and THz images. This analysis shows that despite even with THz incident on the sample at 30°, by accounting for the known angle of incidence the algorithm still operates successfully.

Using this algorithm, it is possible to combine OCT and THz images to overlay the complementary information they provide. OCT in the past has assumed to be able to extract thickness, with the refractive index assumed from the literature or measured for the material. Potentially this algorithm could be used with OCT alone at two different wavelengths; however, the advantage of combining with THz is that the hydration properties may be obtained that informs the functionality not just structure, and this may be of value to clinicians. Using this algorithm, it may be possible to confirm these thickness measurements. This may be of value, for example, in corneal applications, confirming the thickness of the cornea as well as providing enhanced information from the THz that is sensitive to pathology. This may also advance the combined OCT/THz modality in applications for the skin and breast cancer margin detection, as well as in industry, for non-destructive testing and quality control purposes. The next step to achieving this would be examining how the algorithm operates with multiple layered structures and with layered gradients such as that found in the skin.

5. Conclusion

OCT and THz give different and complementary information that may improve diagnostic and characterization capabilities beyond just using one of the modalities alone. We showed through the use of a silicone test object with TiO2-silicone embedded layers that the THz and OCT images may be used to extract structural information of the object, as well as a refractive index using an iterative algorithm. This method overcomes the challenges of unknown refractive index values and the impact these have on the uncertainty in calculating layer depth.

The algorithm we presented solves four simultaneous equations for four unknowns using OCT and THz measurements at two locations. The outcome of the algorithm is the depth of the two locations as well as the refractive index for the OCT wavelength and THz frequencies. Quantitatively, the thicknesses and depths determined from the image agreed very closely to within 3.3% with the measured values. Refractive index values obtained from the algorithm agreed to within 1.5% of the measured values. The refractive index calculation was more dependent on the initial estimate of the refractive index of OCT wavelengths than THz possibly because of the more significant role this plays in the algorithm.

We presented overlaid THz and OCT images based on these refractive index values that aligned well and showed, for the first time, that THz and OCT images could be automatically scaled and co-registered with unsupervised methods, even without knowing the refractive index in either regime. This allows the complementary information from the OCT structure and THz hydration sensitivity to be overlaid to potentially enhance diagnostic capabilities in applications where the two show promise.

Funding

Australian Research Council (DP1096514, DP150100635, FT180100683); National Health and Medical Research Council (1141552); National Research Foundation of Korea funded by the Ministry of Science, ICT and Future PlanningMinistry of Science and ICT (2019H1D3A2A02101784).

Acknowledgements

We thank Qingyun Li for advising on the operation of OCT, sharing processing code and discussing general OCT questions.

This research was supported by the Australian Government through the Australian Research Council's Discovery Projects funding scheme (DP1096514, DP150100635). VPW is supported by the Australian Research Council via a Future Fellowship (project number FT180100683) funded by the Australian Government. MH was funded through NHMRC grant 1141552, http://purl.org/au-research/grants/nhmrc/1141552. BC acknowledges funding through the Brain Pool Program from the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (2019H1D3A2A02101784).

Disclosures

This work is the subject of an Australian Provisional Patent Application (No. 2019903145).

References

1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef]  

2. D. Grischkowsky, S. Keiding, M. van Exter, and C. Fattinger, “Far-infrared time-domain spectroscopy with terahertz beams of dielectrics and semiconductors,” J. Opt. Soc. Am. B 7(10), 2006–2015 (1990). [CrossRef]  

3. M. E. Brezinski, G. J. Tearney, B. E. Bouma, J. A. Izatt, M. R. Hee, E. A. Swanson, J. F. Southern, and J. G. Fujimoto, “Optical Coherence Tomography for Optical Biopsy,” Circulation 93(6), 1206–1213 (1996). [CrossRef]  

4. Y. Pan, R. Birngruber, and R. Engelhardt, “Contrast limits of coherence-gated imaging in scattering media,” Appl. Opt. 36(13), 2979–2983 (1997). [CrossRef]  

5. J. T. Kindt and C. A. Schmuttenmaer, “Far-Infrared Dielectric Properties of Polar Liquids Probed by Femtosecond Terahertz Pulse Spectroscopy,” J. Phys. Chem. 100(24), 10373–10379 (1996). [CrossRef]  

6. E. Pickwell and V. P. Wallace, “Biomedical applications of terahertz technology,” J. Phys. D: Appl. Phys. 39(17), R301–R310 (2006). [CrossRef]  

7. M. Heyden, S. Ebbinghaus, and M Havenith, “Terahertz Spectroscopy as a Tool to Study Hydration Dynamics,” in Encyclopedia of Analytical Chemistry, R. A. Meyers, ed. (Wiley, 2010).

8. Z. D. Taylor, R. S. Singh, D. B. Bennett, P. Tewari, C. P. Kealey, N. Bajwa, M. O. Culjat, A. Stojadinovic, H. Lee, J. Hubschman, E. R. Brown, and W. S. Grundfest, “THz Medical Imaging: in vivo Hydration Sensing,” IEEE Trans. Terahertz Sci. Technol. 1(1), 201–219 (2011). [CrossRef]  

9. Q. Sun, E. P. J. Parrott, Y. He, and E. Pickwell-MacPherson, “In vivo THz imaging of human skin: Accounting for occlusion effects,” J. Biophotonics 11(2), e201700111 (2018). [CrossRef]  

10. E. Berry, J. W. Handley, A. J. Fitzgerald, W. J. Merchant, R. D. Boyle, N. N. Zinov’ev, R. E. Miles, J. M. Chamberlain, and M. A. Smith, “Multispectral classification techniques for terahertz pulsed imaging: an example in histopathology,” Med. Eng. Phys. 26(5), 423–430 (2004). [CrossRef]  

11. A. J. Fitzgerald, V. P. Wallace, S. E. Pinder, A. D. Purushotham, P. O’Kelly, and P. C. Ashworth, “Classification of terahertz-pulsed imaging data from excised breast tissue,” J. Biomed. Opt. 17(1), 016005 (2012). [CrossRef]  

12. A. J. Fitzgerald, V. P. Wallace, M. Jimenez-Linan, L. Bobrow, R. J. Pye, A. D. Purushotham, and D. D. Arnone, “Terahertz Pulsed Imaging of Human Breast Tumors,” Radiology 239(2), 533–540 (2006). [CrossRef]  

13. 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]  

14. N. V. Chernomyrdin, I. Dolganova, S.-I. Beshplav, G. Komandin, I. V. Reshetov, A. A. Potapov, V. Tuchin, K. Zaytsev, P. V. Alexandrova, K. Malakhov, P. V. Nikitin, A. V. Kosyr’kova, and G. Musina, Differentiation of healthy and malignant brain tissues using terahertz pulsed spectroscopy and optical coherence tomography (2019), p. 5.

15. B. J. Vakoc, D. Fukumura, R. K. Jain, and B. E. Bouma, “Cancer imaging by optical coherence tomography: preclinical progress and clinical potential,” Nat. Rev. Cancer 12(5), 363–368 (2012). [CrossRef]  

16. J. Wang, Q. Sun, R. I. Stantchev, T.-W. Chiu, A. T. Ahuja, and E. Pickwell-MacPherson, “In vivo terahertz imaging to evaluate scar treatment strategies: silicone gel sheeting,” Biomed. Opt. Express 10(7), 3584–3590 (2019). [CrossRef]  

17. J. M. Crowther, A. Sieg, P. Blenkiron, C. Marcott, P. J. Matts, J. R. Kaczvinsky, and A. V. Rawlings, “Measuring the effects of topical moisturizers on changes in stratum corneum thickness, water gradients and hydration in vivo,” Br. J. Dermatol. 159, 567–577 (2008). [CrossRef]  

18. D. B. Bennett, Z. D. Taylor, P. Tewari, R. S. Singh, M. O. Culjat, W. S. Grundfest, D. J. Sassoon, R. D. Johnson, J.-P. Hubschman, and E. R. Brown, “Terahertz sensing in corneal tissues,” J. Biomed. Opt. 16(5), 057003 (2011). [CrossRef]  

19. M. Pircher, E. Götzinger, R. Leitgeb, A. F. Fercher, and C. K. Hitzenberger, “Measurement and imaging of water concentration in human cornea with differential absorption optical coherence tomography,” Opt. Express 11(18), 2190–2197 (2003). [CrossRef]  

20. A. J. Fitzgerald, B. E. Cole, and P. F. Taday, “Nondestructive analysis of tablet coating thicknesses using terahertz pulsed imaging,” J. Pharm. Sci. 94(1), 177–183 (2005). [CrossRef]  

21. H. Lin, Y. Dong, D. Markl, B. M. Williams, Y. Zheng, Y. Shen, and J. A. Zeitler, “Measurement of the Intertablet Coating Uniformity of a Pharmaceutical Pan Coating Process With Combined Terahertz and Optical Coherence Tomography In-Line Sensing,” J. Pharm. Sci. 106(4), 1075–1084 (2017). [CrossRef]  

22. A. J. L. Adam, P. C. M. Planken, S. Meloni, and J. Dik, “TeraHertz imaging of hidden paint layers on canvas,” Opt. Express 17(5), 3407–3416 (2009). [CrossRef]  

23. C. Seco-Martorell, V. Lopez-Dominguez, G. Arauz-Garofalo, A. Redo-Sanchez, J. Palacios, and J. Tejada, Goya’s artwork imaging with Terahertz waves (2013), Vol. 21, pp. 17800–17805.

24. C. S. Cheung, M. Spring, and H. Liang, “Ultra-high resolution Fourier domain optical coherence tomography for old master paintings,” Opt. Express 23(8), 10145–10157 (2015). [CrossRef]  

25. M. Iwanicka, E. A. Kwiatkowska, M. Sylwestrzak, and P. Targowski, Application of optical coherence tomography (OCT) for real time monitoring of consolidation of the paint layer in Hinterglasmalerei objects, SPIE Optical Metrology (SPIE, 2011), Vol. 8084.

26. P. Targowski, B. Rouba, M. Góra, L. Tymińska-Widmer, J. Marczak, and A. Kowalczyk, “Optical coherence tomography in art diagnostics and restoration,” Appl. Phys. A 92(1), 1–9 (2008). [CrossRef]  

27. A. Gibson, K. E. Piquette, U. Bergmann, W. Christens-Barry, G. Davis, M. Endrizzi, S. Fan, S. Farsiu, A. Fitzgerald, J. Griffiths, C. Jones, G. Li, P. L. Manning, C. Maughan Jones, R. Mazza, D. Mills, P. Modregger, P. R. T. Munro, A. Olivo, A. Stevenson, B. Venugopal, V. Wallace, R. A. Wogelius, M. B. Toth, and M. Terras, “An assessment of multimodal imaging of subsurface text in mummy cartonnage using surrogate papyrus phantoms,” Heritage Sci. 6(1), 7 (2018). [CrossRef]  

28. A. O. Eghrari and J. D. Gottsch, “Fuchs’ corneal dystrophy,” Expert Rev. Ophthalmol. 5(2), 147–159 (2010). [CrossRef]  

29. J. L. Güell, M. A. El Husseiny, F. Manero, O. Gris, and D. Elies, “Historical Review and Update of Surgical Treatment for Corneal Endothelial Diseases,” Ophthalmol. Ther. 3(1-2), 1–15 (2014). [CrossRef]  

30. R. M. Woodward, V. P. Wallace, R. J. Pye, B. E. Cole, D. D. Arnone, E. H. Linfield, and M. Pepper, “Terahertz pulse imaging of ex vivo basal cell carcinoma,” J. Invest. Dermatol. 120(1), 72–78 (2003). [CrossRef]  

31. Z. D. Taylor, J. Garritano, S. Sung, N. Bajwa, D. B. Bennett, B. Nowroozi, P. Tewari, J. W. Sayre, J. P. Hubschman, S. X. Deng, E. R. Brown, and W. S. Grundfest, “THz and mm-Wave Sensing of Corneal Tissue Water Content: In Vivo Sensing and Imaging Results,” IEEE Trans. Terahertz Sci. Technol. 5(2), 184–196 (2015). [CrossRef]  

32. G. Lamouche, B. F. Kennedy, K. M. Kennedy, C.-E. Bisaillon, A. Curatolo, G. Campbell, V. Pazos, and D. D. Sampson, “Review of tissue simulating phantoms with controllable optical, mechanical and structural properties for use in optical coherence tomography,” Biomed. Opt. Express 3(6), 1381–1398 (2012). [CrossRef]  

33. D. M. M. D. Bruin, R. H. Bremmer, V. M. Kodach, R. D. Kinkelder, J. V. Marle, T. G. V. Leeuwen, and D. J. Faber, “Optical phantoms of varying geometry based on thin building blocks with controlled optical properties,” J. Biomed. Opt. 15(01), 1–10 (2010). [CrossRef]  

34. B. F. Kennedy, T. R. Hillman, R. A. McLaughlin, B. C. Quirk, and D. D. Sampson, “In vivo dynamic optical coherence elastography using a ring actuator,” Opt. Express 17(24), 21762–21772 (2009). [CrossRef]  

35. A. Curatolo, B. F. Kennedy, and D. D. Sampson, “Structured three-dimensional optical phantom for optical coherence tomography,” Opt. Express 19(20), 19480–19485 (2011). [CrossRef]  

36. F. D. J. Brunner, A. Schneider, and P. Günter, “A terahertz time-domain spectrometer for simultaneous transmission and reflection measurements at normal incidence,” Opt. Express 17(23), 20684–20693 (2009). [CrossRef]  

37. Q. Cassar, A. Al-Ibadi, L. Mavarani, P. Hillger, J. Grzyb, G. MacGrogan, T. Zimmer, U. R. Pfeiffer, J.-P. Guillet, and P. Mounaix, “Pilot study of freshly excised breast tissue response in the 300–600 GHz range,” Biomed. Opt. Express 9(7), 2930–2942 (2018). [CrossRef]  

38. Q. Li, K. Karnowski, P. B. Noble, A. Cairncross, A. James, M. Villiger, and D. D. Sampson, “Robust reconstruction of local optic axis orientation with fiber-based polarization-sensitive optical coherence tomography,” Biomed. Opt. Express 9(11), 5437–5455 (2018). [CrossRef]  

39. J. Walther, Q. Li, M. Villiger, C. S. Farah, E. Koch, K. Karnowski, and D. D. Sampson, “Depth-resolved birefringence imaging of collagen fiber organization in the human oral mucosa in vivo,” Biomed. Opt. Express 10(4), 1942–1956 (2019). [CrossRef]  

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 (7)

Fig. 1.
Fig. 1. (a) Photograph of the front surface of the silicone main test object. The white areas are the embedded highly-scattering silicone-TiO2 contrast objects. Areas A, B and C were the locations where the depth of the silicone back interface was calculated and represent the full thickness of the test object. Areas D and E were the locations where the depth of the embedded square was calculated. The red dotted line represents the profile region along which the depth profile of the rectangle was calculated. (b) An en-face THz image, which takes the absolute value of the integral of the reflected THz pulse over a time range of 12 ps for each time-domain waveform and assigns that value to the pixel to form the 2D image and (c) several of the OCT en-face images manually stitched together. The image property used for each pixel in the OCT image was the maximum of the log intensity of the signal returned from the embedded objects.
Fig. 2.
Fig. 2. Illustration of the THz pulse incident on the front surface of the silicone at 30° and refracted due to the refractive index nTHz before being reflected from the embedded TiO2 layer at location 1 and silicone back interface at location 2. The angle ϑ2 is unknown in this case since it depends on nTHz, and nTHz is assumed to be unknown. The time-of-flight for a normally incident THz pulse is shown as tTHz, while the measured time-of-flight of the THz pulse, incident at 30° is tTHz30.
Fig. 3.
Fig. 3. Schematic flowchart of the iterative algorithm to determine the depth and refractive index values from two locations using OCT and THz images, accounting for the angle of incidence of the THz pulse if it is non-normal. The algorithm assumes initial values for the refractive index at OCT wavelengths and THz frequencies and successively updates these and determines the depth of the layers simultaneously.
Fig. 4.
Fig. 4. Plots showing the convergence of the algorithm over 9 iterations for (a) the THz time-of-flight for region A (circles) and region E (triangles). A full line is plotted at tTHz = 15.665 ps and dotted line at tTHz = 6.647 ps showing the convergence value. (b) The depth determined from each iteration of the algorithm for region A (circles) and region E (triangles) with the full line at d1 = 3.116 mm, and dotted line at d2 = 1.296 mm. (c) The refractive index from the algorithm after each iteration for nTHz (circles) and nOCT (triangles), with a full line at nTHz = 1.517 and nOCT = 1.418 showing convergence.
Fig. 5.
Fig. 5. (a) Side-view photograph of the test object showing the embedded square silicone-TiO2 test object. THz and OCT were incident on the top surface. (b) B-scan from the OCT image of the test object, showing a scattered signal from the top surface of the test object and increased scatter from the embedded square. (c) B-scan THz image showing reflections from the top surface of the test object, the top surface of the embedded square, and the bottom surface of the test object. (d) The co-registered overlay of the THz and OCT images showing the agreement between the two imaging methods, where the OCT image is in ‘hot’ colourmap for contrast.
Fig. 6.
Fig. 6. (a) Side-view photograph of the test object showing the embedded rectangular test object. THz and OCT were incident on the top surface. (b) OCT B-scan image, along the profile shown by the red dotted line in Fig. 1, showing a reflection from the top surface of the test object and increased scatter from the embedded rectangle. (c) B-scan THz image showing reflections from the top surface of the test object, the top surface of the embedded rectangle, and the bottom surface of the test object. (d) The co-registered overlay of the THz and OCT B-scan images showing the agreement between the two imaging methods, where the OCT image is in ‘hot’ colourmap for contrast.
Fig. 7.
Fig. 7. (a) Graph showing the depth of the embedded rectangle along the profile marked in Fig. 1 with a red dotted line. The depth of the layer identified from the THz image is shown as blue triangles, while the depth obtained with OCT images is shown as red circles. ‘OCT 1’ represents the data from the left hand OCT image in Fig. 6(b), while ‘OCT 2’ represents the data from the right hand OCT image

Tables (3)

Tables Icon

Table 1. Comparison of the OCT and THz imaging systems

Tables Icon

Table 2. Values of the depth and refractive index calculated from the iterative algorithm for each of the ten unique pair combinations for the 5 locations A, B, C, D and E. The mean of each value is calculated and presented with the standard deviation. These calculated values are compared to the measured values and the percentage difference between the measured and calculated is given in the bottom row of the table.

Tables Icon

Table 3. The impact of the initial starting values of the refractive index for OCT (nOCT) and THz (nTHz) on the outcome of the algorithm for the case example of pair combination [A,E]. A percentage difference was calculated for each value compared to the measured value. Note that the measured values for d1 and d2 are 3.05 mm and 1.27 mm respectively for regions A and E. The measured refractive index for nOCT is 1.43, and measured refractive index for THz, nTHz, is 1.53.

Equations (12)

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

t O C T 1 = d 1 n O C T / c
t O C T 2 = d 2 n O C T / c
t T H z 1 = d 1 n T H z / c
t T H z 2 = d 2 n T H z / c
n O C T = c ( t O C T 1 + t O C T 2 ) d 1 + d 2
n T H z = c ( t T H z 1 + t T H z 2 ) d 1 + d 2
d 1 = c ( t O C T 1 + t T H z 1 ) n O C T + n T H z
d 2 = c ( t O C T 2 + t T H z 2 ) n O C T + n T H z
ϑ 2 = sin 1 ( sin ( 30 ) / sin ( 30 ) n T H z n T H z )
t T H z = t T H z 30 cos ( sin 1 ( sin ( 30 ) / sin ( 30 ) n T H z n T H z ) )
n T H z = n T H z ( t T H z / t T H z t O C T t O C T )
t T H z ( i + 1 ) = t T H z 30 . cos ( sin 1 ( ( sin ( 30 ) . t O C T ) / ( sin ( 30 ) . t O C T ) ( t T H z ( i ) . n O C T ) ( t T H z ( i ) . n O C T ) ) )
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.