ESTIMATING OPTICAL PROPERTIES OF SURGICAL TISSUE

Abstract

A method of assessing optical properties of surgical tissue by observing the tissue thermal response to laser irradiation includes exposing the tissue in a surgical region to a short laser pulse for triggering a localized temperature increase, and observing the triggered temperature increase by a thermal sensor. Optical properties of the tissue are then estimated, for determining laser surgery parameters for an eminent procedure, based on the observed temperature increase. In an example configuration, the thermal sensor is a thermal infrared camera, and estimation includes using an Ensemble Kalman Filter (EnKF) for comparing the temperature sensor data with the output of a computational laser-tissue interaction model.

Claims

1. A method for predicting tissue properties indicative of a thermal response of irradiated tissue, comprising: directing an irradiation signal at tissue for a therapeutic effect; receiving, from an optical sensor, thermal images indicative of a temperature response of the irradiated tissue based on optical properties of the irradiated tissue; applying the temperature response to a thermal laser-tissue interaction model for computing irradiation response coefficients indicative of an ablative effect on the irradiated tissue from the irradiation; and controlling the irradiation signal based on the irradiation response coefficients.

2. The method of claim 1 wherein the tissue properties include: an absorption coefficient indicative of a fraction of photonic energy transferred to the tissue in the form of heat; and and a scattering coefficient indicative of a direction change of a photon upon transferring the fraction of photonic energy.

3. The method of claim 1 wherein the ablative effect causes heating for inducing removal or cutting of the tissue.

4. The method of claim 1 wherein the irradiation signal is a laser and directing the irradiation signal includes actuating the laser for aiming and traversing the laser across the tissue.

5. The method of claim 4 further comprising receiving thermal images from an infrared camera focused on an irradiated region of the tissue defined by the laser.

6. The method of claim 1 wherein the thermal laser-tissue interaction model is a filter for computing optical penetration of the irradiation signal into the tissue and a fraction of energy of the irradiation signal absorbed by the irradiated tissue and manifested as heat by the irradiated tissue receiving the irradiation signal.

7. The method of claim 6 further comprising: determining an initial value for one or more of the optical properties of the tissue; computing the irradiation response coefficients based on the one or more initial values; and iteratively revising the computed irradiation response coefficients based on the received thermal images.

8. The method of claim 7 wherein iteratively revising the irradiation response coefficients includes comparing the temperature response based on the received thermal images with the computed irradiation response coefficients from a previous iteration.

9. The method of claim 1 wherein the filter performs tracking of time varying parameters based on a progression of the irradiation response coefficients in response to changes in the optical properties of the irradiated tissue resulting from a pulsing of the irradiation signal.

10. A computer guided laser device for performing surgical manipulations based on method for predicting tissue properties indicative of a thermal response of irradiated tissue, comprising: a laser operative for directing an irradiation signal at tissue for a therapeutic effect; an infrared camera for transmitting thermal images indicative of a temperature response of the irradiated tissue based on optical properties of the irradiated tissue; a computation circuit for applying the temperature response to a thermal laser-tissue interaction model for computing irradiation response coefficients indicative of an ablative effect on the irradiated tissue from the irradiation; a memory for storing the thermal laser-tissue interaction model; and a laser control responsive to the computation circuit for controlling the irradiation signal based on the irradiation response coefficients.

11. The method of claim 10 wherein the tissue properties include: an absorption coefficient indicative of a fraction of photonic energy transferred to the tissue in the form of heat; and and a scattering coefficient indicative of a direction change of a photon upon transferring the fraction of photonic energy.

12. The method of claim 10 wherein the ablative effect causes heating for inducing removal or cutting of the tissue.

13. The method of claim 10 wherein the irradiation signal is a laser beam and the laser control is configured to actuate the laser for aiming and traversing the laser across the tissue.

14. The method of claim 13 wherein the thermal images contain temperature information from an irradiated region of the tissue defined by the laser.

15. The method of claim 10 wherein the thermal laser-tissue interaction model is a filter for computing optical penetration of the irradiation signal into the tissue and a fraction of energy of the irradiation signal absorbed by the irradiated tissue and manifested as heat by the irradiated tissue receiving the irradiation signal.

16. The method of claim 10 wherein the filter is configured for tracking of time varying parameters based on a progression of the irradiation response coefficients in response to changes in the optical properties of the irradiated tissue resulting from the irradiation signal.

17. In an irradiation treatment environment for ablating tissue through controlled laser energy, a method for predicting a thermal response for absorption and scattering resulting from projected light, comprising: identifying a thermal response of an irradiated tissue receiving, from an optical sensor, thermal images indicative of a temperature response of the irradiated tissue; applying the temperature response to a thermal laser-tissue interaction model for computing irradiation response coefficients indicative of an ablative effect on the irradiated tissue; and controlling an irradiation source based on the irradiation response coefficients.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0010] The foregoing and other objects, features and advantages of the invention will be apparent from the following description of particular embodiments of the invention, as illustrated in the accompanying drawings in which like reference characters refer to the same parts throughout the different views. The drawings are not necessarily to scale, emphasis instead being placed upon illustrating the principles of the invention.

[0011] FIG. 1 is a context view of an irradiation or laser based system suitable for use with configurations herein;

[0012] FIG. 2 is a block diagram of the computation device or process of FIG. 1;

[0013] FIG. 3 shows an irradiation signal as in FIG. 1 intersecting a tissue region or sample;

[0014] FIGS. 4A-4E show graphs of constant parameter estimation using the filter of FIG. 2; and

[0015] FIGS. 5A-5E show graphs of time varying parameter estimation using the filter of FIG. 2.

DETAILED DESCRIPTION

[0016] Configurations below describe implementation of an irradiation method, system and apparatus for laser based treatment of tissue. Alternate irradiation sources and radiation targets other than tissue may be utilized with the disclosed approach. The configurations below depict a particular approach and example using surgical approaches with human tissue, and exhibit differences in laser irradiation cutting, ablation and/or removal techniques with varied tissue types. The disclosed approach is easily deployable in a robotic mechanism for controlling and guiding a robotic actuator for surgical intervention responsive to the computed irradiation parameters as described below.

[0017] FIG. 1 is a context view of an irradiation or laser based system suitable for use with configurations herein. Referring to FIG. 1, in a medical irradiation environment 100, a method for predicting tissue properties indicative of a thermal response of irradiated tissue includes directing an irradiation signal 102 at tissue 110 for a therapeutic effect. An optical sensor 115 focused on the tissue 110 receives thermal images 112 indicative of a temperature response 114 of the irradiated tissue based on optical properties of the irradiated tissue 110. A computation circuit 120 applies the received temperature response 114 to a thermal laser-tissue interaction model 150 for computing irradiation response coefficients indicative of an ablative effect on the irradiated tissue from the irradiation. The irradiation signal 102 from a laser light source 130 has a degrading effect on the tissue 110. Intensity and movement of the light source 130, such as by a laser control 132, causes ablative effects including cutting, removal and eradication of the subject tissue. Precision direction from the laser control 132 is therefore significant in ensuring the resulting therapy is administered correctly, as excessive irradiation and cause tissue charring and insufficient energy will fall short of the desired goal. Accordingly, the laser control 132 directs the irradiation signal based on the computed irradiation response coefficients.

[0018] Conventional robotic systems for laser surgery cannot monitor and control the complex interactions that occur between the laser and the tissue. Configurations depicted below propose a method to identify the physical tissue properties that govern the dynamics of laser-tissue interactions. The disclosed approach focuses on the estimation of two quantities, namely, the absorption and scattering coefficients, whose knowledge is necessary to model the optical-thermal response of laser-irradiated tissue. Giving robots the ability to identify these coefficients is key to enable the automatic control of laser actions. The claimed approach estimates absorption and scattering using a system identification approach, i.e., by measuring the tissue temperature elevation created by a prescribed set of laser inputs, such as laser aiming, spot size, power and exposure time. Other laser inputs for control and intensity of the laser may also be programmed and directed; this set is exemplary.

[0019] Conventional uses of the absorption and scattering coefficients of tissue rely on the use of specialized and/or bulky optical equipment (e.g., integrating spheres) that is cumbersome to use in a surgical setting. Configurations herein, in contrast, characterize absorption and scattering based on the observation of the tissue thermal response to laser exposure. The disclosed approach operates as follows: given an initial guess of the unknown coefficients, first attempt to predict the tissue thermal response using a thermal laser-tissue interaction model; then update the tissue optical properties to minimize the difference between the predicted and observed tissue temperature. To perform the update, an Ensemble Kalman Filter (EnKF) may be employed. The EnKF is a type of Bayesian filtering algorithm for data assimilation.

[0020] The sequential nature of EnKFs enables the implementation of an online, or continual feedback estimation process, i.e., the absorption and scattering coefficients are progressively refined as more and more temperature observations become available over time. It should be noted that the implementation of the approach in the workflow of a robotic laser procedure only imposes the introduction of a thermal sensor to monitor the tissue, such as a miniaturized infrared thermocamera. No additional input or invasive diagnostics are imposed on the laser surgery patient.

[0021] FIG. 2 is a block diagram of the computation device or process of FIG. 1. Referring to FIGS. 1 and 2, the laser control 132 directing an irradiation signal 102 at tissue 110 for a therapeutic effect from the laser light source 130. A thermal infrared camera 121 or other optical sensor 115 receives thermal images 112 indicative of a temperature response of the irradiated tissue based on optical properties of the irradiated tissue. An initial value or guess 162 is determined based on the tissue type, such as from Table I discussed further below. Other mechanisms for determining an initial value for one or more of the optical properties of the tissue may also be invoked. The computation circuit 120 applies the temperature response 114 to a thermal laser-tissue interaction model 160 for computing irradiation response coefficients indicative of an ablative effect on the irradiated tissue from the irradiation 102. The laser control 132 controls the irradiation signal 102 based on the irradiation response coefficients and other laser parameters 134 sent from the computation circuit 120.

[0022] The computation circuit 120 therefore receives thermal images 112 from the infrared camera 121 focused on an irradiated region of the tissue defined by the laser. In a particular configuration, the irradiation signal 102 is a laser and directing the irradiation signal includes actuating the laser for aiming and traversing the laser across the tissue 110.

[0023] The tissue properties include an absorption coefficient indicative of a fraction of photonic energy transferred to the tissue in the form of heat, and a scattering coefficient indicative of a direction change of a photon upon transferring the fraction of photonic energy. The resulting ablative effect causes heating for inducing removal or cutting of the tissue, including cutting, ablating, abrading and/or outright removal from the heat energy of the laser. As underscored above, careful application of the laser is needed to channel the potentially damaging effect in the proper manner to avoid damaging healthy or unexcised tissue in the surgical region.

[0024] The proposed approach identifies the tissue absorption and scattering coefficients based on Ensemble Kalman Filtering. In the EnKF framework, the unknown parameters of a system are modeled as stochastic variables whose probability density functions are encoded by a set of random realizations called an ensemble. Each time new sensor data becomes available, the ensemble is manipulated through a set of update rules to reflect the new probability distributions of the unknown parameters conditioned on the observed data. A more detailed formulation of the filter and a description of the ensemble update rules are described further below. Configurations herein use a version of EnKF capable of tracking time-varying parameters, which allows monitoring of the shift in the tissue optical properties that may occur during a laser procedure. Since these heating effects occur quickly in response to an ongoing laser surgical procedure, prompt response is important to identify appropriate laser control parameters in a timely manner.

[0025] Continuing to refer to FIG. 2, the tissue 110 is irradiated with a laser pulse or light beam 102, triggering a localized temperature increase which is observed by a thermal sensor 115 such as camera 121. For each observation, the EnKF compares the sensor data (temperature response 114 with the output of a laser-tissue interaction model 160, and iteratively updates the ensemble 164 in such a way to minimize the error between the model output and the measured tissue temperature. At any given time, the ensemble mean is used as an estimate of the unknown absorption and scattering coefficients, while the standard deviation provides a measure of uncertainty.

[0026] The thermal laser-tissue interaction model 160 defines a filter for computing optical penetration of the irradiation signal 102 into the tissue 110 and also the fraction of energy of the irradiation signal absorbed and manifested as heat by the irradiated tissue receiving the irradiation signal. It is important to define the role that the absorption and scattering coefficients play in the thermal response of laser-irradiated tissue. FIG. 2 shows the proposed approach for the identification of the tissue absorption and scattering coefficients, Ua and Us, respectively, during robot-assisted laser surgery. The laser light delivered to the tissue is absorbed under the form of heat, and the corresponding temperature increase is observed with a thermal sensor, i.e., an infrared thermal camera. The Ensemble Kalman Filter (EnKF) 150 is used to estimate the unknown coefficients based on the observed temperature dynamics.

[0027] At appropriate intervals, the computation circuit 120 iteratively revises the computed irradiation response coefficients based on the received thermal images. This includes comparing the temperature response based on the received thermal images with the computed irradiation response coefficients from a previous iteration. The filter performs tracking of time varying parameters based on a progression of the irradiation response coefficients in response to changes in the optical properties of the irradiated tissue 110 resulting from the irradiation signal 102.

[0028] FIG. 3 shows an irradiation signal as in FIG. 1 intersecting a tissue region or sample. Referring to FIGS. 1-3, tissue geometry as applied herein considers a rectangular block of tissue 110 exposed in air to an irradiation signal 102 such as a laser beam 102′. For simplicity, assume the tissue surface to be flat, and the laser beam to be perpendicular to the tissue. A Cartesian reference frame is established on the surface of the tissue, with the Z axis 170 coinciding with the optical axis 172 of the laser beam. Define the tissue temperature as a function T (x, y, z, t), where x, y, z are spatial coordinates, and t represents time. The tissue temperature can be calculated by solving the following differential equation:

[00001]$\begin{array}{cc}{c}_v\frac{\partial T}{\partial t}=k\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}\right)T+S,& \left(1\right)\end{array}$

where c.sub.v is the volumetric heat capacity of the tissue (J cm.sup.−3 K.sup.−1), k is the tissue thermal conductivity (W cm.sup.−1 K.sup.−1), and S is the volumetric power density (W cm.sup.−3). This latter term models the heat created by the laser in the tissue, and it is given by the product between the beam power (W) and the light absorption map A(x, y, z) (cm.sup.−3), which represents the fraction of light captured at any given location within the tissue volume.

[0029] To calculate the absorption map A, it is necessary to model the diffusion of light into the tissue. It is in these calculations that the coefficients of absorption μ.sub.a and scattering μ.sub.s (both having units of cm.sup.−1) appear. Obtaining a closed form solution for the absorption map can be challenging, and this quantity is frequently calculated with a Monte Carlo method instead; the idea is to simulate the optical path of a large number of photons as a discrete random walk, and to keep track of where the photons deposit energy. The length of each step of the walk is sampled from a logarithmic distribution, i.e.,

[00002]$\begin{array}{cc}s=\frac{-\ln \left(\zeta \right)}{{\mu }_a+{\mu }_s}& \left(2\right)\end{array}$

where ζ is a computer-generated number sampled uniformly at random between 0 and 1. When a photon moves from one step to the next, its direction of travel will change due to scattering. This change in direction is modeled by means of an azimuthal component, sampled uniformly at random between 0 and 2π, combined with a deflection angle θ, which is typically modeled using the Henyey-Greenstein function, i.e.,

[00003]$\begin{array}{cc}p\left(\cos \left(\theta \right)\right)=\frac{1-g^2}{2{\left(1+g^2-2g\cos \left(\theta \right)\right)}^{3/2}},& \left(3\right)\end{array}$

with g being the expected value of cos(θ). This parameter is also known as the anisotropy factor, and for most biological tissues, its value has experimentally been determined to range between 0.7 and 0.99.

[0030] At each step of the walk, a photon loses a fraction of its energy due to absorption. A photon is terminated either when it escapes the tissue volume or when its residual energy level falls below some arbitrary small positive value €.

[0031] To estimate the absorption and scattering coefficients, μ.sub.a and μ.sub.s in (2), the computation circuit 120 uses an augmented EnKF 150. Given the observed sensor data, the goal is to formulate an approximation of the joint probability density function π(T, μ.sub.a, μ.sub.s) using a discrete sample. For conciseness of notation, we introduce a vector θ=(μ.sub.a, μ.sub.s), so that the probability density function can simply be written as π(T, θ).

[0032] Assume that we have a set of measurements d.sub.j of the obtained sequentially by the thermal sensor at discrete times t.sub.j, with j=1, . . . , M. Further assume that the data are corrupted by measurement errors. Let T.sub.j denote the temperature predicted by the laser-tissue interaction model at time j, and let θ.sub.j=(μ.sub.a,j, μ.sub.s,j) be a vector containing the parameter estimates at time j. The filtering process begins by drawing a random sample of size N from the prior distribution π(T.sub.0, θ.sub.0), which encodes any prior knowledge on the unknown coefficients. This forms the initial ensemble at time j=0. The filter then proceeds in a two-step updating scheme from time j to j+1, as described below.

1) Prediction Step: Given the Current Ensemble

[0033]
S.sub.j={(T.sub.j.sup.n,θ.sub.j.sup.n)}.sub.n=1.sup.N

at time j, the prediction step of the EnKF updates the temperature values using a model approximation; i.e.,

T.sub.j+1|j.sup.n=F(T.sub.j.sup.n,θ.sub.j.sup.n)+v.sub.j+1.sup.n,n=1, . . . ,N  (4)

Where F(T.sub.j.sup.n, θ.sub.j.sup.n) represents the numerical solution to (1) at time j+1, stored as a column vector, and the innovation

v.sub.j+1.sup.n˜(0,C)

accounts for uncertainty in the forward prediction. The parameter values

θ.sub.j.sup.n.

are propagated forward using a random walk model of the form

θ.sub.j+1|j.sup.n=θ.sub.j.sup.n+ξ.sub.j+1.sup.n,n=1, . . . ,N  (5)

where

ξ.sub.j+1.sup.n˜(0,E)

with a prescribed covariance matrix E. Note that the parameter forward prediction in (5) is not necessary if the parameters are assumed to be constants; however, it is vital in tracking time-varying parameters. The predicted temperature and parameter values are augmented into vectors of the form

[00004]$\begin{array}{cc}{z}_{j+1|j}^n=\left[\begin{array}{c}{T}_{j+1|j}^n\\ {\theta }_{j+1|j}^n\end{array}\right],n=1,.Math.,N& \left(6\right)\end{array}$

which are used to compute the mean

[00005]$\begin{array}{cc}{\stackrel{\_}{z}}_{j+1|j}=\frac{1}{N}\underset{n=1}{\overset{N}{.Math.}}{z}_{j+1|j}^n& \left(7\right)\end{array}$

and the covariance

[00006]$\begin{array}{cc}{\Gamma }_{j+1|j}=\frac{1}{N-1}\underset{n=1}{\overset{N}{.Math.}}\left(z_{j+1|j}^n-{\stackrel{\_}{z}}_{j+1|j}\right){\left(z_{j+1|j}^n-{\stackrel{\_}{z}}_{j+1|j}\right)}^T& \left(8\right)\end{array}$

of the prediction ensemble.
2) Analysis Step: During the analysis step, the observed data d.sub.j+1 are assimilated in producing the posterior ensemble, which is computed by

z.sub.j+1.sup.n=.sub.j+1|j.sup.n+K.sub.j+1(d.sub.j+1.sup.n−G(z.sub.j+1|j.sup.n))  (9)

for each n=1, . . . , N. Here

d.sub.j+1.sup.n=d.sub.j+1+w.sub.j+1.sup.n,n=1, . . . ,N  (10)

generates an ensemble of fictitious measurements around the observed data d.sub.j+1, with

w.sub.j+1.sup.n˜(0,D)

representing observation error; G is the observation function, which maps the predicted states and parameters to corresponding model observations; and K.sub.j+1 is the Kalman gain matrix, which contains cross-correlation information between the predicted model states and parameters. In this work, the observation function G in (10) is a linear mapping

G(z.sub.j+1|j.sup.n)=P z.sub.j+1|j.sup.n,n=1, . . . ,N  (11)

where the projection matrix P picks out the tissue location at which the temperature is being measured, and the Kalman gain is computed by

K.sub.j+1=Γ.sub.j+1|jP.sup.T(PΓ.sub.j+1|jP.sup.T+D).sup.−1.  (12)

Posterior ensemble statistics are computed similarly as in (7)-(8) using

z.sub.j+1.sup.n.

The process repeats for j<M. Note that here we assume that temperature data are available at each time j; if measurements are not available at a subset of the filter time steps, the observation update can be neglected such that the prediction ensemble serves as the posterior at these steps.

[0034] To study the viability of the proposed method in estimating the unknown tissue absorption and scattering coefficients, we conducted a set of numerical experiments. Simulations were performed using the MATLAB® programming language (The MathWorks, Inc., Natick, Mass., well known in the art). Other computation packages may also be employed. The computation circuit 102 is a computer having a processor and memory suitable for executing and performing the simulations and related calculations, and for receiving the temperature response and rendering the results.

[0035] In each experiment, we simulate the setup previously illustrated in FIG. 1, where a laser pulse is applied to the surface of a sample of biological tissue. The laser beam has a uniform intensity profile, with a radius of 1 mm and a power of 0.5 W. The tissue is assumed to have the physical properties summarized in Table I. The dimensions of the tissue sample are 0.5 cm 0.5 cm 0.25 cm. For computational purposes, the tissue 110 block is discretized into a grid of 100×100×50 cubic elements, or voxels. The interaction between the laser beam and the tissue is monitored with a virtual thermal sensor, which provides superficial temperature measurements at a rate of 10 Hz. To simulate the presence of sensor noise, temperature measurements are altered with the addition of an error term e˜N (0, 0.01).

[0036] As noted, the goal is to estimate the absorption and scattering coefficients of the tissue given the thermal sensor data. At the beginning of each experiment, we initialize the EnKF with an ensemble of size N=50. The initial tissue temperature is set to 0° C., and the initial values for the unknown tissue coefficients are drawn from uniformly distributed priors, i.e., U (0.5θ{circumflex over ( )}, 2θ{circumflex over ( )}), where θ{circumflex over ( )} are the true coefficient values listed in Table I. To perform the temperature prediction step (§ II-B), the filter internally runs its own implementation of the laser-tissue interaction model described in section II-A with a coarser tissue grid (i.e. 20×20×10 voxels) than the one used to simulate the sensor data. The use of a coarser forward model helps to limit the computational complexity of the filter. While the absorption and scattering coefficients are unknown, we assume that the EnKF has knowledge of the other two tissue parameters in the thermal model, namely, the volumetric heat capacity, c.sub.v, and thermal conductivity, k. We note that these additional parameters may not always be known a priori in a realistic setting, but reasonable approximations can generally be obtained using empirical models available in the laser tissue interactions literature.

TABLE-US-00001 TABLE I TISSUE PHYSICAL PROPERTIES USED IN SIMULATION Symbol Physical Variable Units Value Used μa Absorption Coefficient cm−.sup.1 1 μs Scattering Coefficient cm−.sup.1 100 cv Volumetric Heat Capacity J/(cm−.sup.3 ° C.) 3.76 k Thermal Conductivity W/(cm ° C.) 0.0037

[0037] FIGS. 4A-E and 5A-E respectively describe two sets of numerical experiments: one in which the absorption and scattering coefficients remain constant over time; another in which the absorption and scattering coefficients shift during laser exposure (as may occur during a laser procedure). Each graph depicts a solid line 190 representing the EnKF mean, upper and lower dashed lines 192-1 . . . 192-2 showing an uncertainty bound flanking a true value shown by a dashed line 194.

[0038] FIGS. 4A-4E show graphs of constant parameter estimation using the filter of FIG. 2, depicted in a scenario in which the tissue optical properties remain constant during laser irradiation. FIGS. 4A-4W show application of a laser pulse for 5 seconds, then continue to observe the tissue temperature for 10 more seconds. We utilize the EnKF with the aim of estimating the true tissue absorption and scattering coefficients listed in Table I.

[0039] FIGS. 4A-4B show the resulting estimates of μ.sub.a and μ.sub.s, respectively, along with the corresponding tissue temperature estimates at three locations on and below the surface of the tissue; namely, we estimate the temperature profiles at the sensor location (x, y, z)=(0 cm, 0 cm, 0 cm) and at depths of 0.1 cm and 0.2 cm directly below (FIGS. 4A-4C, respectively). Note that the EnKF estimate in each plot is the ensemble mean, with uncertainty bounds given by ±2 standard deviations around the mean.

[0040] As seen in FIGS. 4A-4E, the EnKF provides an accurate estimate of the absorption coefficient, with uncertainty bounds shrinking over time. The estimate of the scattering coefficient drifts to a slightly higher value after the laser is turned off, with wider uncertainty bounds that do however contain the true parameter value. Further, while only observing data at the surface sensor location, the filter is able to well estimate the tissue temperature at the two locations tracked below the tissue surface, with wider uncertainty bounds as the location becomes farther from the surface.

[0041] In this experiment depicted in FIGS. 5A-5E, the procedure described above is repeated, but simulating a scenario where the tissue optical properties change during laser exposure. Laser induced alterations in the absorption and scattering coefficients have been documented in prior literature. To simulate these shifts, we model both the absorption and scattering coefficients as continuous piecewise functions that increase linearly during the laser pulse and remain constant when the laser is turned off. More specifically, we let

[00007]$\begin{array}{cc}{\mu }_a\left(t\right)=\{\begin{array}{cc}0.6t+1& \mathrm{if}0\le t\le 5\\ 4& \mathrm{if}5<t\le 15\end{array}\mathrm{and}& \left(13\right)\end{array}$$\begin{array}{cc}{\mu }_s\left(t\right)=\{\begin{array}{cc}80t+100& \mathrm{if}0\le t\le 5\\ 500& \mathrm{if}5<t\le 15\end{array}& \left(14\right)\end{array}$

respectively. It should be noted that the EnKF does not assume any knowledge of the relations above, and that the goal of this experiment is precisely to verify if the filter is able to track the absorption and scattering coefficients as their values change over time.

[0042] FIG. 5A-5E displays the resulting estimates of μ.sub.a(t) and μ.sub.s(t) (5D and 5E), along with the corresponding tissue temperature estimates at the three aforementioned locations (5A-5C, respectively). These results show that the filter is able to well track the change in absorption throughout the duration of the experiment. The increase in the scattering coefficient is more difficult to track during the laser pulse, but the filter is able to capture its behavior shortly after the laser is turned off. In both cases, the uncertainty bounds for the time-varying coefficients become increasingly wider once the laser is off. The estimates of the tissue temperature at and below the surface remain accurate throughout the experiment.

[0043] Those skilled in the art should readily appreciate that the programs and methods defined herein are deliverable to a user processing and rendering device in many forms, including but not limited to a) information permanently stored on non-writeable storage media such as ROM devices, b) information alterably stored on writeable non-transitory storage media such as floppy disks, magnetic tapes, CDs, RAM devices, and other magnetic and optical media, or c) information conveyed to a computer through communication media, as in an electronic network such as the Internet or telephone modem lines. The operations and methods may be implemented in a software executable object or as a set of encoded instructions for execution by a processor responsive to the instructions. Alternatively, the operations and methods disclosed herein may be embodied in whole or in part using hardware components, such as Application Specific Integrated Circuits (ASICs), Field Programmable Gate Arrays (FPGAs), state machines, controllers or other hardware components or devices, or a combination of hardware, software, and firmware components.

[0044] While the system and methods defined herein have been particularly shown and described with references to embodiments thereof, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the scope of the invention encompassed by the appended claims.