1.

Introduction

The spatial resolution of optical transillumination imaging systems is severely degraded by light scattering. Several approaches such as confocal microscopy and optical coherence tomography (OCT) have been applied to reduce the degradation in imaging resolution and quality caused by light scattering. However, in skin, for example, (typical ${\mu }_s=100{\mathrm{cm}}^{-1}$, $g=0.9$, absorption coefficient ${\mu }_a=0.1–10{\mathrm{cm}}^{-1}$ for $\lambda \sim 400-700\mathrm{nm}$)¹ light scattering limits the imaging depth to $\sim 500\mu \mathrm{m}$ for confocal microscopy and $\sim 2\mathrm{mm}$ for OCT. In contrast, acoustic waves are more weakly scattered by tissue, and ultrasonic techniques have been widely used in medicine. However, ultrasonic techniques do not provide the optical spectroscopic information, e.g., sensitivity to blood oxygen saturation, provided by optical techniques. Therefore, beyond the penetration depth of confocal microscopy and OCT it is advantageous to combine the functional information provided by optical techniques with the resolution of ultrasound (US). This can be done using photo-acoustic tomography where light is used to generate ultrasonic waves in the tissue. This paper focuses on ultrasound modulated optical tomography (USMOT) where US is used to modulate light within the tissue.

Over the past two decades, the acousto-optic effect in scattering media has received significant attention particularly as a new biomedical imaging modality. This is because it has the potential to overcome the degradation in resolution associated with light scattering. The US changes the optical properties (refractive index, scattering coefficient, motion of scatterers).² These effects cause a phase modulation of the optical field often referred to as US “tagging”; and this light can be used to provide information from a local region within the turbid medium.

One of the first studies in USMOT was conducted by Marks et al.³ They succeeded in detecting the variation of optical intensity in an insonified region. The well studied modulation mechanisms of light traversing through an insonified region are caused by the movement of scatterers and changes in medium refractive index^1,4 under an applied acoustic field. An US transducer driven by continuous wave (CW) was used in many of the early imaging systems.⁵ Several detection methods have been proposed to increase the signal-to-noise ratio (SNR) such as using a single photo-sensor,^6,7 parallel detection using a camera,^8–10 or photorefractive crystal,^11–13 confocal Fabry-Perot,^14,15 or spectral hole burning.^16,17

Although CW systems are better understood and offer narrow bandwidths, they suffer from two main disadvantages. First, low acoustic spatial average intensity is required due to the safety standards in biomedical applications, which limits the modulation depth. Second, there is low axial resolution along the ultrasonic propagation axis as the detected signal is averaged along the focal zone of the US pressure field. In order to improve the axial resolution, a frequency-swept technique was introduced^18,19 where light is tagged with different frequencies at different axial positions; however, the modulation depth remains limited as it employs CW US. With pulsed USMOT techniques,^20,21 the US transducer is driven by a train of high peak voltage pulses. As the modulated optical field is linearly proportional to the acoustic amplitude,^2,22–25 it can provide a stronger modulation effect while the average energy still satisfies safety thresholds. In addition, a pulsed technique also improves the resolution along the US axis because the pulse can be time resolved providing optical information at specific positions due to the low speed of sound relative to light. In CW, frequency swept and pulsed techniques and assuming spherically focused US, the size of the focal zone is also dictated by the focus waist of the US pressure which is proportional to the transducer focal length, and inversely proportional to transducer element diameter and US frequency. Therefore, appropriate modification of these parameters may result in a tighter waist and can provide better lateral resolution for an imaging system. In addition, increasing the acoustic pressure to produce harmonics that have a tighter beam waist may also increase the imaging resolution.²⁶

Some modelling of US modulated light propagation in turbid media has been carried out. Sakadzic and Wang^22–24 proposed a theoretical framework for the mechanisms of producing US modulated light in tissue. The majority of research has involved Monte Carlo simulation of scattered light propagation through a CW US field.^2,23,27,28 Three mechanisms were identified² for US modulation of the optical field: (i) US induced variation of the optical properties (absorption and scattering) due to compression and rarefaction of the medium; (ii) variation of the optical phase due to US induced displacement of scatterers; and (iii) modulation of the optical phase due to US modulation of the index of refraction. (i) is a weak effect which can only be observed with low SNR when incoherent light is used. The contribution of (ii) and (iii) to the modulation depth of the detected signal is comparable when the acoustic wave number $k_a\le 0.559{\mu }_s$, where ${\mu }_s$ is the optical scattering coefficient. As $k_a$ increases above this value (iii) dominates and the modulation depth increases significantly.² Honeysett et al.²⁹ used a Monte Carlo model to investigate USMOT signals when bubbles are introduced in the scattering medium as contrast agents. Yuan et al.³⁰ proposed a diffusion model of fluorescence USMOT.

There have been fewer papers based on image reconstruction of USMOT data. Li and Wang applied the reconstruction method of X-ray computed tomography to an USMOT system where the US beam is scanned linearly and angularly across the sample.^31,32 Allmaras and Bangerth investigated a reconstruction based on a Green’s function approach.³³ Recently, Xu et al.³⁴ investigated time reversal techniques in which the optical properties of the sample at the US focus is reconstructed based on holographic data. Bratchenia et al.³⁵ applied a finite-element method to reconstruct information of optical absorbing objects in USMOT. Therefore, there is a need to develop new algorithms to reconstruct USMOT data.

One candidate is to apply a maximum likelihood (ML) algorithm.^36–41 Given a measured data set and a statistical model for the data, a ML algorithm will provide an estimate to the model parameters. The main advantage of the ML method is that it recognizes that the recorded data are the realization of a random process.^36–39 The ML method also does not require a priori information to be imposed on the solution. This has been successfully applied in many research areas, e.g., dynamic light scattering,³⁹ positron emission tomography,⁴⁰ and X-ray computed tomography.⁴¹

In this paper, we introduce a simple model of the pulsed USMOT signal based on the convolution of the US pulse with the optical field. Although there are limitations to this model, which are discussed later, it provides results that compare well with experimental data. Based on this model, a ML data inversion algorithm for image reconstruction is developed. To our knowledge this is the first time an USMOT image reconstruction technique by the ML method has been investigated. The reconstruction is investigated both theoretically and experimentally using a 10 MHz (central frequency) pulsed USMOT system. The experiments involve non-scattering water samples and scattering water and gel samples.

Section 2 introduces the theoretical description of the pulsed USMOT signal model, together with a description of the ML algorithm for image reconstruction. Section 3 describes the experimental configuration. The validity of the pulsed USMOT model is shown in Sec. 4 where detected signals in many scenarios are compared to the corresponding simulated results. Section 4 also presents results from the application of the ML method to both simulated and experimental data. Discussions and conclusions follow in Secs. 5 and 6, respectively.

2.

Theory

The next two sections describe the simple model that forms the basis for the ML algorithm followed by a description of the ML method itself.

2.1.

Model of Pulsed USMOT Signal

Figure 1 illustrates a model of a pulsed USMOT signal detected at a photo-detector. As an US pulse traverses the sample at a particular time, it introduces a pressure change (compression or rarefaction) at a particular volumetric element of the medium, which contains the US pressure field at a given point in time (defined as a ‘layer’ in this context). This modulates the motion of local scatterers and changes the local sample’s optical properties (scattering coefficient, absorption coefficient, and refractive index). Consequently, light at that layer is phase modulated by the US pulse and produces a phase modulated optical signal. When the US pulse reaches the next layer of the sample, it produces another optical pulse, which is similar to the temporal profile from the previous layer, but with a phase delay due to the time taken for the US pulse to propagate between the two consecutive layers. The speed of sound $v_a$ in water and in agar gel phantoms⁴² has been demonstrated to be approximately $1500\mathrm{m}/\mathrm{s}$ at room temperature. The speed of light $v$ in such media is expressed as $v=c/n$, where $c\approx 3\times {10}^8\mathrm{m}/\mathrm{s}$ is the light velocity in vacuum, and $n\approx 1.33$ is the refractive index of water or tissue. Because the speed of light is much higher than the speed of sound, the time taken for the modulated light to reach the detector is neglected. Hence, the phase difference only depends on the transit time of the US pulse, which produces a slowly changing envelope on the modulated signal. Some preliminary results of imaging small objects embedded in scattering media based on envelope detection method have been shown in Morgan et al.⁴³ Therefore, a pulsed USMOT signal at the detector may be expressed as a function of time, or as a summation of many phase shifted optical pulses.

Fig. 1

Model of a pulsed USMOT signal.

In the absence of absorption, a simple US pulse at the far field in an acoustically-uniform medium has the pressure $U(t)$ at a given layer which can be expressed as,

Eq. (1)

$$U(t)=s(t).\sin ({{\omega }_a}^{}t+{\varphi }_a),$$

where $s(t)$ is the pulse envelope that depends on the bandwidth of a transducer, ${\omega }_a$ is the US central frequency, and ${\varphi }_a$ is the US initial phase. The typical temporal pressure profile of a pulse generated by a 10 MHz US transducer (7.6 cm focal length, Olympus Panametrics-NDT V322) is shown in Fig. 2(b). This was obtained by measuring the amplitude of the acoustic pulse using the configuration shown in Fig. 2(a).

Fig. 2

(a) Pulse-echo setup: a pulser-receiver (Olympus Panametrics-NDT 5800) in pulse-echo mode to receive the reflected acoustic signal from a reflector placed at the transducer focus; (b) A normalized 10 MHz US pulse.

When a vibration source (e.g., US) is applied the small particles suspended within a solution^44,45 move deterministically. For small scatterers suspended in biological and engineered soft tissues, the relatively simple case where US is applied on a viscous flow with no-slip conditions can be considered.^22,44 Based on this assumption, the movement of optical scatterers may be expected to follow closely the movement of the surrounding fluid, but with smaller amplitudes and without any phase delay. These motions cause changes in the medium’s optical properties which then modulate light scattered within this region. The effect is more noticeable within the US focus where the US pressure is much stronger, and hence, the temporal profile of the US pressure is imposed on the temporal modulated optical signal. We define $O(t)$ as the optical profile from a given layer that contributes to the detected AC optical intensity. In this case $O(t)$ is linearly proportional to the US pressure $U(t)$,

Eq. (2)

$$O(t)=M.U(t),$$

where $M$ is a constant ($M=1$ in this paper) representing the efficiency with which the US amplitude is mapped onto $O(t)$ based on the complex interactions between US and optical scatterers. The linear relationship between US pressure and AC optical intensity occurs when there is homodyne detection due to interference of US modulated and un-modulated fields.^10,46

Now consider the US column to be composed of many such layers, each of width $\mathrm{\Delta }z=v_a\mathrm{\Delta }t$, a pulsed USMOT signal (AC signal) at the photo-detector is a summation of many delayed pulses,

Eq. (3)

$$I_{\mathrm{ac}}(t)=\sum _{j=1}^{m}p(z).O(t-j\mathrm{\Delta }t),$$

where $z=v_{a}j\mathrm{\Delta }t$, $\mathrm{\Delta }t$ is a time delay of the acoustic field, related to the number of steps $m$ along the US column. The profile $P(z)$ represents shape of the distribution along the acoustic axis as a result of the combined acoustical and optical characteristics along the US column. We propose a simple expression relating the optical and acoustic properties to the profile $P(z)$ which can be expressed as,

Eq. (4)

$$P(z)=P_{\mathrm{us}}(z).P_{\mathrm{opt}}(z).P_{\mathrm{abs}}(z),$$

where $P_{\mathrm{us}}(z)$ is the axial pressure profile of the US, $P_{\mathrm{opt}}(z)$ the scattered light intensity profile (in the absence of an absorber), and $P_{\mathrm{abs}}(z)$ is the optical absorption profile (indicating position and optical absorption of objects) along the US column. In this simple model, we assume that the optical and absorption profiles can be treated separately⁴⁷ although inevitably there will be some dependence. In the absence of absorption $P_{\mathrm{abs}}=1$ the profile represents the US modulated profile of a homogenous scattering medium.

The total signal $I(t)$ at the detector can be written as,

Eq. (5)

$$I(t)=I_{\mathrm{dc}}+I_{\mathrm{ac}}(t),$$

where $I_{\mathrm{dc}}$ is the unmodulated optical signal.

The pulsed USMOT signal is modelled by convolving a modulated optical pulse $O(t)$ with a profile $P(z)$. Based on the assumption that $O(t)$ initiated at a thin layer has a similar shape in the time domain to that of the US pressure, the US pulse shown in Fig. 2 is used to represent $O(t)$. In these simulations, the constant $M$ in Eq. (2) is set to unity.

Profile $P(z)$ can be modelled as a product of $P_{\mathrm{us}}(z)$, $P_{\mathrm{opt}}(z)$, and $P_{\mathrm{abs}}(z)$ [Eq. (4)]. Because a focused transducer is usually employed, a simplified axial US pressure profile at the far field can be modelled as a Gaussian distribution whose peak is at the US transducer focal zone, together with an attenuation factor depending on the acoustic properties of the medium. According to Challis et al.,⁴⁴ the attenuation of a US pulse travelling along the ultrasonic axis in the medium in the absence of diffraction losses, reflection losses, and focal gain can be expressed as,

Eq. (6)

$$A(z)=A_o{e}^{-\beta .z},$$

where $A_o$ is the un-attenuated magnitude of the propagating wave at an arbitrary location, $A(z)$ is the reduced magnitude after the wave has travelled a distance $z$ from the initial location, and $\beta$ is the attenuation coefficient of the wave travelling in the $z$-direction.

If the US transducer is placed at $z=0$ in the lab coordinates, we assume that $P_{\mathrm{us}}(z)$ can be expressed in a simple form as,

Eq. (7)

$$P_{\mathrm{us}}(z)=A_o{e}^{-\beta .z}{e}^{{-(z-f)}^2/2L_1^2},$$

where $f$ is the US transducer focal length, and $L_1$ depends on the properties of the acoustic lens. The term $L_1$ represents the length of the focal zone of the US column which is inversely proportional to the acoustic lens diameter. For example, a 2.5 cm-diam 10 MHz transducer with focal length of 7.6 cm has a focal zone ($-6\mathrm{dB}$) of 10 mm⁴⁸ hence in this case $L_1=10\mathrm{mm}$.

Assuming that a light beam illuminates the US focal region in an optical homogenous scattering medium, a simple estimate of the optical intensity spatial profile across the US region is a Gaussian distribution,

Eq. (8)

$$P_{\mathrm{opt}}(z)=B.e^{{-(z-f)}^2/2^{}{L}_2^2,}$$

where $B$ is the magnitude of the light intensity, and $L_2$ is the standard deviation of the profile. For example, in the experiments using an expanded beam in clear water described in Sec. 3.2 the beam is expanded to 20 mm so $L_2=20\mathrm{mm}$ in the simulation. The optical absorption profile $P_{\mathrm{abs}}(z)$ can vary from zero to one corresponding to totally-absorbing and non-absorbing objects, respectively. In the presence of an absorbing object, $P_{\mathrm{abs}}(z)$ can be written as,

Eq. (9)

$$P_{\mathrm{abs}}(z)=1-\sum _i^Q{G}_i.\mathrm{rect}\left(\frac{z-z_i}{D_i}\right),$$

where $Q$ is the number of absorbing objects, $G_i$ refers to the light attenuation along the US axis (focal zone), $G_i$ is the size of the object at position $z_i$, and $rect[(z-z_i)/D_i]$ is a rectangular function at position $z_i$ with the width of $D_i$.

There are several assumptions made in this model. For instance, although the detected experimental signal is homodyne detected, diffraction and interference effects (such as speckle) are neglected when estimating the optical intensity profile $P_{\mathrm{opt}}(z)$. In turbid media, light intensity fluctuates (speckles) due to motion (Brownian and deterministic, e.g., US) of scatterers. However, in the experimental system used in this paper, there is some averaging of speckles across the detector plane together with significant temporal averaging, and it is assumed that fluctuations in the optical intensity profile $P_{\mathrm{opt}}(z)$ are averaged out. Therefore, it is assumed that only the average intensity contributes to the optical intensity profile $P_{\mathrm{opt}}(z)$ along the US focal zone. This is also the case for fluorescence USMOT³⁰ where fluorescent light is incoherent and does not produce an observable speckle pattern. Equation (4) also requires that the optical and acoustic profiles in the tissue are known and that the absorbing and optical profiles are independent. In this case, we simply approximate them as Gaussian profiles [Eqs. (7) and (8)] and add absorption according to Eq. (9). This approach has been applied by others,⁴⁰ however, these can be more accurately provided either through measurement on test phantoms or through acoustic and optical simulations. Despite these approximations, encouraging results in Sec. 4.1 demonstrate the potential for using this model as a basis for a ML reconstruction algorithm.

The model suggests that a de-convolution technique can be employed to obtain the profile $P(z)$ from the signal, and then the optical absorption profile $P_{\mathrm{abs}}(z)$ along the US focus can be reconstructed. The next section introduces the ML data inversion algorithm and discusses the possibility of applying the ML technique to the model of pulsed USMOT signals for image reconstruction.

2.2.

Maximum Likelihood Data Inversion Algorithm for Pulsed USMOT Signals

Maximum likelihood estimation is a popular statistical method used for fitting a model to data and providing estimates for the model’s parameters. For a fixed data set and an underlying probability model, a ML method selects the values of the model parameters that can make the data “more likely” than any other values.^36–38 The main advantage of the ML method is its ability to take into account the random nature of noise.³⁹

We apply the ML method to the model of pulsed USMOT signal described in Sec. 2.1 and derive the iterative algorithm in this section. First of all, $I(t)$ in Eq. (3) is written in matrix form as,

Eq. (10)

$$\left[\begin{array}{c}{I}_1\\ I_2\\ ⋮\\ I_u\end{array}\right]=K+\left[\begin{array}{cccc}{H}_{11}& H_{12}& \cdots & H_{1m}\\ H_{21}& H_{22}& \cdots & H_{2m}\\ ⋮& ⋮& ⋮& ⋮\\ H_{u1}& H_{u2}& \cdots & H_{um}\end{array}\right].\left[\begin{array}{c}{P}_1\\ P_2\\ ⋮\\ P_m\end{array}\right]+\left[\begin{array}{c}{N}_1\\ N_2\\ ⋮\\ N_u\end{array}\right]\mathrm{.}$$

Alternatively, $I(t)$ can be expressed as below,

Eq. (11)

$$I_i=K+\sum _n^m{H}_{ij}{p}_j+N_i,$$

where $K$ is a constant (DC baseline), $I_i$ is the $i$th (where $i$ varies from 1 to $u$) data value at time $t_i$, $P_j$ is the value of the profile corresponding to the time delay of the acoustic pulse $j\mathrm{\Delta }t$, $N_i$ represents the noise in the intensity data $I_i$, and $H_{ij}$ is an ($i\times j$) matrix, defined by,

Eq. (12)

$$H_{ij}=O(t_i-j\mathrm{\Delta }t)\mathrm{.}$$

Figure 3 is a visual aid of Eq. (10) (excluding factors $K$ and $N$). The dashed line represents the detected signal $I(t)$, the solid line represents phase-shifted pulses $O(t-j\mathrm{\Delta }t)$, and the dash-dot line symbolizes the profile $P(z)$.

Fig. 3

Visual aid of Eq. (10): $\text{dashed}\text{line}=\text{detected}\text{signal}$ $I(t)$; $\text{solid line}=\text{phase-shifted pulses}$ $O(t-j\mathrm{\Delta }t)$; $\text{dash-dot line}=\text{profile}$ $P(z)$.

The probability that the data have the form $I_i$ with a profile $P_j$ is given by a product over each data value of the Poisson distribution.³⁸

Eq. (13)

$$p(I|P)=\prod _i{e}^{{-(K+H.P)}_i}\frac{{(K+H.P)}_i^{I_i}}{(I_i)!}\mathrm{.}$$

According to Bertero and Boccacci,³⁸ the basic rule in a ML method is to maximize the probability with respect to the values of data $P$. In this maximization procedure it is more convenient to calculate the logarithm of the probability that is usually called the logarithmic likelihood function $l(P)$;³⁹ this step is appropriate because maximizing an arbitrary function also maximizes its logarithm and $P$ is always a non-negative quantity. For Eq. (10), the logarithmic likelihood function has the form

Eq. (14)

$$l(P)=\sum _i\{I_i\ln {[K+H.P]}_i-{[K+H.P]}_i-\ln (I_i!)\}\mathrm{.}$$

The differentiation of $l(P)$ may be written as

Eq. (15)

$$\frac{\partial l(P)}{\partial P_j}=\sum _i\frac{\partial }{\partial P_j}\left\{I_i\ln \right[K+\sum _{j\prime }{H}_{ij\prime }{P}_{j\prime }]-[K+\sum _{j\prime }{H}_{ij\prime }{P}_{j\prime }]-\ln (I_t!)\},$$

which yields

Eq. (16)

$$\frac{\partial (P)}{\partial P_j}P=\left\{\sum _i\right[\frac{I_j{H}_{ij}}{K+{(HP)}_i}-H_{ij}\left]\right\}P\mathrm{.}$$

At the stationary point of $l(P)$, Eq. (16) is equal to zero, the solution gives

Eq. (17)

$$[H^T.(\frac{I}{K+HP}\left)\right]P=\alpha P,$$

where $\alpha =\underset{i}{\sum H_{ij}}$, $H^T$ is the transpose matrix of $H$. Now we define an operator $F(P)$ such that

Eq. (18)

$$F(P)=\frac{[H^T.(\frac{I}{K+HP}\left)\right]P}{\alpha }\mathrm{.}$$

Combining Eqs. (17) and (18),

Eq. (19)

$$P=F(P)\mathrm{.}$$

According to Sun and Walker,³⁹ the solution required is a fixed point of the operator $F$, and where such fixed point conditions exist, an iteration of the form,

Eq. (20)

$$P^{(k+1)}=F(P^{(k)})\mathrm{.}$$

should converge to the fixed point of the operator. An iteration algorithm can be rewritten as follows,

Eq. (21)

$${\hat{P}}^{(k+1)}=\frac{[H^T.(\frac{I}{K+H.{\hat{P}}^{(k)}}\left)\right].{\hat{P}}^{(k)}}{\alpha },$$

where $\hat{P}$ is the estimate of $P$. $\hat{P}$ is the profile corresponding to a maximum of $l(P)$ which corresponds to the most likely light distribution along the acoustic axis.

Also according to Sun and Walker,³⁹ a suitable starting point is to set the first estimate of the profile ${\hat{P}}^{(1)}$ to be a constant so as not to introduce any prior knowledge into the algorithm. Furthermore, when $I\cong K+H.{\hat{P}}^{(k)}$, i.e., if the estimated profile generates a result close to the real data, the bracketed term in Eq. (18) reduces to $\alpha$, and the algorithm will cease to make any further significant changes to the estimated profile.

3.

Experiment Configuration

Initially, the model is validated using clear water as the medium. The experimental set up in this case uses an expanded light beam and an aperture in front of the detector to experimentally control the range of delayed US modulated pulses in Eq. (3) In the case of a scattering medium, the beam expander is removed to increase signal level and improve spatial resolution.

3.1.

Experimental Configuration (Water)

The experimental setup is shown in Fig. 4. An expanded HeNe (632.8 nm, 20 mW) laser illuminates the sample, and a photomultiplier tube (PMT, Hamamatsu H5783-20) is employed as a detector. An adjustable slit is placed 10 mm in front of the PMT in order to experimentally control the range of optical pulses included in the summation described in Eq. (3). The distance between the PMT and the tank is 200 mm. A pulser-receiver (Olympus Panametrics-NDT 5800) was employed to drive a focused US transducer (Olympus Panametrics-NDT V322). The US transducer has a 7.6 cm focal length, a 10 MHz central frequency, 2.54 cm nominal element size, and the $-6\mathrm{dB}$ focal zone is ($0.459\times 10.4{\mathrm{mm}}^2$) ($\text{lateral size}\times \text{axial size}$).⁴⁸ The acoustic pulse is 200 ns long with a maximum pressure of 2 MPa. A $12\times 15\times 5\mathrm{cm}$ ($H\times W\times D$) water tank sits on a computer controlled $XYZ$ motorized stage (Standa 8 MT175-50). The signal from the PMT is fed into an amplifier, before going to an oscilloscope (Tektronix TDS2024B, 8-bit digitizer) and subsequent storage on a PC.

Fig. 4

Experiment setup to investigate pulsed USMOT signals in clear water. The 500 μm wide rod (object 1) is placed close to, but not within, the US focal region to avoid a mismatch in acoustic density (rod and water) which may cause a change in US profile. This does not significantly affect the convolution model in clear water.

The 500 μm wide rod (object 1) is used in the experiment to investigate the disturbance that an optically absorbing object has on the modulated signal. The rod is made of black polyamide which is assumed to be totally optically absorbing. As there is a density mismatch between the rod and water ($\text{polyamide}\sim 1.15\mathrm{g}/{\mathrm{cm}}^3$, $\text{water}\sim 1\mathrm{g}/{\mathrm{cm}}^3$), the US focus is placed at the edge of the rod.

3.2.

Experimental Configuration (Scattering Medium)

The validity of the model is then investigated in a scattering medium. The experimental setup was slightly modified at the detector as shown in Fig. 5. As noted previously^14,15 higher optical intensity in the region of the US focus will increase the peak amplitude of the detected pulse and so the beam expander is removed in this case. A 50 mm focal length plano-convex lens (placed 200 mm away from the tank) is employed to collect more light to the PMT. The number of speckles detected at the detector is approximately 10. Although this is a relatively large number of speckles to be averaged over the detector plane, it is insufficient to cancel the US modulated optical signals in the cases shown here as the scattering is relatively weak ($k_a\gg {\mu }_s$) in both scattering experiments, which means that modulation of the refractive index is the dominant mechanism and the modulation depth of the detected signals is relatively large.² For the experimental results shown in Sec. 4.2, the scattering medium was agarose gel with scatterers (1.6 μm-diam polystyrene microspheres). For the experimental results shown in Sec. 4.4 (2-D image), the scattering medium was water with scatterers (1.6 μm-diam polystyrene microspheres).

Fig. 5

(a) Experiment setup to investigate pulsed USMOT signals in (b) a scattering gel block with an embedded 500 μm black plastic rod (object 1, Fig. 4), the position of the rod is indicated by the dashed rectangle, it sits 8 mm from the top and at the mid plane of the gel, it cannot be seen in the photograph due to multiple scattering of the illumination. The gel block dimensions are $25\times 16\times 32\mathrm{mm}$ ($XYZ$).

4.

Results

The model is initially validated using the set up described in Sec. 3.1. In Sec. 4.1 and 4.2, noise-free simulated data [obtained from Eq. (10)] are presented and compared with experimental data. Section 4.3 employs simulated data with different levels of noise to test the performance of the ML reconstruction. Experimental results obtained using the ML algorithm are shown in the final section.

4.1.

Pulsed USMOT Signal (Water)

The initial experiment is to investigate the model of the pulsed USMOT signal described in Sec. 2. The simple case of a non-scattering medium, in this case water, is used to aid understanding. For the transducer properties described in Sec. 3.1, $P_{\mathrm{us}}(z)$ is fairly uniform within the region of interest (ROI) of the experiments. In addition, the expanded laser beam source ($\sim 20\mathrm{mm}$ standard deviation) also has a wide Gaussian intensity profile compared to the maximum ROI used (2 mm). Hence, parameters $L_1$ and $L_2$ in Eqs. (7) and (8) are set to be larger than the ROI to ensure that this is the case. For example, in the data shown in Fig. 6 (column 1, $L_1=20\mathrm{mm}$ and $L_2=10\mathrm{mm}$), the adjustable slit is used to restrict the region from which the modulated light is detected i.e., the limits of the summation shown in Eq. (3). The attenuation coefficient $\beta$ in Eq. (6) is set to $25\times {10}^{-15}\text{neper}{\mathrm{m}}^{-1}{\mathrm{Hz}}^{-2}$ at 25 °C in distilled water.

Fig. 6

Simulated (column 2) and experimental (column 3) results of pulsed USMOT signals for different profiles (column 1): rows 1 to 4 correspond to detector slit sizes of 200 μm, 1 mm, 2 mm (without an object) and 2 mm (with object 1). All the signals are normalized by their maximum absolute values. Column 1 shows $P(z)$ used in the model to obtain the simulated data. $z$ is defined in Fig. 5(a).

Figure 6 presents the simulated [based on Eq. (10)] and experimental results of the pulsed USMOT signals at the detector when the slit size is adjusted. The first column shows the acousto-optic profile along the acoustic axis, $P(z)$, used within the simulation. The second column shows the noise-free temporal optical signals predicted by the simulation, and the third column shows the experimental temporal optical signals. Each experimental pulse is averaged the maximum 128 times on the oscilloscope and further averaging ($\times 100$) is carried out on the PC. Rows 1 to 3 show slit sizes of 200 μm, 1 mm, 2 mm (without an object). All signals are normalized by their maximum absolute values. There is good agreement between the shape of the simulated and experimental signals.

As the slit size increases, the pulsed USMOT signal becomes more distorted because it is composed of a larger range of phase-shifted pulses. In addition, the contrast of the detected signal is lower when more phase-shifted pulses are detected, and the modulated signals tend to spread out as the slit size increases. Physically, the temporal ultrasonic pressure within the medium maps onto the optical pulse. When the profile set by the slit size is smaller than an acoustic wavelength, there is less destructive interference at the photosensor as most of the optical pulses that reach the detector have a similar phase (Fig. 6, row 1). On the contrary, when the slit size increases, the optical pulses that propagate to the photodetector have a wider range of phase delays and cancel out to produce a DC light level (Fig. 6, rows 2 and 3). There is always a modulated component at the beginning and end of the pulse because first arriving and last arriving modulation does not have a pulse with which to destructively interfere. Figure 6 (row 4) predicts the signal when a 500 μm wide optical absorber (object 1, Fig. 5) is placed in the middle of the ROI. The absorber in this case has disturbed the destructive interference effect which makes a significant change to the shape of the detected optical pulse compared to the case without the object (Fig. 6, row 3).

4.2.

Pulsed USMOT Signal (Scattering Medium)

Figure 7 compares simulation and experiment in the absence (row 1) and presence (row 2) of a 500 μm wide absorbing rod (object 1, Fig. 5). The scattering medium in this case is a $32\times 25\times 16\mathrm{mm}$ ($H\times W\times D$) block of agarose gel mixed with microspheres (scattering coefficient ${\mu }_s=10{\mathrm{cm}}^{-1}$). The first column shows the optical profile $P(z)$ along the acoustic axis $z$ used to generate the simulated temporal profile shown in the second column. In generating the data in the first column of Fig. 7, $L_1=2\mathrm{cm}$ and $L_2=1.5\mathrm{mm}$ ($L_2$ is shorter in this case compared to Fig. 6 as the beam expander has been removed). The simulated signals are generated from an estimated Gaussian optical profile which is clearly approximate as the true optical profile $P(z)$ along the acoustic axis is unknown. However, the results again show a good agreement between the shape of the simulated and experimental detected pulse. This suggests that the proposed model is appropriate for forming the basis of the maximum likelihood algorithm.

Fig. 7

Simulated (column 2) and experimental (column 3) results of pulsed USMOT signals for different profiles (column 1): Without an optical absorber (row 1) and with the 500 μm optical absorber (object 1, row 2). The vertical axes show the signals normalized by their maximum absolute values. Column 1 shows $P(z)$ used in the model to obtain the simulated data.

4.3.

Application of ML to Simulated Pulsed USMOT Signals

The performance of the ML reconstruction method is investigated in this section using simulated data with different noise levels. The SNR in this case is defined as the ratio of the signal energy and the noise energy,

Eq. (22)

$$\mathrm{SNR}=\frac{\sum _i{I}_i^2}{\sum _i{N}_i^2},$$

where $I_i$ and $N_i$ are the amplitudes of the noise-free signal and the imposed noise, respectively. We also introduce an error coefficient to evaluate the reconstruction quality. The error coefficient is defined as,

Eq. (23)

$$\text{Error}=\sqrt{\frac{\sum _i{({\hat{P}}_i-P_i)}^2}{K},}$$

where and $\widehat{P_i}$ and $P_i$ are reconstructed and original profiles, respectively, $K$ is the number of points in the profile.

Two profiles indicating the presence of a single object and two closely spaced objects are used. Figure 8 illustrates the ML image reconstruction of a 500 μm optical absorber placed in the middle of the ROI. In the simulation, the object is assumed to have the same acoustic properties to that of the surrounding medium and the optical attenuation factor $G=0.7$. The number of ML iterations is ${10}^4$; the simulation is repeated 10 times with different randomized noise to produce the error bars in Fig. 8(a). Figure 8(a) shows the error coefficient [Eq. (23)] as a function of SNR. As anticipated, as SNR increases, the error is reduced. Figures 8(b), 8(d), and 8(f) show the simulated pulse with noise added (SNR values of 6.3, 0.8 and 0.1, respectively). Figures 8(c), 8(e), and 8(g) show the original and reconstructed profiles.

Fig. 8

ML reconstruction for a profile with the 500 μm object: (a) Error-SNR plot, Error factor is in arbitrary unit. Simulated optical pulses for different SNRs (b) 6.3, (d) 0.8, (f) 0.1. Reconstructed (dotted) and original (solid) profiles for SNRs (c) 6.3, (e) 0.8, (g) 0.1. The vertical axes show the detected signals normalized by their maximum absolute values.

Figure 9 shows similar results for two closely spaced 250 μm-diam absorbing rods (acoustic properties are matched to the surrounding medium) with different attenuation factors ($G=0.7$ and 0.5). Similarly, ML reconstructs successfully the profile when the SNR is greater than 1. When the SNR is less than 1, ML fits to the noise rather than the signal and the reconstructed profile becomes very noisy.

Fig. 9

ML reconstruction for a profile with two 250 μm wide objects with different optical attenuation factors: (a) Error-SNR plot, the error factor is in arbitrary unit. Simulated optical pulses for different SNRs (b) 8.7, (d) 1.1, (f) 0.1. Reconstructed (dotted) and original (solid) profiles for SNRs (c) 8.7, (e) 1.1, (g) 0.1. Signals are normalized by their maximum absolute values.

4.4.

Application of ML to Experimental Pulsed USMOT Signals

4.4.1.

1-D image results

In this section, the ML data inversion method is applied to the two previous experimental results shown in Fig. 8. Reconstruction is based on the iterative algorithm shown in Eq. (21). The transfer matrix $H$ is constructed by using Eq. (12) with $\mathrm{\Delta }t=25\mathrm{ns}$, corresponding to a 37.5 μm layer in the medium. This is a sufficient resolution as the 10 MHz acoustic wavelength is 150 μm, corresponding to a 150 μm layer in the medium. A shorter $\mathrm{\Delta }t$ will significantly increase the processing time, but does not improve the result.

Figure 10 shows the results after ${10}^4$ iterations of Eq. (21) for a single 500 μm-diam absorbing rod. Figures 10(a) and 10(c) show the detected temporal signals in the absence and presence of the rod, respectively. Figures 10(b) and 10(d) show the reconstructed profiles from the ML algorithm and demonstrate the absence and presence of an absorbing object, respectively.

Fig. 10

ML method for pulsed USMOT image reconstruction in a scattering gel: (a) temporal signal in the absence of an object; (b) reconstructed profile due to the signal from (a); (c) temporal signal in the presence of object 1; (d) reconstructed profile due to the signal from (c).

4.4.2.

2-D image results

The scattering medium used in this case is distilled water with a suspension of scatterers (1.6 μm-diam polystyrene microspheres, ${\mu }_s=6.5{\mathrm{cm}}^{-1}$, $g=0.93$). Although the scattering coefficient in this case is relatively low ($k_a\gg {\mu }_s$) there is no significant ballistic light component. Calculating the attenuation of unscattered light as $\exp (\text{-}{\mu }_{s}d)$ where $d$ is the depth within the medium the proportion of unscattered light is 8.8e−8 at the object plane and 7.7e−15 at the exit plane. This does not produce an observable change in detected signal as the acoustic pulse crosses the optical axis. The optical absorbing object is a 400 μm wide polyamide strip (object 2) which is painted black (assumed to be totally absorbing) at certain positions (1.5 mm absorber, 1.5 mm transparent) and is placed in the middle of a $12\times 15\times 5\mathrm{cm}$ ($H\times W\times D$) tank. The thickness of the strip is $\sim 150\mu \mathrm{m}$. Figure 11 shows the object lying on a ruler (a), and a photograph of the object in a scattering medium contained in a tank (b). In the experiment, the laser and US transducer are stationary. As the US pulse scans along the US axis, it produces a 1-D trace. The stage then scans (0.2 mm scanning step) the water tank laterally to produce a 2-D image of the object. Figure 12 presents the image of the object shown in Fig. 11 using the USMOT system with the ML reconstruction method. The result in Fig. 12 shows only a part ($\sim 4$ transparent and $\sim 3$ absorber) of the object, which requires 40 scans. The total experiment time including image reconstruction (performed as soon as a single data is acquired) is 80 min. This is because SNR is low; hence a long temporal average ($100\times 128$ times) is required to obtain a signal that can be used in reconstruction.

Fig. 11

(a) Optical absorbing object (object 2), and (b) in scattering medium contained in a tank.

Fig. 12

ML method for the USMOT signal from object 2: (a) 2-D image, (b) cross section 1 ($z$-axis), (c) cross section 2 ($x$-axis).

The image provided in Fig. 12(a) from the USMOT system shows a good agreement with the real object. The absorber size is $500\mu \mathrm{m}\times 1.4\mathrm{mm}$ which is similar to the transparent region. When the US column is scanned laterally, only the optical property changes, the acoustic property remains the same. Higher intensity appears when the transparent region is in place, and an intensity drop confirms the presence of an optical absorber. Figures 12(b) and 12(c) present the lateral and axial optical intensity profile. In this experiment, the object and scattering medium have different acoustic properties ($\text{Polyamide}\text{density}=1.15\mathrm{g}/{\mathrm{cm}}^3$, $\text{water}\text{density}=1\mathrm{g}/{\mathrm{cm}}^3$). However, as the object is scanned laterally, the acoustic profile remains relatively unchanged at each measurement. Therefore, the change in the signal is due to the difference in optical properties (absorption and scattering). The absorbing coefficient is larger at the black region, while the scattering coefficient is reduced at the transparent region on the strip. As a result, the low intensity region in Fig. 12(a) corresponds to the black absorber, and the high intensity region is the transparent area on the strip.

We estimate the resolution of the system by using edge spread function (distance it takes for the intensity to change from 10% to 90% of the maximum) to be 600 μm lateral and 160 μm axial. The pure acoustic system using a 10 MHz US transducer has a 459 μm beam waist and 150 μm wavelength.

5.

Discussion

The main aim of this work has been to investigate the use of a ML algorithm for reconstructing USMOT data. To carry out a ML inversion it is necessary to obtain a model for the received data. In this case, the ML reconstruction algorithm is based on a relatively simple model of US and light in tissue. Although the model is simple, experimental results (Figs. 6 and 7) have demonstrated that the trends of USMOT experiments in terms of optical pulse profiles can be predicted in the presence and absence of absorbing objects. The model can also predict the experimental observations of Masaki et al.⁴⁹ where in the absence of an object the optical pulses that propagate to a detector cancel out. In the presence of an object the pulses are perturbed and an optical pulse is detected (simulated results not shown in this paper).

The model neglects the effects that diffraction and laser speckle have on the Gaussian intensity optical profile $P_{\mathrm{opt}}(z)$ along the ROI, however, this is reasonable in the case when there is sufficient spatial or temporal averaging of the speckle pattern or if an US modulated fluorescence signal is detected. The model also assumes a Gaussian profile for the acoustic and optical profiles in the medium, and that the optical profile and absorption profiles can be treated independently [Eq. (4)] which is again approximate. The estimate of the optical profile could be improved by using a more accurate model of light propagation in tissue such as a Monte Carlo or diffusion approximation. Similarly, the axial acoustic pressure can be mapped by using one of many available online simulation tools,⁵⁰ but a more practical way is to characterize the transducer by scanning a hydrophone along the US travelling axis. Various acousto-optic models, both analytical and numerical, have been proposed,^{2,22–25,27–30,51} and could in future work be adopted as the basis for a ML algorithm. This would more accurately model the acousto-optic interactions and provide a better representation of the detected modulated optical signals. For example, a change in the scattering coefficient within a ROI alters both the intensity optical profile and the modulation depth.^22–24,51 Although this could, to a certain extent, be incorporated into the current model by replacing $M$ with a variable $M(z)$ that depends on scattering, it is more accurately dealt with using a more sophisticated model. The model also assumes that the thickness of the absorbing object is infinitesimal and that absorption can be approximated by an attenuation factor $G$ [Eq. (9)]. This could also be refined by a more sophisticated model that takes into account the 3-D nature of absorbing objects with a path-length dependent absorption coefficient ${\mu }_a$. These improvements to the model used as a basis for ML are likely to provide more accurate reconstruction, particularly for geometries that are more complex than those considered here.

Despite these limitations, the model predicts the trends of USMOT experiments and therefore provided a useful basis for initial testing of the performance of the ML algorithm. With noise added to the simulated profiles ML was able to accurately reconstruct the profiles of a single absorber and two absorbers for $\mathrm{SNR}>1$ (Figs. 8 and 9). The ML method has also been shown to be capable of reconstructing experimental profiles $P(z)$ for the case of no object, a single object, and multiple objects. At present the ML algorithm based on Eq. (21) allows us to obtain the profile $P(z)$ which comprises the optical, acoustic, and absorption profiles [Eq. (4)]. More accurate models and deconvolution using a known acoustic profile will further improve the absorption profiles obtained. In an acoustically homogeneous and optically inhomogeneous medium, the problem of obtaining the absorption profile is easier as all the changes occurring in the signal are due to the optical properties. In a medium that is both acoustically and optically inhomogeneous, the method is still capable of obtaining the optical properties of the sample as long as the acoustic properties are known. In this case, purely acoustic imaging can be employed to map the tissue density which can be used as prior knowledge in reconstructing optical properties of a sample via USMOT measurements.

The ML approach could be extended to 3-D reconstruction by making measurements at multiple source and detector positions. A simple extension would be to scan the US across the scattering medium and apply the reconstruction approach described in this paper (Fig. 13) based on multiple line-scans. However, in a more complex case where several optical absorbing or scattering objects are located within a sample, the reconstruction will clearly become more challenging as is the case for most diffuse optical techniques. For example, in a simple configuration where two optically absorbing objects lie on the optical axis ($y$ axis), both objects influence the modulated signal obtained at both positions (Fig. 13).

Fig. 13

Experiment configuration with 2 objects on optical axis.

One approach to extend the line scan approach described in this paper might be to use information obtained from one line scan to predict the optical distribution $P_{\mathrm{opt}}(z)$ at another. For example (Fig. 13), initially if the transducer is placed at scan position (1) and measurements were made at a detector positioned close to the light source, the effect of object 1 will dominate the detected signal. This information could then be used to estimate the optical profile at the position of the second object. However, the most likely route to full 3-D reconstruction would be to develop a more sophisticated 3-D model and apply the ML algorithm to measurements made at multiple source-detector positions. This will depend on how well conditioned the problem is and will require further investigation beyond the scope of this paper.

The main advantage of the ML method is its ability to take into account the random nature of noise. In other words, it recognizes that the recorded data are the realization of a random process.^36–39 Besides, the ML method does not require a priori information to be imposed on the solution. The ML method can be time consuming to achieve the required accuracy. The number of iterations required to obtain a fixed accuracy may differ for different experimental conditions. It is difficult to specify an optimized value for the number of iterations for a general case. However, the ML algorithm always converges to the “most-likely” results after a number of iterations (usually around ${10}^4$ to $5\times {10}^4$ times, approximately 2 min with a 3.4 GHz Pentium4 processor, 1 GB Ram). It then reaches the optimum point before the algorithm starts to fit the noise. An evaluation criterion can be applied to set the stopping point of the method. One of the suggested criteria is the goodness-of-fit $\chi$ (Ref. ³⁹) which is defined as,

Processing speed is important depending on the application area. If real time measurements of tissue are required then a fast processing algorithm is necessary and this will become an even greater problem if the method is extended to 3-D reconstruction. Our main application area is in regenerative medicine and the growth of tissue within scaffolds. In this case, processing speed is less of an issue as the growth of tissue within a scaffold is slow. The reconstruction speed obtained in this work is achieved with a standard m-file MATLAB program which is not optimized for speed. Improvements in coding and hardware could alleviate this issue.

6.

Conclusion

A simple model of a pulsed USMOT signal has been presented which provides a good approximation to experimental pulsed USMOT signals. The model has been used as the basis for a ML reconstruction algorithm. Reconstruction of simulated data in the presence of noise has demonstrated that the algorithm performs well provided the signal to noise ratio is greater than unity. The inversion has also been applied to experimental data and it has been shown to be capable of reconstructing multiple objects in a 2-D imaging experiment. The resolution of the system with a 10 MHz US transducer is estimated to be 160 μm axially and 600 μm laterally.