Identification of Myofascial Trigger Point Using the Combination of Texture Analysis in B-mode Ultrasound with Machine Learning Classifiers

1. Introduction

Nearly one-hundred million adults in the United States are affected by chronic pain (e.g., myofascial pain syndrome (MPS)), with an annual cost between $560 to $635 billion USD¹. MPS is one of the most prevalent musculoskeletal (MSK) pain disorders that occur in every age group and has been associated with primary pain conditions including osteoarthritis, disc syndrome, tendinitis, migraines, spinal dysfunction². MPS is characterized by myofascial trigger points (MTrPs)³. Active MTrPs (A-MTrPs) are spontaneously painful nodules while latent MTrPs (L-MTrPs) are palpable nodules that are painful when palpated.

MTrPs are classically defined as a ‘hyperirritable spot’ in skeletal muscle that is associated with a hypersensitive palpable nodule in a taut band⁴. The diagnostic criteria for MPS involves physical screening but studies have shown that manual detection of MTrPs is unreliable⁵. Techniques that are more quantitative in nature can help improve detection of MTrPs. Ultrasound is an attractive modality for this problem as it has been used to identify MTrPs^5–7, and are relatively cost-effective, non-invasive way to assess muscles, tendons and ligaments^8–10. It has been shown that MTrPs can be visualized and distinguished within a taut band on an affected muscle from normal tissue using US imaging techniques such as Doppler US and elastography US^8,11–14. Unfortunately, not all clinical US machines are equipped with elastography capability and both Doppler and elastography US techniques require comprehensive training sessions. Brightness mode (B-mode) US is readily available and is thus the preferred option for diagnosing and screening musculoskeletal disorders.

Since B-mode US imaging has high variability in echo intensity, texture features have been widely used to discriminate variables on B-mode US images. Texture features plays a vital role in computer vision and radiomics cases providing information such as muscle fiber orientation in the region of the MTrPs, and normal anatomy as well as the extent of adipose, fibrous, and other connective tissues within muscle¹⁵. Previous studies have shown that the muscle fibers within the MTrPs in the affected zone and the muscle fibers in the surrounding regions have different orientation in comparison with normal skeletal muscle.

Although, texture feature analysis of US images has been explored to distinguish MTrPs in affected muscle from normal tissue^9,16,17, there is currently no “gold standard” for the presence of MTrPs in US images. Previous research has used various methods of analyzing texture to tackle this problem such as using entropy characteristics¹⁶, gray level co-occurrence matrices (GLCM), blob analysis, local binary pattern (LBP), and statistical analysis¹⁷. A comprehensive overview of techniques employed for texture analysis or classification, categorize these techniques into four main categories:

• Transform-Based: Transform-based techniques employ a set of predefined filters or kernels to extract texture information from an image. Common filters include Gabor filters, and LBP^18,19. These filters highlight certain frequency components or local variations in pixel values, making them suitable for tasks where texture patterns are characterized by specific spatial frequencies or orientations.
• Structural: Structural techniques focus on describing the spatial arrangement and relationships between different elements in an image. They often involve identifying and characterizing specific patterns or structures within the texture (e.g., GLCM). These methods are valuable for capturing details related to texture regularity, directionality, or organization.
• Statistical: Statistical methods involve the analysis of various statistical properties of pixel intensities within an image or a region of interest (ROI). Common statistical features include entropy, contrast, correlation, homogeneity, energy, mean, and variance. These metrics quantify the distribution and variation of pixel values, providing insights into the texture’s overall properties, such as roughness, homogeneity, or randomness.
• Model-Based: Model-based methods involve fitting mathematical or statistical models to texture patterns in an image. These models can be simple, such as parametric distributions like Gaussian or Markov Random Fields, or more complex, such as deep learning models like convolutional neural networks. Model-based approaches are versatile and can capture intricate texture patterns, making them increasingly popular for texture analysis.

For many clinicians and investigators, the finding of one or more MTrPs is required to assure the diagnosis of MPS. A wide variety of studies describe the MTrPs as a hyperechoic band, hypoechoic elliptical region, or simply a different echo-architecture than the surrounding muscle tissue in clinical examination (e.g., US screening) ^14,20. In addition, MTrPs are typically known to be “knots” in the muscle. However, there remains a gap in knowledge regarding the optimal methods for characterizing these muscle structures (e.g., the associated texture pattern and edge). It was suggested that achieving an objective characterization and quantitative assessment of MTrPs properties has the potential to enhance their localization and diagnosis, while also facilitating the development of clinical outcome measures¹⁶.

Previously, Gabor filters were used to enhance the fiber orientation patterns and detect edges in the US images of the muscle^21,22. Also, LBP feature was used to characterize skeletal muscle (e.g., trapezius) composition in patients with MPS compared with normal healthy participants^9,23. Another proposed technique in which the statistical features were computed over the combination of Edge, GLCM, and LBP features, known as a “SEGL” method, to discriminate stone textures²⁴. It demonstrated the highest performance compared to LBP and GLCM features alone and may be a beneficial way for classifying MTrPs. As an extension of these image processing techniques, machine learning (ML) methods represent a powerful approach for automating the classification of muscle characteristics (e.g., muscle with MTrPs from healthy muscle).

ML techniques offer highly effective means of classification due to their ability to autonomously learn patterns and relationships from data^25–27. Futurists have anticipated that novel autonomous technologies, embedded with ML, will substantially influence healthcare. ML is focused on making predictions as accurate as possible, while traditional statistical models are aimed at inferring relationships between variables²⁷. ML offers advantages in terms of flexibility and scalability when contrasted with conventional statistical methods, allowing its utilization across various tasks like diagnosing, classifying, and predicting survival. Nevertheless, it is crucial to assess and compare the accuracy of muscle characterization through traditional statistical methods and ML within the context of clinical screening²⁸. Supervised ML algorithms (e.g., neural networks, decision trees, etc.) can generalize from training data to make accurate predictions or classifications on new, unseen data. Their adaptability allows them to handle diverse domains and tasks, making them invaluable tools for tasks ranging from image recognition to medical diagnosis, enhancing efficiency and precision in decision-making processes²⁶.

Building upon the advantages of ML and recognizing its potential in healthcare, this study delves into the utilization of various texture feature approaches and ML techniques to classify and characterize MTrPs in US images. We investigate different texture feature approaches (i.e., LBP, Gabor, SEGL method and their combination with texture features) extracted from US images to classify MTrPs. We further employ various ML techniques as well as traditional statistical analysis to explore the effectiveness of the extracted features from the US images to characterize and classify the muscle between A-MTrPs, L-MTrPs and healthy control.

2. Materials and Methods

2.1. Participants

Participants (n = 90) were recruited from the musculoskeletal/pain specialty outpatient clinic at the Toronto Rehabilitation Institute. The upper trapezius muscle of all participants was examined. All participants underwent a physical examination by a trained clinician on our team (BD), who determined the presence or absence of MTrPs (i.e., as A-MTrPs, L-MTrPs) in the upper trapezius muscle according to the standard clinical criteria defined by Travell and Simons⁴, and through visual confirmation on B-mode US. Participants who demonstrated no symptoms or history related to neuromuscular disease, based on diagnostic criteria, were included in this study. Each participant’s muscle was labeled as A-MTrPs, L-MTrPs or healthy control (Table 1).

Informed consent was obtained from all participants involved in the study prior to participating and their upper trapezius muscles were included in our study. This individual protocol was approved by the Institutional Review Board of the University Heath Network (UHN) (protocol code 15-9488). The procedures adhered to the guidelines set out in the Declaration of the World Medical Association of Helsinki.

2.2. Ultrasound Acquisition Protocol & Pre-processing

US videos were acquired using an US system (SonixTouch Q+, Ultrasonix Medical Corporation, British Columbia, Canada) with a linear ultrasonic transducer of 6-15MHz and a depth set to 2.5 cm. The acquisition settings including time gain compensation, depth and sector size were held constant across all participants. Acquisition was performed by the experienced sonographer with the participant sitting upright in a chair with their arms relaxed on their sides, and forearms resting on their thighs. The transducer was placed on the skin in the center of the trapezius muscle (i.e., the midpoint of the muscle belly between the C7 spinous process and the acromioclavicular joint) with enough gel to cover the entire surface (Figure 1). A ten second video (sampling frequency: 30 frames / second) of the trapezius muscle from each side per participant was recorded by moving the transducer towards the acromioclavicular joint (parallel to the orientation of the muscle fibers) at approximately 1 cm/s, generating 300 images per participant for analysis (Figure 1). While recording the video, the researcher manipulated the transducer’s position to reduce artifacts and mitigate muscle distortion caused by the transducer, such as applying downward pressure. From each video, 4 unique frames/images were manually selected out of 300 B-mode images. These selected images captured various section of the muscle (i.e., lateral to medial), and were used to validate the presence or absence of MTrPs evident in the video (Figure 2A). Images from each side of a participant (e.g., left, or right trapezius) were treated as independent sites.

A region of interest (ROI) of the muscle (i.e., the region between the upper trapezius muscle’s superior and inferior fascia) were manually extracted from the acquired images by visual localization. Theses muscle ROIs were further analyzed using the following texture features: LBP, Gabor, SEGL method and combination of them with statistical features.

2.3. Texture Feature Analyses

I. Local Binary Patterns. One of the most popular texture feature analysis operators is LBP, a rotationally invariant texture feature²⁹. It can evaluate the local spatial patterns and contrast of grayscale images. This technique provides eigenvalues for the different patterns in an image, such as edges, and corners within the ROI. LBP was calculated for every B-mode image using the following equation below (Eq. 1).

$$\mathit{L}\mathit{B}{\mathit{P}}_{\mathit{P},\mathit{R}}=\sum _{\mathit{p}=0}^{\mathit{P}-1}\mathit{s}\left({\mathit{g}}_{\mathit{p}}-{\mathit{g}}_{\mathit{c}}\right)2^{\mathit{p}},\mathit{s}\left(\mathit{x}\right)=\left\{\begin{array}{c}1,x\ge 0\\ 0,x<0,\end{array}\right.$$

(1)

Where P is the number of pixels within the neighborhood, $g_p$ represents the p^th neighboring pixel, and c represents the center pixel. In our study, a 3 by 3 neighborhood was used, and its central pixel intensity was compared with its surrounding eight neighbor pixels³⁰. The eight neighbors are labeled using a binary code {0, 1} obtained by comparing their values to the central pixel value. If the tested gray value is below the gray value of the central pixel, then it is labeled 0, otherwise it is assigned the value 1.

$s\left(x\right)$ is the obtained binary code, $s\left(x\right)$ is the original pixel value at position x and $s\left(0\right)$ is the central pixel value. With this technique there is 256 (2⁸) possible patterns (or texture units). The obtained value is then multiplied by weights given to the corresponding pixels. Summing the obtained values gives the measure of the LBP (Eq. 1). The central pixel is replaced by the obtained value (LBP measure). A new LBP image is constructed by processing each pixel and its 3x3 neighbors in the original image.

The LBP is calculated across the entire image and the outer layer of the image (i.e., the border), which does not have eight neighbors, was replaced with the next closest layer’s values (Figure 2C,D).

II. Gabor Feature. Gabor filtering or Gabor filter bank for image textural analysis was introduced by Daugman and used in pattern analysis applications^31–33. Gabor features are obtained by subjecting each filtered image to a nonlinear transformation and computing a measure of “energy” in a window around each pixel in each response image. The Gabor filter is characterized by preferred orientation spatial frequency. The Gabor filter-based features are directly extracted from the gray-level images (i.e., B-mode images). In the spatial domain, a two-dimensional Gabor filter is a Gaussian kernel function modulated by a complex sinusoidal plane wave, defined as (Eq. 2)³⁴:

$$\mathit{G}\left(\mathit{x},\mathit{y}\right)=\mathbf{e}\mathbf{x}\mathbf{p}\left\{\frac{-1}{2}\left(\frac{{\mathit{x}}^{\mathit{\prime }2}}{{\mathit{\sigma }}_1^2}+\frac{{{\mathit{y}}^{\mathit{\prime }}}^2}{{\mathit{\sigma }}_2^2}\right)\right\}\mathit{c}\mathit{o}\mathit{s}\left(2\mathit{\pi }\mathit{ƒ}{\mathit{x}}^{\mathit{\prime }}+\mathit{\phi }\right),$$

(2)

$${\mathit{x}}^{\mathit{\prime }}=\mathit{x}\mathit{s}\mathit{i}\mathit{n}\mathit{\theta }+\mathit{y}\mathit{c}\mathit{o}\mathit{s}\mathit{\theta },{\mathit{y}}^{\mathit{\prime }}=-\mathit{x}\mathit{c}\mathit{o}\mathit{s}\mathit{\theta }+\mathit{y}\mathit{s}\mathit{i}\mathit{n}\mathit{\theta },$$

where ƒ is the spatial frequency of the wave at an angle θ with the x axis, ${\sigma }_1$ and ${\sigma }_2$ are the standard deviations of the 2-D Gaussian envelope, and $\phi$ is the phase.

Gabor features were calculated using the Gabor Feature Extraction function in MATLAB (The MathWorks, Natick, USA) created by Haghighat et al.³⁵; the outputs of this function were used for the Gabor feature. In summary, in their study, forty Gabor filters were calculated in five scales (frequencies) and eight orientations in steps of 45 degrees (i.e., θ: 0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°) for each image (5 x 8). The bank of Gabor filters (40 Gabor filters) extracted the variations in different frequencies and orientations from each image and generated forty ‘Gabor feature image’ corresponding to each B-mode US image (Figure 2B). All the Gabor feature images were normalized to 0 to 255.

III. SEGL Method. This method was proposed by Fekri Ershad S. and used a combinational feature extraction method for textual analysis²⁴. The term SEGL stand for (Statistical, Edge, GLCM, and LBP). This method first uses the LBP image as an input image, described in section I. The GLCM is then calculated on the LBP images, and the edge feature is then calculated from the resulting image.

GLCM was proposed by Haralick and Shanmugam³⁶. GLCM depicts how often different combinations of gray levels co-occur in an image. The GLCM is created by calculating how often a pixel with the intensity value i occurs in a specific spatial relationship to a pixel with the value j. In this study, GLCM were calculated along 8 directions (θ: 0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°) with an empirically determined distance (offset = one pixel).

Edge detection is the process of localizing pixel intensity transitions that has been used to extract information in the image by object recognition, target tracking, segmentation, and etc. It is defined by a discontinuity in gray-level values or a boundary between two regions with relatively distinct gray level values³⁷. Canny Edge detection method was used, as previous literature has shown that Sobel Edge detection method cannot produce smooth and thin edge compared to Canny method³⁸, and this may reduce structural information (e.g., edge) of MTrPs within the muscle.

Finally, the seven statistical features, described below (Section, IV), were calculated over the edge detected images. As described above, the GLCM was calculated in 8 directions; as a result, the total number of extracted attributes is 56 statistical values (8 x 7).

IV. Statistical Feature. Statistical features were used to measure the image variation. In our study, 7 statistical features: entropy, energy, mean, contrast, homogeneity, correlation, and variance were computed over each calculated image from the texture features described above (sections I, II, III) in addition to the B-mode US images using the following equations in Table 2 (Eq. 3 - 9). Therefore, for each US image: 7 (B-mode), 7 (LBP), 280 (Gabor, i.e., 40 Gabor features x 7), and 56 (SEGL: i.e., 8 directions x 7) statistical values were calculated for each approach. Figure 3 shows the summary of the methods that were used for feature extraction of each image in this study.

2.4. Classification Approaches, Evaluation

The statistical features (values) from each image within each method were classified using a variety of supervised ML approaches to discriminate muscle with MTrPs (A-MTrPs and L-MTrPs) from healthy muscle. Classification accuracy was evaluated using following approaches outlined in Table 3. These algorithms were implemented in python using the Scikit Learn library using their common hyperparameters and their associated default values as outlined below in the Table 3.

The NN was a 3-layer network (256 x 512 x 128) with dropout layers between each (70% of nodes dropped). All activation functions were Rectified Linear Units (ReLU). The output layer was a 3-node output with an activation function of SoftMax. The NN was trained for 250 epochs with early stopping criterion of 7 epochs of no improvement in the validation loss. The learning rate was the default set by Keras and was decreased by a factor of 0.1 after 3 epochs of no improvement to a minimum of learning rate = 0.00001.

The inputs to the ML approaches were the 7 statistical features calculated on the B-mode image, the LBP images, the Gabor filtered images or the SEGL output images.

Due to the low number of images a leave-one-out site was used to better evaluate the performance. After the leave-one-out procedure the remaining examples were either used as is for training (i.e., LR, DT, RF, k-NN, NB, SVM) or split into 75% training and 25% validation sets (i.e., NN).

For all approaches, the classification accuracy, F1-score, sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV) was used to evaluate the performance of the methods.

2.5. Ensemble Approaches & feature importance and Statistical Analysis

All the ML approaches with their associated parameters (shown above in Table 3) were used to further investigate the combination of statistical features (e.g., entropy, mean etc.) extracted from the 4 texture features (B-mode LBP, Gabor, SEGL) within each ML algorithm. This was performed to see the effects that statistical feature had on the different classification algorithms. This was done using:

(1)

A single statistical feature from all texture feature approaches on all ML approaches

(2)

A leave-one-feature out approach from all texture feature approaches on all ML approaches

Of the ML approaches investigated, the parameters that provided the best results for each technique were used to investigate ensemble approaches as well as feature importance. These best results associated with selected parameters were shown with asterisks (*) for each ML approach in the Table 3. These techniques were used in a simple majority vote to determine if an ensemble approach could provide improved results on the 4 different images (B-mode, LBP, Gabor, SEGL).

In addition, the mean and standard-deviation (SD) of each statistical feature within each group (A-MTrPs, L-MTrPs, healthy control) and each texture feature (B-mode LBP, Gabor, SEGL) were computed for further traditional statistical analysis (i.e., one-way analysis of variance (ANOVA)). In each texture feature, the ANOVA was computed between each statistical feature (mean values) to compare all 3 groups: A-MTrPs, L-MTrPs, healthy control. For example, ANOVA was performed between the calculated entropy of Gabor filtered image for A-MTrPs, L-MTrPs, and healthy control.

3. Results

Table 4 shows the highest classification accuracy (%), F1-score, sensitivity, specificity, PPV, and NPV for the ML algorithms: between the 7 statistical features for each method: B-mode, LBP, Gabor feature, and SEGL for the majority voting; between each statistical feature of four methods, separately; and between each six-statistical features of four methods, separately, to classify the muscle with MTrPs from healthy muscle.

As shown in Table 4, the ML algorithms for ‘without variance’ with accuracy = 51.67% and F1-score = 0.518, and ‘correlation’ with accuracy = 53.33% and F1-score = 0.4861 showed highest performance. In addition, the ML algorithms for the majority voting approximately performed better than NN between the statistical features for all four methods (B-mode: 49.44% vs. 42.5%, LBP: 47.22% vs. 47.22 %, Gabor: 48.86% vs. 47.78%, and SEGL: 49.44% vs. 43.89).

In addition, the classification accuracies vary between 42.5% (min) and 53.33% (max); there was a low difference (10.83%) between the performance of all approaches between the statistical features and methods (Table 4).

Table 5 shows the traditional statistical analysis (i.e., one-way ANOVA) results between each statistical feature within each method (i.e., B-mode, LBP, Gabor filter, and SEGL) when comparing all 3 groups: A-MTrPs, L-MTrPs, and healthy controls. Statistically significant differences were seen for all statistical features within each method except: the ‘entropy’, ‘contrast’, and ‘energy’ of B-mode, ‘mean’ and ‘correlation’ of Gabor feature (p < 0.05) (Table 5). Nevertheless, Table 5 demonstrated the similar mean and low standard deviation (SD) statistical values between the groups and the four approaches (i.e., Gabor feature, SEGL method, B-mode, LBP).

Figure 3 shows the confusion matrixes of the ML algorithms with the best performance between the 7 statistical features for four methods: B-mode, LBP, Gabor filter, and SEGL for the majority voting; between four methods for each statistical feature, separately; and between four methods for six statistical features, separately. According to the true and predicted classes/labels (e.g., true negative, true positive, false negative, and false positive) for each group, ‘without variance’ approach presents better performance than other approaches for classifying all 3 groups: A-MTrPs, L-MTrPs, and healthy controls.

In addition, as shown in Figure 3, ‘Gabor filter’ approach show better performance than the other three methods between the 7 statistical features for classifying all 3 groups: A-MTrPs, L-MTrPs, and healthy controls.

Figure 3. Confusion matrixes of the learning algorithms with the best performance between the 7 statistical features for four methods: B-mode itself, LBP, Gabor filter, and SEGL for the majority voting and NN, between four methods for each statistical feature, separately, and between four methods for six statistical features, separately, for discriminating the A-MTrPs, L-MTrPs, and healthy controls. This figure shows how many predictions are correct and incorrect per class.

Figure 3. Confusion matrixes of the learning algorithms with the best performance between the 7 statistical features for four methods: B-mode itself, LBP, Gabor filter, and SEGL for the majority voting and NN, between four methods for each statistical feature, separately, and between four methods for six statistical features, separately, for discriminating the A-MTrPs, L-MTrPs, and healthy controls. This figure shows how many predictions are correct and incorrect per class.

4. Discussion

Our study investigated the effectiveness of combining statistical features (specifically focused on edges and spots) derived from US images, utilizing B-mode, LBP, Gabor feature, and SEGL method for the purpose of accurately discriminating muscle with A-MTrPs, L-MTrPs, from healthy muscle (control). Our findings indicate that this combination approach (i.e., statistical-features, + B-mode, LBP, Gabor, and SEGL) did not achieve a high level of accuracy in distinguishing between A-MTrPs, L-MTrPs and healthy. Although, this combination approach showed approximately higher accuracies for B-mode and SEGL method, compared to the LBP and Gabor feature (49.44% and 49.44% vs 47.22% and 48.89%, respectively). We hypothesized that structural, statistical, and combination of them could better classify muscle with MTrPs from healthy muscle. However, the overall accuracies obtained from these combination approaches exhibited a similar range, ranging from 47.22% to 49.40%.

Through the utilization of traditional statistical analysis (i.e., ANOVA) and automated techniques (i.e., seven ML classifier approaches), this study reveals the significant statistical differences among the three groups within the most of statistical features extracted from the four methods (B-mode, LBP, Gabor feature, and SEGL method) (p < 0.05) (Table 5). However, these observed differences do not exhibit sufficient ability to enable high classification accuracy of these three groups, primarily due to the low degree of variation (i.e., SD) and overlapping observed in statistical values between the groups and the four approaches (i.e., Gabor feature, SEGL method, B-mode, LBP), as seen in Table 5.

Notably, the results demonstrated that the ‘correlation’ and ‘mean’ statistical features compared to other statistical features exhibited better discriminatory ability in distinguishing between muscle with MTrPs and healthy muscle (Table 4). These ML algorithms yielded the highest accuracies of 53.33% and 52.5% when utilizing the approaches involving ‘correlation’ and ‘mean’, respectively. Surprisingly, these accuracies decreased to the lowest values of 49.17% and 50.83% when the approaches ‘without correlation’ and ‘without mean’ were employed. It is important to note that despite an increase in the number of training features, which sometimes enhances classification performance, the aforementioned approaches experienced a decline in accuracy. Conversely, other statistical features (i.e., entropy, energy, contrast, variance, homogeneity) exhibited an improvement in performance with an increase in the number of training features. This discrepancy suggests the presence of distinct characteristics among the US images characterizing these three groups (i.e., A-MTrPs, L-MTrPs, and healthy controls) that may be better captured and evaluated by the ‘correlation’ and ‘mean’ statistical features. This suggests differences (different gray-scale values (i.e., echo-intensities)) between the structural patterns of muscle with MTrPs (A- MTrPs and L- MTrPs) and healthy muscle. This is in the agreement with the previous findings that showed the MTrPs as hypoechoic (dark grey) nodules in B-mode US images^39,40.

However, traditional statistical values showed that no statistically significant differences were observed for the ‘correlation’ and ‘mean’ statistical features of the Gabor feature approach (p = 0.0857 and p = 0.2338). This can be attributed to the fact that the Gabor feature measures the gray level of US images⁴¹, and there are similar mean and low standard deviations between and within the correlation and mean values of A-MTrPs and L-MTrPs, as indicated in Table 5. These findings align with previous studies that have reported muscle with MTrPs to exhibit anisotropy⁴².

Overall, MTrPs have been identified and labeled as hypoechoic (dark grey) nodules in previous US studies^39,40. However, recent research has proposed the identification of MTrPs as large hypoechoic contracture knots, which also exhibit smaller hyperechoic “speckles” within the hypoechoic contracture knot⁴³. The presence of these “speckles” can affect the structural information of MTrPs within the muscle ROI and interfere with the characterization of muscle with MTrPs using texture feature analysis. For instance, entropy is capable of describing homogeneity and randomness in the observed patterns in US images, while LBP depict the structural elements (spots, edges, etc.) of US backscatter. Consequently, the presence of different patterns within muscles affected by MTrPs may lead to variations in the values of calculated texture features within each group, thereby reducing the predictive accuracy of the ML techniques. The results seen in this study are in the agreement with the performance found in a recent study using convolutional neural network for discriminating MTrPs from healthy muscle (F1-score = 0.4697 Vs. F1-score = ranging from 0.4582 to 0.4855 in our study)⁴⁴.

Another critical consideration is the relationship between the US image and the clinical scenario. Existing literature has proposed certain clinical criteria for MTrPs; however, these criteria have not been clearly associated with specific US abnormalities. Currently, most researchers in this field concur that MTrPs are a physical entity that exhibits a spherical or elliptical shape, although this assumption is temporary and subject to further investigation. Therefore, it becomes crucial to identify the characteristics of the “center” of this entity, which may appear as hypoechoic and exhibit either homogeneous or heterogeneous features on the US image. Additionally, defining the “border zone” that separates this region from the surrounding normal area is necessary, as preliminary methods suggest that this border or transition zone may provide more valuable information than the lesion (i.e., hypoechoic contracture knot) itself ⁴⁵. Moreover, in cases where a patient experiences pain but does not present a formed MTrPs, it is uncertain if there is an ‘at-risk’ area that later transforms into a visually defined spherical/elliptical MTrPs, and the US characteristics of this at-risk area require further clarification. In our study, the entire muscle region, including the border zone, MTrPs, and normal areas, was considered, which may impact the calculated statistical values. Therefore, this study suggests additional research utilizing quantitative techniques such as texture analysis to analyze the effects of border zones and central regions of visually evident MTrPs within US images.

Here, to the best of our knowledge, this study represents the first attempt to investigate the potential of combining texture features mainly focused on structural information (e.g., spots, edges) within the muscle for discriminating muscle with MTrPs from healthy muscle. Moreover, previous studies have primarily relied on traditional statistical analysis methods, such as ANOVA, as opposed to automated analysis techniques, such as ML classifiers, for distinguishing between these three groups^17,46. Although our study demonstrated significant results using ANOVA in differentiating the groups, the classifier approaches did not show as robust outcomes. Therefore, we propose further exploration and utilization of automated classifier approaches (e.g., ML techniques) with the incorporation of texture features derived from US images, to achieve improved discrimination between MTrPs and healthy controls. Because the key characteristic that sets ML methods apart is their reliance on data-driven processes rather than user-selected approaches, with the aim of achieving precise predictions. This approach helps avoid the mistake of using an inappropriate traditional statistical model on the dataset, which could limit accuracy²⁷. We believe that these imaging techniques have the potential to aid clinicians in confirming the presence of MTrPs and enhancing diagnostic capabilities in the future.

4.1. Limitations

One limitation of this study lies in the proper definition and localization of the region of MTrPs within the US images from the muscle for the analysis. The MTrPs appeared to be 0.05 to 0.5 cm² nodules of varying hypo-echogenicity. However, the consensus was not absolute as there are notable deviations in described presentation, specifically: clusters of hypoechoic structures, and curved muscle form with MTrPs, as well as the hyperechoic presentation. Therefore, in this study, the entire region of muscle was cropped and considered for analysis, limiting the assessment to the textural information of the entire muscle region (ROI) rather than focusing solely on MTrPs themselves. By not confining the ROI to just the MTrPs, which requires more expertise, solving it with just the extraction of muscle would lower the barrier to use the developed technique. However, this approach may reduce the ability to detect required differences between the statistical features (values) within each group.

Furthermore, while hypoechoic images are generally associated with hyperperfused areas and hyperechoic images with hypo perfused areas⁴⁷, it is important to acknowledge the possibility of image artifacts, such as anisotropy, in our patient population. Hypoechoic areas may also arise from acoustic shadowing behind calcifications, lymph nodes, and certain pathological conditions. However, in this study, the manual selection of images aimed to alleviate the presence of any artifacts, and to maintain consistency across all participants, the entire muscle region (ROI) was evaluated.

5. Conclusions

In conclusion, this paper sheds light on the utilization of texture features and their combination (i.e., statistical features with B-mode, Gabor, LBP, and SEGL method) for the classification of A-MTrPs, L-MTrPs, and healthy muscle. The focus was to capture structural information such as edges, spots, and other relevant features. In comparison to traditional statistical analysis methods (e.g., ANOVA), the employed classification approaches did not achieve high classification results, likely due to the significant overlapping observed among the statistical values between the groups (maximum reported accuracy of 53.33%) but were still much higher than chance suggesting the feasibility of discriminating between the three groups; however, some improvements may be required.

Therefore, this study highlights the need for further investigations to explore the potential of extracting advanced texture features in combination with non-traditional statistical analysis, for effectively identifying MTrPs in contrast to healthy muscle. Such endeavors can contribute to the development of more robust diagnostic criteria based on US image characteristics. The findings from these future studies hold promise for the development of improved mechanisms to aid in the accurate identification and diagnosis of MTrPs.