Identification of Controlled-Complexity Thermal Therapy Models Derived from Magnetic Resonance Thermometry Images

Medical imaging provides information valuable in diagnosis, planning, and control of therapies. In this paper, we develop a method that uses a specific type of imaging—the magnetic resonance thermometry—to identify accurate and computationally efficient site and patient-specific computer models for thermal therapies, such as focused ultrasound surgery, hyperthermia, and thermally triggered targeted drug delivery. The developed method uses a sequence of acquired MR thermometry images to identify a treatment model describing the deposition and dissipation of thermal energy in tissues. The proper orthogonal decomposition of thermal images is first used to identify a set of empirical eigenfunctions, which captures spatial correlations in the thermal response of tissues. Using the reduced subset of eigenfunction as a functional basis, low-dimensional thermal response and the ultrasound specific absorption rate models are then identified. Once identified, the treatment models can be used to plan, optimize, and control the treatment. The developed approach is validated experimentally using the results of MR thermal imaging of a tissue phantom during focused ultrasound sonication. The validation demonstrates that our approach produces accurate low-dimensional treatment models and provides a convenient tool for balancing the accuracy of model predictions and the computational complexity of the treatment models.


Introduction
In ultrasound (US) thermal therapies, the goal is to selectively heat the treatment target without excessively elevating the temperature in healthy tissues intervening in the path of transmitted US energy and surrounding the target. In recent years, temperature measurements obtained by MR thermometry has played an increasingly important role in planning and control of thermal therapies. A number of techniques have been developed to obtain MR temperature maps, with water proton resonance frequency (PRF) method being the most widely used in practice [1]. Several studies [2][3][4] have demonstrated that MR-thermometry can be used as a feedback in automatic treatment control systems. If we view images as a collection of pointwise measurements associated with each pixel and use these data to obtain models to plan, optimize and control therapies, the dimension of the resulting models will be high. For example, the pointwise utilization of 512|512 voxels in MR thermal images would lead to a treatment model with over 250,000 states. This presents a problem of a very high computational demand when a model must be used in real time during the therapy, as in the case of a model-based and optimizing treatment control systems, such as the one described in [4].
In this paper, we develop and validate methods that avoid pointwise utilization of imaging data, leading to a highly efficient compression of MR thermometry images and the identification of a low-dimensional dynamic treatment models. Fundamentally, a low-dimensional representation of MR thermal images is possible because the image voxels are spatially correlated, reflecting the dependence of temperature distribution in the treatment target on the specific absorption rate (SAR) of the ultrasound and the heat dissipation by conduction and convection. Furthermore, a time series of MR thermal images are correlated temporally because of the causal dependance of temperatures on the heating history and temperatures at the preceding times. In the developed approach, we apply proper orthogonal decomposition (POD) to identify orthonormal basis functions fw i g N i~1 which capture spatial correlations in the ensemble of N MR thermal images. Each element w i is identified to capture the maximum amount of spatial correlations that have not yet been explained by a subset of previously identified basis elements. We then select the reduced basis fw i g M i~1 , M%N, consisting of the first M basis elements, such that the desired balance between accuracy and complexity is obtained. The selected reduced basis is then used to: N Approximate each image acquired during therapy in the reduced set of basis functions. This step may be viewed as a real time image compression, which exploits spatial correlations of image voxels, as well as the image filtering step since high-frequency spatial valuations, which usually correspond to MRI measurement noise, are removed. N Identify the ultrasound specific absorption rate and the dynamic model that captures temporal correlation in the series of images and provides a prediction the evolution of the ultrasound treatment.
The developed approach is validated experimentally, using an in vitro MRI tissue phantom heated by a focused ultrasound transducer. The results show that the identified low-dimensional models predict the SAR and the temperature response of the target with the expected accuracy, which can be controlled by selecting the order M. Even with a very low-dimensional model, suitable for use in model-predictive treatment control systems, thermal images in the verification set were accurately approximated and the predictions of the SAR and the phantom's thermal response closely agreed with the MRI measurements.

Identification of basis functions using POD method
The proper orthogonal decomposition is a technique often used to extract a set of basis functions for an approximate, modal-like representation of an infinite-dimensional, distributed parameter system (DPS). It was previously shown [5] that POD is the most efficient way to obtain dominant modes of a known infinitedimensional dynamic system, which makes it a popular approach in a variety of applications, including image processing [6]. Often, a numerical solution of a given partial differential model at different times (known as snapshots) is used as the input to the POD algorithm to produce the desired basis functions. The identified basis is then used in combination with the projection methods to obtain a finite-dimensional approximation of the original known infinite-dimensional DPS model [7].
In this paper, we use a time series of images to identify an orthonormal basis of an unknown infinite-dimensional model, characterizing the heat transport in the target due to noninvasive ultrasound heating. The following brief outline gives the exposition to the POD method in the context of our objectives. Further theoretical and algorithmic details are available in [5,[7][8][9].
Consider a set of MR thermal images (snapshots) S~fT m (r,t i ) : 1ƒiƒN,r [ Vg. Here, T m (r,t i ) is an image of the region of interest (ROI), characterized by a vector of coordinates r inside the spatial domain V and acquired at time t i . The problem is to obtain a function w(r), which is the best at capturing the spatial distribution of temperatures in the ensemble S of snapshots T m (r,t i ). Mathematically, the problem is to find w(r), such that the projections of all snapshots T m (r,t i ) onto w(r) are maximized: where Sf ,gT~Ð V f (r)g(r)dr denotes the inner product of square integrable functions f (r) and g(r) defined in V. The normalization condition Sw,wT~1 is imposed to ensure uniqueness of the solution.
The optimization problem (1) is difficult to solve in the general case. The problem is simplified if, following [10], we seek the solution under an additional assumption that w(r) can be expressed as a linear combination of snapshots: and all basis functions w k (r) are orthonormal: In this case, it can be shown [11] that the solution of the optimization problem is reduced to the solution of the following matrix eigenproblem: where C~fc ij g is the covariance matrix of all snapshots with elements For a strictly positive definite matrix C, equation (4) is satisfied by N orthogonal eigenvectors w k~½ w k 1 w k 2 Á Á Á w k N T of dimension N and the corresponding distinct eigenvalues l k , where k~1, . . . ,N. The elements w k i of an eigenvector w k are the coefficients in the sought linear decomposition (2), and thus determine the k-th eigenfunction w k . The requirement that Sw k ,w k T~1 of problem (1) is enforced by normalizing all eigenvectors, such that Assuming that all eigenvalues are ordered (l 1 §l 2 §::: §l N ), the eigenfunction w 1 , corresponding to l 1 , is the desired solution w of the maximization problem (1). The normalized eigenfunctions fwg N i~1~f w 1 ,w 2 , Á Á Á ,w N g form an orthonormal basis of the image ensemble, S. The amount of information, captured by the projection of S on the i-th eigenfunction, is characterized by the corresponding eigenvalue: Consequently, w 1 is the best at explaining spatial correlations in S, followed by w 2 as the next most informative direction, and so on. When the complete set of the identified eigenfunctions is used to represent images in S, there is no loss of information. There is no benefit either since the basis order N is equal to the number of images in S. The question is how to select a reduced basis of order M%N to obtain the desired accuracy of the image approximation in the reduced basis. To answer this question, we start by defining the ''energy'' of an image ensemble as It is easy to show that Therefore, we can use the eigenvalues to guide the selection of the reduced-order POD basis. One approach is to select M such that a predetermined fraction e (ƒ1) of the total energy of the ensemble S is captured. Specifically, we may wish to select the smallest M such that where e is selected by a user. The reduced-order basis fw j (r)g M j~1 can then be used to parsimoniously approximate all images in the ensemble S, as well as all new images acquired during ultrasound treatment. In the sequel, this reduced basis is also used to identify low-dimensional patient-and site-specific ultrasound SAR and thermal treatment models based on the information captured by MR thermometry images.

Image representation in the reduced basis
Consider a thermal image T m (r,t k ), collected at time t k , which reflects the spatial temperature distribution in the region of interest r [ V in response to the sonication history over the time tƒt k . Thermal effect of sonication is described by the power deposition, Q(r,t), which relates to the ultrasound specific adsorption rate (in W/kg) as The problem for approximating a new image can be viewed as the minimization problem of finding projectionsT T m (t) of an image T m (r,t k ) into the manifold, spanned by a reduced basis of empirical eigenfunctions fw j g M j~1 : The solution of (11) is reduced to finding the best solution (in an appropriately selected norm) of the following equation: where W(r)~½w 1 (r) w 2 (r) . . . w M (r) is a row vector of eigenfunction. We will now take into account that r takes only discrete values in V due to image pixelation (finite spatial resolution). Therefore, each basis function, formed according to equation (2), is given by its values in the spatial locations of image voxels, or w j (r)&½w j (r 1 ) w j (r 2 ) . . . w j (r Nvox ) T , where r l is the coordinate of lth voxel and N vox is the total number of voxels. With finite spatial resolution, a vector of function W becomes the N vox |M matrix, which transforms (12) into a system of N vox linear equations in M unknown, and by solving (12) in the least squares sense, obtained. Note that the repeated solution of (12) for each acquired image can be accelerated by parameterizing matrix W (for example, by calculating its LU decomposition).
With the described procedure, the temperature measurements T m (r,t k ) are completely encoded, to the desired accuracy, by a lowdimensional vector of projectionsT T m of dimension M, where M%N vox . Note that only projectionsT T m (t k ) must now be sent from the MRI scanner, which reduces communication traffic. The storage requirements are also reduced because only a Mdimensional vector of projections must be saved for each newly acquired thermal image of N vox (&M) voxels.
Our previous results [11] indicate that by ignoring the contribution of higher order eigenfunctions fwg N i~Mz1 in image representation spatial noises in MR thermometry images are filtered.

Representation of ultrasound power deposition
The temperature distribution measured by MR thermometry depend on the power deposition in tissue, Q (in W : m {3 ), caused by sonication. Therefore, it is reasonable to expect that fwg M i~1 , identified using the acquired images, is also an adequate basis to represent Q(r,t k ). Similarly to equation (12), for known Q(r,t k ), the projection of the power deposition into the reduced basis, u u(t k )~½û u 1 (t k )û u 2 (t k ) Á Á Áû u M (t k ) T , can be found as the least squares solution of the following linear equation: Identification of the projection models of thermal response Our objective is to to identify an ultrasound thermal treatment model in a low-dimensional projection form with the state vector, T T, corresponding to the projection of measurements,T T m . A lowdimensionality is required to enable real time utilization of treatment models for such tasks as intra-operative prediction of specific sonication parameters on temperature distribution and to design efficient model-based treatment controllers. Figure 1 shows a block diagram of a control system that uses projectionsT T m (t k ) (instead of full-dimensional images, T m (r,t k )) in its feedback and determines control inputs,û u, using low-dimensional treatment model in projection form. This figure illustrates that during ultrasound therapy the images from the MRI scanner are sent to the control system in low-dimensional form of image projectionŝ T T m , and the controller uses this information to generateû u d (t k ), which describes Q d (r,t k ) -the desired power deposition which we wish to apply to the patient at the current time;û u d (t k ) and the corresponding Q d (r,t k ) are the control inputs sent to the ultrasound subsystem. The best approximation of Q d (r,t k ) implementable with the available ultrasound transducers or transducer array system, Q(r,t k ), is then applied to the patient. As desired, the described control system entirely avoids highdimensional (pointwise) computations.
We seek to identify the low-dimensional treatment model in the following form:T whereû u is the projection of the ultrasound power deposition, Q(r,t k ); matrices Y and C are unknown, and must be identified from the images T m (r,t) in the training ensemble, S. The state vectorT T is the predicted temperature projection; givenT T, the predictions of the temperature distribution in the region of interest (ROI) can be obtained as It is straight forward to show that all linear PDE models, traditionally used to describe temperature evolution in biological tissues during thermal treatments (such as the convection-diffusion model and its different approximations), can be adequately approximated by model (14). Consider, for example, the following popular Pennes' bioheat transfer equation (BHTE) [12], which is written in terms of the temperature increase T(r,t) above the equilibrium: where r (kg/m 3 ) and k (W/m 0 C) are the tissue density and thermal conductivity, C and C b are the specific heat of tissue and blood (in J/ kg 0 C), respectively. Equation (16) does not require the detailed information on tissue vascularity or blood flow but instead uses an empirical blood perfusion-related parameter W b (kg/m 3 s). After using T(r,t)~W(r)T T(t) and Q(r,t)~W(r)û u(t) in equation (16), and approximating time derivative using backward-difference approximation, the BHTE takes the following form: where Dt is the time discretization step. The weak Galerkin formulation of the Pennes' model is obtained by taking the inner product of equation (17) with the elements w i (r) of the reduced basis, yielding the following system of discrete-time dynamic equations: Dt½Sw i (r),a(r)+ 2 w j (r)T{ where a(r)~k rC , b(r)~W b C b rC , c(r)~1 rC and we took into the account the orthonormality of the basis functions. The compact form of the M-dimensional projection model (18) is exactly the same as equation (14) with matrices Y~fa ij g, C~fb ij g defined by their elements as: where d ij is the Kronecker-delta, and i,j~1, . . . ,M. Note that for a constant tissue density and heat capacity, C is a diagonal matrix cDtI M|M . Matrix Y characterizes the heat dissipation in the target due to conduction and blood flow, while the affine term Cû u(t k ) describes the effect of the ultrasound power deposition Q on the evolution of projected temperature responsesT T(t k ).
The problem is to identify Y, C andû u such that projectionŝ T T(t k ), predicted by the model (14) and used to reconstruct temperature distribution as T(r,t)~W(r)T T(t), give the best agreement with the acquired series of images T m (r,t) [ S.
Identification of matrix Y. In order to decouple the problems of identifying Y, C andû u, the system matrix Y is identified first by utilizing only thermal images T m (r,t k ) acquired during tissue cooling (Q~0) to thermal equilibrium. With zero power inputû u~0, the projection model (14) describes the decay of temperature projections to thermal equilibriumT T~0 from nonzero initial conditions:T This equation is in the form of a first-order multivariate autoregressive (AR) model, which allows us to identify Y using the standard system identification techniques [13]. Identification of the affine term. The specific absorption rate of an ultrasound transducer depends on tissue properties and other variable factors, which can change with treatment site and from patient to patient. At the same time, the knowledge of the patient-and site-specific SAR is critically important for treatment safety and its precise control.
Within the developed approach, the identification of the SAR is accomplished easily. Once the system matrix Y is determined, the time-varying affine term Cû u can be identified at each sampling instant t k as: whereT T m (t k ) andT T m (t kz1 ) are the projections of the consecutive thermal images, collected during sonication of the target. Since C (cf. equation 20) depends only on the MRI scan time Dt, the identified eigenfunctions, and the relatively little-changing tissue density and heat capacity (both often assumed equal to the water values), the matrix C in the projection model can be considered as known. With this assumption,û u(t k ) can be found as a least squares solution of the linear equation (22). The corresponding estimates of Q(r,t k ) and the SAR then be obtained by using equations (13) and (10).

Results
The developed approach to identifying low-dimensional models of ultrasound therapies was tested experimentally using MR measurements of the thermal response of a tissue phantom to focused ultrasound (FUS) sonication. The cubic 11|11|11 cm phantom was prepared following the recipe of Madsen et al. [14]. The T2 relaxation time of the phantom was modified by adding one millimole-per-liter of copper sulfate to the recipe. After preparation, the phantom was allowed to solidify inside of an acrylic box with a Mylar ultrasound treatment window on the bottom surface. The ultrasound power deposition field was created by a single, spherically focused, air backed ultrasound transducer with aperture diameter of 10 cm and radius of curvature of approximately 18 cm, and resonating frequency of 1.5 MHz. The transducer was placed in a bath of degassed and deionized water inside the MR compatible ultrasound positioning system (Figure 2 and [4]), which was designed to move the focal zone in three dimensions. After initial alignment, the transducer position remained fixed in the current experiments. The transducer was driven by a function generator (Stanford Research System, model DS345), amplified using RF amplifier (ENI, model A150). The electrical impedance of the transducer was matched to the output impedance of the amplifier using an external LC matching circuit. The electrical power applied to the transducer (direct and reflected) was measured using power meters (Hewlett-Packard, model 435A/B). Further details of driving circuitry and positioning system are given in [15].
Temperature increase, T m (r,t k ), inside the phantom was measured using Siemens Trio 3T Magnetom scanner. A custom receive-only surface coil was used to improve the temporal and spatial resolution of the acquired thermal images. The coil created a localized sensitivity pattern, which minimized the interferences and improved the signal-to-noise ratio (SNR). Gradient-echo sequence (GRE) with the following acquisition parameters was used to obtain temperature measurements: repetition time TR~30 ms, echo time TE~10 ms, field of view FOV~25:6|25:6 cm and flip angle~25 0 . The voxel size of thermal images was 2|2|3 mm. The scan time was 2.45 s with the phase resolution of fifty percent to increase the sampling rate. The overall image size was N vox~1 28|64 and the k-space data were zero-filled to form a 128|128 data matrix.
Following the PRF shift method [16], the phase difference Dw between the two consecutive phase images was used to calculate the relative temperature change T m as where a~{0:01 ppm C {1 [17] is the coefficient of PRF shift for aqueous tissue, c g is the gyromagnetic ratio, B 0 is the strength of the main magnetic field, and TE is the echo time. Sequential complex MR images S m (r,t k ) and S m (r,t kz1 ) were used to calculate Dw(r,t kz1 ) as the phase of the following product: where Ã denotes the complex conjugate operator. Figure 3 gives a representative temperature image in the US transducer's focal plane (gradient scale is in 0 C). The phantom appears as a rectangular object above the ultrasound positioning system containing a clearly visible 45 0 ultrasound mirror. The selected ROI (r~(x,y) [ V) is the region of an appreciable temperature elevation, which has pixel coordinates 53ƒxƒ58 and 34ƒyƒ59. The number of voxels in the ROI is 6|26; its actual dimension is 1:2 cm|5:2 cm. As expected, the maximum temperature rise is observed on the line of focal symmetry, x~56.

Validation Results
A step increase from zero to either 3.5, 4.8 or 6.5W of total electrical power (direct minus reflected) was applied to the FUS transducer while keeping all other experimental conditions the same. The phantom was allowed to reach thermal equilibrium before changing the power input.
MRI thermal images collected for the case of 6.5W of applied electrical power were used as a training (estimation) dataset to identify the reduced-order basis, and the corresponding models of the thermal response and the ultrasound SAR. The validation datasets, consisting of the images collected during the experiments with the other two power levels, were used to assess the adequacy of the identified basis functions and the accuracy of the temperature predictions obtained with the identified ultrasound treatment model.
A total of 499 MR thermal images were acquired to characterize temperature evolution during each power step test. The estimation dataset included images collected when the power was kept constant at 6.5W (0vtƒ773 seconds) and when the tissue was cooling back to equilibrium (773vtƒ1225) after the power was switched off. Figure 4(a) shows the measured temperature elevation within the ROI at t~773 s when the temperature reaches its peak value. The temporal evolution of temperatures in the selected locations on the line of ultrasound beam symmetry is shown in Figure 4(b). The highest temperature was observed at (x,y)~(56,54), where the maximum temperature increase from ambient temperature was *17:5 0 C.
Identification of the reduced-order basis. All 499 images in the estimation dataset were used to identify the orthonormal basis following the described method. Figure 5 shows the first four identified eigenfunctions and the corresponding eigenvalues, which rapidly decay for the higher order w i 's. Using the criterion (9), it was determined that the first eigenfunction captures approximately 98:93% of the spatial correlations in the collection of 499 images, while w 2 (r) captures only 1:01%. Selecting e~99:9%, we conclude that high accuracy of image approximations in the reduced POD basis fwg M i~1 is achieved with only two basis functions (M~2). Note that the shape of w 1 , identified to maximize the explained spatial correlations in thermal images, is similar to the shape of temperature distribution, as expected. Further note that eigenfunctions w i , i~3,4, . . . capture information at increasingly higher spatial frequencies, and that ignoring their contribution in image representation (15) has the effect of a spatial filtering of imaging data.

Identification
of thermal response and SAR models. Thermal images, collected during tissue cooling after the initial sonication at 6.5W of constant applied power, were used to identify the system matrix Y of the projection model (14). The result of Table 1 was obtained by, first, using Matlab's System Identification Toolbox [13] to obtain the autoregressive model of the cooling process and then converting it into the state space form of equation (21). Handbook values of heat capacity and density (C~4186 J/ (kg 0 C and r = 1000 kg/m 3 ) [18]) were used to calculate c in equation (20). Sequential thermal images were used in equation (22) to estimate the power deposition projectionû u. The identified vectorû u (Table 1) was then used in equations (13) to estimate Q(r); the corresponding SAR(r) was obtained using equation (10). The result (SAR 3 ), scaled with the total applied power of 6.5W, is shown in Figure 6(a). A high degree of correlation between the shapes of the measured temperature distribution (Figure 4a) and the SAR is an expected result for the unperfused phantom.
The prediction of the thermal response to different power inputs is based on the assumption that the shape of the SAR remains constant for the fixed relative transducer-patient position, but the SAR values are scaled with the applied electrical power. This assumption of the linear SAR-power dependence was tested by comparing the predicted SAR pattern for 3.5 and 4.8W of applied power, obtained by scaling SAR 3 by 3.5/6.5 and 4.8/6.5, with the directly identified SAR distribution (SAR 1 and SAR 2 ). The direct identification of SAR 1 and SAR 2 followed the same method used to identify SAR 3 (i.e. a new reduced POD basis and the corresponding thermal model were obtained in each case; equation (22) was used to calculate a newû u; then SAR 1 and SAR 2 were reconstructed using equations (13) and (10)). Figures 6(b) and (c) shows the difference between the predicted and the directly identified SAR 1,2 , both scaled to 1W of the  Validation of the thermal response model. The accuracy of the identified reduced-order treatment model was assessed in projection manifold and in terms of temperature predication errors. In projection space, the prediction of the temperature projection vectorT T(t), obtained with the same two-dimensional state space model of Table 1 and the appropriately scaledû u, was compared with the projection of the actual thermal images T m (r,t), acquired during the experiments with three different power levels. Image projections,T T m (t), were found as the least squares solution of equation (12). Figure 7 shows an excellent agreement between the predictions, T T(t), and measurements,T T m (t). The agreement is the best for 6.5W of applied power -the case used as the estimation dataset. The plot ofT T 1 (t) is similar in shape to the pointwise temperature evolutions (cf. Figure 4b), which indicates that the first component of the vectorT T(t) captures most of the slow temporal variations in the series of thermal images.
When used in equation (15), the model-generatedT T(t) gives the prediction of the temperature elevation, T(r,t), in the ROI, which can be compared with imaging data, T m (r,t). Figure 8 shows the spatial mean and standard deviation (STD) of the temperature prediction errors for all pixels in the ROI. The prediction errors are small, including the two validation cases shown in subplots (a)-(d). The maximum pointwise temperature prediction errors do not exceed 1 0 C, which is of the same order as the measurement noise of MRI thermometry.

Discussion
The developed approach was shown to be effective in identifying low-dimension but accurate models of ultrasound thermal therapies. At a pre-treatment stage, a set of MR thermometry images, characterizing the response of the target and the surrounding normal tissue to thermal excitation, is collected and then used to identify the reduced POD basis, which capture spatial correlations in images. A simple criterion for selecting an appropriate number of basis functions is provided which allows a user to balance the computational complexity of a predictive treatment model with its computational complexity. The selected reduced basis is then used to parsimoniously approximate newly acquired images, thus minimizing storage and data traffic. As an additional benefit, image approximation in the reduced basis filters high-frequency spatial noises in MR images.  The SAR and thermal response models identified following the developed approach are patient-and site-specific and can be used as predictive models in real time (e.g., treatment control) and off line (such as treatment planning) applications. The identification procedures are well suited to perform continuous re-identification of treatment models during the therapy. Such intra-treatment adaptability helps to mitigate the effect of changing tissue properties (such an US absorption) and blood perfusion, caused by elevated temperatures, on the accuracy of model predictions, which is particularly important when the identified low-dimensional models are used to efficiently implement modelbased, optimizing treatment controllers that utilize images in the feedback. A family of related results used to identify continuous-time treatment models is described in reference [19].
The developed methods were validated during in vitro MR experiments with a tissue phantom heated by a focused ultrasound transducer. The experimental results indicate that the SAR and thermal response during the treatment can be accurately predicted by the identified projection models with only two states. The lowdimensionality of the identified models substantially minimizes computational requirements of implementing a model-based treatment control system and communication traffic between the MRI scanner and the treatment controller.   Table 1 was used to make predictions for the case of 6.5W of applied power (plots (a) and (b)). Predictions for the power inputs of 3.5W (plots (c), (d)) and 4.8W are made by scalingû u by 3.5/6.5 and 4.8/6.5, respectively. doi:10.1371/journal.pone.0026830.g007 Though our emphasis is on thermal therapies, the developed approach has a broader applicability in image-based identification and image-guided control of therapies. After straightforward modifications, this approach can be used with images acquired in multiple planes and with three-dimensional MR measurements.