# Non-local Operational Anisotropic Diffusion Filter

Fábio A. M. Cappabianco  Group for Inovation Based on Images and Signals
(GIBIS)
Department of Science and Technology
Federal University of São Paulo
São José dos Campos/SP, Brazil
<EMAIL_ADDRESS>Petrus P. C. E. da Silva  Department of Medical
Phisics
A.C. Camargo Cancer Center
São Paulo/SP, Brazil
<EMAIL_ADDRESS>
###### Abstract

High-frequency noise is present in several modalities of medical images. It
originates from the acquisition process and may be related to the scanner
configurations, the scanned body, or to other external factors. This way,
prospective filters are an important tool to improve the image quality. In
this paper, we propose a non-local weighted operational anisotropic diffusion
filter and evaluate its effect on magnetic resonance images and on kV/CBCT
radiotherapy images. We also provide a detailed analysis of non-local
parameter settings. Results show that the new filter enhances previous local
implementations and has potential application in radiotherapy treatments.

###### Index Terms:

anisotropic diffusion, image filtering, high-frequency noise, radiotherapy, kV
images, cone beam computed tomography

##  I Introduction

Medical imaging is subject to several kinds of noise and artifacts coming from
the acquisition and storage procedure [1]. The nature and appearance of the
noise depend on distinct factors such as the imaging modality, the scanned
body part, scanner parameters and positioning, among others [2, 3].

Despite recent improvements to the acquisition process, high-frequency filters
are essential to improve medical image quality. Nevertheless, their
application is not a trivial task, since some of the original signals may be
removed and new artifacts might be inserted to the image while removing the
noise [4, 5].

There are three generations of high-frequency noise filters: the isotropic
which were proposed since the beginning of digital medical imaging [6, 7]; the
anisotropic – which became popular in the 90s [8] – are still in use for its
robustness to generating artifacts, and is present in some popular medical
imaging tools such as FLS [9] and 3D Slicer [10]; and the non-local
anisotropic filters which started to be developed around 2005 with the Non-
Local Means (NLM) algorithm [11]. Even though the third generation produces
better results in terms of several metrics, in the context of medical imaging,
one should be even more careful since they tend to insert more artifacts to
the filtered image [4].

In this paper we focus on the implementation and testing of an operational
third generation filter which is an extension of the anisotropic diffusion
filter (ADF) [12] for reducing addictive noise present in MRI and CT images.
The proposed implementation is operational with respect to local and non-local
conservativeness parameter. We also conduct an analysis of the non-local
parameters in the models proposed in [13]. Our experiments include synthetic
MRI images to estimate parameter values and also a qualitative evaluation of
kV/CBCT radiotherapy images by specialists.

The reminder of this paper is organized as follows: Section II presents a
review of related works in high-frequency image segmentation; Section III
contains an explanation of ADF and its extensions; Section IV proposes our new
filtering method; Section V shows our quantitative and qualitative
experiments; and Section VI states our conclusions.

##  II Related Works

The first generation of high-frequency noise filters are isotropic, that is,
the same operation is applied to each image pixel despite its contents or
position. Some of the traditional isotropic filters are the arithmetic and
geometrical means, median, and alpha-trimmed mean filter [7].

The second generation of high-frequency noise filters is applied over other
transform domains and/or with local anisotropic functions. These filters apply
distinct operations according to pixel features such as intensity frequency or
local patterns. Wavelet transform and cosine transform based [14, 15],
bilateral and trilateral [16, 17], and the diffusion [18, 19] are among the
most popular filters. While first generation filters were simpleton
implementations, the second generation filters bear several non-trivial
parameters which yielded lots of tuning and evaluation papers [2].

Third generation filters are in most part extensions of the existing second-
generation implemented to take advantage of non-local data. The basic idea is
to use similar patches in order to estimate and remove high-frequency noise.
Examples of bilateral filter extensions are the NLM [11], BM3D [20]. Third
generation filters have also been applied to medical imaging in the context of
MRI and CT [3].

In the literature, there is not much effort in developing a non-local ADF
filter for MRI and CT images [2, 3]. Only two third generation of ADF were
proposed in [21, 22]. These extensions are just a trivial addition of non-
local patches to the algorithm with no mathematical or empirical selection of
optimum parameters.

In [21], the authors use the original edge-stopping functions. In [22], the
Tukey’s biweight robust estimator was used. Either of them estimated similar
patches based on pre-filtered images. The second implementation computed
adjacent patches at each iteration of the filter which is unnecessary and
computationally expensive.

Therefore, there is a lack of studies of a non-local ADF which is the main
subject of this paper.

##  III Anisotropic Diffusion Filter

The ADF was proposed by Perona and Malik [12]. It consists in applying
concepts of fluid thermodynamics to filter images. Just as molecules with
different temperature from its adjacency tends to exchange heat until
achieving homogeneous temperature, there would have a propensity to match
high-frequency noisy pixel intensities with their local neighbours. In this
scenario, image edges are similar to tow large adjacent volumes of fluid with
higher capacitance, demanding more time to homogenize the temperature in
comparison to the noisy smaller volumes. Therefore, ADF is capable of
eliminating the high-frequency noise while preserving image strong edges.

The discrete implementation of ADF is given by Equation 1:

$I_{s}^{t+1}\approx
I_{s}^{t}+\frac{\lambda}{|\eta_{s}|}\sum_{p\in\eta_{s}}{g(\nabla
I_{s,p}^{t},\gamma^{t})\nabla I_{s,p}^{t}}$ (1)

$I^{t}_{s}$ stands for the intensity of pixel s in image $I$ at the instant
$t$, $\lambda$ is a diffusion rate scalar, $\gamma^{t}$ is a positive variable
related to the diffusion or smoothing strength which decreases at each
iteration, $\eta_{s}$ represents the set of pixels adjacent to $s$, $g(\cdot)$
is an edge stopping function (ESF) which controls the diffusion process, and
$\nabla I^{t}_{s,p}$ is the image directional gradient from $s$ to $p$ at
instant $t$. The directional gradient $\nabla I^{t}_{s,p}$ can be approximated
by $I^{t}_{p}-I^{t}_{s}$.

The diffusion process follows the ESF which has been widely studied and tested
[19, 23]. Tukey’s biweight robust estimator ESF is given in Equation 2 and is
an excellent option since it does not affect edges with intensity variation
over a certain threshold [13].

$g(\nabla
I_{s,p}^{t},\gamma^{t})=\left\\{\begin{array}[]{l}\left[1-\left(\frac{\nabla
I_{s,p}^{t}}{\sqrt{5}\gamma^{t}}\right)^{2}\right]^{2},\nabla
I_{s,p}^{t}\leq\gamma^{t}\sqrt{5}\\\
0,\text{otherwise.}\end{array}\color[rgb]{1,1,1}\right\\}$ (2)

With respect to parameter $\lambda$, the authors in [18] proposed the maximum
values in order to keep a monotonic variation of the intensities at each
iteration, so that no artifacts are inserted to the image. In [13], the
authors corrected one of the constants for 26-adjacency. Table I shows the
expected maximum values of $\lambda$ and of the fraction $\lambda/|\eta_{s}|$
according to the number of dimensions of the image ($D$) and to the adjacency
size ($|\eta_{s}|$).

TABLE I: Maximum ADF values for constant $\lambda$ and expression $\lambda/|\eta_{s}|$, given the input image dimensions $D$ and the adjacency size $|\eta_{s}|$. $D$ | $|\eta_{s}|$ | $\lambda$ | $\lambda/|\eta_{s}|$
---|---|---|---
2 | 4 | 4/5 | 1/5
2 | 8 | 8/7 | 1/7
3 | 6 | 6/7 | 1/7
3 | 18 | 18/13 | 1/13
3 | 26 | 78/47 | 3/47

There were several different solutions proposed to estimate the initial
$\gamma^{0}$ for the first iteration and the number of iterations $T$. Some
approaches used global and/or local gradient to compute $\gamma$ [12, 19, 24]
while others use a planar region (i.e. without strong edges) [25].

The approach in [13] used an adaptive method based on both the global gradient
and a planar region in order to estimate the best initial $\gamma^{0}$. While
the gradient indicates an upper bound $\gamma{E}$ to avoid filtering edges,
the planar region provides a lower bound $\gamma{F}$ so that practically all
noise is removed. In case of low signal to noise ratio (SNR)111In the
experiments, Gaussian or Rician noise of around 3% of amplitude is the maximum
SNR which is perfectly filtered without blurring significant edges., it may
happen that $\gamma{F}>\gamma{E}$. In these cases, a value between $\gamma{E}$
and $\gamma{F}$ should be chosen. Lower $\gamma^{0}$ translates into a
conservative ADF which preserves the stronger edges, while higher $\gamma^{0}$
constitutes into an aggressive filtering, removing all noise.

Also, since its proposal [12] and as corroborated by [24], it is very
important to update $\gamma^{t}$ at each iteration. As showed in [13], it is
fundamental to reduce $\gamma^{t}$ at each iteration in order to preserve the
borders and it may be very expensive to recompute $\gamma^{t}$ at each
iteration based on the noise in the planar region and on the strong edges. A
conservative proposed solution was to set $\gamma^{t}=0.25\gamma^{t-1}$. That
is because it is expected that the strongest noise has at most one-fourth of
its adjacent pixels with a similar intensity, while for the edges, it is
reasonable to expect that at least one-fourth of the pixels have distinct
intensity. Therefore, reducing $\gamma$ by $0.16\gamma$ was a conservative
adopted criteria222For $\gamma$ reduction values of other ESFs, please refer
to [13]..

Finally, with respect to the number of iterations $T$, in [25] the authors
verify the intensity of the strongest edges. That solution is not suitable for
decreasing $\gamma^{t}$ with Tukey’s biweight robust estimator given in
Equation 2. In [13], the authors propose to stop as $\gamma^{t}<\gamma_{\cal
E}/7$ which does not depend on the noise intensity itself.

##  IV Non-Local Optimal Anisotropic Diffusion Filter

The non-local ADF (NL-ADF) proposed in [21, 22] do not optimize any filter
parameter. We propose a novel non-local ADF with enhanced parameter
optimization. It has improvements over the previous optimization procedure
proposed in [13]. We also study the non-local patches influence in the
filtering process.

###  IV-A Local Parameter Optimization

Instead of estimating the initial $\gamma^{0}$ based on the standard deviation
of edgy or noisy regions, in this paper we propose a more robust and
straightforward strategy. Noise pixels in flat regions have a distinct
intensity value than all their adjacents. An edge pixel, on the other hand,
has a similar intensity to pixels in the same side of the edge and distinct
intensity from pixels in the opposite side of the edge. Therefore, instead of
computing the standard deviation of pixels to estimate noise, we propose to
use the sum of the directional gradients $G_{s}$ as in Equation 3:

$G_{s}=\sum_{p\in\eta_{s}}{I_{s}-I_{p}}$ (3)

Note that the direction of $G_{s}$ is always computed toward the source pixel
$s$. This way, edge pixels will have opposite edge pixel intensities
neutralized, which noise pixels will display the highest values of $G_{s}$. We
then select standard deviation ($\sigma_{G_{s}}$) of the 5$\%$ highest $G_{s}$
pixels as our $\gamma^{0}$, as the maximum effect of of the ADF with Tukey’s
biweight robust estimator happens for $\gamma^{0}=\sigma_{G_{s}}$.

A second contribution is to revise the $\gamma^{t}$ reduction strategy
proposed in [13] based on 1/4 of the adjacent pixels as explained in Section
III as it was a very conservative choice. In this paper, we evaluated faster
and slower $\gamma^{t}$ reductions based on values in range
$[0.16\gamma^{t-1},0.96\gamma^{t-1}]$.

With respect to the stopping criteria a simple alternative is to stop
filtering the image as $\gamma^{t}\leq 0.01I^{M}$ ($I^{M}$ is the maximum
intensity of an image $I$), as 1$\%$ of maximum intensity is imperceptible by
humans eyes in medical images.

Finally, diffusion of adjacent pixels should also consider the distance
between the source and target pixels. Therefore, we propose the weighted ADF
(WADF) based on Equation 4:

$I_{s}^{t+1}\approx
I_{s}^{t}+\frac{1}{|\eta_{s}|}\sum_{p\in\eta_{s}}{\frac{g(\nabla
I_{s,p}^{t},\gamma^{t})}{d_{s,t}}\nabla I_{s,p}^{t}}$ (4)

where $d(s,t)$ is the distance between pixels $s$ and $t$. Note that,
according to the study proposed in [18, 13] with the WADF for a monotonic
variation of the intensities at each iteration we may set $\lambda=1$,
eliminating this parameter from the original ADF Equation 1.

###  IV-B Non-Local Anisotropic Diffusion Function

An interesting question overlooked in previous publications is the weight or
distance of non-local patches. Most of the papers, including the previous NL-
ADF [21, 22], position the non-local patches in a discrete way in an
additional third or fourth image dimension. That is, the closest patch is one-
pixel-wide distant from the patch of interest. In a similar way, the original
NLM [11] uses the spatial distance as a weighting factor, but this is an
arbitrary solution. BM3D [20] and BM4D [26] also build a structure with an
additional discrete dimension.

As the ADF does not depend on the distance of patches itself, we decided to
evaluate patch distance and quantity over the filtering result. The $n$ more
similar patches will have the same distance $d$ to the filtered pixel $s$ as a
third (or fourth) dimension to a 2D (or 3D) image. Note that $d$ may assume
any float value. We will evaluate multiple patch distance variation in a
future work.

The choice of the distance $d$ also impacts on the computation of constant
$\lambda$. Still, considering the WADF, we may set $\lambda=1$ even for the
non-local weighted anisotropic diffusion filter (NL-WADF) in Equation 5:

$I_{s}^{t+1}\approx
I_{s}^{t}+\frac{1}{|H_{s}|}\sum_{p\in\eta_{s}}{\frac{g(\nabla
I_{s,p}^{t},\gamma^{t})}{d_{s,t}}\nabla I_{s,p}^{t}}$ (5)

$H_{s}$ comprehends all adjacent pixels to $s$ both local and non-local.

In the case of NL-WADF, the number of pixels with similar intensity to noise
and edge pixels should stay proportionally the same. Therefore, we propose
also to evaluate the $\gamma^{t}$ decrease in every iteration using the same
range of $[0.32\gamma^{t-1},0.96\gamma^{t-1}]$.

Also, we determine the similar patches based on filtered images by a median
filter [27]. After this, the patches are used by NL-WADF over the original
image.

##  V Experiments

The first experiment involves estimating the best parameter values for local
and non-local WADF. For this purpose, we used BrainWeb Phantom dataset [28]
changing the noise intensity between 1 and 9$\%$ of the image maximum
intensity. Figure 1 contains slices of the ground-truth image and of an image
corrupted by a noise of 7$\%$ of the maximum intensity. The second experiment
is a qualitative evaluation in the context of radiotherapy imaging as shown in
Figure 2. Images were acquired with Varian TrueBeam and Varian Clinac IX
Radiotherapy Systems.

---
(a)
(b)
Figure 1: Sample slices of BrainWeb Phantom dataset image (a) Ground-truth
with no noise; (b) Image corrupted with gaussian noise with 7$\%$ of maximum
image intensity.
---
(a)
(b)
Figure 2: Sample slices of radiotherapy images (a) 2D kV planar image; (b) 3D
CBCT image.

###  V-A Quantitative Evaluation

We used the following metrics to evaluate the quality of the filtered image in
comparison with the ground-truth image with no noise: structural similarity
index(SSIM) Equation 6, mean square error(MSE) Equation 7, peak signal-to-
noise ratio(PSNR) Equation 8, and image quality index(IQI) Equation 9. We will
not present here results for SSIM, since all images achieved accuracy equal to
$1.0$. Lower scores are better in terms of MSE and higher scores are better
with respect to SSIM and PSNR.

$SSIM(I,J)=\frac{(2\mu_{I}\mu_{J}+c_{1})(2\sigma_{IJ}+c_{2})}{(\mu_{I}^{2}+\mu_{J}^{2}+c_{1})(\sigma_{I}^{2}+\sigma_{J}^{2}+c_{2})}$
(6) $c_{1}=0.0001|I||J|$ $c_{2}=0.0009|I||J|$

$I$ and $J$ are the filtered and ground-truth images, respectively, $\mu_{I}$
and $\sigma_{I}$ are the mean and standard deviation of pixel intensity in
$I$, respectively, and $|I|$ is number of pixels in $I$. $\sigma_{I}J$ is the
covariance between images $I$ and $J$.

$MSE(I,J)=\frac{1}{|I|}\sum_{s}(I_{s}-J_{s})^{2}$ (7)
$PSNR(I,J)=20*\log(\max(I^{M},J^{M}))-10\log(MSE(I,J))$ (8)

$\max(x,y)$ is the maximum value between x and y. $I^{M}$ is the maximum pixel
intensity in $I$.

$IQI(I,J)=\frac{(4\mu_{I}\mu_{J}\sigma_{IJ})}{(\mu_{I}^{2}+\mu_{J}^{2})(\sigma_{I}^{2}+\sigma_{J}^{2})}$
(9)

For the WADF, we tested the $\gamma^{t}$ reduction rate of $\gamma
R=\\{0.16\gamma^{t-1},0.32\gamma^{t-1},0.48\gamma^{t-1},0.64\gamma^{t-1},0.80\gamma^{t-1},0.96\gamma^{t-1}\\}$.
We also evaluated a conservativeness parameter which sets the initial
$\gamma_{0}=0.2\sigma_{G_{s}},0.4\sigma_{G_{s}},0.6\sigma_{G_{s}},0.8\sigma_{G_{s}},and\sigma_{G_{s}}$.

We notice in Table II that the best parameter set for WADF with respect to
$\gamma^{t}$ reduction rate in average to all noise levels is between
0.64$\gamma^{t-1}$ and 0.80$\gamma^{t-1}$ which is much higher than the value
adopted in [13]. It is important to mention that the accuracy given by the
proposed metrics are not always favorable to 0.64$\gamma^{t-1}$ and
0.80$\gamma^{t-1}$ reduction, but 0.32$\gamma^{t-1}$ provides better results
for more conservative filtering(i.e. $\gamma_{0}\leq 0.4\sigma_{G_{s}}$).

TABLE II: Quantitative results of application of WADF over BrainWeb Phantom MRI dataset. $\gamma R$ | IQI | MSE | PSNR
---|---|---|---
16$\gamma^{t-1}$ | 0.9898 | 3,742 | 27.59
32$\gamma^{t-1}$ | 0.9945 | 2,986 | 28.77
48$\gamma^{t-1}$ | 0.9964 | 2,836 | 29.09
64$\gamma^{t-1}$ | 0.9970 | 2,788 | 29.21
80$\gamma^{t-1}$ | 0.9949 | 2,777 | 29.26
96$\gamma^{t-1}$ | 0.9942 | 2,802 | 29.22

For the NL-WADF, we tested the following parameters: $\gamma
R=\\{0.32\gamma^{t-1},0.48\gamma^{t-1},0.64\gamma^{t-1},0.80\gamma^{t-1},0.96\gamma^{t-1}\\}$;
$\gamma_{0}=0.2\sigma_{G_{s}},0.4\sigma_{G_{s}},0.6\sigma_{G_{s}},0.8\sigma_{G_{s}},and\sigma_{G_{s}}$;
search radius (SR = 2.0, 3.0 4.0 pixels), non-local patch radius (PR = 1.1,
1.5, 1.8 pixels), the distance from the central pixel of non-local patches (PD
= 0.5, 1.0, 2.0 pixels), and the number of non-local patches ($\\#$P = 1, 2).

Figure 3 shows sample slices of image showed in Figure 1(b) filtered
conservatively($\gamma^{0}=0.4\sigma_{G_{s}}$) and
aggressively($\gamma^{0}=0.8\sigma_{G_{s}}$) with parameters SR=3.0, PR=1.9,
PD=1.0, $\\#$P=2, and $\gamma$R=0.80$\gamma^{t-1}$.

---
(a)
(b)
Figure 3: Result of filtering image in Figure 1(b) with (a) conservative NL-
WADF and (b) aggressive NL-WADF.

Table III contains the best results for NL-WADF in terms of the proposed
metrics. The best set of parameters in average according to MSE and PSNR is to
use SR=4.0, PR=1.1, PD=0.5, $\\#$P=2, and a $\gamma$R= 0.32$\gamma^{t-1}$.
Note that $\gamma_{0}$ is not related to the quality of the result, but to the
desired conservativeness, which is more related to the goal application.
According to IQI, the best set of parameters is SR=3.0, PR=1.9, PD=0.5,
$\\#$P=1, and $\gamma$R=0.64$\%$.

TABLE III: Quantitative results of application of NL-WADF over BrainWeb Phantom MRI dataset. SR; PR; PD; #R; $\gamma R$ | IQI | MSE | PSNR
---|---|---|---
2.0; 1.1; 0.5; 2; 0.32 | 0.9968 | 2,719 | 29.46
3.0; 1.1; 0.5; 2; 0.32 | 0.9967 | 2,712 | 29.48
4.0; 1.1; 0.5; 2; 0.32 | 0.9967 | 2,710 | 29.49
4.0; 1.9; 1.0; 2; 0.32 | 0.9941 | 3,022 | 28.71
4.0; 1.1; 0.5; 2; 0.48 | 0.9956 | 2,749 | 29.45
2.0; 1.9; 2.0; 1; 0.64 | 0.9972 | 2,793 | 29.20
3.0; 1.9; 0.5; 1; 0.64 | 0.9972 | 2,776 | 29.23
4.0; 1.9; 1.0; 2; 0.64 | 0.9972 | 2,788 | 29.20

Note though, that this is by no means the best parameter set for all noise
levels and desired conservativeness. Based on the results in Tables II and III
we can also notice that NL-WADF may generate better results as compared to
WADF, but this improvement is not statistically significant.

Other important observations are: using one non-local patch produces in
65.7$\%$ of cases a better accuracy than employing two of them and in 69.9$\%$
of the times for more aggressive filtering (i.e. $\gamma_{0}\geq 0.6$). In
59.1$\%$ of the times, for a noise level of 7$\%$, using two non-local patches
produces a better accuracy. That is the only noise level with better results
using 2 non-local patches.

In terms of IQI, $\gamma^{t}$ reduction of 0.96$\gamma^{t-1}$ produces the
higher accuracy in 99$\%$ of the tests, except while dealing with a noise
level of 7$\%$.

The analysis of the size of the non-local patch radius is very similar to the
number of patches. A radius of 1.1 pixels generates higher accuracy in
52.3$\%$ of times. Smaller radius are a better choice in 62.2$\%$ of times for
aggressive filters (i.e $\gamma_{0}=0.6\sigma_{G_{s}}$).

Setting non-local patches at 0.5 pixels from the source patch generates the
best results in 55$\%$ of the times and in 61.8$\%$ of the times for more
aggressive filtering (i.e. $\gamma_{0}\geq 0.6$). This means that non-local
patches should have a stronger impact in removing the noise.

With respect to the size of the similar patches search region, a radius of 2.0
pixels is usually sufficient and produces the best results. Conservative
parameters are an exception (i.e. $\gamma_{0}\leq 0.4$) over high-intensity
noise in which a radius of 4.0 pixels produces better results. This indicates
that a larger search radius is an important parameter to conserve image edges.

###  V-B Qualitative Evaluation

Using the best parameter set for NL-WADF, we performed a qualitative test with
cone beam computed tomography (CBCT) and kV planar radiotherapy images of the
low pelvis region provided by hospital AC Camargo in São Paulo/SP, Brazil. A
physicist specialist evaluated the quality of the filtered images and their
utility for image visualization purposes.

Figures 4 and 5 shows sample slices of filtered planar kV and 3D CBCT images
with conservative and aggressive NL-WADF.

---
(a)
(b)
Figure 4: Resulting image after filtering image in Figure 2(a) with (a)
conservative NL-WADF and (b) aggressive NL-WADF.
---
(a)
(b)
Figure 5: Resulting image after filtering image in Figure 2(b) with (a)
conservative NL-WADF and (b) aggressive NL-WADF.

Five filtering levels for conservative to aggressive of ten images (3 planar
and 7 3D CBCT volumes) were evaluated. The specialist assigned a score from 1
(poor) to 5 (excellent) to each of the filtered images. Table IV show the
scores. As noticed for radiotherapy visualization purposes, the best filter
consists of applying $\gamma^{0}=0.4\sigma_{G_{s}}$. More conservative
parameters do not remove significantly the noise. Using
$\gamma^{0}>0.4\sigma_{G_{s}}$ on the other hand on more aggressive filter
removes important details in the opinion of the specialist.

TABLE IV: Qualitative evaluation of planar kV and 3D CBCT images by physicist specialist. Grades from 1 (poor) to 5 (excellent) were given to each image. Mean values are presented here. $\gamma^{0}$ | Score
---|---
0.2$\sigma_{G_{s}}$ | 4.9
0.4$\sigma_{G_{s}}$ | 5.0
0.6$\sigma_{G_{s}}$ | 4.9
0.8$\sigma_{G_{s}}$ | 4.9
$\sigma_{G_{s}}$ | 4.8

##  VI Conclusions

In this paper, we propose an operational non-local anisotropic diffusion
filter for additive noises present in medical images such as MRI and CT. The
main contributions of the filter were the distance factor never considered
before for non-local patches filtering for ADFs and the weighted contribution
of adjacent pixels. Also, we present an extensive evaluation of the parameter
values in order to define the best parameter set in terms of SSIM, PSNR, MSE,
and DE. Results show that the new non-local version of the ADF is
significantly better than the local filter. Finally, we also present a
qualitative evaluation of the filter in the context of visual quality of
radiotherapy imaging which showed that filtering the image is an important
procedure, but the procedure must be conservative avoiding removing important
anatomical details.

Future works involve proposing more sophisticated strategies for setting the
distance between non-local adjacent patches depending on the source pixel
neighbourhood and also based on the neighbourhood similarity.

## Acknowledgment

We would like to thank AC Camargo and Dr Luiz Juliano Neto for the support
with the radiotherapy images. We also thank the training provided by Varian
Medical Systems in agreement with the Brazilian Ministry of Health.

## References

  * [1] J. Beutel, H. L. Kundel, and R. L. Van Metter, _Handbook of medical imaging_. Spie Press, 2000, vol. 1.
  * [2] J. Mohan, V. Krishnaveni, and Y. Guo, “A survey on the magnetic resonance image denoising methods,” _Biomedical signal processing and control_ , vol. 9, pp. 56–69, 2014.
  * [3] M. Diwakar and M. Kumar, “A review on ct image noise and its denoising,” _Biomedical Signal Processing and Control_ , vol. 42, pp. 73–88, 2018.
  * [4] M. Maggioni, V. Katkovnik, K. Egiazarian, and A. Foi, “Nonlocal transform-domain filter for volumetric data denoising and reconstruction,” _IEEE transactions on image processing_ , vol. 22, no. 1, pp. 119–133, 2013\.
  * [5] H. K. Rafsanjani, M. H. Sedaaghi, and S. Saryazdi, “An adaptive diffusion coefficient selection for image denoising,” _Digital Signal Processing_ , vol. 64, pp. 71–82, 2017.
  * [6] M. Lindenbaum, M. Fischer, and A. Bruckstein, “On gabor’s contribution to image enhancement,” _Pattern Recognition_ , vol. 27, no. 1, pp. 1–8, 1994\.
  * [7] M. Petrou and C. Petrou, _Image processing: the fundamentals_. John Wiley & Sons, 2010.
  * [8] L. Alvarez, P.-L. Lions, and J.-M. Morel, “Image selective smoothing and edge detection by nonlinear diffusion. ii,” _SIAM Journal on numerical analysis_ , vol. 29, no. 3, pp. 845–866, 1992.
  * [9] S. M. Smith and J. M. Brady, “Susan—a new approach to low level image processing,” _International journal of computer vision_ , vol. 23, no. 1, pp. 45–78, 1997.
  * [10] A. Fedorov, R. Beichel, J. Kalpathy-Cramer, J. Finet, J.-C. Fillion-Robin, S. Pujol, C. Bauer, D. Jennings, F. Fennessy, M. Sonka _et al._ , “3d slicer as an image computing platform for the quantitative imaging network,” _Magnetic resonance imaging_ , vol. 30, no. 9, pp. 1323–1341, 2012.
  * [11] A. Buades, B. Coll, and J.-M. Morel, “A non-local algorithm for image denoising,” in _Computer Vision and Pattern Recognition, 2005. CVPR 2005\. IEEE Computer Society Conference on_ , vol. 2. IEEE, 2005, pp. 60–65.
  * [12] P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” _IEEE Transactions on pattern analysis and machine intelligence_ , vol. 12, no. 7, pp. 629–639, 1990.
  * [13] C. A. Palma, F. A. Cappabianco, J. S. Ide, and P. A. Miranda, “Anisotropic diffusion filtering operation and limitations-magnetic resonance imaging evaluation,” _IFAC Proceedings Volumes_ , vol. 47, no. 3, pp. 3887–3892, 2014.
  * [14] C. S. Anand and J. S. Sahambi, “Wavelet domain non-linear filtering for mri denoising,” _Magnetic Resonance Imaging_ , vol. 28, no. 6, pp. 842–861, 2010\.
  * [15] J. V. Manjón, P. Coupé, A. Buades, D. L. Collins, and M. Robles, “New methods for mri denoising based on sparseness and self-similarity,” _Medical image analysis_ , vol. 16, no. 1, pp. 18–27, 2012.
  * [16] W. C. Wong, A. C. Chung, and S. C. Yu, “Trilateral filtering for biomedical images,” in _Biomedical Imaging: Nano to Macro, 2004. IEEE International Symposium on_. Citeseer, 2004, pp. 820–823.
  * [17] J. Xie, P.-A. Heng, and M. Shah, “Image diffusion using saliency bilateral filter,” _IEEE Transactions on Information Technology in Biomedicine_ , vol. 12, no. 6, pp. 768–771, 2008.
  * [18] G. Gerig, O. Kubler, R. Kikinis, and F. A. Jolesz, “Nonlinear anisotropic filtering of mri data,” _IEEE Transactions on medical imaging_ , vol. 11, no. 2, pp. 221–232, 1992.
  * [19] M. J. Black, G. Sapiro, D. H. Marimont, and D. Heeger, “Robust anisotropic diffusion,” _IEEE Transactions on image processing_ , vol. 7, no. 3, pp. 421–432, 1998.
  * [20] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Bm3d image denoising with shape-adaptive principal component analysis,” in _SPARS’09-Signal Processing with Adaptive Sparse Structured Representations_ , 2009.
  * [21] M. Yang, J. Liang, J. Zhang, H. Gao, F. Meng, L. Xingdong, and S.-J. Song, “Non-local means theory based perona–malik model for image denosing,” _Neurocomputing_ , vol. 120, pp. 262–267, 2013.
  * [22] J. Yuan, “Improved anisotropic diffusion equation based on new non-local information scheme for image denoising,” _IET Computer Vision_ , vol. 9, no. 6, pp. 864–870, 2015.
  * [23] V. Kamalaveni, R. A. Rajalakshmi, and K. Narayanankutty, “Image denoising using variations of perona-malik model with different edge stopping functions,” _Procedia Computer Science_ , vol. 58, pp. 673–682, 2015.
  * [24] F. Voci, S. Eiho, N. Sugimoto, and H. Sekibuchi, “Estimating the gradient in the perona-malik equation,” _IEEE Signal Processing Magazine_ , vol. 21, no. 3, pp. 39–65, 2004.
  * [25] C. Tsiotsios and M. Petrou, “On the choice of the parameters for anisotropic diffusion in image processing,” _Pattern recognition_ , vol. 46, no. 5, pp. 1369–1381, 2013.
  * [26] M. Maggioni, V. Katkovnik, K. Egiazarian, and A. Foi, “Nonlocal transform-domain filter for volumetric data denoising and reconstruction,” _IEEE transactions on image processing_ , vol. 22, no. 1, pp. 119–133, 2013\.
  * [27] J. V. Manjón, P. Coupé, and A. Buades, “Mri noise estimation and denoising using non-local pca,” _Medical image analysis_ , vol. 22, no. 1, pp. 35–47, 2015.
  * [28] C. A. Cocosco, V. Kollokian, R. K.-S. Kwan, G. B. Pike, and A. C. Evans, “Brainweb: Online interface to a 3d mri simulated brain database,” in _NeuroImage_. Citeseer, 1997.


# Continuum modeling perspectives of non-Fourier heat conduction in biological
systems

Á. Sudár1,2, G. Futaki1, R. Kovács1,2,3 1Department of Energy Engineering,
Faculty of Mechanical Engineering, BME, Budapest, Hungary
2Department of Theoretical Physics, Wigner Research Centre for Physics,
Institute for Particle and Nuclear Physics, Budapest, Hungary
3Montavid Thermodynamic Research Group

###### Abstract.

The thermal modeling of biological systems has increasing importance in
developing more advanced, more precise techniques such as ultrasound surgery.
One of the primary barriers is the complexity of biological materials: the
geometrical, structural, and material properties vary in a wide range, and
they depend on many factors. Despite these difficulties, there is a tremendous
effort to develop a reliable and implementable thermal model. In the present
paper, we focus on the continuum modeling of heterogeneous materials with
biological origin. There are numerous examples in the literature for non-
Fourier thermal models. However, as we realized, they are associated with a
few common misconceptions. Therefore, we first aim to clarify the basic
concepts of non-Fourier thermal models. These concepts are demonstrated by
revisiting two experiments from the literature in which the Cattaneo-Vernotte
and the dual phase lag models are utilized. Our investigation revealed that
using these non-Fourier models is based on misinterpretations of the measured
data, and the seeming deviation from Fourier’s law originates in the source
terms and boundary conditions.

Keywords: non-equilibrium thermodynamics, heat conduction, biological systems

##  1\. Introduction

Thermodynamics is one of the most crucial fundamental blocks of
interdisciplinary researches. It is proved to be useful in numerous coupled
problems such as thermo-electricity, thermo-diffusion, and thermo-acoustics
[1, 2, 3], thanks to the universal definition of the thermodynamic laws.
Regarding the topic in this paper, we consider only the I. and the II. laws of
thermodynamics for heat conduction in a continuum medium. Moreover, we also
assume the material to be rigid, although it is not characteristic of
biological systems. We want to keep the discussion as simple as possible to
clearly present the essential aspects of non-Fourier heat conduction for such
complex systems. Otherwise, mechanics must be included together with its
consequences: including source and coupling terms in the balance and
constitutive equations, respectively.

These equations serve the basis for continuum models, i.e., the I. law of
thermodynamics describes the balance of internal energy:

$\rho\partial_{t}e+\nabla\cdot\mathbf{q}=0,$ (1)

with $\mathbf{q}$ being the heat flux, $\rho$ is the mass density, $\nabla$ is
the usual nabla operator, $e=cT$ is the internal energy density with $c$
standing for the specific heat, and $T$ is the temperature. Also, the zero on
the right-hand side in (1) expresses the lack of source terms. This equation
requires a ‘closure,’ an expression between the temperature $T$ and the heat
flux $\mathbf{q}$. Such expression is called constitutive equation, describing
a material’s behavior, and strictly restricted by the II. law of
thermodynamics [4, 5]. In the simplest model, one can use Fourier’s law:

$\mathbf{q}=-\lambda\nabla T,$ (2)

in which $\lambda$ is the thermal conductivity. However, as observed in many
cases [6, 7, 8], this relation is restricted on a certain time, and spatial
scale [9], and can lose its validity. Therefore, one must find a suitable,
thermodynamically compatible generalization.

One often referred example in papers on bioheat transfer is the dual phase lag
(DPL) model,

$\mathbf{q}(x,t+\tau_{q})=-\lambda\nabla T(x,t+\tau_{T}),$ (3)

introduced by Tzou [10]. Here, $\tau_{q}$ and $\tau_{T}$ are the corresponding
relaxation times (time lags), vaguely described in many papers. According to
the original argument of Tzou, $\tau_{q}$ can be either smaller or higher than
$\tau_{T}$, depending on the phenomenon. Sadly, numerous works concluded that
these relaxation times are not independent of each other and cannot be
arbitrary. Furthermore, the DPL model can violate the basic physical
principles [11, 12, 13, 14, 15], depending on $\tau_{q}$ and $\tau_{T}$. This
problem originates in the Taylor series expansion of equation (3), considered
up to arbitrary orders, which is not a thermodynamically compatible way to
generalize and derive constitutive relations.

Instead of the DPL model, we propose a similar one with stronger physical
background, called Guyer-Krumhansl (GK) equation. In what follows, we first
shortly gather the thermodynamic background of the GK model. After, we briefly
summarize the known experiences on non-Fourier heat conduction experiments,
restricting ourselves on heterogeneous macro-scale materials and room
temperature conditions. Then, we revisit two experiments in which both Fourier
and non-Fourier models are utilized, showing common misconceptions about the
possible non-Fourier effects. With revisiting these experiments, we aim to
highlight the less-known properties of non-Fourier models and present our
perspectives of modeling heterogeneous materials.

##  2\. Thermodynamics of constitutive relations

As we previously mentioned, the II. law of thermodynamics restricts the
constitutive equations. In fact, it offers a constructive methodology to find
a ‘closure’ and derive the possible set of constitutive equations. The
framework of Classical Irreversible Thermodynamics [1, 16] formulate the
entropy density $s$ as a potential function of the state variables. For heat
conduction, it reads $s=s(e)$, and the Gibbs relation
$\textrm{d}e=T\textrm{d}s$ expresses the functional relationship between $e$
and $s$. The II. law can be formulated similarly to the I. law:

$\rho\partial_{t}s+\nabla\cdot\mathbf{J}_{s}=\sigma_{s}\geq 0,$ (4)

where $\mathbf{J}_{s}$ stands for the entropy flux, and $\sigma_{s}$ is the
entropy production. Eventually, that inequality leads to the set of
constitutive equations. ‘Set’ is intentional: there are infinitely many
solutions of such inequality. For instance, in case of $s=s(e)$ with
$\mathbf{J}_{s}=\mathbf{q}/T$, it yields

$\sigma_{s}=\mathbf{q}\nabla\frac{1}{T}\geq 0,$ (5)

for which the simplest linear solution is the Fourier’s law (2). Although many
other solutions exist in the nonlinear regime (i.e., nonlinear functions
respect to $\nabla(1/T)$), but we keep our focus on the linear ones. Such
linear solution of (4) still allows to have a state variable-dependent
transport coefficient in the constitutive equation such as
$\lambda=\lambda(T)$.

In the non-Fourier case, either the set of state variables can be extended,
e.g., with the heat flux $\mathbf{q}$, or the entropy flux $\mathbf{J}_{s}$
can be more general [17, 18]. Among the various approaches, the non-
equilibrium thermodynamics with internal variables (NET-IV) is applied here
due to its useful generality [19, 20, 21, 22, 23]. In other words, while the
kinetic theory-based Rational Extended Thermodynamics (RET) [24, 25] builds
its heat conduction models on phonon interactions, the continuum formulation
of NET-IV leaves the particular transport mechanism aside. Consequently, in
the RET model, the transport coefficients can be calculated, while in the
framework of NET-IV, these coefficients can be determined through a fitting
procedure. Therefore, while it is easier to determine the transport
coefficients in RET, it is restricted to systems with high Knudsen number (Kn
$\geq 0.1$), characteristic for nano or low-temperature systems. On the
contrary, the NET-IV model has more freedom by having tunable parameters, and
that generality allows us to apply the generalized equation for room
temperature problems. Particularly for the Guyer-Krumhansl equation, in the
NET-IV approach we have $s=s(e,\mathbf{q})$ and
$\mathbf{J}_{s}=\mathbf{B}\mathbf{q}$, in which $\mathbf{B}$ is called current
or Nyíri multiplier [17], helping to introduce the so-called nonlocal terms
(e.g., the $\partial_{xx}q$). For the detailed derivation, explanation, and
background, we refer to [26]. The GK equation for a one-dimensional situation
in Cartesian coordinates is

$\tau_{q}\partial_{t}q+q=-\lambda\partial_{x}T+\kappa^{2}\partial_{xx}q,$ (6)

where $\kappa^{2}$ is a sort of ‘dissipation parameter’, usually connected to
the mean free path in the kinetic theory. Eq. (6) has two special cases. The
first one is the Cattaneo-Vernotte (CV) equation (with $\kappa^{2}=0$), the
other is a less-known one, called Nyíri equation (with $\tau_{q}=0$):

$q=-\lambda\partial_{x}T+\kappa^{2}\partial_{xx}q,$ (7)

which is analogous with the Brinkman equation (the nonlocal generalization of
Darcy’s law). Although the CV equation is well-applicable for dissipative wave
propagation (‘second sound’), it is not observed for macro-scale heterogeneous
materials on room temperature. At that point, we note that many other
procedures are possible, depending on the purpose and on the thermodynamic
framework in which the problem is formulated [27, 28, 29, 30, 31, 32].

###  2.1. Why the Guyer-Krumhansl model?

In our recent heat pulse experiments, we investigated the transient thermal
behavior of several heterogeneous materials such as metal foams, rocks, and
capacitor specimen [33, 34]. Among them, the capacitor sample has the most
regular structure in which good conductor and insulator layers are arranged
periodically, two materials with different characteristic thermal time scales.
The presence of multiple time scales attributes the heterogeneous materials,
in general. In the $T$-representation of the GK model,

$\tau_{q}\partial_{tt}T+\partial_{t}T=\alpha\partial_{xx}T+\kappa^{2}\partial_{xxt}T,$
(8)

the thermal diffusivity $\alpha=\lambda/(\rho c)$ is comparable with the ratio
of $\kappa^{2}/\tau_{q}$. When $\kappa^{2}/\tau_{q}=\alpha$, it is called
Fourier resonance, and it recovers the Fourier behavior exactly. Otherwise,
non-Fourier solutions appear, and these GK parameters represent the existence
of parallel heat conduction time scales. On the contrary to the usual belief
that merely wave propagation is possible when non-Fourier effects come into
the picture, we found the so-called over-diffusion for which
$\kappa^{2}/\tau_{q}>\alpha$, and Figure 1 shows its characteristics.
Remarkably, it seems similar to the classical case, but the Fourier equation
cannot describe the data with acceptable accuracy. Moreover, both the thermal
diffusivity and the non-Fourier effects can be size-dependent [35]. Overall,
based on our experimental findings, the GK equation stands as the next
candidate beyond the Fourier equation. Neither of its special cases provided
acceptable predictions.

Figure 1. The signs of non-Fourier heat conduction for heterogeneous materials
at room temperature.

##  3\. Bioheat models

Roughly speaking, the related heat conduction models belong to two groups,
based on the considered constitutive equation, but all aim to predict the
tissue temperature. In the first group, Fourier’s law is utilized with various
biologically motivated extensions, mostly being source terms in the balance of
internal energy. Particular examples are the internal heat transfer by blood
perfusion, which appeared first in the Pennes’ model [36]. Later, more and
more physiological details are tried to be implemented. A few well-known
examples are related to the coupled two-temperature Chen-Holmes [37] and the
Weinbaum-Jiji-Lemons models [38]. However, the methodology of implementing the
detailed structure (e.g., the artery-vein pairs, various tissue structures)
does not seem to be advantageous due to the more complicated equations, the
unknown variables, and limited validity region. Also, in some cases, the
extensions are physically contradictory: in Wulff’s model [39], the mean blood
velocity $\mathbf{u}_{b}$ explicitly appears:

$\mathbf{q}=-\lambda_{t}\nabla T_{t}+\rho_{b}f_{b}\mathbf{u}_{b}$ (9)

with $f_{b}$ being the blood enthalpy, and the indicies $t$ and $b$ are
referring to the tissue and blood properties. Such coupling between
$\mathbf{q}$ and $\mathbf{u}$ is not possible, unfortunately, due to
objectivity reasons [40, 41, 42]. Similar shortcomings appear in Klinger’s
approach [43, 44].

The second group of biologically motivated heat conduction models consists of
the DPL and the CV equations in the vast majority of problems, with keeping
specific source terms. Sadly, the DPL model suffers from physical and
mathematical shortcomings, and the CV equation predicts waves, which are not
observed so far.

Beyond these models, there are more specific ones in the literature. For
instance, finding the thermal response of the cornea is a long-lasting
question. In this respect, we mention the work of Taflove and Brodwin [45] due
to their remarkable idea: it is more advantageous to use ‘effective’ models.
Effective in the sense that the complex heterogeneous structure is substituted
with a homogeneous one; consequently, the model parameters are also effective,
i.e., ‘averaging’ the complete thermal behavior. This methodology agrees with
the GK model’s perspective: it introduces two new, physically strongly
motivated parameters, which characterize the entire medium, despite its
particular heterogeneities. From a purely thermodynamic point of view, rocks
do not differ significantly from biological materials: they both porous,
having various heterogeneities and irregularities, thus parallel time scales
are present. If such material shows non-Fourier thermal behavior, then it
should be observable when all biological specific effects are excluded, and no
heat source is present. In other words, the difficulties with complicated heat
sources and boundary conditions must be separated from the constitutive model.
It stands as the motivation to revisit two preceding experiments in which
these aspects appear and could lead to the misinterpretation of the observed
data.

##  4\. Experiment I.: processed meat samples

The first experiment we re-evaluate is performed by Tang et al. [46] on
processed meat samples, including fat as inclusions. They prepared three
samples with $2$, $3$ and $4$ mm thickness, and all are having the same $10$
mm diameter. They applied a heat pulse on the front side uniformly, and it
lasts $1$ s while they registered the rear side temperature history for $90$
s. Figure 2 shows their results partially; for their complete investigation,
we refer to [46]. According to their interpretation, convection cannot be
responsible for the measured deviation from Fourier’s law. They investigated
two approaches to explain this phenomenon. On one hand, they proposed to use
the DPL model with adiabatic (zero flux) boundaries, hence neglecting the
possible heat transfer to the environment. On the other hand, they introduced
the size of the fat inclusion as a variable parameter and performed a detailed
numerical analysis using the Fourier equation.

In the end, they concluded that the Fourier equation seems to be a better
choice together with these new parameters than the DPL model. These models
differ mostly at the end of the temperature history. Despite the valuable
insight they provided by these experiments, there are a few
misinterpretations, which we found to be important presenting here.

Figure 2. The observed temperature history for the thinnest sample ($2$ mm).
A) presents the solutions for Fourier equation with various heat transfer
coefficients. B) shows the outcomes of other approaches. The complete figures
and for further information, we refer to [46].

We re-evaluate that data using a novel algorithm developed explicitly for the
GK equation using its analytical solution [47]. Shortly, this evaluation
procedure first finds the heat transfer coefficient and the thermal
diffusivity. If needed, it continues to find the GK parameters ($\tau_{q}$ and
$\kappa^{2}$). We found that omitting the heat transfer to the environment
yields significant errors in the modeling, i.e., Fourier’s law and convection
boundary condition provide a suitable model to adequately describe the
observed data. In other words, no generalized heat conduction equation is
needed, merely a proper optimization on the effective thermal diffusivity.
Consequently, despite that they found the relaxation times in the
mathematically valid domain ($\tau_{q}>\tau_{T}$), it remains misleading to
interpret the temperature history as proof for the wave nature of heat
conduction under such conditions.

Figure 3. Comparing the Fourier model to the sample with $4$ mm thickness. The
measured data is recovered using a plot digitizer. Figure 4. Comparing the
Fourier model to the sample with $2$ mm thickness. The measured data is
recovered using a plot digitizer.

Nevertheless, the size-dependent behavior is indeed present within this
spatial scale ($2-4$ mm) and requires further investigation, and probably new
experiments would be helpful for the literature. Table 1 summarizes the
parameters, and Figures 3-4 show the fitting of the Fourier equation.

Sample thickness |  | Thermal diffusivity
---
$10^{-7}$ [m2/s]
$2$ mm | $0.66$
$3$ mm | $0.73$
$4$ mm | $1.04$
Table 1. The found thermal diffusivities with the Fourier equation.

##  5\. Experiment II.: skin tissue samples

In this section, we revisit the experiment of Jaunich et al. [48], in which
they prepared various cylindrical samples made from the mixture of araldite,
DDSA resin, and DMP-30 epoxies in order to mimic the behavior of real tissue,
moreover, they added titanium dioxide to this mixture as a ‘scatterer’ for the
laser light. They assumed that this mixture also thermally substitutes the
real skin tissue, therefore utilizing the same optical and thermal parameters,
detailed in [48]. A skin tissue has a promising structure in regard to non-
Fourier heat conduction due to the parallel time scales (the presence of
various layers and heat transfer mechanisms). Nevertheless, its observation
might also require the heat flux to be parallel with these layers [33].
Contrary to the previous experiment, the thermal excitation now lasts $10$ s
instead of $1$ s, and registering the temperature history for $10$ s.

Regarding the aspects of non-Fourier heat conduction, all these are extremely
important. The material structure has no sharp interfaces inside due to the
relatively homogenous mixture of epoxy components, therefore it is different
from real heterogeneous porous materials. Furthermore, the excitation is
unusual in the sense that it is not ‘pulse-like’ but much longer,
significantly absorbed under the surface, and the temperature history is
measured only in the period of laser heating, and there is no available data
after. Consequently, the entire experiment is heat source dominated, which
hides heat conduction and, eventually, almost any other thermal transport
effects.

According to their findings, Fourier’s law is seemingly inadequate for this
problem since it predicts the temperature field quite far from the measured
one. Figure 5 presents the measured and fitted results partially, only the
data, which we also tried to reproduce. The assumption that other constitutive
equations should be used instead of Fourier’s law requires accurate knowledge
of any other influencing property, most importantly, on the internal heat
generation caused by the laser light. In their study, the CV (hyperbolic)
equation is tested and found to be better than the Fourier equation, providing
more accurate predictions. Unfortunately, even the CV equation fails to offer
the proper values, and at some points, it is still the same as the Fourier
equation.

Figure 5. The partial results of Jaunich et al., for their detailed
description and for the complete figures, we refer to [48]. Triangles denote
the measured temperature increase. A) presents the temperature history at $2$
mm depth on the axis of the cylinder. B) and C) show the radial and axial
temperature distribution at specific points.

We also tried to model this experiment by solving the Fourier, CV, and GK
equations in cylindrical coordinates based on the scheme published in [9]. We
found that the problem indeed originates in the heat source since neither the
CV nor the GK model could improve the fitting, as it was supposed from the
beginning. Therefore we tried various volumetric heat generation models from
which the following one provides the closest values with the Fourier equation:

$Q_{v}=\left\\{\begin{array}[]{ll}\displaystyle\frac{q_{m}}{2}\left(1+\cos\big{(}\pi
r_{e}/r_{p}\big{)}\right)&\ r_{{e}}<r_{{p}}\\\ 0&\
\textrm{else}\end{array}\right.,$ (10)

in which

$r_{e}=\sqrt{r^{2}+(z-z_{c})^{2}}$ (11)

is the effective radius where heat generation $Q_{v}$ occurs, $r_{p}$ stands
for the pulse radius; $z_{c}$ appears in the case when the center of heat
generation is under the surface, in our case, it is found to be $z_{c}=0$;
$q_{m}$ is the maximum heat generation. The value of $q_{m}$ is assumed to be
the only adjustable parameter due to the unknown amount of absorbed heat, and
it is found to be $q_{m}=1.7\cdot 10^{7}$ W/m3, which means that around the
half of the irradiated energy is turned into heat generation. Furthermore, $z$
and $r$ are the axial and radial coordinates, respectively. Although this is
an empirical model, it is motivated by statistical physics: photons propagate
with equal probability in every direction after enough number of scattering,
resulting in a spherically symmetric absorption. Figures 6 and 7 present the
best outcome we managed to achieve.

Figure 6. Comparing the Fourier model to the measured values in both axial and
radial directions. The measured data is recovered using a plot digitizer.
Figure 7. Comparing the Fourier model to the temperature history recorded at
$z=2$ mm. The measured data is recovered using a plot digitizer.

As it is apparent from Fig. 6, we also applied a convection boundary on the
surface, which was irradiated by the laser pulse, the heat transfer
coefficient is assumed to be $10$ W/(m2K). Unfortunately, the measured
transient characteristic at $z=2$ mm is still uncovered, and non-Fourier
models cannot offer a better description, in agreement with [48].

##  6\. Discussion

The thermal modeling of biological materials is a challenging task and has
increasing attention due to its numerous applications. In the present paper,
we aimed to offer a brief overview of the essential aspects of non-Fourier
equations in continuum thermodynamics, focusing on their observability in room
temperature experiments.

First, we want to emphasize that one must pursue to separate the heat sources
from the heat conduction effects as much as possible. Although volumetric heat
generation is inevitable in many practical situations and also has a
significant influence in medical diagnostics [49, 50, 51], it could dominate
the time evolution of the temperature field. Therefore, the experiments, which
aim to measure biological material’s thermal properties, should be designed to
let the heat conduction be the dominant heat transport mechanism instead of
including complex heat source and heat generation mechanisms. In such a
preferable situation, it is easier to decide whether we observed a non-Fourier
phenomenon or not; consequently, the role of a constitutive equation becomes
apparent in the thermal model. Moreover, the non-Fourier thermal parameters
are more comfortable to determine. Otherwise, one can easily misinterpret the
measurements and concluding in a different outcome, and it is recommended to
always test all the possibilities in a model before implementing a generalized
constitutive equation.

One must be careful since the observation of non-Fourier propagation depends
both on the boundary conditions and material properties. In other words, it
requires the coexistence of suitable spatial and time scales, probably having
a size-dependent behavior as well. This field is still under development, and
there are countless questions, which wait to be solved in the future. Also, as
our re-evaluation attempts are standing here, the Fourier equation can explain
most of the observations, and most of the seeming deviation disappeared.

Overall, we aim to collect and call attention to these modeling aspects. The
research field of non-Fourier heat conduction is vast, continuously growing,
and increasingly challenging to keep up with the up-to-date results, mainly
because these constitutive equations require a deeper understanding of
boundary conditions as well. The classical definitions established for the
Fourier equation do not work, and even for a thermodynamically strongly
motivated model (such as the GK equation) can lead to unphysical, false
solutions with inadequate boundary conditions. Therefore, despite the
attractive results of non-equilibrium thermodynamics, it is suggested to be
careful with generalized models.

## Acknowledgement

We dedicate our paper to the memory of our respected colleague, Prof. José
Casas-Vázquez, who passed away recently.
The research reported in this paper and carried out at BME has been supported
by the grants National Research, Development and Innovation Office-NKFIH FK
134277, and by the NRDI Fund (TKP2020 NC, Grant No. BME-NC) based on the
charter of bolster issued by the NRDI Office under the auspices of the
Ministry for Innovation and Technology. This paper was supported by the János
Bolyai Research Scholarship of the Hungarian Academy of Sciences.

## References

  * [1] S. R. de Groot and P. Mazur.  Non-Equilibrium Thermodynamics.  Dover Publications, 1963.
  * [2] L. Onsager.  Reciprocal relations in irreversible processes. I.  Physical Review, 37(4):405, 1931.
  * [3] L. Onsager.  Reciprocal relations in irreversible processes. II.  Physical Review, 38(12):2265, 1931.
  * [4] T. Matolcsi.  Ordinary Thermodynamics.  Akadémiai Kiadó, 2004.
  * [5] A. Berezovski and Ván P.  Internal Variables in Thermoelasticity.  Springer, 2017.
  * [6] H. E. Jackson, C. T. Walker, and T. F. McNelly.  Second sound in NaF.  Physical Review Letters, 25(1):26–28, 1970.
  * [7] V. Narayanamurti, R. C. Dynes, and K. Andres.  Propagation of sound and second sound using heat pulses.  Physical Review B, 11(7):2500–2524, 1975.
  * [8] T. F. McNelly.  Second Sound and Anharmonic Processes in Isotopically Pure Alkali-Halides.  1974\.  Ph.D. Thesis, Cornell University.
  * [9] V. Józsa and R. Kovács.  Solving Problems in Thermal Engineering: A Toolbox for Engineers.  Springer, 2020.
  * [10] D. Y. Tzou.  A unified field approach for heat conduction from macro- to micro-scales.  Journal of Heat Transfer, 117(1):8–16, 1995.
  * [11] M. Fabrizio and B. Lazzari.  Stability and second law of thermodynamics in dual-phase-lag heat conduction.  International Journal of Heat and Mass Transfer, 74:484–489, 2014\.
  * [12] R. Quintanilla and R. Racke.  Qualitative aspects in dual-phase-lag heat conduction.  Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 463(2079):659–674, 2007.
  * [13] R. Kovács and P. Ván.  Thermodynamical consistency of the Dual Phase Lag heat conduction equation.  Continuum Mechanics and Thermodynamics, pages 1–8, 2017.
  * [14] S. A. Rukolaine.  Unphysical effects of the dual-phase-lag model of heat conduction.  International Journal of Heat and Mass Transfer, 78:58–63, 2014\.
  * [15] S. A. Rukolaine.  Unphysical effects of the dual-phase-lag model of heat conduction: higher-order approximations.  International Journal of Thermal Sciences, 113:83–88, 2017.
  * [16] G. Lebon.  From classical irreversible thermodynamics to extended thermodynamics.  Acta Physica Hungarica, 66(1-4):241–249, 1989.
  * [17] B. Nyíri.  On the entropy current.  Journal of Non-Equilibrium Thermodynamics, 16(2):179–186, 1991\.
  * [18] M. Szücs, R. Kovács, and S. Simić.  Open mathematical aspects of continuum thermodynamics: Hyperbolicity, boundaries and nonlinearities.  Symmetry, 12:1469, 2020.
  * [19] T. Fülöp, Cs. Asszonyi, and P. Ván.  Distinguished rheological models in the framework of a thermodynamical internal variable theory.  Continuum Mechanics and Thermodynamics, 27(6):971–986, 2015.
  * [20] A. Berezovski, J. Engelbrecht, and G. A. Maugin.  Thermoelasticity with dual internal variables.  Journal of Thermal Stresses, 34(5-6):413–430, 2011.
  * [21] P. Ván, A. Berezovski, and J. Engelbrecht.  Internal variables and dynamic degrees of freedom.  Journal of Non-Equilibrium Thermodynamics, 33(3):235–254, 2008\.
  * [22] G. A. Maugin and W. Muschik.  Thermodynamics with internal variables. Part I. General concepts.  Journal of Non-Equilibrium Thermodynamics, 19(3):217–249, 1994\.
  * [23] G. A. Maugin and W. Muschik.  Thermodynamics with internal variables. Part II. Applications.  Journal of Non-Equilibrium Thermodynamics, 19(3):250–289, 1994\.
  * [24] I. Müller and T. Ruggeri.  Rational Extended Thermodynamics.  Springer, 1998.
  * [25] T. Ruggeri and M. Sugiyama.  Rational extended thermodynamics beyond the monatomic gas.  Springer, 2015.
  * [26] P. Ván.  Weakly nonlocal irreversible thermodynamics—the Guyer–Krumhansl and the Cahn–Hilliard equations.  Physics Letters A, 290(1-2):88–92, 2001.
  * [27] V. A. Cimmelli.  Different thermodynamic theories and different heat conduction laws.  Journal of Non-Equilibrium Thermodynamics, 34(4):299–333, 2009\.
  * [28] I. Carlomagno, A. Sellitto, and V. A. Cimmelli.  Dynamical temperature and generalized heat-conduction equation.  International Journal of Non-Linear Mechanics, 79:76–82, 2016.
  * [29] M. Grmela.  Generic guide to the multiscale dynamics and thermodynamics.  Journal of Physics Communications, 2(3):032001, 2018.
  * [30] M. Grmela, G. Lebon, and C. Dubois.  Multiscale thermodynamics and mechanics of heat.  Physical Review E, 83(6):061134, 2011.
  * [31] D. Jou, J. Casas-Vazquez, and G. Lebon.  Extended irreversible thermodynamics revisited (1988-98).  Reports on Progress in Physics, 62(7):1035, 1999.
  * [32] M. Sauermoser, S. Kjelstrup, N. Kizilova, B. G. Pollet, and E. G. Flekkøy.  Seeking minimum entropy production for a tree-like flow-field in a fuel cell.  Phys. Chem. Chem. Phys., 22(13):6993-7003, 2020.
  * [33] S. Both, B. Czél, T. Fülöp, Gy. Gróf, Á. Gyenis, R. Kovács, P. Ván, and J. Verhás.  Deviation from the Fourier law in room-temperature heat pulse experiments.  Journal of Non-Equilibrium Thermodynamics, 41(1):41–48, 2016.
  * [34] P. Ván, A. Berezovski, T. Fülöp, Gy. Gróf, R. Kovács, Á. Lovas, and J. Verhás.  Guyer-Krumhansl-type heat conduction at room temperature.  EPL, 118(5):50005, 2017.  arXiv:1704.00341v1.
  * [35] T. Fülöp, R. Kovács, Á. Lovas, Á. Rieth, T. Fodor, M. Szücs, P. Ván, and Gy. Gróf.  Emergence of non-Fourier hierarchies.  Entropy, 20(11):832, 2018.  ArXiv: 1808.06858.
  * [36] H. H. Pennes.  Analysis of tissue and arterial blood temperatures in the resting human forearm.  Journal of Applied Physiology, 1(2):93–122, 1948.
  * [37] M. M. Chen and K. R. Holmes.  Microvascular contributions in tissue heat transfer.  Annals of the New York Academy of Sciences, 335(1):137–150, 1980\.
  * [38] S. Weinbaum, L. M. Jiji, and D. E. Lemons.  Theory and experiment for the effect of vascular microstructure on surface tissue heat transfer—Part I: Anatomical foundation and model conceptualization.  Journal of Biomechanical Engineering, 106(4):321–330, 1984.
  * [39] W. Wulff.  The energy conservation equation for living tissue.  IEEE Transactions on Biomedical Engineering, 6(BME-21):494–495, 1974.
  * [40] W. Muschik.  Objectivity and frame indifference, revisited.  Archives of Mechanics, 50(3):541–547, 1998.
  * [41] T. Fülöp.  Objective thermomechanics.  arXiv preprint arXiv:1510.08038, 2015.
  * [42] T. Matolcsi and P. Ván.  Can material time derivative be objective?  Physics Letters A, 353(2):109–112, 2006.
  * [43] H. G. Klinger.  Heat transfer in perfused biological tissue—I: General theory.  Bulletin of Mathematical Biology, 36:403–415, 1974.
  * [44] H. G. Klinger.  Heat transfer in perfused biological tissue—II: The “macroscopic” temperature distribution in tissue.  Bulletin of Mathematical Biology, 40(2):183–199, 1978.
  * [45] A. Taflove and M. E. Brodwin.  Computation of the electromagnetic fields and induced temperatures within a model of the microwave-irradiated human eye.  IEEE Transactions on Microwave Theory and Techniques, 23(11):888–896, 1975.
  * [46] D. Tang, N. Araki, and N. Yamagishi.  Transient temperature responses in biological materials under pulsed IR irradiation.  Heat and Mass Transfer, 43(6):579–585, 2007.
  * [47] A. Fehér and R. Kovács.  Novel evaluation method for non-Fourier effects in heat pulse experiments.  2021\.  arXiv:2101.01123.
  * [48] M. Jaunich, S. Raje, K. Kim, K. Mitra, and Z. Guo.  Bio-heat transfer analysis during short pulse laser irradiation of tissues.  International Journal of Heat and Mass Transfer, 51(23):5511–5521, 2008.
  * [49] G. Tambave, J. Alme, G.G. Barnaföldi, R. Barthel, A. van den Brink, et al.  Characterization of monolithic cmos pixel sensor chip with ion beams for application in particle computed tomography.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 958:162626, 2020\.  Proceedings of the Vienna Conference on Instrumentation 2019.
  * [50] J. Alme, G. G. Barnaföldi, R. Barthel, et al.  A high-granularity digital tracking calorimeter optimized for proton ct.  Frontiers in Physics, 8:460, 2020.
  * [51] A. Andreozzi, L. Brunese, M. Iasiello, C. Tucci, and G. P. Vanoli.  Bioheat transfer in a spherical biological tissue: a comparison among various models.  Journal of Physics: Conference Series, 1224(1):012001, 2019.