Optical imaging of soft biological tissues is highly desirable since it is nonionizing and provides sensitive contrast information which enables detection of physiological functions and abnormalities, including potentially early cancer detection. However, due to the diffusion of light, it is dificult to achieve simultaneously both good spatial resolution and good imaging depth with the pure optical imaging modalities. This work focuses on the ultrasound-modulated optical tomography - a hybrid technique which combines advantages of ultrasonic resolution and optical contrast. In this technique, focused ultrasound and optical radiation of high temporal co-herence are simultaneously applied to soft biological tissue, and the intensity of the ultrasound-modulated light is measured. This provides information about the optical properties of the tissue, spatially localized at the interaction region of the ultrasonic and electromagnetic waves. In experimental part of this work we present a novel implementation of high-resolution ultrasound-modulated optical tomography that, based on optical contrast, can image several millimeters deep into soft biological tissues. A long-cavity confocal Fabry-Perot interferometer was used to detect the ultrasound-modulated coherent light that traversed the scattering biological tissue. Using 15-MHz ultrasound, we imaged with high contrast light absorbing structures placed 3 mm below the surface of chicken breast tissue. The resolution along the axial and the lateral directions with respect to the ultrasound propagation direction was better than 70 and 120¹m, respectively. This technology is complementary to other imaging technologies, such as confocal microscopy and optical-coherence tomography, and has potential for broad biomedical applications. In the theoretical part we present various methods to model interaction be-tween the ultrasonic and electromagnetic waves in optically scattering media. We first extend the existing theoretical model based on the diffusing-wave spectroscopy approach to account for anisotropic optical scattering, Brownian motion, pulsed ul-trasound, and strong correlations between the ultrasound-induced optical phase in-crements. Based on the Bethe-Salpeter equation, we further develop a more general correlation transfer equation, and subsequently a correlation diffusion equation, for ultrasound-modulated multiply scattered light. We expect these equations to be applicable to a wide spectrum of conditions in the ultrasound-modulated optical tomography of soft biological tissues.

Optical imaging of soft biological tissue is highly desirable since it is nonionizing and provides sensitive contrast information which enables the detection of physiological functions and abnormalities, including potentially early cancer detection. However, due to the diffusive nature of light in soft biological tissue, it is difficult to achieve simultaneously good spatial resolution and good imaging depth with pure optical imaging modalities. This work focuses on the ultrasound-modulated optical tomography (UOT): a hybrid technique which combines the advantages of ultrasonic resolution and optical contrast. In this technique, focused ultrasound and optical radiation of high temporal coherence are simultaneously applied to soft biological tissue. The intensity of the sideband, or ultrasound ‗tagged‘ photons depends on the optical absorption in the region of interest where the ultrasound is focused. Demodulation of the optical speckle pattern yields the intensity of tagged photons for each location of the ultrasonic focal spot. Thus UOT yields an image with spatial resolution of the focused ultrasound — typically submillimeter — whose contrast is related to local optical absorption and the diffusive properties of light in the organ. Thus it extends all the advantages of optical imaging deep into highly scattering tissue. However lack of efficient tagged light detection techniques has so far prevented ultrasound-modulated optical tomography from achieving maturity. The signal-to-noise ratio (SNR) and imaging speed are two of the most important figures of merit and need further improvement for UOT to become widely applicable. In the first part of this work, nonlinear optics detection methods have been implemented to demodulate the ―tagged‖ photons. The most common of these is photorefractive (PR) two wave mixing (TWM) interferometry, which is a time-domain filtering technique. When used for UOT, it is found that this approach extracts not only optical properties but also mechanical properties for the area of interest. To improve on TWM, PR four wave mixing (FWM) experiments were performed to read out only the modulated light and at the same time strongly suppressing the ‗untagged‘ light. Spectral-hole burning (SHB) in a rare-earth-ion-doped crystal has been developed for UOT more recently. Experiments in Tm3 :Y3Al5O12 (Tm:YAG) show the outstanding features of SHB: large angle acceptance (etendue), light speckle processing in parallel (insensitive to the diffusive light nature) and real-time signal collection (immune to light speckle decorrelation). With the help of advanced laser stabilization techniques, two orders of magnitude improvement of SNR have been achieved in a persistent SHB material (Pr^3 :Y2SiO5) compared to Tm:YAG. Also slow light with PSHB further reduces noise in Pr:YSO UOT that is caused by polarization leakage by performing time-domain filtering.

In this dissertation, two recent problems from tomographic imaging are studied, and results from numerical simulations with synthetic data are presented. The first part deals with ultrasound modulated optical tomography, a method for imaging interior optical properties of partially translucent media that combines optical contrast with ultrasound resolution. The primary application is the optical imaging of soft tissue, for which scattering and absorption rates contain important functional and structural information about the physiological state of tissue cells. We developed a mathematical model based on the diffusion approximation for photon propagation in highly scattering media. Simple reconstruction schemes for recovering optical absorption rates from boundary measurements with focused ultrasound are presented. We show numerical reconstructions from synthetic data generated for mathematical absorption phantoms. The results indicate that high resolution imaging with quantitatively correct values of absorption is possible. Synthetic focusing techniques are suggested that allow reconstruction from measurements with certain types of non-focused ultrasound signals. A preliminary stability analysis for a linearized model is given that provides an initial explanation for the observed stability of reconstruction. In the second part, backprojection schemes are proposed for the detection of small amounts of highly enriched nuclear material inside 3D volumes. These schemes rely on the geometrically singular structure that small radioactive sources represent, compared to natural background radiation. The details of the detection problem are explained, and two types of measurements, collimated and Compton-type measurements, are discussed. Computationally, we implemented backprojection by counting the number of particle trajectories intersecting each voxel of a regular rectangular grid covering the domain of detection. For collimated measurements, we derived confidence estimates indicating when voxel trajectory counts are deviating significantly from what is expected from background radiation. Monte Carlo simulations of random background radiation confirm the estimated confidence values. Numerical results for backprojection applied to synthetic measurements are shown that indicate that small sources can be detected for signal-to-noise ratios as low as 0.1%.

Stable and computationally efficient reconstruction methodologies are developed to solve two important medical imaging problems which use near-infrared (NIR) light as the source of interrogation, namely, diffuse optical tomography (DOT) and one of its variations, ultrasound-modulated optical tomography (UMOT). Since in both these imaging modalities the system matrices are ill-conditioned owing to insufficient and noisy data, the emphasis in this work is to develop robust stochastic filtering algorithms which can handle measurement noise and also account for inaccuracies in forward models through an appropriate assignment of a process noise. However, we start with demonstration of speeding of a Gauss-Newton (GN) algorithm for DOT so that a video-rate reconstruction from data recorded on a CCD camera is rendered feasible. Towards this, a computationally efficient linear iterative scheme is proposed to invert the normal equation of a Gauss-Newton scheme in the context of recovery of absorption coefficient distribution from DOT data, which involved the singular value decomposition (SVD) of the Jacobian matrix appearing in the update equation. This has sufficiently speeded up the inversion that a video rate recovery of time evolving absorption coefficient distribution is demonstrated from experimental data. The SVD-based algorithm has made the number of operations in image reconstruction to be rather than. 2()ONN3()ONN The rest of the algorithms are based on different forms of stochastic filtering wherein we arrive at a mean-square estimate of the parameters through computing their joint probability distributions conditioned on the measurement up to the current instant. Under this, the first algorithm developed uses a Bootstrap particle filter which also uses a quasi-Newton direction within. Since keeping track of the Newton direction necessitates repetitive computation of the Jacobian, for all particle locations and for all time steps, to make the recovery computationally feasible, we devised a faster update of the Jacobian. It is demonstrated, through analytical reasoning and numerical simulations, that the proposed scheme, not only accelerates convergence but also yields substantially reduced sample variance in the estimates vis-à-vis the conventional BS filter. Both accelerated convergence and reduced sample variance in the estimates are demonstrated in DOT optical parameter recovery using simulated and experimental data. In the next demonstration a derivative free variant of the pseudo-dynamic ensemble Kalman filter (PD-EnKF) is developed for DOT wherein the size of the unknown parameter is reduced by representing of the inhomogeneities through simple geometrical shapes. Also the optical parameter fields within the inhomogeneities are approximated via an expansion based on the circular harmonics (CH) (Fourier basis functions). The EnKF is then used to recover the coefficients in the expansion with both simulated and experimentally obtained photon fluence data on phantoms with inhomogeneous inclusions. The process and measurement equations in the Pseudo-Dynamic EnKF (PD-EnKF) presently yield a parsimonious representation of the filter variables, which consist of only the Fourier coefficients and the constant scalar parameter value within the inclusion. Using fictitious, low-intensity Wiener noise processes in suitably constructed ‘measurement’ equations, the filter variables are treated as pseudo-stochastic processes so that their recovery within a stochastic filtering framework is made possible. In our numerical simulations we have considered both elliptical inclusions (two inhomogeneities) and those with more complex shapes ( such as an annular ring and a dumbbell) in 2-D objects which are cross-sections of a cylinder with background absorption and (reduced) scattering coefficient chosen as = 0.01 mm-1 and = 1.0 mm-1respectively. We also assume=0.02 mm-1 within the inhomogeneity (for the single inhomogeneity case) and=0.02 and 0.03 mm-1 (for the two inhomogeneities case). The reconstruction results by the PD-EnKF are shown to be consistently superior to those through a deterministic and explicitly regularized Gauss-Newton algorithm. We have also estimated the unknown from experimentally gathered fluence data and verified the reconstruction by matching the experimental data with the computed one. The superiority of a modified version of the PD-EnKF, which uses an ensemble square root filter, is also demonstrated in the context of UMOT by recovering the distribution of mean-squared amplitude of vibration, related to the Young’s modulus, in the ultrasound focal volume. Since the ability of a coherent light probe to pick-up the overall optical path-length change is limited to modulo an optical wavelength, the individual displacements suffered owing to the US forcing should be very small, say within a few angstroms. The sensitivity of modulation depth to changes in these small displacements could be very small, especially when the ROI is far removed from the source and detector. The contrast recovery of the unknown distribution in such cases could be seriously impaired whilst using a quasi-Newton scheme (e.g. the GN scheme) which crucially makes use of the derivative information. The derivative-free gain-based Monte Carlo filter not only remedies this deficiency, but also provides a regularization insensitive and computationally competitive alternative to the GN scheme. The inherent ability of a stochastic filter in accommodating the model error owing to a diffusion approximation of the correlation transport may be cited as an added advantage in the context of the UMOT inverse problem. Finally to speed up forward solve of the partial differential equation (PDE) modeling photon transport in the context of UMOT for which the PDE has time as a parameter, a spectral decomposition of the PDE operator is demonstrated. This allows the computation of the time dependent forward solution in terms of the eigen functions of the PDE operator which has speeded up the forward solution, which in turn has rendered the UMOT parameter recovery computationally efficient.

This thesis describes the development and applications of an extension of DWS which enables the recovery of ‘localized’ mechanical properties, in a specified region of a complex jelly-like object which is inhomogeneous, marked out by the focal volume of an ultrasound transducer, also called the region-of-interest (ROI). Introduction of the sinusoidal forcing creates a sinusoidal phase variation in the detected light in a DWS experiment which modulates the measured intensity autocorrelation, g2 (τ ). Decay in the modulation depth with τ is used to recover the visco-elastic spectrum of the material in the ROI. En route to this, growth of the mean-squared dis- placement (MSD) with time is extracted from the modulation depth decay, which was verified first by the usual DWS experimental data from an homogeneous object with properties matching those in the ROI of the inhomogeneous object and then those obtained by solving the generalized Langevin equation (GLE) modelling the dynamics of a typical scattering centre in the ROI. A region-specific visco-elastic spectral map was obtained by scanning the inhomogeneous object by the ultrasound focal volume. Further, the resonant modes of the vibrating ROI were measured by locating the peaks of the modulation depth variation in g2(τ ) with respect to the ultrasound frequency. These resonant modes were made use of to recover elasticity of the material of the object in the ROI. Using a similar strategy, it was also shown that flow in pipe can be detected and flow rate computed by ‘tagging’ the photons passing through the pipe with a focussed ultrasound beam. It is demonstrated, both through experiments and simulations that the ultrasound-assisted technique devel- oped is better suited to both detect and quantitatively assess flow in a background of Brownian dynamics than the usual DWS. In particular, the MSD of particles in the flow, which shows forth a super-diffusive dynamics with MSD growing following τ α with α < 2, is captured over larger intervals of τ than was possible using existing methods. On the theoretical front, the main contribution is the derivation of the GLE, with multiplicative noise modulating the interaction ‘spring constant’. The noise is derived as an average effect of the micropolar rotations suffered by the ‘bath’ particles on the ‘system’ particle modelled. It has been shown that the ‘local’ dynamics of the system particle is nontrivially influenced by the dynamics, both translation and rotation, of ‘nonlocal’ bath particles.

The focus of the thesis is on the development of a few stochastic search schemes for inverse problems and their applications in medical imaging. After the introduction in Chapter 1 that motivates and puts in perspective the work done in later chapters, the main body of the thesis may be viewed as composed of two parts: while the first part concerns the development of stochastic search algorithms for inverse problems (Chapters 2 and 3), the second part elucidates on the applicability of search schemes to inverse problems of interest in tomographic imaging (Chapters 4 and 5). The chapter-wise contributions of the thesis are summarized below. Chapter 2 proposes a Monte Carlo stochastic filtering algorithm for the recursive estimation of diffusive processes in linear/nonlinear dynamical systems that modulate the instantaneous rates of Poisson measurements. The same scheme is applicable when the set of partial and noisy measurements are of a diffusive nature. A key aspect of our development here is the filter-update scheme, derived from an ensemble approximation of the time-discretized nonlinear Kushner Stratonovich equation, that is modified to account for Poisson-type measurements. Specifically, the additive update through a gain-like correction term, empirically approximated from the innovation integral in the filtering equation, eliminates the problem of particle collapse encountered in many conventional particle filters that adopt weight-based updates. Through a few numerical demonstrations, the versatility of the proposed filter is brought forth, first with application to filtering problems with diffusive or Poisson-type measurements and then to an automatic control problem wherein the exterminations of the associated cost functional is achieved simply by an appropriate redefinition of the innovation process. The aim of one of the numerical examples in Chapter 2 is to minimize the structural response of a duffing oscillator under external forcing. We pose this problem of active control within a filtering framework wherein the goal is to estimate the control force that minimizes an appropriately chosen performance index. We employ the proposed filtering algorithm to estimate the control force and the oscillator displacements and velocities that are minimized as a result of the application of the control force. While Fig. 1 shows the time histories of the uncontrolled and controlled displacements and velocities of the oscillator, a plot of the estimated control force against the external force applied is given in Fig. 2. (a) (b) Fig. 1. A plot of the time histories of the uncontrolled and controlled (a) displacements and (b) velocities. Fig. 2. A plot of the time histories of the external force and the estimated control force Stochastic filtering, despite its numerous applications, amounts only to a directed search and is best suited for inverse problems and optimization problems with unimodal solutions. In view of general optimization problems involving multimodal objective functions with a priori unknown optima, filtering, similar to a regularized Gauss-Newton (GN) method, may only serve as a local (or quasi-local) search. In Chapter 3, therefore, we propose a stochastic search (SS) scheme that whilst maintaining the basic structure of a filtered martingale problem, also incorporates randomization techniques such as scrambling and blending, which are meant to aid in avoiding the so-called local traps. The key contribution of this chapter is the introduction of yet another technique, termed as the state space splitting (3S) which is a paradigm based on the principle of divide-and-conquer. The 3S technique, incorporated within the optimization scheme, offers a better assimilation of measurements and is found to outperform filtering in the context of quantitative photoacoustic tomography (PAT) to recover the optical absorption field from sparsely available PAT data using a bare minimum ensemble. Other than that, the proposed scheme is numerically shown to be better than or at least as good as CMA-ES (covariance matrix adaptation evolution strategies), one of the best performing optimization schemes in minimizing a set of benchmark functions. Table 1 gives the comparative performance of the proposed scheme and CMA-ES in minimizing a set of 40-dimensional functions (F1-F20), all of which have their global minimum at 0, using an ensemble size of 20. Here, 10 5 is the tolerance limit to be attained for the objective function value and MAX is the maximum number of iterations permissible to the optimization scheme to arrive at the global minimum. Table 1. Performance of the SS scheme and Chapter 4 gathers numerical and experimental evidence to support our conjecture in the previous chapters that even a quasi-local search (afforded, for instance, by the filtered martingale problem) is generally superior to a regularized GN method in solving inverse problems. Specifically, in this chapter, we solve the inverse problems of ultrasound modulated optical tomography (UMOT) and diffraction tomography (DT). In UMOT, we perform a spatially resolved recovery of the mean-squared displacements, p r of the scattering centres in a diffusive object by measuring the modulation depth in the decaying autocorrelation of the incident coherent light. This modulation is induced by the input ultrasound focussed to a specific region referred to as the region of interest (ROI) in the object. Since the ultrasound-induced displacements are a measure of the material stiffness, in principle, UMOT can be applied for the early diagnosis of cancer in soft tissues. In DT, on the other hand, we recover the real refractive index distribution, n r of an optical fiber from experimentally acquired transmitted intensity of light traversing through it. In both cases, the filtering step encoded within the optimization scheme recovers superior reconstruction images vis-à-vis the GN method in terms of quantitative accuracies. Fig. 3 gives a comparative cross-sectional plot through the centre of the reference and reconstructed p r images in UMOT when the ROI is at the centre of the object. Here, the anomaly is presented as an increase in the displacements and is at the centre of the ROI. Fig. 4 shows the comparative cross-sectional plot of the reference and reconstructed refractive index distributions, n r of the optical fiber in DT. Fig. 3. Cross-sectional plot through the center of the reference and reconstructed p r images. Fig. 4. Cross-sectional plot through the center of the reference and reconstructed n r distributions. In Chapter 5, the SS scheme is applied to our main application, viz. photoacoustic tomography (PAT) for the recovery of the absorbed energy map, the optical absorption coefficient and the chromophore concentrations in soft tissues. Nevertheless, the main contribution of this chapter is to provide a single-step method for the recovery of the optical absorption field from both simulated and experimental time-domain PAT data. A single-step direct recovery is shown to yield better reconstruction than the generally adopted two-step method for quantitative PAT. Such a quantitative reconstruction maybe converted to a functional image through a linear map. Alternatively, one could also perform a one-step recovery of the chromophore concentrations from the boundary pressure, as shown using simulated data in this chapter. Being a Monte Carlo scheme, the SS scheme is highly parallelizable and the availability of such a machine-ready inversion scheme should finally enable PAT to emerge as a clinical tool in medical diagnostics. Given below in Fig. 5 is a comparison of the optical absorption map of the Shepp-Logan phantom with the reconstruction obtained as a result of a direct (1-step) recovery. Fig. 5. The (a) exact and (b) reconstructed optical absorption maps of the Shepp-Logan phantom. The x- and y-axes are in m and the colormap is in mm-1. Chapter 6 concludes the work with a brief summary of the results obtained and suggestions for future exploration of some of the schemes and applications described in this thesis.