Almpion Ratsakou

Anastassios Skarlatos

Christophe Reboud

Dominique Lesselier

Shape reconstruction of delamination defects using thermographic infrared signals based on an enhanced Canny approach

Keywords: infrared testing, delamination, image reconstruction, Canny algorithm, TSR, semi-analytical solution

The Canny shape-reconstruction algorithm is combined with the thermographic signal reconstruction (TSR) technique for the characterisation of delamination defects in planar multilayer structures. The adopted approach consists in fitting the raw thermographic images to a low-degree polynomial in the log-log representation of the time axis, according to TSR technique, and then apply the Canny algorithm to the reconstruction of the original signals. The proposed juxtaposition of the two treatments benefits from the inherent signal de-noising and compression of the TSR algorithm, thus enhancing the performance of the overall reconstruction process. The proposed approach has been applied using synthetic signals obtained via a 3D semi-analytical simulation code and the obtained reconstructed images have been compared with the (known) delamination geometry demonstrating satisfactory results.

Introduction

Infrared non-destructive testing (IR) techniques have received growing attention in recent years thanks to their advantages, that is, fast, high resolution and contactless control, and benefiting from the technological progress in the domain of the hardware, namely infrared cameras and data acquisition equipment. The acquired infrared images are then analysed for the detection and characterization of flaws. Subsurface irregularities will affect the heat diffusion rate leading to a thermal contrast on the surface of the homogeneous work-piece. By analysing alterations or the contrast in the thermal pattern of the material surface, one can obtain information about subsurface flaws.

The material under testing may be heated using various sources as optical excitation, electromagnetic excitation, acoustic excitation, stress excitation. The most common form of active thermography [START_REF] Maldague | Introduction to NDT by active infrared thermog-857 raphy[END_REF] for material evaluation consists in using sources as flash lamps [START_REF] Omar | A quantitative review of three flash ther-859 mography processing routines[END_REF] or lasers [START_REF] Li | Crack imaging by scanning 862 pulsed laser spot thermography[END_REF][START_REF] Burrows | Ther-865 mographic detection of surface breaking defects using a scan-866 ning laser source[END_REF] where a pulse of light instantaneously heats a surface and the resulting temperature change is observed with a thermal camera. This technique is also known as pulsed thermography (PT) [START_REF] Maldague | Advances in pulsed phase 869 thermography[END_REF] or square pulse thermography (SPT) [START_REF] Feuillet | De-871 fect detection and characterization in composite materials using 872 square pulse thermography coupled with singular value decom-873 position analysis and thermal quadrupole modeling[END_REF] and has been extensively used as inspection technology for composites and layered structures. According to the heating function, other common types of active thermography are stepped thermography (ST) [START_REF] Badghaish | Non-destructive inspection of 876 composites using step heating thermography[END_REF], lock-in thermography (LT) [START_REF] Maldague | Theory and Practice of Infrared Technology for 879 Nondestructive Testing[END_REF][START_REF] Sakagami | Applications of pulse heating thermog-881 raphy and lock-in thermography to quantitative nondestructive 882 evaluations[END_REF], pulsed phase thermography (PPT) [START_REF] Maldague | Pulse phase infrared thermography[END_REF], and frequency modulated thermography (FMT) [START_REF] Tabatabaei | Thermal-wave radar: A novel 887 subsurface imaging modality with extended depth-resolution 888 dynamic range[END_REF].

The detection of the defects can be on the same side as the heat source or on the opposite side, depending on the type of access to the sample. Typical TNDT procedure re-sults in a sequence of infrared (IR) images, obtained via an 32 IR camera, that reflects the evolution of temperature over Analysing the thermal images may be an intense and 39 not an easy task to perform. The detection and character-40 ization of the defects highly depends on the quality of the itly in the reconstructed images with a different intensity.

Morphological image processing techniques can be used to detect these abnormalities. Almond et al. here [START_REF] Almond | Defect sizing by transient thermogra-934 phy. I. An analytical treatment[END_REF] suggested the identification procedure called Full Width Half Maximum (FWHM) which is based on the measurement of the distance between the points which correspond to half of the maximum temperature signal. Thermal edge effects for crack-like defects have been calculated using the Wiener-Hopf technique. Lugin and Netzelmann [START_REF] Lugin | A defect shape reconstruction algo-937 rithm for pulsed thermography[END_REF] proposed a model-based iterative algorithm for the reconstruction of two-and three-dimensional defect shapes.

More generic image processing techniques can be used for the detection of the defects' edges. A suitable edgedetection technique which can be used for accurate recognition and shape reconstruction is the Canny algorithm [START_REF] Canny | A computational approach to edge detection[END_REF].

Canny edge detector is a multi-step technique that processes an input image and detects the edges of any objects.

Because of its great accuracy, the algorithm is used in computer vision systems, yet not only. In the NDT field the applications of the algorithm are numerous and its usage can be found in several inspection techniques. Alaknanda et al. [START_REF] Alaknanda | Flaw detection in radio-943 graphic weld images using morphological approach[END_REF] used the Canny algorithm for the flaw detection in radiographic weld images. In a recent work by Yan et al. [START_REF] Yan | Recognition and evaluation of corrosion profile via pulse-947 modulation eddy current inspection in conjunction with im-948 proved Canny algorithm[END_REF] the algorithm has been used for the detection of corrosion via pulse-modulation eddy current inspection.

In thermography, Sham et al. [START_REF] Sham | Surface crack detection by flash 951 thermography on concrete surface[END_REF] used the algorithm for the detection of cracks on concrete surfaces heated by a flash lamp. This paper proposes a new approach for the real-time shape reconstruction of delamination defects in planar multilayer structures based on a combination of the TSR method and the Canny algorithm. To assess the performance of the approach, the proposed algorithm is applied to synthetic thermograms obtained by corrupting simulation results with different levels of additive Gaussian noise.

The numerical simulation of IR procedures involves the solution of the heat conduction problem in the considered work-piece, with and without defects, in order to obtain the temperature distribution at the piece interfaces, which constitutes the measurement. The full solution to this problem can be obtained using a numerical technique like the finite elements method (FEM). Nevertheless, in practical situations it is often meaningful to renounce the detailed information of the complete numerical solution in favour of fast analytical or semi-analytical approximations, which hold the essence of the thermal flow behaviour. A powerful, well-established approach is the so-called thermal quadrupole method, where the original three dimensional problem is approximated as a multilayer one dimensional problem (by ignoring the heat flow in the lateral layers directions) and modelled as a cascade of "quadrupoles" in analogy with the electrical network theory [START_REF] Maillet | Thermal quadrupoles: soving the heat equation through inte-956 gral transforms[END_REF][START_REF] Fudym | A seminumerical approach 958 for heat diffusion in heterogeneous media: one extension of the 959 analytical quadrupole method[END_REF][START_REF] Fudym | Heat diffusion at the 962 boundary of stratified media: Homogenized temperature field 963 and thermal constriction[END_REF].

Should the lateral propagation become no longer negligible, as in the case of thermal image processing at late times considered in this work, the quadrupole approach is not valid anymore, and a more general approach properly accounting for the 3D diffusion should be followed. (eddy-current) problems [START_REF] Theodoulidis | Eddy Current Canonical 966 Problems (with applications to nondestructive evaluation[END_REF][START_REF] Theodoulidis | Interaction of an eddy-current 969 coil with a right-angled conductive wedge[END_REF][START_REF] Skarlatos | Semi-analytical calculation 972 of the low-frequency electromagnetic scattering from a 973 near-surface spherical inclusion in a conducting half-space[END_REF][START_REF] Skarlatos | Calculation of the eddy-977 current flow around a cylindrical through-hole in a finite-978 thickness plate[END_REF][START_REF] Vafeas | Semianalytical method for the identification of inclusions by aircored coil interaction in ferromagnetic media[END_REF]. Defects can 130 be modelled either directly as part of the geometry like 131 in [START_REF] Skarlatos | Semi-analytical calculation 972 of the low-frequency electromagnetic scattering from a 973 near-surface spherical inclusion in a conducting half-space[END_REF][START_REF] Skarlatos | Calculation of the eddy-977 current flow around a cylindrical through-hole in a finite-978 thickness plate[END_REF][START_REF] Vafeas | Semianalytical method for the identification of inclusions by aircored coil interaction in ferromagnetic media[END_REF], or indirectly by first applying the TREE 132 method to construct the Green's function accounting for 133 the geometry of the flawless piece and treating the defect 134 as a perturbation by solving the appropriate integral equa-135 tion [START_REF] Bowler | Eddy current probe signals due to a crack at a right-angled corner[END_REF][START_REF] Pipis | ECTsignal calculation of cracks near fastener holes using an integral equation formalism with dedicated Green's kernel[END_REF].

136

The TREE approach has already been tested for the 137 inversion of thermographic signals in order to estimate the 138 thickness of delamination defects [START_REF] Ratsakou | Model based characterisation of delamination by means of thermographic inspection[END_REF]. In this work, this 139 approach is extended to multilayer specimens and is com- Let T (i) (r, t) be the temperature distribution inside the ith layer, as shown in Fig. 1, expressed as a function of position r and time t. Since no heat sources exist inside the layer by hypothesis, T (i) (r, t) satisfies the homogeneous heat equation in this layer

T (0) T (2) T (N ) T (N +1)
∇ 2 - 1 α i ∂ ∂t T (i) (r, t) = 0, (1)
where α i = κ i /c pi ρ i stands for the diffusivity of the layer.

187

Since we are interested in the transient solution of the heat equation, and we are particularly keen to the early time response, the established (and also most efficient) approach is to consider the problem in Laplace domain. Application of the Laplace transform to (1) yields

∇ 2 - s α i T (i) (r, s) = 0, (2)
s being the Laplace variable.

188

In media with flaws the total solution for the tempera- 205

∂T (i) ∂x x=0 = ∂T (i) ∂x x=Lx = 0, (3)
206

∂T (i) ∂y y=0 = ∂T (i) ∂y y=Ly = 0, (4
T (i) x=0 = T (i) x=Lx = 0, ( 5
)
222

T (i) y=0 = T (i) y=Ly = 0. ( 6
)
223

Notice here that re-gauging the temperature field to be zero at the ambient conditions, we can apply the homogeneous Dirichlet condition at the domain boundaries. The formal solution for the temperature field in the ith layer reads as

T (i) = ∞ m=1 ∞ n=1 C (i) mn e η (i) mn (z-zi-1) w mn (x, y) + ∞ m=1 ∞ n=1 D (i) mn e -η (i) mn (z-zi) w mn (x, y) (7)
for i = 1 . . . , N , with

w mn (x, y) = sin(κ m x) sin(λ n y) , (8)
the discrete eigenvalues κ m and λ n being determined by the truncation conditions in the x and y directions (5),(6), namely,

sin(κ m L x ) = 0, κ m = mπ L x , m = 1 . . . , ∞, (9) sin
(λ n L y ) = 0, λ n = nπ L y , n = 1 . . . , ∞, (10) and η (i)
mn being calculated using the dispersion equation

η (i) mn = κ 2 m + λ 2 n + s a i . (11)
Note that the development coefficients

C (i)
mn and D (i) mn standing for the upwards and downwards evanescent modes have been normalised with respect to the values of the corresponding exponential terms at the lower z i and upper z i-1 interface, respectively. There are different ways to normalise the modes. The one chosen here assures that all the corresponding exponents are negative, resulting in values lower than one, and thus overflows are avoided.

In the air regions above and underneath the piece, only the terms vanishing to infinity must be kept, reducing the general expansion to the following expressions

T (0) = ∞ m=1 ∞ n=1
D (1) mn e -η (0) mn z w mn (x, y)

and

T (N +1) = ∞ m=1 ∞ n=1 C (N +1) mn e η (0) mn (z-z N ) w mn (x, y) (13)
for the regions above and below the piece, respectively, with η (0) mn = κ 2 m + λ 2 n + s/a 0 , a 0 standing for the air diffusivity.

The expansion coefficients in [START_REF] Badghaish | Non-destructive inspection of 876 composites using step heating thermography[END_REF] and ( 12),(13) are determined according to the standard procedure by applying the continuity of the temperature and the normal to the interface heat flux across the piece interfaces.

We treat the problem in two steps. First, the solution for the flawless medium is developed, assuming a perfect thermal contact between all layers. The effect of delamination effects is accounted for in a second step by modification of the temperature continuity at the respective interfaces.

Equation system for the flawless medium

Considering the interface between the i and the (i+1)th layers, these two relations are written as

T (i) -T (i+1) z=zi = 0 (14)
and

J (i) n -J (i+1) n z=zi = J e , (15)
respectively, with the normal heat flux in the ith layer being given by the Fourier's law

J (i) n (x, y) = -κ i ∂T (i) ∂z . (16)
J e stands for the incoming flux of the heat source, which for the class of problems which we are dealing admits a non-zero value only at the upper interface, where the thermal excitation is applied. Assuming an impinging beam of constant intensity J e and of arbitrary support S(x, y), the source term can be developed in terms of the modal functions w mn (x, y) as

J e (x, y) = ∞ m=1 ∞ n=1
C (e) mn w mn (x, y) ,

where the development coefficients are given by the integrals

C (e) mn = J e Lx 0 Ly 0 S(x, y) w mn (x, y) dx dy. ( 18
)
The above integrals admit closed form expressions for the 246 most usual shapes such as rectangular and circular.

247 Substituting ( 7),( 12),( 13) and ( 17) in ( 14),(15) and weighting with w mn we obtain

C (i) mn e -η (i) mn di + D (i) mn -C (i+1) mn -D (i+1) mn e -η (i+1) mn di+1 = 0 (19)
and

-η (i) mn κ i C (i) mn e -η (i) mn di -D (i) mn η (i+1) mn κ i+1 C (i+1) mn -D (i+1) mn e -η (i+1) mn di+1 = C (e) mn δ i,0 , (20)
where d i = z i-1 -z i and δ i,0 is the Kronecker delta. No-

248

tice that for the first and last interface, C

mn and D

(N +1) mn 249 vanish.

250

Since the development basis is infinite, we need to truncate the sums ( 7),( 12) and ( 13) in order to arrive at a solvable system of equations. We are allowed to do that since the series terms are decreasing for increasing mode numbers m and n. Thus, assuming an exact arithmetic, the error can be made arbitrarily small by judicious choice of the number of modes taken into account. Let N m and N n be this number for the modes along the x and y directions. The previous result for the continuity relations at the ith interface can be conveniently represented in a matrix form as follows

A (i) X (i) = B (i) , (21)
with

A (i) = e -η (i) di 1 -1 -e -η (i+1) di+1 -β i e -η (i) di β i β i+1 -β i+1 e -η (i+1) di+1 (22)
X (i) = C (i) D (i) C (i+1) D (i+1) T (23)
and

B (i) = 0 δ i,0 C (e) T . (24)
The bold symbols in the system matrix ( 22) are understood as diagonal matrices, namely

η (i) = diag η (i)
mn , m = 0, . . . , N m , n = 0, . . . , N n and β (i) = κ i η (i) , and the coefficient vectors such as C (e) are column vectors, i.e.,

C (i) = C (i) mn (i)
, m = 0, . . . , N m , n = 0, . . . , N n , etc.

Assembling the relations [START_REF] Gao | Unsupervised Sparse 930 Pattern Diagnostic of Defects With Inductive Thermography 931 Imaging System[END_REF] for the N + 1 interfaces of the medium, we arrive at a linear system of equations, which relates the 2(N +1)N m N n sought development coefficients, i.e., the system rank is 2(N + 1)N m N n . It should be noticed here that the inversion of ( 21) using algorithms for sparse matrices is equivalent to the explicit iterative calculation of ( 7), ( 12),( 13) or the reflection/transmission coefficients approach adopted elsewhere, both from theoretical and computational point of view. The matrix form representation adopted in this work allows a more compact notation, which is the reason for choosing it.

Introduction of delamination defects

Let us consider now the multilayer medium of the previous paragraph but this time we assume an imperfect contact between the ith and the (i + 1)th layer, which can be attributed to the existence of a thin void inclusion between the two layers. This kind of imperfect contacts between two adjacent layers is practically met in media affected by delamination defects. The effect of the defect on the heat propagation inside the medium can be taken into account by locally modifying the continuity relation between the two layers i and i + 1. Consider the general case of a thin void inclusion depicted in Fig. 2, where D is void domain and d(x, y) stands for its thickness, which in the general case is a function of the x, y position. By applying the Fourier's law across the void volume, one has

T (i) T (i+) d(x, y) T (i) κ d
J z = -κ d ∇T, (25)
which for relatively small values of d(x, z) reduces to the finite difference relation

J z (x, y) ≈ -κ d T (i) (x, y) -T (i+1) (x, y) d(x, y) . ( 26
)
By introducing the local thermal resistance as [START_REF] Maillet | Thermal quadrupoles: soving the heat equation through inte-956 gral transforms[END_REF] R

(x, y) = d(x, y) κ d , (27)
the previous relation can be written as follows:

T (i) -T (i+1) = -R(x, y) J z . (28)
Substitution of ( 28) into [START_REF] Liang | Low en-922 ergy impact damage detection in CFRP using eddy current 923 pulsed thermography[END_REF] and taking (16

) into ac- count yields Nm m =0 Nn n =0 δ nm,n m -R m n κ i η (i) m n C (i) m n e -η (i) m n di + Nm m =0 Nn n =0 δ nm,n m + R m n κ i η (i) m n D (i) m n e -η (i) m n di -C (i+1) mn -D (i+1) mn e -η (i+1) mn di+1 = 0, (29)
where R mn stands for the weighted thermal resistance function

R mn = Lx -Lx Ly -Ly R(x, y) w(x, y) dxdy. (30)
Adopting again the matrix notation presented above, the continuity relations for the ith layer become

A (i) + δA (i) X (i) = B (i) (31)
with A (i) , X (i) and B (i) being given by ( 22), ( 23) and ( 24) respectively, and

δA (i) = R -β i β i 0 0 0 0 0 0 , (32)
where this time R is a full matrix that couples all modes There are several methods reported in the literature for performing this inversion. In this work we have chosen to carry out the inversion numerically using the Stehfest's algorithm [START_REF] Stehfest | Algorithm 368. Numerical inversion of the Laplace transforms[END_REF][START_REF] Stehfest | Numerical inversion of the Laplace transforms[END_REF], which is an improved variant of Gaver's method [START_REF] Gaver | Observing stochastic process, and approximate transform inversion[END_REF]. This method is well adapted to the problem of impulse thermography studied here [START_REF] Maillet | Thermal quadrupoles: soving the heat equation through inte-956 gral transforms[END_REF].

If F (s) is the known Laplace transform of the function f (t), evaluated at s = a j /t where a j = j ln(2), then an approximate value of this function at time t can be calculated as

f (t) ∼ = ln(2) t N J=1 V j F j ln(2) t . (33)
The coefficients V j are given by the following expression for an even value of N :

V j = (-1) j+N/2 × M in(N/2,j) k=Int((j+1)/2) k N/2(2k)! (N/2 -k)!k!(k -1)!(j -k)!(2k -j)! . (34)
In this equation 'Int' designates the integer part of a real number and 'Min' the minimum of two numbers.

Applying [START_REF] Skarlatos | Calculation of the eddy-977 current flow around a cylindrical through-hole in a finite-978 thickness plate[END_REF] at a number of preselected instances t,

we then obtain the temporal profile of our observable at those points. The optimal choice of the points, itself, will depend on the thermal signal processing that we wish to apply. In the following, the TSR approach is used for the representation of the thermograms at the plate interfaces, and consequently the choice of the suitable sampling points will be based on this approach.

371

Especially, in a defect-free sample, the thermal diffusion in the lateral direction could be neglected compared to the normal direction, so one can consider the one-dimensional heat equation for the thermography. In principle, the method exploits the well-known observation that in a semi-infinite flawless sample, or in a very thick slab, the surface temperature response to instantaneous uniform heating, ideal impulsive heat flux, is described by the one-dimensional heat diffusion equation

∂T ∂t = α ∂ 2 T ∂z 2 , ( 35
)
where α = κ ρc is the thermal diffusivity of the material with the solution

T (z, t) = Q √ πt e z 2 4αt (36)
at distance z from its surface, where = √ κρc is the heat effusivity and Q is the quantity of energy absorbed by the surface. Since thermal imaging is only applicable to surface temperature, from Eq. ( 36) one has the temperature increase ∆T (thermogram) as a function of time t at the surface:

∆T (t) = Q √ πt . (37)
The one-dimensional approximation of Eq. ( 37) assumes that the lateral diffusion components more or less cancel in a defect-free sample. However, in the presence of an adiabatic subsurface boundary such as a void, a wall, a buried defect, or in the case of a multi-layer sample the incident heat flow from the sample surface is impeded, and this solution no longer applies locally. The effect of a wall is shown in the left plot in Fig. 3 where the surface temperature for a semi-infinite sample is compared with the case of the presence of an adiabatic wall at three different depths. The effect of a buried defect in the thermograms is shown in the right plot in Fig. 3 where the separation of the temperature response at the surface of a solid between a sound area and a defected area should be simple when these defects are large or very close to the surface. However, when one attempts to detect thin buried defects, the effects of infrared camera noise as well as the complexity found in many samples limit the ability to discriminate between sound and flawed areas. Considering the time evolution of the surface temperature in the logarithmic domain, where Eq. ( 37) can be written in the logarithmic scale, additional insight into the surface temperature response during the cooling period of the sample is gained:

ln (∆T ) = ln Q √ π - 1 2 ln (t) . (38)
In Eq [START_REF] Ratsakou | Model based characterisation of delamination by means of thermographic inspection[END_REF] the time dependence has been separated from the input energy and material properties. Only the offset of the response will change as the sample material and the input energy vary. In Fig. 4 the thermograms have been normalised based 406 on Eq. [START_REF] Ratsakou | Model based characterisation of delamination by means of thermographic inspection[END_REF] where the term which describes the source and 407 the material effect have been removed.

408

For a given pixel, the response given by Eq. ( 38) can be approximated, in general, by a function or set of orthogonal functions. The TSR method uses a polynomial series to fit the experimental data in log-log space:

ln (∆T ) = N n=0 a n [ln(t)] n . ( 39
)
The fitting of the log-log thermogram, for each recorded pixel (i, j), by the logarithmic polynomial replaces the full sequence of the temperature response images T (i, j, t) by the series of (n + 1) images of the polynomial coefficients a 0 (i, j), • • • , a n (i, j). Once the time evolution of each pixel has been approximated by Eq. ( 39), the original data can be reconstructed as

∆T = exp N n=0 a n [ln(t)] n . (40)
Thus, it is only necessary to store the polynomial coef- The use of TSR effectively removes the temporal noise from the recorded raw signals. However, the reduction of temporal noise does not necessarily increase the flaw detectability. This can be partially addressed by the computation of the first and second time derivatives, using Eq. ( 39), without additional noise contributions, which leads to equations

d ln (∆T ) d [ln(t)] = N n=1 na n [ln(t)] n-1 , (41)
d 2 ln (∆T ) d [ln(t)] 2 = N n=2 n(n -1)a n [ln(t)] n-2 , (42)
for the first and the second derivative, respectively. The temperature in the logarithmic scale is compared with the first and the second time derivative for a sound area and a flawed area in Fig. 6.

The denoising provided by TSR is not limited to the time signals but has been partially transmitted to the spatial signals also, this is an indirect effect of the method.

The derivatives, by definition, are much more sensitive to In its classical form, the Canny algorithm, in its first step, processes the image smoothly through Gaussian convolution and obtains the gradient image through differential operation on the image which is processed via Gaussian convolution. Consider the two-dimensional Gaussian function

G(x, y) = 1 2πσ 2 exp - x 2 + y 2 2σ 2 , ( 43
)
with mean µ = 0 and standard deviation σ. The parameter σ is the width of the Gaussian filter and directly determines the effect of filtering where larger σ provides smoother images. The filtered image is derived from the convolution

I (x, y) = G(x, y) * I (x, y), (44)
where I (x, y) is the original image matrix. Making use of 457 Gaussian function's separability, ∇G can be decomposed 458 into two one-dimension filters:

459 ∂G ∂x = kx exp - x 2 2σ 2 exp - y 2 2σ 2 , ( 45a
) 460 ∂G ∂y = ky exp - y 2 2σ 2 exp - x 2 2σ 2 . ( 45b
)
461

By convolving these equations with the image we obtain: 476 483

462 L x = ∂G ∂x * I (x, y), (46a)
K x = 1 -1 1 -1 , ( 48a
) 477 K y = 1 1 -1 -1 . ( 48b
)
K x =   1 0 -1 2 1 -2 1 0 -1   , (49a)
484

K y =   1 2 1 0 1 0 -1 -2 -1   . ( 49b
)
485

After convolving the image with the kernels, the gradientcomponent intensity of the image is derived from

M (i, j) = L 2 x (i, j) + L 2 y (i, j), (50)
and its normal vector direction at the pixel (i, j) is defined as

θ(i, j) = arctan L x (i, j) L y (i, j) . (51)
To follow, after smoothing the image using Gaussian those two classes so that their combined spread is minimal.

541

Suppose that G = [0, L -1] is the range of greyscale in image F and P i is the probability of every greyscale and the threshold value τ has splitted the image in two classes which are C 0 = [0, τ ] and

C 1 = [τ + 1, L -1]. The two classes probabilities are α 0 = τ i=0 P i and α 1 = 1 -α 0 ,
respectively. The average grey values of the two classes are

µ 0 = τ i=0 iP i α 0 = µ τ α 0 and µ 1 = L-1 i=τ +1 iP i α 1 = µ -µ τ 1 -α 0 , respectively, therein µ = L-1 i=0 iP i , µ τ = τ i=0 iP i .
The criterion function has been defined as variance be-542 tween the two classes, expressed as

543 η 2 (τ ) = α 0 (µ 0 -µ) 2 + α 1 (µ 1 -µ) 2 (52) 544 = α 0 α 1 (µ 0 -µ 1 ) 2 . ( 53
)
545

The optimal threshold value τ is given by

η 2 (τ ) = max 0≤τ ≤L η 2 (τ ). ( 54
)
Threshold τ will be used as the high threshold parameter 546 τ h . The value of the low threshold τ l , usually, is set to

547 be τ l = τ h 2 .
548

The last step of the algorithm is the connection of into the plate with their larger faces set to be parallel to the plate surface. As an excitation term, a flash lamp, set above the plate and parallel to its surface, depositing a heating power density of Q = 10 4 W/m 2 at the surface of the plate, has been modelled as a Dirac's delta function in time, whereas its spatial distribution is considered to be uniform and covers all the domain of interest, as already done before.

TSR and noise reduction

The described configuration is used in this part with the addition that the flaws are located in the middle of the plate along the z-axis. From the illustrated results the beneficial effect of the

33 time.

 33 The recorded thermal response usually is degraded 34 because of several factors. Uneven heating and variations 35 of emissivity on the observed surface are only few of the 36 factors which can cause signal-to-noise-ratio (SNR) prob-37 lems and limit the potential sensitivity of any method.

38

 38

  Semi-analytical formulations are still tractable in this case, 121 provided that the considered geometry consists of canon-122 ical pieces and making the assumption that the solution 123 is negligible (or stationary) at sufficient distance from the 124 region of interest. This approach is based on the artificial 125 truncation of the computational domain, referred to in the 126 literature as the truncated region eigenfunction expansion 127 (TREE), and has been successfully applied in electromag-128 netics for the solution of magnetostatic and low-frequency 129

  140 bined with the TSR method for providing the simulated 141 signals to the reconstruction algorithm. 142 The paper is organised as follows. The considered in-143 spection problem is posed and a semi-analytical solution 144 of the heat conduction problem using the TREE approach 145 is presented in section 2. The defect reconstruction al-146 gorithm is then discussed in section 3, where a detailed 147 presentation of the constituting bricks, namely the TSR 148 and the Canny algorithm, is provided. The application of 149 the proposed approach in the case of a two-layer medium 150 embedding a number of delamination patches between the 151 two layers is examined in section 4, and the article is con-152 cluded with a brief discussion in section 5.

153 2 .

 2 Heat propagation problem in multilayer pieces 154 with delamination defects 155 2.1. Problem position 156 Let us consider a planar piece composed of N stacked 157 sheets perfectly joined with another. The medium above 158 and underneath the considered piece is air. The piece is 159 thermally excited at its upper interface by means of an im-160 pinging heat flux (flash lamp, laser source, air flow). The 161 exact form of the excitation is irrelevant for the hereafter 162 developed mathematical analysis. Solely its geometrical 163 features, i.e., the shape and the dimensions of the flux 164 spot as well as its intensity, that is, the heat power per 165 unit of surface, are of interest here. We also assume that 166 the considered medium is affected by a delamination flaw 167 located between the ith and the (i + 1)th layer. As such, 168 it is understood a local loss of contact between the two 169 layers, which physically takes the form of a very thin cav-170 ity between the two layers, filled with air. The considered 171 configuration is depicted in Fig. 1. 172 Each layer i is characterised by its mass density ρ i , its 173 heat capacity at constant (atmospheric) pressure c pi and

Figure 1 :

 1 Figure 1: Problem configuration: a planar multilayer medium affected by a number of delamination flaws is thermally excited by an impinging heat flux.

189

  ture field can be decomposed into two terms: the solution 190 for the unperturbed medium under the excitation of the 191 considered heat source, and a perturbation term produced 192 by the interaction of the incoming heat flux with the mate-193 rial flaw. Given that we are dealing with a diffusion prob-194 lem, the latter contribution can be considered as negligible 195 at sufficiently remote distance from the flaw, in the same 196 way we treat the flaw field in eddy current problems. As a 197 consequence, the temperature field at these distances will 198 be dominated by the solution for the unperturbed medium, 199 which for an infinite planar piece depends solely on the di-200 mensions of the incoming heat flux. Hence, for extended 201 (practically infinite) fluxes, the computational domain can 202 be closed (truncated) by a Neumann boundary condition, 203 which assures zero escaping heat flux from the domain 204 boundaries (otherwise stated, an adiabatic condition).

Figure 2 :

 2 Figure 2: Void inclusion between two adjacent layers. The void is assumed homogeneous and much thiner than the embedding layers.

278together.

  The introduction of the defect, in other words, 279 comes at the cost of the partial loss of the system spar-280 sity. An interesting limiting case, though, is the one where 281 the support of the delamination defect exceeds the com-282 putational domain, i.e. R(x, y) = c, for -L x ≤ x ≤ L x 283 and -L y ≤ y ≤ L y , with c constant. In this case R be-284 comes diagonal and the sparsity of the original system is 285 restored. This is the quasi-1D case in the sense that the 286 considered geometry is invariant along x and y directions, 287 yet excited via an arbitrary flux profile. An intermediate 288 case would be also having an infinite source along one of 289 the x or y directions. In this latter case, the R matrix be-290 comes block-diagonal, i.e., it has the form of a Kronecker 291 product of a diagonal and a full submatrix.

292 2 . 5 .

 25 Inverse Laplace transform 293 For the development of the semi-analytical solution we 294 have worked so far entirely in the Laplace domain. We are 295 interested though in the transient response of a number of 296 observables such as the temperature evolution at specific 297 points of the geometry and especially on the two surfaces 298 of the medium, which is actually the only accessible exper-299 imental information of the method. In order to obtain the 300 temporal evolution of the temperature field the solution in 301 the Laplace domain should be transformed back into the time domain.

3 . 3 . 1 .

 331 Defect shape reconstruction from the thermal images In order to properly characterize subsurface defects in laminates, the signals recorded by an infrared camera should be processed by some algorithms. Data processing algorithms in TNDT are either one-dimensional or two-dimensional, with the first ones being applied to pixel-based temperature evolution in time and the second ones applied to single images. The time series of each recorded pixel is separately evaluated, usually to find any alterations of material' properties or the depth of those alterations. Single IR images are normally filtered or segmented to reduce random noise or to analyse geometrical features of the areas of interest. Since its introduction [42], the thermographic signal reconstruction (TSR) method has emerged as one of the most widely used methods for enhancement, analysis and compression of raw thermographic sequences. The technique was originally developed for pulse thermography to improve contrast results and is one of the pixel-based methods. One of the well-established edge-detection algorithms 341 is the so-called Canny algorithm [24]. The methods pro-342 cesses an input image and detects the edges of any present 343 objects (defects). Because of a better signal to noise ratio 344 and detection accuracy, the Canny operator becomes the 345 evaluation criterion of other methods. The algorithm con-346 sists of Gauss filtering, gradient calculation, non-maximum 347 suppression, double thresholding, checking the edges and 348 connecting the edges. 349 In this section, the theoretical background of the two 350 mentioned methods is given and numerical results from 351 their application for the detection of delamination-type 352 defect in planar media shown. 353 Thermographic signal reconstruction 354 The technique is based on the pixel-based evaluation 355 of the time series which represent the surface temperature 356 variation over time. The technique consists of two basic 357 steps. The first step is the fitting of the recorded time 358 series in the log-log space by a logarithmic polynomial of 359 degree n. This step provides a significant compression of 360 the raw data. In a later step, the reconstruction of the 361 temperature signals in the logarithmic domain using the 362 polynomial has to be performed, providing noise-reduced 363 copy of each pixel time series. Fitting to the thermograms 364 highly depends on the time window chosen and the poly-365 nomial degree. Thus, the choice of the time window has 366 to be made with the objective to consider only the part of 367 the thermograms influenced by the physical phenomena to 368 characterize, and the polynomial degree has to be defined 369 considering that this will be a trade-off between accuracy 370 of signal reproduction and denoising.

Figure 3 :

 3 Figure 3: Surface temperature decay curves for a single point of a steel plate. Left: Comparison of temperature time plot of a semiinfinite sample with three adiabatic samples of different thickness. Right: Comparison of temperature time plot of a sound sample with two defected samples at different depth.

For aFig. 4 .

 4 Fig.[START_REF] Burrows | Ther-865 mographic detection of surface breaking defects using a scan-866 ning laser source[END_REF]. This particular time is correlated with the thickness of the plate. In the presence of a subsurface defect in a plate, or in a semi-infinite solid, the time evolution plot of the temperature corresponding to those pixels depart from that behaviour in a particular time but in a different way

Figure 4 :

 4 Figure 4: Surface temperature decay curves for a single point of a steel plate in the logarithmic scale. Left: Comparison of temperature time plot of a semi-infinite sample with three adiabatic samples of different thickness in the logarithmic space. Right: Comparison of temperature time plot of a sound sample with two defected samples at different depth in the logarithmic space.

409

  ficients a n and reconstruct the images stack at required 410 time samples. This approach provides a significant degree 411 of data compression.412A low-degree polynomial expansion is applied usually 413 to serve as a low-pass filter. Thus, this approximation pre-414 serves the essential thermal response, while rejecting nonthermal noise contributions. The use of higher-order polynomials reproduces the original data and replicates also part of the noise that appears in the later, low-amplitude data. In Fig.5reconstructed noisy signals by using polynomial degrees n = 7 and 17, with and without the presence of a defect, are compared, and the trade-off between accuracy of signal reproduction and signal denoising is clearly exhibited.

Figure 5 :

 5 Figure 5: Regression by a logarithmic polynomial of degree n equal to 7 and to 17, for a sound area, in the top plot, and a flawed area, in the bottom plot.

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Plots of the logarithmic thermogram, 1st and 2nd logarithmic time derivatives for a thick steel plate. Top: Homogeneous sound plate. Bottom: Homogeneous plate with an embedded flaw.

3. 2 .

 2 Shape reconstruction using the Canny edge detection 455 algorithm 456

  463L y = ∂G ∂y * I (x, y), (46b)464where in a matrix representation this can be written as:465 L x = K x * I (i, j), (47a) 466 L y = K y * I (i, j), (47b)467where i and j are the coordinates of a pixel in the image. 468 This is a way to compute the gradient, but not the 469 only one. Different kernels can be used to calculate the 470 image gradient. This could separate the denoising part of 471 the algorithm from the computation of the gradient. For 472 the traditional Canny algorithm, two 2 × 2 convolution 473 operators K x and K y are deployed to calculate the image 474 gradient in the x and y directions, respectively. These 475 operators are written as:

478

  In general, 2 × 2 or 3 × 3 kernels based on finite dif-479 ferences are used for the computation of the gradient. In 480 the literature one can find many other derivative kernels 481 which provide also some smoothing in the data, with the 482 most used ones being the Sobel operator:

486

  smoothing and convolving it with derivative kernels, one 487 ends up with the gradient magnitude image M (i, j) which 488 reflects the edge intensity at the pixel (i, j), and θ(i, j), 489 which reflects the normal vector at the pixel (i, j) in the 490 image. Edges of objects can be extracted from the gradient 491 component intensity image, but they will be quite blurry. 492 In this step, the algorithm aims at thinning those edges 493 by setting the pixels around local maxima in the gradient 494 image M to 0. 495 Firstly, the direction angle is rounded to 0 • , 45

549 1 .Figure 9 :

 19 Figure 9: Sketch of the configuration in the (x, y) and (x, z)-plane, left and right respectively.

  Fig. 10 the reconstructed signals versus time are compared 624

Figure 10 :

 10 Figure 10: Reconstruction of temperature noisy signals with the TSR technique for a sound area and a flawed area. Comparison of the reconstructed signals with the reference signals and the noisy signals. Left: Time signals corresponding to a sound area of the plate. Right: Time signals corresponding to a flawed area of the plate.An additional important detrimental factor for the raw

650 4 . 3 .Figure 11 :Figure 12 :

 431112 Figure 11: Comparison of the thermal signal obtained using zero and finite round-off error for a horizontal line through the centre row of delamination. Up: Raw signal. Down: Signal after application of the TSR algorithm. The TSR seems to remove the finite precision effect yielding an SNR comparable with the one obtained for zero round-off.

Figure 13 :

 13 Figure 13: Comparison of a raw image with the reconstructed image corresponding to the same optimal frame number. Left: Raw image, I . Right: Reconstructed image, I * .

Figure 14 :

 14 Figure 14: RGB(I * , I * 1 , I * 2 ) image corresponding to the optimal frame number.A similar RGB image can be constructed from the

Figure 15 :

 15 Figure 15: Reconstructed RGB image from the projection of three monomials after applying TSR on signals with SNR = 40 dB using polynomials of 5 th degree, left, and 7 th degree, right.At this point, the image I * will be provided to the Canny algorithm for the detection and the reconstruction of defects' shape. The first step of the algorithm consists of a Gaussian smoothing of which the image will be a subject. Thus the parameter σ that defines the amount of smoothing will be provided to the Canny algorithm since the smoothing algorithm has been integrated into the Canny algorithm. In our case, we perform a slight smoothing by using σ = 5. A comparison in greyscale of the input image and the smoothed one is shown in Fig.16, where the two images have been rescaled.

Figure 16 :

 16 Figure 16: The impact of the Gaussian smoothing of the input image to the Canny algorithm shown in grey scale. Left: The noisy input image I * . Right: The smoothed image I * . A gradient magnitude image M , Fig. 17, left, of the smoothed image I * , Fig. 16, right, will be derived after its convolution with the derivative kernel, and in this case, the Sobel kernel has been used. By applying non-maximum suppression to the gradient image, all values along the line of the gradient that are not peak values of the ridge have been suppressed. This leads to the image M , Fig. 17, right, which contains one-pixel-wide edges. Due to noise in the original image, M contains pixels depicting false edges and to mitigate these spurious edges, hysteresis tracking is performed using dual thresholding.

Figure 17 :

 17 Figure 17: The impact of the non-maximum suppression on the gradient intensity image. Left: The gradient intensity image before applying non-maximum suppression, M (i, j). Right: The gradient intensity image after applying non-maximum suppression, M (i, j).

Figure 18 :

 18 Figure 18: The stages of the edge detection algorithm illustrated with different colours.

Figure 19 :

 19 Figure 19: The final binary image.

Figure 20 :

 20 Figure 20: Gaussian smoothing of the input image to the Canny algorithm shown in grey scale. Left: The noisy input image I * . Right: The smoothed image I * .

765 4 . 4 .Fig. 21 .

 4421 Fig. 21. The corresponding results of the Canny algorithm 782

Figure 21 :

 21 Figure 21: Best TSR reconstructed frame with different heat-front gradient. Left: φx = 5. Right: φx = 30.

Figure 22 :

 22 Figure 22: Results for the edge detection obtained using the frames of Fig. 21 after application of non-maximum suppression. Left: φx = 5. Right: φx = 30.

797 4 . 5 .

 45 Depth effect 798The reconstruction results are expected to be also de-799 pendent from the depth of the defects, since the primary 800 field of the heat source becomes more diffuse. This ef-801 fect is important, when the inspection is carried out from 802 the opposite side of the excitation (transmission thermog-803 raphy). In order to test the algorithm performance with 804 increasing defect depth, the simulations have been car-805 ried out for three different depths, namely z = 0.75mm, 806 z = 1.5mm and z = 2.25mm, and the reconstruction al-807 gorithm has been applied using the corresponding noise-808 corrupted images. The raw and TSR reconstructed results 809 for the optimal-contrast frame are shown in Fig. 23. The 810 SNR has been set equal to 10 dB.

Figure 23 :

 23 Figure 23: Comparison of a raw image with the reconstructed image corresponding to the same optimal frame number for three different defect depths z = 0.75mm, z = 1.5mm and z = 2.25mm. Left Column: Raw images. Right Column: Reconstructed images.

  812TSR-induced denoising becomes evident and the resulting 813 enhancement of the edge detection algorithm is obvious.

814 4 . 6 .

 46 Complex defect shapes 815 For a better judgement of the reconstruction capabili-816 ties in the case of defects with more complex shapes, the 817 proposed algorithm has been applied in the case of a cir-818 cular and an irregular shaped defect, whose geometries are 819 given in Fig. 24. In the same figure are illustrated the edges 820 forms obtained after the application of the TSR-Canny ap-821 proach with non-maximum suppression. Both defects are 822 embedded in the same depth with the one of the nomi-823 nal case, namely 1.5 mm under the upper surface and the 824 SNR has been set to 10 dB. It appears that the applied 825 approach yields overall satisfactory results for both cases.

Figure 24 :

 24 Figure 24: Defect geometries in the (x, y)-plane and the corresponding edge-detection results after application of the TSR-Canny approach. Left column: A disc of radius 1.5 mm included in a ring of internal radius 8.5 mm and external radius 10 mm. Right column: A complex shape of dimensions 20 × 18 mm.

827 5 .

 5 Discussion 828 The TSR method has been combined with the Canny 829 algorithm for the detection and the dimensioning of de-830 lamination defects in thermograms obtained using pulsed 831 thermography. In this combined approach TSR is used 832 for compression and denoising the raw signals, and the 833 thus proceeded images are passed to the Canny algorithm.

834

  It is interesting to note the partial spatial denoising of 835 the images (beside the temporal filtering) achieved by the 836 TSR, which contributes to a beneficial improvement to the 837 overall processing thus enhancing the performance of the 838 Canny algorithm. 839 TSR has been also used as an acceleration technique 840 for the simulation of the heat propagation problem, in the 841 sense that the temperature field needs to be calculated 842 at a limited number of time instances, used for the poly-843 nomial fitting when the TSR approximation is computed. 844 A powerful tool for the enhancement of transient thermo-896 graphic images, Biosybern. Biomed. Eng. 35 (1) (2015) 1-9.

  897

  )

	207
	208	with the computational domain extending from 0 to L x in
	209	the x direction and from 0 to L y in the y direction. This
	210	condition is in accordance with the physics of the solution.
	211	If, instead, the consider source has a finite (with respect
	212	to the characteristic dimensions of the problem) support,
	213	the domain can be truncated either by a Dirichlet (zero
	214	temperature) or Neumann (zero flux) condition, as it is
	215	the usual practice with the TREE method 1 .
	216	In this work we deliberately choose to work with finite
	217	support sources since this case leads to simpler mathemat-
	218	ical treatment (absence of a zero-order term). It does not
	219	lack practical significance either, since narrow lamp spots
	220	and point sources like laser beams are also of interest in
	221	real world applications. We thus impose

  doi:10.1016/j.bbe.2014.07.002. Wang, Y.-H. Wang, W. Peng, J.-Q. Ji, J.-Y. Liu, P. Xiao, Current Pulsed Thermography Through Integration of PCA 910 and ICA, IEEE Sensors J. 14 (5) (2014) 1655-1663. 911 doi:10.1109/JSEN.2014.2301168.

		912
	898
	899	[14] N. Rajic, Principal component thermography for flaw con-
	900	trast enhancement and flaw depth characterisation in com-
	901	posite structures, Compos. Struct. 58 (4) (2002) 521-528.
		doi:10.1016/S0263-8223(02)00161-7.

902

[15] F. 903 Y. Wang, Independent component analysis enhanced pulse ther-904 mography for high silicon oxygen phenolic resin (HSOPR) sheet 905 with subsurface defects, Infrared Phys. Technol. 92 (2018) 345-906 349. doi:10.1016/j.infrared.2018.06.027. 907 [16] L. Cheng, B. Gao, G. Y. Tian, W. L. Woo, G. Berthiau, 908 Impact Damage Detection and Identification Using Eddy

909

It should be underlined here the difference with the eddy-current problems, where the TREE formulation has been extensively studied. In these problems, the excitation term is always localised around the coil and hence one can use either Dirichlet or Neumann conditions in almost every case. The basic criterion in the final choice is mathematical convenience.

Combined with a fast semi-analytical solution based on the TREE approach, TSR provides a powerful tool.

846

This work is restricted to the characterisation of delam-847 ination defects in homogeneous and isotropic media with 848 the thermal excitation being applied upon one of the piece 849 interfaces. Volumetric heat production via ultrasonic of 850 eddy current excitation [START_REF] Ratsakou | Fast models dedicated to simulation of eddy current thermography[END_REF][START_REF] Siakavellas | The influence of the heating rate and thermal energy on crack detection by eddy current thermography[END_REF] are of great practical in-851 terest as well, and work is under way in order to extend 852 the herein developed approach in order to treat such kind 853 of inspection scenarios as well.


MD Julie Morel
email: julie.morel@chu-brest.fr

MD, PhD Armelle Gentric

MD, PhD Cédric Annweiler

MD Morel Julie

Lower-limb hypopallesthesia: a risk factor for falls in cognitively intact non-diabetic older adults

Keywords: Peripheral neuropathies, hypopallesthesia, fall history, older adults

Peripheral neuropathies (PN) in older adults often involve altered vibrational perception, i.e. hypopallesthesia. The main objective of this cross-sectional study was to determine whether

INTRODUCTION

Falls are a major public health issue in older adults because of their high prevalence, their consequences in terms of morbidity and mortality, and their high cost [START_REF] Pfortmueller | Reducing fall risk in the elderly: risk factors and fall prevention, a systematic review[END_REF]. It is estimated that one third of those over 65 years, and half over 80, fall at least once a year [START_REF] Soriano | Falls in the community-dwelling older adult: A review for primary-care providers[END_REF]. To prevent falls, it is crucial to identify and correct the risk factors as early as possible. However some risk factors remain unknown due to changes in the expression of diseases or in the clinical interpretation with aging (3).

Among the factors with modified expression are peripheral neuropathies (PN). PN in older adults are predominantly sensitive, marked by altered vibrational perception i.e. hypopallesthesia [START_REF] Barrell | Peripheral neuropathy[END_REF]. The proportion of idiopathic PN increases with advancing age, ranging from 9-20% before age 80 to 40% beyond [START_REF] Forestier | Peripheral neuropathy in the elderly[END_REF]. While PNs are recognized as increasing the risk of falls in diabetic patients, their possible role in falling has not yet been explored among non-diabetic older adults to our knowledge (4).

We hypothesized that age-related hypopallesthesia could be an independent risk factor for falls. The objective of this cross-sectional observational study was to determine whether lower-limb hypopallesthesia is associated with the history of falls in cognitively intact nondiabetic older adults.

METHODS

The data were extracted from the French MERE cohort (ClinicalTrials.gov number, NCT01315704), which included older community-dwellers aged 60 and over consulting at the Memory Clinic of the University Hospital of Angers, France. Sampling and data collection procedures have been described and published (6). Participants were included in the present analysis when they had no history of diabetes mellitus and no cognitive impairment (Mini-Mental Status Examination score, MMSE > 24/30), and when they had a tuning fork test during the consultation. Participants were included after having given their informed consent.

The study was conducted in accordance with the ethical standards set forth in the Helsinki Declaration (1983) and protocol was approved by the local Ethical Committee (2009/15).

An operator made vibrate a dedicated tuning fork placed successively on each tibial tuberosity in order to measure the vibratory sensitivity, with a score graded between 8 (vibratory sensation until the end of the vibration = no disorder) and 0 (no sensation = severe disorder).

Hypopallesthesia was defined as a tuning fork score ≤2 bilaterally, likely to reflect PN. History of falls assessed using a standardized questionnaire by interviewing participants and their families, was retained if at least one fall was reported during the preceding 12 months. A fall was defined as an event resulting in a person coming to rest unintentionally on the ground or another lower level, not as the result of a major intrinsic event or an overwhelming hazard.

The following covariates were used as potential confounders: age, sex, polypharmacy (i.e. A descriptive analysis of the participants' characteristics was firstly performed using effectives (n) and frequencies for qualitative variables, and means and standard deviations (SD) for quantitative variables. Comparisons were performed according to those with and without hypopallesthesia using chi² or Fisher exact test for qualitative variables, and Student t-test or Mann-Whitney U test for quantitative variables according to the normal distribution assumption, as appropriate. Multivariate models were obtained by backward selection, using a P-value>0.2 for removal from the model. Sex, age, Hb level, and vitamin D supplementation were forced into the model. P-values<0.05 were considered significant. Statistical analyzes were performed using SAS software version 9.4.

RESULTS

Among 912 participants in the cohort (mean age 72years, 45.3% female), the results of the tuning fork test were available for 400 participants, including 39 with diabetes mellitus, who were excluded. Of the 361 non-diabetic people, 157 had no cognitive impairment and were included in the analysis. Participants' characteristics are summarized in Table 1 (mean age±SD, 71.5±3.8years; 35.7% female; 19.7% fallers). Fallers exhibited more often hypopallesthesia than non-fallers (13.3% versus 1.6%, P=0.04). In contrast, performance on postural tests did not differ between groups. The multiple logistic regression showed that hypopallesthesia was associated with the history of falls (odds ratio (OR)=19.5 [95% confidence interval (CI): 2.7-143.7], P=0.004). Female sex was also associated with the history of falls (OR=5.0 [95CI: 1.9-13.0], P=0.001).

DISCUSSION

Our results show that the finding of lower-limb hypopallesthesia, suggestive of PN, was associated with the history of falls in cognitively-intact non-diabetic older adults, irrespective of all studied potential confounders.

Clinical PN is a known risk factor for disability and mobility limitations in those with diabetes mellitus (7). However, advancing age is also a risk factor for PN, regardless of diabetes and our results suggest that this age-related PN is also associated to falls. Thus our findings encourage systematic testing of lower-limb hypopallesthesia in all subjects at risk of falls, even in the absence of diabetes mellitus.

Literature has reported that somatosensory information, including both tactile and proprioceptive systems, is responsible for maintaining balance. While standing, humans continuously oscillate and neurological integrity is needed to detect these oscillations to ensure efficient postural control. In the case of NP, more oscillations are needed, hence an increased fall risk (8).

Our study has some limitations. First, our study is observational with a cross-sectional design, which prevents any causal inference. Second, the study was monocentric with a limited number of participants who may be unrepresentative of all older adults. Third, falls were self-reported, and up to 50% of self-reported falls may be forgotten. Fourth, balance disorders were assessed with clinical tests that are less accurate than spatial temporal gait measures.

In conclusion, we found a positive independent and possibly early association between lowerlimb hypopallesthesia and history of falls in cognitively-intact non-diabetic older adults in whom postural tests were still unaltered. This result should be confirmed in larger and preferentially prospective studies, in order to determine with a higher level of evidence the potential role of PN as a risk factor for falls in older adults. If these results are confirmed, the measurement of hypopallestesia could be systematized in ambulatory medicine to set up early fall prevention programs.

Contributors

Julie Morel contributed to study concept and design, acquisition of data, analysis and interpretation of data, and drafting of the manuscript, and has full access to the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analyses.

Armelle Gentric contributed to study concept and design, analysis and interpretation of data, drafting of the manuscript, critical revision of the manuscript for important intellectual content, and study supervision.

Cédric Annweiler contributed to study concept and design, acquisition of data, and critical revision of the manuscript for important intellectual content.

Highlights

• The main objective of this cross-sectional study was to determine whether age-related lower-limb hypopallesthesia is associated with a history of falls in cognitively intact non-diabetic older adults.

• This association was found in the sample of 157 people.

• Systematization of the measurement of hypopallestesia would be to the benefit of early fall prevention programmes.

  using >4 different drugs per day), usual vitamin D supplementation, MMSE score, depression (i.e. 4-item Geriatric Depression Scale [GDS] score≥1), anemia defined as hemoglobin (Hb) level<13g/dL for men and 12g/L for women, and performance on postural tests (one-leg balance test [OLB; normal>5seconds] and five time sit-to-stand test [FTSS; normal<15seconds]).

  6. G.T. Duval, A. Brangier, J. Barré, C.P. Launay, O. Beauchet, C. Annweiler C. Vitamin D deficiency and incident onset of orthostatic hypotension in older adults: preliminary results from the 'MERE' Study. J. Am. Geriatr. Soc. 63 (2015) 1245-1247. 7. P. Hewston, N. Deshpande. Falls and balance impairments in older adults with type 2 diabetes: thinking beyond diabetic peripheral neuropathy. Can. J. Diabetes. 40 (2016) 6-9. 8. L.A. Lipsitz, B. Manor, D. Habtemariam, I. Iloputaife, J. Zhou, T.G. Travison. The pace and prognosis of peripheral sensory loss in advanced age: association with gait speed and falls. BMC Geriatr. 18 (2018): 274.

Table 1 . Characteristics of participants (n=157).

 1

	Whole sample	No fall history	Fall history	P-value
	n = 157	n = 126	n = 31

ACKNOWLEDGMENTS

We are grateful to the participants for their cooperation.

We also thank Maelys Consigny from Clinical Research and Investigation Center, University Hospital, Brest, France, for daily assistance.

Funding

No sources of funding were used to assist in the preparation of this article.

Ethical approval

The study was conducted in accordance with the ethical standards set forth in the Helsinki Declaration (1983). Ethics approval was obtained from the Ethics Board of the University Hospital of Angers, France (2009/15). Informed consent was obtained at enrolment according to the protocol approved by the local institutional review board.

Provenance and peer review

This article was not commissioned. Peer review was directed by Teresa Liu-Ambrose independently of Cédric Annweiler, an author and Maturitas editor, who was blinded to the process.

Research data (data sharing and collaboration)

Conflict of interest

The authors declare that they have no conflict of interest.