Detection of time-varying heat sources using an analytic forward model

TL;DR

This paper introduces an analytic model for detecting static and dynamic heat sources by solving the 3D heat equation, with algorithms tested via simulations, applicable in material testing and medical imaging.

Contribution

It provides a novel simple analytic point source model for time-varying heat sources and algorithms for source detection and spectral analysis.

Findings

Algorithms successfully detect source location and spectral content.

Numerical tests demonstrate the model's effectiveness.

Framework shows promise for practical inverse heat source problems.

Abstract

We present a simple, analytic point source model for both static and time-varying point-like heat sources and the resulting temperature profile that solves the heat equation in dimension three. Simple algorithms to detect the location and spectral content of these sources are developed and numerically tested using Finite Element Mesh simulations. The resulting framework for heat source reconstruction problems, which are ill-posed inverse problems, seems promising and warrants for future research. Possible fields of application for our work are material testing, to detect manufacturing defects, and medical imaging to detect abnormal health conditions.

Tables

Table 1. Table 1. Heat source parameters for the simulations.

	static, R=10	dynamic, f=0.2Hz
center	$Q_1{Q}_2{Q}_3{Q}_4$	$S_1{S}_2$
A [40 40 50]	$1.00.80.60.4$	$2.21.8$
B [35 35 45]	$1.00.80.60.4$	$2.21.8$
C [30 30 40]	$1.00.80.60.4$

Table 2. Table 2. Static reconstructions for the three source locations, four power levels and three simulated noise levels.

	location error	$Q_1$	$Q_2$	$Q_3$	$Q_4$
center A
no noise	5.0	1	0.83	0.56	0.39
noise 1%	5.2	1.1	1.0	0.79	0.48
noise 5%	10.2	1.6	1.5	1.1	0.67
noise 10%	10.5	1.4	1.3	1.0	0.76
center B
no noise	6.6	0.90	0.48	0.38	0.34
noise 1%	7.0	1.0	0.68	0.52	0.57
noise 5%	19	1.0	0.68	0.87	0.52
noise 10%	14	2.7	0.80	0.66	0.75
center C
no noise	11	0.61	0.41	0.33	0.34
noise 1%	8.7	0.90	0.56	0.52	0.80
noise 5%	19	1.6	0.50	0.59	0.44
noise 10%	19	3.0	0.84	0.68	0.69

Table 3. Table 3. Dynamic reconstructions for the source location A, two modulation levels 2.2 and 1.8 and two simulated noise levels.

	location error	spectrum magnitude S
center A	(amp)     (phase)	(amp)     (phase)
no noise(A1)	6.1     5.0	2.2     2.2
no noise(A2)	6.1     5.0	1.8     1.8
noise 1% (A1)	6.1     8.2	2.2     4.0
noise 3% (A2)	6.1     8.2	1.8     3.2

Table 4. Table 4. Dynamic reconstructions for the source location B, two modulation levels 2.2 and 1.8 and two simulated noise levels.

	location error	spectrum magnitude S
center B	(amp)     (phase)	(amp)     (phase)
no noise(B1)	3.8     1.5	5.2     438
no noise(B2)	3.8     2.3	4.2     2409
noise 1% (B1)	5.3     6.8	13.4     411
noise 3% (B2)	5.3     8.1	11.0     14945

Equations

$$q(x,t)=s(t)\delta (x_0)\text{ in }\Omega ,$$

$$s(t)=\sum _{n=0}^N{S}_n\cos (2\pi f_{n}t+{\varphi }_n).$$

$$c(x)\frac{\partial T}{\partial t}=\nabla \cdot (\kappa (x)\nabla T)+q(x,t),$$

$$q(x,t)=\frac{1}{c}s(t)\delta (x_0),$$

$$\frac{\partial T}{\partial t}=\alpha \Delta T+q(x,t).$$

$$G(x,t)=\frac{1}{(4\pi \alpha t{)}^{n/2}}{e}^{-\frac{x\cdot x}{4\alpha t}},$$

$$T(x,t)=\frac{1}{(4\pi \alpha {)}^{n/2}}{\int }_0^t\frac{e^{-x\cdot x/(4\alpha (t-y))}}{(t-y{)}^{n/2}}s(y)\mathbf{d}y,$$

$$T(x,t)$$

$${\partial }_t{T}_1$$

$${\partial }_t{T}_2$$

$${\partial }_{t}T=\alpha \Delta T+q_1+q_2.$$

$$q(x,t)=q_t(x,t)+q_c(x),$$

$$\rho c\frac{\partial T}{\partial t}=\nabla (k\cdot \nabla T)+w_b{\rho }_b{c}_b(T_a-T)+Q_m,$$

$$\kappa \frac{\partial T}{\partial \nu }=h(T_{atm}-T)\text{on }\partial \Omega ,$$

$$\Delta T=-\frac{Q}{\alpha }\cdot \delta (r),$$

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2\frac{\partial T}{\partial r}\right)=-\frac{Q}{\alpha }\cdot \delta (r).$$

$$T(r)=\frac{Q}{4\pi \alpha r}+C,$$

$$q(x,t)=S_n\cos (2\pi f_{n}t+{\varphi }_n)\delta (x_0).$$

$$f(w)=\frac{e^{-x\cdot x/w}}{w^{n/2}},\lambda =2\pi f_{n}t+{\varphi }_n,\beta =\frac{2\pi f_n}{4\alpha }A=\frac{S_n}{4\alpha {\pi }^{n/2}},$$

$$T(x,t)$$

$$\hat{f}(\omega )=(\sqrt{2+\sqrt{2}}+\mathrm{i}\sqrt{2-\sqrt{2}})\frac{\sqrt[4]{\omega }}{\sqrt{\left|x\right|}}{K}_{1/2}(2\sqrt{-\mathrm{i}\omega }\left|x\right|),$$

$$K_{1/2}(z)=\sqrt{\frac{\pi }{2}}\frac{e^{-z}}{\sqrt{z}}.$$

$$\hat{f}(\omega )=\frac{\sqrt{\pi }}{\sqrt{2}\left|x\right|}(1+\mathrm{i})e^{-\sqrt{2\omega }|x|}\left(\cos (\sqrt{2\omega }\left|x\right|)+\mathrm{i}\sin (\sqrt{2\omega }\left|x\right|)\right),$$

$$T(x,t)$$

$$T(x,t)\approx \frac{1}{4\alpha \pi \left|x\right|}\sum _{n=1}^N{\tilde{S}}_n\cos (2\pi f_{n}t+\widetilde{{\varphi }_n}),$$

$${\tilde{S}}_n$$

$${\tilde{\varphi }}_n$$

$$s(t)=\cos (\omega t+B\sin ({\omega }_{m}t))=\sum _{k=-\infty }^{\infty }{J}_k(B)\cos ((\omega +k{\omega }_m)t),$$

$$s(t)$$

$$\hat{M}$$

$${\hat{M}}_0$$

$$F(x)$$

Full text

1. 1 Introduction

2. 2 Mathematical model

3. 2.1 About Pennes’ bio-heat model

4. 2.2 Static full-space solution

5. 2.3 Dynamic full-space solution

6. 2.4 Some examples

7. 3 Literature review

8. 4 Algorithms

9. 4.1 General algorithm descriptions

10. 4.2 Source power reconstruction

11. 4.3 Source spectrum reconstruction

12. 4.4 Static reconstruction algorithm

13. 4.5 Dynamic reconstruction algorithm

14. 5 Simulation

15. 6 Numerical tests

16. 6.1 Static reconstructions

17. 6.2 Dynamic reconstructions

18. 7 Conclusions

1. Introduction

This work stems from the field of medical infrared (IR) imaging, where the boundary temperature of an object is measured using an IR camera and subsequently analysed using e.g. statistical methods or thermodynamic modeling. More specifically, the starting point of this paper is the method of dynamic infrared imaging (DIRI), where time-varying temperature phenomena are recorded with framerates of 30Hz and higher. Thus we are dealing with a time-varying temperature distribution on the boundary of the object, from where we are interested in computing some thermal properties inside the object e.g. heat conductivity, blood perfusion or the heat source distribution. In this paper we deal with the last case in 3D: given the temperature measurement $T(x,t)|_{\partial\Omega}$ on the boundary of the domain $\Omega\subset{\mathbb{R}}^{3}$ in a given time interval, reconstruct the heat source

[TABLE]

where $\delta(x_{0})$ refers to a point source (delta-distribution) at the spatial location $x_{0}$ and $s(t)$ is a signal of the form

[TABLE]

We specifically investigate two problems:

1. (1)

Given the static measurement $T(x)|_{\partial\Omega}$, reconstruct the location $x_{0}$ and magnitude $Q$ of a static heat source $q(x)=Q\delta(x_{0})$. 2. (2)

Given the dynamic temperature measurement $T(x,t)|_{\partial\Omega}$, reconstruct the location $x_{0}$, the amplitudes $S_{n}$ (and the phases $\phi_{n}$ should we need them).

Mathematically speaking, reconstructing the heat source from boundary measurements is an ill-posed inverse problem, which is very sensitive to noise in the measurement and thus needs a regularized solution method. Inverse problem algorithms can be categorized, roughly speaking, in three: analytic or direct methods, iterative optimization methods, and statistical methods. In this paper we develop an analytical model which in principle gives a direct analytical reconstruction formula, but which will be too unstable to use directly. So we will use the analytical model to compute a penalty functional from the noisy measurement data in a parameter space. We do not need optimization algorithms since the space can be fully computed and the minimizer even visually inspected.

This particular problem setting with time-dynamic heat sources (1) in 3D may not be very common in scientific literature, in 1D similar ideas are used in the field of material testing using Pulsed Phase Tomography algorithms [1, 2]. In our case it is loosely based on the application of breast cancer detection using Fast Fourier Trasform (FFT) as described in [25, 26] and in [27, 28, 29, 30]. In those papers, by considering the measurement $T(\cdot,t)|_{\partial\Omega}$ a stationary signal and computing its FFT it was hypothized that the frequency spectrum would then reveil abnormal thermodynamic behaviour, for example because of increased blood pulsation and vasodilation. In our paper we build on this idea and investigate how the heat source variation presents itself on the boundary. It turns out that each frequency component $S_{n}\cos(2\pi f_{n}+\phi_{n})$ of a point heat source results to a temperature signal of the same frequency, but attenuated and phase-shifted according the the frequency.

The main novelty in our approach is the possibility to reconstruct a seemingly complicated heat source with low computational cost, consequently we can plot the complete penalty functional in several dimensions and see how the nature of the inverse problem presents itself there. Another option with potentially less approximation error would be to use Finite Element Method (FEM) or other numerical solvers to compute the forward model for any heat source. Since this can be computationally expensive, to optimize our penalty functional we would use some iterative scheme and hope that the obtained minumum is not local, but global. The present work serves as a proof of concept of the idea of using point source models as building blocks; in the future we aim to investigate linear forward models built from an array of point sources and then the robust inversion, for example using truncated singular value decomposition, to reconstruct source information from boundary temperature measurements.

Unfortunately it would be too much to investigate all of these matters thoroughly in one simple scientific article, and thus we have to make approximations and simplifications, some of which seem severe. We assume that the object is homogenic in material, so that the specific heat capacity, density and the heat conductivity is constant. We don’t simulate the IR-sensor technology and their specific noise profiles, different in space and time dimensions, which will always be present in real world measurements, but instead add simulated gaussian (white) noise to the simulated measurements. Also we make use of analytical solutions of the heat equation in the full space, whereas in reality we are always dealing with finite objects and subsequently with certain boundary conditions. Moreover, we do not investigate the actual bio-physical models nor the phenomena they formulate.

We continue with the presentation of our analytical model in section 2. Then we give a brief review of relevant scientific work in section 3. After, we describe the reconstruction algorithm in section 4. The simulation of data is described in section 5. In section 6 we present our numerical results, and finally conclude the paper in the final section 7.

2. Mathematical model

Define the heat equation by

[TABLE]

where $c(x)$ is the volumetric heat capacity, $\kappa(x)$ is the heat conductivity, $q(x,t)$ is the volumetric heat source and $T(x,t)$ is the temperature distribution ${\mathbb{R}}^{3}\times{\mathbb{R}}\to{\mathbb{R}}$. In this document the source has the form (1) and we first aim to (analytically) compute $T(x,t)$, later in the detection algorithms we aim to compute the location and signal content of $q(x,t)$ from the knowledge of $T(x,t)$. Note that in real life $T(x,t)$ emerges automatically as it is the result of the natural heat diffusion process.

We might not have an analytical solution for the temperature distribution $T(x,t)$ in the most general case. To simplify, we assume that the coefficients $\kappa,c$ and thus the diffusivity $\alpha=\kappa/c$ are constant, and abuse notation a bit to write the heat source

[TABLE]

and consider

[TABLE]

Then, using standard analytical tools such as the fundamental solution (Gaussian kernel)

[TABLE]

where $n=3$ for 3D problems, and Duhamel’s principle, we get a solution of (5) with zero initial condition,

[TABLE]

where $s(y)$ is the time-varying signal of the heat source (4). We may also change the order of convolution, to get the later useful formula

[TABLE]

An important property of heat equations with constant diffusivity $\alpha$, such as (5), is the following: let $T_{1},T_{2}$ be solutions to

[TABLE]

Then it is clear that $T=T_{1}+T_{2}$ solves

[TABLE]

This means that by solving analytical solutions $T_{n}$ for heat sources $q_{n}$ we also get the solution corresponding to the source $\sum_{n}q_{n}$ by a simple summation of the solutions. In particular this principle can be used to combine analytical solutions of multiple separate point sources, and also each heat signal source can be separated into its harmonics, each of which constitutes as a single (time varying) heat source. Moreover this means that given an added time-constant source,

[TABLE]

the resulting time-dynamic temperature profile can be computed by subtracting the constant temperature solution.

Of course all of the above mentioned solutions exist only in certain function spaces and with certain theoretical considerations not explained here. What now follows is some extra considerations on how the simple equation (5) differs from more realistic models in section 2.1 and what kind of full-space solution formulas derived from (7) exist in sections 2.2 and 2.3.

2.1. About Pennes’ bio-heat model

According to the well known Pennes’ bio-heat transfer model the temperature $T$ in biological tissue satisfies

[TABLE]

where $\rho,c$ and $k$ are the density, thermal capacity and heat conductivity of the tissue respectively, $w_{b},\rho_{b}$ and $c_{b}$ are the perfusion, density and thermal capacity of the blood respectively, $T_{a}$ is the temperature of arterial blood and $Q_{m}$ is the metabolic heat source. When compared to (3), the difference is the form of the source term. Of course in any specific situation a model using a source of the general type $q(x,t)$ might be a good approximation, but in general the convection term $w_{b}\rho_{b}c_{b}(T_{a}-T)$ is needed to model the blood circulation and its effect on the heat dynamics. The Pennes’ bio-heat model and its modified versions are well studied, see for example [31, 32, 33, 34].

We should also consider finite objects, say domain $\Omega\subset{\mathbb{R}}^{3}$ with a boundary $\partial\Omega$, and the full boundary value problem with Neumann-type boundary condition, so that the temperature distribution solves (9) inside $\Omega$ and satisfies

[TABLE]

where $h$ is the heat transfer coefficient and $\nu$ is the outwards unit normal vector. This type of boundary condition is usually found in conjunction with the Pennes’ model. It is not clear in what way this type of condition will ”affect” the full-space solution, or in other words how the solution in a finite domain satisfying (10) differs from the full space solution of (3) or (5).

In the simulations and numerical tests of this paper we use the model (5); analysis of the error caused by the omission of the convection term and the boundary condition in real world applications is left for future works.

2.2. Static full-space solution

Assume a static heat source (4) with $s(t)=Q\delta(x_{0})$ (and abuse the notation again to infuse the term $c$ into $Q$). We now follow the method described in [35, 36]. Take the equation (5) in steady-state and in spherical coordinates, with the origin set at the point source location $x_{0}$, to get

[TABLE]

where $r$ is the distance from $x_{0}$. The solution will be spherically symmetric. The Laplace-operator simplifies to only include $r$-terms, and we get

[TABLE]

This differential equation is easily solved for its general solution

[TABLE]

where $C$ is a constant. In [35] this analysis is followed by a method to then compute the unknown $Q$ and the depth $d$, defined as the shortest distance between $x_{0}$ and the boundary, from the boundary temperature distribution $T(x)|_{\partial\Omega}$. In [36] another method is presented, using curve fitting to plot the ”Q-r” -curve which according to their clinical study has different profile for different diagnostic cases. Neither of these methods are using the usual, regularized mathematical machinery from the field of inverse problems, and so we continue with our own implementation in section 4. One aspect we take care of is to not only search for the depth of the heat source, but the location as well, since it might not be clear where the point of maximum temperature is located exactly in a noisy measurement.

It should be noted here that in [35] there is also a step to filter the temperature distribution to account for non-singular heat sources. We omit this step, and just simply test our resulting algorithm in a simulation where the heat source is not as small as possible, but a decently sized sphere.

In figure 1 we plot the distinguishability of an abnormal heat source using different values of diffusivity. The heat source is spherical with values $Q=29000W/m^{3}$ and radius $R=1cm$, mimiking the metabolic rate of a cancerous tumor. The distinguishability here means the difference of temperature between the hottest spot and 5cm away on a plane. In real applications the distinguishability will be less because of noise and boundary conditions such as (10). Moreover the distinguishability here is a crude approximation of what it would be like to actually perform an analysis on the boundary temperature data, which entails the full temperature data between the hottest spot and the nearby area, up to 5cm radius. Fat and muscle tissue have diffusivity of 0.98E-7 and 1.31E-7 respectively, which will be suitable for the static detection of heat sources, based on figure 1. In the next section we will consider dynamic heat sources which have larger optimal diffusivities depending on the frequency of the source.

2.3. Dynamic full-space solution

The idea of what follows can be found in the book [38] and is also available in the lecture note [37]. Assume the time-variable heat source (4) with one harmonic component,

[TABLE]

Then, by writing

[TABLE]

from (8) and trigonometric identities, as $t\to\infty$, we obtain

[TABLE]

where $\hat{f}$ is the half-line Fourier-transform of $f$. In 3D we have

[TABLE]

where $K_{1/2}$ is a modified Bessel function

[TABLE]

Thus we get

[TABLE]

and subsequently, using again some trigonometric identities,

[TABLE]

Note that the original heat signal source will inhibit phase change and amplitude attenuation depending on the frequency and distance from the source; in particular the amplitude attenuation is exponential.

Take the full heat source of (1). By using the summation property and (13), the corresponding analytical solution of the temperature distribution is

[TABLE]

where

[TABLE]

Many specific heat sources can be divided into its components in a similar way, and their respective temperature distributions deduced using the harmonics. Some of these are presented in the next section.

Note the different attenuation terms in (13). The term $1/\left|{x}\right|$ is the same for all frequencies as in the static case (12). The exponential term depends also on the frequency. Both terms depend on the diffusivity. As in figure 1, it is interesting to see the distinguishability for different diffusivity coefficients and frequencies. The scale of the depth and diffusivity in figures 2, 3 and 4 are different from the static case because the attenuation is much more severe because of the exponential term. In this dynamic case, it is more difficult to estimate what the distinguishability, computed in the same way as in the static case, means. In principle the amplitude information is computed using FFT and if the phenomena are perodic and we have long enough acquisition time, we are able to detect very small temperature variation amplitudes.

2.4. Some examples

Consider a phase-modulated source and its series expansion

[TABLE]

where $\omega_{m}$ is the frequency of modulation and $J_{k}(x)$ is a Bessel function of the first kind, of order $k$. Following earlier reasoning, each component in the sum (17) contributes to a solution given by (13); each component will be attenuated and phase-shifted according to their frequency, higher frequencies more than lower. The source term itself has main frequency $\omega$ and symmetrical sidelobes of frequencies $n\omega_{m}$ with amplitudes $J_{n}(B)$. The smaller the modulating amplitude $B$ is, the smaller these sidelobes are; when $B\geq 1$ the sidelobes become significant and will tail off as a sinusoid in their amplitude. The important thing is that for larger amplitude, say $B=2$, the amplitude of the main frequency will be lower than the first sidelobes, and in general for even larger $B$ the largest harmonic components are (symmetrically) further away from $\omega_{c}$. Thus without knowing the attenuation behaviour described by equation (13), we could register a lower frequency from the temperature variation than the actual source frequency.

Consider an amplitude modulated source and its trigonometric expansion

[TABLE]

The first term can be thought of as our signal of interest, the second introduces sidelobes situated $\omega_{a}$ away from the original, with amplitudes $M/2$. Again as with the phase-modulated source, each of these components contributes to a solution given by (13).

3. Literature review

The heat source reconstruction problem has been studied in many aspects: in spatial dimensions 1,2 and 3, with various forms of the heat source, with constant and variable heat conductivity and with various assumptions on the measurement, both Dirichlet and Neumann conditions on either partial or full boundary. Many papers exist on the stability of the problem, that is an estimate on how the error between the reconstructed solution and the true solution behaves as a function of the measurement error. Some papers concentrate on the mathematical theory, others include some numerical tests and even numerical algorithms.

In [3] analytical solutions of the three dimensional heat equation, using a series expansion, are computed for a multi-layered sphere. In [4] analytical solutions are used to model moving heat sources i.e. welding process. Both of these papers could be the inspiration for an even more advanced analytic modeling and subsequent heat source detection.

The papers [5, 6, 7] deal with one dimensional numerical heat source reconstruction, but they do not have any regularization even though the problem is ill-posed. A more mathematical background dealing with uniqueness and establishing the problem as an inverse problem is given by J. R. Cannon in [8, 9] and in numerous other papers by the same author. Continuing the analytical approach we have the theoretical stability analysis of papers [16, 17, 18, 19, 20, 21]. Examples of regularization strategies to battle against the ill-posedness are found in [10, 11, 12, 13, 14, 15].

Finally, in [22] we have a regularized, numerical algorithm to reconstruct time-dynamic heat source in 2D, for the purposes of material testing and the analysis of thermomechanical effects. The method is further developed in [23].

Please see also the references in the publications above, as the list presented here is only a very limited and somewhat arbitrary subset of all of the work done in this field.

4. Algorithms

We can think of several classes of heat source reconstruction algorithms, from the most general case to special cases using a priori information. In order to describe the algorithms we need more notation. The boundary of the object is $\partial\Omega=B(x)$. The ”full” measurement is $M(x,t)=T(x,t)|_{B(x)}$, where $t\in[0,t_{max}]$. From $M$ we may manually choose a more interesting neighbourhood $M_{0}(x,t)$ where $x\in\mathcal{N}\subset B$, containing a visible hot spot we want to analyze in more detail. The FFT of these measurements are $\hat{M}$ and $\hat{M}_{0}$, for each spatial point $x$, more specifically we write

[TABLE]

where the amplitudes $A_{n},A^{0}_{n}$ and phases $\theta_{n},\theta^{0}_{n}$ for each $n$ describe the signal component corresponding to the frequency $f_{n}=n/N_{s}\cdot f_{s}$, where $f_{s}$ is the sampling rate and $N_{s}$ is the number of samples. The index $n=1,\ldots,N$ is meant to go through the interesting frequencies, if not all available from FFT, in which case we would have roughly $N=N_{s}/2$ because of Shannon-Nyquist sampling theorem.

For each point $x_{0}$ and frequency $f_{n}$ we can write

[TABLE]

and subsequently we can express the diffused source component of (13) as

[TABLE]

where

[TABLE]

Also the static solution (12) can then be written as

[TABLE]

There are many ways to arrange the detection algorithm in terms of whether it is model or data driven. In our algorithms, we take a grid of source location candidates $x_{0}$ and output the one having the smallest penalty functional. But to compute the penalty, we first use the model and the data together to compute some intermediate steps. In the static case, we compute a maximum and minimum temperature value on $M_{0}$ from which we get the model complient parameters $Q,C$ for that particular source location. Then, the penalty functional is

[TABLE]

In the dynamic case, for each candidate point source $x_{0}$ we compute its nearest boundary point $C$ on $M_{0}$ and then the values $\{A^{0}_{k},\theta^{0}_{k}\}$ for each frequency component $k$. Using the dynamic model (25) these are then turned to reconstructions $\{S_{k},\phi_{k}\}$. However, we both cut away impossible spectral content and find the best fit from the source locations. By defining truncation parameters $A_{t}$ and $M_{t}$, we define the set $M_{0}^{k}\subset M_{0}$ as the set for which has spectral amplitude over $A_{t}$ (cutting off too small data) and which has the size of at least $M_{t}$ (cutting off too little pieces of measurement data). The lower $A_{t}$ and higher $M_{t}$ is, the more robust we need elsewhere in the algorithm. For each source candidate, for the amplitude data, we use the penalty functional

[TABLE]

where the profile $\tilde{S}_{k}$ is given by the source candidate $x_{0}$, the computed value $A^{0}_{k}$ and the dynamic model (25). The total amplitude penalty functional is

[TABLE]

For the phase penalty functional, we use

[TABLE]

where $\tilde{\phi}_{k}^{1}$ is the model computed profile (at each source location $x_{0}$) normalized to have the same scale as the measurement $\theta_{k}$ on $M_{0}^{k}$. Similarly to the amplitude penalty we write

[TABLE]

We were not able to combine the information given by the amplitude and phase profiles in a good way, so instead we get two reconstructions using either information.

Lastly, we denote the normalized collection of harmonics by $\{S^{1}_{k},\phi_{k}\}$, satisfying $S_{k}=AS_{k}^{1}$ for each $k=1,\ldots,N$, where $A$ is the amplitude of the signal with harmonic content $\{S^{1}_{k},\phi_{k}\}$.

4.1. General algorithm descriptions

There are many possibilities how to arrange the algorithm to find the location $x_{0}$ and the spectral content $\{S_{n},\phi_{n}\}$ from the measurement, depending also on what we expect the situation to be. We can think of the following examples:

• (A1)

Given $M$, reconstruct $q(x,t)$.

• (A2)

Given $M$, reconstruct a collection of point sources $\{x_{0}\}_{m}$ with signals $\{S_{n},\phi_{n}\}_{m}$.

• (A3a)

Given $M_{0}$, reconstruct $Q,C,x_{0}$.

• (A3b)

Given $M_{0}$, reconstruct $x_{0}$ and $\{S_{n},\phi_{n}\}$.

• (A4)

Given $M_{0}$ and $\{S^{1}_{n},\phi_{n}\}$, reconstruct $x_{0}$ and $A$.

The algorithm (A1) will be difficult to realize and is not in the scope of this paper. Our model only considers point-like sources. In principle by using the linearity of the heat equation with constant diffusion constant one might be able to approximate different geometries of heat sources by using a collection of point sources. A moving pulse-like heat source can be approximated by a series of point sources, but then an issue of phase cancellation emerges, when two different points sources produce two temperature variations in different phase resulting to reduced temperature signal. Part of the algorithm (A1) would be to reconstruct the geometry of moving sources from the phase information.

To understand the difficulty of algorithm (A2), consider the fact that each point source will be visible in the part of $M$ closest to them, as a hot spot from which we want reconstruct the parameters of the point source. The algorithm (A2) would need an automatic scan of the full measurement $M$ in such a way that these visible hot spots are first scanned, after which they are analyzed separately. Again, we do not attempt this, but instead assume that the user is able to choose an interesting neighbourhood which contains a hot spot, to reconstruct each heat source $x_{0}$ of interest.

The algorithm (A3a) uses the static model (25) and is described in more detail in the next section. The algorithm (A3b) uses the dynamic model (22) and is described in section 4.5.

The algorithm (A4) is thought of as a simple way of using a priori information about the signal, for example we could have another means of measuring the heat source signal content in conjunction with the temperature profile on the boundary after which the algorithm would search for this signal content in the temperature data. The details for this is left in future works.

4.2. Source power reconstruction

Note that a finite domain has a boundary and boundary conditions. If there are positive heat sources and perfect insulation at the boundary, so that the domain would not lose heat to the surroundings, then the temperature of the domain would increase infinitely. A more realistic scenario, which would be true for all objects in the physical world, is a boundary condition such as equation (10) which would result to a stable state after a sufficient stabilization time. The stable temperature of the domain would be dependent on the values of the heat source power and the heat transfer coefficient on the boundary. Having a larger power heat source and larger heat transfer coefficient could result to approximately the same average temperature as smaller power heat source and smaller heat transfer coefficient; thus our analytical model which takes none of this into account can never perform an absolute reconstruction of the heat power.

However, using the equation (25) we may still compute an absolute value for the power amplitude $Q$, even if it won’t be a physically correct value. Subsequently our model can still compare different powers of heat sources since they can be reconstructed having different values of $Q$. To do this we use the following sub-algorithm (AQ) for each candidate location $x_{0}$:

1. (1)

Sort the values of $M$ or $M_{0}$ and take an average of the first and last $\mathrm{Nav}$ values to get the minimum temperature $T_{1}$ and maximum temperature $T_{2}$ of the measurement patch. Besides the sorted vector of measurements, use their locations to compute the corresponding average locations $x_{1}$ and $x_{2}$ and their distances $r_{1}$,$r_{2}$ from $x_{0}$. 2. (2)

Since we seek for values $Q,C$ for the model (25), from there we get the formula

[TABLE] 3. (3)

Compute the values $C_{1},C_{2}$ again from the equation (25) using $T_{1},T_{2},r_{1},r_{2}$ and $Q$, and finally compute

[TABLE]

In this way we can compute for any candidate $x_{0}$ the power $Q$ and the parameter $C$. If $M$ or $M_{0}$ is very noisy the above will perform less well, but the integer $\mathrm{Nav}$ can be used to alleviate this problem a little.

4.3. Source spectrum reconstruction

The case of a dynamic heat source becomes more complicated as we have one more dimension in data and source parameters. The biggest problem is the sensible use of the phase information; the model could easily predict a total phase change of $6\pi$, for example, but the FFT processing of any simulated or real data will not be accurate throughout the measurement patch, but instead only inside a small radius. In addition to this the overall level of the phase error norm will be different form the amplitude level, thus making it difficult to combine the phase and amplitude information. The algorithm presented here is only one somewhat simple suggestion and should be investigated more in a subsequent paper.

One of the features is that we automatically cut away frequency information that is not large enough and only reconstruct the best source fit for the remaining data. In some future work when tackling more noisy data this approach will probably be too limiting. Other feature that we had to reluctantly settle with is that the algorithm finds the best fit separately using amplitude and phase information since we did not find a good way to combine the differently behaving penalty functionals.

In the dynamic case the measurement $M$ or $M_{0}$ has been processed by FFT to get the collection $\hat{M}$ or $\hat{M_{0}}$ of amplitudes and phases. Define an amplitude limit $A_{t}$ below which we discard amplitude data. Then, given any source location candidate $x_{0}$ we use the following sub-algorithm (AA):

1. (1)

Compute the boundary point $C$ closest to the source and write $A_{n}^{C}$ for the amplitude values on a neighbourhood of $C$, for each frequency bin $n$. Write $r=\left|{x-x_{0}}\right|$. 2. (2)

For each $n$, use (22) and $A_{n}^{C}$ to compute the model predicted amplitude profile $\tilde{A}_{n}=A_{n}^{C}F(r)e^{-D_{n}(r)}$. Compute the sub-measurement patch $\tilde{M}_{n}=\{x\in B(x):\quad\tilde{A}_{n}>A_{t}\}$. 3. (3)

If $\tilde{M}_{n}=\emptyset$ or a very small set, put $S_{n}=\phi_{n}=0$ as the reconstruction. Otherwise use $\tilde{M}_{n}$ as the domain where the phase error is computed.

The phase profile $\phi_{n}$ needs special processing as its raw values are always inside $[0,2\pi)$. We use Matlab’s unwrap-function to first make the phase profile continuous, then put the phase value on $C$ to zero by subtracting the value computed earlier in conjunction with the amplitudes $A_{n}^{C}$. The analytical model will produce a reference phase profile that is already continuous, and we process it the same way so that on the boundary point $C$ the phase will be zero. The phase information can be averaged by first changing to unit vectors in 2D, then averaging the vectors and finally changing back to the angle of the averaged vector.

Since we only use $\tilde{M}_{n}$ to compute the phase error, two problems arise. First, the smaller the area of $\tilde{M}_{n}$, the smaller the magnitude of the phase error. Second, the phase information outside of $\tilde{M}_{n}$ will not contribute to the error, which is clearly not optimal. The first issue is somewhat solved by dividing the error norm by the area of $\tilde{M}_{n}$. The second issue is not properly addressed, but instead we compute the optimal source location and the accompanying spectral content using the amplitude and the phase information separately.

4.4. Static reconstruction algorithm

The algorithm (3a) in more detail is the following:

1. (1)

From $M$, pick an interesting neighbourhood $M_{0}$ manually. 2. (2)

Take a collection of source candidates $x_{0}$ and for each of them use (AQ) to compute the parameters $Q$ and $C$. For each triplet $(x_{0},Q,C)$ compute $\epsilon$ from (26). 3. (3)

The reconstruction is the triplet having the smallest error $\epsilon$.

4.5. Dynamic reconstruction algorithm

The algorithm (3b) in more detail is the following:

1. (1)

From $M$, pick an interesting neighbourhood $M_{0}$ manually. 2. (2)

Compute FFT of $M_{0}$ to get the amplitude and phase data. 3. (3)

Take a grid of source candidates $x_{0}$ and for each use the sub-algorithm (AAn) to get the model-conforming amplitudes $A_{n}^{C}$ and usable phase areas $\tilde{M}_{n}$. 4. (4)

For each $x_{0}$ and each $k$ that is not cut away, compute amplitude and phase penalty functionals(27) and (29). 5. (5)

The reconstruction is the location $x_{0}$ and the accompanying signal $\{S_{k},\phi_{k}\}$ using either (28) or (30).

5. Simulation

For medical imaging, ideally we would want to simulate data on a realistic 3D geometry, using realistic values of heat conductivity and capacity for anatomically correct tissues. The discrete computing mesh should enable high accuracy in the parts of the domain where we measure physical properties, but should be adaptively course to keep down computational demands. Unfortunately in this paper we needed to keep things simple. The domain is a half ball of radius 100, the computational mesh is similar everywhere in the domain and has maximum edge length of 10 for the static examples, 14 for dynamic examples, the diffusivity is constant 2 and the values of heat generation and heat conduction outwards are chosen so that we get somewhat realistic (between 20 and 40 degree Celsius) temperature values inside the object.

We investigate static and dynamic ball-shaped heat sources of radius 10. For the static cases, there are three different locations of the heat source each with four different volumetric power. For the dynamic case we only compute two locations closest to the boundary and two even larger amplitudes of power modulation compared to the static case, because the effect of attenuation is severe.

We use Finite Element Method (FEM) to compute the temperature profile in three spatial dimensions in addition to the time dimension. For the static cases we only use the temperature profile at the last available discrete time position. For the dynamic cases we use a time-dynamic profile of the last ten seconds of the simulation. In general we computed several rounds of FEM solutions in a serial manner, so that the temperature profile at the end of one simulation was used as the initial value for the next simulation. In this way we could use a course time dimensional computation for the stabilization period, and finer time dimensional computation at the very end after the system had stabilized.

The dynamic solution is transformed to amplitude and phase information using FFT with frame rate of 10Hz without windowing since it was not needed in our particular choices of frequency and acquisition time. Before FFT we also removed linear trend of warming from the simulated measurement data.

Accuracy problems were encountered because of the four dimensional computations and subsequent need for course FEM mesh in spatial dimensions. We interpolated the temperature profile at radius 90 in the static examples and at radius 85 in the dynamic examples, inside the domain, as the boundary measurement. The dynamic simulations were problematic especially since the dynamic effects are exponentially attenuated, and thus very small in general.

We simulate noise by adding Gaussian white noise at several levels. One percent noise means the variance of the Gaussian distribution is one percent of the maximum temperature difference in the measurement patch $M_{0}$.

The parameters of the heat sources are given in table 1. In figures 5 and 6 we have an example of manually picked temperature measurement patch $M_{0}$ near the boundary, non-noisy and with three different simulated noise levels, from static cases A1 and C2. The two figures use the same temperature scale and colormap and shows how much smaller and more difficult to detect the ”hot spot” becomes as the heat source is located deeper in the domain and has a smaller power.

6. Numerical tests

We used truncation values $A_{t}=0.02$ and $M_{t}=900$ pixels. In the static case we took $40^{3}$ different source candidate locations, 40 in each spherical dimension to cover the whole domain $\Omega$ near the measurement patch $M$. In the dynamic case we used $30^{3}$ different source candidate locations, giving inherently a bit smaller detection resolution.

6.1. Static reconstructions

In table 2 are the results of the static simulations. The location column has the distance from the reconstructed source center to the true source center in same units as the simulations. The power columns give the reconstructed power compared to the non-noisy reconstruction of source A1, which should be the easiest to reconstruct. In figure 7 we see example reconstructions, or the best fit temperature profiles, to show how close the point source model computed measurements are to the simulated measurements.

6.2. Dynamic reconstructions

In tables 3 and 4 are the results of the dynamic simulations. The location columns have the distance from the reconstructed source center to the true source center in same units as the simulations, both using the amplitude and phase information. The spectrum columns give the reconstructed source modulation amplitude compared to the non-noisy reconstruction of source A1, which again should be the easiest to reconstruct. In figures 8 and 9 we see examples of how close the model computed FFT amplitudes are to the simulated FFT amplitudes. What is not presented here is the full reconstructed spectra; the cutoff parameter $At$ was efficiently chosen so that all other frequency components other than the third FFT bin corresponding to exactly $f=0.2Hz$ were zero, in all reconstructions.

7. Conclusions

In this article we have presented a new approach to detecting static and time-varying point-like heat sources. The numerical work here can be thought of as a proof of concept of this approach, even if there is a lot of work to do in order to develop the algorithms, especially as the data changes from simple and simulated to realistic and noisy. This simple model can in principle be used to approximate more complicated cases.

The numerical simulations turned out to be very challenging as we would need high numerical accuracy in 4D. In the future we hope to revisit the simulations with realistic parameters. As can be seen in figures 8 and 9 on the left, the simulations are not looking smooth - this could be improved by using a better and finer FEM mesh and possibly by adjusting FEM solver parameters. This is a challenging and laborous task in itself.

On the other hand, the example model fits of figures 7,8 and 9 are good enough to label this approach as promising. The simulation was done on a finite half ball with a semi-realistic boundary condition used in real bio-heat models; still the analytical point-source model using a solution in the infinite 3D space fits reasonably well. Another future project would be to further investigate ways to either modify our model to accommodate realistic features such as, for example, blood convection, or to analytically compute the different errors it induces when used in realistic scenarios. For example, the boundary condition (10) could perhaps be addressed analytically as well as different source shapes. In fact, any source could be thought of as a collection of point sources, each contributing as a simple model solution of (25) and (22).

One very interesting future research subject is the efficient use of phase information. As the model (25) predicts, the phase profile will be not attenuated exponentially, but linearly, although the phase is related to the amplitude and we will not be able to compute it more accurately from the data.

The results of tables 2, 3 and 4 have to be taken critically because of the accuracy problems of the simulations and inherent complexity of the algorithms. Still, one could argue that the basic trends of an inverse problem are present: the more deeply situated source, the more difficult it is to reconstruct, and the more noise, the results get worse quickly. The dynamic results might look perplexing, but we think that the following explains them:

• •

Using the amplitude information, the location detection is almost unchanged as noise is increased, in all of the simulations. What is not shown here is the results for higher noise levels, from around 5% onwards the reconstruction break down very suddenly. Also the noise presents itself differently as we compute FFT from the noisy data. This is another issue completely and is not covered here, but in general it is possible to detect very small periodic phenomena from noisy data if the sampling time is large. In our simulation the sampling time was quite average ten seconds.

• •

For some of the simulations, especially using the phase information, the reconstructed spectrum amplitude value change radically. These values are explained by a best fitting source location around 10 units towards or away from the true location, and so the source amplitude is then computed with or without this extra attenuating layer. This behaviour really highlights the ill-posedness of the problem, and to make the algorithms more usable we would need more robust handling of data.

• •

The dynamic results of the more deeply located source B are better than of source A and even the static cases (based on the location error). This can be explained by the fact that the more deeply located source is actually more like a point source and thus the analytic model could work better in that case. Also, the dynamic cases use a lot more data for the source detection, and so if the FFT can be computed properly like in our simulations, it makes sense that the location is computed more accurately.

• •

The reconstructions using the phase information might suffer from the more difficult algorithm implementation and also from larger deviation between the simulation and the analytical point source model, as the phase profile is not so attenuated.