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 timevarying temperature phenomena are recorded with framerates of 30Hz and higher. Thus we are dealing with a timevarying 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 on the boundary of the domain in a given time interval, reconstruct the heat source
(1) 
where refers to a point source (deltadistribution) at the spatial location and is a signal of the form
(2) 
We specifically investigate two problems:

Given the static measurement , reconstruct the location and magnitude of a static heat source .

Given the dynamic temperature measurement , reconstruct the location , the amplitudes (and the phases should we need them).
Mathematically speaking, reconstructing the heat source from boundary measurements is an illposed 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 timedynamic 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 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 of a point heat source results to a temperature signal of the same frequency, but attenuated and phaseshifted 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 IRsensor 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 biophysical 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
(3) 
where is the volumetric heat capacity, is the heat conductivity, is the volumetric heat source and is the temperature distribution . In this document the source has the form (1) and we first aim to (analytically) compute , later in the detection algorithms we aim to compute the location and signal content of from the knowledge of . Note that in real life emerges automatically as it is the result of the natural heat diffusion process.
We might not have an analytical solution for the temperature distribution in the most general case. To simplify, we assume that the coefficients and thus the diffusivity are constant, and abuse notation a bit to write the heat source
(4) 
and consider
(5) 
Then, using standard analytical tools such as the fundamental solution (Gaussian kernel)
(6) 
where for 3D problems, and Duhamel’s principle, we get a solution of (5) with zero initial condition,
(7) 
where is the timevarying signal of the heat source (4). We may also change the order of convolution, to get the later useful formula
(8)  
An important property of heat equations with constant diffusivity , such as (5), is the following: let be solutions to
Then it is clear that solves
This means that by solving analytical solutions for heat sources we also get the solution corresponding to the source 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 timeconstant source,
the resulting timedynamic 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 fullspace solution formulas derived from (7) exist in sections 2.2 and 2.3.
2.1. About Pennes’ bioheat model
According to the well known Pennes’ bioheat transfer model the temperature in biological tissue satisfies
(9) 
where and are the density, thermal capacity and heat conductivity of the tissue respectively, and are the perfusion, density and thermal capacity of the blood respectively, is the temperature of arterial blood and 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 might be a good approximation, but in general the convection term is needed to model the blood circulation and its effect on the heat dynamics. The Pennes’ bioheat model and its modified versions are well studied, see for example [31, 32, 33, 34].
We should also consider finite objects, say domain with a boundary , and the full boundary value problem with Neumanntype boundary condition, so that the temperature distribution solves (9) inside and satisfies
(10) 
where is the heat transfer coefficient and
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 fullspace 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 fullspace solution
Assume a static heat source (4) with (and abuse the notation again to infuse the term into ). We now follow the method described in [35, 36]. Take the equation (5) in steadystate and in spherical coordinates, with the origin set at the point source location , to get
(11) 
where is the distance from . The solution will be spherically symmetric. The Laplaceoperator simplifies to only include terms, and we get
This differential equation is easily solved for its general solution
(12) 
where is a constant. In [35] this analysis is followed by a method to then compute the unknown and the depth , defined as the shortest distance between and the boundary, from the boundary temperature distribution . In [36] another method is presented, using curve fitting to plot the ”Qr” 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 nonsingular 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 and radius , 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.98E7 and 1.31E7 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 fullspace solution
The idea of what follows can be found in the book [38] and is also available in the lecture note [37]. Assume the timevariable heat source (4) with one harmonic component,
Then, by writing
from (8) and trigonometric identities, as , we obtain
where
is the halfline Fouriertransform of
. In 3D we havewhere is a modified Bessel function
Thus we get
and subsequently, using again some trigonometric identities,
(13)  
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
(14) 
where
(15)  
(16) 
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 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 phasemodulated source and its series expansion
(17) 
where is the frequency of modulation and is a Bessel function of the first kind, of order . Following earlier reasoning, each component in the sum (17) contributes to a solution given by (13); each component will be attenuated and phaseshifted according to their frequency, higher frequencies more than lower. The source term itself has main frequency and symmetrical sidelobes of frequencies with amplitudes . The smaller the modulating amplitude is, the smaller these sidelobes are; when the sidelobes become significant and will tail off as a sinusoid in their amplitude. The important thing is that for larger amplitude, say , the amplitude of the main frequency will be lower than the first sidelobes, and in general for even larger the largest harmonic components are (symmetrically) further away from . 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
The first term can be thought of as our signal of interest, the second introduces sidelobes situated away from the original, with amplitudes . Again as with the phasemodulated 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 multilayered 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 illposed. 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 illposedness are found in [10, 11, 12, 13, 14, 15].
Finally, in [22] we have a regularized, numerical algorithm to reconstruct timedynamic 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 . The ”full” measurement is , where . From we may manually choose a more interesting neighbourhood where , containing a visible hot spot we want to analyze in more detail. The FFT of these measurements are and , for each spatial point , more specifically we write
(18)  
(19) 
where the amplitudes and phases for each describe the signal component corresponding to the frequency , where is the sampling rate and is the number of samples. The index is meant to go through the interesting frequencies, if not all available from FFT, in which case we would have roughly because of ShannonNyquist sampling theorem.
For each point and frequency we can write
(20)  
(21) 
and subsequently we can express the diffused source component of (13) as
(22) 
where
(23)  
(24) 
Also the static solution (12) can then be written as
(25) 
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 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 from which we get the model complient parameters for that particular source location. Then, the penalty functional is
(26) 
In the dynamic case, for each candidate point source we compute its nearest boundary point on and then the values for each frequency component . Using the dynamic model (25) these are then turned to reconstructions . However, we both cut away impossible spectral content and find the best fit from the source locations. By defining truncation parameters and , we define the set as the set for which has spectral amplitude over (cutting off too small data) and which has the size of at least (cutting off too little pieces of measurement data). The lower and higher is, the more robust we need elsewhere in the algorithm. For each source candidate, for the amplitude data, we use the penalty functional
(27) 
where the profile is given by the source candidate , the computed value and the dynamic model (25). The total amplitude penalty functional is
(28) 
For the phase penalty functional, we use
(29) 
where is the model computed profile (at each source location ) normalized to have the same scale as the measurement on . Similarly to the amplitude penalty we write
(30) 
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 , satisfying for each , where is the amplitude of the signal with harmonic content .
4.1. General algorithm descriptions
There are many possibilities how to arrange the algorithm to find the location and the spectral content from the measurement, depending also on what we expect the situation to be. We can think of the following examples:

Given , reconstruct .

Given , reconstruct a collection of point sources with signals .

Given , reconstruct .

Given , reconstruct and .

Given and , reconstruct and .
The algorithm (A1) will be difficult to realize and is not in the scope of this paper. Our model only considers pointlike 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 pulselike 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 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 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 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 , 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 . To do this we use the following subalgorithm (AQ) for each candidate location :

Sort the values of or and take an average of the first and last values to get the minimum temperature and maximum temperature of the measurement patch. Besides the sorted vector of measurements, use their locations to compute the corresponding average locations and and their distances , from .

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

Compute the values again from the equation (25) using and , and finally compute
In this way we can compute for any candidate the power and the parameter . If or is very noisy the above will perform less well, but the integer 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 , 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 or has been processed by FFT to get the collection or of amplitudes and phases. Define an amplitude limit below which we discard amplitude data. Then, given any source location candidate we use the following subalgorithm (AA):

Compute the boundary point closest to the source and write for the amplitude values on a neighbourhood of , for each frequency bin . Write .

For each , use (22) and to compute the model predicted amplitude profile . Compute the submeasurement patch .

If or a very small set, put as the reconstruction. Otherwise use as the domain where the phase error is computed.
The phase profile needs special processing as its raw values are always inside . We use Matlab’s unwrapfunction to first make the phase profile continuous, then put the phase value on to zero by subtracting the value computed earlier in conjunction with the amplitudes . 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 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 to compute the phase error, two problems arise. First, the smaller the area of , the smaller the magnitude of the phase error. Second, the phase information outside of 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 . 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:

From , pick an interesting neighbourhood manually.

Take a collection of source candidates and for each of them use (AQ) to compute the parameters and . For each triplet compute from (26).

The reconstruction is the triplet having the smallest error .
4.5. Dynamic reconstruction algorithm
The algorithm (3b) in more detail is the following:

From , pick an interesting neighbourhood manually.

Compute FFT of to get the amplitude and phase data.

Take a grid of source candidates and for each use the subalgorithm (AAn) to get the modelconforming amplitudes and usable phase areas .
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 ballshaped 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 timedynamic 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
.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 near the boundary, nonnoisy 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.
static, R=10  dynamic, f=0.2Hz  

center  
A [40 40 50]  
B [35 35 45]  
C [30 30 40] 
6. Numerical tests
We used truncation values and pixels. In the static case we took different source candidate locations, 40 in each spherical dimension to cover the whole domain near the measurement patch . In the dynamic case we used 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 nonnoisy 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.
location error  

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 
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 nonnoisy 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 was efficiently chosen so that all other frequency components other than the third FFT bin corresponding to exactly were zero, in all reconstructions.
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 
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 
7. Conclusions
In this article we have presented a new approach to detecting static and timevarying pointlike 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 semirealistic boundary condition used in real bioheat models; still the analytical pointsource 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 illposedness 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.
References
 [1] W. Ben Lardi, C. IbarraCastanedo, M. Klein, A. Bendada, X. Maldague, Experimental Comparison of Lockin and Pulsed Thermography for the Nondestructive Evaluation of Aerospace Materials, Electrical and Computing Engineering Department, Université Laval, Quebec City (2010).
 [2] S. Marinetti, Y. A. Plotnikov, W. P. Winfree, A. Braggiotti, Pulse phase thermography for defect detection and visualization, Proc. SPIE 3586, Nondestructive Evaluation of Aging Aircraft, Airports, and Aerospace Hardware III (1999).

[3]
N. Dalir,
Exact Analytical Solution for 3D TimeDependent Heat Conduction in a Multilayer Sphere with Heat Sources Using Eigenfunction Expansion Method
, Hindawi Publishing Corporation, International Scholarly Research Notices (2014).  [4] T. F. Flint, J. A. Francis, M. C. Smith, A. N. Vasileiou, Semianalytical solutions for the transient temperature fields induced by a moving heat source in an orthogonal domain, International Journal of Thermal Sciences 123 (2018), 140 – 150.
 [5] F. Davoodi, A. Abbas Nejad, A. Shahrezaee, M. J. Maghrebi, Control parameter estimation in a semilinear parabolic inverse problem using a high accurate method, Applied Mathematics and Computation 218 (2011), 1798 – 1804.

[6]
W. Liao, M. Dehghan, A. Mohebbi,
Direct numerical method for an inverse problem of a parabolic partial differential equation
, Journal of Computational and Applied Mathematics 232 (2009), 351 – 360.  [7] J. Liu, B. Wang, Z. Liu, Determination of a source term in a heat equation, International Journal of Computer Mathematics 87(5) (2010), 969 – 975.
 [8] J. R. Cannon, Determination of an unknown heat source from overspecified boundary data, SIAM Journal of Numerical Analysis 5(2) (1968), 275 – 286.
 [9] J. R. Cannon, D. Zachmann, Parameter determination in parabolic partial differential equations from overspecified boundary data, Int. J. Engng Sci. 20(6) (1982), 779 – 788.
 [10] W. Cheng, LL. Zhao, CL. Fu, Source term identification for an axisymmetric inverse heat conduction problem, Computers and Mathematics with Applications 59 (2010), 142 – 148.
 [11] W. Cheng, YJ. Ma, CL. Fu, Identifying an unknown source term in radial heat conduction, Inverse Problems in Science and Engineering 20(3) (2012), 335 – 349.
 [12] FF. Dou, CL. Fu, Determining an unknown source in the heat equation by a wavelet dual least squares method, Applied Mathematics Letters 22 (2009), 661 – 667.
 [13] FF. Dou, WaveletGalerkin Method for Identifying an Unknown Source Term in a Heat Equation, Mathematical Problems in Engineering (2012).
 [14] D. N. Hao, B. V. Huong, N. T. N. Oanh, P. X. Thanh, Determination of a term in the righthand side of parabolic equations, Journal of Computational and Applied Mathematics 309(2017), 28 – 43.
 [15] P. Wang, K. Zheng, Reconstruction of spatial heat sources in heat conduction problems, Applicable Analysis 85(5)(2006), 459 – 465.
 [16] J. R. Cannon, R. E. Ewing, Determination of a Source Term in a Linear Parabolic Partial Differential Equation, Journal of Applied Mathematics and Physics 27 (1976), 393 – 401.
 [17] J. R. Cannon, S. PerezEsteva, An inverse problem for the heat equation, Inverse Problems 2 (1986), 395 – 403.
 [18] J. R. Cannon, S. PerezEsteva, Uniqueness and stability of 3D heat sources, Inverse Problems 7 (1991), 57 – 62.
 [19] G. Talenti, S. Vessella, A note on an illposed problem for the heat equation, J. Austral. Math. Soc. (Series A) 32 (1982), 358 – 368.
 [20] D. D. Trong, T. T. Tuyen, P. T. Nam, A. P. N Dinh, Determine the special term of a twodimensional heat source, Applicable Analysis 88(3) (2010), 457 – 474.
 [21] M. Yamamoto, Conditional Stability in Determination of Force Terms of Heat Equations in a Rectangle, Mathl. Comput. Modelling 18(1) (1993), 79 – 88.
 [22] N. Renault, S. Andre, D. Maillet, C. Cunat, A twostep regularized inverse solution for 2D heat source reconstruction, International Journal of Thermal Sciences 47 (2008), 834–847.
 [23] N. Renault, S. Andre, D. Maillet, C. Cunat, A spectral method for the estimation of a thermomechanical heat source from infrared temperature measurements, International Journal of Thermal Sciences 49 (2010), 1394–1406.
 [24] Y. X. Zhang, C. L. Fu, Y. J. Ma, An a posteriori parameter choice rule for the truncation regularization method for solving backward parabolic problems, Journal of Computational and Applied Mathematics 255 (2014), 150 – 160.
 [25] M. Anbar, L. Milescu, A. Naumov, C. Brown, T. Button, C. Carty, K. AlDulaimi, Detection of Cancerous Breasts by Dynamic Area Telethermometry, IEEE Engineering in Medicine and Biology 20(5) (2001), 80–91.
 [26] T. M. Button, H. Li, P. Fisher, R. Rosenblatt, K. Dulaimy, S. Li, B. O’Hea, M. Salvitti, V. Geronimo, C. Geronimo, S. Jambawalikar, P. Carvelli, R. Weiss, Dynamic infrared imaging for the detection of malignancy, Phys. Med. Biol. 49 (2004), 3105–3116.
 [27] R. Joro, A.L. Lääperi, P. Dastidar, S. Soimakallio, T. Kuukasjärvi, T. Toivonen, R. Saaristo, R. Järvenpää, Imaging of breast cancer with mid and longwave infrared camera, Journal of Medical Engineering & Technology 32(3) (2008), 189–197.
 [28] R. Joro, A.L. Lääperi, S. Soimakallio, R. Järvenpää, T. Kuukasjärvi, T. Toivonen, R. Saaristo, P. Dastidar, Dynamic infrared imaging in identification of breast cancer tissue with combined image processing and frequency analysis, Journal of Medical Engineering & Technology 32(4) (2008), 325–334.
 [29] R. Joro, A.L. Lääperi, P. Dastidar, R. Järvenpää, T. Kuukasjärvi, T. Toivonen, R. Saaristo, S. Soimakallio, A Dynamic Infrared ImagingBased Diagnostic Process for Breast Cancer, Acta Radiol 8 (2009), 860–869.
 [30] R. Joro, P. Dastidar, V. Iivonen, H. Ylänen, S. Soimakallio, NADINE: new approaches to detecting breast cancer by sequential wavelength imaging with the aid of novel frequency analysis techniques, Journal of Medical Engineering & Technology 36(5) (2012), 251–260.
 [31] H. Ahmadikia, R. Fazlali, A. Moradi, Analytical solution of the parabolic and hyperbolic heat transfer equations with constant and transient heat flux conditions on skin tissue, International Communications in Heat and Mass Transfer 39 (2012), 121 – 130.
 [32] A. Lakhssassi, E. Kengne, H. Semmaoui, Modifed pennes’ equation modelling bioheat transfer in living tissues: analytical and numerical analysis , Natural Science 2 (2010), 1375 – 1385.
 [33] K. Das, R. Singh, S. C. Mishra, Numerical analysis for determination of the presence of a tumor and estimation of its size and location in a tissue, Journal of Thermal Biology 38 (2013), 32 – 40.
 [34] K. Das, S. C. Mishra, Study of thermal behavior of a biological tissue: An equivalence of Pennes bioheat equation and Wulff continuum model, Journal of Thermal Biology 45 (2014), 103 – 109.
 [35] D. Arthur, Towards Application of Thermal Infrared Imaging in Medical Diagnosis: Protocols and Investigations, PhD thesis, Curtin University (2014).
 [36] G. Shi, F. Han, L. Wang, C. Liang and K. Li, Qr curve of thermal tomography and its clinical application on breast tumor diagnosis, Optical Society of America 6(4) (2015), 1109–1123.
 [37] D. Gurarie, M445: Heat equation with sources, Lecture note on ”MATH 445. Introduction to Partial Differential Equations”, Case Western Reserve University.
 [38] H. S. Carslaw, J. C. Jaeger, Conduction of Heat in Solids, Oxford University Press, second edition (1959).
Comments
There are no comments yet.