ON THE APPLICATION OF THE FOURIER METHOD TO SOLVE THE PROBLEM OF CORRECTION OF THERMOGRAPHIC IMAGES Текст научной статьи по специальности «Математика»

Discrete & Continuous Models

#& Applied Computational Science 2022, 30 (3) 205-216

ISSN 2658-7149 (online), 2658-4670 (print) http://journals-rudn-ru/miph

Research article

UDC 519.6

PACS 07.05.Tp, 02.60.Pn, 02.70.Bf DOI: 10.22363/2658-4670-2022-30-3-205-216

On the application of the Fourier method to solve the problem of correction of thermographic images

Obaida Baaj

Peoples' Friendship University of Russia (RUDN University), 6, Miklukho-Maklaya St., Moscow, 117198, Russian Federation

(received: March 18, 2022; revised: March 28, 2022; accepted: August 8, 2022)

Abstract. The work is devoted to the construction of computational algorithms implementing the method of correction of thermographic images. The correction is carried out on the basis of solving some ill-posed mixed problem for the Laplace equation in a cylindrical region of rectangular cross-section. This problem corresponds to the problem of the analytical continuation of the stationary temperature distribution as a harmonic function from the surface of the object under study towards the heat sources. The cylindrical region is bounded by an arbitrary surface and plane. On an arbitrary surface, a temperature distribution is measured (and thus is known). It is called a thermogram and reproduces an image of the internal heat-generating structure. On this surface, which is the boundary of the object under study, convec-tive heat exchange with the external environment of a given temperature takes place, which is described by Newton's law. This is the third boundary condition, which together with the first boundary condition corresponds to the Cauchy conditions — the boundary values of the desired function and its normal derivative. The problem is ill-posed. In this paper, using the Tikhonov regularization method, an approximate solution of the problem was obtained, stable with respect to the error in the Cauchy data, and which can be used to build effective computational algorithms. The paper considers algorithms that can significantly reduce the amount of calculations.

Key words and phrases: thermogram, ill-posed problem, Cauchy problem for the Laplace equation, integral equation of the first kind, Tikhonov regularization method

1. Introduction

Improving the quality and information content of images obtained by thermal imaging methods using a thermal imager that registers thermal electromagnetic radiation from the surface of the object under study in the infrared range by their mathematical (digital) processing is an urgent problem. In particular, in medicine, thermal imaging has become an effective diagnostic tool [1-4]. The image on the thermogram, which is a visualization of the temperature distribution on the surface of the patient's body, makes it possible

© BaajO., 2022

This work is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.org/licenses/by-nc/4.0/legalcode

@ CD®

to assess functional anomalies in the state of his internal organs. At the same time, the image on the thermogram in some cases turns out to be somewhat distorted due to the processes of thermal conductivity and heat exchange. The paper proposes a method of image correction on a thermogram within a certain mathematical model. As an adjusted thermogram, the image of the temperature field on the plane near the density of heat sources is considered as more accurately transmitting the image of heat sources. It is proposed to obtain this field as a result of the continuation (similar to the continuation of gravitational fields in geophysics problems [5]) of the temperature distribution from the surface from which the initial thermogram is taken. The problem under consideration is ill-posed, since small errors in the initial data (the initial thermogram) may correspond to significant errors in solving the inverse problem. To construct its stable approximate solution, the Tikhonov regularization method [6] is used.

2. Mathematical model and problem statement

Let's consider a physical and mathematical model, in which we set the task of continuing from the boundary of the stationary temperature distribution.

The physical model is a homogeneous heat-conducting body in the form of a rectangular cylinder, bounded by the surface S and containing heat sources with a time-independent density function that create a stationary temperature distribution in the body. We associate the density function of heat sources with the object under study. We assume that a given temperature distribution (equal to zero) is maintained on the lateral faces of the cylinder, and on the surface S there is convective heat exchange with the external environment of temperature U0, described by Newton's law, according to which the density of the heat flux at the point of the surface S is directly proportional to the temperature difference inside and outside.

Let's move on to the mathematical model. In a rectangular cylinder

D^ = {(x,y,z) ■ 0 < x <lx, 0 <y <ly, -to < z < to} C K3 (1) consider a cylindrical region

D(F, to) = {(x, y,z) ■ 0 <x <lx ,0<y<ly, F(x, y) < z < to}, (2) limited by the surface

Let r be the set of side faces of the domain D(F,to). In the domain D(F, to) we consider the following mixed boundary value problem for the Laplace equation

S = {(x,y,z) ■ 0 < x < lx, 0 < y < ly, z = F(x, y) < H}.

(3)

We'll assume that we also know that

ai < F(x, y) <a2 < H, (x, y) e n, n = {(x,y) ■ 0<x<lx,0<y< ly}.

(4)

(5)

fAu(M) = p(M), M e D(F, œ) _u limited at z ^ œ.

(6)

The problem (6) corresponds to the steady-state temperature distribution created by heat sources with the distribution density function p, on the surface S — a third boundary condition is set corresponding to convective heat exchange with a medium of temperature U0 with a coefficient h, zero temperature is set at the boundary r.

We assume that the function p is such that the solution of the problem (6) exists in C2(D(F,o))^CX(D(F,<)). In particular, the solution of the problem (6) allows us to find , i.e. the temperature distribution of u on the surface S, which we will call a thermogram.

Now let the thermogram be obtained as a result of measurements. Let us now set the inverse problem. We set the problem of continuation of the temperature distribution from the surface towards the sources in order to obtain an adjusted thermogram as the temperature distribution u\z=H on the plane z = H, closer to the density carrier than the surface S.

We assume that the carrier of the function p is located in the domain z > H, then the solution of the problem (6) in the domain

D(F, H) = {(x,y,z)-. 0<x<lx,0<y<ly, F(x, y) < z < H} (7)

satisfies the Laplace equation. The set of side faces of the domain D(F, H) is denoted by rH.

Inverse problem. Let the function be given within the framework of the model (6)

It is required to find u\_ . Since the value of H sufficiently arbitrarily

defines the plane between the support of p and the surface S, then in fact the inverse problem consists in obtaining a solution u in the domain D(F, H) of the boundary value problem

(8)

fAu(M) = 0, MeD(F,H)

U\o = f,

(9)

We assume that the function f in (8), (9) is taken from the set of solutions to the direct problem (6), so the solution to the inverse problem exists in

C 2(D(F,H))^C 1(D(F,H)).

Note that in the problem (9) on the surface S of the form (3), Cauchy conditions are set, that is, the boundary values f of the desired function u and the values of its normal derivative are set, so the problem (9) has a unique solution. The boundary z = H of the domain D(F, H) is free and, thus, the problem (9) is unstable with respect to errors in the data, i.e. ill-posed.

The function u\z=H will be considered as an adjusted thermogram. Since the plane z = H is located closer to the support of density p, it should be expected that the corrected thermogram more accurately conveys information about the distribution of heat sources than the original thermogram.

Further we give an explicit representation of the exact solution of the problem (9).

3. Exact solution of the inverse problem

Based on the [7] scheme, an exact solution of the problem (9) is constructed in [8].

Let p(M, P) be the source function of the Dirichlet problem in the cylinder D™ :

Au(P) = -p(P), P e D™,

u\ = 0, u\ = 0, (10)

]x=0,lx , \y=o,iy , (10)

u ^ 0 at Izl ^ œ. In the domain zM < H in the cylinder (1), we introduce the notation

$(M) =

h(Uo - f(P)MM, P) - f(P)^-(M, P)

daP. (11)

In [8], the following representation of the solution of the problem is obtained (9)

u(M) = v(M) + $(M), MeD(F,H), (12)

where the function $ is calculated on the known functions f and fx, and the function v has the form:

nnx nmy -a)} sin —:— sin

v(M) = - ^ $nm(a) exp {knm(z-a)}

n,m=1 ^ ly (13)

M(x,y,z) e D(-œ,H),

where

n2 m2 \ ^ A .

Km =n[-2 + ) (14)

V lx 'Jy J

and &nm(a) — Fourier coefficients of the function &(M)

£ r \ 4 [ Knx . nmy

®nm (a) = T~T I Q(x,y,a) sin ——sin—— dxdy (15)

'JX Jy J Jx 'J'U

n(a) *

on the auxiliary plane:

n(a) = {(x, y,z) ■ 0 < x < lx, 0 < y < ly, z = a}, a < a1. (16)

For a $ function of the form (11) considering that daP = (xP, yp)dxpdyp, where the normal n1 to the surface S is calculated by the formula

n = grad (F(x,y) — z) = VxyF- k, ^ = |nI, (17)

we will use the representation

$(M) = i [h(Uo -f(xp,yp))<p(M,P)\p€Sni(xP,yP)-

n

-f(xP,yP)(ni,VMM,P))\peS dxpdyp. (18)

When calculating the function $(M)IMeU(a) on the rectangle n(a) for the source function tp(M,P), you can use the formula

<p(M,P) =

2 e-knm\zM-zp\ nnxM KmyM nnxp nmyP , ,

= —— > -;-sin—--sin—--sin—-—sin—--, (19)

1 1 Z_/ b 1 1 1 J

'Jx by n m=1 nm 'Jx by 'Jx by

which for zM = a and P G S takes the form

2 ^ e-knm{Fixp,yp)-a)

*(M,p) = TT E --1-x

hx hy n,m=1 nrn

nnxM KmyM nnxp nmyP , x x sin —-—— sin yM sin —-—— sin —. (20)

ix Ly ix Ly

The series converges uniformly, since the exponent is estimated by exp{-knm(a1 — a)}. When calculating the function $ in (12), the source function at a2 < zM < H and P G S takes the form

2 ^ {ZM — P{xp,yp})

*(M,p) = TT E --i-x

hx hy n,m=1 nrn

nnxM KmyM nnxp nmyP , x x sin —-—— sin yM sin —-—— sin —. (21)

Lx by ix Ly

The series converges uniformly on any fixed plane zM = const, since the exponent is estimated by exp{—A;nm(zM — a2)}, that is important for applications. At the points zM < a2, the source function can be calculated by the reflection method.

4. Construction of an approximate solution to the problem

Let the function f in the problem (9) be given with an error, that is, instead of f, the function fs is given, so that

IfS -fhap) (22)

In this case, the function (11) is calculated approximately

$(M)= i [h(U0 -f5(xP,yP))^(M,P)\ppsn1 (xP, yP)-

n

-fs(xP,yP)(n1 ,VMM,P))\peS dxpdyp. (23)

The approximate solution to the problem (9) is constructed using the Tikhonov regularization method [6] and in accordance with (12) has the form

usa (M) = v5a (M) + ¥ (M), MeD(F,H), (24)

where is a function of the form (23) and

VS (M) = - V (a) exp{knm(zM - a)} sin ™xM sin KmyM

a ( ) ^ 1 + a exp{2knm (H-a)} lx ly . ()

Note that the members of the series (25) differs from the members of the series (13) by the regularizing factor (1 + aexp{2knm(H - a)})-1, ensuring the convergence of the series.

In the numerical solution, the bulk of the calculations is related to the calculation of the Fourier coefficients of the function by the formula (15). The next section is devoted to the calculation of Fourier coefficients with a significant reduction in the amount of calculations.

5. Calculation of Fourier coefficients

As follows from the formulas (15), (23), (20), when calculating the Fourier coefficient for each pair of indices n and m, a superposition of the following calculations is required: summation of the series for p, integration on the surface S, integration on the rectangle n(a). Thus, when discretizing [9] the problem (Nx points on the variable x, Ny points on the variable y) when calculating Fourier coefficients, about 0(NXNy)4 operations are required. This is the largest volume of operations when constructing a solution to the problem (9), during which, in addition to time, there is a loss of accuracy and an additional error is formed in calculating the Fourier coefficients and solving the problem as a whole.

It seems advisable to carry out some of these operations analytically, reducing the subsequent amount of calculations, namely. Let us carry out the

integration in the formula for calculating the Fourier coefficients (15) under the sign of the integral in (23) and under the sign of the sum in (20), and use the orthogonality of the complete system of functions

I nnx nmy I , x

{sin——sin —-— } . (26)

I 'Jx 'Jy J ,

Calculate the Fourier coefficient from the first term in (23)

£ t \ 4 f irnx . irmy

^i,nm(a) = y~r I ®i(x,y,a) sin — sin -—dxdy =

JX Jy J Jx 'Jy

n(a]

4 i' tTUX irmy

sin —-— sin —-— dxdyx

xJ

n

^x^y J ^x ly

"ma)

h(U0 -fd(xP,yp)MM,P)\

Dr-C^l (XP,

yP)\dxpdyp. (27)

By integrating on the rectangle n(a) under the sign of the integral on the rectangle n, using the representation (20), we calculate the value

4 r . nnx . irmy , ,

sin —-— sin —-— dxdy~p(M, P)\

11 1 1 ' np^s

hx ly J hx Ly

nnx nmy 2 sin —:—sin—— dxdy—— y -;-x

i il "" i — i -dxdyrr ^ k

'Jxby J 'Jx by uxuy n/ m' = i "'n'm

n(a) ,

. nn x . nm y . nn xP . nm yP , x sin —-— sin —-— sin —--sin —--. (28)

ix Ly Lx by

By performing integration under the sign of the sum of uniformly convergent series and using the orthogonality of the system (26), we obtain

4 r . nnx . irmy , ,

sin —-— sin —-— dxdy~p(M, P)\

ill i 7 wt-v * ''Pes

hx ljy J hx ljy

4 2 ^ e-k^>™>(F(*p,VP)-°) ljy . nn'xp . nm'yp

/ ï Qnn' ^sin ~ sin ■

1111 Z._/ U A nn mm; J J

'x'y 'x'y n/ m' = i n'm' 4 '-r L

2 e k^^(F(xp,yp) a) sin 7TnxP sin tvmyP (29)

Using (29), for the Fourier coefficients (27), replacing integration variables xp and yP with x and y, we get

2

$l,nm (a) = TTT X

JxJy "'nm.

X

n

h(U0 — fs(x,y))e knm(F(x,y) a)ni(x,y) sin ^^ sin

^X ly

dxdy. (30)

From the formula (30) it follows that to calculate the Fourier coefficient of the function $ on the rectangle n(a) there is no need to calculate the function itself. You can use the formula (30), which formally coincides with the formula (15) for the Fourier coefficients on the system (26) of some function depending on the Fourier indices and including information about the surface S in the form of a function F and the normal n1 calculated by the formula

n1 (x,y) = ^(Fli^^yj^^+jF^i^,^))2^!.

In this case, the number of operations has the order of 0(NXNy)2, that is, the second order in terms of the number of points, which is two orders of magnitude less than the direct calculation of the Fourier coefficients by the formulas (15), (23), (20).

Similarly, the Fourier coefficient of the second term is calculated in the formula (23)

£ i \ 4 f nnx . nmy , ,

$2,nm(a) = YY I $2(x,y,a) sin — sin -—dxdy =

n(a)

AT nnx nmy

sin —-— sin —-— dxdyx

IX ly J IX ly

X

n

fb(xP,yP)(ni,VMM,P))\Pes dxpdyp. (31)

Using the representation (20), we calculate the value

(n1 Pp(M,P))\ =

2 ^ e-(F(x p, yp)-a) nn'x

= TT ^ -k-sin^x

hxljy n/ m' = i n'm' hx

. ■nm'y nn'xP . Km'yPnn'., , X sin —-— cos —--sin —----—bx (xP, yP) +

by ix by bx

2 ^ e-( F(x p, yp)-a) nn' x

+ ^ / -;-sin—— X

J J Z_/ U J

ux°y n' m' = i n'mv 'Jx

. Km' y . nn'xp Km' yP nm',, , x sin —— sin —--cos —----—Fy (xP, yP)+

by bx by by

'

+ _ v e-k(F(^p,yp)-^) sin ^J^x I I ^ I

x y n',m' = 1 x

Km' y ku' xP Km' yP . x x sin—-—sin—--sin—--. (32)

y Jx Jy

By integrating on the rectangle n(a) under the sign of the integral on the rectangle n, performing integration under the sign of the sum of a uniformly convergent series and using the orthogonality of the system (26), we obtain

4 r . kux . Kmy , , . _ . ^

sin —-— sin —-— dxdy (n1, vPp(M,P))\

11 1 1 \ n rr \ 1 JJ\p^_s

hx ly J hx by

2 e-knm{F(xp,yp)-a)

j j b

'Jx by nnm,

Knxp . KmyP Kn ^.. ,

cos —-— sin —--—Fx (Xp, yP) +

'x 'y 'x

Knxp KmyP Km kuxp KmyP + sin —-— cos —----—by (xp, yP) + knm sin —-— sin —--

ill II

(33)

Hence and from (31) follows

= i f5(x,y)e-k™Wx,y)~°)x

'Jxynm J

n

^, , kux . Kmy , , x bx (x, y) cos —-— sin —-— dxdy+

'x 'y

+ i f5(x, y)e-knm(F(x,y)-a)F'y(x, y) sin ^ cos ^dxdy+

(¡x^y""nm J Jx '"y

n

+ JT [f 5(x,y)e-knm{F{x,y)-a) sin ^ sin ^dxdy. (34)

x y x y

n

Thus, the Fourier coefficient nm(a) is calculated as the sum of formally calculated Fourier coefficients over orthogonal systems

to f to

kux Kmy sin —:-sin

r c

I kux Kmy I

{ cos —-— sin —-— } I 'Jx Ly J

Ix ly J , y "x "y y -,

a n,m=1 a n,m=1

^ CO (35)

kux Kmy sin —:— cos

Ix lys ,

a n,m=\

of functions depending, among other things, on the indices of the Fourier coefficients. In this case, as well as when calculating the Fourier coefficient from the first term, the number of operations has the order of 0(NXNy)2, that is, the second order in terms of the number of points, which is two orders

of magnitude less than the direct calculation of the Fourier coefficients by the formulas (15), (23), (20).

Summing (30), (34), we get the Fourier coefficient

$nm(a) = $1,nm (a) + $2,nm(a). (36)

According to the remarks to the formulas (30), (34) in general, the number of operations when calculating the Fourier coefficients using these formulas relative to the number of NxNy points on the thermogram has the order of

0(NX Ny )2.

To calculate the Fourier coefficients using the formulas (30), (34), the Hamming method [10] is used.

6. Conclusion and discussion

Stable solution of the inverse problem (9) can be used for mathematical processing of thermograms taken with a thermal imager, in particular, in medicine [4], in order to correct the image on the thermogram. Note that taking into account the blood flow leads to the need to use the metaharmonic equation [11, 12] in problem (9). As already mentioned, a thermogram, with one or another reliability, convey an image of the structure of heat sources inside the body. Refinement of the image on the thermogram can be carried out within the framework of the problem (9). In this case, the function f is associated with the original thermogram, and the function uH is considered the result of processing the thermogram. Since the function uH represents the temperature distribution on a plane closer to the studied heat sources than the original surface S, we can expect a more accurate reproduction of the image of the sources on the calculated thermogram uH. The results of calculations carried out on a model example show the effectiveness of the proposed method and algorithm based on the formulas (24), (25), (23), (36), which can be used to process thermographic images.

References

[1] E. F. J. Ring, "Progress in the measurement of human body temperature," IEEE Engineering in Medicine and Biology Magazine, vol. 17, no. 4, pp. 19-24, 1998. DOI: 10.1109/51.687959.

[2] E. Y. K. Ng and N. M. Sudarshan, "Numerical computation as a tool to aid thermographic interpretation," Journal of Medical Engineering and Technology, vol. 25, no. 2, pp. 53-60, 2001. DOI: 10.1080/ 03091900110043621.

[3] B. F. Jones and P. Plassmann, "Digital infrared thermal imaging of human skin," IEEE Eng. in Med. Biol. Mag., vol. 21, no. 6, pp. 41-48, 2002. DOI: 10.1109/memb.2002.1175137.

[4] G. R. Ivanitskii, "Thermovision in medicine [Teplovideniye v meditsine]," Vestnik RAN, vol. 76, no. 1, pp. 44-53, 2006, in Russian.

[5] A. N. Tikhonov, V. B. Glasko, O. K. Litvinenko, and V. R. Melihov, "On the continuation of the potential towards disturbing masses based on the regularization method [O prodolzhenii potentsiala v storonu vozmushchayushchih mass na osnove metoda regulyarizatsii]," Izvestiya AN SSSR. Fizika Zemli, no. 1, pp. 30-48, 1968, in Russian.

[6] A. N. Tikhonov and V. J. Arsenin, Methods for solving ill-posed problems [Metody resheniya nekorrektnyh zadach]. Moscow: Nauka, 1979, in Russian.

[7] E. B. Laneev, "Construction of a Carleman function based on the Tikhonov regularization method in an ill-posed problem for the Laplace equation," Differential Equations, vol. 54, no. 4, pp. 476-485, 2018. DOI: 10.1134/S0012266118040055.

[8] E. B. Laneev, N. Y. Chernikova, and O. Baaj, "Application of the minimum principle of a Tikhonov smoothing functional in the problem of processing thermographic data," Advances in Systems Science and Applications, vol. 1, pp. 139-149, 2021. DOI: 10.25728/assa.2021.21. 1.1055.

[9] E. B. Laneev, M. N. Mouratov, and E. P. Zhidkov, "Discretization and its proof for numerical solution of a Cauchy problem for Laplace equation with inaccurately given Cauchy conditions on an inaccurately defined arbitrary surface," Physics of Particles and Nuclei Letters, vol. 5, no. 3, pp. 164-167, 2002. DOI: 10.1134/S1547477108030059.

[10] R. W. Hamming, Numerical methods for scientists and engineers. New York: McGraw-Hill Book Company, 1962.

[11] H. Pennes, "Analysis of tissue and arterial blood temperature in the resting human forearm," J. Appl. Physiol., no. 1, pp. 93-122, 1948.

[12] J. P. Agnelli, A. A. Barrea, and C. V. Turner, "Tumor location and parameter estimation by thermography," Mathematical and Computer Modelling, vol. 53, no. 7-8, pp. 1527-1534, 2011. DOI: 10.1016/j.mcm. 2010.04.003.

For citation:

O. Baaj, On the application of the Fourier method to solve the problem of correction of thermographic images, Discrete and Continuous Models and Applied Computational Science 30 (3) (2022) 205-216. DOI: 10.22363/26584670-2022-30-3-205-216.

Information about the authors:

Baaj, Obaida — postgraduate student of Nikolskiy Mathematical Institute of Peoples' Friendship University of Russia (RUDN University) (e-mail: [email protected], phone: +7(916)6890863, ORCID: https://orcid.org/0000-0003-4813-7981)

УДК 519.6

PACS 07.05.Tp, 02.60.Pn, 02.70.Bf DOI: 10.22363/2658-4670-2022-30-3-205-216

О применении метода Фурье для решения задачи коррекции термографических изображений

Обаида Бааж

Российский университет дружбы народов, ул. Миклухо-Маклая, д. 6, Москва, 117198, Россия

Аннотация. Работа посвящена построению вычислительных алгоритмов, реализующих метод коррекции термографических изображений. Коррекция осуществляется на основе решения некоторой некорректно поставленной смешанной задачи для уравнения Лапласа в цилиндрической области прямоугольного сечения. Эта задача соответствует задаче аналитического продолжения стационарного распределения температуры как гармонической функции с поверхности исследуемого объекта в сторону источников тепла. Цилиндрическая область ограничена произвольной поверхностью и плоскостью. На произвольной поверхности измеряется (и таким образом, задано) распределение температуры, называемое термограммой и воспроизводящее изображение внутренней тепловыделяющей структуры. На этой поверхности — границе исследуемого объекта — имеет место конвективный теплообмен с внешней средой заданной температуры, который описывается законом Ньютона. Это третье краевое условие, которое в совокупности с первым краевым условием соответствует заданию условий Коши — граничным значениям искомой функции и ее нормальной производной. Задача некорректно поставлена. В статье применением метода регуляризации Тихонова получено приближённое решение поставленной задачи, устойчивое по отношению к погрешности к данным Коши, и которое может быть использовано для построения эффективных вычислительных алгоритмов. В работе рассматриваются алгоритмы, позволяющие существенно уменьшить объем вычислений.

Ключевые слова: термограмма, некорректная задача, задача Коши для уравнения Лапласа, интегральное уравнение первого рода, метод регуляризации Тихонова