M Marmonier
email: maxance.marmonier@cea.fr

S Robert

J Laurent

C Prada

Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing

Keywords:

Real-time 3D ultrasound imaging with matrix arrays remains a challenge in Non-Destructive Testing (NDT) due to the time-consuming reconstruction algorithms based on delay-and-sum operations. Other algorithms operating in the Fourier domain have lower algorithmic complexities and therefore higher frame rates at the cost of more storage space, which may limit the number of reconstruction points. In this paper, we present an implementation for real-time 3D imaging of the Total Focusing Method (TFM) and the Plane Wave Imaging (PWI), as well as of their Fourier-domain counterparts, referred to as k-TFM and k-PWI. For both acquisition types, the Fourier-domain algorithms are used to increase frame rates, and they are compared to the time-domain TFM and PWI in terms of image quality, frame rates and memory requirements. In order to greatly reduce their memory requirements, a new implementation of k-TFM and k-PWI is proposed. The four imaging methods are then evaluated by imaging in real time a block of stainless steel containing a 3D network of spherical porosities produced by additive layer manufacturing using a powder bed laser fusion process.

Introduction

Nowadays, real-time 2D imaging methods with linear or matrix arrays are available in most industrial Non-Destructive Testing (NDT) systems [START_REF] Sutcliffe | Real-time full matrix capture for ultrasonic non-destructive testing with acceleration of post-processing through graphic hardware[END_REF][START_REF] Lewandowski | A GPU-Based Portable Phased-Array System with Full-Matrix Capture[END_REF]. The most popular is the Total Focusing Method (TFM), which requires the inter-element response acquisition often referred to as Full Matrix Capture (FMC) [START_REF] Drinkwater | Ultrasonic arrays for non-destructive evaluation: A review[END_REF]. The set of FMC signals are focused numerically both in transmit and receive mode everywhere in the region of interest, which greatly improves the spatial resolution compared to focused B-scan or S-scan [START_REF] Holmes | Post-processing of the full matrix of ultrasonic transmit-receive array data for nondestructive evaluation[END_REF]. 2D matrix arrays offer the great advantage of allowing 3D imaging without any displacement [START_REF] Lane | 3D ultrasonic inspection of anisotropic aerospace components[END_REF][START_REF] Kitazawa | Three-dimensional visualisation and evaluation techniques for volumetrically scanned data of ultrasonic phased arrays[END_REF][START_REF] Brizuela | Improving elevation resolution in phased-array inspections for NDT[END_REF]. However, the 3D TFM remains a challenge for real-time applications due to time-consuming reconstructions. Indeed, the number of signals to be processed is huge as it corresponds to the square of the number of array elements, which is typically equal to 256 in NDT. In addition, the number of points to be reconstructed is generally two orders of magnitude greater than in 2D imaging. The problem is therefore to process this large amount of data in a short period of time. To this end, the Plane Wave Imaging (PWI) method, initially developed for elastography in medical imaging [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF], has been adapted to various NDT inspection configurations, such as multi-mode imaging of crack-like defects [START_REF] Le Jeune | Plane Wave Imaging for ultrasonic non-destructive testing: Generalization to multimodal imaging[END_REF][START_REF] Merabet | The multi-mode plane wave imaging in the Fourier domain: Theory and applications to fast ultrasound imaging of cracks[END_REF] or imaging under complex surfaces [START_REF] Le Jeune | Plane wave imaging for ultrasonic inspection of irregular structures with high frame rates[END_REF][START_REF] Rachev | Plane Wave Imaging Techniques for Immersion Testing of Components With Nonplanar Surfaces[END_REF]. In most cases, the transmission of plane waves reduces the amount of data to be processed while maintaining a similar image quality compared to the TFM. Today the PWI method is already available in some NDT acquisition systems to speed up imaging in 2D [START_REF] Ducousso | Real-time imaging of microcracks on metallic surface using total focusing method and plane wave imaging with Rayleigh waves[END_REF], but the principle of plane wave emissions re-mains still insufficient to achieve 3D real-time imaging with current NDT systems.

The TFM and PWI reconstructions take place in the time-domain and are based on the Delay-And-Sum (DAS) principle, where the signals are delayed for each divergent or plane wave emission and then coherently summed. Another way to compute images is to use Fourier-domain algorithms that lead to a shorter computation times than the DAS algorithms. The two best-known Fourier-domain imaging methods are the Stolt's migration based on the exploding reflector model [START_REF] Stolt | Migration by Fourier transform[END_REF][START_REF] Garcia | Stolt's f-k migration for plane wave ultrasound imaging[END_REF] and the wavenumber algorithm. The second is more suitable in many NDT configurations because it remains valid when images are wider than the array aperture [START_REF] Merabet | 2-D and 3-D Reconstruction Algorithms in the Fourier Domain for Plane-Wave Imaging in Nondestructive Testing[END_REF] or when images of crack type defects involve multi-mode ultrasonic paths [START_REF] Merabet | The multi-mode plane wave imaging in the Fourier domain: Theory and applications to fast ultrasound imaging of cracks[END_REF]. The wavenumber algorithm consists in computing the spectrum of the set of signals on a specific wavenumber grid that depends on the chosen reconstruction region, and then in obtaining the image by inverse Fourier transform. It was first used in Synthetic Aperture Radar (SAR) imaging with a monostatic configuration [START_REF] Franceschetti | A SAR processor based on twodimensional FFT codes[END_REF] and later with a bistatic one [START_REF] Callow | Wavenumber domain reconstruction of SAR/SAS imagery using single transmitter and multiple-receiver geometry[END_REF]. In the field of NDT, Hunter et al. [START_REF] Hunter | The wavenumber algorithm for full-matrix imaging using an ultrasonic array[END_REF] propose to extend the method to FMC data in order to reduce computation times in 2D imaging compared to TFM. Almost at the same time, Cheng and Lu [START_REF] Cheng | Extended high-frame rate imaging method with limited-diffraction beams[END_REF] proposed a similar algorithm dedicated to 2D imaging with plane wave emissions. In NDT, the performances of the Fourier-domain algorithms are generally evaluated in 2D cases [START_REF] Zhuang | Comparison of Time Domain and Frequency-Wavenumber Domain Ultrasonic Array Imaging Algorithms for Non-Destructive Evaluation[END_REF] where real-time frame rates can be achieved with approaches that are more usual. Here, we propose to study the efficiency of these algorithms for 3D applications where real-time imaging remains a difficult task whatever the method.

In this paper, we present an implementation of both TFM and PWI in 3D as well as their frequency-domain counterparts, referred to as k-TFM and k-PWI in the following.

Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing

In both cases, the aforementioned frequency-domain algorithms are used in order to increase frame rates. They are compared to the time-domain TFM and PWI in terms of image quality, frame rates and memory requirements. For the k-TFM, the wavenumber algorithm is generalized to 3D imaging with matrix arrays, while for the k-PWI, we follow the same approach as Merabet et al. [START_REF] Merabet | 2-D and 3-D Reconstruction Algorithms in the Fourier Domain for Plane-Wave Imaging in Nondestructive Testing[END_REF] based on the theory developed by Cheng and Lu [START_REF] Cheng | Extended high-frame rate imaging method with limited-diffraction beams[END_REF]. The four imaging algorithms are implemented in a versatile Vantage system allowing to perform acquisitions and reconstructions directly with Matlab (Mathworks). The Matlab parallel computing toolbox is used in order to launch all computations on a graphic card. Finally, the four real-time imaging algorithms are evaluated on an additively manufactured stainless steel sample containing a 3D network of equispaced spherical porosities.

The paper is organized as follows. In Section 2, the theoretical background of 3D imaging is presented. We first recall the mathematical formulation of the time-domain TFM and PWI based on the DAS algorithms, and of their frequency-domain counterparts k-TFM and k-PWI. The relationship between the image and the signal spectra is derived for both types of acquisition. In Section 3, some technical details and tricks to improve real-time imaging are provided.

In particular, we propose an implementation of the Fourierdomain algorithms that greatly reduces the required memory space. Minimum storage capacity and algorithmic complexities for k-TFM and k-PWI are derived and compared to those obtained with TFM and PWI. In Section 4, we present the experimental setup and the results obtained from the realtime imaging of the steel sample. The experimental results are discussed and compared with theoretical and numerical predictions.

Theoretical background in 3D imaging

In this section, the 3D imaging algorithms are overviewed. For sake of clarity, we only present the theory for a square matrix array and an inspection configuration where elements are in direct contact with the sample. However, given that inspections through a water layer are quite common in NDT, in Appendix A we derive the method accounting for the refraction at the fluid-solid interface.

Time-domain methods: TFM and PWI

For the FMC, the array elements are excited one by one to transmit divergent waves, where is the total number of elements. As signals are recorded for each transmission, the dataset to handle is composed of 2 samples with the number of time samples. As indicated in Fig. 1(a), = ( 1 , 2 , 0) and = ( 1 , 2 , 0) are the vectors that locate respectively the transmitting and receiving element centers. The origin of the Cartesian coordinate system is placed at the center of the matrix and the sets of transmitter and receiver positions are denoted  and  .

The TFM consists in focusing in transmission and reception at every point in the region of interest. If ( , , ) is the signal received by the element at when the element at is excited by an impulse signal, the TFM image is built as the coherent sum of 2 analytic signals ̂ ( , , ) = ( , , )+ { ( , , )} with appropriate time delays, where  denotes the Hilbert transform. Using these notations, the image amplitude at a given voxel located by the position vector for a Dirac excitation is defined as

TFM ( ) = | | | | | ∑ ∈ ∑ ∈ ̂ , , ‖ -‖ + ‖ -‖ | | | | | , (1)
where is the phase velocity of longitudinal or transverse ultrasonics waves depending on the NDT application.

For the PWI method, we assume that the solid is insonified by plane waves propagating in different directions. As shown in Fig. 1(b), each plane wave is indexed by where ∈  = [1, ... , ] and have a wave vector equal to where = sin cos , sin sin , cos ,

with ( , ) the pair of angles defining the direction of propagation of the -th plane wave and the wavenumber. If all elements are used as receivers, then the number of samples to be stored after transmission is . Similarly to the TFM, if ( , ) is the signal received by the element at for the -th emission, the PWI image results from a coherent sum of analytic signals ̂ ( , ) = ( , ) + { ( , )}. The image amplitude at a point is defined as

PWI ( ) = | | | | | ∑ ∈ ∑ ∈ ̂ , . + ‖ -‖ + | | | | | , (3)
where is the emission delay of one element taken as reference.

In practice, the signals and are known only for times in a set of recorded time samples  . Thus, the acquisition gate has to be carefully chosen to allow the image computation at all points in the region of interest. Furthermore, in Eqs. ( 2) and (3), ̂ and ̂ must be interpolated to accurately extract the amplitudes corresponding to the theoretical arrival times. A nearest neighbour interpolation is a sufficiently good approximation if the sampling frequency is large with respect to the central frequency of the array.

Frequency-domain methods: k-TFM and k-PWI

The general idea of the Fourier-domain imaging methods is to use the mathematical relation between the Fourier transform of the dataset and the Fourier transform of the reflector distribution in the medium. Here, the theory of the Fourierdomain imaging is presented in 3D: first k-TFM where spherical waves insonify the medium, then k-PWI with plane wave transmissions. 3D k-TFM is an extension of the work of Hunter et al. [START_REF] Hunter | The wavenumber algorithm for full-matrix imaging using an ultrasonic array[END_REF] to matrix arrays, while 3D k-PWI was recently developed by Merabet et al. [START_REF] Merabet | 2-D and 3-D Reconstruction Algorithms in the Fourier Domain for Plane-Wave Imaging in Nondestructive Testing[END_REF].

3D k-TFM method

Let us consider a spherical wave transmitted by an element located at , backscattered by a distribution of isotropic point-like reflectors, and received by an element located at . Assuming that there is no attenuation and that the elements have omnidirectional radiation patterns, the time Fourier transform of may be written as

( , , ) = 1 (4 ) 2 ∫ e -i ‖ -‖ ‖ -‖ ( ) e -i ‖ -‖ ‖ -‖ d , ( 4
)
where is the distribution of point-like reflectors and d = d d d . For sake of brevity, the factor in front of the integral in Eq. ( 4) is omitted in the following. We denote by ̂ the 5-D Fourier transform of , = ( , , ) , = ( 1 , 2 , 0) and = ( 1 , 2 , 0) the wave vectors corresponding to , and , respectively. Following the same approach than in 2D [START_REF] Hunter | The wavenumber algorithm for full-matrix imaging using an ultrasonic array[END_REF], we use the Weyl identity to decomposed each Green function in Eq. (4) as continuous sum of plane waves [START_REF] Chew | Waves and Fields in Inhomogeneous Media[END_REF], and the 3D Fourier transform of can be expressed as

̂ ( ) = ̂ ( , , ), (5)
provided that , , and satisfy the relation

= + (6)
where, noting = or , = ( 1 , 2 , ) and

= √ 2 -2 1 -2 2 .
To compute an image, we need to express , and as a function of . As the problem is ill-posed with three equations for five unknowns, we solve it by fixing in Eq. ( 6). The solution is expressed as

⎧ ⎪ ⎨ ⎪ ⎩ 1 = -1 2 = -2 = 1 2 √ 2 + ‖ ‖ 2 -‖ ‖ 2 2 + 4 2 ‖ ‖ 2 . ( 7
)
Thus, we can compute a spectrum ̂ from the spectrum ̂ at a fixed value of that will be denoted ̂ in the following. Then, the spectra are added for all , and the final image is obtained by the 3D inverse Fourier Transform over dimensions 1 , 2 and k, noted  -1

, . Finally, the imaging equation is

k-TFM ( ) = | | | | | |  -1 , ∑ ∈  ̂ ( , ) | | | | | | , (8)
where Eq. ( 7) provides the required frequencies of the spectra ̂ and where the set   is the wavenumber counterpart of .

3D k-PWI method

The imaging equation for k-PWI is derived in the same way as previously. First, we express the time Fourier transform of as

( , ) = 1 4 ∫ e -i . ( ) e -i ‖ -‖ ‖ -‖ d , (9)
Here, the Weyl identity needs to be used only once to replace the Green function in Eq. ( 9), which leads to

̂ ( ) = ̂ ( , ), (10)
with ̂ the 3D Fourier transform of and where , and must satisfy the relation

= + . ( 11
)
From Eq. ( 11), and can be expressed as a function of . We obtain

⎧ ⎪ ⎨ ⎪ ⎩ 1 = -sin cos 2 = -sin sin = ‖ ‖ 2 2 . . ( 12
)
Thus, a spectrum ̂ is calculated for each incident plane wave. The spectra corresponding to plane waves are added and the final image is obtained by 3D inverse Fourier transform. The imaging equation is

k-PWI ( ) = | | | | | | -1 , ∑ ∈ ̂ ( , ) | | | | | | , ( 13
)
where Eq. ( 12) provides the required frequencies of the spectra ̂ .

Implementation of imaging algorithms

In this section, the implementation of the time-domain and Fourier-domain methods is presented in detail. First, the strategy adopted to compute the signal spectra is described.

In particular, we propose a method to compute the signal spectra without resorting to any interpolation along the two lateral spatial dimensions. Then, the imaging algorithms in the time and frequency domains are compared in terms of required memory space and complexities.

Computation of the signal spectra

The frequencies ( , ) required for k-TFM and k-PWI depend on the desired image frequencies and are defined by Eqs. ( 7) and [START_REF] Rachev | Plane Wave Imaging Techniques for Immersion Testing of Components With Nonplanar Surfaces[END_REF]. However, in practice, the spectra ̂ and ̂ are computed with Discrete Fourier Transforms (DFTs) on the wavenumber grid   ×   given by

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩   = -1 2 + 1, ..., 1 2 2 1 × -2 2 + 1, ..., 2 2 2 2 × {0}   = -2 + 1, ..., 2 2 , ( 14
)
where is the pitch of the matrix array and to the number of time samples . The problem is to obtain the spectra on the required wavenumber grid instead, and a classical way to do this is to interpolate the initial spectra known on   ×   . In this paper, only the dimension will be treated in this way, while a method is proposed to compute the signal spectra without any interpolation along the lateral dimensions 1 and 2 .

In this method, the shifted signal spectra ̂ ′ and ̂ ′ are used instead of the original signal spectra ̂ and ̂ . These spectra are calculated applying a phase shift before the spatial Fourier transform, which is written as

̂ ′ ( , ) =  ( , )e i . ( 15
)
and

̂ ′ ( , ) =  ( , )e i ′ . , ( 16
)
where  is the Fourier transform along the dimensions 1 and 2 . If ̂ and ̂ are replaced respectively by ̂ ′ and ̂ ′ , the expressions of the required wavenumbers 1 and 2 in Eqs. ( 7) and ( 12) simplify and become 1 = and 2 = . Satisfying these two relations means having equality between the first two dimensions of   and the desired lateral grid   ×   , where

⎧ ⎪ ⎨ ⎪ ⎩   = -2 + 1, ..., 2 2   = -2 + 1, ..., 2 2 . ( 17
)
In other words, it is necessary to have equal sampling steps as well as equal number of points for   ×   and the first two dimensions of   . To this end, the method comprises two steps. The first is zero-padded Fourier transform along the 1 and 2 dimensions of the above shifted signal spectra ′ and ′ . The number of points 1 and 2 are chosen so that the wavenumber pitch is identical for the dimensions 1 and as well as for the dimensions 2 and . For a volume with lateral sizes and , the wavenumber sampling step along the dimensions and must be equal to 2 ∕ and 2 ∕ . Furthermore, the pitch of the matrix probe along the  and  axes imposes that the computed spectra vary between -∕ and ∕ . Thus, we choose

1 = (18)
and

2 = . ( 19
)
The second step is to extend the obtained discrete spectra. For a given spectrum known between -∕ and ∕ , an extension by periodicity over any interval is possible, as illustrated in Fig. 2 for the dimension 1 . The lower and upper bounds of the wavenumber grid are chosen equal to -∕ and ∕ for the dimension 1 and equal to -∕ and ∕ for the dimension 2 . In this way, the amplitudes of the spectra ̂ ′ and ̂ ′ are directly obtained on the desired wavenumber grid   ×   , without any interpolation. Note that we do not directly apply an inverse DFT to this extended spectrum, it is only used for the interpolation step. This makes it possible to obtain an image spectrum after interpolation having a wider support than if zeros were added instead of extending the signal spectra, which improves the quality of the reconstructed images.

The expression of in Eqs. ( 7) and ( 12) is non-linear with respect to and so, the required wavenumber grid for is non-uniform. Thus, the method described above is not applicable for this dimension and, as previously mentioned, an interpolation of the signal spectra along the dimension is performed. A nearest neighbour interpolation is chosen but, knowing that the interpolation accuracy is a preponderant factor for the quality of the Fourier-domain reconstructions, zero-padding is used for the time dimension before DFT computation and is chosen large enough to obtain a satisfactory image quality. Note that different ways to control the quality of images obtained with Fourier methods are proposed in the literature. For example, by improving the quality of the interpolation step [START_REF] Harlan | Avoiding interpolation artifacts in Stolt migration[END_REF], by using a non-uniform FFT [START_REF] Kruizinga | Plane-wave ultrasound beamforming using a nonuniform fast fourier transform[END_REF][START_REF] Li | A Fast Radial Scanned Near-Field 3-D SAR Imaging System and the Reconstruction Method[END_REF] or by using the interpolation-free Stolt mapping method developed by Li et al. [START_REF] Li | Interpolation-Free Stolt Mapping for SAR Imaging[END_REF]. However, the amount of memory required for these algorithms will be greater and the obtained frame rate lower. For real-time 3D imaging, the use of zero-padding and nearest neighbour interpolation seems to be a good compromise to obtain good performance without excessively degrading the image quality.

Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing

Required memory space

Regarding the real-time implementation of imaging algorithms, high frame rates can be reached if all matrices needed to reconstruct an image is computed before starting the imaging loop. However, large enough memory space is required to store the reconstruction algorithm matrices, in addition to the recorded signals. This is the main constraint for real-time imaging. In this paper, we do not make a tradeoff between frame rates and memory requirements and we only present the fastest possible implementation and therefore the most memory intensive.

For TFM and PWI, the computation of × and × arrival times is needed for each reconstruction point in the region of interest. For each of these times, the index of the closest value in  is stored in the graphic card. In the same way, for the Fourier-domain algorithms, the indices corresponding to the closest values in   of the required values of are stored in the memory space of the graphic card. Using shifted signal spectra reduces the number of indexes to be stored by a factor of for the k-TFM and k-PWI methods because the required in Eqs. [START_REF] Brizuela | Improving elevation resolution in phased-array inspections for NDT[END_REF] and [START_REF] Rachev | Plane Wave Imaging Techniques for Immersion Testing of Components With Nonplanar Surfaces[END_REF] no longer depend on 1 and 2 which are imposed equal to and in this case. The skew factors and , as well as the phase shifts in Eqs. ( 15) and ( 16) and those used for taking into account a water/solid interface (see Appendix A), are also computed and stored.

The sizes of the different stored matrices are summarized in Table 1 and the memory requirements corresponding to our experimental settings are presented in Fig. 3(a). It appears that our implementation of the Fourier-domain algorithms provides lower memory costs compared to those in the time domain. This is especially valuable in the context of 3D real-time imaging where the amount of available memory in NDT systems is a huge constraint if large volumes have to be reconstructed with high resolutions. In our experiment, where = 256, = 33, = 237, = 711, = 9.9 mm, = 9.9 mm and = = ∕2, k-TFM requires approximately 80 times less memory space than TFM for 10 5 < < 10 8 where = is the total number of voxels. In the case of plane wave emission and for the same range of values, k-PWI has memory requirement between 2 and 60 times lower than PWI.

Algorithmic complexity analysis

For the TFM and PWI methods, the signal amplitudes in the dataset are extracted using the indices stored in the graphic card. This operation is a nearest neighbour interpolation of the signals that comprises only memory accesses. The results of the interpolation step is stored in a 5-D double array of size × × × × for the TFM method and × × × × for the PWI method. Each component of those 5-D variable corresponds to an amplitude extracted in the experimental dataset at a corresponding arrival time. For both methods, the calculation of a single voxel requires a double summation (see Eqs. ( 1) and ( 3)). Thus, for voxels, the algorithmic complexities for TFM and PWI are respectively

 TFM = 2 (20)
and

 PWI = . ( 21
)
For the k-TFM algorithm, the 5-D Fourier transform of the dataset must be computed before interpolation. As explained before, zero-padding is used for the DFT computations along the dimensions 1 , 2 and with the number of points and . Thus, this 5-D Fourier transform comprises the following computation steps:

1. Zero-padded DFT along the time dimension with suppression of the useless negative frequencies; 2. DFTs along the dimensions 1 and 2 ; 3. Application of phase shifts before zero-padded DFTs along the dimensions 1 and 2 .

Direct and inverse DFTs are computed using the Fast Fourier Transform (FFT) and the algorithmic complexity of the 5-D FFT is therefore equal to

1 2 log 2 √ 1 2
. The application of the phase shift corresponds to a term-by-term multiplication of two 5-D arrays of size 2 × √ × √ × √ × √ which has therefore an algorithmic complexity equal to Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing with complexity equal to log 2 . All these above operations lead to an algorithmic complexity for k-TFM equal to

 k-TFM = 1 2 log 2 √ 1 2 + 2 2 + 2 + log 2 . ( 22
)
k-PWI includes the same computation steps, which leads to a similar algorithmic complexity :

 k-PWI = 1 2 log 2 √ 1 2 + 2 + 2 + log 2 . ( 23
)
In this expression, the first term corresponds to the calculation of 3D FFT, one for each transmission. It is the main difference with the algorithmic complexity expression of k-TFM where the first term corresponds to a single 5-D FFT of the entire recorded data set. The other differences are simply due to the different numbers of transmissions for both methods. Each algorithmic complexity consists in two mains contributions: the first corresponds to the Fourier transform of the dataset and mainly depends on the number of array elements; and the second corresponds to the inverse Fourier transform of the image spectrum and mainly depends on the number of voxels . The algorithmic complexities of the four imaging methods were estimated considering the matrix probe and the imaging parameters used in our experiment, and are plotted as a function of the number of voxels in Fig. 3(b). Using the Fourier-domain algorithms leads to lower algorithmic complexities, but only for a number of reconstructed voxels greater than 5.10 3 in our configuration. As the desired number of voxels in 3D imaging typically varies between 10 5 and 10 8 , the Fourier-domain algorithms are therefore useful to speed up real-time 3D imaging. However, if a matrix array with more elements is used, the imaging algorithms will be more efficient in the time domain rather than in the frequency domain for a larger range of number of voxels. Indeed, the terms that mainly depend of the number of elements in the algorithmic complexities of the Fourier-domain methods will stay preponderant for a larger number of reconstruction points.

An analysis of ( 22) and [START_REF] Harlan | Avoiding interpolation artifacts in Stolt migration[END_REF] shows that the gains in terms of numbers of operations for the Fourier-domain algorithms Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing for 10 5 < < 10 8 are given by

 k-TFM  TFM ∼ 2 + log 2 2 (24)
and

 k-PWI  PWI ∼ 2 + log 2 . ( 25
)
For example, these expressions lead to approximatively 125 times less operations for k-TFM over TFM and of approximatively 95 for k-PWI over PWI. Although this does not directly give us the gain in terms of improved frame rate, it indicates that Fourier domain algorithms are a good tool for performing real-time imaging.

Real-time 3D imaging results

In the present section, an additively manufactured 316L stainless steel block containing a set of spherical porosities is imaged with the four imaging methods. First, the experimental setup and the manufacturing additive sample are presented. Then, the obtained results are analyzed and the performance of the methods are compared.

Experimental setup

The experimental setup is shown in Fig. 4. The matrix array manufactured and distributed by Imasonic (Voray-surl'Ognon, France) is composed of 16 × 16 square elements of 1 mm 2 arranged with a pitch of 1.1 mm in the  and  directions. The elements operate at the central frequency of 5 MHz. The sample to be imaged is a 316L stainless steel block elaborated by additive manufacturing with a laser powder-bed fusion process. The block includes a network of 27 porosities of 0.8 diameter organized in a cubic mesh of 4 mm period. The block is 20-mm high and is assumed to be isotropic with longitudinal wave velocity equal to 5.74 mm. s -1 [START_REF] Kim | Ultrasonic nondestructive evaluation of laser powder bed fusion 316L stainless steel[END_REF]. The experiment was carried out in a water tank and the matrix array placed 10 mm above the sample surface so that its central axis coincides with the central porosities of the 3D network. It should be pointed out that the printed sample was controlled by X-ray radiography to ensure that all the defects are identical spheres and that they are not completely filled with steel powder.

The data acquisition with divergent or plane waves were performed using a 256-channel Vantage (Verasonics) system and the 3D imaging codes were implemented on Matlab 2021a and were launch on a graphic card (Nvidia Geforce RTX 3090) with 24 GB of available memory. The recorded time window is set between 16.25 and 20.05 µs, so the echoes from the specimen surfaces are not digitized. This avoids imaging artifacts possibly caused by aliasing for the Fourier-domain methods [START_REF] Hunter | The wavenumber algorithm for full-matrix imaging using an ultrasonic array[END_REF]. The Vantage system has a programmable sample rate of up to 62.5 MHz. The maximum available sample rate is used, so each recorded signal contains 237 time samples.

For the PWI and k-PWI methods, a satisfactory quality image of the printed sample can be obtained with only 33 plane waves provided that the transmission angles , are defined as

, = (0, 0) (2.5 , 45 ) where ∈ [1, ... ,4] ∈ [0, ... ,7] . ( 26
)
It should be noted that our objective is not to obtain the best possible image by increasing the number of plane waves, but to evaluate the benefit of a reconstruction in the Fourier domain for a given number of plane waves. Note also that no apodization was applied at the time of acquisition in order to maintain a satisfactory spatial resolution.

Comparison of obtained image quality

The 3D images corresponding to TFM, k-TFM, PWI and k-PWI are displayed in Fig. 5 in the form of isosurfaces. The region of interest is a 10 mm 3 cubic volume centered on the porosity of coordinates (0, 0, 22.1) located at the center of the 3D network. The number of voxels is × × = 24 × 24 × 50 which corresponds to voxels of sizes close to ∕3 in the two lateral directions and close to ∕6 in the axial direction. For each volume, different isovalues are used depending on the depth of the porosities so that all the echoes remain visible with similar sizes.

A first observation is that all porosities are properly imaged. The different used isovalues used suggest that the reconstructed porosity amplitudes decrease faster as a function of the distance from the array for both Fourier-domain reconstructions. For a quantitative comparison, the amplitudes along the vertical central axis of the matrix are plotted in Fig. 6. The amplitude gaps between the top porosity at  = 20 mm and the bottom porosity at  = 28 mm are 7.9 dB for TFM and 13 dB for k-TFM. For PWI and k-PWI, the amplitude gaps correspond to 5.4 and 11.5 dB, respectively. As expected, the amplitudes decrease more rapidly for TFM and k-TFM compared to PWI and k-PWI due to the attenuation by spherical spreading.

The horizontal cross sections are displayed in Fig. 7 for depths  = 18.1 mm ,  = 22.1 mm and  = 26.1 mm where the image amplitude is maximum. They correspond to the top of the spherical defects. For the Fourierdomain methods, we observe that the spatial resolution of the upper plane is better and the image quality degrades more rapidly with depth than the time-domain methods. This is particularly clear for k-TFM. In order to compare the quality of the reconstructions, two metrics are used: the Signalto-Noise Ratio (SNR) and the Array Performance Indicator (API) [START_REF] Mondal | Design and Evaluation of Two Dimensional Phased Array Ultrasonic Transducers[END_REF]. The SNR measures the detectability of a defect in the presence of noise and can be expressed as

SNR = echo noise , ( 27
)
where echo is the maximum amplitude of a defect echo and noise is the standard deviation of the noise at the same depth. In our experiment, noise was measured by moving the probe over a defect-free region. The API is introduced to quantify the spatial resolution of an imaging method and can be expressed as

API = -6dB 3 , (28)
where , and are the dimensions of the voxels along the ,  and  axes, and -6dB is the number of voxels with amplitudes greater than the -6 dB threshold where 0 dB corresponds to the maximum reconstructed amplitude of the defect. The SNR and API values measured around the positions of the central porosities are presented in Table 2. The Fourier-domain imaging results in lower SNR values compared to the time-domain reconstructions, and this SNR decay with depth is similar for the TFM and PWI methods. For the upper porosity, Fourier-domain and time-domain imaging provide similar API. However, the API degrades more rapidly with depth for the Fourier-domain imaging.

Following the approach taken by Velichko et al. [START_REF] Velichko | An analytical comparison of ultrasonic array imaging algorithms[END_REF] and Merabet et al. [START_REF] Merabet | 2-D and 3-D Reconstruction Algorithms in the Fourier Domain for Plane-Wave Imaging in Nondestructive Testing[END_REF], the differences observed experimentally between the 3D images computed in the time and frequency domains can be justified theoretically using asymptotic expansions of the integral expressions of the reflector distributions. The analytical demonstration leading to the imaging equations is summarized in Appendix B. In the simplified case where the array elements are directly in contact with the sample, the imaging equations for TFM and k-TFM can be put into the following integral forms:

TFM ( ) = ∭ TFM ( , , , ) d d d (29)
and

k-TFM ( ) ∼ →∞ 2 ∭ 4 cos 2 ( , ) cos 2 ( , ) ‖ -‖‖ -‖ TFM ( , , , ) d d d , ( 30
)
where TFM ( , , , ) = ( , , )e i (‖ -‖+‖ -‖) , is a complex constant and, noting = or , ( , ) is the angle at which a point located at is seen from the element positioned at . For the PWI method, the expressions obtained for a single plane wave emission are where PWI ( , , ) = ( , )e i ( . +‖ -‖) . For the frequency-domain methods k-TFM and k-PWI, the terms TFM ( , , , ) and PWI ( , , ) are multiplied by a function that acts as a space-time filter. The effects of this filters can be summarized in three points:

PWI ( ) = ∬ PWI ( , , ) d d (31)
1. The factors ‖ -‖ -1 theoretically confirm that the echo amplitudes decrease faster with the distance to the array when images are computed in the frequency domain; 2. The factors ( , ) show that the Fourier-domain imaging penalizes the detection of defects located outside of the matrix aperture; 3. The factors in the form of powers of shift the spectral content towards the higher frequencies, which inproves the resolution of images in the Fourier domain but can also increases the noise sensitivity.

For the first point, the experimental observations are in perfect agreement with the theory since the porosity echoes in the k-TFM and k-PWI images decrease more rapidly with depth compared to what is observed in the TFM and PWI Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing Depth (mm) TFM k-TFM PWI k-PWI images. As the porosities are all located below the matrix aperture, we cannot make a comparison in relation to the second point. Concerning the third point, our previous study in 3D imaging with simulated echoes confirmed that a better spatial resolution is obtained for the reconstructed images in the Fourier domain in the absence of noise [START_REF] Merabet | 2-D and 3-D Reconstruction Algorithms in the Fourier Domain for Plane-Wave Imaging in Nondestructive Testing[END_REF]. Looking at the API values, it seems that the greater noise sensitivity of these methods does not allow to find the same results with the experimental data except for the highest central porosity where the k-PWI obtains a better API than the PWI. Furthermore, in the presence of experimental noise, the SNR obtained with the Fourier domain algorithms is lower than for the time domain.

Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing

Comparison of obtained frame rates

The frame rates obtained with the four real-time imaging methods are given in Table 3 for voxel numbers ranging from 3.5 × 10 3 to 1.3 × 10 7 . Due to the limited available memory in the graphic card, some algorithms cannot be applied to the larger number of voxels. Our implementation of the algorithms in the Fourier domain allows us to gain a factor 46 on the maximum number of voxels for the TFM and a factor 59 for the PWI. Furthermore, as expected from the analysis of algorithmic complexities, better frame rates are obtained with the Fourier-domain imaging beyond a certain number of voxels, approximatively 10 4 voxels for the k-TFM method and 10 5 voxels for the k-PWI method. For an image containing 24 × 24 × 48 voxels, the obtained frame rate is approximately 4 times greater with k-TFM compared to TFM, and for an image formed of 48 × 48 × 96 voxels the frame rate is approximately 11 times higher with k-PWI than with PWI. The larger the number of voxels, the more useful the Fourier-domain algorithms. Thus, a low memory cost implementation of these algorithms makes it possible to overcome the difficulties of 3D real-time imaging.

Conclusion

In this paper, we have presented a way to perform fast and memory-efficient 3D imaging with matrix arrays. Two types of acquisition have been investigated: the first consists in the transmission of spherical waves by exciting the elements one by one, and the second uses incident plane waves where the elements are all excited in parallel with delay laws. In each case, we have presented the theory of the imaging method in the Fourier domain and demonstrated a significant decrease of the algorithmic complexity compared to the equivalent method operating in the time domain. In addition, the Fourier-domain methods have been improved in terms of required memory space, which makes it possible to compute 3D images with many more voxels in comparison with the implementation usually adopted in the literature.

The four methods, namely TFM, k-TFM, PWI and k-PWI, have been evaluated experimentally using a 256channel Vantage (Verasonics) system and thanks to parallel computations on graphic card. The solid sample is a stainless steel block elaborated by additive manufacturing and containing a 3D network of porosities of 0.8 mm diameter, and 3D images are produced with a 16 × 16 matrix array at 5 MHz. In terms of required memory space, it has been shown that for an image composed of 10 8 voxels, the k-TFM algorithm requires 80 times less memory than TFM while k-PWI 60 times less than PWI.

The gain in frame rate depends strongly on the number of voxels. For an image with 2.2×10 5 voxels, it has been shown that the frame rate is 10 times higher for k-PWI compared to PWI. For k-TFM, the frame rate is 4 times higher compared to TFM for an image with 2.8 × 10 4 voxels. These results show that the use of imaging algorithms in the Fourier domain is a relevant approach for real-time 3D imaging. In particular, the performances obtained with the k-PWI method in terms of frame rate and maximum number of voxels greatly exceed those of the industrial imaging systems used nowadays for on-site inspections.

From a qualitative point of view, the experimental comparison of the 3D images formed with time and frequencydomain algorithms confirmed the theory: the faster amplitude decay with depth and the greater sensitivity to noise of the Fourier-domain algorithms lead to lower values of SNR compared to time-domain methods, in particular for defects far from the array. However in case of electronic noise, it Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing Table 3 Frame rates expressed in frames-per-second reached with the dierent imaging algorithms as a function of the number of voxels. A cross indicates that reconstruction is not feasible due to lack of memory in the graphic card. A star indicates that the associated reconstruction is performed with the maximum possible number of voxels.

× ×

Total voxels TFM k-TFM PWI k-PWI should be noted that since the frame rates are higher, the images obtained with k-TFM and k-PWI could be improved by averaging several reconstructed images by choosing a compromise between frame rate and image quality. Provided acquisition time is small compared to image computation time, is should be also possible to increase signal to noise ratio by averaging the acquired signals. This is the case when the number of emissions is reduced by using plane wave emissions, for example, and when the region of interest is close to the probe, which allows for a short signal digitisation time.

Reducing noise in the data in this way specifically improves the quality of images reconstructed with k-TFM and k-PWI being more noise sensitive methods than TFM and PWI. Further theoretical and experimental work is needed to better understand the origin of this attenuation and sensitivity to noise, which is specific to Fourier domain algorithms. Future work will focus on finding methods to correct these effects using specific spatio-temporal filters for example. For the Fourier-domain algorithms, we use an idea developed by Skjelvareid et al. [START_REF] Skjelvareid | Synthetic aperture focusing of ultrasonic data from multilayered media using an omega-K algorithm[END_REF]. The spectra which can be obtained with a control in direct contact are extrapolated from the spectra obtained in immersion. This allows to apply the Fourier-domain algorithms presented in this paper without any modification. The idea is to apply a phase shift to the spectra of recorded signals corresponding to the inverse propagation in water. For k-TFM, we multiply ̂ by the phase factor exp iℎ √ ̃ -+ √ ̃ -, where ̃ = ∕ ̃ is the wavenumber corresponding to the propagation in water. The term exp iℎ √ ̃corresponds to the first propagation in water, just after

Figure 1 :

 1 Figure 1: Geometries and notations used for 3D imaging with divergent or plane waves: (a) transmission of a spherical wave from a single element centered at = ( 1 , 2 , 0) ; (b) transmission of a plane wavefront of normal vector by time delayed excitations of all array elements. For both transmission types, all the receivers of coordinates = ( 1 , 2 , 0) are used to record signals.

1 , 2

 12 and are respectively the lengths of the Fourier transforms along the dimensions 1 , 2 and . Provided that DFTs are computed without zero-padding,

1

 1

1 Figure 2 :

 12 Figure 2: Diagram of spectral extension principle for an image of length along the  axis: (a) initial spectrum ; (b) spectrum after extension. The trapezoidal area in gray represents the support of propagative waves.

2 2 .

 2 After the extension of the spectrum along the dimensions 1 and 2 , a nearest neighbour interpolation is performed along the dimension . We thus obtain a 5-D array of size × × × √ × √ that we multiply by and sum over the two last dimensions. Each of these two operations has a complexity equal to and finally, an image spectrum is stored in a 3D array of size × × . The last operation is an inverse 3D FFT

Figure 3 :

 3 Figure 3: Required memories (a) and algorithmic complexities (b) plotted as a function of the number of voxels for TFM, PWI, k-TFM and k-PWI. The results correspond to the acquisition and reconstruction parameters used in the experiments: = 256, = 33, = 237, = 711,

Figure 4 :

 4 Figure 4: Experimental setup with a matrix array and a stainless steel block immersed in water (a). The block was elaborated by additive manufacturing and contains 27 porosities of 0.8 mm diameter arranged in a 3D network of 4 mm period: side view with a water column heigth of 10 mm (b); 3D view (c).

  PWI ( , , , ) d d ,[START_REF] Hörmander | The Analysis of Linear Partial Differential Operators I: Distribution Theory and Fourier Analysis[END_REF]

Figure 5 :

 5 Figure 5: Experimental 3D images of the network of 27 spherical porosities shown in the form of isosurfaces: (a) TFM; (b) k-TFM; (c) PWI; (d) k-PWI. The 0 dB reference level corresponds to the maximum over amplitude measured in the whole volume, and the isovalues chosen to display the three horizontal networks of 9 porosities are indicated in dB in the table below. (a) TFM. (b) k-TFM. (c) PWI. (d) k-PWI.

Figure 6 :

 6 Figure 6: Amplitudes of the 3D images extracted along the vertical central axis of the matrix array: (a) TFM and k-TFM; (b) PWI and k-PWI.

Z = 18

 18

Figure 7 :

 7 Figure 7: Horizontal cross sections of reconstructed volumes at  = 18.1 mm,  = 22.1 mm and  = 26.1 mm. (a) TFM. (b) k-TFM. (c) PWI. (d) k-PWI. Level limits are adapted to noise and therefore TFM images are displayed between 0 and -30 dB and PWI ones between 0 and -20 dB.

  point on the surface of the incident plane wave and̃ = sin ̃ cos , sin ̃ sin , cos ̃ with ̃ = arcsin ̃ sinis the normal vector of the plane wavefront in water.

Table 1

 1 Types and sizes of the variables precomputed and stored in the graphic card.

	Variable Interpolation indices	Type uint32	TFM	2	k-TFM	Size	PWI	k-PWI
	Skew factors Signals phase shift single complex single complex Water phase shift single complex			2	∕2		∕2 ∕2

Table 2

 2 SNR and API values measured around the positions of the central porosities at  = 18.1 mm,  = 22.1 mm and  = 26.1 mm.

	Depth (mm) 18.1 22.1 26.1	SNR (dB) TFM k-TFM PWI k-PWI TFM k-TFM PWI k-PWI API 39 33 16 11 0.43 0.53 0.81 0.7 41 33 15 8 0.48 0.86 1.2 1.34 35 28 18 9 1.2 1.9 1.3 2.1

Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing

Acknowledgement

The authors would like to thank the team in charge of the Additive Manufacturing Hub (AFH) platform hosted at CEA for the elaboration of the stainless steel sample with spherical porosities.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A Imaging equation for a water-steel interface

The theoretical background presented in the paper considers a matrix array directly in contact with the specimen. It is common to carry out an "immersion" control where the probe is placed in water at a certain distance above the sample. The water height, noted ℎ, is taken into account differently depending on the type of algorithm used.

For the time-domain algorithms, the calculation of the delay corresponding to the direct path of the wave is adapted. For TFM, the image amplitude can be expressed as

where ̃ is the phase velocity in water and the vectors and locate the impact points on the surface for the incident path and outgoing path, respectively [START_REF] Le Jeune | Plane Wave Imaging for ultrasonic non-destructive testing: Generalization to multimodal imaging[END_REF]. The vectors and are computed by minimising the time of flight while respecting Fermat's principle. Likewise, the image amplitude for PWI is expressed as

Real-time 3D imaging with Fourier-domain algorithms and matrix arrays applied to non-destructive testing emission, and the term exp iℎ √ ̃ -corresponds to the second propagation in water of form the surface to the receivers. For k-PWI, we multiply ̂ by the phase factor exp iℎ ̃ cos ̂ + √ ̃ -where the term exp iℎ ̃ cos ̂ corresponds to the propagation in water of the plane wave. These operations provide a good approximation of the recorded data spectra that we would have obtained with a contact control.

Appendix B Derivation of the reflector distributions

In the following, the expressions of the distribution of point-like reflectors are derived for the four imaging methods following the idea proposed in [START_REF] Velichko | An analytical comparison of ultrasonic array imaging algorithms[END_REF] for the divergent waves emissions and adapted in [START_REF] Merabet | 2-D and 3-D Reconstruction Algorithms in the Fourier Domain for Plane-Wave Imaging in Nondestructive Testing[END_REF] to plane waves emissions. The objective is to obtain expressions in the form of integral over the same variables for both types of acquisition in order to compare the theoretical image amplitudes.

The isotropic point-like reflector distribution obtained with TFM at the point is expressed as

where TFM ( , , , ) = ( , , )e i (‖ -‖+‖ -‖) . Using Eq. ( 5), the reflector distribution at the point obtained for k-TFM can be expressed as To obtain an expression comparable to Eq. ( 35), the integral ( , , ) = ∫ e i( .( -)+ ) d (37) must be calculated. An asymptotic expression of the integrals of the form

can be derived for → ∞ using the stationary phase method [START_REF] Bender | Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory[END_REF][START_REF] Hörmander | The Analysis of Linear Partial Differential Operators I: Distribution Theory and Fourier Analysis[END_REF]. The result is given by

where * 1 , * 2 is a stationary point of , i.e., satisfying ∇ * 1 , * 2 = 0, Q is the Hessian matrix of and

. Thus, by setting = ,

1 , 2 = and

an asymptotic expression of integral in Eq. ( 37) can be obtained for → ∞, which is given by

with a complex constant. The integrals over and in Eq. ( 36) can be evaluated for → ∞ using the result given in Eq. ( 41). The corresponding asymptotic expression is given by ,,,) d d d (42) with ( , ) = arccos ‖ -‖ the angle at which the point is seen from the element positioned at . This expression is valid for the points sufficiently far away from the array, i.e. for ≫ .

The reflector distribution obtained with PWI for one plane wave emitted in the direction at the point is expressed as PWI ( ) = ∬ PWI ( , , ) d d .

(

where PWI ( , , ) = ( , )e i ( . +‖ -‖) . Using Eq. ( 10), the reflector distribution obtained for k-PWI at the point is expressed as


Sogo Amagai
email: samagai@s.h.k.u-tokyo.ac.jp

Pierre Maret
email: pierre.maret@univ-st-etienne.fr

Kamal Singh
email: kamal.singh@univ-st-etienne.fr

Improving maintainability of calendar-based IoT-driven office with knowledge graph and rule-based reasoning

Keywords: Internet of Things, Knowledge graph, rule-based reasoning

Designing smart Internet of Things (IoT)-driven offices has recently gained considerable attention. Such smart offices are expected to enhance worker productivity and economise electricity for sustainability. However, there are several open issues related to the maintenance of IoT applications.

In this study, we implement a calendar-based IoT-driven office exploiting knowledge graphs and rule-based reasoning. The smart IoT-driven office consists of micro-controllers (MCUs) with sensors and actuators, and a controller running on PC. These smart nodes send sensor data which in turn is represented as a Knowledge graph. This has different advantages in terms of interoperability, ease of integration with other systems and improved device compatibility. Decision-making is carried out by a rule-based reasoner. The reasoning process of a rule-based reasoner is more transparent than some machine learning approaches. This allows for easy configuration and adaptation of rules.

CCS CONCEPTS

• Information systems → Graph-based database models.

INTRODUCTION

The Internet of Things (IoT) is a network of physical devices which are embedded with sensors, software, and network connectivity, allowing them to collect and exchange data. In recent years, different solutions for smart offices based on IoT have been proposed [START_REF] Horch | Why Should Only Your Home Be Smart? -A Vision for the Office of Tomorrow[END_REF]. An IoT-driven office offers several benefits, including improved employee productivity [START_REF] Hossain | Toward the Prediction of Environmental Thermal Comfort Sensation Using Wearables[END_REF] and reduced energy consumption [START_REF] Röcker | Services and Applications for Smart Office Environments -A Survey of State-of-the-Art Usage Scenarios[END_REF]. However, it is well known that IoT has several challenges to overcome [START_REF] Sobin | A Survey on Architecture, Protocols and Challenges in IoT[END_REF]. Some of the challenges are interoperability, integration and maintenance which needs to be overcome for enabling automation using IoT in large systems such as Smart Cities [START_REF] Röcker | Services and Applications for Smart Office Environments -A Survey of State-of-the-Art Usage Scenarios[END_REF]. In this paper, we argue that the maintainability of IoT applications can be improved through the use of Knowledge graph (KG). This will bring benefits such as interoperability, ease of changing hardware devices, ease of integrating with other systems, and expressivity of rule based reasoning. IoT systems with higher maintainability are expected to better address the multiple challenges faced by the IoT domain.

This work proposes a calendar-based IoT-driven office that utilizes KG as well as rule-based reasoning. A calendar-based IoT is suitable for office workers to easily manage conditions of many rooms. However, as the number of rooms increase, the amount of data will increase and managing data will become complicated. Thus, a KG which represents knowledge in a graph format [START_REF] Hogan | Knowledge Graphs[END_REF] is more suitable for such scenarios. KG can connect different types of data in meaningful and flexible ways. The decision-making process is then carried out by a rule-based reasoner that operates on explicit rules and the knowledge base, such as the KG [START_REF] Chen | A review: Knowledge reasoning over knowledge graph[END_REF][START_REF] Dutta | Integrating case-and rule-based reasoning[END_REF], as illustrated in Figure 1.

Through the use of KG and rule-based reasoners, a prototype of a calendar-based IoT-driven office with higher maintainability has been developed. Our work proposes easily maintainable IoT systems that can be generalised to other IoT-based systems. This system serves as a basis for the development of IoT-based systems on a larger scale. The problems identified during the actual development of this system will be valuable for future IoT development.

This paper is organized as follows. Section 1 provides the background and overview of this work. Section 2 outlines how we implemented our proposed system. Section 3 presents the tasks

MAINTAINABILITY OF IoT APPLICATIONS 2.1 Graph-based Data Model

Graph-based data models, such as KGs, are a popular way of modeling and analyzing data that contains complex relationships between entities. This type of data model represents entities as nodes and their relationships as edges, creating a network of interconnected nodes that can be queried and analyzed in a variety of ways [START_REF] Hogan | Knowledge Graphs[END_REF].

Figure 2 is a simplified illustration of the KG used in this system. This KG is based on Semantic Sensor Network (SSN) Ontology [START_REF] Haller | The modular SSN ontology: A joint W3C and OGC standard specifying the semantics of sensors, observations, sampling, and actuation[END_REF], Sensor, Observation, Sample, and Actuator Ontology (SOSA) [START_REF] Janowicz | SOSA: A lightweight ontology for sensors, observations, samples, and actuators[END_REF] and Time Ontology [START_REF] Cox | Time ontology in OWL[END_REF]. The entity called Room is associated with an entity of physical quantities observed by sensors and information extracted from several APIs. The room entity also has entities for objects in the room, such as windows and actuators. In addition, room's automation is kept off when not needed in order to reduce the power consumption of the room. This on or off state is determined by calendar event information. KG can meaningfully represent different types of data and knowledge, such as sensor output and calendar event information. The room entities hold information about the state of the room derived from the reasoning results. Figure 3 gives snippets of Resource Description Framework (RDF) data to create a knowledge graph. The first few lines are used for the triple that describes the windows, actuators, etc. in the room. In the following lines, another triple describes the information about the room. Information such as the state of the room and the time is stored. The entity named as Observation stores information related to temperatures observed by sensors. The entity calendarEvent holds the event's start time, end time, and event title retrieved from Google Calendar API. Google Calendar can also be easily used by other systems. That would be helpful during the integration with other systems.

Such knowledge graphs in general can easily integrate with other data models and are a powerful and flexible way to represent knowledge. Indeed, this integration is expected not only for smart offices but also for large-scale systems such as smart cities. Another advantage is high compatibility with hardware devices for the maintenance of IoT systems: by representing device capabilities and relationships in a standardised way, knowledge graphs can facilitate interoperability between different devices and systems, making it easier to integrate new devices into IoT applications [START_REF] Janowicz | SOSA: A lightweight ontology for sensors, observations, samples, and actuators[END_REF]. Consequently, knowledge graph data models will be one of the necessary elements of IoT applications. This will become even more prevalent in the future as IoT applications become larger.

Rule-based reasoning

Rule-based reasoners are a type of reasoning engine which infer new knowledge from existing knowledge using a set of logical rules. We used Apache Jena as the reasoning engine. These reasoning engines work by inputting a set of rules and a set of explicit facts, which are RDF triples in our case, and generating new implicit facts derived by using the rules. The main advantage of rule-based reasoners is their flexibility. Rules can be defined in many different ways, allowing for a wide range of reasoning tasks. In addition, rules can be updated or modified as new information becomes available, allowing continuous improvement and refinement of the reasoning Table 1: General rule expressions setting the levels of temperature, air quality, humidity, and energy consumption respectively Table 2: Rule expressions related to alerts Table 3: Rule expressions related to calendar events process. Figure 4 shows the example rules used in our application. These rules are written in Notation 3 syntax. These rules are written manually by the developers. Here we describe four rules among several. The first rule is to turn on room automation if the calendar event comes within 30 minutes. Thus, the event start time is used in the rule to compare it with the current time. The lessThan and greaterThan operators used for comparison are builtin primitives from Apache Jena. The second rule is to determine if the room temperature is higher than the threshold value. The third rule is to open the room window when the room needs ventilation. The fourth rule is to operate the actuator when the room window is closed and the room is hot.

Table 1 shows the meaning of these rule expressions used in this application. These rules relate to physical quantities such as temperature and humidity in the room, as observed by the sensors. These rules set a threshold and add new information, such as alerts, to the room entities in the knowledge graph when they are above or below that threshold. The threshold can be changed by the users.

Table 2 shows the rule expressions that are based on alerts or alarms. The actual actions to be done by the actuators and the material to be described in the user interface are obtained from these rule expressions via room alarms derived from the general rule expressions. Furthermore, regulations are defined not just for a single alarm but also for many alarms and window conditions. This method of expressing rules allows for workplace comfort as well as energy savings.

Today, many workers use virtual calendars such as Google Calendar. Knowledge of room occupancy information from the calendar allows more sophisticated inferences to be made to improve workspace conditions. Table 3 shows the rule expressions for reasoning based on calendar events. As shown in the first row of the table, the information from the calendar can be used to condition the room in advance. A pre-conditioned room makes it easier for workers to start their work. As shown in the second row of the table, energy savings can also be expected by switching off the automation to make the room comfortable ahead of the leaving time. Conversely, there is the problem that these functions cause the room automation to be switched off regardless of whether people are in the room. We used PIR sensors for presence detection and wrote rules as shown in the fourth row of the table. With KG, it's easy to utilize sensor data and calendar information simultaneously, enabling sophisticated reasoning.

Implementation

Figure 5 depicts a detailed system diagram of this application. The main circuit comprises several sensors surrounding the ESP32 that measure temperature, humidity, CO 2 , human presence, and power consumption. The measured values are sent to the PC via Message Queuing Telemetry Transport (MQTT). MQTT communication is bidirectional, and it has many advantages, such as being lightweight and battery-friendly, which makes it widely used in machine-tomachine systems. Furthermore, information on calendar events and weather forecasts is sent to the PC via HyperText Transfer Protocol (HTTP). The knowledge graph is constructed from this information, and the rule base is developed by considering use cases and scenarios. The rule base could include more rules by considering more detailed situations. Based on the rule base and the knowledge graph, reasoning is performed to derive the workspace conditions and based on that the actuators can take actions. On the PC, the code for MQTT communication, API usage, knowledge graph composition, and rule engine execution is written in Java due to the availability of the library. The new information derived by reasoning is then sent to a LED display, smartphone, and actuator. The LED display shows basic information, such as the room's temperature and humidity, as well as tasks that the user needs to perform. For many people, room conditions may not be interesting enough to be constantly monitored. Hence, a to-do list is used to encourage users to take actions regularly, and it is less likely to be forgotten than other methods. It also does not reduce productivity by avoiding high number of notifications that are usually sent by smartphones and laptops. The user interface, which can be viewed on smartphones and laptops, was developed using the Node-RED visual programming language. Electromagnetic relays are used to switch devices on and off for room automation.

DISCUSSION

As explained so far, we have developed a calendar-based IoT-driven office that is expected to be actively expanded and integrated into larger IoT applications in the future. Through this development, IoT applications with knowledge graphs do not need to redesign the data model when adding new sensors. In other words, the addition of new sensors, etc., is easier as compared to IoT applications which lack knowledge graphs. Additionally, the feature of knowledge graphs, which can be easily integrated with other data models, can be expected to be integrated with larger-scale systems.

The use of a rule-based reasoner makes it easy to improve and add new rules. Furthermore, the open-box nature of the system allows developers and others to understand the reasoning process. On the other hand, a disadvantage of this system is the developer effort required to write rules for possible scenarios. There are quite a lot of possible scenarios in a room, and it can be quite difficult and labor-intensive to write rules that cover all of them. Another potential major problem is rule conflict. It is possible that several logically correct but incompatible rules could be created, causing the system to behave unexpectedly.

Furthermore, in the calendar-based IoT, using KG and rule-based reasoner can be suitable because it can operate a lot of data in a human-readable way and can easily integrate with other IoT systems.

CONCLUSION

In conclusion, we have developed a calendar-based IoT-driven office system with a knowledge graph and a rule-based reasoner to improve the maintainability of IoT systems. The use of knowledge graphs in this system has several advantages, such as high device compatibility and an open-box rule base. Thus, this study has contributed to improving the maintainability of IoT systems that can perform sophisticated reasoning by combining sensor and calendar information.

However, it has become clear that there are still more challenges that need to be addressed to further improve maintainability. Future challenges include the following.

• Generalisation of the knowledge graph • Simplifying rule editing and compatibility checking • Experimental validation and analysis on larger systems • Running the reasoning system itself on the embedded devices The IoT's high maintainability is likely to lead to further adoption, and we intend to address these challenges in the future.

Figure 1 :

 1 Figure 1: Overview of rule-based reasoner for IoT

Figure 2 :

 2 Figure 2: knowledge graph for this application

Figure 3 :

 3 Figure 3: Example of RDF data

Figure 4 :

 4 Figure 4: Snippets of rules

Figure 5 :

 5 Figure 5: System configuration

ACKNOWLEDGMENTS

This work was partially supported by grant ANR-19-CE23-0012 from the Agence Nationale de la Recherche, France, for the CoSWoT project.