Figures
Abstract
Knowledge of the location of muscle Innervation Zones (IZs) is important in many applications, e.g. for minimizing the quantity of injected botulinum toxin for the treatment of spasticity or for deciding on the type of episiotomy during child delivery. Surface EMG (sEMG) can be noninvasively recorded to assess physiological and morphological characteristics of contracting muscles. However, it is not often possible to record signals of high quality. Moreover, muscles could have multiple IZs, which should all be identified. We designed a fully-automatic algorithm based on the enhanced image Graph-Cut segmentation and morphological image processing methods to identify up to five IZs in 60-ms intervals of very-low to moderate quality sEMG signal detected with multi-channel electrodes (20 bipolar channels with Inter Electrode Distance (IED) of 5 mm). An anisotropic multilayered cylinder model was used to simulate 750 sEMG signals with signal-to-noise ratio ranging from -5 to 15 dB (using Gaussian noise) and in each 60-ms signal frame, 1 to 5 IZs were included. The micro- and macro- averaged performance indices were then reported for the proposed IZ detection algorithm. In the micro-averaging procedure, the number of True Positives, False Positives and False Negatives in each frame were summed up to generate cumulative measures. In the macro-averaging, on the other hand, precision and recall were calculated for each frame and their averages are used to determine F1-score. Overall, the micro (macro)-averaged sensitivity, precision and F1-score of the algorithm for IZ channel identification were 82.7% (87.5%), 92.9% (94.0%) and 87.5% (90.6%), respectively. For the correctly identified IZ locations, the average bias error was of 0.02±0.10 IED ratio. Also, the average absolute conduction velocity estimation error was 0.41±0.40 m/s for such frames. The sensitivity analysis including increasing IED and reducing interpolation coefficient for time samples was performed. Meanwhile, the effect of adding power-line interference and using other image interpolation methods on the deterioration of the performance of the proposed algorithm was investigated. The average running time of the proposed algorithm on each 60-ms sEMG frame was 25.5±8.9 (s) on an Intel dual-core 1.83 GHz CPU with 2 GB of RAM. The proposed algorithm correctly and precisely identified multiple IZs in each signal epoch in a wide range of signal quality and is thus a promising new offline tool for electrophysiological studies.
Citation: Marateb HR, Farahi M, Rojas M, Mañanas MA, Farina D (2016) Detection of Multiple Innervation Zones from Multi-Channel Surface EMG Recordings with Low Signal-to-Noise Ratio Using Graph-Cut Segmentation. PLoS ONE 11(12): e0167954. https://doi.org/10.1371/journal.pone.0167954
Editor: William D. Phillips, University of Sydney, AUSTRALIA
Received: July 7, 2016; Accepted: November 23, 2016; Published: December 15, 2016
Copyright: © 2016 Marateb et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All the raw EMG files are available from the figshare database (DOI: 10.6084/m9.figshare.4212123).
Funding: The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the Seventh Framework Programme of the European Union (FP7/2007-2013) under REA grant agreement no. 600388 (TECNIOspring programme), from the Agency for Business Competitiveness of the Government of Catalonia, ACCIÓ, and from Spanish Ministry of Economy and Competitiveness- Spain (project DPI2014-59049-R).
Competing interests: The authors have declared that no competing interests exist.
Introduction
The electromyographic signal (EMG) is the electrical activity associated to contracting muscles. The EMG signal is generated by the summed electrical activity of the muscle fibers stimulated by motoneurons [1]. This signal is recorded either invasively, known as intramuscular EMG, or noninvasively, called surface EMG (sEMG). sEMG is usually applied to asses physiological and morphological characteristics of contracting muscles and their neural strategies [2]. The use of sEMG spans from neurophysiological and medical research (aging, gait and posture analysis), to rehabilitation (biofeedback), ergonomics, sports and movement sciences (biomechanics) [3].
sEMG can be recorded using electrodes of different shapes, sizes, and arrangements [4]. Traditionally, sEMG is recorded as the potential difference between two electrodes placed on the skin, known as bipolar derivation. 1D or 2D electrode arrays including several electrodes along the lines parallel to the muscle fiber orientation are also used [5]. Such electrode systems could be used to investigate the processes of generation, propagation, and extinction of the action potentials in fusiform muscles with fibers parallel to the skin [6]. It is also possible to estimate the muscle fiber Conduction Velocity (CV) and to identify the location of the muscle Innervation Zones (IZs).
Knowledge of the location of muscle IZs in fusiform muscles with fibers parallel to the skin is important for many reasons [6]. The correct estimation of sEMG variables for the monitoring of muscular activity and the detection of muscle fatigue implies electrode positions with knowledge on the IZ location [7, 8]. Moreover, the quantity of injected botulinum toxin for the treatment of spasticity can be minimized if injected close to the IZ [9]. Optimization of neuromuscular electrical stimulation [10], motor point biopsy for diagnosis of neuromuscular diseases [7, 11], and decisions on the type of episiotomy during child delivery [12, 13] are other examples of uses of the information on IZ locations.
Manual identification of the location of muscle IZs by visual analysis [14] is a time-consuming procedure. Several methods have been proposed in the literature to automate this procedure [15–21]. Some of these methods can detect at most one IZ location in each signal interval [15, 17, 18, 20]. There are also various methods to identify multiple IZs. Masuda et al developed a technique based on correlation analysis [16]; the developed program was able to identify up to two IZs in good-quality recordings of biceps brachii muscle. Cescon et al. developed a technique based on Radon transform [21] which could detect multiple IZs in high-SNR sEMG recordings of the gracilis muscles. Mesin et al, proposed a method based on the matched filter [19] which was validated on 20-dB simulated sEMG signals with a maximum of two IZs. All these previous studies have focused on detecting one or two IZs on high-quality sEMG signals. However, it is not often possible to record good quality signals, particularly when using linear arrays of electrodes [3, 22, 23]. In addition, muscles have usually multiple IZs, e.g., the external anal sphincter, brachioradialis, and biceps brachii [3, 16, 24], muscle fibers in partially denervated muscles [25] and fasciculating Motor Units (MUs) in Amyotrophic Lateral Sclerosis (ALS) patients [26] may have more than one IZ as well. The aim of this paper is to automatically identify multiple IZs in medium and low-SNR sEMG signals. A method based on Graph-cut image segmentation and morphological image processing is presented. Part of this work was presented in abstract form [27].
Materials and Methods
Simulation
The model proposed by Farina et al was used to generate surface EMG signals [28]. This model is more complete in comparison with previous approaches [29]. In this model, the volume conductor was described as an anisotropic multilayered cylinder and the source was a spatio-temporal function describing the generation, propagation, and extinction of the intracellular action potential at the end-plate, along the fiber, and at the tendons, respectively. The Inter-Electrode-Distance (IED) was set to 5 mm as recommended in [30] to locate IZs. The remainder of the model parameters used in our study were in principle the same as those used by Mesin et al [19]. Finally, the number of active MUs in each 60-ms simulated signal interval was between 1 and 5. Signals were zero-phase digitally band-pass filtered [31] using an overall eighth-order Butterworth filter with cut-off frequencies 20 and 500 Hz.
For each MU number category (1 to 5), sEMG signals with SNR values of -5, 0, 5, 10 and 15 dB were simulated to include very low to moderate quality sEMG signals. Twenty Single-Differential (SD) channels were simulated along the muscle fiber direction and sampling frequency was 4096 Hz. Thirty frames (or images) with up to 5 IZs were simulated for each SNR value. The temporal location of the IZs was created randomly in each frame. The signal SNR for each simulated 60-ms epoch was defined as the RMS of the raw sEMG divided by the standard deviation of the added Gaussian noise, expressed in dB [32]. Thus, a total of 750 1-D linear array sEMG signals were simulated, considering five SNR values and maximum five MUs.
Methods
Fig 1 shows the flowchart of the proposed IZ detection method. The whole task breaks into 5 main parts. Detailed descriptions will be presented in the following sections; however, each part is briefly introduced as follows:
- Pre-processing: images are generated from spatio-temporal epochs of sEMG, consisting of 1) transfer coefficient calculation to tune the relation between physical units and the pixel data needed to calculate IZ information (such as location and CV) in the time domain, and 2) initializing graph-cut image segmentation algorithm.
- Image segmentation: this decomposes images into the specific regions related to propagations consisting of: 1) Fuzzy clustering (used as a common kernel function) and 2) Graph-cut segmentation.
- Pruning section: a morphological technique is used as a post-processing operation to remove possible outliers.
- Region identification: this extracts required parameters (such as coordinates of the region’s center/edge) as basic information of the IZ detection procedure.
- IZ detection: this locates the IZ location and MU Action Potential (MUAP) CV when the image includes suitable propagating regions in which a propagation swing could be seen.
EMG matrix preparation extracts an appropriate epoch of data for image conversion. Graph-Cut algorithm was used for image segmentation. Parameters (Slope, center/edge coordinates) in step IV were estimated to consider the interaction between regions in the image.
The proposed algorithm was implemented with Matlab® [MathWorks™ Inc. Natick (MA), release 7.12.0].
Image generation from sEMG signals.
A digital image could be defined by a 2-D N by M matrix [33], where 1 ≤ n ≤ N and 1 ≤ m ≤ M are the spatial coordinates of each pixel. The sEMG signal recorded from a linear-array of electrodes, could be considered as an image where N is the number of channels in the recording array, and M is the recorded samples in time. The activity of the muscle (amplitude of the sEMG), was normalized with respect to the maximum value of sEMG signals in the entire recording channels, to represent the intensity of each pixel as a fraction of the maximum value.
Let E be the original N by M sEMG data matrix containing samples sufficiently close so that both the sampling frequencies in space and in time satisfy Nyquist theorem and aliasing is absent in both directions. Interpolation between samples in space and/or in time can then be applied to obtain a new matrix I with dimensions αN and βM, where α and β are interpolation coefficients. Cubic spline interpolation [34] by a factor of 10 (α = β = 10) was used for both dimensions.
In fusiform muscles with fibers approximately parallel to the skin, the amplitude of bipolar EMG signal is relatively low near the IZ locations because of potential cancellations. Propagation of MUAPs occurs from the IZ toward the muscle-tendon junctions. Thus, it creates special patterns in the corresponding regions in the image (Fig 2). By extracting such regions, it is possible to detect the IZ location as well as the MUAP conduction velocity and the direction of propagation.
(a) The image was generated from 60-ms epochs of linear-electrode arrays (spatial interpolation factor of 100 for visualization and the inter-electrode-distance of 5 mm) with the sampling frequency of 4096 Hz and (b) the result of image segmentation via kernel-based Graph-Cut. Propagated potentials from each IZ towards both directions (up and down in the Figure) can be seen.
In this approach, the location of IZ’s was obtained in pixel units. However it is necessary to report the results (IZ and MUAP CV) in the space and time domain. Thus a transfer coefficient was required from E to I. The calibration was performed by the algorithm to report the results in the suitable domain.
Image segmentation using kernel Graph-Cut.
Many problems in computer vision such as segmentation, can be put into labeling problems where an undirected graph is used for data representation. It is in fact a summary representation of locations and their neighborhood structure. The solution of the labeling problem is the assignment of a label to each vertex in the graph based on the criteria in the problem’s requirements. Such criteria are represented with the energy, a cost function evaluating how good the given labeling is, where the smaller energy the better solution. Thus, the labeling problem is considered as an energy minimization problem. This function usually has many local minimums. Using graph cuts, it is possible to find the minimum of such energy functions [35]. The information about the graph-cut image segmentation method used in our study, was provided in S1 Appendix.
Let I: p ∈ Ω ⊂ R2 → Ip = I(P) ∈ Ψ be the image in the R2 space [36] defined by Ω transferred to Ψ feature space. For segmenting I into maximum G regions, each pixel must be labeled as l ≤ G from the finite set of labels Γ. Each region is defined as a set of pixels that have the same label. Now, λ indexing function is introduced to assign a label from the finite Γ set to each point in the image as:
(1)
In computer vision, the segmentation is usually performed via minimizing an energy function (F) [37, 38]. The energy function (F) can then be written as:
(2)
where D is the data energy (that encodes the constraints of the data by measuring the disagreement between the segment representative intensity and the intensity of the pixels in each segment i.e. similar to the objective or cost function of a classic clustering algorithm), R is the prior energy (that encodes the constraints provided by prior knowledge about the optimal labeling by measuring the extent to which labeling is not piece-wise smooth i.e. extracted segments must be smooth) and the constant γ is a positive regularized coefficient that controls relative importance of data and prior term [39]. This form of objective function is also used in Regularized Least Squares (RLS), where the first term is supposed to minimize. However, the knowledge about the possible solutions are also expressed as a penalized problem in which a regularized function (the second term) is added to the objective function [40].
In our study the regularized coefficient was adaptively estimated based on the quality of the image frame (Quality factor: QF) estimated by the method proposed by Zhang et al [41]. We fitted a linear regression model between QF (input) and the regularized coefficient (output) in a set of image frames with variety of SNR values. Such frames were not used for validation. The regression equation was output = 0.801×QF+0.888 (R2 = 0.8). The regularized coefficient ranged from 0.226 to 0.639, in our study.
According to the maximum a posterior (MAP) formulation and piecewise constant segmentation model, the data term could be defined as follows [36, 39, 42]:
(3)
where μl is the piecewise constant model of Rl (i.e. the segment representative intensity). The term DP(λ(P)) measures how much assigning λ(P) to pixel P disagrees with the original pixels’ intensity. The prior energy term expressing the smoothness constraints is defined as follows:
(4)
The entire neighboring pairs {p,q} are stored in a neighborhood set H. The term r is the smoothness regularization or neighbor interaction function, gives penalties to neighboring pixels when they have different labels. Thus the smoothness energy prior is the sum of interaction value between every pair in the neighborhood set [39]. Briefly, it measures how different the intensities of pixels are from each other in each segment.
The smoothness prior function could be obtained by the truncated squared absolute difference principle i.e. min(const2,|μλ(p) − μλ(q)|2) [42]. Unlike the other functions such as mean square error and mean absolute difference, the truncated quadratic function is robust and thus it limits the influence of mismatches. This assigns higher value to the larger difference between labels and vice versa. The “const” term is a constant value that is used to make the function more robust and prevent possible outliers [39, 42, 43]. In the pre-processing procedure, the weight of smoothness constraint was estimated (Graph-Cut initialization). Accordingly, any abrupt changes of the image were detected and the corresponding smoothness weights were calculated.
In image segmentation, an objective function measures the goodness of the possible solutions (i.e. the label of pixels) according to a set of constrains (i.e. the prior term). The lower the value of the objective function, the better goodness of fit is reached. Thus the global minimum of the objective function gives the optimized solution[39]. In order to minimize the objective functionF(λ), it is useful to transform image data to the kernel space via a mapping function φ and to write the segmentation function in the kernel space. In machine learning, the ‘kernel trick’ is using a linear classifier to solve a nonlinear problem by mapping the original nonlinear data into a higher dimensional space. It is particularly useful to consider the complexity of image data and insufficiency of computationally efficient models such as piecewise Gaussian to partition nonlinear separable data [36]. Meanwhile, data are not handled directly and the dot product in the kernel space suffices to write the kernel induced function based on the Mercer’s theorem. The kernel function must be continuous, symmetric, and positive semidefinite. In fact, the function φ is not needed to be computed explicitly [36, 44, 45]. The kernel induced segmentation function could be defined as the following:
(5)
where JK is non-Euclidean distance measure defined as:
(6)
According to the mercer’s theorem [45], the kernel function K(x,y) can be written as:
(7)
In our study, the Radial basis Function (RBF) was used as the kernel function (K). Accordingly, Eq 2 was transformed into Eq 5, in which JK was calculated using Eq 6. Solving image segmentation in a kernel-induced space with graph cuts consists of finding the labeling to minimize Eq 5. In fact, the kernel induced segmentation function is used via the following optimization strategy:
We can optimize the kernel segmentation function in two steps. First, we use RBF kernel function to fix some points as region centers, by using the Fuzzy C-Mean (FCM) as to initiate the segmentation. Second, the locations of the points are optimized and upgraded by the graph-cut technique [36]. The result of such segmentations on an epoch of EMG signal is shown in Fig 2.
Pruning and region identification.
Removing outliers (i.e. segments not related to propagation) is necessary prior to extracting features from the extracted regions. It is thus necessary to separate relevant regions and to remove overlapped/irrelevant regions. The pruning procedure consists of four steps performed by morphological processing. The detailed tutorial about the morphological image processing techniques used in our study, was provided in S2 Appendix.
Morphological processing is used as a tool for extracting image components useful in a variety of machine learning problems, such as pruning [33]. We used several morphological techniques, such as dilation, by which the discontinuity of the extracted regions is compensated. The breakpoints i.e. the areas in which the propagation pattern is strongly influenced by noise are recovered via the pair-point structure (the logical OR operation between the pair structures with the offset values of [-3,3] and [3,3]) (Fig 3). Opening is used to separate the overlapped regions in order to extract symmetric regions via disk structure (the radius of 6 pixels). Then, Erosion operation is used to suppress any irrelevant regions not related to symmetric propagation via the line structure in both directions of the propagation (with the length of 20 pixels and degrees of 45° and -45° for the upper and lower V-shaped propagations, respectively). The erosion of the image was performed using these two line elements, in parallel, and the logical OR operator was used to combine these two images. Finally, such regions are reconstructed using 4-connectivity morphological reconstruction where the marker was the opened image and the mask was the resulting image after erosion. These procedures are shown in Fig 4. The extracted regions were classified as unsuccessful (i.e. no clear propagation), or successful at this step. Candidate regions are then processed at the next step. The MUAP propagation has positive and negative swings (bright and dark portions). The image segmentation is performed on both swings. For unsuccessful regions in each swing, the other swing is also analyzed and the results are then combined when successful regions are formed. The parameters of morphological structuring elements were tuned using trial and error. For the detailed information about the morphological image processing, the reader is invited to read the review on mathematical morphology by Heijmans [46].
(a) the segmented image with discontinuous regions (some breakpoints are shown in rectangular area). (b) The compensated image with fewer breakpoints and intensified regions.
In such procedures, the irrelevant structures i.e. non-propagating parts were suppressed while the detectability of the propagating regions close to the innervation zones was improved.
IZ detection and feature extraction.
The images generated in the previous section (Fig 4d), are directly related to the IZ and MUAP propagations. The parameters slope and center/edge coordinate are extracted for each region to perform paired-region labeling and offset calculation (Fig 5). Regions not paired in the algorithm, were excluded from further analysis. The offset of paired-region is calculated by averaging two edge coordinates of the regions with the closest distance.
The slope parameters found by center/edge coordinates are shown by triangles and pentagons, respectively. Bold triangles show the closest distance of edges. The center is defined as the center of each propagating region. The edges are the upper and lower boundaries of such regions. The slope is calculated based on the angle between a virtual line representing the propagation region (bold line) and the horizontal line.
The IZ detection procedure is then performed in three steps. First, labels are assigned to each “paired-region” defined as the symmetric opposite propagation. This region is then represented with a line (Fig 6a). Second, the offset of the paired-regions and intersection of lines are calculated (Fig 6b). The intersection points and offsets of nearby paired-regions with minimum distance are linked together. This operation iterates until all points are checked. (Fig 6c; the labels with min. distances are shown by pentagon stars). Finally, MUAP conduction velocity is estimated based on the slope of the extracted propagating areas.
The feature extraction procedure including (a) assigning labels to the paired-regions (paired regions are shown with the same labels). Each region is represented by a straight line to calculate the intersection points. (b) Calculation of the offset of the paired-region and intersection of lines, offsets and intersection points are shown by hexagons and stars inside black circles respectively. (c) The estimation of the innervation zone by choosing appropriate hexagon according to its distance with corresponding intersection point. The detected innervation zones are shown in pentagon stars.
Validation
Each 60-msec sEMG frame was automatically analyzed by the proposed algorithm. In each frame, if one IZ was found within ±0.5 IED (and ±0.5 msec) of a simulated IZ, it was counted as a True Positive (TP). When no IZ was identified within the above threshold, it was considered as a missed IZ (False Negative: FN). Those detected IZs in a frame not considered as TPs, were counted as False Positives (FPs). No post-processing was performed on the propagating areas to reject FPs if the CV was not within the normal range [47, 48]. In each SNR category, the total TPs, FNs and FPs were calculated in the entire 150 signals and the performance of the proposed method was assessed. The following performance measures are reported: Sensitivity (= Recall) (Se) as the proportion of simulated IZs correctly identified, Precision (= PPV) (Pr) as the proportion of identified IZs that are correct, % Missed IZs, % Erroneous IZs, and F1-score, defined as micro- or macro-averaging measures[49], defined as:
(8)
(9)
(10)
(11)
(12)
In the micro-averaging procedure, the number of TPs, FPs and FNs in each frame are summed up to generate cumulative measures and precision, recall and F1-score are then evaluated. In the macro-averaging, on the other hand, precision and recall are calculated for each frame and their averages are used to determine F1-score. For each correctly identified IZ, the location error measured as the absolute channel error (IED ratio), and the absolute and relative muscle fiber conduction velocity error are also reported.
To examine the effect of spatial distribution of the IZ locations on the performance of the proposed algorithm, the Average Nearest Neighbor (ANN) distance [50] was calculated in the image space for different distributions of the IZ channels. ANN was defined for each analyzed image frame as follows:
(13)
where n is the total number of IZ locations simulated in that frame, A is the total number of pixels in the image frame with the size of αN by βM (i.e.αβ×NM), and di is the Euclidean distance between the IZ location i and the closest neighboring IZ location in pixels. The lower the ANN value, the closer the IZ locations. In our study, the values of parameters M and N are (60 ms × 4096; rounded) and 20, respectively. Parameters α and β are identical to 10. Thus, the value of A is 492000.
On the other hand, to identify whether the closeness of the IZ locations to the frame center could affect the performance of the proposed algorithm, the spatial feature Center Closeness (CC) in the range of [0,1] was defined as follows:
(14)
where di is the Euclidean distance between the IZ location i and the center of the image frame in pixels. The higher the CC value, the closer the IZ locations to the frame borders.
Statistical analysis
We aimed to identify whether the complexity of the recorded signal and the spatial information of the IZ locations could affect the performance of the proposed IZ detection algorithm. Thus, the Poisson regression model [51] was used to determine the significance of the relationship between the count data number of missed (FNs) and erroneous (FPs) (dependent variable) and the number of IZ locations in each frame, frame SNR, and spatial parameters ANN and CC (independent variables). The level of significance was set to P<0.05. Data was analyzed using STATA 10 [52].
Results
Some of simulated sEMG signals are shown in Fig 7 as examples. Image frames A to D contained 2 IZs (-5 dB SNR), 3 IZs (0 dB), 4 IZs (5 dB), and 5 IZs (10 dB), respectively. The proposed algorithm could accurately locate all of the IZ locations with the average channel error of 0.07±0.05 IED.
The image frames A to D contained 2 IZs (-5 dB SNR), 3 IZs (0 dB), 4 IZs (5 dB) and 5 IZs (10 dB), respectively. The location of the simulated IZs is shown by circles. The developed program automatically identified the location of IZs as the crossing of the ‘v’ shape propagation lines (upper lines in blue and lower lines in red color). The CV of the identified propagation pattern was then estimated by the proposed algorithm.
The spatial distribution of the simulated IZ locations is shown in Table 1 in terms of ANN and CC parameters for each IZ number category. Note that ANN parameter is not defined for frames with 1 IZ.
The average number of missed (FN), erroneous (FP) and correctly identified (TP) IZs in each frame is shown in Table 2 for each number of IZs and SNR category.
The overall micro-averaged performance of the proposed IZ location detection algorithm for each SNR category is shown in Table 3. Overall, sensitivity, precision and F1-score of the algorithm on the sEMG frames with positive SNR were 92.5%, 98.3% and 95.3%, respectively. Also, the macro-averaged sensitivity, precision and F1-score of the proposed algorithm on the entire frames were 87.5%, 94.0% and 90.6%, respectively.
For the correctly identified IZ locations, the bias error for the innervation zone detection algorithm for each IZ number and SNR category is listed in Table 4. Overall the average bias for IZ channel selection was 0.02±0.10 IED ratio.
The absolute and relative absolute muscle fiber CV error of the proposed algorithm for each IZ number and SNR category is listed in Table 5. Overall the average absolute and relative CV estimation errors were 0.41±0.40 and 11±10%, respectively.
The computational complexity of the proposed algorithm was assessed in terms of the program running time on each image frame in average for each IZ number category (Table 6). Overall, the average running time of the proposed algorithm on each 60-ms sEMG frame was 25.5±8.9 (s). The program was analyzed on an Intel dual-core 1.83 GHz CPU with 2 GB of RAM.
The regression models produced moderate fits (Pearson dispersion ratio ≈ 0.8 in both FN and FP models). The regression results showed that, all the considered parameters except CC were significant in the FN model (i.e. the number of missed IZs) while the SNR category and the number of IZs in the frame were significant in the FP model (i.e. the number of erroneous IZs) (p < 0.05). Thus, whether the IZ locations are closer to the frame border or not, does not affect the performance of the algorithm.
Discussion
Settings of the algorithm
For the series-fibered muscles, such as brachioradialis, up to 6 distinct endplate zones (4 zones on average) have been observed [24]. In a study performed by Cescon et al, about 5 to 6 different IZs were identified and further analyzed in the external anal sphincter muscles [12]. Thus, an automatic detection algorithm should address the issue of multiple IZs (Fig 7).
Among different image segmentation algorithms, Graph-cut was used in our study. It was used in low-level vision problems in the literature [53]. There is a link between this method and diffusion and it enforces piecewise smoothness while preserving relevant sharp discontinuities [54]. Watershed image segmentation, on the other hand, was used for segmentation of sEMG images to improve muscle activity estimation [55]. It is possible to combine these two techniques to improve image segmentation [56, 57], which is the focus of our future activity.
In morphological image processing steps of our algorithm, dilation was first used to fill the gaps between propagation regions. The pair-point structure was used in this step, since it reduces the probability of hits. Thus, the possibility of merging propagation regions of different IZs is reduced. Meanwhile, the detected propagation lines must span within three recording electrodes to accept the related propagation region in our algorithm.
Utilizing the IED of 5 mm in our simulations, which was also used in similar studies [16, 18], the spatial aliasing is negligible [30, 58]. In our study, the sampling frequency of 4096 Hz was used to simulate sEMG signals. Then, an interpolation factor of 10 was used to obtain the image frame. Thus, another sampling frequency could be used changing the factor accordingly.
Twenty SD channels were simulated in each array in our study. The detection algorithm does not require a long propagation line (Fig 7A) and therefore many channels. On the other hand, the number of channels influences the number of True Negatives (TNs), defined as the number of recording channels without any IZs that were not identified as IZ channels by the algorithm in each image frame. Thus, the greater the number of channels is, the greater the number of TNs. Moreover, the number of TNs is usually much higher than that of FNs or FPs in each frame. This is the reason to use performance indexes which do not depend on TNs (Eqs 8–12). They are objective performance measures in the unbalanced datasets [49]. Conversely, accuracy, which is also defined in information theory and is not used in our study, overestimates the performance of the algorithm tending to 100% with a very high number of TNs.
Finally, the SD derivation was used because it is commonly applied to detect IZs [16, 17]. In SD EMG signals, propagating potentials along the electrode array are enhanced while the non-propagating components related to power-line-interference and end-of-fiber effect are reduced[3]. Also, it is known that monopolar EMG signals contain so-called far-field potentials [59].
Detection performance
The proposed algorithm could perfectly identify up to two IZs in frames with the whole range of SNR (Table 2). For frames with more IZs, the performance of the detection algorithm was very good under medium levels of SNR (Tables 2 and 3). The average IZ channel bias error was less than 1 mm (0.2 IED) in more than 90% of cases and only in less than 5% of the frames, the error was slightly greater than 1.5 mm (0.32 IED).
The morphological Image processing technique used in our algorithm relies on the slope of propagation. Thus, it cannot be used in pathological cases (such as ALS [26]) where the slope of propagation is close to 90°. Finally, in fusiform muscles with multiple IZs, the CV cannot be usually reliably estimated [3]. However, in the proposed algorithm the average absolute CV estimation error was only 0.41 m/s which is in the range of the expected error in real signals [60].
Comparison with state-of-the-art
The performance of the designed algorithm was compared with that of other methods, proposed in the literature [15–17]. Such methods identified the IZ channel in the SD signal array based on the lowest RMS amplitude or highest mean frequency values (MNF), or the lowest correlation coefficient among the channels with opposite propagations.
Although such methods required a rather long stationary sEMG frame to estimate the time or frequency domain features (e.g. 250 msec) [61], much shorter epochs (e.g., 4 msec) were used to identify the IZ channels in signals with multiple IZs [16]. This epoch length is slightly less than the action potential duration in sEMG signals. Such epochs were analyzed based on the RMS and correlation analysis to find the IZ channels in the array. However, due to the low frequency resolution of the analyzed 60-msec epoch and also its lower accuracy in comparison with other methods [15], MNF criterion was not used for comparison. The value of (micro) F1-score measure was 17% and 25% for the RMS and correlation analysis on the positive SNR data, respectively compared with what obtained from the proposed algorithm (F1-score = 95%).
The matched filter method proposed by Mesin et al [19] had an average error of 0.3 IED (ratio) reaching values higher than 1 IED for the simulated sEMG signals with up to two IZs in each frame and SNR value of 20 dB [20]. This method is in fact the first step (i.e., template matching) of the sEMG decomposition method proposed by Gazzoni et al [62]. Since the approach proposed by Gazzoni et al relied on MUAP segmentation in each recording channel, it failed in not only very high-SNR data to some extent [20], but in low-SNR data more frequently. Image processing methods are preferred because no detection threshold is required for the MUAPs and the whole image frame is analyzed, compared with each signal.
Östlund et al, developed an image processing method based on Particle Image Velocimetry (PIV) to identify the location of muscle IZ [18]. Although a wide SNR range of -15 to 30 dB was analyzed in their study, only one IZ channel was simulated and then identified. The average channel error was about 0.60 IED ratio for the PIV method. Our algorithm correctly identified all of IZs in similar cases (Table 2; one IZ category) and the error was 0.13 IED ratio.
Ullah et el, developed an image processing method based on bi-dimensional cross correlation between the interpolated image frame and its flipped transformation [20]. That method was applied on the MUAP templates obtained by spike-triggered averaging, thus simulated sEMG frames had a very high SNR, no noise was added, and only one IZ using the results of sEMG decomposition was considered. [20]. Although this method was virtually bias free (average IZ channel error of 1.3% IED) compared with that of our method i.e. 1.8%, the inter-quartile-ratio (IQR) was 0.23 IED ratio whereas this was only 0.18 IED in our study for our entire dataset (from one to 5 IZs, SNR range of -5 to 15 dB) considering not only one IZ or very high SNR.
Soares et al proposed an algorithm based on morphological image processing to estimate the CV of the sEMG frames evaluated on the simulated signals with the SNR values of 12, 16, 20, 30, and ∞ (i.e. no noise) in cases with only one IZ channel yielding an average RMS error (RMSE) for the CV of 0.07 m/s [63]. In addition, Farina et al proposed an approach based on the modified multi-channel maximum likelihood estimation and beamforming for sEMG frames with only one IZ [60], yielded an average CV RMSE of 0.13 m/s on 12 dB SNR simulated signals[63].
On the other hand, the CV RMSE in our algorithm was 0.20 m/s for sEMG frames including one IZ with positive SNR (5, 10 and 15 dB), that is, with lower SNR considered in the previous studies.
Additional analysis
The performance of the proposed IZ detection algorithm was further assessed considering different simulation conditions and analysis methods. In each analysis, only a simulation condition or analysis method was changed in comparison with those of what presented in this manuscript (baseline) and the performance of the algorithm was presented.
An IED of 5 mm was used in our study. However, arrays with IED of 10 mm are often used in practice. Thus, sEMG signals with 10 mm IED were simulated and the performance of the proposed algorithm was further assessed (Table 7; Scenario 1). The performance is rather lower than that of the baseline that could be because of the spatial aliasing present in sEMG signals with 10 mm IED but not 5 mm IED [30, 58].
In our method, the interpolation factor of 10 was used for the time dimension (β = 10). Also, the sampling frequency of 4096 Hz was used. The performance of the algorithm was assessed when no interpolation was performed in the time dimension (β = 1) (Table 7; Scenario 2). The performance of the algorithm was deteriorated in comparison with that of baseline. Thus, it is necessary to interpolate sEMG signals in the time dimension.
In our study, cubic spline interpolation was used to increase the spatial and temporal resolution of the sEMG image frames. The performance of the proposed IZ detection method was further assessed using another image interpolation method, namely as sinc-based (Table 7; Scenario 3). From sampling theory, it follows that the ideal interpolation kernel is the sinc function, i.e. ideal low-pass filter, which is of infinite extent (unbounded support). In practical applications its support is limited by using truncation. However, there is excessive ripple in the pass band and poor attenuation in the stop band (c.f. Gibbs effect) [64],resulting in poor performance in image interpolation [65]. Windowing is used to reduce the abruptness of the truncated ends and thus improve the frequency response. Among windowed-sinc filters, Lanczos filter with 2 spatial grid points was shown to be best compromise in terms of reduction of aliasing, sharpness, and minimal ringing[66].
Cubic spline, on the other hand, has only positive values which is very interesting for image processing to guarantee positive interpolation values. It is a linear, shift-invariant filter [67]. In our study, the performance of the cubic spline was better than that of Lanczos-2 filter. One reason could be that cubic spline has a favorable stop-band response [68]. It was also indicated in the literature that it is sufficient for several practical applications [69] or conditions e.g. sampling at Nyquist rate [70].
Gaussian noise is not the only disturbance present in recorded sEMG signals. Power line interferences were further simulated and added to the signal based on the method proposed by Allen [71].
The power-line interference was simulated at the fundamental power-line frequency (f1 = 50 Hz) and all of its harmonics up to f10 = 500 Hz. The frequency, amplitude and phases of each interference signal for a given sEMG were uniformly distributed to model the interference variations. The simulated main power line frequency was within 49.5–50.5 Hz. The phase distribution for the entire interference frequencies was within ±pi/2 radians. The magnitude of the fundamental frequency was determined randomly in the range 0–50% of the sEMG Root-Mean-Square (RMS) amplitude. Such ratios ranged from 25% to 50%, and 12% to 25% for frequencies f2, and f3, respectively. For the other frequencies (f4 to f10), it ranged from 6% to 12%. The Gaussian noise was then added to have the SNR range of -5 dB to 10 dB. The performance of the proposed algorithm was assessed (Table 7; Scenario 4). Its performance is rather lower than that of the baseline method. However, in practice, the power line interferences are significantly reduced when recording in SD mode. In our simulations, such interferences were directly added to the SD signals to resemble poor recording conditions. Some of simulated sEMG signals using power-line inferences in addition to Gaussian noise are shown in Fig 8 as examples. Image frames A to D contained 2 IZs (7.2 dB average SNR), 3 IZs (8.5 dB), 4 IZs (7.5 dB), and 5 IZs (6.7 dB), respectively. The proposed algorithm could accurately locate all of the IZ locations with the average channel error of 0.08±0.07 IED.
The image frames A to D contained 2 IZs (-5 dB SNR), 3 IZs (0 dB), 4 IZs (5 dB) and 5 IZs (10 dB), respectively. The location of the simulated IZs is shown by circles. The developed program automatically identified the location of IZs as the crossing of the ‘v’ shape propagation lines (upper lines in blue and lower lines in red color). The CV of the identified propagation pattern was then estimated by the proposed algorithm.
Practically, it is also possible to have intermittent contacts, whose sEMG signals are known as ‘outliers’. It is necessary to identify outliers and reconstruct sEMG time samples, prior to using our proposed IZ detection algorithm. The algorithms proposed by Marateb et al [23] and Ghaderi and Marateb [72] could be used. However, they require recording HDsEMG signals using 2D electrode arrays. Our algorithm could then be used for each electrode column.
Final considerations
The type of muscle is an important factor for the IZ arrangement [18]. In fact, human muscles could have single, multiple and scattered IZs [73]. The proposed detection algorithm can be used in muscles of the first two types, since the locations of IZs are identified based on the crossing of the isolated propagation lines. In our study, simulated sEMG signals were used to evaluate the performance of the algorithm and a large variety of signal complexities was considered. Such findings could not be easily generalized to recorded sEMG signals because such signals lack the common sense of actual EMG signals [74]. However, the simulation parameters were set as to resemble fair recording conditions. Moreover, the proposed IZ detection algorithm could be used for sEMG signals with low to moderate complexity where isolated MUAPs could be seen. It requires bi-directional propagation for morphological signal processing. Also, our proposed algorithm in its current form is applicable for sEMG signals recorded with 1D electrode arrays. Using HDsEMG (High Density sEMG) signals, recorded using 2D electrode arrays, implies combining 2D spatial information and time propagation which is the focus of our future work.
Conclusions
In conclusion, we developed a fully-automatic offline algorithm to identify multiple IZs in low-SNR multi-channel sEMG signals. In addition, CV was also estimated. This algorithm is robust and accurate and thus a promising new tool for non-invasive IZ identification. Future work will focus on the application of the algorithm to experimental signals.
Supporting Information
S1 Appendix. Image segmentation using kernel Graph-Cut.
https://doi.org/10.1371/journal.pone.0167954.s001
(PDF)
Acknowledgments
The authors would like to thank Roberto Merletti for reviewing a draft of this paper and Marjan Mansourian for valuable discussions on the statistical methods.
Author Contributions
- Conceptualization: HRM MF.
- Data curation: HRM MF.
- Formal analysis: HRM MF MR MAM.
- Funding acquisition: HRM MAM.
- Investigation: HRM MF DF.
- Methodology: HRM MF DF.
- Project administration: HRM MAM.
- Resources: DF.
- Software: HRM MF.
- Supervision: MAM.
- Validation: HRM MR MAM DF.
- Visualization: MF.
- Writing – original draft: HRM MF.
- Writing – review & editing: HRM MF MR MAM DF.
References
- 1. De Luca CJ. Physiology and mathematics of myoelectric signals. IEEE transactions on bio-medical engineering. 1979;26(6):313–25. pmid:468280
- 2. Farina D, Merletti R, Enoka RM. The extraction of neural strategies from the surface EMG: an update. Journal of applied physiology (Bethesda, Md: 1985). 2014;117(11):1215–30.
- 3.
Barbero M, Merletti R, Rainoldi A. Atlas of muscle innervation zones: understanding surface electromyography and its applications: Springer Science & Business Media; 2012.
- 4. Hogrel JY. Clinical applications of surface electromyography in neuromuscular disorders. Clinical neurophysiology. 2005;35(2–3):59–71. pmid:16087069
- 5. Zwarts MJ, Stegeman DF. Multichannel surface EMG: basic aspects and clinical utility. Muscle & nerve. 2003;28(1):1–17.
- 6. Merletti R, Farina D, Gazzoni M. The linear electrode array: a useful tool with many applications. Journal of Electromyography and Kinesiology. 2003;13(1):37–47. pmid:12488085
- 7. Saitou K, Masuda T, Michikami D, Kojima R, Okada M. INNERVATION ZONES OF THE UPPER AND LOWER LIMB MUSCLES ESTIMATED BY USING MULTICHANNEL SURFACE EMG. Journal of Human Ergology. 2000;29(1–2):35–52. pmid:12696320
- 8. Merletti R, Rainoldi A, Farina D. Surface electromyography for noninvasive characterization of muscle. Exercise and sport sciences reviews. 2001;29(1):20–5. pmid:11210442
- 9. Lapatki BG, van Dijk JP, van de Warrenburg BPC, Zwarts MJ. Botulinum toxin has an increased effect when targeted toward the muscle's endplate zone: A high-density surface EMG guided study. Clinical Neurophysiology. 122(8):1611–6. pmid:21195024
- 10. Gobbo M, Maffiuletti NA, Orizio C, Minetto MA. Muscle motor point identification is essential for optimizing neuromuscular electrical stimulation use. Journal of neuroengineering and rehabilitation. 2014;11(1):1.
- 11. Brown R, Dissanayake KN, Skehel PA, Ribchester RR. Endomicroscopy and electromyography of neuromuscular junctions in situ. Annals of clinical and translational neurology. 2014;1(11):867–83. pmid:25540801
- 12. Cescon C, Raimondi EE, Zacesta V, Drusany-Staric K, Martsidis K, Merletti R. Characterization of the motor units of the external anal sphincter in pregnant women with multichannel surface EMG. International urogynecology journal. 2014;25(8):1097–103. pmid:24623258
- 13. Cescon C, Riva D, Zacesta V, Drusany-Staric K, Martsidis K, Protsepko O, et al. Effect of vaginal delivery on the external anal sphincter muscle innervation pattern evaluated by multichannel surface EMG: results of the multicentre study TASI-2. International urogynecology journal. 2014;25(11):1491–9. pmid:24687366
- 14. Gallina A, Merletti R, Gazzoni M. Innervation zone of the vastus medialis muscle: position and effect on surface EMG variables. Physiological Measurement. 2013;34(11):1411. pmid:24081116
- 15. Beck TW, DeFreitas JM, Stock MS. Accuracy of three different techniques for automatically estimating innervation zone location. Computer Methods and Programs in Biomedicine. 2012;105(1):13–21. pmid:20692064
- 16. Masuda T, Miyano H, Sadoyama T. The Position of Innervation Zones in the Biceps Brachii Investigated by Surface Electromyography. Biomedical Engineering, IEEE Transactions on. 1985;BME-32(1):36–42.
- 17. Rainoldi A, Melchiorri G, Caruso I. A method for positioning electrodes during surface EMG recordings in lower limb muscles. Journal of Neuroscience Methods. 2004;134(1):37–43. pmid:15102501
- 18. Östlund N, Gerdle B, Stefan Karlsson J. Location of innervation zone determined with multichannel surface electromyography using an optical flow technique. Journal of Electromyography and Kinesiology. 2007;17(5):549–55. pmid:16890457
- 19. Mesin L, Gazzoni M, Merletti R. Automatic localisation of innervation zones: A simulation study of the external anal sphincter. Journal of Electromyography and Kinesiology. 2009;19(6):e413–e21. pmid:19269857
- 20. Ullah K, Cescon C, Afsharipour B, Merletti R. Automatic detection of motor unit innervation zones of the external anal sphincter by multichannel surface EMG. Journal of electromyography and kinesiology: official journal of the International Society of Electrophysiological Kinesiology. 2014;24(6):860–7.
- 21.
Cescon C. Automatic Location Of Muscle Innervation Zones From Multi-Channel Surface EMG Signals. Medical Measurement and Applications IEEE International Workshop on; Benevento, Italy2006. p. 87–90.
- 22. Chowdhury RH, Reaz MBI, Ali MABM, Bakar AAA, Chellappan K, Chang TG. Surface Electromyography Signal Processing and Classification Techniques. Sensors (Basel, Switzerland). 2013;13(9):12431–66.
- 23. Marateb HR, Rojas-Martínez M, Mansourian M, Merletti R, Mañanas Villanueva MA. Outlier detection in high-density surface electromyographic signals. Medical & Biological Engineering & Computing. 2011;50(1):79–89.
- 24. Lateva ZC, McGill KC, Johanson ME. The innervation and organization of motor units in a series-fibered human muscle: the brachioradialis. Journal of Applied Physiology. 2010;108(6):1530–41. pmid:20360433
- 25.
MacIntosh BR, Gardiner PF, McComas AJ. Skeletal muscle: form and function. 2nd ed. Champaign, IL: Human Kinetics; 2006. viii, 423 p. p.
- 26. Jahanmiri-Nezhad F, Barkhaus PE, Rymer WZ, Zhou P. Innervation zones of fasciculating motor units: observations by a linear electrode array. Frontiers in human neuroscience. 2015;9:239. pmid:26029076
- 27.
Farahi M, Marateb HR, Mansourian M, Muceli S, Farina D. DETECTION OF MULTI-INNERVATION ZONES FROM HIGH-DENSITY SURFACE EMG USING THE GRAPH-CUT SEGMENTATION. In: Felici F, editor. XX conference of the international society of electrophysiology and kinesiology; Rome: ISEK; 2014. p. 432.
- 28. Farina D, Mesin L, Martina S, Merletti R. A surface EMG generation model with multilayer cylindrical description of the volume conductor. IEEE transactions on bio-medical engineering. 2004;51(3):415–26. pmid:15000373
- 29. Farina D, Merletti R. A novel approach for precise simulation of the EMG signal detected by surface electrodes. IEEE transactions on bio-medical engineering. 2001;48(6):637–46. pmid:11396594
- 30.
Afsharipour B, Ullah K, Merletti R. Spatial Aliasing and EMG Amplitude in Time and Space: Simulated Action Potential Maps. In: Roa Romero ML, editor. XIII Mediterranean Conference on Medical and Biological Engineering and Computing 2013: MEDICON 2013, 25–28 September 2013, Seville, Spain. Cham: Springer International Publishing; 2014. p. 293–6.
- 31. Gustafsson F. Determining the initial states in forward-backward filtering. IEEE Transactions on Signal Processing. 1996;44(4):988–92.
- 32.
Kay SM. Fundamentals of statistical signal processing. Englewood Cliffs, N.J.: Prentice-Hall PTR; 1993.
- 33.
Gonzalez RC, Woods RE. Digital Image Processing (2ed Edition). Secend ed. New jersey: Prentice-Hall, Inc.; 2001. 793 p.
- 34. Hsieh H, Andrews H. Cubic splines for image interpolation and digital filtering. IEEE Transactions on Acoustics, Speech, and Signal Processing. 1978;26(6):508–17.
- 35.
Lézoray O, Grady L. Image processing and analysis with graphs: theory and practice. Boca Raton, FL: Taylor & Francis; 2012. xxiv, 537 p. p.
- 36. Salah MB, Mitiche A, Ayed IB. Multiregion Image Segmentation by Parametric Kernel Graph Cuts. Image Processing, IEEE Transactions on. 2011;20(2):545–57.
- 37. Hofmann T, Puzicha J, Buhmann JM. Unsupervised texture segmentation in a deterministic annealing framework. Pattern Analysis and Machine Intelligence, IEEE Transactions on. 1998;20(8):803–18.
- 38. Poggio T, Torre V, Koch C. Computational vision and regularization theory. Nature. 1985;317(6035):314–9. pmid:2413361
- 39.
Veksler O. Efficient Graph-Based Energy Minimization Methods in Computer Vision: Cornell University; 1999.
- 40.
Beck A. Introduction to nonlinear optimization: theory, algorithms, and applications with MATLAB: SIAM; 2014. xii, 282 pages p.
- 41. Zhang L, Zhang L, Mou X, Zhang D. FSIM: a feature similarity index for image quality assessment. IEEE transactions on Image Processing. 2011;20(8):2378–86. pmid:21292594
- 42. Boykov Y, Veksler O, Zabih R. Fast approximate energy minimization via graph cuts. Pattern Analysis and Machine Intelligence, IEEE Transactions on. 2001;23(11):1222–39.
- 43.
Salah M, Mitiche A, Ayed I. A Continuous Labeling for Multiphase Graph Cut Image Partitioning. In: Bebis G, Boyle R, Parvin B, Koracin D, Remagnino P, Porikli F, et al., editors. Advances in Visual Computing. Lecture Notes in Computer Science. 5358: Springer Berlin Heidelberg; 2008. p. 268–77.
- 44.
Oliveira JVd, Pedrycz W. Advances in Fuzzy Clustering and its Applications: John Wiley& Sons, Inc.; 2007.
- 45. Muller K, Mika S, Ratsch G, Tsuda K, Scholkopf B. An introduction to kernel-based learning algorithms. Neural Networks, IEEE Transactions on. 2001;12(2):181–201.
- 46. Heijmans HJAM. Mathematical Morphology: A Modern Approach in Image Processing Based on Algebra and Geometry. SIAM Review. 1995;37(1):1–36.
- 47. Lynn PA. Direct on-line estimation of muscle fiber conduction velocity by surface electromyography. IEEE transactions on bio-medical engineering. 1979;26(10):564–71. pmid:540906
- 48.
Moritani T, Stegeman D, Merletti R. Basic Physiology and Biophysics of EMG Signal Generation. Electromyography: John Wiley & Sons, Inc.; 2005. p. 1–25.
- 49. Sokolova M, Lapalme G. A systematic analysis of performance measures for classification tasks. Information Processing & Management. 2009;45(4):427–37.
- 50.
Ebdon D. Statistics in geography. 2nd ed. Oxford Oxfordshire; New York, NY, USA: B. Blackwell; 1985. viii, 232 p. p.
- 51.
Hilbe JM. Negative binomial regression. Cambridge; New York: Cambridge University Press; 2007. xii, 251 p. p.
- 52.
StataCorp. Stata Statistical Software. College Station. TX: Stata Corporation; 2007.
- 53. Felzenszwalb PF, Zabih R. Dynamic Programming and Graph Algorithms in Computer Vision. IEEE transactions on pattern analysis and machine intelligence. 2011;33(4):721–40. pmid:20660950
- 54.
Boykov Y, Veksler O. Graph Cuts in Vision and Graphics: Theories and Applications. In: Paragios N, Chen Y, Faugeras O, editors. Handbook of Mathematical Models in Computer Vision. Boston, MA: Springer US; 2006. p. 79–96.
- 55. Vieira TM, Merletti R, Mesin L. Automatic segmentation of surface EMG images: Improving the estimation of neuromuscular activity. Journal of biomechanics. 2010;43(11):2149–58. pmid:20444452
- 56. Al-Kofahi Y, Lassoued W, Lee W, Roysam B. Improved automatic detection and segmentation of cell nuclei in histopathology images. IEEE transactions on bio-medical engineering. 2010;57(4):841–52. pmid:19884070
- 57. Couprie C, Grady L, Najman L, Talbot H. Power Watershed: A Unifying Graph-Based Optimization Framework. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2011;33(7):1384–99. pmid:21079274
- 58.
Merletti R, Afsharipour B, Piervirgili G. High Density Surface EMG Technology. In: Pons LJ, Torricelli D, Pajaro M, editors. Converging Clinical and Engineering Research on Neurorehabilitation. Berlin, Heidelberg: Springer Berlin Heidelberg; 2013. p. 1205–9.
- 59. Stegeman DF, Dumitru D, King JC, Roeleveld K. Near- and far-fields: source characteristics and the conducting medium in neurophysiology. J Clin Neurophysiol. 1997;14(5):429–42. pmid:9415390
- 60. Farina D, Muhammad W, Fortunato E, Meste O, Merletti R, Rix H. Estimation of single motor unit conduction velocity from surface electromyogram signals detected with linear electrode arrays. Medical and Biological Engineering and Computing. 39(2):225–36. pmid:11361250
- 61. Bilodeau M, Cincera M, Arsenault AB, Gravel D. Normality and stationarity of EMG signals of elbow flexor muscles during ramp and step isometric contractions. Journal of electromyography and kinesiology: official journal of the International Society of Electrophysiological Kinesiology. 1997;7(2):87–96.
- 62. Gazzoni M, Farina D, Merletti R. A new method for the extraction and classification of single motor unit action potentials from surface EMG signals. J Neurosci Methods. 2004;136(2):165–77. pmid:15183268
- 63. Soares FA, Carvalho JLA, Miosso CJ, de Andrade MM, da Rocha AF. Motor unit action potential conduction velocity estimated from surface electromyographic signals using image processing techniques. BioMedical Engineering OnLine. 2015;14(1):1–17.
- 64.
Smith SW. The scientist and engineer's guide to digital signal processing. 1st ed. San Diego, Calif.: California Technical Pub.; 1997. xiv, 626 p. p.
- 65.
Meijering EHW, Niessen WJ, Pluim JPW, Viergever MA. Quantitative Comparison of Sinc-Approximating Kernels for Medical Image Interpolation. In: Taylor C, Colchester A, editors. Medical Image Computing and Computer-Assisted Intervention—MICCAI’99: Second International Conference, Cambridge, UK, September 19–22, 1999 Proceedings. Berlin, Heidelberg: Springer Berlin Heidelberg; 1999. p. 210–7.
- 66.
Turkowski K. Filters for common resampling tasks. In: Andrew SG, editor. Graphics gems: Academic Press Professional, Inc.; 1990. p. 147–65.
- 67.
Mihalik J, Zavacky J, Kuba I. Spline interpolation of image. Radioengineering. 1995.
- 68. Lehmann TM, Gonner C, Spitzer K. Survey: interpolation methods in medical image processing. IEEE transactions on medical imaging. 1999;18(11):1049–75. pmid:10661324
- 69. Unser M, Aldroubi A, Eden M. Fast B-spline transforms for continuous image representation and interpolation. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1991;13(3):277–85.
- 70. Unser M. Splines: a perfect fit for signal and image processing. IEEE Signal Processing Magazine. 1999;16(6):22–38.
- 71. Allen DP. A frequency domain Hampel filter for blind rejection of sinusoidal interference from electromyograms. J Neurosci Methods. 2009;177(2):303–10. pmid:19010353
- 72. Ghaderi P, Marateb HR. Muscle Activity Map Reconstruction from High Density Surface EMG Signals with Missing Channels Using Image Inpainting and Surface Reconstruction Methods. IEEE transactions on bio-medical engineering. 2016.
- 73. Childers MK. Targeting the neuromuscular junction in skeletal muscles. American journal of physical medicine & rehabilitation / Association of Academic Physiatrists. 2004;83(10 Suppl):S38–44.
- 74. Nawab SH, Chang S-S, De Luca CJ. High-Yield Decomposition of Surface EMG Signals. Clinical neurophysiology: official journal of the International Federation of Clinical Neurophysiology. 2010;121(10):1602–15.