Terahertz in-line digital holography of human hepatocellular carcinoma tissue

Terahertz waves provide a better contrast in imaging soft biomedical tissues than X-rays, and unlike X-rays, they cause no ionisation damage, making them a good option for biomedical imaging. Terahertz absorption imaging has conventionally been used for cancer diagnosis. However, the absorption properties of a cancerous sample are influenced by two opposing factors: an increase in absorption due to a higher degree of hydration and a decrease in absorption due to structural changes. It is therefore difficult to diagnose cancer from an absorption image. Phase imaging can thus be critical for diagnostics. We demonstrate imaging of the absorption and phase-shift distributions of 3.2 mm × 2.3 mm × 30-μm-thick human hepatocellular carcinoma tissue by continuous-wave terahertz digital in-line holography. The acquisition time of a few seconds for a single in-line hologram is much shorter than that of other terahertz diagnostic techniques, and future detectors will allow acquisition of meaningful holograms without sample dehydration. The resolution of the reconstructions was enhanced by sub-pixel shifting and extrapolation. Another advantage of this technique is its relaxed minimal sample size limitation. The fibrosis indicated in the phase distribution demonstrates the potential of terahertz holographic imaging to obtain a more objective, early diagnosis of cancer.

O n the electromagnetic spectrum, terahertz waves are between infrared radiation and microwaves, and they share some of the properties of each of these spectra. Similar to infrared radiation, terahertz waves are sensitive to infrared-active modes and provide a better contrast in the imaging of soft biomedical tissues than do X-rays. Similar to microwave radiation, and in contrast to X-rays, terahertz radiation can penetrate non-conducting materials, such as clothes, without causing ionising radiation damage. These properties make terahertz waves suitable for non-invasive biomedical imaging and cancer diagnosis [1][2][3][4][5][6][7][8][9][10][11][12][13] .
For most cancers, confident diagnosis at an early stage is a challenging task, and several cross-checking tests are typically performed. Diseased or cancerous tissues contain more interstitial water than healthy tissues as a result of oedema or increased vascularity 3 , and the analysis of water content can be performed by measuring the absorption of terahertz radiation [1][2][3][4][5][6][7][8][9][10][11][12] . However, the absorption contrast is affected not only by the water content but also by the structural changes caused by cancer 4 . Intentionally dehydrated cancerous tissue or massively inflamed liver tissue has a lower absorption than normal tissue at certain terahertz frequencies 5,6 . Thus, there are two competing influences on the absorption of cancerous tissues: an increase in absorption due to higher water content 3 and a decrease in absorption due to the altered cell density and protein content 2 . It is therefore difficult to evaluate the degree of cancerous damage from absorption imaging alone, and phase-shift imaging can be critical for diagnostics.
Terahertz time domain spectroscopy (THz-TDS) 14,15 is capable of extracting both the amplitude and phase information of the pulsed terahertz signal, enabling the absorption coefficient and the refractive index of the object to be obtained simultaneously [3][4][5] . To further determine the histological features of the tissue, the image is mapped out pixel by pixel by raster scanning the object based on THz-TDS in the reflective mode [6][7][8][9] , also referred to as terahertz pulsed imaging (TPI) 16 . However, TPI has the following disadvantages: the sample size must be larger than 5 mm 3 5 mm to avoid diffraction effects at the sample's edges, and dehydration must be prevented during the slow acquisition time, which varies from 10 min for a 26 mm 3 20 mm sample 6 to 6 hours for an 18.7 mm 3 18.7 mm 9 sample depending on the hardware and the scanning step. Continuous-wave (CW) terahertz imaging of cancerous tissues has become more popular due to the ready availability of commercial high-power terahertz laser systems. However, the acquisition time of the CW terahertz fibre-scanning near-field imaging technique is even longer than that of TPI and is approximately 50 min for 10 mm 3 10 mm samples 10-12 . Holography 17 , by contrast, allows recording of both the amplitude and the phase distribution of the scattered waves in a single holographic record, which can then be used to retrieve the absorption and phase-shifting properties of the object by numerical reconstruction 18 . Terahertz digital holography 13,[19][20][21][22][23] has been applied to achieve the simultaneous reconstruction of absorption and phase-shift distributions 13,21 , and a resolution of approximately 0.2 mm was achieved 23 , although this technique has not been applied to image cancer tissues.

Results
Experimental setup. In the present study, we imaged cancerous tissues using in-line holography. A Gabor-type in-line holographic scheme 17 does not require optical elements between the sample and the detector, which allows aberrations to be avoided and makes this approach a good choice for terahertz imaging 13,23 . An experimental set-up for terahertz in-line digital holography is depicted in Fig. 1. A far-infrared radiation laser system pumped by a CO 2 laser was used to provide a 2.52 THz CW beam with a wavelength of 118.83 mm. The in-line hologram was formed by interference between the wave scattered by the object and the unscattered beam, and it was recorded by a pyroelectric array detector with 124 3 124 pixels, a pixel size of 85 mm 3 85 mm, and a pitch of 100 mm 3 100 mm. The conditions for realising holography have been outlined by Dennis Gabor 24 : for the successful reconstruction of an object from its in-line hologram, the area occupied by the object must be approximately 1% of the entire area illuminated by the incident wave. Thus, in in-line holography, the smaller the object size, the more reference wave information there is in the hologram and the more reliable is the object reconstruction. This condition relaxes the requirement for the minimal size of the sample. However, the sample must not be smaller than the wavelength of the probing waves.
As a calibration object to test the absorption and phase imaging of our terahertz in-line holographic setup, we selected a metallic object -a steel screwdriver blade, a photo of which is shown in Fig. 2(a). The sample was placed at a distance of 19 mm from the detector. To enable the simultaneous reconstruction of the absorption and phase distributions from a single hologram, the hologram must be normalised with respect to the background intensity 25,26 . Thus, images with and without the sample were recorded. To increase the sampling of the recorded images, the detector was laterally shifted in a stepwise fashion (as indicated in Fig. 1) with sub-pixel-size steps for a total of 25 steps. At each detector position, 1,000 frames of both the hologram and background images were recorded and averaged (see Methods: Statistical averaging of individual frames). The sampling required 1/48 sec per single frame and 20.83 sec per 1,000 frames; thus, the total acquisition time for all 25 steps was approximately 9 min. The resulting 25 averaged holograms and 25 background images were each upsampled from 124 3 124 pixels to 992 3 992 pixels and aligned using a sub-pixel registration method 27 . The normalised hologram (see Methods: Hologram normalisation) is shown in Fig. 2(b). The reconstructed absorption and phase distributions are shown in Fig. 2(c)-(f) (see Methods: Hologram reconstruction). The metal blade of the screwdriver is relatively thin and exhibits finite absorption of terahertz waves, and some phase contrast is also observed. Moreover, after applying an iterative reconstruction routine 25 (see Methods: Hologram reconstruction), the twin image in the reconstructions is suppressed and the pronounced fringe structure becomes apparent in the phase distribution.
Cancerous liver sample. A cancerous liver sample was obtained from a patient who underwent a liver resection. The section of liver was diagnosed as containing a hepatocellular carcinoma (HCC) tumour by histology. The tissue was stored immediately after resection at 280uC. While still frozen, it was cut into 30-mmthick slices with a microtome, and the slices were placed in pHbalanced saline (pH 7.4), where they were kept for two days prior to the holographic experiments. For the holographic imaging, a patch of approximately 3.2 mm 3 2.3 mm was excised from one of the slices and mounted on a 3-mm-thick quartz slide. The slide was placed approximately 25 mm in front of the detector and tilted at an angle of approximately 12u relative to the incident beam to  Hologram acquisition and reconstruction. A digital photo of the sample is shown in Fig. 3(a). A single normalised hologram 124 3 124 pixels in size, resulting from the averaging of 1,000 frames recorded at a selected detector position, is shown in Fig. 3(b). The reconstructed absorption and phase-shift distributions from the normalised hologram are shown in Fig. 3(c) and Fig. 3(d), respectively. Thus, a single hologram acquired during a 20.83-sec interval is sufficient to arrive at low-resolution reconstructions of the absorption and phase-shifting properties of the sample. To retrieve these distributions at a higher resolution, the holograms acquired at different detector positions were upsampled, aligned and averaged in the same manner as described above for the test sample. The reconstructed absorption and phase-shift distributions from the resulting normalised hologram of 992 3 992 pixels size are shown in Fig. 4(b) and Fig. 4 (c). Whereas the reconstructed absorption of the sample is quite low, which can be explained by the dehydration during data acquisition, the reconstructed phaseshift distribution demonstrates good contrast, matches well with the contour of the tissue sample, and exhibits detailed features not observed in the absorption image. For example, the hole in the top left corner of the sample (indicated by a green arrow in Fig. 3(a)), which could be either a cut across a vessel or the result of damage caused by the freezing procedure, is pronounced in the phase reconstruction shown in Fig. 4 (c).
The small numerical aperture (N.A. 5 0.24) of the set-up leads to a limited resolution in the reconstructed images 28 . It is, however, possible to increase the N.A. and therefore the resolution by a numerical extrapolation of the experimental hologram 29,30 . The hologram was extrapolated by an iterative reconstruction based on the error-reduction algorithm 31 , wherein the propagation between the detector and the object planes was calculated using the angular spectrum   approach 32 (see Methods: Extrapolation procedure). The resulting extrapolated hologram is shown in Fig. 4(f), where the edges of the original hologram are still visible. In comparison, the amplitude of the complex-valued field at the detector after the last iteration shown in Fig. 4(g) has a smoother appearance than the extrapolated normalised hologram. The discrepancy between Fig. 4(f) and Fig. 4(g) can be explained by the fact that the extrapolated hologram is created by substantial, and therefore not artefact-free, data analysis, whereas the amplitude at the hologram plane after each iteration approaches that of an ''ideal'' hologram.
The reconstructed absorption and phase-shift distributions from the extrapolated hologram are shown in Fig. 4(h) and Fig. 4(i). The phase-shift distribution reconstructed from the extrapolated hologram ( Fig. 4(i)) reveals additional details of the inner structure that are not observed in the reconstruction of the normalised hologram ( Fig. 4(d) and Fig. 4(e)). For example, the vertical line that can be seen in the middle of the tissue (indicated by a blue arrow in Fig. 4(i)) is apparent in the reconstruction of the extrapolated hologram and can be attributed to the fibrosis in the tissue. Fibrosis is often an indication of cirrhosis, which when untreated can develop into HCC. Thus, in our sample, the phase image contains more structural information, which can be analysed for indications of cancer.
To cross-validate our results, we also performed imaging of healthy liver tissue, which was placed at a distance of 36.6 mm from the detector. The imaging and reconstruction was done in exactly the same manner as for the cancerous tissue. The results are shown in Fig. 5. The reconstructed phase-shift image of the healthy liver tissue (Fig. 5 (h)) has a more uniform appearance than that of the cancerous tissue ( Fig. 4(i)).
Estimation of resolution. Quantitatively, the increase in resolution can be estimated as follows. The intrinsic resolution of in-line holography is given by R~l d ND 18 , where l is the wavelength, d is the distance between the sample and the detector, N is number of pixels, D is the pixel size, and ND 3 ND is the size of the hologram. For terahertz in-line holography at sample-to-detector distances greater than 3 mm, the resolution is given by the intrinsic resolution 33 . In the formula for calculating the resolution 18 , the hologram size is limited to the area where the interference pattern is observed. For the single holographic record shown in Fig. 3, as well as for the upsampled hologram shown in Fig. 4, the size of the hologram is 12.4 3 12.4 mm 2 and the resolution is R 5 240 mm. After extrapolation, the effective size of the hologram has been increased to 18.75 3 18.75 mm 2 with a resolution of R 5 158 mm.
Acquisition time and hydration level. To investigate the change in the sample hydration during the acquisition of in-line holograms as a function of acquisition time, we recorded holograms of a sample of 60-mm-thick healthy mouse liver. The sample was twice as thick as the human samples; thus, the absorption was expected to be stronger. The entire period of data acquisition was less than 1 min, including retrieval of the sample from the saline. In less than 5 min (starting from the time the sample was mounted on the slide), the samples completely rolled up. For timing reasons, we thus recorded 1,000 frames of holograms without shifting the detector. The first 100 (1-100) and the last 100 (901-1,000) frames were averaged into two holograms with good signal-to noise ratio, which were then reconstructed. The results are shown in Fig. 6. It is obvious from a comparison of Fig. 6(b) and Fig. 6(e) that the maximum absorption for the first hologram (frames 1-100) is higher than the maximum absorption for the second hologram (frames 901-1,000). This finding demonstrates that the sample quickly dehydrated during the acquisition procedure. Once a better detector with better dynamics and larger number of pixels is available, the number of frames, and therefore the acquisition time, can be further reduced.
This efficiency, in turn, can enable better preservation of the hydration level of samples.

Discussion
With the method proposed here, the absorption and the phase-shift distributions were simultaneously reconstructed from a normalised holographic record, and both distributions can contribute to the diagnosis of early-stage cancer, for example, by imaging an indication of fibrosis in the liver. The reconstructed absorption is low, which indicates that the sample is transparent to the terahertz waves and that imaging of the internal structure of the sample is possible.  Thus, in terahertz in-line holography, the phase image exhibits more structure and can be used for diagnosis. Further development of terahertz array detectors with larger areas and smaller pixel sizes will eliminate the need for recording and averaging over a sequence of holograms, and a single-shot hologram will be sufficient for highresolution reconstruction of the sample, which will ultimately reduce the acquisition time to just a few seconds and thus allow the preservation of the natural hydration level. Conventional cancer diagnosis by biopsy can be complemented by the method proposed here because the tissues extracted for biopsy have the optimal size for terahertz digital holography. Overall, we believe that these results will encourage the development of terahertz digital holography for further biomedical applications.

Methods
Statistical averaging of individual frames. Owing to the thermal working principle of the pyroelectric detector, zero or negative values appeared at some pixels in some frames. For each pixel, a histogram of its intensity values (bin size 5 1,000 1/3 5 10) was plotted and fitted with a Gaussian distribution function, and the mean of the Gaussian distribution function was considered to be the pixel value 13 . The pixels whose values were zero for all 1,000 frames were marked as ''dead'' pixels; they were set free of constraints during the phase retrieval, which allowed their values to be retrieved.
Hologram normalisation. The averaged hologram H 0 was divided by the averaged background B 0 to obtain the normalised hologram. Owing to the inhomogeneous absorption of the quartz slide and the slightly unstable intensity of the terahertz beam, the intensity of H 0 /B 0 did not vary from a value of approximately 1 at the edges of the hologram, where the amplitude of the object wave was very weak and only the reference wave was detected. To solve this problem, H 0 /B 0 was smoothed 1,000 times using a low-pass filter, which gave B9, and was then divided by B9 to give H Normalised 5 H 0 /(B 0 B9). Fourier domain filtering was applied to eliminate the horizontal fringes created by the detector scanning procedure.
Hologram reconstruction. The complex-valued transmission function of the sample is described by the function t(x,y) 5 exp[2a(x,y)2iQ(x,y)], where a(x,y) is the absorption and Q(x,y) is the phase shift introduced by the object in the passing wave 25 .
The transmission function of the sample t(x,y) was reconstructed by the propagation of the optical field from the detector plane backwards to the object plane using the angular spectrum approach 32 : where FT and FT 21 are the Fourier transform and inverse Fourier transform, respectively, and f x and f y are the spatial frequencies. The distribution of the reconstructed transmission function was refined by applying iterative reconstruction 25 , where the constraint of positive absorption values was applied.
Extrapolation procedure. Extrapolation is similar to the procedure explained elsewhere 29 . For the first iteration, the complex-valued distribution in the recording plane was created as follows. The amplitude was obtained by padding the square root of the normalised hologram H Normalised (size of 992 3 992 pixels) up to a size of 1,500 3 1,500 pixels with noise. The phase distribution was assumed to be randomly distributed in the range [2p, p]. The phase retrieval procedure was based on the error-reduction algorithm 31 applied for 200 iterations, and the field propagation between the detector and the object plane was calculated by using the angular spectrum method 18 . The following constraints were applied in the detector plane: the amplitude of the central 992 3 992 pixel region was replaced by the square root of the hologram. The amplitude values at the remaining pixels and the phase values at all pixels were updated after each iteration. In the object domain, the constraints of a positive absorption 25 and a supporting mask were applied.