Two gesture-computing approaches by using electromagnetic waves

We are concerned with a novel sensor-based gesture input/instruction technology which enables human beings to interact with computers conveniently. The human being wears an emitter on the finger or holds a digital pen that generates a time harmonic point charge. The inputs/instructions are performed through moving the finger or the digital pen. The computer recognizes the instruction by determining the motion trajectory of the dynamic point charge from the collected electromagnetic field measurement data. The identification process is mathematically modelled as a dynamic inverse source problem for time-dependent Maxwell's equations. From a practical point of view, the point source should be assumed to move in an unknown inhomogeneous background medium, which models the human body and the surroundings. Moreover, a salient feature is that the electromagnetic radiated data are only collected in a limited aperture. For the inverse problem, we develop, from the respectively deterministic and stochastic viewpoints, a dynamic direct sampling method and a modified particle filter method. Both approaches can effectively recover the motion trajectory. Rigorous theoretical justifications are presented for the mathematical modelling and the proposed recovery methods. Extensive numerical experiments are conducted to illustrate the promising features of the two proposed recognition approaches.


Introduction
Among various human-computer interaction technologies, most people prefer to interact with computers in a more personal way, e.g. by using voice, touch and gesture rather Γ z 0 Ω Figure 1: Schematic illustration of the proposed input/instruction technology using a moving emitter. than a mouse or a keyboard. In particular, people take an interest in gesture-computing technology since it enables people to communicate with the machine and interact more naturally without the help of any mechanical devices. Some existing technologies have been developed using cameras to capture the human body gesture, and then using imaging processing algorithms to interpret the body language in order to understand the instructions or inputs; see [24] and the references therein. We also refer the interested readers to a recent article by Liu et al. [19], where they first attempted to achieve gesture recognition by using inverse scattering techniques. In fact in [19], instead of using of a camera, it is proposed that one uses wave probing to identify the body gesture. However, the proposed method in [19] is mainly suitable for static instructions for computers, but not suitable for dynamic text inputs. In order to develop a novel gesture-computing method that is suitable for both instruction and input, Guo et al. [13] propose to use a moving emitter generating an acoustic point wave and then to identify the motion trajectory of the emitter that carries the intended instruction/input. The purpose of the present article is to further develop the idea in [13] to produce much more practical, effective and efficient gesture-computing methods. Fig. 1 provides a schematic illustration of the proposed gesture-computing technology. In order to give the specific input or instruction to the computer, a human being wears an emitter on his/her finger or holds a digital pen that generates a time harmonic point charge. The input/instruction is performed through moving the finger or the digital pen. There are sensors installed on the computer that timely collect the electromagnetic wave data generated by the motion of the point charge. Then the computer recognizes the instruction/input by determining the motion trajectory of the dynamic point charge from the collected electromagnetic field data. Mathematically, the motion trajectory identification can be modelled as an inverse problem where one intends to identify a moving emitter from the measurement of the electromagnetic wave fields generated by the emitter. There are several practically important issues that should be incorporated into the design and modelling. First, the choice of the point charge is a critical ingredient in our design. The physical principle of using a moving charge to generate electromagnetic fields can be found in [12]. We choose to utilize a time-harmonic point charge for two considerations. On the one hand, the time-harmonic point charge will generate electromagnetic waves with a fix frequency. This enables us to distinguish the signals due to the motion of the emitter from the possibly various background signals such as those from radios, telephones and televisions. On the other hand, in a certain practical scenario, we can show that by properly choosing the frequency of the waves, one can eliminate the scattering influence from certain inhomogeneous background scatterers such as the body of the human being who performs the input/instruction. Second, the emitter should be modelled as moving within an unknown inhomogeneous background medium as described above, and the trajectory identification should be independent of the background medium. Third, the measurement data should be collected only in a limited aperture on a surface. Indeed, as described earlier, the sensors used to collect the wave data are installed on the computer. Hence, one would only have limited-view data for the inverse problem. Moreover, we shall see in the subsequent study that the location of the measurement surface is also an important ingredient in our design. Finally, the trajectory identification should be conducted in a timely manner. All of those challenging issues distinguish our study from the existing ones in the literature on inverse source problems associated with electromagnetic wave probings. We are aware of some existing theoretical and computational developments on the identification of an unknown source from the measurement of the generated electromagnetic wave data away from the source. Those problems arise in various practical applications; see e.g. [1,2,4,9,10,14,15,20] and the references therein. Finally, as mentioned earlier that in [13], the similar input/instruction technology has been proposed using a moving acoustic emitter. Clearly, using electromagnetic waves for the proposed technology is more realistic and practical. Indeed, according to the mathematical and theoretical analysis in [13], the applicability of the proposed technology is limited if acoustic wave is used due to the low-speed of propagation. The use of electromagnetic wave shall overcome this problematic issue. Nevertheless, we would like to emphasize that our study is mainly concerned with the conceptual design and theoretical analysis, and the issue on engineering realization of the proposed technology is beyond the scope of this paper.
For the dynamic inverse electromagnetic problem described above, we develop two methods for the trajectory identification: one is a direct imaging method motivated by our theoretical analysis and the other one is a modified particle filter method motived by the particle filter methods developed in [6,16,23] for various dynamic inverse problems. Both methods are shown to be effective and efficient by extensive numerical examples. Finally, we also note some existing results on recovering moving targets in different contexts [3,11,21].
The rest of the paper is organized as follows. Section 2 is devoted to the mathematical modelling of the moving-emitter-based gesture recognition technique as well as some necessary theoretical analysis. In Section 3, we develop two recovery methods for the trajectory identification. In Section 4, we conduct numerical experiments to illustrate the proposed instruction/input technology as well as the effectiveness and efficiency of the methods developed for the identification.

Mathematical modelling and theoretical analysis
In this section, we build up the mathematical modelling and present some theoretical analysis for the proposed sensor-based gesture recognition technology.

Mathematical modelling
Assume that the point charge emanates a causal sinusoidal signal where t denotes the temporal variable, and ω 0 ∈ R + denotes the frequency. Suppose that the electric charge density ρ has the following form where x denotes the spatial variable, and z : R + → R 3 is a C 2 smooth function that signifies the instantaneous position of the point charge at time t and δ is the Dirac's delta distribution. Let the moving trajectory of the point charge be where T ∈ R + is the terminal time of the motion. As the point charge moves, the electric current density satisfies where v(t) = dz(t)/dt is the instantaneous velocity. Suppose a human being is interacting with the gesture computing device, his/her body could be modelled by a bounded moving domain Ω(t) ⊂ R 3 (0 < t ≤ T ) such that its boundary is time-varying. In what follows, we set The receivers are fixed and located on a static open surface Γ ⊂ R 3 \Ω such that Γ∩Λ z = ∅; see Figure 1 for the schematic illustration of the problem setting. We set ǫ(x, t) ∈ L ∞ (R 3 × [0, T ]) and µ(x, t) ∈ L ∞ (R 3 × [0, T ]) be positive functions to, respectively, signify the electric permittivity and magnetic permeability at the space point x and time point t; and σ(x, t) ∈ L ∞ (R 3 × [0, T ]) be a non-negative function to signify the electric conductivity at (x, t). We assume that the medium is homogeneous in the background space R 3 \ Ω(t), that is, where ǫ 0 ∈ R + and µ 0 ∈ R + are, respectively, the constant permittivity and permeability of the background space. Biologically, it is reasonable to assume that µ(x, t) = µ 0 in Ω(t).
In the setup described above, the electric current J(x, t) generates electromagnetic waves that propagate in the space. We denote by E(x, t) and H(x, t) the electric field and the magnetic field in R 3 × (0, T ], respectively. The electromagnetic field satisfies the following Maxwell system, 5) and the initial condition E| t=0 = 0, H| t=0 = 0. (2.6) The Maxwell system (2.5)-(2.6) is well understood, and we refer to [18] for the study on its well-posedness and especially the unique existence of a pair of solutions (E, H) ∈ . The trajectory identification associated with the prosed input/instruction technology can be formulated as follows, That is, by monitoring the change of the electromagnetic field on the surface Γ generated by the emitter, we intend to recover the motion trajectory of the emitter. Here, we would like to emphasize that the recovery should be independent of the inhomogeneity (Ω; ǫ, µ, σ).

A theoretical model approximation
As emphasized at the end of Section 2.1, the recovery for the inverse problem (2.7) should be independent of the background inhomogeneity (Ω; ǫ, µ, σ). Clearly, in the measurement data E Γ×(0,T ] , there are both scattering data due to the moving emitter and the background inhomogeneity that are coupled together. Furthermore, the background inhomogeneity could be changed due to the change of the human being that performs the input or instruction. Next, we present some practical conditions that the scattering influence due to the background inhomogeneity can be eliminated. To that end, we introduce the following Maxwell system for where J is given in (2.4). The system (2.8) describes the electromagnetic wave propagation generated by the moving emitter in the free space without any inhomogeneity presented. In the following, we shall show that under certain practical conditions, the difference between the two fields, E in (2.5) and E 0 in (2.8) can be small; that is If (2.9) holds true, then clearly the inverse problem (3.1) can be approximately replaced by the following one Here, we note that in (2.10), there is no inhomogeneity presented. Next, we shall show that (2.9) can indeed hold under certain conditions. By (2.1), (2.2) and (2.4), we know that the electromagnetic field is time-harmonic with frequency ω 0 . Hence, we can assume that (cf. [22]) and Furthermore, in order to show (2.9), we could only consider a fixed instant, say t 0 . In the following, at t 0 , we still use the notations ǫ(x), µ(x), σ(x) and Ω to denote the optical parameters of the inhomogeneity and its support respectively, and this should be clear from the context. By plugging (2.11) into (2.5), one then has where J(x) signifies the electric current density at the instant t 0 . It can be shown that J(x) takes the following form 14) It is known that E and H in (2.13) satisfy the following so-called Silver-Müller radiation condition (cf. [7,22] We can show that be respectively the solutions to the Maxwell systems (2.13) and (2.16). Then we have the following results.
Remark 2.1. In order to give a completely rigorous justification of (2.9), one should make use of the Fourier transform to convert the time-domain problems (2.5) and (2.8) into their frequency-domain counterparts. However, since the optical parameters ǫ and σ may also vary according to time, one may meet difficulties in such a Fourier argument. Hence in Theorem 2.1, we simplify our study by using the time-harmonic assumption (2.11) as well as by considering a fixed instant t 0 . Nevertheless, it is unobjectionably to claim that under the conditions specified in Theorem 2.1, (2.9) should also hold true as long as the emitter is moving within a certain bounded domain.
Remark 2.2. By Theorem 2.1, (i), if the inhomogeneity is away from the measurement surface by a reasonable distance, then the scattering influence from the inhomogeneity can be neglected under low frequency emission. We note that |E 0 (x)| ∼ 1/|x − z 0 | for |x − z 0 | sufficiently large, and hence (2.17) indicates that the difference between E and E 0 on Γ is indeed nearly negligible. The inhomogeneity is mainly used to model the body of the human being who performs the instruction/input, and this means that one can keep a reasonable distance away from the computing device when performing the instruction/input. Next, in Theorem 2.1, (ii), we note that the assumption on z 0 is obviously a reasonable one. The assumption on ǫ and µ being constant in Ω can be relaxed as being with small variations. In such a case, one can choose a small frequency ω 0 , and in doing so, the scattering influence from the inhomogeneity can also be eliminated. All in all, Theorem 2.1 indicates that one should make use of lowfrequency emission, and the human being that performs the input/instruction should keep a reasonable distance away from the computing device in the design of the proposed gesture-computing device. This has been confirmed by our numerical experiments in what follows under practical and mild conditions on ω 0 and L, not as restrictive as the theoretical assumptions in Theorem 2.1.

Motion trajectory recovery
We are now in a position to present two imaging schemes for qualitatively determining the motion trajectory Λ z by knowledge of |E(x, t; Λ z )| , namely, We would like to point out that we are using the data without phase information, i. e., only the strength of the wave field is available. For inverse problems with phaseless data, we refer to the recent work [17].

Imaging via the direct sampling method
where D ⊂ R 3 is the static compact sampling region, such that Λ z ⊂ D. In the present study, for (ẑ, t) ∈ D × (0, T ], we propose the indicator function as follows: where α > 0 is a parameter and |Γ| denotes the surface area of Γ. Theorem 3.1. Let E(x, t; Λ z ) be the measurement data for (x, t) ∈ Γ × (0, T ], corresponding to a moving point charge described in (2.5) and let I(y, t) be defined in (3.3). Define Let the parameters be chosen such that (2.17) and (2.18) in Theorem 2.1 hold, then for each fixed t 0 ∈ (0, T ] and any given ǫ > 0, there exists an ω 0 > 0, such that |I 0 (y, t 0 ) − I(y, t 0 )| < Cε, ∀y ∈ D, where C > 0 is a constant depending on t 0 , α and |Γ|.
Remark 3.1. Theorem 3.1 illustrates that I(y, t) is quite close to I 0 (y, t). Thus the indicator I should inherit I 0 's maximum indicating behavior. For any fixed t 0 ∈ (0, T ], it can be easily deduced that I 0 (y, t 0 ) attains its maximum when y = z(t 0 ). Therefore, I(y, t 0 ) attains its maximum when y ≈ z(t 0 ).
Based on Theorem 3.1 and Remark 3.1, we proposed the first trajectory reconstruction scheme for the inverse problem (3.1), see Algorithm 1. It is remarked that the proposed sampling scheme could be able to recover the moving trajectory if the sampling grid is sufficiently fine. However, the computational cost would be relatively high when the indicator functional is evaluated over a very fine mesh in Step 3. To speed up the reconstruction process, in the next subsection, we shall develop an alternative meshless approach, namely, the modified particle filter method.

Imaging via the modified particle filter method
In this subsection, we are going to reformulate the problem of trajectory reconstruction as a stochastic inverse problem based on a probability space. The method provided here could be considered as a modified version of the classical particle filter method. For the ease of the readers, we recall some essential rudiments of the particle filter method in Appendix A.
The underlying assumption of the stochastic inverse problem is that the time-discrete reconstruction is a sequence of states of random variables, which sample some particular probability distribution. Assume that the measured data were collected at N t discrete Algorithm 1: Reconstruction of the trajectory via direct sampling Step 1 Properly choose a low frequency ω 0 and a point charge of the form (2.1).
Step 2 Set the emitter in motion following a specific path Λ z , depending on the desired input/instruction. The sensors collect the electric field data |E(x m , t n ; Λ z )| at the measurement points {x m } ∈ Γ and a sequence of discrete time points {t n } ∈ (0, T ]. Step 3 Select a sampling mesh T h in a region D such that Λ z ⊂ D. For each time point t n , evaluate the imaging functional I(y, t n ) defined in (3.3) for each z ∈ T h . Step 4 Locate the global maximum point z n of I(y, t n ) for y ∈ T h , which is an approximation to z(t n ).
Step 5 The ordered chain {z n } Nt n=1 forms a discrete version of the reconstruction of Λ z . times t n = n∆t, n = 1, 2, · · · , N t , with time step ∆t = T /N t . Now let us consider a time-discrete Markov chain of states {ξ n } corresponding to time t n for n = 1, · · · , N t . The transition probability of the Markov chain is denoted by p(ξ n | ξ n−1 ), i.e., the probability to draw ξ n at time t n when ξ n−1 was drawn in its previous step t n−1 .
To specifically quantify the trajectory of the point emitter by the strategy of Markov chains, we consider the following Gaussian random walk model {ξ n } satisfying ξ n = ξ n−1 + G n , n = 1, 2, · · · .
Here G n denotes a Gaussian distribution with zero mean and covariance matrix Υ := diag(γ 2 , γ 2 , γ 2 ) where γ > 0 controls the step size in the spatial evolution. Correspondingly, the transition density distribution is Considering the observation model in light of (3.3), we define the density function of W n by p(|E(x, t n ; Λ z )| | ξ n ) = I(ξ n , t n ) where I is defined in (3.3).
By Appendix A, the posterior probability density function is approximate to where ξ (i) n ∈ D, 1 ≤ i ≤ N s , are N s random sample points at the instance t n . Based on the above discussion, we next present the trajectory reconstruction scheme via the modified particle filter method in Algorithm 2.
Algorithm 2: Reconstruction of the trajectory via particle filter.
Step 1 Properly choose a low frequency ω 0 and a point charge of the form (2.1).
Step 2 Set the emitter in motion following a specific path Λ z , depending on the desired input/instruction. The sensors collect the electric field data |E(x m , t n ; Λ z )| at the measurement points {x m } ∈ Γ and a sequence of discrete time points {t n } ∈ (0, T ]. Step 3 Initialization: for n = 0, draw an initial random sample {ξ (i) n } Ns i=1 from a uniform distribution in a region D such that Λ z ⊂ D.
Step 5 Resample: let U be the uniformly distribution density function, draw a random number q ∼ U ([0, 1]). For every i, set ξ Step 6 If (n + 1)∆t ≥ N t , then the reconstruction is finished. Otherwise, set n = n + 1 and repeat from Step 4.
Step 7 The ordered chain {z n } Nt n=1 forms a discrete version of the reconstruction of Λ z .
Remark 3.2. The critical feature of the modified particle filter method is that we use multi-measurements to identify the trajectory at discrete instants. Using the discrete observation points, the discrete version of the indicator function in (3.3) could be written as

Numerical examples
In this section, we will present several numerical examples to demonstrate the feasibility and effectiveness of the proposed method.
All the following numerical experiments are carried out using MATLAB R2016a on a Lenovo workstation with 2.3GHz Intel Xeon E5-2670 v3 processor and 128GB RAM. Synthetic electromagnetic field data are generated by solving direct problem (2.5) by using the quadratic finite elements on a truncated spherical domain enclosed by absorbing boundary condition. The mesh of the forward solver is successively refined till the relative error of the successive measured electromagnetic wave data is below 0.1%. To test the stability of the proposed reconstruction algorithm, Gaussian noise was point-wisely added to the synthetic data, that is,  Figure 2 for the geometry setting. At last, let N m = 400 and N t = T /∆t, where the time step ∆t is set to be 0.1 s. Assume that the average velocity of the point emitter is 0.5 m/s, hence the step size parameter γ is chosen as 0.5 m/s × ∆t = 0.05 m. Heuristically, we would like to remark that in our proposed imaging methods, the parameter α is such chosen that α ∼ c 2 0 , then they can produce fine reconstructions. Thus α = 10 16 is used. Finally, as a convention in the following figures, the blue points denote the measurement points and the blue cubic domain denotes the person's motion domain.
To show the accuracy of the proposed methods, we also define the discrete relative L 2 error between the exact trajectory and the reconstruction by . Example 1. In the first example, we compare the direct sampling method and the particle filter technique in reconstructing a moving trajectory. Consider a simplified scenario that a person is wearing an emitter on one of his/her finger and moving the finger to write an Arabic number "3", which is modelled by the trajectory: The indicating scatter plots of both direct sampling method and particle filter method are shown in Figure 3. In the following figures concerning the problem geometry, some 2D projections (shadows with gray color) are also added in order to facilitate the 3D visualizations. Both methods are able to identify the moving trajectory even if the measured data has 10% noise. Through comparing the error between the exact trajectory and the reconstructed trajectory in discrete instants, one can find that the recovery results of both methods have similar error fluctuations, see Figure 3(d). In order to exhibit the accuracy and computational cost of the reconstructions, the relative error in discrete L 2 norm and the CPU time are listed in Table 1. Reconstruction of a conical spiral. In the previous example, we only study moving trajectory in a plane. This example is intended to demonstrate the performance of both methods for reconstructing a complicated 3D trajectory. Here the moving trajectory is like an upward conical spiral, z(t) = ( 0.05t cos 2t + 1, 0.05t sin 2t + 1, 0.1t + 0.5 ) , t ∈ (0, 10 s]. From equation (4.2), it is clearly that the velocity |v(t)| is monotone increasing. As shown in Figure 4(b)-(d), both methods could produce satisfactory reconstructions. In addition, Figure 4(e) shows that the particle filter method works better with a higher sampling density, namely, the number of particles. However, the computational cost increases as the number of particles increases, so a proper number of particles should be considered.  moving path is as follows , 3 5 cos t + 1 , t ∈ (0, 6 s], This example is challenging since there exists a sudden skip in the moving trajectory, where the skip distance is far longer than the normal transition step size. The dotted lines in Figure 5(a) represent the skip trajectory. As discussed in Example 1, the standard deviation γ usually depends on normal transition step size. We note that the particle  Figure 6: Reconstruct the moving trajectory by the particle filter method and snapshots at different instants, where the black points denotes the particles. (a) t = 3.0 s, (b) t = 6.6 s, (c) t = 9.0 s, (d) t = 9.1 s, (e) t = 9.2 s, (f) t = 12.0 s filter method will stop tracking the real trajectory path when there is a sudden skip. An efficient method is to redistribute the particles such as Step 3 in Algorithm 2, see Figure 6(d). Figure 6 illustrates that the particle filter method could recover the trajectory after redistribution. However, in Figure 5(d), there exists large perturbation at the skip instants, this is due to the shortcoming of the particle filter method.
Remark 4.1. In this paper, the reconstructed trajectories are only in the discrete form. To obtain a smooth curve as the reconstructed trajectory, one could resort to the truncated Fourier approximation for post-processing the discrete points {z n } Nt n=1 [13]. The implementation of the direct sampling method could also be accelerated by incorporating the sequential or parallel local sampling strategies, see [13] for more details.

Conclusions
In this work, we develop a conceptual framework of the novel gesture recognition techniques using electromagnetic waves. The model problem is formulated as an inverse source problem of determining the moving trajectory of a moving point charge. Two methods are proposed to deal with the inverse problem, respectively, in the deterministic and statistical manners. Mathematical justifications are presented and extensive numerical examples are provided to validate the effectiveness and efficiency of the methods.
We would like to remark that this work is mainly from a theoretical point of view and the physical and engineering realizations are beyond the scope of our current study.
where f : R nx → R nx and g : R ny → R ny are two measurable functions. V t ∈ R nx denotes the process noise and W t ∈ R ny denotes the measurement noise. Let p V be a density distribution of V t and p W be a density distribution of W t and define p(x t | x t−1 ) := p V (x t − f (x t−1 )), p(y t | x t ) := p W (y t − g(x t )). (A.1) The optimal nonlinear filtering problem is to find the posterior probability density function p(x t | y 1:t ) at the state x t from the observation data y 1:t := {y 1 , · · · , y t }. Here the posterior probability density function is given by Bayes' formula [8]: p(x t | y 1:t ) = p(y t | x t )p(x t | y 1:t−1 ) The key idea behind the particle filter method is to use a set of samples with weight to approximate the posterior probability density function in (A.2). Given N random particles {x , correspondingly, the posterior probability density function with associated weights {w where δ x denotes the delta-Dirac mass located at x and .
(A.4) Equation (A.4) indicates that it is possible that only one particle has a significant weight value after several recursive steps. A resampling stage [8] allows to replace the samples with low weights by copies of the samples with high weights, which ensures more particles in statistically significant areas. The classical particle filter algorithm proceeds in three main steps, see Algorithm PF.
Algorithm PF: The classical particle filter method.
Step 3 Resampling. Sample where the posterior probability density function is