Photothermal tomography for the functional and structural evaluation, and early mineral loss monitoring in bones

Salient features of a new non-ionizing bone diagnostics technique, truncated-correlation photothermal coherence tomography (TC-PCT), exhibiting optical-grade contrast and capable of resolving the trabecular network in three dimensions through the cortical region with and without a soft-tissue overlayer are presented. The absolute nature and early demineralization-detection capability of a marker called thermal wave occupation index, estimated using the proposed modality, have been established. Selective imaging of regions of a specific mineral density range has been demonstrated in a mouse femur. The method is maximum-permissible-exposure compatible. In a matrix of bone and soft-tissue a depth range of ~3.8 mm has been achieved, which can be increased through instrumental and modulation waveform optimization. Furthermore, photoacoustic microscopy, a comparable modality with TC-PCT, has been used to resolve the trabecular structure and for comparison with the photothermal tomography.


Introduction
The history of skeletal imaging extends back to the discovery of X-rays in 1896, the year in which Wilhelm Conrad Röntgen published an article incorporating 'a shadow of the bones of the hand' back projection by 'a new kind of rays' on a photographic plate [1,2]. With the advent of digital computing and computed tomography, planar X-ray projections in two dimensions acquired three-dimensional form and progress in diagnostic skeletal radiology has been phenomenal ever since [3,4]. In the context of bone mineral density (BMD) estimation and general skeletal health screening, dual energy X-ray absorptiometry (DEXA) continues to be the gold standard while highly detailed information is provided by its variants: quantitative computed tomography (QCT) and micro-computed tomography (µCT) [5][6][7]. X-ray methods are limited to structural evaluation only and this limitation has been addressed by proposing mixed modal approaches, which couple this with other modalities, like positron emission tomography (PET), single-photon emission computed tomography (SPECT), and magnetic resonance imaging (MRI) [8][9][10][11]. In these cases, detailed functional and quantitative assessment is possible for bone tissues. Diagnosis of angiogenesis, tumor growth and metastasis, enzyme activities etc., could thus be effectively accomplished in addition to the BMD estimation and osteoporosis (OP) tracking. Despite the vast potential for skeletal diagnosis, X-ray-based ionizing methods are criticized for their adverse effects on tissues and poor sensitivity especially in the case of osteoporosis diagnosis. Furthermore, X-rays are attenuated only by the mineral component of bones while the organic component remains unaccounted for. Since the true bone strength depends on both these components, X-ray BMD analysis yields an incomplete picture [12]. Another demanding situation is the BMD monitoring of astronauts under microgravity environments, as disuse osteoporosis puts them at a higher risk of fracture on return to Earth [13,14]. Portable and non-ionizing methods suitable for continuous/frequent bone strength monitoring are preferred in this situation. Quantitative ultrasound (QUS) and mechanical response tissue analysis (MRTA), though non-ionizing counterparts of X-ray methods for BMD estimation, are short on reliability and reproducibility [15][16][17].
The past few decades have witnessed significant developments in the medical applications of light, both in diagnosis and in therapy. Since the sensitivity and specificity of many of these methods are inherently high or can be improved through the use of proper contrast agents, the advantage of early detection of diseases has been a biophotonic highlight [18]. Bioluminescence and biofluorescence are successful agents for tracing bone metastases, tumor growth, angiogenesis, pro-angiogenic signaling, etc [19][20][21]. However, poor spatial resolution, inability for discriminating the signal from cortical and trabecular sections, and the depth integrated rather than depth resolved nature, are some of the shortcomings of these optical techniques. Optical coherence tomography (OCT), although it offers high resolution and depth resolved visualization, has its probe range limited to <1 mm in bone [22,23]. Another attempt made to resolve bone structure was by using diffuse optical tomography, however, its resolution was too small to observe the trabecular network [24]. Using direct measurement of scattered light although a correlation was established between the scattered image density and BMD, depth-resolved analysis was not feasible [25]. Since the extracellular bone matrix contains protein-rich molecular vibrational modes, mostly due to collagen components, Raman and infrared (IR) spectroscopies have been used to study collagen quality, crystallinity, and mineral-to-matrix and carbonate-to-phosphate ratios, among other bone properties. This in turn leads to the assessment of the variation of bone composition with aging or diseases [12,[26][27][28].
Apart from these methodologies, no effort to develop an optical equivalent of QCT or µCT has yet been reported, capable of resolving the trabecular network in three-dimensions through the cortical region with and/or without soft-tissue overlayers. Detecting bone mineral loss in a depth resolved manner is another desired competence. In this article, we discuss some significant features of the first hybrid-optical tomography of bones featuring opticalgrade contrast and photothermal (PT) depth profilometry capable of resolving the trabecular structure through the cortical overlayer in a depth-selective manner [29]. Furthermore, a nonionizing quantitative volumetric marker has been defined, validating its absolute nature, to monitor the variation in BMD. The potential of this marker for the early detection of mineral loss has been investigated and the results are compared with the µCT-measured BMD variations. Discussions have been extended to the response in bone, with a thin soft tissue overlayer, containing intentionally injected blood within the trabecular pores to mimic the presence of marrow. Yet another useful feature, the feasibility of BMD-scaled visualization, has been demonstrated in a mouse femur. The method is fast with an analysis time of a few tens of seconds and does not require any kind of relative motion between the bone and detector.

The photothermal effect
When light modulated at an angular frequency ω impinges on a material, there will be periodic absorption of energy by molecules or atoms. The non-radiative relaxation of molecules/atoms results in the periodic heating of the material which induces diffusive heat flow from the region of optical absorption into the colder bulk of the material or into its surroundings. This periodic thermal diffusion is called a thermal wave (TW), it is spatially damped, and it further attenuates with distance from the source due to optical absorption and scattering. The depth to which a TW travels until its amplitude is attenuated by exp(−2π) with respect to its surface value is called thermal diffusion length, l d (ω) [29]. The thermal diffusion length depends on the excitation frequency, density (ρ), specific heat capacity (c) and thermal conductivity (K) of the material and is given by It follows that, for a given material, a TW at a lower frequency travels longer than at a higher frequency. Due to their diffusive nature, the TW amplitude and phase measured at the sample surface will be functions of sub-surface thermo-optical properties and discontinuities inasmuch as they disrupt or perturb the flow of optically generated heat. For a homogeneous scattering medium excited by a pulsed laser beam the PT relaxation signal recorded by an IR detector, assuming 1-dimensional thermal diffusion, is given by [30] where, μ a , μ IR , μ tr and μ eff are the optical absorption coefficient at the excitation wavelength, IR absorption coefficient, transport coefficient and effective attenuation coefficient 30. Although bone is a highly heterogeneous medium, the general features of its PT response can be derived from Eq. (1). All or most of the parameters appearing in this equation vary due to demineralization with more or less sensitivity and one can expect corresponding changes in the PT signal as follows.
(i) Laser absorption and scattering coefficients: Both these parameters decrease with demineralization over 520-960 nm wavelength range [31]. Decrease in these parameters leads to the reduction of PT amplitude [30]. In this article, we deal with only the conductive coupling of PT signal and hence the variation of IR absorption coefficient with demineralization is neglected.
(ii) Density: Demineralization leads to the removal of bone material. Assuming the rate of mass reduction is larger than that of volume, a fall in the physical density would result. So, the PT amplitude should increase with mineral loss. No reliable data on the variation of bone specific heat capacity with mineral loss is available.
(iii) Bone volume: Prolonged artificial demineralization may lead to trabecular collapse leading to bone volume shrinking. If the rate of this volume compression dominates the mass-loss rate then the density may increase, which would reduce the PT amplitude.
(iv) Void fraction: For low enough demineralization that doesn't lead to trabecular collapse, the void fraction may increase. Since ρc, the volumetric specific heat capacity, of air ( 3 1.2 10 × ) is much lower than that of bone ( 5 8.36 10 × ) the PT signal will be amplified by the increasing porosity. So, the PT signal strength for a specific demineralization level is controlled by the complex interplay of the above mechanisms, the relative contributions being variable.

Photothermal radar and truncated-correlation photothermal coherence tomography (TC-PCT)
In contrast to electromagnetic or acoustic waves which are described in terms of hyperbolic differential equations, photothermal diffusion waves lack wavefronts and are reflection-and refraction-less in nature [32]. Detection of distance-integrated, rather than localized, distributions of energy is a unique characteristic of these parabolic diffusion-wave fields. This results in poor axial resolution which deteriorates with time and distance from the source and is a serious limiting factor when depth-resolved or 3D visualization of the probed medium is considered [33,34]. Several efforts have been made to achieve energy localization in a thermal wave field. The photothermal radar [35,36] is a prominent class of these methodologies. If R(τ) and S(t + τ) are the excitation (reference) and response (signal), respectively, the latter being measured at the delay time t then the cross-correlation defines the radar output. Here, ( ) * R τ is the complex conjugate of the excitation signal, τ being the dummy variable. Fundamentally, these energy localization approaches utilize matched filtering using pulse compression, which is a traditional radar technique for enhancing the range resolution and signal to-noise ratio (SNR) in a hyperbolic wave field [37]. Here, the coded signal can be described by the frequency response H(ω) of the coding filter. The frequency response of the matched filter which receives the signal is the complex conjugate H*(ω) of the coding filter response. The output of the matched filter is the inverse Fourier transform of the product of the signal spectrum and the matched filter response: Recently, we introduced a new kind of photothermal imaging called the truncated-correlation photothermal coherence tomography (TC-PCT) capable of resolving the depth distribution of thermo-optical parameters of an optically excited medium in a layer-by-layer fashion along the direction of light propagation [38]. This modality enables the reconstruction of a threedimensional picture of the sub-surface thermo-optical features of the medium. A pulsed diode laser triggered by a custom developed pulse generator is the excitation source ( Fig. 1). The current through the laser is controlled by a laser driver. The laser output is coupled to an optical fiber and is intensity homogenized before irradiating the bone sample. A 2-5 μm spectral band IR camera (Cedip 520M) captures the evolution of the laser-induced thermal signature from the bone surface at a rate of 370 frames/second. The entire system is under the control of a computer, which executes the TC-PCT algorithm in LabView environment.

Principle of TC-PCT and volumetric reconstruction of the optothermal structure of bones
The recently introduced chirped pulse PT radar exhibits high depth-resolving and penetrating capabilities along with excellent sensitivity [39]. Therefore, it is capable of acceptable SNR using low enough laser intensity to meet the laser safety requirement for maximum permissible exposure (MPE) in biological media [40]. The TC-PCT concept was derived from this pulsed PT radar. A linear frequency modulated (LFM) train of short optical pulses (E) is used to excite the sample ( Fig. 2(a)). In TC-PCT, truncated reference chirps of pulsewidth W T (W T < excitation pulsewidth) with continuously variable delay/initial-phase (d) are generated. Each of these chirps is cross-correlated with the PT image (S) captured by the camera in a pixel-by-pixel manner to produce two-dimensional axially resolved layer-by-layer images of the sample. For example, cross-correlation of the reference R1, which is in phase with E, with the PT image chirp, will give a surface picture of the sample. In the case of R2, the resulting image will be that of a subsurface thin layer. As d increases, the depth of the layer below the surface also increases. Furthermore, the depth resolution (thinness of an image layer) will increase as W T shortens. The TC-PCT is thus analogous to a bandpass filter whose operating frequency and quality (Q) factor are variable. The PT relaxation signal, following a pulsed excitation, is continuous in the frequency domain and during TC-PCT execution the frequency contributing the planar image keeps on decreasing with the increase in d. Parameters controlling the maximum penetration depth are the starting LFM frequency and period/length of the chirp.
The algorithm is implemented in the frequency domain in LabView environment ( Fig.  2(b)): The excitation chirp is first converted to a square wave chirp and then passed through the delay-sweep module. The delay-incremented square wave chirp (C 1f,0 ) is frequency doubled (C 2f,0 ). Subsequently, C 1f,0 and C 2f,0 are subjected to the binary exclusive-or (EX-OR) operation to generate a quadrature square chirp C 1f,90 . Truncated references R 0 and R 90 are synthesized from C 1f,0 and C 1f,90 , respectively. The complex conjugate of Fourier transformed truncated-references, both in-phase and quadrature, is multiplied with the Fourier transform of the PT chirp and the products are inverse-Fourier transformed to generate the in-phase and quadrature TC-PCT output from which the amplitude and phase can be evaluated.
We chose an 808 nm laser (Jenoptic JOLD-120-QPXF-2P) for its deep penetrating capability in tissues. It has been proven that for the highest SNR, a short optical pulse of high peak power should be used in TC-PCT [39]. Unless otherwise specified, a pulsewidth of 10 ms with ~50 W peak-power was considered for excitation. The starting and ending frequencies for LFM were 0.2 and 0.6 Hz, respectively with a period of 12 s (henceforth termed chirp-1 waveform). With the above parameters and ~3 cm 2 illumination area the chirp energy density is ~0.8 Jcm −2 ; the MPE ceiling for soft tissue is ~3 Jcm −2 and is even higher for bones [40]. The truncated chirp pulse width was fixed to 2 ms, which was the temporal resolution of our camera.

Planar-image sequence formation and three-dimensional reconstruction
A sequence of 2D truncated-correlation amplitude images at increasing depths below the surface can be generated by incrementing the value of d. Goat ribs were subjected to natural meat decomposition by keeping in water. Clean ribs were washed, dried and cut into suitable lengths. Typical 2D images obtained for a rib (~5x1.6x0.45 cm 3 ; average cortical thickness ~0.7 mm) with chirp-1 excitation are shown (Fig. 3). Imaged dimensions are 4 and 3.2 mm along the X and Y directions, respectively. W T was 2 ms. The zero-delay (0 ms) corresponds to the surface layer and a delay of 128 ms defines a plane that is ~1.6 mm below the surface. These sequences can be stacked to generate a three-dimensional visualization and the stack can be binarized to reconstruct a binary volume tomogram (BVT), Fig. 4 [38,41]. BVTs are independent of the absolute PT amplitude and they depict the presence or absence (binary 1 or 0) of thermal wave power within the sample. ImageJ software was used for all image processing tasks [42]. For binarization, unless otherwise specified, we used Huang algorithm [41-43], which offered the optimum noise immunity and sensitivity to demineralization. BVTs were recorded for demineralized ribs to investigate the potential of this modality for monitoring early mineral loss noninvasively. For that study, the sample was etched in ethylenediaminetetraacetic acid (EDTA) solution, which was prepared by diluting 0.5 M acid with equal volume of distilled water [44]. The time of demineralization was considered as the independent variable. Bone is composed of an organic phase that contains (i) bone forming cells (osteoblasts and osteocytes), (ii) bone resorbing cells (osteoclasts) and (iii) a matrix of collagen and non-collagenous proteins (osteoid), and a mineral phase of inorganic salts. Osteoid contains type-I collagen (~95%) and ground substances (~5%). The bone mineral, calcium hydroxyapatite [Ca 10 (PO 4 ) 6 (OH) 2 ] of low crystallinity, is tightly embedded within the organic matrix. Normal calcified bone contains ~25% of collagen matrix, ~5% noncollagenous protein and ~70% mineral. On optimum EDTA treatment, the mineral content dissolves in the solution leaving the organic matrix and the cellular components undisturbed. In addition, traces of fat, hemoglobin and water may also be present in the demineralized sample. None of these constituents has specific absorption peak at 808 nm and hence the photothermal excitation at this wavelength is due to the collective absorption of these components.

Axial and lateral resolution of bone TC-PCT
The PT signal, in general, carries information about the sample's optical and thermal properties in a composite manner. So, the nature of TC-PCT image formation cannot be attributed to the optical inhomogeneity alone. Instead, the coupled opto-thermal response determined by the relative magnitudes of the thermal diffusion length, optical penetration depth and sample thickness characterizes this process. Hence, TC-PCT is a functional imaging modality rather than a structural imaging platform. Since l d is inversely correlated with the square root of ω, higher resolution would be possible at higher frequencies that demand shorter excitation pulsewidth. However, this will limit the depth range. Due to the cross-coupled optothermal signal generation mechanism, the axial and lateral resolution of TC-PCT were determined experimentally, the easiest approach [38]. An absorbing graphite mark was made on the back surface of a mechanically polished, ~100 μm thick, cortical goat femur slab (30x10 mm 2 ). It was wetted to prevent direct coupling of IR photons. TC-PCT was performed by irradiating (chirp-1, 808 nm) the front surface for different values of W T : 2, 4, 6, 8, 10 and 12 ms at different locations. On an average, the back absorber mark appeared in the TC-PCT image for W T ≥ 8 ms. For 200 μm thickness this occurred at ~16 ms. Extrapolating this response a bone axial resolution of ~25 μm for 2-ms (camera resolution) truncated-reference was estimated. The lateral resolution is ~100 μm. For a given material, imaging area and chirp parameters, the lateral resolution is controlled by the camera resolution.

Binary-amplitude volume TC-PCT of intact and artificially demineralized bones
A rib sample (~50x16x4.2 mm 3 ; ~0.7 mm cortical thickness) was demineralized for known time intervals (0.5, 1, 2.5, 5, 10, 20 and 30 hours) and chirp-1 TC-PCT analysis was performed after each treatment (W T = 2 ms) to generate 120 planar images. Imaged area was 4x3.2 mm 2 . The first 30 images were used for reconstructing the BVT of the cortical layer and the next 90 sequences for the trabecular region. This ratio of sequences was determined by visually examining the onset of the trabecular signature in 2D images. The combined cortical-trabecular tomogram was made using the entire 120 sequences. Results obtained for the intact (A1), 0.5-(A2) and 20-(A3) hour demineralized samples are shown in Fig. 4. The irradiated bottom cortical layer (C1-3) and the approx. 2-mm thick trabecular section above it (B1-3) were reconstructed independently and in tandem. The combined tomogram was sliced along PQ (marked on the trabecular section) in the (X-Z) plane and the rear half is projected in D1-3 to see the interior. Tomograms E1-3 are the counterparts of D1-3, binarized using mean-threshold algorithm instead of Huang's. Images F1-3 are obtained by adjusting the threshold, to see the cortical geometry clearly, in absolute amplitude TC-PCT from which D1-3 and E1-3 are derived. The relevance of the last two series will be discussed in Sect. 4. Cross-section of the combined cortical-trabecular regions obtained by slicing B1-3 along PQ in the (X-Z) plane (see text). Laser irradiation was shone from the bottom cortical surface, which was the plane of imaging. B, C and D series were generated using Huang binarization algorithm. E1-3: Counterpart of D1-3, binarized using mean-threshold algorithm. F1-3: Threshold-adjusted absolute amplitude TC-PCT from which D1-3 and E1-3 were derived.

Thermal wave occupation index (TWOI): A new marker for non ionizing demineralization monitoring and its absolute nature
Since the thermal wave field properties are strongly thermal diffusivity dependent which is a function of density, a measure of the volume occupied by the thermal wave power should be indicative of the mineral density distribution (or in general the degree of demineralization) within the bone. The TWOI is defined for a BVT of arbitrary volume. It is estimated for a regular volume envelope within the bone specimen ( Fig. 5(a)). It is defined as the ratio of the volume occupied by the binarized thermal wave field (u) to the total volume of the quadrangular envelope (V) surrounding the photothermally activated region of the bone: In the illustration, the bluish box is the volume enveloping the binarized thermal wave distribution (yellow). For a given excitation conditions, the binarized TC-PCT properties of a sample are controlled by W T and the total delay (d max ). d max refers to the phase difference between the first and the last truncated-references and it determines the probed depth. The influence of W T and d max on the estimation of the TWOI of a sample has been investigated by reconstructing the BVTs of a rib sample, excited using chirp-1, for W T varying in the range 2-100 ms and d max fixed at 200 ms. In all cases, to ensure sweep (axial) continuity, the delay increment was the same as the pulsewidth. d max = 200 ms corresponds to a depth of ~2.5 mm this kind is naturally expected as TC-PCT is analogous to a continuously tunable bandpass filter: the quality factor improves as the truncated-pulsewidth shortens and central or operating frequency decreases with the increase in delay [38]. TWOI is analogous to the filter output observed in the time domain after excitation. Output settling time for low frequencies would be longer than that for high frequencies for a given filter time constant. d max should be increased proportionally with W T to enable stable (steady or non-oscillatory) evaluation of TWOI. For the response in Fig. 5(b), the stability is better than 98.5% if W T ≤ 6 ms. Here, stability refers to the constancy observed in the estimated value of TWOI over the interval 0-W T in the plot. This leads to the conclusion that, for a given excitation chirp, the TWOI is independent of W T if d max >> W T . In other words, the axial resolution should be made much higher than the depth-of-interest for the reliable estimation of this marker. We investigated the influence of excitation pulsewidth and pump intensity also and found that if the SNR remains unaltered, TWOI will be independent of these parameters as well. All these observations point towards the absolute nature of TWOI. We have evaluated the TWOI of a set of goat rib samples as a function of demineralization time (1, 2.5, 5, 10, 20, 30 and 40 hours in EDTA solution). Based on the fact that natural early mineral loss (osteoporosis) mainly occurs in the trabecular region [5], BVTs were constructed for ~1 mm thick trabecular region. TWOI was estimated for a volume of ~3x2x1 mm 3 before (TWOI B ) and after (TWOI A ) demineralization, and the fractional loss, ( ) TWOI TWOI TWOI − × , was calculated. Furthermore, μCT (Skyscan 1173, Bruker, USA) analysis was carried out to evaluate the BMD of these samples, intact and demineralized. Figure 7 shows the variation of fractional loss in TWOI and BMD with demineralization time. BMD follows a monotonic response up to 30 hours of demineralization while for TWOI this behavior exists up to 20 hours. However, for low demineralization levels (up to 5 hours) that is likely to occur naturally in osteoporotic bones, the sensitivity of TWOI is ~5.6 times the BMD. The TC-PCT sensitivity to the changes in organic phase as well can be a reason for this advantage. For very high demineralization both the fractional losses undergo a reversal in direction. In the case of BMD this is due to a failure in the μCT procedure: the protocol defined for analyzing samples of mild demineralization did not converge for the 40-hour sample. Fractional loss in TWOI reverses sign for the 30and 40-hour samples. Fig. 7. Variation of fractional loss in μCT-estimated BMD and TWOI with demineralizationtime for the goat rib sample (trabecular region). TC-PCT is ~5.6 times more sensitive than μCT at low demineralization levels.

Intensity/BMD-scaled bone tomography
A unique characteristic of absolute amplitude TC-PCT is the feasibility of intensity scaling by which the tomogram can be modified to show the regions of a specific band of PT amplitudes. Since the PT signal amplitude is directly correlated to the BMD, this feature turns out to be a means for the selective mapping of regions of a specific BMD range. We have examined this possibility using a mouse femur (Fig. 8). The illuminated cortical surface was the plane of imaging (BF). The sequential two-dimensional images are stacked, without binary conversion, to form the amplitude tomogram, TF. Images TR1-5 are the rear views of the amplitude Fig. 8. Variable-threshold TC-PCT amplitude tomograms for a mouse femur. Selective mapping of regions of a specific BMD range is possible by adjusting the threshold. Amplitude tomogram (TF) is developed by irradiating the cortical surface in BF. TF is rotated to view from the rear side (TR1-5). BR is the rear view photograph.
tomogram TF for different amplitude thresholds: 0, 20, 40, 60 and 80% of the maximum amplitude, respectively. BR is the corresponding rear view photograph. Threshold selection enables the blocking of PT sources below it from contributing to the image. As the threshold approaches 60% of the maximum signal only the hard cortical portions appear in the image (TR4). The bluish background is due to the ambient IR signal.

Discussion
Limitations of traditional bone diagnostics/imaging techniques such as tissue ionizing disadvantage that forbids continuous monitoring, poor sensitivity and the impracticality to monitor disuse osteoporosis in space missions are the key motivators to search for alternate methods in this area. A related issue is the lack of sensitivity of traditional techniques to the biological conditions of bones resulting from changes in the organic phase. To address all these challenges together, our present efforts have been geared to demonstrate a photothermal 3D imaging modality featuring optical grade contrast for general skeletal diagnosis with the specific axial and lateral spatial resolution abilities for BMD monitoring and early mineral loss detection in micro-gravity environments as well as on Earth.
With TC-PCT, the irradiation energy is within the MPE ceiling and typical analysis duration is a few tens of seconds. It has the advantage of 3D tomographic reconstruction without the need for relative motion between the object and detector, which is the case with QCT or µCT. The depth-resolved imaging capability enables computational sectioning of cortical and trabecular regions for independent analysis, which is a remarkable feature also offered by µCT. Since the attenuation of thermal waves is stronger than that of X-rays in tissues, TC-PCT has a depth range limited to a few millimeters for excitation frequencies in the range 0.1-1 Hz. This range could be extended by using lower chirp frequencies, therefore slowing down the analysis with a compromise on the SNR and tissue dehydration. In a bonetissue matrix, with 24 s long chirp, we have observed a depth range of ~4 mm.
The potential of TC-PCT for depth selective imaging has found immediate application in the reconstruction of cortical and trabecular bone regions separately, and in the analysis of the relative influence of artificial demineralization in those regions. It is obvious from Figs. 4(B1) and 4(B2) that the modification to trabecular structure following very short demineralization time (0.5 hours) is marginal, yet observable with TC-PCT. So is the case for cortical layers. However, after 20-hour demineralization, there appear complex modifications to both cortical and trabecular regions and partial collapse of the trabecular network ( Fig. 4(A3)). The TC-PCT volume appears expanded for both regions in this case. If this were true then the corresponding TWOI would have remained larger than its intact value according to Eq. (4). But what Fig. 7 reveals is that the TWOI keeps on decreasing with demineralization up to 20 hours. To address these contradictory results, first we analyzed the absolute PT amplitude for the intact and demineralized conditions in a layer-by-layer manner taking the sum over every layer. It was found that the intact state has the highest PT amplitude and it decreases with the demineralization time up to 20 hours. This ensures the accuracy of the TWOI response in Fig.  7 and leaves questions on the image processing aspects of BVT. To further explore this issue, we generated BVTs for the intact, 0.5-and 20-hour samples using mean-threshold algorithm (another binarization technique) for the same surface-mesh parameter used for Huang method, Fig. 4(E1-3). On comparing D3 and E3, one can see that the apparent cortical expansion in D3 is due to the inferior resolution of the Huang binarization algorithm, compared to mean-thresholding, for TC-PCT reconstruction Although the mean-threshold algorithm offers better resolution, its noise immunity and sensitivity to demineralization are lower than Huang's. We further examined this discrepancy of apparent volume expansion by considering the absolute amplitude tomograms, Fig. 4(F1-3). Here, the threshold was adjusted to improve the cortical contrast. The cortical thickness appears almost the same in all the three cases. Fortunately, the volume expansion or filling appears in the image only, not in the binary data, and it is assumed to be due to some artifacts at the surface-mesh generating stage: The algorithm first generates a set of 3-dimensional binarized data corresponding to the absolute amplitude TC-PCT. Using these data, with user defined mesh parameters, 3D surfaces are constructed to render volume visualization. For the successful reconstruction of a BVT the spatial separation between adjacent binary data should be greater than the resolution of the algorithm. For the 20-hour sample, the collapsed trabeculae become more closely packed due to bone volume compression so that their binary spatial separation in all directions falls below the algorithm resolution at many locations. Now, adjacent structures may get joined together to form a single object filling up the inter-data space, thus leading to an apparent volume expansion in the BVT. Since TWOI is calculated using the binary data, not the reconstructed volume, it remains unaffected by this discrepancy. However, the efficient use of TC-PCT for structural and functional imaging, and demineralization monitoring in bones demands optimum binarization algorithm and reconstruction protocols to be developed.
As per Fig. 7, the PT signal power keeps on decreasing with demineralization time up to 20 hours. The primary reason is the fall in the laser absorption and scattering coefficients with mineral loss [31]. As another reason there can be increase in the physical density of bone, in case, if the rate of mass loss is lower than that of volume compression. For 30-and 40-hour demineralization, the fractional loss in TWOI reverses sign indicating the postdemineralization signal amplitude is larger than the intact value. A possible reason is, following prolonged demineralization, excessive removal of bone material may enlarge the void fraction, especially within the trabecular region. The decreased geometric dimensions of the trabecular network will result in increased thermal confinement and higher thermal-wave amplitudes under similar optical absorption conditions. This effect can be amplified by the poor thermal conductivity of air: the increased void volume may have impeded the conductive loss rate, thereby sustaining a stronger thermal wave field within the trabecular pores. Also, the loss in mineral content may lead to a reduction in the physical density if the rate of mass loss is higher than that of volume loss. These physical mechanisms set an upper limit on the mineral loss that can be monitored by TC-PCT, however, this is not an impediment to clinical applications, as the physiological osteoporotic bone loss never reaches the demineralization level of Fig. 4(B3). This hypothesis is further corroborated by the TWOI response, Fig. 7.
The ~3.8 mm depth range achieved with our present instrumentation enables the in-vitro analysis of small animal bones, for example, mouse bones. However, for in-vivo imaging, the strong IR photon absorption of soft tissue and skin overlayers must be overcome. We mimicked such a situation by placing a thin (~0.8 mm) layer of chicken breast tissue over a goat rib bone sample (~5x2x0.5 cm 3 ; cortical thickness ~0.8 mm) and the matrix was subjected to TC-PCT analysis with chirp-1 waveform (Fig. 9). The result was promising as the three layers (soft-tissue, cortical and trabecular sub-layers) were resolvable in the BVT. The thermophysical properties of soft-tissue and bone differ considerably and hence there exists thermal wave impedance mismatch at the boundary of these two-media [29]. The discontinuities or empty spaces appearing at the interface of bone and soft-tissue are assumed a consequence of inhomogeneous thermal wave impedance mismatch. We further investigated this feature by placing a tissue layer (~0.8 mm thick) on a thermally thick (~1.5 cm) tissue slab, both extracted from the same chicken breast piece. No such severe discontinuities were observed in the TC-PCT images for this sample configuration. In general, these tomograms reveal the high dynamic-range advantage of TC-PCT that enables the depiction of soft and hard tissues features simultaneously without sacrificing image contrast.  Because of its non-ionizing and non-invasive nature, and the depth-resolved 3-dimensional visualization advantage in biological media, photoacoustic (PA) imaging [46][47][48][49] is a comparable modality with TC-PCT. According to previous reports [50,51], and despite its high resolution, contrast, deep penetration and capability for depth-resolved imaging in soft tissues, photoacoustic tomography (PAT) performs poorly with bone tissues due to the multiple reflections and severe attenuation of ultrasound signals within the trabecular pores where high impedance-mismatch exists, especially at higher frequencies necessary for enhanced spatial and axial resolution.

Comparison with photoacoustic microscopy imaging
This, in turn, makes the existing reconstruction methods, which are not valid for velocity discontinuities, unsuccessful in resolving the trabecular pores. As an alternative we considered photoacoustic microscopy (PAM), which does not involve any reconstructionmediated image formation. The experimental setup consisted of a 690-nm, 8-ns pulsewidth OPO with 10-Hz pulse repetition frequency. The exposure energy density was ~0.03 mJ/mm 2 per pulse. PA detection was carried out at two frequencies, using 50-and 20-MHz transducers raster-scanned over the sample with a step size of 0.02 mm. Figure 10(a) depicts the PAM obtained at 50 MHz for the goat rib sample whose cross-sectional photograph is shown below it. For averaging, data were collected for 10 excitations at each point. Due to the strong attenuation of photoacoustically generated ultrasound signal at this frequency, neither the bone-water nor the cortical-trabecular interface is resolved clearly. PAM obtained for the more penetrating frequency of 20 MHz is shown in Fig. 10(b). Here the trabecular structure appears resolved along with the water-bone and cortical-trabecular interface.

Conclusions
We demonstrated a non-ionizing equivalent of μCT/QCT for bone diagnostics. This hybridoptical tomography of bones, truncated-correlation photothermal coherence tomography (TC-PCT), features optical-grade contrast and photothermal depth profilometry capable of resolving the trabecular structure through the cortical overlayer, with and without an additional soft tissue overlayer, in a depth-selective manner involving 2-dimensional image slice stacks forming a 3-dimensional image in a manner similar to OCT. This enables image slicing (sectioning) of cortical and trabecular regions for depth-specific analysis. Furthermore, a new non-ionizing quantitative volumetric marker, TWOI, has been defined to monitor BMD variation. If the truncated-correlation pulsewidth is << total reference phase increment then TWOI is independent of the pulsewidth. This marker was found to be more sensitive than the μCT-measured BMD for monitoring low demineralization levels in animal rib bones. This is very promising toward the application of TC-PCT to early osteoporosis diagnosis. Another notable feature of the bone TC-PCT is the intensity scaled bone tomography that enables the mapping of regions within a specific BMD range. Axial and lateral resolutions are found to be ~25 and 100 μm, respectively, which can be enhanced through instrumental improvements. The current depth range (~3.8 mm) is the highest in the field of photothermal imaging to-date and is sufficient to meet the requirements of small animal diagnosis. Deep-penetrating bone TC-PCT is expected to find applications in regular clinical diagnosis as well as in space missions for BMD monitoring of astronauts. As a comparable modality, we examined the potential of photoacoustic microscopy (PAM) to resolve trabecular features through the cortical layer. The results, though early, are promising. A major advantage of TC-PCT is the availability of advanced lock-in compatible mid-IR infrared cameras which can be used to implement the full imaging/diagnostic technology in a non-contacting manner, as opposed to transducer arrays and their fluid coupling requirement required for PAM.