We conducted a retrospective collection of T2-weighted mid-sagittal cervical spine MRI scans from patients of varying ages across ten hospitals within Mainland China. Subsequently, we segregated tese scans into testing datasets. The participating hospitals comprised Longhua Hospital Affiliated to Shanghai University of Traditional Chinese Medicine, Shanghai Municipal Hospital of Traditional Chinese Medicine, Dongzhimen Hospital Affiliated to Beijing University of Chinese Medicine, The Second Affiliated Hospital of Guangzhou University of Chinese Medicine, Shenzhen Pingle Orthopedic Hospital Affiliated to Guangzhou University of Chinese Medicine, Ningxia Hui Autonomous Region TCM Hospital and TCM Research Institute, Luoyang Orthopedic-Traumatological Hospital Of Henan Province, Affiliated Hospital of Shaanxi University of Chinese Medicine, Shenzhen Traditional Chinese Medicine Hospital, and Suzhou TCM Hospital Affiliated to Nanjing University of Chinese Medicine. Scanning sequences, parameters, and scanners (including strength and vendor) were not standardized across sites to evaluate the generalization of PDB-SLIC in a heterogenous, real-world clinical sample.

This study adhered strictly to the principles outlined in the Declaration of Helsinki and received approval from the Institutional Review Board (IRB) of Longhua Hospital Shanghai University of Traditional Chinese Medicine (approval number 2023LCST006). The IRB thoroughly reviewed and approved the research protocol, ensuring the protection of patient privacy and the confidentiality of the clinical data utilized. Given the retrospective nature of the study and the de-identified nature of the data employed, patient consent was waived. All methodologies were conducted in accordance with relevant guidelines and regulations. Prior to analysis, data were anonymized and de-identified uphold patient confidentiality. Furthermore, all researchers involved in data analysis underwent training and adhered to stringent data protection practices.

To address the issue of inconsistent brightness levels in different images, (Figure 1A) we employed CLAHE equalization to enhance the overall contrast of the images (Figure 1B). And bilateral filtering was used for preprocessing, which effectively addresses issues such as noise interference, unclear image edges, and image blurring, thus further improving the stability of our algorithm (Figure 1C). The calculation formula for bilateral filtering is shown below. G σ s ∥ p − q ∥ = ⁡ exp   − ∥ p − q ∥ 2 2 σ s 2 (1) G σ r ∣ I p − I q ∣ = ⁡ exp   − ∣ I p − I q ∣ 2 2 σ r 2 (2) I ¯ p = 1 W p ∑ q ∈ S G σ s ∥ p − q ∥ G σ r ( ∣ I p − I q ∣ I q (3)

Where W p = ∑ q ∈ S G σ s ∥ p − q ∥ G σ r ∣ I p − I q ∣ is the normalization factor, G σ s represents the spatial domain kernel, G σ r represents the pixel domain kernel, I ¯ p represents the image after bilateral filtering, I q represents the input image, p 、 q represents the pixel position, I p 、 I q represents the corresponding pixel value.

Superpixels refer to irregular pixel blocks composed of adjacent pixels with similar texture, color, brightness, and other features that have certain visual significance. By grouping pixels based on the similarity of features between them, a small number of superpixels can be used to represent the characteristics of a large number of pixels, achieving low-complexity and efficient segmentation.

The superpixel segmentation method designed in this study generates superpixels through multiple K-means clustering on the image. In order to make the generated superpixels more regular, a distributed hexagonal distribution-based initial clustering center method is designed specifically for vertebral features in this study (Figure 2).

After bilateral filtering, in order to improve the clustering speed and obtain more global image information, a downsampling process is adopted to remove even rows and even columns of the upper-level pyramid image with K-means clustering.

The objective function D is calculated for each clustering center and all pixels within its range. A pixel is assigned to the class whose objective function D with that clustering center is the minimum. The objective function D is defined as follows: d c = L j − L i 2 (4) d s = x j − x i 2 + y j − y i 2 (5) D = d c + m S c 2 × d s (6)

Where m is the weighting factor between brightness and spatial differences. In this experiment, m is set to 15. S c represents the column sampling distance of the clustering centers.

After the initial K-means clustering, the positions and brightness information of the clustering centers are updated using the average values of the position and brightness features of the superpixels, respectively. A new clustering center in the original image corresponds to one updated clustering center in the downsampled image.

Next, iterative clustering with the objective function D is performed within their respective ranges based on the new clustering centers. Figure 3 shows the different superpixel segmentation performance under different numbers of superpixels.

The integration of Pyramid-based K-means clustering methodology with multi-resolution images serves to broaden the scope of exploration for clustering centroids, thereby mitigating a constraint inherent in the SLIC algorithm pertaining to its capacity to encapsulate the overarching characteristics of images. Such augmentation markedly amplifies the algorithm’s efficacy in delineating the boundaries of cervical vertebral bodies with heightened precision.

In the stage of superpixel merging, we use the brightness distance metric function Dc to measure the similarity between the seed superpixel and its neighboring superpixels. The definition of the brightness distance metric function Dc is as follows: D c = L n − L m 2 (7) where m represents the seed superpixel, and n represents the neighboring superpixel.

Next, a label set, C, and a candidate set, S, are defined. All superpixels are initially stored in the candidate set. The top-left superpixel of the image is set as the first seed superpixel and assigned a label. The newly formed superpixel resulting from the fusion of the superpixel in the label set with the seed superpixel continues to search for neighboring superpixels within the candidate set until no neighboring superpixels are found in the candidate set. The seed superpixel must be an unlabeled superpixel until all superpixels are labeled, resulting in a complete region segmentation. At this point, doctors can manually select the vertebral regions they consider appropriate to generate the corresponding segmentation results.

After obtaining the segmented region of the vertebrae, the vertebral centroid is calculated by the ratio of the horizontal and vertical centroids of the vertebrae to the sum of the pixels in that region (Figure 4). The centroid calculation formula is shown below. x a i = ∑ x = ⁡ min ⁡ x i max ⁡   x i ∑ y = ⁡ min ⁡ y i max ⁡ y i x   * P xy ∑ x = ⁡ min ⁡ x i max ⁡   x i ∑ y = ⁡ min ⁡ y i max ⁡ y i P xy (8) y a i = ∑ x = ⁡ min ⁡ x i max ⁡   x i ∑ y = ⁡ min ⁡ y i max ⁡ y i y   * P xy ∑ x = ⁡ min ⁡ x i max ⁡   x i ∑ y = ⁡ min ⁡ y i max ⁡ y i P xy (9)

Where i is the ith vertebra. For CCL-A, i ranges from 1 to 4, while for CCL-B, i ranges from 1 to 3.   x i denotes the set of horizontal coordinates of the ith vertebra, y i represents the set of vertical coordinates of the ith vertebra, and P xy represents the pixel value when the height coordinate is x and the width coordinate is y, with P xy value of 0 or 1.

CCL-A: After obtaining the centroid for each vertebra, the acute angle ( θ a ) formed by the line connecting the centroids of C2 and C3 and the line connecting the centroids of C6 and C7 is calculated.

CCL-B: After obtaining the centroid for each vertebra, the acute angle ( θ b ) formed by the line connecting the centroids of C2 and C4 and the line connecting the centroids of C4 and C7 is calculated.

We also propose a calculation method based on the minimum bounding rectangle fitting to achieve efficient and accurate automatic extraction of the tangent line on the upper edge of the vertebra and measure the corresponding slope. By sequentially locating the pixels in the segmented region of the vertebra to determine four extreme points, the minimum perimeter fitting rectangle is obtained using the rotating calipers method. The upper edge of the rectangle is approximated as the tangent line on the upper edge of the vertebra. The calculation formula for C7 slope is shown below (Figure 5). θ slope i = ⁡ arccos   w a P x min i P x max i → × 180 π (10)

Where i represents the ith vertebra, with i being 6. P x min i and P x max i represent the pixels corresponding to the minimum and maximum horizontal coordinates in the segmented region of the vertebra. w a represents the horizontal distance between these two points.

One of the most important requirements for superpixels is to preserve adherence to object boundaries (Yuan et al., 2021). In order to quantitatively compare the performance of SLIC(Achanta et al., 2012), VBseg (Barbieri et al., 2015), and PDB-SLIC in vertebral segmentation quality across different numbers of superpixels, precision, recall, and the Jaccard index are utilized to assess the similarity between algorithmic segmentation and true segmentation of the vertebral region. Additionally, the segmentation performance is compared with classical deep learning networks such as FCN (Shelhamer et al., 2017), DeeplabV3 (Chen et al., 2017), and U-Net (Ronneberger et al., 2015) to identify the best segmentation method.

The comparison of performance parameters among disparate algorithms across numerous superpixel counts, utilizing one-way analysis of variance (ANOVA). The test results are initially subjected to the Shapiro-Wilk test for normality and the Levene’s test for homogeneity of variances, the outcomes of which dictate the appropriate method for variance analysis.

To compare the performance of the semi-automatic measurement method with experienced spinal surgeons, a subset was randomly extracted from the cervical spine MRI dataset, consisting of 464 T2-weighted cervical MR images. Two spinal surgeons with 5–15 years of clinical experience were recruited to participate in this performance comparison study. Cervical curvature and C7 slope were manually measured using ImageJ2 software (National Institutes of Health, Bethesda, Maryland, United States, https://imagej.net/software/imagej2/) (Rueden et al., 2017). Each image was measured twice independently by both physicians, with the average of the measurements used as the final result. The manual results were then compared with the machine measurements for consistency assessment. The consistency between manual measurements and semi-automatic measurements was assessed using intraclass correlation coefficient (ICC) analysis.

To validate the performance of PDB-SLIC across different hospital settings, MRI devices, and age demographics, we categorized the cervical MRI dataset based on age group. PDB-SLIC was utilized to measure cervical curvature, enabling the analysis of parameter distribution trends across various age cohorts. Our focus was primariy on cervical curvature, C7 slope, and their ratio, with correlation analysis conducted across the datasets.

The distribution of cervical curvature was summarized using means ± standard deviations (SD). To assess the effects of different hospitals, MRI machines, and gneders on cervical curvature measurements, we examined data for normality and equal variances. Multivariate ANOVA or Kruskal-Wallis tests were applied accordingly. Additionally, for comparisons among age groups, we employed one-way ANOVA with Welch’s adjustment and Games-Howell post hoc analysis. Pearson’s correlation analysis was utilized to evaluate the correlations among cervical curvature, C7 slope, and their ratio.

IBM SPSS Statistics for Windows, version 26.0 (IBM Corp., Armonk, NY, United States) was applied to complete the statistical analysis. A p-value less than 0.05 was considered statistically significant.

A total of 4,258 cervical spine MR images from ten hospitals were included in the study, with patient ages ranging from 20 to 90 years. The MR manufacturers included GE, Philips, and Siemens, and the MRI field strength was either 1.5T or 3.0T. The baseline characteristics of the cervical spine MRI image dataset are shown in Table 1.

Figure 6 illustrated the performance evaluation of the SLIC algorithm, showcasing its inferiority compared to the other algorithms. Notably, significant deficiencies in vertebra segmentation are observed, particularly evident whem the number of superpixels is set to 100, resulting in substantial loss of are in vertebra segmentation. While the Vbseg algorithm achieves superior vertebral body boundary delineation, challenges arise in segmenting local extreme value regions within the vertebrae. This difficulty stems from the algorithm’s utilization of the Otsu method for secondary segmentation of superpixels, leading to enhanced superpixel quality but complicating the segmentation of local extreme value region and subsequent superpixel merging. Consequently, the Vbseg algorithm fails to extract the local extreme value regions within the vertebral body. In contrast, PDB-SLIC demonstrates consistent excellence in both the vertebral body region and boundary segmentation, surpassing the performance of the other algorithms overall.

In Figure 7, precision, recall and Jaccard index evaluation of vertebral body segmentation in cervical MR images are depicted across varying superpixel numbers. A one-way ANOVA, utilizing Fisher’s method, was conducted to analyze the results. The findings reveal that both PDB-SLIC and Vbseg algorithms show superior performance compared to the SLIC algorithms in terms of precision (F = 22.517; p < 0.001) and Jaccard index (F = 14.842; p < 0.001), while the disparity in UE values is statistically insignificant (F = 0.990; p = 0.391). Notably, the application of the Otsu method in the Vbseg algorithm contributes to the inability to separate the local extreme value regions within the vertebral body, consequently leading to decreased accuracy and Jaccard index. Upon analysis of the vertebral body segmentation results, both PDB-SLIC and SLIC algorithms demonstrate optimal performance when the number of superpixels is set of 200, whereas the Vbseg algorithm performs best with 500 superpixels. Therefore, PDB-SLIC exhibits superior vertebrae segmentation performance compared to the other two algorithms.

Our algorithm presents a significant improvement in accuracy when compared to traditional unsupervised methods, as shown in Table 2. Our results highlight the efficacy of the PDB-SLIC algorithm, showcasing a substantial enhancement of 4.9% in the Jaccard index and 5.7% in precision relative to SLIC. Furthermore, compared to VBSeg, the PDB-SLIC algorithm demonstrates a 2.7% increase in the Jaccard index and a 3.5% improvement in precision. Notably, it even outperforms established supervised neural network algorithms such as FCN and DeeplabV3 networks. The difference in Jaccard value between the PDB-SLIC algorithm and the U-Net network is a mere 0.2%, while the precision value differs by only 0.3%. The variations in recall values are not statistically significant (F = 2.691; p = 0.127). These slight disparities compared to classical deep learning networks highlight the exceptional accuracy achieved by our algorithm, eliminating the need for manually annotated datasets by medical professionals.

For the CCL-A method, PDB-SLIC’s yielded a mean absolute error of 2.980°, with a standard deviation of 3.167°, and an intraclass correlation coefficient (ICC) of 95.8%. These results closely mirrored those obtained through manual measurements. Conversely, for the CCL-B method, PDB-SLIC exhibited a mean absolute error of 1.604°, with a standard deviation of 1.674°, and an ICC of 97.9%, surpassing the accuracy of manual measurements. These findings indicate the reliability of PDB-SLIC in producing measurement results, thereby validating the algorithm’s efficacy and acceptability.

PDB-SLIC was employed to measure the cervical curvature across the dataset, as summarized in Table 3. A one-way ANOVA, employing Welch’s method, was conducted, followed by a Games-Howell post hoc test for pairwise comparisons among groups. The results revealed a significant increases in cervical curvature with age (F = 83.097; p < 0.001) (as shown in Figure 8). While no statistical significant differences were observed among the 20–30, 31–40, 41–50, 51–60, and 61–70 age groups (p > 0.05), all remaining between-group comparisons yielded statistically significant differences (p < 0.05).

Non-parametric analysis using the Kruskal-Wallis test was employed to assess the impact of various factors, including different hospitals, MRI machines and genders, on cervical curvature across different age groups. The results showed that there were no statistically significant differences in cervical curvature among different hospitals in the 20–30 and 81–90 age groups (p > 0.05), while statistically significant differences were observed in all other age groups (p < 0.05); Similarly, in the 20–30, 31–40, and 41–50 age groups, no statistically significant differences were found in the effect of different MRI machines on cervical curvature (p > 0.05), whereas significant differences were observed in all other age groups (p < 0.05). In the 81–90 age group, sex did not have a statistically significant effect on cervical curvature (p > 0.05), while significant differences were noted in all other age groups (p < 0.05) (refer to Table 4 and Figure 9).

A strong positive correlation was observed between cervical curvature and C7 slope (r = 0.45; 95% CI [0.38, 0.52]; p < 0.001). Conversely, the correlation between cervical curvature and the curvature/C7 slope ratio was positive but weaker (r = 0.29; 95% CI [0.21, 0.37]; p < 0.001). Additionally, there was a weak negative correlation between C7 slope and the curvature/C7 slope ratio (r = −0.34; 95% CI [−0.42, −0.25]; p < 0.001) (as shown in Figure 10).