Multi-Channel Surface EMG Spatio-Temporal Image Enhancement Using Multi-Scale Hessian-Based Filters

Surface electromyography (sEMG) signals acquired with linear electrode array are useful in analyzing muscle anatomy and physiology. Most algorithms for signal processing, detection, and estimation require adequate quality of the input signals, however, multi-channel sEMG signals are commonly contaminated due to several noise sources. The sEMG signal needs to be enhanced prior to the digital signal and image processing to achieve the best results. This study is using spatio-temporal images to represent surface EMG signals. The motor unit action potential (MUAP) in these images looks like a linear structure, making certain angles with the x-axis, depending on the conduction velocity of the MU. A multi-scale Hessian-based filter is used to enhance the linear structure, i.e., the MUAP region, and to suppress the background noise. The proposed framework is compared with some of the existing algorithms using synthetic, simulated, and experimental sEMG signals. Results show improved detection accuracy of the motor unit action potential after the proposed enhancement as a preprocessing step.


Introduction
Human muscle is made up of a certain number of motor units (MU) which are comprised of a set of muscle fibers innervated by the motor neuron. The electrical activity of these motor units is called motor unit action potential (MUAP) [1]. The sum of these potentials for a muscle is called electromyograph (EMG). Surface EMG (sEMG) signals are recorded by placing the electrodes on muscle surface [2]. Multi-channel EMG signals are often recorded using a linear array of electrodes [3]. These multi-channel sEMG signals can be represented using spatio-temporal two-dimensional (2D) images. Rows and columns of these images represent the sEMG electrodes and the time samples respectively, while the gray level corresponds to the sEMG amplitude [3][4][5]. Representing sEMG signals as 2D images can enable the use of state-of-the-art digital image processing techniques for further processing and features extraction.
High-density and multi-channel sEMG signals are useful in muscle geometry estimation, detection of muscle fibers alignment, motor unit length, and innervation zones [3]. Accurate detection and estimation of these muscle's parameters depend on the precise localization of MUAPs in sEMG images. The shape of the MUAP in sEMG images is distorted by noise and interference due to nearby MUAPs, which make it difficult to detect and estimate these parameters [3,4].
Due to interference of the neighbor MUAPs, the MUAP structure has different gray levels and it is difficult to be distinguished from the background noise (see Figure 1a,b) [2]. sEMG signals have varying amplitudes due to variations in the size of mucosa in case of external anal sphincter (EAS) muscle or pinnation of the fibers in case of skeletal muscles [3]. Moreover, different noise sources can contaminate the sEMG signals, which ultimately affect the MUAP structure (see Figure 1b) [6,7]. The spatio-temporal 2D sEMG images of the MUAP propagation structure resembles a line-like structure with either a single line or a v-like shape. The MUAP structure needs to be detected before estimating the muscle's physiological and anatomical parameters accurately. The sEMG signals are affected by different factors such as sampling frequency in space and time, the electrode size, interpolation, and truncation effects, and the interelectrode distance [8]. Moreover, sEMG signals are also affected due to severe crosstalk between deep and superficial muscle activation [9]. Due to these limitations, preprocessing of sEMG is needed before further processing. The visualization of the MUAP structure is significantly improved in the preprocessed enhanced spatio-temporal images. This facilitates the segmentation of the MUAP regions due to increased difference between the background and foreground. It also improves the binarization and edge detection. A review of the literature reveals significant research activity in applying digital image processing techniques to sEMG spatio-temporal images. These techniques are used to enhance the sEMG spatio-temporal images and detect the MUAP structure automatically. Authors in Reference [3] used the 2D cross-correlation technique to detect the MUAP and motor unit innervation zone (MU IZ) from sEMG images. In References [10,11], the authors used digital image processing techniques to automatically detect MU IZ and to estimate muscle fiber conduction velocity. Various digital image processing techniques are also used for enhancing and detecting line-like structures in 2D images [12][13][14]. Most of these methods are based on fixed scale techniques. These techniques are not suitable for MUAP enhancement and detection in sEMG images due to the non-uniform width of the MUAP propagation pattern. In a recent study [4], the authors have analyzed several line-enhancing filters for enhancing and detecting the MUAP propagation pattern in sEMG spatio-temporal images. Some authors have used the Gabor filter [15,16], Hermite Shape filter [17], and steerable filters [18,19] to enhance the sEMG images. From the results of these studies, it is evident that the Hermite shape filter is the best candidate to represent and enhance the MUAP propagation pattern. However, it was found that the Hermite shape filter and other filters analyzed in Reference [4] are unable to retain the polarities of the EMG signals. In some recent studies [20][21][22], multi-scale methods based on Hessian are proposed for enhancing line-like structure in digital images. However, none of these methods are utilized for the MUAP enhancement in the EMG images. In Reference [5], the authors have proposed that sEMG signals can be enhanced using the Hessian-based filter. In this preliminary work, only the bright portion of the MUAP is enhanced. Moreover, the analytical justification of the Hessian filter and its performance analysis on synthetic, simulated, and experimental sEMG was not studied. This proposed study is an extension of the work in Reference [5]. It presents a similar framework based on the Hessian matrix that enhances a Gaussian line-like structure (MUAP propagation structure) in sEMG spatio-temporal images with analytical justification. A novel approach is presented based on the multi-scale Hessian matrix for detecting and enhancing the MUAP propagation structure in the spatio-temporal sEMG images. The proposed work is an extension of the initial work done by the authors in Reference [20] and then by the authors of References [22,23] for enhancement of linear and tubular structure in 2D and three-dimensional (3D) images. The same approach is also successfully used in Reference [21] to enhance and segment vessels in retinal digital images. The proposed method is compared with different linear structure enhancement methods using synthetic, simulated, and experimental sEMG signals. The proposed method has an average accuracy of 98% in enhancing and detecting the MUAP propagation region in sEMG spatio-temporal images.
This manuscript is organized into 4 sections. Section 2 describes the analytical justification of using the Hessian matrix for MUAP enhancement. It provides a detailed mathematical background for the design of the proposed multi-scale filter. Section 3 analyzes the results of the proposed filter for the enhancement of MUAP in the synthetic, simulated, and experimental signals in detail. Section 4 provides the conclusion of the proposed work.

Model of the MUAP Propagation
The MUAP propagation, in spatio-temporal sEMG images, has bright (positive) and dark (negative) portions over gray (zero) background, and in general, these two portions can be modelled separately based on a second-order structure, i.e., Hessian matrix of the pixels of the sEMG image. We assume that in spatio-temporal sEMG images, the brightness of the bright portion is high in the middle and gradually decreases toward the end (Figure 2a). Therefore, the bright portion of the MUAP can be modelled as a line-like Gaussian structure transversal to its axis, i.e., x-axis, as given in the following equation (see Figure 2a): where is the variance of the Gaussian which is proportional to the width of the MUAP propagation profile shown in Figure 2a. Similarly, the dark region is assumed to be darkest in the middle and then decreases gradually as approaching towards the borders (Figure 2b).
The Hessian matrix for the line-like inclined structure, like a MUAP ( Figure 2a) or any other image, can be computed as: As from Equation (1), ( , ) = , thus the derivatives in Equation (2) can be obtained analytically as follows: Replacing all these expressions in Equation (2), we get: The eigenvalues and eigenvectors for the Hessian matrix in Equation (3) are computed analytically as: It is clear form Equation (4) that for a line-like Gaussian structure, λ is zero and λ is negative for the pixels inside the Gaussian region, i.e., < . Based on these eigenvalues, each pixel can be classified as belonging either to the MUAP propagation region or background region. These regions can then be enhanced using the multi-scale filter described in the next subsection. A synthetic 15- pixels-wide image with Gaussian profile is shown in Figure 3. For the line-like Gaussian structure, the second eigenvalue of the Hessian for each pixel is negative, and zero for the background pixels (see Figure 3c). This agrees with the above analytical proof.

Filter for Enhancing the MUAP Propagation Region
In single differential sEMG spatial-temporal images, the bipolar (positive and negative) MUAP structure looks like a bright and dark region on a gray background, respectively (see Figure 1b) [3,4]. The propagation of the positive amplitude appears as a bright line-like structure, while the negative EMG amplitude appears as a dark line-like structure in sEMG images, as shown in Figure 1b. These bright and dark regions in spatio-temporal sEMG images can be approximated by a Gaussian waveform with EMG amplitude maximum at the middle and decreasing gradually towards the end. Second-order derivative of the image with respect to the x-and y-axis, i.e., Hessian, is then used to design a filter whose response is maximum for the MUAP region, while minimum for the background, resulting in enhancing the MUAP structure. The filter is based on the eigenvalues of the Hessian matrix of each pixel of sEMG spatio-temporal sEMG images. The same method is applied in References [20][21][22][23] to design a filter that enhances line-like Gaussian structures in digital images, like retinal images. Some of the methods use fixed scale and the (nonlinear) combinations of finite difference operators applied in a set of orientations for enhancement of digital images [24]. All these methods have shown great performance in enhancing linear Gaussian structures with different scales but none of these is applied to sEMG spatio-temporal image enhancement. In sEMG images, the width of the MUAP propagation region is not fixed due to the interference of other MUAPs and noises; thus, a fixed scale approach with a fixed sigma will not be able to enhance the EMG signals due to variations in the MUAP [3]. Therefore, the response of the enhancement filter may not be the highest for every pixel at a fixed scale. On the other hand, with the use of a multi-scale approach where a range of sigma values are used, whenever the sigma is equal to the width of the MUAP structure, the response of the enhancing filter is highest.
In the proposed method, first, to achieve the multi-scale Hessian matrix at each pixel, the input sEMG spatio-temporal image is convolved with the first-and second-order derivative of a Gaussian at various scales (that is with multiple σ values) to estimate the first-and second-order derivative of each pixel with respect to the x-and y-axis. The second-order derivative approximates the Hessian. Eigenvalues λ1 and λ2 of the Hessian matrix are then computed for each pixel to find whether it is a background or foreground pixel. For a pixel in the MUAP region in a sEMG image, the first eigenvalue of the Hessian matrix is almost zero, while the second eigenvalue is negative and has a big absolute value for the bright portion of the MUAP. The second eigenvalue is positive and bigger for the dark portion of the MUAP. Norm of the eigenvalues is then computed, and the filter response depends on the value of the norm. In the MUAP region, the value of the norm is big because one of the eigenvalues is high, thus the filter response is high.
For the spatio-temporal sEMG image I(x,y), the Taylor series up to the second order can be written as (Equation (5) [20]. The bright line-like structure on grey background for each pixel is represented by a high negative magnitude of second eigenvalue and a low absolute first eigenvalue. Based on these eigenvalues, a higher response filter for the line-like region can be defined as having almost zero response in the background. To track the MUAP propagation for varying grey levels and varying widths, the Hessian is approximated using a Gaussian kernel with varying scales (scale space theory) [20,25]. The Gaussian Kernel can be represented as: where σ is the width of the Gaussian vertical to the axis of the MUAP structure. The Hessian matrix of each pixel of an image is analytically given as: In scale space theory, the Hessian is computed by the convolution of the sEMG spatio-temporal image with the second-order derivatives of the Gaussian, described in Equation (6) for a chosen scale σk. Mathematically, this approximation of the Hessian is as follows [18,19]: where γ is a parameter used for normalization [25]. The σ k is the scaling factor proportional to the width of the MUAP structure in sEMG images. Using the eigenvalue decomposition, the Hessian of each pixel of the sEMG image is decomposed to eigenvectors and eigenvalues. Let these eigenvalues for the 2D sEMG image be λ and λ , such that |λ | < |λ |, with associated eigenvectors V1 along the MUAP structure while V2 is normal to the MUAP structure. As shown in Section 2.1., the λ is zero while the λ is negative inside the bright tubular structure, thus classifying a pixel as either belonging to MUAP or the background region. A filter is now designed to enhance the MUAP structure by first calculating the norm (N) of the Hessian, as given below: ( , ) = (λ 1 + λ 2 ) The response of the proposed filter which enhances the bright portion of the sEMG image is as follows: where R B =λ 1 /λ 2 . This ratio is nearly zero for the MUAP region with smaller value of λ and larger absolute value of λ . The value of parameters α and c are adjusted to the scale response of the filter. The multiplication operation (Equation (10)) ensures that the pixel corresponding to the MUAP region is enhanced only. Both factors of the product are nearly equal to one for the pixel belonging to MUAP region, and otherwise equal to zero for the background. This filter response is calculated for σk of the Gaussian kernel with varying pixels ranging between 1 and 10. This range is chosen according to the possible widths of the MUAP [3]. When the width of the MUAP is like the σk, then the filter output (Equation (10)) is the highest value, which is the desired output. Thus, the final response is: Similarly, the Hessian matrix for the dark region (negative signal amplitude) of the MUAP propagation is estimated and the corresponding eigenvalues are obtained. The eigenvalue λ is zero while the eigenvalue λ is positive for the dark region. The filter output for enhancing the dark region of the MUAP propagation is: After choosing the negative minimum (which is the highest response for some scale, σk), the final response is: ( , ) = min ( ( ( , ), )) (13) Enhancing this dark region of the MUAP separately holds a negative polarity. The overall response of the system for enhancing both the bright and dark portions of the MUAP can be obtained by adding Equation (11) and Equation (13): One advantage of enhancing the bright and the dark portions separately is that noise in the background is suppressed twice, which ensures a clearer image of the MUAP.

Results and Discussions
sEMG images have poor MUAP visualization due to different noises and destructive interference of nearby MUAPs. Sometimes, the width of the MUAP propagation region is nonuniform due to a mix of destructive and constructive interferences and because of the presence and absence of noises at different electrodes [3]. Thus, a multi-scale filter is needed to preserve and enhance the MUAP regions at different scales proportional to the MUAP width and to attenuate the background noise.
The performance of the multi-scale Hessian-based filter is studied and analyzed in different scenarios, such as synthetic images with gaussian propagating profile, simulated EMG signals using the model proposed by the authors in Reference [26], and the experimental EMG signals recorded from different human muscles. As discussed in the previous section, the MUAP propagation across the electrodes with respect to time appears as a Gaussian line-like structure. The second eigenvalue of the Hessian matrix is negative for the MUAP propagation region (see Equation (4) and Figure 3c). To verify this, a synthetic multi-channel signal with a propagating Gaussian profile is generated and contaminated with additive white gaussian noise. In the presence of this high-level noise with SNR of 6.996 dB, the signal is passed through the proposed multi-scale filter.
The results show that the multi-scale filter has preserved and enhanced the region with Gaussian profile and has suppressed the background noise. This is because the second eigenvalue of Hessian is negative and its magnitude is higher in the Gaussian region, and the output of the proposed filter is also high, while in the background region, both the RB and N parameters of the eigenvalues of the Hessian matrix are very small; thus, the filter output is also nearly zero. To further quantify the performance of the proposed filter for various noise levels, noises are added to the synthetic signal of Figure 4a and the input signal to interference ratio ( ) is computed. The is the ratio between the power of the clean input signal and the noise added to the signal. The signal is then enhanced using the proposed filter and the remaining interference, Ri(t), after enhancement for each channel is calculated, as follows, by subtracting the clean simulated sEMG V(t) from the enhanced signal ( , , ), as follows: . Using this ( ), the performance of the proposed filter is evaluated by calculating the output signal to the interference ratio ( ) of the enhanced signals. This is named as the output SIR. The output SIR is the ratio of the power of the simulated signal (clean original signal) to the power of the b) a)

d) c)
remaining interference (noise left in the signal after enhancement, i.e., ( )). The output SIR is computed as: (16) where ( ) is the power of the original simulated signal and ( ) is the power of the remaining interference left in the signal after filtering. The result of this analysis for an input SIR varying from 3 to 25 dB for the signal in Figure 4a is shown in Figure 5. It can be seen from the result that initially, at low-input SIR, the proposed filter is able to enhance the signal with less improvement in the output SIR; however, as the input noise level decreases, the performance of the proposed filter is better and the output SIR increases. The higher output SIR shows that this filter may be a better choice as a preprocessing step for EMG signal analysis. The results of the proposed filter for a multi-channel sEMG signal acquired from the external anal sphincter muscle acquired using an anal probe developed in the LISiN politecnico di Torino Italy [3] are shown in Figure 6. The signals were acquired from the EAS muscle of 478 pregnant women at clinical partners from five European Countries (Germany, Italy, Latvia, Slovenia, Ukraine). Each clinical partner obtained the approval from the local ethical committee. Each subject was informed about the study protocol and signed an informed consent prior to the experiment [27]. It can be seen in this example signal that the MUAP propagation pattern in the original signal is non-uniform at different locations due to varying size of the mucosa at different locations of the muscle [3]. The proposed filter has successfully enhanced the MUAP regions and the background noise is suppressed, as is evident from Figure 6c. The proposed multi-scale filter is finally evaluated as a preprocessing step for the segmentation of the MUAP propagation region. Graph cut segmentation as recommended in Reference [28] and discussed in detail in Reference [29] is applied to the enhanced image using the proposed filter. The proposed multi-scale filter is also compared with some of the line-like structure enhancement algorithms, namely Hermite Filter, Gabor filter, and Steerable filters. The results of these algorithms for a multi-channel sEMG of the right trapezius muscle are shown in Figure 7.  From the results of Figure 7, it can be summarized that the output of the proposed multi-scale filter is better because it enhanced both the bright and dark portions of MUAP and retained its negative and positive polarity. The Gabor filter added distortion to the image while other algorithms were unable to retain the same polarity of the EMG values, except the proposed filter. From Figure  7f, it is evident that the proposed filter enhanced both the dark and bright regions of the MUAP and suppressed the background pixels. This is because the second eigenvalues of the Hessian matrix for both the dark and bright MUAP regions are high, thus producing a higher filter output. In case of background, both the eigenvalues of Hessian are smaller, thus producing almost a zero-filter output.
To analyze the results of the of the MUAP enhancement algorithms quantitatively, the enhanced image by each algorithm was segmented using the graph cut segmentation. True Positive (TP), True Negative (TN), False Positive (FP), and False Negative (FN) for computer-simulated sEMG signals were computed by labeling the pixels and counting the correctly labeled pixels. For these simulated signals, the skin layer thickness was 1 mm, fat layer thickness was 3 mm, conductivity ratio between skin and fat layers was 20, and conductivity ratio between fat and muscle layers was 0.5. For simulated signals, the ground truth (MUAP location) is already known. On the basis of TP, TN, FP, and FN values, the sensitivity, specificity, and accuracy of these algorithms was calculated for EMG signal, with a signal-to-noise ratio in the range of 50 dB (good-quality EMG) to 5 dB (low-quality EMG). The results of this analysis, for a set of simulated EMG signals, for each algorithm analyzed in this study, are shown in Figure 8a-c. Accuracy is also calculated for another group of sEMG signals of 12 sets, with a signal-to-noise ratio in the range of 5 to 50 dB, as shown in Figure 9a-d From the results, it is evident that the multi-scale filter has the highest accuracy, and has little variations as compare to the other methods, so it should be preferred over other methods. As the Hermite shape filter has almost the same shape as that of the MUAP, it more resembles the MUAP in the sEMG signals and shows good performance at higher noises. However, it has a drawback that it cannot retain the polarity of the MUAP, and also, in case of bidirectional propagation of the sEMG, it creates a distortion to the MUAP and some of the non-propagating noise is also detected as MUAP, which reduces its performance at a low SNR. These results show that the proposed filter can enhance the sEMG images, which will facilitate its further processing for studying muscle anatomy and physiology. Figure 8. (a) Accuracy, (b) specificity, and (c) sensitivity of the proposed multi-scale Hessian-based filter compared with various algorithms analyzed in this study for the enhancement of a set of simulated EMG signals. From the results, it is evident that the multi-scale Hessian filter has a higher accuracy, sensitivity, and specificity for a signal with singal to noise ratio (SNR) greater than 10 dB, so it should be preferred over other methods.

Conclusions
In this study, a novel method based on the Hessian matrix was presented for sEMG spatiotemporal image enhancement. The proposed algorithm was also compared with some of the available digital line-like structure enhancement algorithms using experimental and simulated signals. We proposed multi-scale-based filter enhancement that detects and models the MUAP pattern in spatiotemporal EMG images. Results on various noisy signals showed the precision of the proposed filter for detection and enhancement of the MUAP pattern. The proposed multi-scale filter first modeled the MUAP propagation as a Gaussian line-like structure, and then enhanced it, using the Hessian eigenvalues. The proposed multi-scale filter enhanced both the bright and dark portions of the MUAP propagation structure in the sEMG images, as opposed to the other methods, which can only enhance the bright region. The results for noisy signals at a very low SNR show that the proposed filter is the best method for enhancing the EMG and suppressing the noise.

Conflicts of Interest:
The authors declare no conflict of interest.