A method and device for automatically extracting a skeletal mechanical framework

Abstract

The application provides a kind of skeleton mechanical skeleton automatic extraction method and device, it is related to medical data processing technical field.The method comprises: obtaining multiple groups of bone image data, grouping according to clinical features, obtaining statistical average bone section;On each average bone section, extract the initial bone region boundary, sequentially execute the operation of internal shrinkage and the operation of differential translation based on anatomical partition, generate the optimized region of interest boundary;In the optimized region of interest, the high-density bone sub-region is extracted using the adaptive threshold method, and the corresponding binary mask is generated;Skeletonization processing is performed on the binary mask, and the centerline network is extracted as the skeleton mechanical skeleton of the section;The topological characteristics of the skeleton mechanical skeleton are calculated, and the bone density distribution in the region of interest is combined to evaluate the mechanical properties of the skeleton.By statistical average modeling independent of anatomical feature points, individual, equipment and parameter differences are eliminated, standardized bone sections are constructed, and the degree of automation and robustness is significantly improved.

Description

Technical Field

[0001] This invention relates to the field of medical data processing technology, and in particular to an automatic method and apparatus for extracting skeletal biomechanics. Background Technology

[0002] Skeleton is the primary biomechanical component of the human musculoskeletal system. Describing the spatial structure of the main load-bearing components of bone helps to understand the distribution and transmission of forces in space, providing a theoretical basis for osteoporosis assessment, fracture risk prediction, and personalized implant design. In recent years, data mining based on Hounsfield Units (HUs) in computed tomography (CT) images has become one of the core methods for non-invasive assessment of local bone quality. Multiple studies have confirmed a significant positive correlation between HUs and bone mineral density, enabling spatial modeling of bone microstructure through quantitative analysis.

[0003] However, existing methods for skeletal biomechanics analysis and skeleton extraction still have the following technical limitations: Existing averaging model methods rely on predefined anatomical feature points (such as the femoral bulb center, the cervical isthmus center, etc.) for registration. When there are anatomical variations, pathological changes, or poor image quality in the bones, feature point recognition is difficult or fails, resulting in low automation and poor robustness in model construction. Summary of the Invention

[0004] To address the technical problem that existing averaging models rely on predefined anatomical feature points (such as the femoral bulb center and the cervical isthmus center) for registration, which leads to difficulties or failures in feature point recognition when there are anatomical variations, pathological changes, or poor image quality, resulting in low automation and poor robustness in model construction, this invention provides an automatic skeletal biomechanical skeleton extraction method and apparatus. The technical solution is as follows: On the one hand, an automatic method for extracting the skeletal biomechanics skeleton is provided, including the following steps: Step S1: Acquire multiple sets of bone image data, group them according to clinical characteristics, select multiple two-dimensional cross sections along the long axis of the bone in each set of data, acquire grayscale information of multiple sampling positions on each cross section, and generate the statistical average bone cross section corresponding to each cross section based on the statistical average of multiple sets of data. Step S2: Extract the initial bone region boundary on each average bone section, and sequentially perform the inward shrinkage operation and the differential translation operation based on anatomical partitions to generate the optimized region of interest boundary; Step S3: Within the optimized region of interest, a high-density sub-region is extracted using an adaptive thresholding method to generate a corresponding binary mask. The adaptive thresholding method determines the threshold based on the density distribution within the current region of interest. Step S4: Perform skeletonization processing on the binary mask to extract the centerline network as the skeletal mechanical skeleton of the cross section. The skeletal mechanical skeleton reflects the topological connectivity of the load-bearing trabecular network. Step S5: Calculate the topological features of the skeletal mechanical framework and evaluate the mechanical properties of the bone by combining the bone density distribution in the region of interest.

[0005] Optionally, the step S1 of generating the statistical average bone cross section corresponding to each cross section based on the statistical average of multiple sets of data includes: The geometric coordinates of each sampling location are statistically averaged, and the grayscale information of each sampling location is interpolated and standardized before calculating the statistical average grayscale.

[0006] Optionally, the clinical characteristics include age and sex; The plurality of two-dimensional sections include a lower section, a middle section, and an upper section distributed along the long axis of the bone; The process of obtaining grayscale information from multiple sampling locations includes drawing multiple equally spaced sampling lines on each cross section, collecting the grayscale values ​​of pixels on the sampling lines, and performing geometric correction on the sampling coordinates.

[0007] Optionally, the retraction operation includes: Iteratively move the boundary vertices inward with a preset step size, and calculate the image feature value at the current vertex after each iteration; When the image feature value meets the preset condition, the vertex stops moving; The image feature values ​​include image gradient intensity and / or density values.

[0008] Optionally, the translation operation includes: classifying the recessed boundary vertices into at least two regions according to their anatomical positions, and applying different translation directions and / or translation amounts to each region; the translation directions include inward translation, outward translation, and / or outward expansion.

[0009] Optionally, the anatomical location classification is based on the coordinate values ​​of the boundary vertices and is divided into left-side region, middle region, and right-side region; The vertices of the left region shift to the right, the vertices of the right region shift to the left, and the vertices of the middle region expand outward along the ray direction; The coordinate range and translation amount of each region are determined based on the statistical average of the skeletal anatomical divisions.

[0010] Optionally, the translation operation includes: The translation amount is determined based on the distance between the vertex and its adjacent vertices or the distance to the boundary, and each has a preset lower threshold or upper threshold. After translation, the number of vertices of all polygons in the region of interest is unified to a preset number.

[0011] Optionally, the skeletonization process in step S4 includes: extracting the morphological skeleton of the region to obtain a centerline network with a single pixel width.

[0012] Optionally, it also includes: comparing the extracted skeletal features with a reference database of the same age and sex to generate a skeletal biomechanical status assessment report; and / or, combining the bone mineral density histogram distribution with skeletal topology indicators to construct an osteoporosis risk prediction model.

[0013] On the one hand, an automatic skeletal biomechanical skeleton extraction device is provided, comprising: Data input module: used to receive multiple sets of skeletal image data; Preprocessing module: used to perform section selection, sampling, correction, and average section generation; ROI Extraction and Optimization Module: Used to perform boundary extraction and polygon optimization; Skeleton extraction and analysis module: used to perform skeleton generation and feature calculation; Output and Reporting Module: Used to generate visualizations and mechanical evaluation reports.

[0014] The beneficial effects of the technical solutions provided in the embodiments of the present invention include at least the following: The method provided by this invention generates a statistically average bone cross-section based on the statistical average of multiple sets of data, eliminating data differences caused by different individuals, different CT equipment, and different scanning parameters, and constructing a statistically significant standard bone cross-section model (ABS). This makes subsequent analysis independent of the randomness of specific samples, significantly improving the standardization of the skeleton extraction process and the reproducibility of the results. This invention employs a statistical averaging method based on a large number of samples. The coordinates and grayscale values ​​of each sampling location are derived from the arithmetic mean of multiple sets of data. Anatomical variations or pathological changes in individual samples do not significantly affect the overall morphology and density distribution of the average cross-section, thus ensuring the stability and representativeness of the model. Grayscale information of all pixels is collected along the sampling line, and random noise and individual differences are suppressed by averaging the grayscale values ​​of multiple sets of data. Even if a single image has local quality defects, the statistical averaging process can effectively reduce its impact, improving the robustness of model construction. This method significantly improves the automation and robustness of average bone cross-section generation, enabling it to adapt to a wider range of clinical real-world data. Attached Figure Description

[0015] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1 This is a flowchart of an automatic skeletal biomechanical skeleton extraction method provided in an embodiment of the present invention; Figure 2 This is a data extraction diagram of an automatic skeletal biomechanics skeleton extraction method provided in an embodiment of the present invention; Figure 3 This is a grayscale extraction and prediction image of an automatic skeletal biomechanical skeleton extraction method provided in an embodiment of the present invention; Figure 4 This is a schematic diagram of the ROI on a male ABS (anatomical skeleton) according to an embodiment of the present invention for an automatic extraction method of skeletal biomechanics skeleton; Figure 5 This is a cross-sectional morphological analysis diagram of an automatic skeletal biomechanics extraction method provided in an embodiment of the present invention; Figure 6 This is a block diagram of an automatic skeletal biomechanical skeleton extraction device provided in an embodiment of the present invention. Detailed Implementation

[0017] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0018] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0019] In this embodiment of the invention, sometimes a subscript such as W1 may be mistakenly written as a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

[0020] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0021] This invention provides an automatic method for extracting the skeletal biomechanical skeleton. Figure 1 This is a flowchart of an automatic skeletal biomechanical skeleton extraction method provided by an embodiment of the present invention. Please refer to [link / reference]. Figure 1 This includes the following steps: 101: Acquire multiple sets of bone image data, group them according to clinical characteristics, select multiple two-dimensional sections along the long axis of the bone in each set of data, acquire grayscale information at multiple sampling positions on each section, and generate the statistical average bone section corresponding to each section based on the statistical average of multiple sets of data.

[0022] Skeletal imaging data: refers to DICOM format images acquired through clinical computed tomography (CT) equipment, containing grayscale information of the bones (CT values, in HU).

[0023] Clinical characteristics: These include demographic parameters such as age and sex, used to group data to reduce the impact of individual differences on the statistical mean.

[0024] Two-dimensional cross-section: Representative positions selected along the long axis of the bone (such as the direction from the articular surface to the shaft of the distal femur), including the lower cross-section (close to the articular surface), the middle cross-section (between the shaft and the articular surface), and the upper cross-section (close to the shaft).

[0025] Sampling location: The set of pixels where equally spaced sampling lines drawn on each two-dimensional section intersect with bone tissue.

[0026] Statistical average bone section: A standardized cross-sectional model generated by statistically averaging the geometric coordinates and gray values ​​of multiple samples (e.g., 215 sets of CT data) at the same anatomical location, denoted as ABS (Average Bone Section).

[0027] For details, please see Figure 2 Multiple sets (e.g., 215 sets) of distal femoral CT images were collected and grouped by age and gender. Within each set, three sections—lower, middle, and upper—were selected along the long axis of the bone. Twenty equally spaced horizontal sampling lines were drawn on each section, and the coordinates and grayscale values ​​of each pixel on the lines were collected. Since the original number of sampling points on each line varied (8–25 points for the lower section, 38–62 points for the upper section), linear interpolation was used to unify the number of points on all sampling lines to 1000. Then, the geometric mean of the coordinates and the arithmetic mean of the grayscale values ​​were calculated for all samples at the same section and sampling line position within each set of data, generating the "average bone section" (ABS) for the three positions. The Alpha Shape algorithm was used to generate the envelope boundary of the ABS, and a uniformly distributed set of points (e.g., 10,000 points) was generated within the boundary. The grayscale value of each point was predicted using Delaunay triangulation and linear interpolation.

[0028] This step resolves the incomparability issues caused by differences in scanning parameters and individual anatomy among multi-center, multi-device CT data, and provides a standardized analytical benchmark.

[0029] 102: Extract the initial bone region boundary on each average bone section, and sequentially perform inward contraction operation and differential translation operation based on anatomical partitions to generate the optimized region of interest boundary.

[0030] Initial bone region boundary: The outer contour of the bone cross section extracted through image segmentation (such as thresholding segmentation, connected component analysis).

[0031] Inward contraction: The initial boundary is moved inward to avoid artifacts and partial volume effects at the cortical bone edge, so that the boundary stops at the true boundary between cortical and cancellous bone. In one embodiment, the stopping criterion for the inward contraction operation is based on the mean and standard deviation of the gradient intensity of all vertices of the initial boundary (i.e., the boundary before the inward contraction begins), and stops when the gradient intensity exceeds the sum of the mean and standard deviation.

[0032] Differential translation operation based on anatomical regions: According to the functional differences of different anatomical regions of the bone (such as the medial condyle, lateral condyle, and intercondylar fossa), displacements of different directions and amplitudes are applied to the boundary vertices to make the region of interest (ROI) accurately cover the core weight-bearing area of ​​cancellous bone.

[0033] Region of Interest (ROI): The bone cross-sectional area ultimately used for bone density analysis and skeleton extraction, i.e., the polygonal region enclosed by the optimized boundary.

[0034] For details, please see Figure 4 Extract the largest connected component from the ABS as the initial bone region boundary, and resample the boundary to a fixed number of vertices (e.g., 80).

[0035] This step, through a two-step strategy of inward retraction and translation, automatically and precisely focuses on the core weight-bearing area of ​​cancellous bone, eliminating the interference of cortical bone edge artifacts and significantly improving the repeatability and accuracy of ROI localization.

[0036] 103: Within the optimized region of interest, an adaptive thresholding method is used to extract high-density sub-regions and generate corresponding binary masks. The adaptive thresholding method determines the threshold based on the density distribution within the current region of interest.

[0037] Adaptive thresholding: A method that automatically determines the binarization threshold based on the density distribution within the current ROI, without relying on a fixed absolute HU value.

[0038] High-density trabecular region: The highest density set of voxels within the ROI, corresponding to the load-bearing trabecular network that bears the main load.

[0039] Binary mask: A binary image in which high-density bone regions are marked as 1 (foreground) and the remaining regions are marked as 0 (background).

[0040] Specifically, within the optimized ROI, the HU values ​​of all voxels are sorted from high to low, and the HU values ​​corresponding to the top 30% of the cumulative pixel count (i.e., the top 30% of the highest density voxels) are taken as the threshold to generate a binary mask (1 for high-density areas and 0 for the rest). This threshold does not depend on any fixed absolute HU value and is automatically determined by the density distribution of the current ROI.

[0041] This step enables fully data-driven extraction of high-density bone subregions, unaffected by HU value drift from different CT devices, and can sensitively capture the preferential loss of high-density trabeculae during the osteoporosis process.

[0042] 104: Perform skeletonization on the binary mask and extract the centerline network as the skeletal mechanical skeleton of the cross section. The skeletal mechanical skeleton reflects the topological connectivity of the load-bearing trabecular network.

[0043] Skeletonization: Compresses a binary object into a single-pixel-wide centerline network while preserving the topological connectivity of the original object.

[0044] Skeletal biomechanical framework: The topological structure reflecting the central path of the trabecular load-bearing network, distinct from the geometric central axis of the bone.

[0045] For details, please see Figure 5 Skeletonization can be achieved using morphological thinning algorithms (such as iterative boundary stripping) or ridge extraction algorithms based on distance transform, extracting a network of connected centerlines with a single-pixel width. This network constitutes the skeletal mechanical skeleton, representing the central path of high-density load-bearing trabeculae. The number of branches and their connection methods reflect the integrity and connectivity of the trabeculae.

[0046] This step generates the topological skeleton of the load-bearing trabecular bone network, rather than a simple geometric centerline, providing a direct quantitative indicator for assessing the connectivity of bone microstructures.

[0047] 105: Calculate the topological features of the skeletal biomechanics framework and evaluate the mechanical properties of the skeleton by combining the bone density distribution in the region of interest.

[0048] Topological features include quantitative indicators such as the number of branches in the skeleton, the number of endpoints, and the length distribution.

[0049] Number of branches: The number of independent paths derived from all branch points (degree ≥ 3) in the skeleton.

[0050] Endpoint count: The number of pixels with a skeleton density of 1.

[0051] Length distribution: Statistical values ​​(mean, standard deviation, maximum value, etc.) of the pixel length or physical length (millimeters) of each branch.

[0052] Bone mineral density distribution: The distribution of bone mineral density (BMD) calculated from the HU values ​​of each voxel within the ROI after phantom calibration.

[0053] Specifically, the calculation of the topological features of the mechanical skeleton, such as the number of branches, the number of endpoints, and the total length, is performed. The individual's skeletal features are compared with a normal reference database of the same age and sex, and the deviation value is calculated. A deviation exceeding a preset threshold (e.g., 1.96 times the standard deviation) is considered an abnormal bone biomechanical state. Simultaneously, the average bone mineral density and skeletal topological indicators within the ROI can be combined and input into an osteoporosis risk prediction model (e.g., logistic regression) to output the fracture risk probability.

[0054] This step combines bone mineral density with skeletal topology information, providing a more comprehensive assessment of bone biomechanical properties than a single bone mineral density measurement, which can assist in fracture risk stratification and clinical decision-making.

[0055] The method provided by this invention generates a statistically average bone cross-section based on the statistical average of multiple sets of data, eliminating data differences caused by different individuals, different CT equipment, and different scanning parameters, and constructing a statistically significant standard bone cross-section model (ABS). This makes subsequent analysis independent of the randomness of specific samples, significantly improving the standardization of the skeleton extraction process and the reproducibility of the results. This invention employs a statistical averaging method based on a large number of samples. The coordinates and grayscale values ​​of each sampling location are derived from the arithmetic mean of multiple sets of data. Anatomical variations or pathological changes in individual samples do not significantly affect the overall morphology and density distribution of the average cross-section, thus ensuring the stability and representativeness of the model. Grayscale information of all pixels is collected along the sampling line, and random noise and individual differences are suppressed by averaging the grayscale values ​​of multiple sets of data. Even if a single image has local quality defects, the statistical averaging process can effectively reduce its impact, improving the robustness of model construction. This method significantly improves the automation and robustness of average bone cross-section generation, enabling it to adapt to a wider range of clinical real-world data.

[0056] For further details, please see Figure 3 The statistical average bone cross section generated based on the statistical average of multiple sets of data in section 101 includes: The geometric coordinates of each sampling location are statistically averaged, and the grayscale information of each sampling location is interpolated and standardized before calculating the statistical average grayscale.

[0057] Wherein, geometric coordinates are the planar coordinates (x, y) of each pixel on each sampling line.

[0058] Statistical average: Calculate the arithmetic mean of the coordinates of the same location (same cross section, same sampling line index) of multiple sets of data.

[0059] formula: ,in Let be the coordinate sequence of the k-th sampling line in the m-th data set.

[0060] Interpolation standardization: Since the number of pixels on each original sampling line is inconsistent (8-25 in the lower section, 38-62 in the upper section), it is impossible to directly average them. Linear interpolation is used to unify the number of pixels on all sampling lines to a fixed length (e.g., 1000 points), so that the grayscale sequences between different samples and different sections can be matched point by point.

[0061] Statistical mean gray level: For the standardized gray level sequence, calculate the arithmetic mean of the gray levels of all samples at the same sampling line position and the same interpolation point index. .

[0062] In this step, for each sampling line of each section in each set of data, the original pixel coordinate sequence and grayscale sequence are recorded.

[0063] For each sampling line's coordinate sequence, sort them in spatial order. Then, for sampling lines with the same index across all samples, calculate the average coordinates point by point to obtain the average coordinate sequence.

[0064] For each sampling line's grayscale sequence, determine the minimum and maximum values ​​of the original coordinates. Then, uniformly insert 1000 interpolation points within the interval. For each interpolation point, use linear interpolation of the grayscale values ​​of adjacent original points to obtain the grayscale value of that point. This generates a grayscale sequence of length 1000 for each sampling line.

[0065] For all samples with the same cross section, the same sampling line index, and the same interpolation point position (the p-th interpolation point), calculate the arithmetic mean to obtain the average gray value at that position.

[0066] In this step, geometric coordinate averaging eliminates the rotation and translational deviations of bone cross-sections caused by different patient scanning positions, aligning the anatomical positions of all samples. Interpolation standardization resolves the inconsistency in the number of sampling points caused by different bone cross-sectional widths, providing a data foundation for subsequent point-by-point statistical averaging. Gray-level averaging suppresses random noise in the image and physiological differences in bone density between individuals, generating an "average gray-level cross-section" that represents the typical bone density distribution pattern of this population.

[0067] In one embodiment of the present invention, clinical characteristics include age and sex; Multiple two-dimensional sections include a lower section, a middle section, and an upper section distributed along the long axis of the skeleton; Obtaining grayscale information from multiple sampling locations involves drawing multiple equally spaced sampling lines on each cross section, collecting the grayscale values ​​of pixels on the sampling lines, and performing geometric correction on the sampling coordinates.

[0068] The age and sex mentioned above are used to divide the samples into different subgroups (e.g., men aged 50-60, women aged 60-70, etc.), because bone mineral density and bone microstructure vary significantly with age and sex. Constructing the average bone cross-section for specific subgroups can improve the personalization accuracy of the assessment.

[0069] Lower, middle, and upper sections: These are three key locations selected along the long axis of the distal femur, from the articular surface (lowest point) towards the diaphysis (highest point). The lower section is close to the articular surface, has a small cross-sectional area, and an irregular shape; the middle section is located between the diaphysis and the articular surface; and the upper section is close to the diaphysis, has a large cross-sectional area, and is nearly circular. Selecting these three sections can indirectly reflect the three-dimensional morphological changes of the skeleton.

[0070] Equally spaced sampling lines: A set of horizontal (or vertical) lines drawn at equal intervals from top to bottom or from left to right on a two-dimensional cross-sectional image, with a spacing of 5-10 pixels, to ensure that the sampling covers the entire skeletal cross-section.

[0071] Geometric correction: Because the patient's limbs may be tilted during CT scans, the cross-sections are not strictly horizontal or vertical. By calculating the average slope of the sampling lines, a rotation matrix is ​​used to rotate the coordinates to the horizontal direction, and then the coordinate system is translated to unify the origin of the coordinate system to the upper left corner of the image.

[0072] In this step, the collected CT data are divided into multiple subgroups based on age (e.g., every 10 years) and sex. Each subgroup undergoes independent processing. Within each data group, for each distal femur sample, the distal end (lowest point of the articular surface) and the superior end (a point on the diaphysis) of the femur are first located, and two-dimensional images are cropped at three heights along the long axis: lower, middle, and upper. On each cross-section, starting from the upper boundary of the image, 20 horizontal lines are drawn at fixed steps (e.g., 10 pixels), covering the entire bone cross-section. For each horizontal line, the coordinates and grayscale values ​​of all pixels intersecting the line with the bone region (pixels identified as bone tissue through thresholding) are extracted. For the coordinate sequence extracted from each line, the slope of the fitted straight line is calculated, and all points are rotated to the horizontal using a rotation matrix, then translated to the origin of a unified coordinate system (e.g., the x-coordinate of the leftmost point is set to zero).

[0073] The age and sex groupings described above ensure that the subsequently constructed ABS (analytical bone sample) is representative of a specific population and can be used for accurate comparison between individuals and normal reference values. The selection of multiple cross-sections overcomes the limitation of two-dimensional analysis in reflecting three-dimensional morphology; the lower, middle, and upper cross-sections represent different functional areas of the skeleton (articular surfaces, the transition zone between the diaphysis and articular surfaces, and the diaphysis). Equal-interval sampling ensures the uniformity of sampling and statistical comparability. Geometric correction eliminates systematic errors introduced by scanning tilt, ensuring accurate alignment of cross-sections from different samples.

[0074] Furthermore, the indentation operation includes a variety of different methods. In one embodiment, a fixed-distance indentation is used. For example, the boundary can be offset inward by a certain distance (e.g., 5–10 pixels) based on experience.

[0075] In another embodiment, adaptive iterative shrinkage is employed, specifically: Iteratively move the boundary vertices inward with a preset step size, and calculate the image feature value at the current vertex after each iteration; The vertex stops moving when the image feature value meets the preset condition; Image feature values ​​include image gradient intensity and / or density values.

[0076] The preset step size mentioned above refers to the distance moved in each iteration, typically 1 pixel, to ensure a smooth transition of the boundary. Inward iterative movement involves gradually moving the boundary vertices along a direction pointing inwards from the polygon (e.g., the direction of the angle bisector between adjacent sides).

[0077] Image feature values: Computable quantities reflecting the local tissue characteristics of a vertex, including: Image gradient intensity: the rate of change of gray level at this location, with a sharp increase in gradient at the junction of cortical and cancellous bone.

[0078] Density value: CT value (HU) at this location, with high cortical bone density (usually >700HU) and low cancellous bone density.

[0079] Preset conditions: A set of criteria for determining whether a vertex has reached the ideal stopping position, such as gradient strength exceeding a threshold or density value exceeding a threshold.

[0080] In this step, starting from the initial bone region boundary (polygon), each boundary vertex is assigned an initial position. The following steps are repeated until all vertices stop or the maximum number of iterations is reached: for each vertex that has not yet stopped, move it 1 pixel along the shrinking direction.

[0081] At the new location, calculate the image gradient intensity (e.g., Sobel operator magnitude) and / or CT value for that point. Compare the calculated feature values ​​with preset conditions. If the conditions are met (e.g., gradient intensity > 50), the vertex is marked as "stopped" and will not move further. When all vertices have stopped or the maximum number of iterations (e.g., 15 steps) has been reached, the shrinking operation ends, yielding the shrunken boundary.

[0082] The iterative movement method described above can finely adjust the boundary position, avoiding over- or under-adjustment caused by improper one-time inward distance. It automatically identifies the true boundary between cortical and cancellous bone using image features (gradient, density), replacing fixed-distance inward adjustment, and is robust to skeletal anatomical variations and image artifacts. Different vertices can stop at different positions due to local tissue differences, making the final boundary more closely resemble the actual bone tissue boundary.

[0083] In one embodiment of the present invention, the translation operation includes: classifying the inwardly recessed boundary vertices into at least two regions according to their anatomical positions, and applying different translation directions and / or translation amounts to each region; the translation directions include inward translation, outward translation, and / or outward expansion.

[0084] Among them, the anatomical location classification is based on the functional division of the boundary vertices in the skeletal anatomy (such as the medial condyle, lateral condyle, and intercondylar fossa). The distribution of cancellous bone and the load-bearing requirements are different in different regions.

[0085] Translation direction: Lateral translation: moving towards the center line of the bone (e.g., moving the left femur to the right).

[0086] Lateral displacement: moving outwards from the bone (e.g., moving the left femur to the left).

[0087] Expand outward: Move outward along the line connecting the vertex and the center of the bone cross section to increase the area of ​​the region.

[0088] Translation amount: The distance each vertex moves, which can be a fixed value, a dynamic value based on local geometric features, or a value based on statistical priors.

[0089] In this step, the boundary polygon after the inward retraction is obtained. Each vertex is assigned to an anatomical region based on its coordinates (or position relative to anatomical landmarks). For example, the distal femur is divided into the medial condyle, lateral condyle, and intercondylar fossa.

[0090] Preset translation strategy for each partition: Medial condyle apex: translated medially (so that the ROI covers the medial cancellous bone core area).

[0091] Lateral condyle vertex: translated outwards.

[0092] The apex of the intercondylar fossa extends outward (expanding its coverage area).

[0093] Calculate the translation vector for each vertex according to the strategy and move the vertex to obtain the translated boundary.

[0094] In this step, differential translation based on anatomical regions allows the ROI to accurately match the distribution morphology of cancellous bone in different areas, avoiding incomplete coverage or excessive cortical bone coverage in some areas due to uniform translation. By extending the intercondylar fossa region outward, the trabecular bone structure in this region can be included, improving the integrity of the skeletal extraction. The adjustability of the translation direction and amount allows this method to adapt to the anatomical characteristics of different bones (such as the femur and tibia).

[0095] Furthermore, the anatomical location classification is based on the coordinate values ​​of the boundary vertices, and is divided into the left side region, the middle region, and the right side region; The vertices of the left region shift to the right, the vertices of the right region shift to the left, and the vertices of the middle region expand outward along the ray direction; The coordinate range and translation amount of each region are determined based on the statistical average of the skeletal anatomical regions.

[0096] The left, middle, and right regions correspond to the projections of the medial condyle, intercondylar fossa, and lateral condyle of the distal femur. On the two-dimensional cross-section, the X-coordinates, from smallest to largest, correspond to the movement from medial to lateral. Rightward translation: movement in the direction of increasing X-coordinate. Leftward translation: movement in the direction of decreasing X-coordinate. The ray direction is from the vertex to the geometric center of the cross-section. Outward expansion along the ray direction: rays emanating from the center of the bone cross-section (e.g., the centroid) towards the vertex, with the vertex moving outward along the ray. Statistical average: empirical values ​​of typical coordinate ranges and translations for each region calculated through measurements of a large number of distal femoral CT images from healthy individuals.

[0097] In this step, the coordinates (xi, yi) of the shrunken boundary vertices are obtained. The region of each vertex is determined based on a preset X-coordinate threshold. If xi ≤ X1, then it belongs to the left region; If X1 < xi < X2, then it belongs to the middle region; If xi ≥ X2, then it belongs to the right-hand region.

[0098] X1 and X2 are determined by statistical averaging of the boundaries of the medial and lateral femoral condyles in the training samples (e.g., X1 = 120 pixels, X2 = 470 pixels).

[0099] Left region vertex: Move this vertex to the right by the greater of the first horizontal distance and 35 pixels. The first horizontal distance refers to the minimum horizontal distance between this vertex and its left and right adjacent vertices. If adjacent vertices are very close (horizontal distance less than 35), force a movement of at least 35 pixels to ensure adequate coverage of the medial cancellous bone; if adjacent vertices are already far apart (greater than 35), move them by the actual distance to maintain the local shape.

[0100] Right-side region vertex: Move this vertex to the left by the larger of the second horizontal distance and 50 pixels. The second horizontal distance refers to the minimum horizontal distance between this vertex and its left and right adjacent vertices.

[0101] Intermediate region vertex: Move this vertex outward along the ray by the smaller of the distance from the vertex to the boundary and 75 pixels. Intermediate regions (such as the intercondylar fossa) need to be expanded outward to include more trabeculae, but not too much to avoid entering the cortical bone or exceeding the bone boundary. As a rule of thumb, limiting the maximum to 75 pixels is a safe upper limit.

[0102] After moving all vertices, we obtain the translated boundary.

[0103] For details, please refer to Table 1: Table 1 Translation Operations

[0104] This coordinate-based classification method is simple and effective, requiring no identification of complex anatomical landmarks and offering a high degree of automation. The translation amount combines local geometric features (distances between adjacent vertices and distances to boundaries) with empirical statistical lower limits (35, 75, 50), ensuring individual adaptability while avoiding over- or under-translation. This approach has been validated on distal femoral data, enabling precise alignment of the ROI with the cancellous bone distribution area.

[0105] Furthermore, translation operations include: The translation amount is determined based on the distance between the vertex and its adjacent vertices or the distance to the boundary, and each has a preset lower threshold or upper threshold. After translation, the number of vertices of all polygons in the region of interest is unified to a preset number.

[0106] Distance between adjacent vertices: For the left / right region, the translation amount is referenced to the horizontal distance between adjacent vertices on the same horizontal line. Specifically, it is the absolute value of the horizontal distance between the vertex and its left (or right) adjacent vertex, taken as the maximum value among all adjacent pairs (as in formula (3-1)). ).

[0107] Distance to boundary: For intermediate regions, the translation amount is referenced to the distance from the vertex to the current ROI boundary along the ray direction (i.e., the maximum space that the vertex can expand in that direction).

[0108] Preset lower threshold: the left area should be shifted by at least 35 pixels and the right area by at least 50 pixels to ensure coverage of the core cancellous bone.

[0109] Preset upper limit threshold: The expansion amount of the middle area shall not exceed 75 pixels to prevent excessive expansion into the cortical bone.

[0110] The vertex count is standardized to a preset number: the number of vertices in the translated polygon may change due to vertex movement (such as some vertices overlapping). By uniform downsampling (if there are too many vertices) or edge interpolation (if there are too few vertices), the vertex count of all polygons is standardized to 60. This ensures that all ROIs have the same parameterized representation, facilitating subsequent batch processing.

[0111] In this embodiment, for each vertex in the left region, the horizontal distance between it and its adjacent vertex on the left (with the same or similar y-coordinate) is calculated, and the minimum of all such distances is taken as the reference distance. The translation amount is the larger of this reference distance and 35 pixels.

[0112] For the vertices in the right region, similarly calculate the minimum horizontal distance to the adjacent vertices on the right, and take the larger of this minimum value and 50 pixels for the translation.

[0113] For a vertex in the middle region, calculate the distance from that vertex along the ray direction to the intersection point of the current ROI boundary (i.e., the pixel length traversed by moving along the ray from the vertex to the boundary). The extension is taken as the smaller of this distance and 75 pixels.

[0114] After performing the translation, a new polygon is obtained. Count the number of vertices N_new.

[0115] If N_new > 60, then select 60 indices evenly from the original vertex sequence (e.g., using the linspace function), keep these vertices, and delete the rest.

[0116] If N_new < 60, then calculate the number of vertices that need to be added on each edge, and add vertices at equal intervals on the edges to make the total number reach 60.

[0117] The adaptive calculation of the translation amount in this step (based on the distance between adjacent vertices or the distance to the boundary) ensures that the translation amplitude conforms to local geometric features, avoiding boundary deformation caused by rigid fixed values. Lower / upper thresholds provide a safety guarantee, ensuring that the translation is not too small (ineffective) or too large (exceeding the bone tissue range). A uniform vertex count of 60 ensures that all samples have the same ROI representation, facilitating subsequent batch feature extraction and statistical analysis.

[0118] Furthermore, the skeletonization process in step 104 includes: extracting the morphological skeleton of the region to obtain a centerline network with a single pixel width.

[0119] Binarization: Voxels in the high-density bone region are labeled as foreground (1), and the rest are labeled as background (0) to generate a binary image.

[0120] Morphological skeletonization: By performing repeated boundary stripping (thinning) operations on a binary image, the topological connectivity of the original objects is preserved, and regions are compressed into lines with a width of one pixel. Common algorithms include Hilditch thinning and Zhang-Suen thinning.

[0121] Centerline network with a single pixel width: Each pixel in the skeleton has at most two adjacent pixels (except for endpoints), and the entire skeleton is composed of curve segments and branch points.

[0122] Specifically, within the optimized ROI, a binary mask (from 103) of the high-density skeleton region is obtained. A morphological thinning algorithm is applied to this binary mask, such as iteratively deleting boundary pixels while retaining those points whose deletion would alter connectivity, until no further deletions are possible. This results in a connected network composed of single pixels. Topological analysis is performed on the skeleton image: each skeleton pixel is scanned, and the number of foreground pixels in its 8-neighborhood is counted to determine whether the pixel is an endpoint (neighborhood = 1), a connector (neighborhood = 2), or a branch (neighborhood ≥ 3). A depth-first search (DFS) is used to traverse the skeleton, identifying all branches (segments between two branch points or between an endpoint and a branch point), and recording the length of each branch (the physical length is obtained by multiplying the number of pixels by the pixel spacing). The number of branches, the number of endpoints, and the length distribution (e.g., mean, standard deviation, maximum length) are output.

[0123] In this step, morphological skeleton extraction preserves the original connectivity of the trabecular bone network, reflecting the integrity of load-bearing pathways (e.g., whether branches are missing or broken). The number of branches and endpoints are sensitive indicators of bone microstructural degradation: in osteoporosis, high-density trabeculae gradually disappear, branches decrease, and endpoints increase. Length distribution can assess the average length and continuity of load-bearing trabeculae, providing a quantitative basis for fracture risk prediction.

[0124] Furthermore, the method also includes: comparing the extracted skeletal features with a reference database of the same age and gender to generate a skeletal biomechanical status assessment report; and / or, combining the bone mineral density histogram distribution with skeletal topology indicators to construct an osteoporosis risk prediction model.

[0125] The reference database consists of statistical reference values ​​constructed from the skeletal topological features extracted by this method from skeletal images of a large number of healthy individuals (different ages and genders), including the mean and standard deviation of the number of branches, the number of endpoints, and the total length of the skeleton for each age and gender group.

[0126] Comparison: Compare the skeletal feature values ​​of the individual to be evaluated with the mean values ​​of the corresponding age and gender groups in the reference database, and calculate the deviation (e.g., Z-score = (individual value - mean) / standard deviation).

[0127] Skeletal biomechanical status assessment report: includes individual skeletal characteristic values, percentage deviation from normal range, and judgment conclusions based on preset thresholds (such as Z-score < -1.96) (such as "reduced trabecular network connectivity, further assessment of fracture risk is recommended").

[0128] Bone mineral density histogram distribution: A statistical histogram of the HU values ​​(or converted BMD) of all voxels within the ROI, reflecting the overall distribution of bone mineral density (such as peak value, skewness, kurtosis, etc.).

[0129] Osteoporosis risk prediction model: Employs machine learning (such as logistic regression, random forest) or multivariate regression models. Input features include skeletal topology indicators (number of branches, number of endpoints, total length) and bone mineral density statistics (mean BMD, peak BMD position, etc.). Outputs osteoporosis risk probability or fracture risk score.

[0130] In this step: Distal femoral CT data were collected from at least several hundred healthy volunteers (covering all ages and sexes). The skeletal topological features and bone mineral density histograms for each case were extracted using the methods described above, thus constructing a reference database. The mean, standard deviation, percentile, and other statistics were calculated for each age-sex group. Skeletal features and bone mineral density distribution were extracted for the individuals to be evaluated using the same methods.

[0131] Calculate the Z-score of an individual's skeletal features relative to the mean of the same age and sex population. For example, if the Z-score for the number of branches is below -1.96 (P<0.05, bilateral), it is judged as "significant degradation of the trabecular network"; if it is between -1.96 and -1, it is "marginal abnormality"; if it is above -1, it is "normal".

[0132] Generate an evaluation report, outputting it in text and chart formats (such as skeleton visualizations and comparison curves with the reference range).

[0133] To construct a risk prediction model: use the skeletal features and bone mineral density features of the training set as input, and known osteoporosis diagnoses (verified by DXA bone mineral density T-score ≤ -2.5) or fracture outcomes as labels to train the classification model. After obtaining the model, for a new individual, substitute its feature values ​​to output the probability of disease.

[0134] This step, comparing with a reference database of the same age and sex, eliminates the natural influence of age and sex on bone microstructure, making the assessment more personalized. Statistical bias assessment provides objective criteria for abnormality detection, avoiding subjective thresholds. Combining bone mineral density (BMD) and skeletal topology information to build a predictive model is more accurate in assessing fracture risk than using BMD alone, because the effect of trabecular connectivity on bone strength is independent of BMD.

[0135] The following describes a specific embodiment provided by the present invention: 201. Data Acquisition and Data Preprocessing In the quantitative analysis of tibial anatomy, data preprocessing is a crucial step in ensuring the reliability of the results, and its quality directly determines the accuracy and clinical applicability of subsequent bone mineral density distribution modeling. This study extracted tibial anatomical features from 215 sets of clinical CT scan data, covering different gender groups aged 50 to 80 years, stratified into 50-60, 60-70, 70-80, and over 80 years old. The specific preprocessing workflow is as follows: First, the distal femur region was located from the aforementioned CT scan data, and three sections—lower, middle, and upper—representing different anatomical positions were selected. For each section, 20 equally spaced horizontal lines were drawn as sampling paths, and the HU (Hounsfield Unit) values ​​of the sections were collected along these lines. Experimental data showed a significant inconsistency in the number of sampling points along each horizontal line: Lower section (lowest point of the distal femur): the number of sampling points on each line is 8 to 25; Mid-section (middle of the distal femur): 23 to 40 sampling points on each line; Top section (highest position of the distal femur): The number of sampling points on each line is 38 to 62.

[0136] This fluctuation in the number of sampling points along a single line results in significant differences in the total number of sampling points (the sum of sampling points across 20 lines) across different cross-sections: the lower cross-section has approximately 160–500 sampling points (8×20–25×20), the middle cross-section has approximately 460–800 sampling points (23×20–40×20), and the upper cross-section has approximately 760–1240 sampling points (38×20–62×20). This imbalance in the total number of sampling points across cross-sections may affect the accuracy and clinical applicability of subsequent bone mineral density distribution modeling, necessitating optimization strategies to improve preprocessing quality.

[0137] 202. Generation of average bone cross section

[0138] 1) Scanning and sorting

[0139] Given Group CT scan data ( Each group contains There are three sections (bottom, middle, and top), each section has... There are 10 sampling lines. The coordinate sequence of each sampling line is as follows: ( (Original number of pixels), grayscale value sequence is The goal is to generate an average bone cross-section in a unified coordinate system that satisfies: Geometric consistency: All cross-sectional coordinate systems are aligned with the tibial long axis (eliminating tilt deviation). Statistical significance: Mean cross-section based on statistical cumulative data from multiple sets. Point standardization: Each line is standardized to 1000 equally spaced points. To eliminate cross-sectional offsets caused by differences in scanning angle or equipment, the coordinates of each sampling line must be strictly sorted according to the cross-sectional position. Specific steps: Cross-sectional location encoding: 60 sampling lines The code based on section type is: 1 20 is the lower section, 21 40 is the mid-section, 41 60 is the upper section; Sorting operation: Number the sampling lines for each group of data. Sort the coordinates in ascending order and rearrange them synchronously: (1-1) in, This is a list of sorted sample line numbers; This is a sorted array of indices; This is a list of coordinates arranged in sorted order. This ensures that the cross-sectional order is consistent across all data groups, allowing for proper alignment of cross-sectional data between different data groups and providing a unified benchmark for subsequent averaging calculations.

[0140] 2) Calculation of geometric mean coordinates

[0141] Geometric mean is applied to the coordinates of each sampling line (position k) to eliminate geometric bias in individual data sets. Specifically: Coordinate vector definition: The coordinate vector of the k-th line This represents the coordinates of the line in the m-th data set. ).

[0142] Use the arithmetic mean to obtain the new coordinates: (1-2) in, The number of valid data sets containing the k-th line ( (If a line is missing, skip it). These are the normalized average coordinates.

[0143] 3) Gray value standardization and interpolation

[0144] Original grayscale sequence To ensure consistent processing, the following standardized interpolation is performed: Normalized coordinate axes: Map the original coordinates to the interval [1, 1000] and define equally spaced points: (1-3) Linear interpolation formula: for any The grayscale value is calculated through linear interpolation between adjacent points: (1-4) in, .when Falling on the boundary ( or When ), directly take the endpoint value.

[0145] Standardized output: Generates a grayscale sequence of length 1000, providing a uniform input for averaging calculations.

[0146] 4) Based on the standardized grayscale sequence, calculate the statistical average grayscale value of each line to suppress noise and enhance significance: (1-5) in, This is the standard grayscale sequence (length 1000) of the k-th line in the m-th data set. To be the number of valid data sets, obviously Greater than or equal to 1; represents the average grayscale distribution of the k-th line, used for subsequent envelope generation and grayscale prediction.

[0147] 203. Envelope generation and uniform distribution of point sets

[0148] 1) The Alpha Shape algorithm is used to generate the cross-sectional envelope, the core of which is to calculate the convex hull boundary of the point set: (1-5) in For shape parameters (0 Controlling the tightness of the envelope (in the experiment, take...) ).

[0149] 2) Generate N=10000 uniformly distributed points inside the envelope, achieved through bounding box mesh generation and filtering, as follows: Boundary box generation (calculating the coordinate range of the point set): (1-6) Grid resolution calculation (to approximate uniform distribution, determine the grid spacing d to satisfy the target number of points): (1-7) Filter and retain points inside the envelope: final retention The internal points are used to ensure uniformity.

[0150] 204. Grayscale Value Prediction and Visualization

[0151] For predicting grayscale values ​​of points inside the envelope, linear interpolation based on Delaunay triangulation is used: (1-8) in For point The coordinates of the centroid in the triangle (satisfying) ); This represents the grayscale value of the triangle's vertices.

[0152] 205. Bone mineral density distribution modeling and stratification analysis

[0153] A high-precision bone mineral density distribution model was constructed, and based on this, stratified analysis tailored to clinical needs was conducted. The core content is as follows: Load ABS data grouped by gender-age; Construct a standardized data structure and perform HU value range calibration; Define the anatomical ROI (Region of Interest), and generate ROI boundaries with a uniform number of vertices through image processing and polygon optimization; Extract the bone density (HU) distribution within the ROI, generate a hierarchical histogram, and save the visualization results.

[0154] 206. ROI Definition and Polygon Optimization

[0155] It can be seen that the ABS data contains 24 sections (gender (2) × section location (3) × age group (4)). After integrating the original data into a unified structure, the normalized gray values ​​[0,255] are linearly mapped to the typical HU range [226,1428] of clinical CT: (1-9) This conversion ensures that the grayscale values ​​conform to the physical characteristics of bone tissue CT scans (cortical bone HU ≈ 700–1500).

[0156] Perform the following operations for each gender-age-section combination: 1) Read and visualize PNG images; 2) Extract low grayscale regions (Igray ≤ 120) as candidate regions for bone mineral density; 3) Calculate the envelope boundary: In order to accurately define the anatomical contour of the bone section and construct the region of interest (ROI) required for subsequent analysis, this study uses a method based on image morphology and geometric resampling to automatically extract the envelope boundary of bone tissue from the original CT image and optimize it into a polygon with a standard number of vertices for subsequent density analysis and morphological evaluation.

[0157] Obtaining the largest connected component boundary: Since the input image (as shown in the figure) has a white background, while bone tissue appears as a dark, low-grayscale area, background information can be used to achieve effective segmentation: First, the image is converted into a binary image: a threshold is set to 0 for non-bone areas (i.e., white background) and 1 for bone tissue areas; then, connected component analysis is used to identify all independent bone tissue areas; finally, the connected component with the largest area is extracted as the main bone cross-sectional area to eliminate possible noise spots or small artifacts.

[0158] Boundary extraction and closure processing: a) Perform boundary tracing on the largest connected component (e.g., using algorithms such as bwboundaries or imfindcircles) to obtain its outer contour curve; b) Since the image may have slight discontinuities or jagged edges, the boundaries need to be smoothed (e.g., using the smooth function or morphological opening operation) to ensure that the contour is closed and the topological structure is reasonable.

[0159] Resample the boundary to 80 uniform vertices: Resample the original boundary curve at equal intervals to include 80 uniformly distributed vertices.

[0160] 3) ROI Construction: Shrinking and Translation

[0161] After obtaining the envelope boundary of 80 vertices, the final ROI polygon is generated through the following steps: Inflation / Offset: Offset the boundary inward by a certain distance (e.g., 5–10 pixels) to avoid potential artifacts or high gradient areas at the edge of the cortical bone and focus on the internal cancellous bone structure.

[0162] Translation: Fine-tune the position of the ROI center based on clinical needs or anatomical symmetry (e.g., shift towards areas with dense trabecular bone) to ensure that the ROI covers key bone quality areas.

[0163] The specific translation operation is as follows: The boundary vertices are divided into three categories according to their X-coordinates (as shown in Table 1 above), and their positions are adjusted accordingly to conform to the anatomical structure. The translation formula for the left-side region is as follows: (3-1) 4) Since the number of polygon vertices may change after resampling, it is necessary to forcibly unify the number of vertices of all ROI polygons to 60: Too many vertices: Uniform downsampling (linspace index selection); Too few vertices: Add vertices by interpolation on the edges.

[0164] 207. Morphological Analysis and Verification

[0165] To further evaluate the microstructural integrity and topological connectivity of bone tissue, this study introduced a series of morphological operations based on binary images. This process not only generates auxiliary visualization results but, more importantly, provides a closed-loop verification mechanism for the anatomical rationality of the Region of Interest (ROI) and extracts potential bone quality biomarkers (such as convex hull fullness and skeletal branching complexity). The entire process consists of four core steps: ROI clipping, high-density sub-region extraction, convex hull and skeletal generation, and multimodal visualization verification.

[0166] High-density bone sub-region extraction: Osteoporosis is not only manifested as a decrease in overall density, but more importantly, it is characterized by the preferential loss of high-density load-bearing structures (such as cortical bone and thick trabecular bone). Therefore, focusing on the subset with the highest HU value can more sensitively reflect microstructural degradation.

[0167] Convex Hull Analysis: The convex hull is the smallest convex polygon that contains all foreground pixels. In bone analysis, it reflects the overall "fullness" and continuity of the high-density bone structure.

[0168] Skeleton extraction: The skeleton is a central axis representation of a binary object, preserving its topological structure (connectivity, number of branches) while compressing it to a single pixel width. In trabecular bone analysis, the skeleton directly corresponds to the central path of the load-bearing network.

[0169] Multimodal visualization verification: To systematically display the morphological analysis results, this study generated a four-panel composite image for each age group (1. Original image + ROI boundary; 2. Cropped image; 3. Convex hull image; 4. Skeleton image).

[0170] This embodiment can achieve the following: (1) Fully automated process: From inputting raw CT data, coordinate correction, sampling standardization, to generating average bone cross section, optimizing ROI extraction, generating skeleton and outputting feature analysis and visualization reports, no manual intervention is required throughout the process, effectively reducing the dependence on operator experience and making it suitable for opportunistic screening and batch processing of large-scale clinical data.

[0171] (2) High robustness: By eliminating scanning angle deviation and position offset through coordinate rotation and translation correction, the modeling of multiple sets of data averaging suppresses individual differences and image noise, and by using morphological smoothing and connected component analysis to eliminate artifact interference, it can effectively resist image quality problems caused by factors such as bundle hardening effect, metal artifacts, and patient movement, and adapt to the data processing needs of different CT manufacturers and scanning parameters.

[0172] (3) Indirect representation of three-dimensional information: Although the analysis is based on two-dimensional cross-sectional data, by selecting three key cross-sections (lower, middle and upper) along the long axis of the distal femur and combining the bone density distribution and skeletal topological features of multiple cross-sections for comprehensive analysis, the three-dimensional topological characteristics and spatial distribution of bone microstructure can be indirectly reflected, thus making up for the limitations of two-dimensional analysis.

[0173] (4) Strong clinical compatibility: It can directly process routine clinical CT archive data without special scanning protocols or phantom calibration. It achieves standardized calibration through linear mapping of gray values ​​and HU values, adapts to DICOM format data widely used in clinical practice, requires no additional equipment investment, and is easy to promote and apply in medical institutions at all levels.

[0174] (5) High interpretability: The output bone density histogram, ROI overlay image, convex hull analysis image and skeleton visualization image are intuitive and clear, which can directly reflect the "connection status" of bone trabeculae, the distribution pattern of bone density and the degree of microstructure degeneration. The relevant quantitative indicators (such as the proportion of high density area and the number of skeleton branches) are easy for clinicians to understand and apply, providing clear support for diagnostic decisions.

[0175] (6) High precision: The system accurately focuses on the cancellous bone region through ROI shrinkage and intelligent translation strategy to avoid cortical bone interference; linear interpolation is used to standardize sampling points to ensure data consistency; and Delaunay triangulation is used to accurately predict gray values. Multiple steps ensure the accuracy and reliability of the analysis results.

[0176] Figure 6 This is a block diagram of an automatic skeletal biomechanical skeleton extraction device provided in an embodiment of the present invention. The present invention also provides an automatic skeletal biomechanical skeleton extraction device, which is used in an automatic skeletal biomechanical skeleton extraction method. (Refer to...) Figure 6The device includes a data input module 610, a preprocessing module 620, an ROI extraction and optimization module 630, a skeleton extraction and analysis module 640, and an output and reporting module 650. Among them: Data input module 610: Used to receive multiple sets of skeletal image data. This module serves as the input interface for the software system, supporting the reading of DICOM format CT image files and organizing the image data into data structures (such as three-dimensional arrays) that can be processed by subsequent modules.

[0177] Preprocessing module 620: Used to perform section selection, sampling, correction, and average section generation. This module is a software unit that implements function 101 in the method, specifically including: Section selection: Automatically locates the lower, middle, and upper sections based on the long axis of the bone.

[0178] Sampling: Draw equally spaced sampling lines on each cross section and extract the coordinates and grayscale values ​​of the pixels on the lines.

[0179] Correction: Perform geometric rotation and translation on the sampled coordinates to eliminate tilt error.

[0180] Average cross section generation: ABS is generated based on the statistical average of multiple sets of data, including coordinate average and grayscale interpolation average.

[0181] ROI Extraction and Optimization Module 630: This module performs boundary extraction and polygon optimization. It is a software unit that implements function 102 in the method, specifically including: Boundary extraction: Extract the largest connected component from the ABS as the initial bone region boundary.

[0182] Polygon optimization: Perform shrinkage operations (based on gradient / density adaptive stopping) and translation operations (based on anatomical partition differential movement) on the boundary, and unify the number of vertices to a preset number.

[0183] Skeleton Extraction and Analysis Module 640: Used to perform skeleton generation and feature calculation. This module is a software unit that implements functions 103 and 104 in the method, specifically including: High-density sub-region extraction (using an adaptive thresholding method to generate a binary mask).

[0184] Morphological skeletonization (generating a single-pixel wide centerline network).

[0185] Calculate topological features (number of branches, number of endpoints, length distribution).

[0186] Output and Reporting Module 650: Used to generate visualization results and mechanical evaluation reports. This module is a software unit that implements function 105 in the method, specifically including: Visualization results: Generate ROI boundary overlay map, skeleton map, bone density histogram, etc.

[0187] Assessment Report: Output a skeletal biomechanical status assessment report (including results compared with a reference database and risk predictions).

[0188] Specifically, the user drags a batch of CT images (such as a DICOM folder) into the system interface. The data input module reads and parses the image metadata (pixel spacing, patient age, gender, etc.). The data is then passed to the preprocessing module, which automatically identifies the distal femur region, extracts cross-sections at preset heights (lower, middle, and upper), draws 20 equally spaced horizontal lines on each cross-section, and collects the pixel grayscale values ​​on each line. The collected data undergoes geometric correction and interpolation standardization. Then, the coordinates and grayscale values ​​of all samples in the same group (e.g., same age, gender) are statistically averaged to generate three ABS (Abstract Areas). The ROI extraction and optimization module receives the ABS, first extracting the largest connected component as the initial boundary, and resampling it to 80 vertices. Then, a shrinking operation is performed: iteratively moving the vertices, adaptively stopping based on gradient strength and HU (Human Hull Number). Next, a translation operation is performed: the vertices are divided into left, middle, and right sides, and translated to the right, outward, and left respectively (the translation amount is based on the distance between adjacent vertices and a preset threshold). Finally, the number of vertices is unified to 60. The optimized ROI boundary is output.

[0189] The skeleton extraction and analysis module applies an adaptive percentile thresholding method (taking the top 30% of the highest density voxels) within the ROI to obtain a binary mask of the high-density skeleton region. This mask is then morphologically refined to extract the skeleton. The skeleton topology is analyzed, and features such as the number of branches, endpoints, and total length are calculated. Simultaneously, the average bone mineral density and bone mineral density histogram within the ROI are calculated.

[0190] The output and reporting module overlays images of ROI boundaries and skeletons onto the original ABS, generating a visual comparison chart. Simultaneously, it compares the individual's skeletal features with a built-in reference database (pre-existing in the system), generating a PDF assessment report including: individual feature values, the degree of deviation from the normal range, and an osteoporosis risk score based on a logistic regression model.

[0191] This system achieves a fully automated workflow from raw CT images to bone mechanics assessment reports, requiring no manual intervention and suitable for opportunistic screening of large-scale clinical data. Its modular design allows for independent development and maintenance of each function, facilitating future upgrades (such as replacing different skeletonization algorithms or updating the reference database). It can directly interface with hospital PACS systems, automatically capturing patient CT data, generating reports, and assisting radiologists and orthopedic surgeons in diagnosing osteoporosis and stratifying fracture risk.

[0192] It should be understood that the term "and / or" in this article is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. A and B can be singular or plural. Additionally, the character " / " in this article generally indicates an "or" relationship between the preceding and following related objects, but it can also represent an "and / or" relationship. Please refer to the context for a more accurate understanding.

[0193] In this invention, "at least one" means one or more, and "more than one" means two or more. "At least one of the following" or similar expressions refer to any combination of these items, including any combination of a single item or a plurality of items. For example, at least one of a, b, or c can represent: a, b, c, ab, ac, bc, or abc, where a, b, and c can be a single item or multiple items.

[0194] It should be understood that, in various embodiments of the present invention, the order of the above-mentioned process numbers does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present invention.

[0195] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. An automatic method for extracting a skeletal mechanical skeleton, characterized in that, Includes the following steps: Step S1: Acquire multiple sets of bone image data, group them according to clinical characteristics, select multiple two-dimensional cross sections along the long axis of the bone in each set of data, acquire grayscale information of multiple sampling positions on each cross section, and generate the statistical average bone cross section corresponding to each cross section based on the statistical average of multiple sets of data. Step S2: Extract the initial bone region boundary on each average bone section, and sequentially perform the inward shrinkage operation and the differential translation operation based on anatomical partitions to generate the optimized region of interest boundary; Step S3: Within the optimized region of interest, a high-density sub-region is extracted using an adaptive thresholding method to generate a corresponding binary mask. The adaptive thresholding method determines the threshold based on the density distribution within the current region of interest. Step S4: Perform skeletonization processing on the binary mask to extract the centerline network as the skeletal mechanical skeleton of the cross section. The skeletal mechanical skeleton reflects the topological connectivity of the load-bearing trabecular network. Step S5: Calculate the topological features of the skeletal mechanical framework and evaluate the mechanical properties of the bone by combining the bone density distribution in the region of interest.

2. The method according to claim 1, characterized in that, The step S1, which generates the statistical average bone cross section corresponding to each cross section based on the statistical average of multiple sets of data, includes: The geometric coordinates of each sampling location are statistically averaged, and the grayscale information of each sampling location is interpolated and standardized before calculating the statistical average grayscale.

3. The method according to claim 1, characterized in that, The clinical characteristics include age and sex; The plurality of two-dimensional sections include a lower section, a middle section, and an upper section distributed along the long axis of the bone; The process of obtaining grayscale information from multiple sampling locations includes drawing multiple equally spaced sampling lines on each cross section, collecting the grayscale values ​​of pixels on the sampling lines, and performing geometric correction on the sampling coordinates.

4. The method according to claim 1, characterized in that, The inward retraction operation includes: Iteratively move the boundary vertices inward with a preset step size, and calculate the image feature value at the current vertex after each iteration; When the image feature value meets the preset condition, the vertex stops moving; The image feature values ​​include image gradient intensity and / or density values.

5. The method according to claim 1, characterized in that, The translation operation includes: classifying the folded boundary vertices into at least two regions according to their anatomical positions, and applying different translation directions and / or translation amounts to each region; the translation directions include inward translation, outward translation, and / or outward expansion.

6. The method according to claim 5, characterized in that, The anatomical location classification is based on the coordinate values ​​of the boundary vertices and is divided into left, middle, and right regions. The vertices of the left region shift to the right, the vertices of the right region shift to the left, and the vertices of the middle region expand outward along the ray direction; The coordinate range and translation amount of each region are determined based on the statistical average of the skeletal anatomical divisions.

7. The method according to claim 6, characterized in that, The translation operation includes: The translation amount is determined based on the distance between the vertex and its adjacent vertices or the distance to the boundary, and each has a preset lower threshold or upper threshold. After translation, the number of vertices of all polygons in the region of interest is unified to a preset number.

8. The method according to claim 1, characterized in that, The skeletonization process in step S4 includes: extracting the morphological skeleton of the region to obtain a centerline network with a single pixel width.

9. The method according to claim 1, characterized in that, Also includes: The extracted skeletal features are compared with a reference database of the same age and gender to generate a skeletal biomechanical status assessment report; and / or, by combining the bone mineral density histogram distribution with skeletal topology indicators, an osteoporosis risk prediction model is constructed.

10. An automatic skeletal mechanical skeleton extraction device, characterized in that, include: Data input module: used to receive multiple sets of skeletal image data; Preprocessing module: used to perform section selection, sampling, correction, and average section generation; ROI Extraction and Optimization Module: Used to perform boundary extraction and polygon optimization; Skeleton extraction and analysis module: used to perform skeleton generation and feature calculation; Output and Reporting Module: Used to generate visualizations and mechanical evaluation reports.

Images