Automatic analysis method for CT images of femurs of rats and mice

Abstract

This application discloses an automatic analysis method for CT images of long bones in mice and rats, relating to the field of image processing technology. The method includes: acquiring a CT image sequence of the target long bone, extracting the three-dimensional spatial distribution of bone tissue, and performing pose alignment and three-dimensional cropping to obtain a standardized bone image sequence; then constructing candidate depth intervals based on butterfly-shaped anatomical features and key points of growth plate strip-shaped gaps, respectively, and obtaining a fine search range through intersection operation; subsequently identifying candidate slices of triangular fractures and determining the anatomical zero point; finally, mapping the target analysis region based on the anatomical zero point, completing cortical bone segmentation, medullary cavity segmentation, and trabecular bone extraction under closed mask constraints, and outputting quantitative parameters of bone tissue. This application achieves accurate positioning of the growth plate and unified analysis region mapping in long bone CT images, improving the accuracy of trabecular bone extraction and the consistency of bone tissue quantification results.

Description

[0001] This application relates to the field of image processing technology, and more specifically, to an automatic analysis method for CT images of long bones in mice and rats. Background Technology

[0002] CT imaging analysis of long bones in mice and rats is widely used in bone development research, bone metabolism assessment, drug efficacy validation, and characterization of bone disease models. By performing tomographic imaging on long bones such as the femur and tibia and obtaining continuous slice sequences, morphological changes in bone tissue near the growth plate, the distribution of trabecular bone in the metaphysis, and the structural relationship between cortical bone and medullary cavity can be observed, thus providing basic data for monitoring long bone development and quantitative analysis of bone microstructure. With the continuous increase in the number of experimental animal samples, automatic localization, region mapping, and parameter extraction of target long bones based on CT images have become a key technical direction for improving analysis efficiency, reducing manual workload, and enhancing the consistency of experimental results. Especially in research tasks that use the growth plate as an anatomical reference, how to accurately determine a unified reference layer from continuous slices and construct a stable target analysis region based on the reference layer is directly related to the reliability of subsequent bone tissue quantification results.

[0003] In existing technologies, long bone CT image analysis typically relies on manual experience or single-viewpoint features to determine the location of the growth plate, and then selects a fixed area near the experienced slice for bone tissue quantification. Due to differences in scan positioning, bone axis tilt, slice initiation slice, and background range among different samples, existing methods are prone to problems such as distorted cross-sectional morphology, inconsistent longitudinal slice correspondence, and mapping offset of the analysis area, making it difficult to establish a unified anatomical reference benchmark between different individuals. Furthermore, when locating based solely on a single cross-sectional morphology or a single longitudinal gap, it is easily affected by local noise, artifacts, and similar morphological layers, resulting in an excessively large search range, boundary drift, and the inclusion of false positive layers. Regarding the growth plate initiation boundary itself, existing methods typically lack a mechanism for identifying the critical layer where bone tissue transitions from a contact state to a disconnected state, making it difficult to form a stable and repeatable anatomical zero point, thus affecting the consistency of subsequent target analysis areas. Simultaneously, in the trabecular bone extraction stage, without the constraint of a closed space within the medullary cavity, high-grayscale boundaries and local noise in cortical bone are easily misidentified as trabecular bone structures, leading to a decrease in the accuracy and comparability of bone tissue quantification parameters.

[0004] Therefore, there is an urgent need for an automated analysis method that can achieve standardized processing of long bone CT images, precise localization of growth plates, unified determination of anatomical zero points, and quantification of bone tissue within closed constraint spaces, in order to solve the technical problems of inconsistent localization, unstable regional mapping, and susceptibility to interference in quantification results in existing technologies.

[0005] Application content To overcome the problems in existing long bone CT image analysis techniques, such as inconsistent scanning postures, insufficient growth plate localization accuracy, unstable anatomical reference levels, and susceptibility of trabecular bone quantification to interference from cortical bone boundaries and local noise, this application proposes an automated analysis method for long bone CT images of mice and rats. This method performs spatial standardization on the target long bone CT image sequence, combines cross-sectional and longitudinal slice anatomical feature localization to determine the fine search range corresponding to the growth plate, and further locks the anatomical zero point. Based on this, a unified target analysis region is constructed, completing cortical bone segmentation, medullary cavity segmentation, trabecular bone extraction under closed mask constraints, and outputting bone tissue quantification parameters, thereby improving the consistency and stability of analysis results among different samples. This application provides the following technical solution: an automated analysis method for long bone CT images of mice and rats, including: Acquire the target long bone CT image sequence, extract the three-dimensional spatial distribution of bone tissue, and perform pose alignment and three-dimensional cropping based on the main direction of bone tissue to obtain a standardized bone image sequence; A first candidate depth interval is constructed based on the butterfly-shaped anatomical features in the cross-sectional slices, and a second candidate depth interval is constructed based on the depth coordinates of key points in the growth plate-like gaps in the longitudinal slices. An intersection operation is performed on the first candidate depth interval and the second candidate depth interval to obtain a fine search range. Within a fine search range, candidate slice sequences of triangular fractures are identified where epiphyseal and metaphyseal bone tissues transition from a contact state to a disconnected state. The geometric gap distribution sequence of the epiphyseal and metaphyseal edges in each candidate slice of triangular fracture is calculated, and the slice layer number with the smallest variance of the geometric gap distribution sequence is determined as the anatomical zero point. The target analysis region is obtained by offsetting the anatomical zero point along the depth axis by a preset number of layers and truncating the preset number of sampling layers. The target analysis region is then segmented into cortical bone and medullary cavity. The medullary cavity segmentation results are used to form a medullary cavity closure mask. Trabecular bone tissue is extracted within the medullary cavity closure mask, and quantitative parameters of bone tissue are output.

[0006] Furthermore, methods for obtaining the target analysis region include: The anatomical zero-point layer number, as well as the offset layer number and sampling layer number pre-written into the parameter configuration table, are read. The anatomical zero-point layer number is added to the offset layer number to determine the starting layer number of the target analysis region. The starting layer number is added to the sampling layer number minus 1 to determine the ending layer number of the target analysis region. When the ending layer number exceeds the maximum layer number of the standardized skeletal image sequence, the maximum layer number is used as the ending layer number, and the sampling layer number is pushed back to determine the starting layer number. When the pushed-back starting layer number is less than 0, the starting layer number is set to 0.

[0007] Furthermore, the target analysis region is segmented into cortical bone and medullary cavity. The medullary cavity segmentation results are used to form a medullary cavity closure mask, including: Each cross-sectional slice in the target analysis region is input into the bone tissue segmentation network, which is implemented using a shared encoder and a dual segmentation output head structure. The encoder is used to extract the texture and contour features of bone tissue in the slices. The first segmentation output head is used to generate a cortical bone probability map, and the second segmentation output head is used to generate a medullary cavity probability map. Threshold binarization is performed on the cortical bone probability map and the medullary cavity probability map respectively to obtain the initial mask of cortical bone and the initial mask of medullary cavity. Connected component filtering, closing operation, hole filling and boundary smoothing are performed on the initial mask of medullary cavity to obtain the closed mask of medullary cavity. After deleting the overlapping area with the initial mask of cortical bone, boundary contraction processing is performed based on the inner boundary of cortical bone.

[0008] Furthermore, trabecular bone tissue was extracted within the closed mask of the medullary cavity, and quantitative parameters of the bone tissue were output, including: Only the pixel grayscale values ​​inside the closed mask of the medullary cavity are retained. A dual threshold fusion method is used to determine the trabecular segmentation threshold. The first threshold is the maximum inter-class variance threshold calculated based on the pixel grayscale histogram inside the closed mask. The second threshold is the sum of the mean pixel grayscale value inside the closed mask and a preset multiple multiplied by the standard deviation of pixel grayscale value inside the closed mask. The larger value between the first threshold and the second threshold is taken as the trabecular segmentation threshold. Pixels with gray values ​​greater than or equal to the trabecular segmentation threshold are marked as trabecular pixels. Trabecular pixels from each slice within the target analysis region are stacked along the depth axis to form a three-dimensional trabecular voxel set. A 26-neighborhood connected component analysis is performed on the three-dimensional trabecular voxel set. Connected components with a voxel count not less than a preset threshold are retained as the trabecular extraction results. The trabecular volume, medullary cavity volume, and trabecular volume fraction are calculated by combining the pixel spacing and depth layer spacing.

[0009] Furthermore, the method for obtaining the candidate slice sequence of the triangular fracture includes: Each cross-sectional slice within the fine search range is input into the critical feature detection neural network model, which outputs candidate target boxes, state category indexes, and state confidence scores. The state category indexes include contact state and disconnected state. Non-maximum suppression is performed on the candidate target boxes. Based on the state confidence scores and preset state thresholds, each cross-sectional slice is marked as a contact state slice, a disconnected state slice, or a missed slice. The interlayer position where the slice transitions from a contact state slice to a disconnected state slice is searched along the depth layer number from small to large. Consecutive adjacent disconnected state slices are identified as a triangular fracture candidate slice sequence.

[0010] Furthermore, methods for obtaining the geometric gap distribution sequence include: Based on the effective target boxes corresponding to each candidate slice in the triangular fracture candidate slice sequence, the local analysis region is determined. The lowest boundary point of bone tissue above the local low-density gap is extracted to form the epiphyseal edge curve, and the highest boundary point of bone tissue below the local low-density gap is extracted to form the metaphyseal edge curve. With the central gap position as the center, the preset sampling width and preset number of sampling points are set, and the edge spacing between the epiphyseal edge and the metaphyseal edge at each sampling position is calculated to form a geometric gap distribution sequence.

[0011] Furthermore, the method for constructing the first candidate depth interval includes: A preset retrieval step size is set for cross-sectional slices. Starting from the initial depth layer of the standardized skeletal image sequence, cross-sectional slices are extracted sequentially according to the preset retrieval step size. For each cross-sectional slice, size uniformity, grayscale normalization, and butterfly-shaped anatomical feature recognition are performed. Confidence ranking and non-maximum suppression are performed on the candidate detection boxes in the recognition results. The highest confidence score in the valid detection results is read. When the highest confidence score is greater than or equal to the preset judgment threshold, the corresponding depth layer number is recorded to form a butterfly-shaped feature hit layer number set. The first candidate depth interval is determined based on the minimum hit layer number, the maximum hit layer number, and the preset boundary compensation layer number.

[0012] Furthermore, the method for constructing the second candidate depth interval includes: In a standardized skeletal image sequence, longitudinal slices are reconstructed by fixing the horizontal or vertical coordinates. Keypoint identification is performed on the growth plate-like gaps in the longitudinal slices. The coordinate values ​​of keypoints with confidence greater than or equal to a preset keypoint threshold in the depth direction are extracted and mapped to depth layer numbers to form a set of keypoint depth layer numbers. The second candidate depth interval is determined based on the minimum keypoint layer number, the maximum keypoint layer number, and the preset keypoint boundary compensation layer number.

[0013] Furthermore, methods for acquiring standardized skeletal image sequences include: A three-dimensional spatial coordinate matrix is ​​constructed based on the voxel coordinates of bone tissue in the three-dimensional spatial distribution of bone tissue. The coordinates of the spatial center are determined. Covariance analysis and eigenvalue decomposition are performed on the centered three-dimensional spatial coordinate matrix to extract the principal direction vector of bone tissue. A dot product operation is performed on the principal direction vector of bone tissue and the unit direction vector of the preset depth axis, and the skew angle of the principal direction is determined by the inverse cosine operation. A cross product operation is performed on the principal direction vector of bone tissue and the unit direction vector of the preset depth axis to determine the rotation axis direction vector. A rotation matrix is ​​generated based on the rotation axis direction vector, the skew angle of the principal direction, and the Rodriguez rotation rule. Based on the rotation matrix, rotation transformation and interpolation resampling are performed on the voxel coordinates to generate pose-aligned 3D volume data. Based on the pose-aligned 3D volume data, 3D clipping is performed to obtain a standardized bone image sequence.

[0014] Furthermore, methods for generating rotation matrices based on rotation axis direction vectors, principal direction skew angles, and Rodriguez rotation rules include: The rotation axis direction vector is normalized to obtain the first, second, and third components of the unit rotation axis. Calculate the magnitude of the rotation axis direction vector; when the magnitude is less than a preset minimum threshold, set the rotation matrix as an identity matrix; when the magnitude of the rotation axis is not less than the preset minimum threshold, calculate the cosine and sine values ​​of the principal direction skew angle, and based on the first, second, and third components of the unit rotation axis, the cosine value of the principal direction skew angle, and the sine value of the principal direction skew angle, determine the elements at each position of the three-row, three-column rotation matrix according to the Rodriguez rotation rule.

[0015] Compared with existing technologies, the automatic analysis method for long bone CT images of mice and rats proposed in this application has the following technical effects and advantages: This application first addresses the issue of differences in scanning posture, bone axis tilt direction, and background range among different samples. It extracts the main direction of bone tissue based on the three-dimensional spatial distribution of bone tissue, and obtains standardized bone image sequences through posture alignment and three-dimensional cropping. This allows for the construction of subsequent cross-sectional and longitudinal slices under a unified spatial coordinate system, avoiding anatomical deformation and torsional shift caused by inconsistent slice orientation, and improving the basic consistency of subsequent localization analysis.

[0016] This application does not rely directly on a single image view to determine the location of the growth plate. Instead, it constructs a first candidate depth interval and a second candidate depth interval based on the butterfly-shaped anatomical features in the cross-sectional slice and the key point depth coordinates of the strip-shaped gaps in the growth plate in the longitudinal slice. The fine search range is obtained through intersection operation, thus forming a dual-path localization mechanism that combines cross-sectional morphological constraints and longitudinal spatial constraints. This can effectively compress the search depth range and reduce the impact of single-path recognition errors on the localization results.

[0017] This application addresses the technical problems of difficulty in accurately unifying the initial boundary of the growth plate and the easy drift of anatomical reference points between different individuals. Within a fine search range, it identifies candidate slice sequences of triangular fractures where epiphyseal and metaphyseal bone tissues transition from a contact state to a disconnected state. It also compares the variance of the geometric gap distribution sequence between the epiphyseal and metaphyseal edges in each candidate slice and determines the slice number of the candidate slice with the smallest variance as the anatomical zero point. This establishes the zero point positioning on the critical level where the bone tissue connection relationship changes abruptly, overcoming the inconsistency of reference caused by relying on experience to select slices, manually estimating boundaries, or coarse positioning over wide intervals in existing technologies. This significantly improves the stability and repeatability of the anatomical benchmark.

[0018] Meanwhile, this application also constructs a target analysis region by offsetting the anatomical zero point along the depth axis by a preset number of layers and truncating a preset number of sampling layers. Then, through a processing chain of cortical bone segmentation, medullary cavity segmentation, medullary cavity closure mask construction, and trabecular bone extraction within the closure mask, the identification of trabecular bone is strictly limited to the closed space inside the medullary cavity. This can effectively reduce the interference of high grayscale boundaries of cortical bone, local noise, artifacts, and interlayer discontinuities on the trabecular bone extraction results. Furthermore, by combining dual-threshold fusion segmentation and cross-layer connectivity consistency filtering, the authenticity of trabecular bone identification and interlayer continuity can be improved, thereby improving the stability and comparability of quantitative parameters such as trabecular bone volume, medullary cavity volume, and trabecular bone volume fraction.

[0019] In summary, this application, by constructing a complete technical chain of spatial standardization, dual candidate interval convergence, zero-point locking of the critical layer of the growth plate, and quantification of closed spatial constraints, not only solves the technical problems of inconsistent growth plate positioning, inconsistent analysis area diameter, and easy boundary interference in trabecular quantification in long bone CT images, but also provides unified, stable, and reproducible automatic analysis results for batches of mouse and rat samples. Attached Figure Description

[0020] Figure 1 The flowchart of the automatic analysis method for long bone CT images of mice and rats in this application is shown. Figure 2 This is a flowchart illustrating the method for acquiring standardized skeletal image sequences according to this application. Figure 3 Flowchart of the method for obtaining candidate slice sequences of triangular fracture; Figure 4 The flowchart shows the method for obtaining the geometric gap distribution sequence. Figure 5 This diagram illustrates the identification of butterfly-shaped anatomical features and the localization of the analysis area. Detailed Implementation

[0021] The technical solutions of this application will be described in detail, clearly, and completely below with reference to the accompanying drawings of the embodiments. It should be particularly noted that the specific embodiments described below are only used to better illustrate and explain the technical solutions of this application, and are intended to enable those skilled in the art to better understand and implement this application, and should not be construed as limiting the scope of protection of this application. Without departing from the spirit and substance of this application, those skilled in the art can modify, adjust, or make equivalent substitutions based on the content disclosed in this application, and these modifications, adjustments, or equivalent substitutions should all be considered within the scope of protection of this application.

[0022] Example 1: Please refer to Figure 1 As shown in the figure, this embodiment discloses an automatic analysis method for CT images of the leg bones of rats and mice, including: Step S1: Obtain the target long bone CT image sequence, extract the three-dimensional spatial distribution of bone tissue, and perform pose alignment and three-dimensional cropping based on the main direction of bone tissue to obtain a standardized bone image sequence.

[0023] It should be noted that the flowchart for obtaining standardized skeletal image sequences is as follows: Figure 2 As shown.

[0024] Furthermore, a more specific implementation of step S1 is as follows: The target long bone CT image sequence was acquired, and the three-dimensional spatial distribution of bone tissue was extracted. Specifically, the original CT scan results of the target long bone were first sequentially organized, and the slice number, inter-slice physical distance, and single-slice image resolution information of the consecutive slices were read. The slices were then sorted from smallest to largest according to their physical location to construct the original CT image sequence. Taking the femur of a mouse or rat as an example of the target long bone, the original CT image sequence can consist of 1800 consecutive slices, with a single slice resolution of 1024×1024 and a physical distance of 0.01mm between adjacent slices. This resulted in a three-dimensional grayscale volume data with a size of 1024×1024×1800, where each voxel location stores the corresponding CT grayscale value.

[0025] To improve the processing efficiency of subsequent extraction of the three-dimensional spatial distribution of bone tissue, resolution compression is performed on the original CT image sequence. Specifically, the horizontal and vertical dimensions of the original CT image sequence are resampled proportionally according to a preset compression ratio to generate resolution-compressed three-dimensional volume data. In one embodiment, the preset compression ratio can be set to 0.25, meaning that grayscale resampling is performed once for every 4×4 adjacent pixel regions, generating one compressed pixel. Bilinear interpolation can be used for the resampling method. After compression, the original CT image sequence is transformed from 1024×1024×1800 to 256×256×1800. Subsequently, threshold segmentation is performed on the resolution-compressed three-dimensional volume data to distinguish bone tissue regions from background regions, obtaining a three-dimensional binary mask of the bone tissue. The threshold segmentation threshold can be automatically determined using the maximum inter-class variance method. Taking one example, when traversing the grayscale threshold range of 0 to 255, the optimal segmentation threshold corresponding to the largest inter-class variance is calculated to be 148. Voxels with a grayscale value greater than or equal to 148 are marked as bone tissue voxels, and voxels with a grayscale value less than 148 are marked as background voxels.

[0026] The three-dimensional spatial distribution of bone tissue is extracted based on a three-dimensional binary mask. Specifically, the three-dimensional binary mask is traversed layer by layer and searched point by point to extract the voxel coordinates of all voxels marked as bone tissue, and written into the three-dimensional spatial coordinate matrix in the order of horizontal coordinate, vertical coordinate, and depth coordinate. Taking one embodiment as an example, when traversing to the position of depth coordinate 12, horizontal coordinate 86, vertical coordinate 143, if the corresponding mask value is 1, then the coordinates (86, 143, 12) are written into the three-dimensional spatial coordinate matrix; when traversing to the position of depth coordinate 12, horizontal coordinate 87, vertical coordinate 143, if the corresponding mask value is also 1, then the coordinates (87, 143, 12) are continued to be written. After all traversals, a three-column, 352481-row three-dimensional spatial coordinate matrix is ​​formed, which is used to characterize the overall spatial distribution range and spatial extension direction of the target long bone in the three-dimensional volume data.

[0027] After obtaining the three-dimensional spatial distribution of bone tissue, pose alignment is performed based on the principal direction of the bone tissue. Specifically, the geometric center of the three-dimensional spatial coordinate matrix is ​​first calculated to obtain the spatial center coordinates of the target long bone; then, using the spatial center coordinates as a reference, the coordinates of all bone tissue voxels are centered; subsequently, covariance analysis and eigenvalue decomposition are performed on the centered coordinate matrix to extract the eigenvector corresponding to the largest eigenvalue, which is determined as the principal direction vector of the bone tissue. The unit direction vector of the preset depth axis is set to (0, 0, 1). A dot product operation is performed on the principal direction vector of the bone tissue and the unit direction vector of the preset depth axis to obtain the cosine value of the angle; then, an inverse cosine operation is performed on the cosine value of the angle to obtain the skew angle of the principal direction of the bone tissue relative to the principal direction of the preset depth axis. Further, an outer product operation is performed on the principal direction vector of the bone tissue and the unit direction vector of the preset depth axis to obtain the rotation axis direction vector perpendicular to both.

[0028] After obtaining the rotation axis direction vector and the principal direction skew angle, the rotation axis direction vector is first normalized. Then, based on the normalized rotation axis direction vector, the principal direction skew angle, and the Rodriguez rotation rule, a 3x3 rotation matrix is ​​generated for attitude alignment to rotate the principal direction of the bone tissue to align with the preset depth axis. Once the rotation matrix is ​​determined, the following processing is performed on each voxel coordinate in the original CT image sequence: first, the spatial center coordinates are subtracted, and the voxel coordinates are translated to a local coordinate system with the spatial center coordinates as the origin; then, matrix multiplication is performed on the translated local coordinates using the rotation matrix to obtain the rotated local coordinates; finally, the rotated local coordinates are added back to the spatial center coordinates to obtain the transformed voxel coordinates. Since the transformed voxel coordinates usually contain a decimal part and no longer correspond to integer grid positions, trilinear interpolation resampling is performed on the rotated grayscale values ​​to generate the attitude-aligned 3D volume data.

[0029] Taking one embodiment as an example, the calculated spatial center coordinates are (508.4, 497.2, 266.8), and the extracted principal direction vector of the bone tissue is (0.118, 0.094, 0.989). A dot product operation is performed on the principal direction vector of the bone tissue and the preset depth axis unit direction vector (0, 0, 1), yielding a cosine value of approximately 0.9890. Then, an inverse cosine operation is performed on 0.9890, yielding a principal direction skew angle of approximately 8.5°. Finally, an outer product operation is performed on the principal direction vector of the bone tissue and the preset depth axis unit direction vector, yielding a rotation axis direction vector of approximately (0.623, -0.782, 0). Based on the normalized rotation axis direction vector, the principal direction skew angle of 8.5°, and the Rodriguez rotation rule, a rotation matrix is ​​constructed. The first row contains 0.9933, -0.0054, and -0.1156; the second row contains -0.0054, 0.9957, and -0.0921; and the third row contains 0.1156, 0.0921, and 0.9890. Taking the original voxel coordinates (620, 510, 300) as an example, after spatial center translation, rotation matrix transformation, and coordinate shift, the transformed voxel coordinates are approximately (615.34, 506.29, 313.72).

[0030] After pose alignment, 3D cropping is performed based on the pose-aligned 3D volume data to obtain a standardized skeletal image sequence. Specifically, bone tissue regions are re-identified in the pose-aligned 3D volume data, and the extreme coordinate values ​​of all bone tissue voxels in the lateral, longitudinal, and depth directions are extracted. The 3D cropping boundary surrounding the target long bone is determined based on the extreme values ​​in each direction, and the background region outside the cropping boundary is removed. For example, after alignment, the minimum and maximum coordinates of the bone tissue voxels in the lateral direction are 398 and 618, respectively; in the longitudinal direction, 421 and 657, respectively; and in the depth direction, 52 and 448, respectively. A 3D cropping region of size 220×236×396 can be constructed accordingly. The pose-aligned 3D volume data is then cropped according to this 3D cropping region to obtain a standardized skeletal image sequence that closely follows the edge of the target long bone and whose long axis direction is consistent with the preset depth axis.

[0031] Through the above processing, the standardized skeletal image sequence output in step S1 simultaneously possesses a unified spatial orientation and a compact analysis range. The standardized skeletal image sequence is three-dimensional volume data arranged along the horizontal, vertical, and depth coordinates. In subsequent steps, the depth coordinate is fixed to obtain cross-sectional slices, and the horizontal or vertical coordinate is fixed to obtain longitudinal slices, providing a consistent spatial reference for subsequent cross-sectional butterfly-shaped anatomical feature recognition, extraction of key points in longitudinal strip-shaped gaps, and anatomical zero-point locking. Taking one embodiment as an example, after rotation matrix pose alignment and three-dimensional cropping, the principal axis skew angle of the target long bone in different samples of the same batch can converge from the original state of 6.8° to 12.4° to within 0.5°. The subsequently extracted cross-sectional slices can more stably reflect the true cross-sectional topology near the growth plate, thereby reducing the impact of oblique cutting deformation on subsequent recognition results.

[0032] Methods for generating rotation matrices based on rotation axis direction vectors, principal direction skew angles, and Rodriguez rotation rules include: After obtaining the rotation axis direction vector and the principal direction skew angle, the rotation axis direction vector is first normalized. Let the three directional components of the rotation axis direction vector be the first, second, and third components of the rotation axis. First, the sum of squares and the square root of the three directional components are performed to obtain the rotation axis modulus. Then, the first, second, and third components of the rotation axis are divided by the rotation axis modulus to obtain the first, second, and third components of the unit rotation axis. If the rotation axis modulus is less than a preset minimum threshold, it is determined that the principal direction vector of the bone tissue is already in a state of overlap or approximately overlap with the preset depth axis. In this case, the rotation matrix is ​​directly set to the identity matrix, and no further rotation transformation is performed. For example, when the rotation axis direction vector is (0.623, -0.782, 0), the calculated rotation axis modulus is approximately 1.000. After normalization, the unit rotation axis component is still approximately (0.623, -0.782, 0).

[0033] After normalization, a 3x3 rotation matrix is ​​generated according to the Rodriguez rotation rule. Specifically, the cosine and sine values ​​of the skew angle of the principal directions are first calculated; then, based on the first, second, and third components of the unit rotation axis, the elements at each position of the rotation matrix are determined sequentially.

[0034] The elements in the first row and first column are equal to the product of the cosine of the skew angle plus the square of the first component of the unit rotation axis and "1 minus the cosine of the skew angle". The elements in the first row and second column are equal to the product of the first and second components of the unit rotation axis multiplied by "1 minus the cosine of the skew angle", minus the product of the third component of the unit rotation axis and the sine of the skew angle. The elements in the first row and third column are equal to the product of the first and third components of the unit rotation axis multiplied by "1 minus the cosine of the skew angle", plus the product of the second component of the unit rotation axis and the sine of the skew angle. The elements in the second row and first column are equal to the product of the second and first components of the unit rotation axis multiplied by "1 minus the cosine of the skew angle", plus the product of the third component of the unit rotation axis and the sine of the skew angle. The elements in the second row and second column are equal to the product of the cosine of the skew angle plus the square of the first component of the unit rotation axis and the square of the skew angle. The element in the second row and third column is equal to the product of the square of the second component of the rotation axis and "1 minus the cosine of the skew angle". The element in the third row and first column is equal to the product of the second component of the unit rotation axis and the third component of the unit rotation axis, multiplied by "1 minus the cosine of the skew angle", and minus the product of the first component of the unit rotation axis and the sine of the skew angle. The element in the third row and second column is equal to the product of the third component of the unit rotation axis and the first component of the unit rotation axis, multiplied by "1 minus the cosine of the skew angle", and minus the product of the second component of the unit rotation axis and the sine of the skew angle. The element in the third row and second column is equal to the product of the third component of the unit rotation axis and the second component of the unit rotation axis, multiplied by "1 minus the cosine of the skew angle", and added to the product of the first component of the unit rotation axis and the sine of the skew angle. The element in the third row and third column is equal to the cosine of the skew angle plus the product of the square of the third component of the unit rotation axis and "1 minus the cosine of the skew angle".

[0035] Taking one embodiment as an example, when the skew angle in the principal direction is 8.5°, the calculated cosine of the skew angle is approximately 0.9890, and the sine of the skew angle is approximately 0.1478. Substituting the unit rotation axis component (0.623, -0.782, 0), the cosine of the skew angle (0.9890), and the sine of the skew angle (0.1478) into the calculation rules of the matrix elements above, we can obtain the rotation matrix with the following values: first row approximately 0.9933, -0.0054, -0.1156; second row approximately -0.0054, 0.9957, -0.0921; and third row approximately 0.1156, 0.0921, 0.9890. The rotation matrix generated in this way can be directly used for subsequent voxel coordinate rotation transformations.

[0036] When the principal direction vector of the bone tissue is opposite to the preset depth axis, their dot product is close to -1. In this case, the cross product result is close to zero, making it impossible to directly determine a stable rotation axis from the cross product result. For this situation, an auxiliary direction vector that is not collinear with the principal direction vector of the bone tissue can be preset. The auxiliary direction vector and the principal direction vector of the bone tissue can then be cross-producted to obtain a substitute rotation axis. The skew angle of the principal direction is then set to 180°, and the rotation matrix is ​​constructed according to the aforementioned rotation matrix generation rules.

[0037] Furthermore, the settings for the preset depth axis, preset minimum threshold, and preset auxiliary direction vector are illustrated in the following example: The preset depth axis is used as the alignment reference for the main direction of the target long bone and is preset according to the coordinate definition method of the standardized skeletal image sequence. Specifically, in step S1, after sorting the original CT image sequence according to the physical position of the slices in ascending order and defining the layer number increment direction as the positive direction of the depth coordinate, the preset depth axis is set as a unit direction vector along the positive direction of the depth coordinate. To avoid inconsistencies in the definition of the depth direction between different samples, the slice sorting direction is unified before entering the pose alignment calculation, and then the same unit direction vector is fixed as the preset depth axis. Taking one embodiment as an example, when the standardized skeletal image sequence is arranged according to the horizontal coordinate, vertical coordinate, and depth coordinate, and the depth layer number increases from layer 0 to layer 395, the preset depth axis can be set as the unit direction vector (0, 0, 1). If the slice storage order of the original scan result is opposite to the predetermined depth direction, the slice order is first rearranged to the layer number increment direction, and then (0, 0, 1) is continued to be used as the preset depth axis, thereby ensuring that the subsequent dot product, outer product, and rotation matrix construction uses a unified direction reference.

[0038] The preset minimum threshold is used to determine whether the magnitude of the rotation axis direction vector is close to 0, in order to distinguish between the case where "the principal direction of bone tissue coincides with or is approximately coincident with the preset depth axis" and the case where "normal rotation still needs to be performed". Specifically, the preset minimum threshold is set based on the coordinate normalization accuracy, the floating-point operation error range, and the voxel coordinate perturbation level, so that the minimum magnitude caused by numerical error will not be misjudged as a valid rotation axis. To ensure that the threshold setting is clear and reproducible, the principal direction vector of bone tissue and the unit direction vector of the preset depth axis can be normalized first, and then the preset minimum threshold can be determined according to the number of decimal places retained by the normalized vector components. For example, in one embodiment, when each component of the normalized direction vector retains 4 decimal places, the preset minimum threshold can be set to 0.001. If the calculated rotation axis direction vector is (0.0002, -0.0003, 0.0001), its magnitude is approximately 0.0004, which is less than 0.001. At this point, it is determined that the principal direction of the bone tissue approximately coincides with the preset depth axis, and the rotation matrix is ​​set as the identity matrix. In another embodiment, when the rotation axis direction vector is (0.0200, -0.0150, 0.0000), its magnitude is approximately 0.0250, which is greater than 0.001. Therefore, the normalization process and rotation matrix construction continue. Furthermore, when the magnitude of the rotation axis direction vector is less than a preset minimum threshold and the dot product of the bone tissue principal direction vector and the preset depth axis unit direction vector is greater than 0, it is determined to be in the same direction or approximately in the same direction; when the magnitude of the rotation axis direction vector is less than a preset minimum threshold and the dot product is less than 0, it is determined to be in the opposite direction or approximately in the opposite direction, and enters the alternative rotation axis determination process involving auxiliary direction vectors.

[0039] The preset auxiliary direction vector is used to handle situations where the main direction vector of bone tissue coincides or is approximately coincident with the unit direction vector of the preset depth axis in opposite directions, and the cross product result is close to zero, making it impossible to directly obtain a stable rotation axis. Specifically, the preset auxiliary direction vector is selected from three pre-defined coordinate axis basis vectors according to fixed rules to ensure that the selected auxiliary direction vector is not collinear with the main direction vector of bone tissue. Preferably, the first auxiliary candidate vector is preset as (1, 0, 0), the second auxiliary candidate vector is (0, 1, 0), and the third auxiliary candidate vector is (0, 0, 1); the absolute value of the dot product between the main direction vector of bone tissue and the three auxiliary candidate vectors is calculated respectively, and the auxiliary candidate vector with the smallest absolute value of the dot product is selected as the preset auxiliary direction vector; when two or more auxiliary candidate vectors have the same absolute value of the dot product, the preset auxiliary direction vector is determined according to the priority order of the first auxiliary candidate vector, the second auxiliary candidate vector, and the third auxiliary candidate vector. Taking one embodiment as an example, when the principal direction vector of bone tissue is (0.0000, 0.0000, -1.0000), the absolute values ​​of its dot product with (1, 0, 0), (0, 1, 0), and (0, 0, 1) are 0, 0, and 1, respectively. At this time, (1, 0, 0) is selected as the preset auxiliary direction vector according to the priority order. Then, the outer product operation is performed on the preset auxiliary direction vector (1, 0, 0) and the principal direction vector of bone tissue (0.0000, 0.0000, -1.0000) to obtain the replacement rotation axis direction vector (0, 1, 0), and the skew angle of the principal direction is set to 180°. Based on this, the rotation matrix is ​​constructed. In another example, when the principal direction vector of bone tissue is (0.1200, -0.9800, 0.1500), the absolute values ​​of its dot product with (1, 0, 0), (0, 1, 0), and (0, 0, 1) are 0.1200, 0.9800, and 0.1500, respectively. Therefore, (1, 0, 0) is selected as the preset auxiliary direction vector.

[0040] Step S2: Construct a first candidate depth interval based on the butterfly-shaped anatomical features in the cross-sectional slice, and construct a second candidate depth interval based on the key point depth coordinates of the growth plate-like gaps in the longitudinal slice. Perform an intersection operation on the first candidate depth interval and the second candidate depth interval to obtain the fine search range.

[0041] Furthermore, a more specific implementation of step S2 is as follows: Obtain the standardized skeletal image sequence output in step S1. The standardized skeletal image sequence is three-dimensional volumetric data arranged along horizontal, vertical, and depth coordinates. The depth coordinates correspond to the layer number order of the standardized skeletal image sequence, and the depth coordinate direction is consistent with a preset depth axis. Based on the standardized skeletal image sequence, slices perpendicular to the preset depth axis are defined as cross-sectional slices. These cross-sectional slices are obtained by fixing the depth coordinates and extracting the corresponding two-dimensional image. Slices containing the preset depth axis are defined as longitudinal slices. These longitudinal slices are obtained by fixing either the horizontal or vertical coordinates and reconstructing along the remaining coordinate directions and the depth coordinate direction. The longitudinal slices reconstructed by fixing the horizontal coordinates are coronal slices, and the longitudinal slices reconstructed by fixing the vertical coordinates are sagittal slices. For example, if the size of the standardized bone image sequence is 220×236×396, then fixing the depth coordinate 150 can obtain a cross-sectional slice with a size of 220×236, fixing the transverse coordinate 110 can obtain a coronal slice with a size of 236×396, and fixing the longitudinal coordinate 118 can obtain a sagittal slice with a size of 220×396.

[0042] The first candidate depth interval is constructed based on the butterfly-shaped anatomical features in cross-sectional slices. Specifically, a preset retrieval step size is set for the cross-sectional slices. Starting from the initial depth layer of the standardized skeletal image sequence, cross-sectional slices are extracted sequentially according to the preset retrieval step size, and butterfly-shaped anatomical feature recognition is performed on each extracted cross-sectional slice. For cross-sectional slices with a confidence level greater than or equal to a preset judgment threshold in the recognition results, the corresponding depth layer number is recorded to form a set of butterfly-shaped feature hit layer numbers. All layer numbers in the butterfly-shaped feature hit layer number set are compared to determine the minimum and maximum hit layer numbers. Then, based on a preset boundary compensation layer number, the minimum and maximum hit layer numbers are expanded forward and backward, respectively, to obtain the starting and ending layer numbers of the first candidate depth interval. The preset boundary compensation layer is used to compensate for boundary omissions caused by step size extraction. It can be set to the number of layers corresponding to the preset retrieval step size, or a part of the number of layers corresponding to the preset retrieval step size. When the expanded starting layer number is less than the starting layer number of the standardized skeletal image sequence, the starting layer number of the standardized skeletal image sequence is used as the starting layer number of the first candidate depth interval. When the expanded ending layer number is greater than the ending layer number of the standardized skeletal image sequence, the ending layer number of the standardized skeletal image sequence is used as the ending layer number of the first candidate depth interval.

[0043] It should be noted that the epiphyseal tissue near the joint end of long bones exhibits a symmetrical, bilaterally expanded morphology in cross-section. This morphology appears as a stable butterfly-shaped anatomical feature in standardized skeletal imaging sequences. Since this feature is continuously distributed in the depth direction and located near the metaphyseal region containing the growth plate, it can be used to narrow the search area from the entire long bone to a local depth range near the growth plate, thus providing coarse localization results in the cross-sectional direction for subsequent fine localization.

[0044] Taking one embodiment as an example, the standardized skeletal image sequence has a total depth of 396 layers, a depth interlayer spacing of 0.01 mm, a preset retrieval step size of 8 layers, a preset boundary compensation layer size of 8 layers, and a confidence threshold of 0.55 for butterfly-shaped anatomical features. After extracting cross-sectional slices according to the 8-layer step size, butterfly-shaped anatomical features are identified in the cross-sectional slices corresponding to depth layer numbers 120, 128, 136, 144, 152, 160, 168, 176, 184, 192, and 200, forming a butterfly-shaped feature hit layer number set {120, 128, 136, 144, 152, 160, 168, 176, 184, 192, 200}. Comparison of this set yields a minimum hit layer number of 120 and a maximum hit layer number of 200; then, compensation is applied 8 layers forward and backward respectively, resulting in the first candidate depth range of layers 112 to 208.

[0045] In one embodiment, butterfly-shaped anatomical feature recognition is implemented using an object detection network, which can be implemented using the YOLOv8-based detect task. Specifically, cross-sectional slices are first extracted from a standardized skeletal image sequence according to a preset retrieval step size. Size unification and grayscale normalization are then performed on each cross-sectional slice to generate the input image for the object detection network. Size unification transforms cross-sectional slices of different sizes into a fixed input size, while grayscale normalization maps the original grayscale values ​​to a preset numerical range. For example, in one embodiment, each cross-sectional slice is uniformly adjusted to 640×640 pixels, and the grayscale values ​​are linearly mapped to the 0-1 range.

[0046] The object detection network takes a single cross-sectional slice that has been sized and grayscale normalized as input, and outputs a set of candidate detection boxes within the current cross-sectional slice. Each candidate detection box includes at least the horizontal coordinates of its center, the vertical coordinates of its center, its width, its height, and a confidence score for its butterfly-shaped anatomical features. For multiple candidate detection boxes output from the same cross-sectional slice, they are first sorted according to their confidence scores, and then non-maximum suppression is used to remove redundant candidate detection boxes whose overlap exceeds a preset overlap threshold. The remaining candidate detection boxes are retained as valid detection results for the current cross-sectional slice. In one embodiment, the preset overlap threshold can be set to 0.45.

[0047] After obtaining the valid detection results for the current cross-sectional slice, layer hit determination is performed. Specifically, the highest confidence score among all valid detection results for the current cross-sectional slice is read. When the highest confidence score is greater than or equal to a preset determination threshold, the depth layer number corresponding to the current cross-sectional slice is recorded in the butterfly feature hit layer number set. When the highest confidence score is less than the preset determination threshold, or when there are no valid detection results for the current cross-sectional slice, the depth layer number corresponding to the current cross-sectional slice is not recorded in the butterfly feature hit layer number set. Taking one embodiment as an example, the preset determination threshold is set to 0.55. If the cross-sectional slice corresponding to depth layer number 144 outputs two candidate detection boxes after detection, and one valid detection box is retained after non-maximum suppression processing with a confidence score of 0.81, then depth layer number 144 is recorded in the butterfly feature hit layer number set. If the highest confidence score of the cross-sectional slice corresponding to depth layer number 152 is 0.42, then depth layer number 152 is not recorded in the butterfly feature hit layer number set.

[0048] To ensure the object detection network can stably identify butterfly-shaped anatomical features in cross-sectional slices, a corresponding training sample set is first constructed. Specifically, cross-sectional slices are extracted from multiple standardized skeletal image sequences that have undergone pose alignment and 3D cropping. Annotators draw rectangular bounding boxes around slices containing butterfly-shaped anatomical features, delineating the circumscribed region of the butterfly-shaped anatomical features within the cross-sectional slice; slices without butterfly-shaped anatomical features are marked as empty target slices. For example, in one embodiment, 4800 cross-sectional slices can be extracted from 120 sets of standardized skeletal image sequences, including 2600 positive sample slices and 2200 negative sample slices. After rectangular bounding box annotation of the butterfly-shaped anatomical features in the positive sample slices, the data is divided into training, validation, and test sets in an 8:1:1 ratio for object detection network training, parameter validation, and recognition performance evaluation.

[0049] In one embodiment, the object detection network employs a single-stage object detection structure. During training, labeled cross-sectional slices are input into the object detection network. A loss function is constructed using the positional deviation between the labeled bounding box and the network's predicted bounding box, as well as the object category prediction deviation. The network parameters are updated through iterative training. To enhance the network's adaptability to different bone morphologies, imaging brightness, and local noise conditions, random flipping, random scaling, random cropping, and brightness perturbation processing can be performed on the input image during the training phase. For example, in one embodiment, the training epochs are set to 300, the batch size to 16, the initial learning rate to 0.001, and the weight decay coefficient to 0.0005.

[0050] In the application phase, for each cross-sectional slice extracted according to a preset retrieval step size, the processing flow of input image generation, candidate detection box output, non-maximum suppression, and layer hit determination is repeatedly executed to form a butterfly-shaped feature hit layer number set. Then, all layer numbers in the butterfly-shaped feature hit layer number set are compared to determine the minimum and maximum hit layer numbers, and combined with the preset boundary compensation layer number to obtain the first candidate depth interval. Taking one embodiment as an example, if the butterfly-shaped feature hit layer number set is {120, 128, 136, 144, 152, 160, 168, 176, 184, 192, 200}, and the preset boundary compensation layer number is 8 layers, then the minimum hit layer number 120 and the maximum hit layer number 200 are taken, and 8 layers are compensated forward and backward respectively to obtain the first candidate depth interval of layers 112 to 208.

[0051] A second candidate depth interval is constructed based on the depth coordinates of key points in the growth plate-like gaps within longitudinal slices. Specifically, longitudinal slices are reconstructed in a standardized skeletal image sequence, and key point identification is performed on the growth plate-like gaps in each longitudinal slice. For key points with a confidence level greater than or equal to a preset key point threshold in the key point identification results, the coordinate values ​​of the key points in the depth direction of the longitudinal slice are extracted, and these coordinate values ​​are mapped to depth layer numbers in the standardized skeletal image sequence to form a set of key point depth layer numbers. For the depth layer numbers corresponding to all valid key points, a comparison is performed to determine the minimum and maximum key point layer numbers; then, based on a preset key point boundary compensation layer number, the minimum and maximum key point layer numbers are expanded forward and backward, respectively, to obtain the starting and ending layer numbers of the second candidate depth interval. The preset key point boundary compensation layers are used to compensate for key point recognition errors and boundary fluctuations caused by differences in the vertical slice reconstruction position. When the expanded starting layer number is less than the starting layer number of the standardized skeletal image sequence, the starting layer number of the standardized skeletal image sequence is used as the starting layer number of the second candidate depth interval. When the expanded ending layer number is greater than the ending layer number of the standardized skeletal image sequence, the ending layer number of the standardized skeletal image sequence is used as the ending layer number of the second candidate depth interval.

[0052] Taking one embodiment as an example, five coronal slices are reconstructed near the lateral center of a standardized skeletal image sequence. The preset number of keypoints is set to 10, the preset keypoint threshold is set to 0.50, and the preset number of keypoint boundary compensation layers is set to 2. A total of 13 valid keypoints are identified in the five coronal slices. The corresponding depth coordinates are rounded to obtain the keypoint depth layer number set {145, 148, 150, 152, 156, 160, 164, 169, 173, 176, 180, 184, 186}. Comparison of this set yields the minimum keypoint layer number as 145 and the maximum keypoint layer number as 186. Then, two layers are compensated forward and backward respectively to obtain the second candidate depth range from layer 143 to layer 188.

[0053] After obtaining the first and second candidate depth intervals, an intersection operation is performed on the two candidate depth intervals to obtain a refined search range. Specifically, the first and second candidate depth intervals are first represented as sets of continuous depth layer numbers containing a start layer number and an end layer number, respectively. The start and end layer numbers of the first and second candidate depth intervals are read; the larger of the two start layer numbers is determined as the undetermined start layer number, and the smaller of the two end layer numbers is determined as the undetermined end layer number. When the undetermined start layer number is less than or equal to the undetermined end layer number, the set of continuous depth layer numbers including the undetermined start and end layer numbers is determined as the refined search range; when the undetermined start layer number is greater than the undetermined end layer number, it is determined that the first and second candidate depth intervals do not have a direct intersection, and the first and second candidate depth intervals are expanded outward by a preset conflict compensation layer number, and the intersection operation is performed again until a refined search range is obtained that satisfies the condition that the undetermined start layer number is less than or equal to the undetermined end layer number, or the preset recalculation limit is reached. The preset conflict compensation layers are used to compensate for errors in cross-sectional slice step size extraction, errors in longitudinal slice key point identification, and interval misalignment caused by insufficient boundary compensation; the preset maximum number of recalculations can be set to 1 to 3 times. After the fine search range is determined, all continuous depth layer numbers in the fine search range are used as the set of search layer numbers for subsequent triangular fracture candidate slices; if it is necessary to output the physical depth range, the starting layer number and ending layer number of the fine search range are multiplied by the depth layer spacing to obtain the starting physical depth and ending physical depth corresponding to the fine search range.

[0054] Taking one embodiment as an example, if the first candidate depth range is from layer 112 to layer 208, and the second candidate depth range is from layer 143 to layer 188, then the larger value of the starting layer number, 143, is determined as the undetermined starting layer number, and the smaller value of the ending layer number, 188, is determined as the undetermined ending layer number. Since 143 is less than 188, the set of consecutive depth layer numbers corresponding to layers 143 to 188 is determined as the fine search range. If the depth layer spacing is 0.01 mm, then the physical depth range corresponding to the fine search range is 1.43 mm to 1.88 mm. Taking another example, if the first candidate depth range is from 112 to 156 layers and the second candidate depth range is from 160 to 188 layers, then the intersection is empty. When the preset conflict compensation layer is set to 4 layers, the first candidate depth range is expanded to 108 to 160 layers and the second candidate depth range is expanded to 156 to 192 layers. After performing the intersection operation again, the fine search range is obtained as 156 to 160 layers.

[0055] Through the above processing, the fine search range output in step S2 is simultaneously constrained by both the cross-sectional butterfly-shaped anatomical features and the depth coordinates of key points in the longitudinal strip-shaped gaps. Specifically, the first candidate depth interval provides a coarse localization result based on the cross-sectional morphology, the second candidate depth interval provides a verification result based on the spatial distribution of the longitudinal strip-shaped gaps, and the fine search range obtained through the intersection operation is used to limit the retrieval depth range of subsequent triangular fracture candidate slices, thereby reducing boundary drift and false positive interference caused by a single identification result.

[0056] Step S3: Identify candidate slice sequences of triangular fractures where epiphyseal and metaphyseal bone tissues transition from a contact state to a disconnected state within a fine search range. Calculate the geometric gap distribution sequence of the epiphyseal and metaphyseal edges in each candidate slice of triangular fractures. Determine the slice layer number with the smallest variance of the geometric gap distribution sequence as the anatomical zero point.

[0057] The flowchart for obtaining the candidate slice sequence of the triangular fracture is as follows: Figure 3 As shown, the flowchart for obtaining the geometric gap distribution sequence is as follows: Figure 4 As shown.

[0058] It should be noted that as the cross-sectional slices progress along the depth axis from small to large, the relative relationship between the epiphyseal and metaphyseal bone tissues near the growth plate gradually transitions from mutual contact to mutual separation. Near the initial boundary of the growth plate, the inner edges of both the upper and lower epiphyseal and metaphyseal bone tissues converge towards the center, forming corresponding upper and lower apical contours. When the two apical edges change from contact to separation, local low-density gaps appear on the cross-section, forming a triangular fracture morphology. This triangular fracture morphology corresponds to the critical layer where the epiphyseal and metaphyseal bone tissues transition from a contact state to a disconnected state. Based on a critical feature detection neural network model, state determination and candidate region screening are performed on the cross-sectional slices within a fine search range, extracting a sequence of candidate slices for triangular fractures transitioning from a contact state to a disconnected state. Then, based on the geometric gap distribution sequence of the epiphyseal and metaphyseal edges corresponding to each candidate slice in the triangular fracture candidate slice sequence, the variance of each candidate slice is compared, and the candidate slice layer number with the most uniform geometric gap distribution is determined as the anatomical zero point. By first screening candidate slices of triangular fracture and then performing a comparison of the stability of geometric gap distribution, the determination of the anatomical zero point can be limited to the critical level where the bone tissue connection relationship changes abruptly, reducing the interference of morphological differences in non-critical levels on the zero point localization results, thereby providing a unified spatial reference benchmark for subsequent target analysis region mapping.

[0059] Furthermore, a more specific implementation of step S3 is as follows: Obtain the fine search range determined in step S2, and extract all cross-sectional slices within the fine search range. Specifically, each layer number in the set of consecutive depth layer numbers corresponding to the fine search range corresponds to one cross-sectional slice. Sort all cross-sectional slices in ascending order of depth layer number to form a sequence of slices to be judged. Taking one embodiment as an example, if the fine search range is from layer 143 to layer 188, then the sequence of slices to be judged includes 46 cross-sectional slices corresponding to depth layer numbers 143, 144, 145 up to 188.

[0060] In one implementation, a critical feature detection neural network model is used to identify candidate triangular fracture slices transitioning from a contact state to a disconnected state. This model performs local target localization and state determination on cross-sectional slices within a fine search range, and outputs the local region localization results required for subsequent geometric analysis. Specifically, each cross-sectional slice in the slice sequence to be determined is input into the critical feature detection neural network model, which outputs a set of candidate bounding boxes, a state category index, and a state confidence score for the current cross-sectional slice. The state category index includes at least contact and disconnected states. The candidate bounding boxes delineate the local analysis region where contact or disconnection occurs between epiphyseal and metaphyseal bone tissue, and the state category index characterizes whether the current cross-sectional slice is in a contact or disconnected state within that local analysis region. To improve the distinguishability of bone tissue edges and local low-density gaps, grayscale stretching, size unification, and grayscale normalization can be performed on the cross-sectional slices before inputting them into the critical feature detection neural network model. For example, in one embodiment, the cross-sectional slices can be uniformly adjusted to 640×640 pixels, and the grayscale values ​​can be linearly mapped to the 0-1 range.

[0061] In one optional implementation, the critical feature detection neural network model is implemented using the YOLOv8-based detect task. Specifically, the labeled cross-sectional slice is input into the network corresponding to the YOLOv8-based detect task, and a set of candidate bounding boxes is output. Each candidate bounding box includes at least the horizontal coordinates of the box center, the vertical coordinates of the box center, the box width, the box height, the state category index, and the state confidence score. For multiple candidate bounding boxes output from the same cross-sectional slice, they are first sorted according to their state confidence scores, and then non-maximum suppression is used to delete redundant candidate bounding boxes whose overlap exceeds a preset overlap threshold, retaining the remaining candidate bounding boxes as the valid detection results of the current cross-sectional slice. For example, in one embodiment, the preset overlap threshold can be set to 0.45.

[0062] Specifically, a training sample set for the critical feature detection neural network model is first constructed. The sample source is cross-sectional slices located near the growth plate in the standardized bone image sequence output in step S1. Taking one embodiment as an example, 5400 cross-sectional slices can be extracted from 180 sets of standardized bone image sequences, including 2200 slices in contact, 1900 slices in disconnection, and 1300 slices that do not meet the criteria. The slices in contact correspond to slices where the inner edges of the epiphyseal bone tissue and the inner edges of the metaphyseal bone tissue are continuously connected and the minimum edge spacing is no greater than 0.5 pixels; the slices in disconnection correspond to slices where there is a local low-density gap between the inner edges of the epiphyseal bone tissue and the inner edges of the metaphyseal bone tissue, the minimum edge spacing is greater than 0.5 pixels and no greater than 6.0 pixels, and the upper and lower edges form a triangular apex outline near the center gap position; the slices that do not meet the above two conditions are slices.

[0063] During the annotation process, rectangular bounding boxes are drawn for both contact and disconnected slices. These boxes define the local analysis area where the epiphyseal and metaphyseal bone tissues are in contact or disconnected, and are assigned the category labels "Contact State" and "Disconnected State," respectively. Slices that do not match the annotation are not drawn and are used as empty annotation samples for training. To avoid the same sample data appearing in both the training and validation sets simultaneously, the dataset is divided according to individual samples. For example, in one embodiment, the training, validation, and test sets can be divided in an 8:1:1 ratio, corresponding to 4320 images in the training set, 540 images in the validation set, and 540 images in the test set.

[0064] For input construction, each cross-sectional slice is first subjected to grayscale stretching, then uniformly adjusted to 640×640 pixels, and the grayscale values ​​are linearly mapped to the 0 to 1 range. For single-channel grayscale images, a 3-channel duplication method can be used to generate network input, that is, the same grayscale image is copied to 3 input channels. In order to maintain the anatomical consistency of the contact state and the disconnected state in the vertical direction, vertical flip enhancement is not enabled during the training phase; random horizontal flip, random translation, random scaling, brightness perturbation, and mild noise perturbation can be enabled. Taking one embodiment as an example, the probability of random horizontal flip is set to 0.5, the random translation ratio is set to -0.10 to 0.10, the random scaling ratio is set to 0.85 to 1.15, the brightness perturbation ratio is set to -0.15 to 0.15, and the standard deviation of Gaussian noise is set to 0.01.

[0065] Regarding network structure and training parameters, the YOLOv8-s detect task model is used as an example of a critical feature detection neural network model. The training epochs are set to 300, the batch size to 16, and the AdamW optimizer is used. The initial learning rate is set to 0.001, the minimum learning rate to 0.00001, the weight decay coefficient to 0.0005, and the warm-up epochs to 3. During the warm-up phase, the learning rate linearly increases from 0.0001 to 0.001. After the warm-up, a cosine annealing strategy is used to gradually decrease the learning rate. The first-order momentum parameter of the AdamW optimizer is set to 0.9, and the second-order momentum parameter is set to 0.999. The loss function includes bounding box regression loss, class classification loss, and distribution focus loss. The weight of the bounding box regression loss is set to 7.5, the weight of the class classification loss is set to 0.5, and the weight of the distribution focus loss is set to 1.5.

[0066] During training, in each iteration, cross-sectional slices from the training set are input into the critical feature detection neural network model to obtain candidate bounding boxes, category prediction results, and localization results. Then, the loss function is calculated based on the bounding boxes and category labels, and backpropagation is used to update the model parameters. In the validation phase, after every 5 training iterations, the validation set is used to evaluate model performance, and the average precision of the contact state class, the average precision of the disconnected state class, and the recall rate of the disconnected state class are recorded. To improve the stability of the disconnected state slice selection in step S3, the disconnected state class recall rate can be used as a priority metric. For example, in one embodiment, when the disconnected state class recall rate reaches 0.95 and the overall average precision reaches 0.90, the model parameters for the current iteration are saved as a candidate model; when no better validation results appear after 30 consecutive training iterations, training is terminated, and the model with the best validation metrics is selected as the final model.

[0067] During the application phase, each cross-sectional slice within the fine search range is input into the trained critical feature detection neural network model, which outputs a set of candidate bounding boxes, a state category index, and a state confidence score. Non-maximum suppression is performed on the candidate bounding box set, with a preset overlap threshold set to 0.45. The highest state confidence score among the retained candidate bounding boxes is then read. If the highest state confidence score is greater than or equal to 0.60 and the corresponding category is the contact state, the current cross-sectional slice is marked as a contact state slice. If the highest state confidence score is greater than or equal to 0.60 and the corresponding category is the disconnected state, the current cross-sectional slice is marked as a disconnected state slice. If there are no valid candidate bounding boxes, or the highest state confidence score is less than 0.60, the current cross-sectional slice is marked as a missed slice.

[0068] After determining the state of all cross-sectional slices, the sequence of slices to be determined is traversed in ascending order of depth layer number to search for the interlayer position where the state transitions from contact to disconnection. Specifically, when the cross-sectional slice with depth layer number i is marked as a contact slice and the cross-sectional slice with depth layer number i+1 is marked as a disconnection slice, the cross-sectional slice corresponding to depth layer number i+1 is determined as the first candidate slice for triangular fracture. The traversal continues, and cross-sectional slices that are continuously adjacent to the first candidate slice for triangular fracture and are continuously marked as disconnection slices are added to the candidate slice for triangular fracture. When subsequent cross-sectional slices are no longer marked as disconnection slices, the expansion of the current candidate slice for triangular fracture is terminated. Taking one embodiment as an example, within the fine search range of layers 143 to 188, the cross-sectional slices corresponding to depth layers 143 to 149 are all marked as contact state slices by the critical feature detection neural network model, the cross-sectional slices corresponding to depth layers 150, 151 and 152 are all marked as disconnected state slices by the critical feature detection neural network model, and the cross-sectional slice corresponding to depth layer 153 is marked as a missed slice. Thus, the triangular break candidate slice sequence is obtained as 3 cross-sectional slices corresponding to layers 150 to 152.

[0069] For each candidate slice in the triangular fracture candidate slice sequence, a local analysis region is determined based on the effective target box output by the critical feature detection neural network model, and the epiphyseal edge and metaphyseal edge are extracted within the local analysis region. Specifically, the region enclosed by the effective target box is used as the local analysis region, and bone tissue mask extraction is performed on the local analysis region to obtain a binary image of bone tissue; then, the region is scanned column by column in the horizontal direction, and the lowest boundary point of bone tissue above the local low-density gap is recorded in each column to form the epiphyseal edge curve; the highest boundary point of bone tissue below the local low-density gap is recorded to form the metaphyseal edge curve. To improve the edge localization accuracy, interpolation smoothing processing can be performed on the epiphyseal edge curve and the metaphyseal edge curve; when the sampling position falls at a non-integer pixel coordinate, linear interpolation can be performed based on the coordinates of adjacent edge points to obtain the corresponding longitudinal coordinate value of the edge.

[0070] It should be noted that, Figure 5 This is a schematic diagram of butterfly-shaped anatomical feature recognition and local analysis region localization in an embodiment of this application. The outer bounding box is used to represent the candidate detection region of butterfly-shaped anatomical features in the cross-sectional slice, and the inner bounding box is used to represent the local analysis region further determined based on the candidate detection region.

[0071] For each candidate slice, a geometric gap distribution sequence is calculated. Specifically, the central gap position is first determined in the current candidate slice. The central gap position is preferably determined by the lateral coordinates of the center of the effective target box output by the critical feature detection neural network model, or it can be determined by the average of the lateral coordinates of the epiphyseal edge tip and the metaphyseal edge tip. Using the central gap position as the center, a preset sampling width is set on both sides in the lateral direction, and multiple sampling positions are evenly selected within the sampling width according to a preset number of sampling points. For each sampling position, the longitudinal coordinates of the epiphyseal edge and the metaphyseal edge at the corresponding position are read, and the difference between the two is determined as the edge spacing of the current sampling position. All edge spacings are arranged sequentially according to the sampling position order to form the geometric gap distribution sequence of the current candidate slice. If it is necessary to unify the physical scale between different samples, the edge spacing can be multiplied by the pixel spacing corresponding to the cross-sectional slice to obtain a geometric gap distribution sequence in the form of physical distance.

[0072] Taking one embodiment as an example, in the candidate slice corresponding to depth layer number 150, the effective target box output by the critical feature detection neural network model has a horizontal center coordinate of 110.0, a vertical center coordinate of 152.0, a box width of 24.0, a box height of 18.0, a state category index of "open state", and a state confidence score of 0.87. Using the horizontal center coordinate of the effective target box, 110.0, as the center gap position, 5 pixels are taken on each side in the horizontal direction as the sampling width, and 5 sampling positions are evenly selected. The horizontal coordinates of the 5 sampling positions are 104.5, 107.0, 109.5, 112.0, and 114.5, respectively. At the five sampling locations, the longitudinal coordinates of the epiphyseal edge are 150.1, 150.4, 150.3, 150.5, and 150.2, respectively, while the longitudinal coordinates of the metaphyseal edge are 153.5, 153.7, 153.5, 153.8, and 153.7, respectively. Subtracting these two values ​​yields a geometric gap distribution sequence of 3.4, 3.3, 3.2, 3.3, and 3.5 pixels. Taking the candidate slice corresponding to depth layer 151 as an example, the geometric gap distribution sequence obtained using the same sampling rule is 3.2, 3.2, 3.2, 3.3, and 3.2 pixels; and taking the candidate slice corresponding to depth layer 152 as an example, the geometric gap distribution sequence obtained using the same sampling rule is 3.1, 3.5, 3.2, 3.8, and 3.0 pixels.

[0073] The anatomical zero point is determined based on the variance of the geometric gap distribution sequence. Specifically, the mean of the geometric gap distribution sequence for each candidate slice is calculated. Then, the sum of the squares of the differences between the edge spacing at each sampling location and the mean is calculated and divided by the number of sampling points to obtain the variance of the geometric gap distribution sequence for the current candidate slice. The smaller the variance value, the more uniform the spacing between the epiphyseal and metaphyseal edges in the current candidate slice, indicating that the current layer is closer to the regular critical layer where the growth plate transitions from contact to separation. The variance values ​​corresponding to all candidate slices are compared, and the candidate slice layer number with the smallest variance value is taken as the anatomical zero point.

[0074] Taking one embodiment as an example, variance calculation was performed on the geometric gap distribution sequences 3.4, 3.3, 3.2, 3.3, and 3.5 corresponding to depth slice number 150, yielding a variance value of 0.01; variance calculation was performed on the geometric gap distribution sequences 3.2, 3.2, 3.2, 3.3, and 3.2 corresponding to depth slice number 151, yielding a variance value of 0.002; and variance calculation was performed on the geometric gap distribution sequences 3.1, 3.5, 3.2, 3.8, and 3.0 corresponding to depth slice number 152, yielding a variance value of 0.10. Comparing the variance values ​​of the three candidate slices, the minimum value of 0.002 corresponds to depth slice number 151; therefore, depth slice number 151 was determined as the anatomical zero point.

[0075] Through the above processing, the anatomical zero-point layer number output in step S3 directly corresponds to the cross-sectional slice position with the most uniform geometric gap distribution among the triangular fracture candidate slices selected by the critical feature detection neural network model within the fine search range. This anatomical zero point provides a unified spatial reference benchmark for subsequent target analysis region mapping, thereby reducing quantitative analysis errors caused by growth plate positioning offsets between different samples.

[0076] It should be noted that identifying candidate triangular fracture slices transitioning from a contact state to a disconnected state within the fine search range serves to further converge the depth interval-level localization results obtained in step S2 to the set of critical layers directly corresponding to the growth plate's initial boundary. Specifically, while the intersection of the first and second candidate depth intervals can limit the search range to the region near the growth plate, the fine search range still simultaneously includes bone tissue layers in a contact state, bone tissue layers in a disconnected state, and layers that do not meet the critical morphological conditions. If anatomical zero-point determination is performed directly within the entire fine search range, it is easily affected by morphological differences in non-critical layers, leading to inter-layer drift in the zero-point localization results. Therefore, this invention invokes a critical feature detection neural network model within the fine search range to perform local target localization and state determination on each cross-sectional slice, and filters out candidate triangular fracture slice sequences transitioning from a contact state to a disconnected state, thereby limiting subsequent geometric analysis to the critical layers where the bone tissue connectivity changes abruptly.

[0077] Furthermore, the triangular fracture candidate slices correspond to the transitional stage where the epiphyseal and metaphyseal bone tissues change from continuous connection to gap separation. This transitional stage simultaneously features localized low-density gaps and corresponding triangular apex contours, distinguishing it from both the continuous bone tissue morphology in the contact state and the wide gap morphology in the clearly disconnected state. Therefore, it possesses strong structural exclusivity and anatomical indicativeness. Based on the triangular fracture candidate slice sequence screened by the critical feature detection neural network model, the geometric gap distribution sequence of the epiphyseal and metaphyseal edges in each candidate slice is calculated, and the variance of the geometric gap distribution sequence corresponding to each candidate slice is compared. This allows for further identification of the layer with the most uniform geometric gap distribution from the critical transitional layer, and the corresponding layer number is determined as the anatomical zero point that best matches the growth plate's initial boundary. Thus, this invention improves the consistency and anti-interference capability of anatomical zero point localization through a step-by-step convergence method of "fine search range limitation—critical feature detection neural network model screening of triangular fracture candidate slices—variance comparison of geometric gap distribution sequence," providing a stable spatial reference benchmark for subsequent target analysis region mapping.

[0078] Step S4: Based on the anatomical zero point, offset along the depth axis by a preset number of layers and extract the preset number of sampling layers to obtain the target analysis area; perform cortical bone segmentation and medullary cavity segmentation on the target analysis area, use the medullary cavity segmentation results to form a medullary cavity closure mask, extract trabecular bone tissue within the medullary cavity closure mask, and output the quantitative parameters of bone tissue.

[0079] It should be noted that the anatomical zero-point layer number output in step S3 corresponds to the unified reference position of the growth plate's initial boundary in the standardized skeletal image sequence. By offsetting the anatomical zero point along the depth axis by a preset number of layers and truncating a preset number of sampling layers, the target analysis region can be stably mapped to the metaphysis analysis area below the growth plate, avoiding the direct inclusion of adjacent layers in the rapid reconstruction phase into the quantitative analysis. Furthermore, by performing cortical bone segmentation and medullary cavity segmentation on the target analysis region, and using the medullary cavity segmentation results to form a medullary cavity closure mask, the extraction of trabeculae is confined to the closed space inside the medullary cavity. This reduces the interference of high-grayscale boundaries of cortical bone on trabeculae identification, thereby improving the comparability and stability of subsequent bone tissue quantitative parameters.

[0080] Furthermore, a more specific implementation of step S4 is as follows: Obtain the anatomical zero-point layer number output in step S3 and map it along the depth axis to the target analysis region. Specifically, the offset layer number and sampling layer number are pre-written in the parameter configuration table, and the direction of depth coordinate increment is pre-defined as the direction from the epiphysis to the metaphysis. After reading the anatomical zero-point layer number, add the offset layer number to the anatomical zero-point layer number to obtain the starting layer number of the target analysis region; then add the sampling layer number minus 1 to the starting layer number to obtain the ending layer number of the target analysis region. If the ending layer number of the target analysis region exceeds the maximum layer number of the standardized skeletal image sequence, the maximum layer number is used as the ending layer number, and the sampling layer number is pushed back to determine the starting layer number; if the pushed-back starting layer number is less than 0, the starting layer number is set to 0. Taking one embodiment as an example, if the anatomical zero-point layer number is 151, the offset layer number is set to 50, the sampling layer number is set to 100, and the maximum layer number of the standardized bone image sequence is 395, then the starting layer number of the target analysis area is 201, the ending layer number is 300, corresponding to 100 consecutive cross-sectional slices; when the depth interlayer spacing is 0.01mm, the corresponding physical depth range is 2.01mm to 3.00mm.

[0081] Cortical bone segmentation and medullary cavity segmentation are performed on each cross-sectional slice in the target analysis region. Specifically, each cross-sectional slice in the target analysis region is input into a bone tissue segmentation network, which outputs a cortical bone probability map and a medullary cavity probability map. The bone tissue segmentation network is implemented using a shared encoder and a dual-segmentation output head structure. The encoder is used to extract bone tissue texture and contour features from the slices, the first segmentation output head is used to generate the cortical bone probability map, and the second segmentation output head is used to generate the medullary cavity probability map. As an example, in one implementation, the bone tissue segmentation network can be implemented using a dual-output semantic segmentation network based on a U-shaped encoder-decoder structure. The input image is first subjected to size unification and grayscale normalization processing before being input into the bone tissue segmentation network. In one embodiment, each cross-sectional slice is uniformly adjusted to 256×256 pixels, and the grayscale values ​​are linearly mapped to the 0-1 range. The cortical bone probability map and the medullary cavity probability map are binarized using 0.50 as the segmentation threshold to obtain the initial masks for cortical bone and medullary cavity.

[0082] To ensure that those skilled in the art can perform cortical bone segmentation and medullary cavity segmentation according to the specifications, a corresponding training sample set is first constructed. Specifically, cross-sectional slices located within the target analysis region are extracted from multiple standardized skeletal image sequences that have undergone pose alignment and 3D cropping. Annotators then perform pixel-level delineation of the cortical bone and medullary cavity regions in each slice, forming cortical bone and medullary cavity annotation masks. For example, in one embodiment, 9600 cross-sectional slices can be extracted from 120 sets of standardized skeletal image sequences, with 7680 slices in the training set, 960 in the validation set, and 960 in the test set. During the training phase, the input slices and the corresponding cortical bone and medullary cavity annotation masks are simultaneously input into the bone tissue segmentation network. Parameter updates are performed using a weighted sum of the cross-entropy loss function and the Dice loss function; the cortical bone segmentation loss weight is set to 1.2, and the medullary cavity segmentation loss weight is set to 1.0. The network was trained for a total of 250 epochs, with a batch size of 16. The AdamW optimizer was used, with an initial learning rate of 0.001, a minimum learning rate of 0.00001, and a weight decay coefficient of 0.0005. Linear warm-up was used for the first three epochs, followed by cosine annealing for learning rate scheduling. During training, random horizontal flipping, random scaling, and brightness perturbation could be enabled, with random scaling ratios set from 0.90 to 1.10 and brightness perturbation ratios set from -0.12 to 0.12. In the application phase, the trained bone segmentation network was directly used to generate initial masks for cortical bone and medullary cavity.

[0083] A closed mask for the medullary cavity is formed using the segmentation results. Specifically, connected component filtering, closing operations, hole filling, and boundary smoothing are performed on the initial medullary cavity mask to obtain the closed mask. Connected component filtering is used to remove isolated noise regions, and the connected component with the largest retained area is selected as the main connected component of the medullary cavity. Closing operations are used to repair local gaps at the boundary of the medullary cavity. Hole filling is used to fill isolated cavities outside the medullary cavity within the main connected component. Boundary smoothing is used to reduce the jaggedness of the mask boundary. Taking one embodiment as an example, a circular structuring element with a radius of 2 pixels is used to perform a closing operation on the initial medullary cavity mask, and then the largest connected component with a pixel area greater than 80 is retained. Hole filling of 8 neighboring areas is then performed to obtain the closed mask for the medullary cavity. If there is an overlapping area between the closed mask for the medullary cavity and the initial mask for the cortical bone, the overlapping area is deleted, and a boundary contraction of 1 pixel is performed again with the inner boundary of the cortical bone as a reference to ensure that the closed mask is located in the inner space of the cortical bone. Through the above processing, the closed mask for the medullary cavity forms a continuous, closed, and topologically complete internally constrained region.

[0084] Trabecular bone tissue was extracted within the closed mask of the medullary cavity. Specifically, for each slice of the target analysis region, only the pixel grayscale values ​​inside the closed mask were retained, while pixels outside the closed mask were masked; then, the trabecular bone segmentation threshold was calculated inside the closed mask. To improve the adaptability of trabecular bone extraction to grayscale fluctuations in different samples, a dual-threshold fusion method was used to determine the trabecular bone segmentation threshold: first, a first threshold (maximum inter-class variance threshold) was calculated based on the pixel grayscale histogram inside the closed mask using the maximum inter-class variance method; then, a second threshold was calculated based on the mean and standard deviation of pixel grayscale values ​​inside the closed mask, which is the sum of the mean pixel grayscale values ​​inside the closed mask multiplied by a preset factor and the standard deviation of pixel grayscale values ​​inside the closed mask; finally, the larger of the first and second thresholds was taken as the trabecular bone segmentation threshold for the current slice.

[0085] It should be noted that the pixels inside the closed mask mainly include bone marrow background pixels and trabecular bone pixels. Trabecular bone pixels are typically fewer in number than bone marrow background pixels, and their boundaries are susceptible to local noise, grayscale fluctuations, and cortical bone highlight artifacts. Therefore, this embodiment sets the second threshold to the average grayscale value of the pixels inside the closed mask plus 0.8 times the grayscale standard deviation. The average grayscale value represents the center of the background distribution inside the closed mask, and the grayscale standard deviation represents the degree of grayscale dispersion. By shifting the average value towards higher grayscale values ​​by a preset multiple of the standard deviation, the second threshold is preferentially placed outside the upper edge of the background grayscale distribution, thereby reducing the misidentification of high-grayscale background noise and boundary artifacts as trabecular bone. Furthermore, 0.8 times is a preferred coefficient, achieving a balance between suppressing false positives and avoiding the omission of small trabecular bone. When this coefficient is too small, the threshold rise is insufficient, easily introducing high-grayscale noise from non-trabecular bone; when this coefficient is too large, the threshold rise is too high, easily missing true trabecular bone with moderate grayscale. Therefore, setting the second threshold to the average gray level of pixels inside the closed mask plus 0.8 times the gray level standard deviation threshold can, together with the maximum inter-class variance threshold, form a dual-threshold fusion method, improving the stability of trabecular bone extraction. The multiplier can also be set within the range of 0.6 to 1.0 based on the dispersion of gray level distribution inside the closed mask.

[0086] Pixels with gray values ​​greater than or equal to the trabecular segmentation threshold inside the closed mask are marked as trabecular pixels; pixels with gray values ​​less than the trabecular segmentation threshold inside the closed mask are marked as non-trabecular pixels. For example, in a cross-sectional slice corresponding to depth layer 210, the gray-level histogram inside the closed mask yields a first threshold of 128, a mean gray value of 96, and a standard deviation of 48. Based on this, a second threshold of 134.4 is obtained. 134.4 is then rounded down to 135 and used as the trabecular segmentation threshold for the current slice. When a pixel within the closed mask has a gray value of 142, it is marked as a trabecular pixel; when another pixel has a gray value of 118, it is marked as a non-trabecular pixel.

[0087] To reduce the interference of noise points and single-layer artifacts on the trabecular bone extraction results, cross-layer connectivity consistency filtering was performed on the initial trabecular bone results in all slices. Specifically, the initial trabecular bone results of each slice in the target analysis region were stacked along the depth axis into a three-dimensional trabecular bone voxel set. A 26-neighborhood connected component analysis was performed on the three-dimensional trabecular bone voxel set, and isolated connected components with fewer than 20 voxels were deleted, retaining the remaining connected components as the final trabecular bone extraction results. This combined processing of "double-threshold fusion segmentation within a closed mask + cross-layer connectivity consistency filtering" not only limited the physical distribution space of the trabecular bone but also enhanced the continuity of the trabecular bone recognition results in the inter-layer direction.

[0088] Output bone tissue quantification parameters. Specifically, based on the initial cortical bone mask, the medullary cavity closure mask, and the final extraction results of trabeculae, the number of pixels in cortical bone, the number of pixels in the medullary cavity closure mask, and the number of pixels in trabeculae are counted in each slice, and the area and volume parameters are calculated by combining the pixel spacing and depth layer spacing. For example, when the horizontal and vertical pixel spacing are both 0.01 mm, and the depth layer spacing is 0.01 mm, the volume corresponding to a single voxel is 0.000001 mm. 3 In the cross-sectional slice corresponding to depth layer 210, if the number of pixels in the cortical bone is 4300, the number of pixels in the medullary canal closure mask is 8200, and the number of pixels in the trabecular bone is 1650, then the area of ​​the cortical bone is 0.43 mm². 2 The area of ​​the medullary cavity is 0.82 mm². 2 The trabecular area is 0.165 mm². 2 After statistical analysis of all 100 slices in the target analysis area, if the total number of voxels in the bone marrow cavity closure mask is 820,000 and the total number of voxels in the trabecular bone is 165,000, then the total volume of the bone marrow cavity is 0.82 mm. 3 The total volume of the trabecular bone is 0.165 mm. 3The trabecular bone volume fraction was 20.12%. Further, the average cortical bone area of ​​all sections was calculated, and the average medullary cavity area of ​​all sections was calculated, yielding the average medullary cavity area. The ratio of the total trabecular bone volume to the total medullary cavity volume was output as one of the quantitative parameters of bone tissue in the target analysis region. Finally, the anatomical zero-point layer number, the starting layer number of the target analysis region, the ending layer number of the target analysis region, the average cortical bone area, the average medullary cavity area, the total trabecular bone volume, and the trabecular bone volume fraction were entered into the quantitative analysis results table to form the quantitative parameters of bone tissue for the current sample.

[0089] Through the above processing, step S4 achieves a quantitative workflow for the target analysis region, using the anatomical zero point as a unified reference, the medullary cavity closure mask as the internal constraint boundary, and the extraction of bone trabeculae within the closure mask as the core. This workflow sequentially connects target analysis region mapping, bone tissue segmentation, closure mask construction, bone trabeculae extraction, and quantitative parameter output, providing a unified analysis region scope and stable quantitative parameter output results across different samples. This reduces statistical errors caused by positioning bias, cortical bone interference, and inter-slice noise.

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

[0091] Finally: The above description is only a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. An automatic analysis method for CT images of long bones in rats and mice, characterized in that, include: Acquire the target long bone CT image sequence, extract the three-dimensional spatial distribution of bone tissue, and perform pose alignment and three-dimensional cropping based on the main direction of bone tissue to obtain a standardized bone image sequence; A first candidate depth interval is constructed based on the butterfly-shaped anatomical features in the cross-sectional slices, and a second candidate depth interval is constructed based on the depth coordinates of key points in the growth plate-like gaps in the longitudinal slices. An intersection operation is performed on the first candidate depth interval and the second candidate depth interval to obtain a fine search range. Within a fine search range, candidate slice sequences of triangular fractures are identified where epiphyseal and metaphyseal bone tissues transition from a contact state to a disconnected state. The geometric gap distribution sequence of the epiphyseal and metaphyseal edges in each candidate slice of triangular fracture is calculated, and the slice layer number with the smallest variance of the geometric gap distribution sequence is determined as the anatomical zero point. The target analysis area is obtained by offsetting the anatomical zero point along the depth axis by a preset number of layers and extracting a preset number of sampling layers. The target analysis region is segmented into cortical bone and medullary cavity. The medullary cavity segmentation results are used to form a medullary cavity closure mask. Bone trabecular tissue is extracted within the medullary cavity closure mask, and bone tissue quantitative parameters are output.

2. The automatic analysis method for long bone CT images of mice and rats according to claim 1, characterized in that, Methods for obtaining the target analysis region include: The anatomical zero-point layer number, as well as the offset layer number and sampling layer number pre-written into the parameter configuration table, are read. The anatomical zero-point layer number is added to the offset layer number to determine the starting layer number of the target analysis region. The starting layer number is added to the sampling layer number minus 1 to determine the ending layer number of the target analysis region. When the ending layer number exceeds the maximum layer number of the standardized skeletal image sequence, the maximum layer number is used as the ending layer number, and the sampling layer number is pushed back to determine the starting layer number. When the pushed-back starting layer number is less than 0, the starting layer number is set to 0.

3. The automatic analysis method for long bone CT images of mice and rats according to claim 1, characterized in that, The target analysis region is segmented into cortical bone and medullary cavity. The medullary cavity segmentation results are used to form a medullary cavity closure mask, including: Each cross-sectional slice in the target analysis region is input into the bone tissue segmentation network, which is implemented using a shared encoder and a dual segmentation output head structure. The encoder is used to extract the texture and contour features of bone tissue in the slices. The first segmentation output head is used to generate a cortical bone probability map, and the second segmentation output head is used to generate a medullary cavity probability map. Threshold binarization is performed on the cortical bone probability map and the medullary cavity probability map respectively to obtain the initial mask of cortical bone and the initial mask of medullary cavity. Connected component filtering, closing operation, hole filling and boundary smoothing are performed on the initial mask of medullary cavity to obtain the closed mask of medullary cavity. After deleting the overlapping area with the initial mask of cortical bone, boundary contraction processing is performed based on the inner boundary of cortical bone.

4. The automatic analysis method for long bone CT images of mice and rats according to claim 1, characterized in that, Trabecular bone tissue was extracted within the closed mask of the medullary cavity, and quantitative parameters of the bone tissue were output, including: Only the pixel grayscale values ​​inside the closed mask of the medullary cavity are retained. A dual threshold fusion method is used to determine the trabecular segmentation threshold. The first threshold is the maximum inter-class variance threshold calculated based on the pixel grayscale histogram inside the closed mask. The second threshold is the sum of the mean pixel grayscale value inside the closed mask and a preset multiple multiplied by the standard deviation of pixel grayscale value inside the closed mask. The larger value between the first threshold and the second threshold is taken as the trabecular segmentation threshold. Pixels with gray values ​​greater than or equal to the trabecular segmentation threshold are marked as trabecular pixels. Trabecular pixels from each slice within the target analysis region are stacked along the depth axis to form a three-dimensional trabecular voxel set. A 26-neighborhood connected component analysis is performed on the three-dimensional trabecular voxel set. Connected components with a voxel count not less than a preset threshold are retained as the trabecular extraction results. The trabecular volume, medullary cavity volume, and trabecular volume fraction are calculated by combining the pixel spacing and depth layer spacing.

5. The automatic analysis method for long bone CT images of mice and rats according to claim 1, characterized in that, Methods for obtaining candidate slice sequences of triangular fracture include: Each cross-sectional slice within the fine search range is input into the critical feature detection neural network model, which outputs candidate target boxes, state category indexes, and state confidence scores. The state category indexes include contact state and disconnected state. Non-maximum suppression is performed on the candidate target boxes. Based on the state confidence scores and preset state thresholds, each cross-sectional slice is marked as a contact state slice, a disconnected state slice, or a missed slice. The interlayer position where the slice transitions from a contact state slice to a disconnected state slice is searched along the depth layer number from small to large. Consecutive adjacent disconnected state slices are identified as a triangular fracture candidate slice sequence.

6. The automatic analysis method for long bone CT images of mice and rats according to claim 5, characterized in that, Methods for obtaining the geometric gap distribution sequence include: Based on the effective target boxes corresponding to each candidate slice in the triangular fracture candidate slice sequence, the local analysis region is determined. The lowest boundary point of bone tissue above the local low-density gap is extracted to form the epiphyseal edge curve, and the highest boundary point of bone tissue below the local low-density gap is extracted to form the metaphyseal edge curve. With the central gap position as the center, the preset sampling width and preset number of sampling points are set, and the edge spacing between the epiphyseal edge and the metaphyseal edge at each sampling position is calculated to form a geometric gap distribution sequence.

7. The automatic analysis method for long bone CT images of mice and rats according to claim 1, characterized in that, The methods for constructing the first candidate depth interval include: A preset retrieval step size is set for cross-sectional slices. Starting from the initial depth layer of the standardized skeletal image sequence, cross-sectional slices are extracted sequentially according to the preset retrieval step size. For each cross-sectional slice, size uniformity, grayscale normalization, and butterfly-shaped anatomical feature recognition are performed. Confidence ranking and non-maximum suppression are performed on the candidate detection boxes in the recognition results. The highest confidence score in the valid detection results is read. When the highest confidence score is greater than or equal to the preset judgment threshold, the corresponding depth layer number is recorded to form a butterfly-shaped feature hit layer number set. The first candidate depth interval is determined based on the minimum hit layer number, the maximum hit layer number, and the preset boundary compensation layer number.

8. The automatic analysis method for long bone CT images of mice and rats according to claim 1, characterized in that, The methods for constructing the second candidate depth interval include: In a standardized skeletal image sequence, longitudinal slices are reconstructed by fixing the horizontal or vertical coordinates. Keypoint identification is performed on the growth plate-like gaps in the longitudinal slices. The coordinate values ​​of keypoints with confidence greater than or equal to a preset keypoint threshold in the depth direction are extracted and mapped to depth layer numbers to form a set of keypoint depth layer numbers. The second candidate depth interval is determined based on the minimum keypoint layer number, the maximum keypoint layer number, and the preset keypoint boundary compensation layer number.

9. The automatic analysis method for long bone CT images of mice and rats according to claim 1, characterized in that, Methods for obtaining standardized skeletal image sequences include: A three-dimensional spatial coordinate matrix is ​​constructed based on the voxel coordinates of bone tissue in the three-dimensional spatial distribution of bone tissue. The coordinates of the spatial center are determined. Covariance analysis and eigenvalue decomposition are performed on the centered three-dimensional spatial coordinate matrix to extract the principal direction vector of bone tissue. A dot product operation is performed on the principal direction vector of bone tissue and the unit direction vector of the preset depth axis, and the skew angle of the principal direction is determined by the inverse cosine operation. A cross product operation is performed on the principal direction vector of bone tissue and the unit direction vector of the preset depth axis to determine the rotation axis direction vector. A rotation matrix is ​​generated based on the rotation axis direction vector, the skew angle of the principal direction, and the Rodriguez rotation rule. Based on the rotation matrix, rotation transformation and interpolation resampling are performed on the voxel coordinates to generate pose-aligned 3D volume data. Based on the pose-aligned 3D volume data, 3D clipping is performed to obtain a standardized bone image sequence.

10. The automatic analysis method for CT images of long bones in mice and rats according to claim 9, characterized in that, Methods for generating rotation matrices based on rotation axis direction vectors, principal direction skew angles, and Rodriguez rotation rules include: The rotation axis direction vector is normalized to obtain the first, second, and third components of the unit rotation axis. Calculate the magnitude of the rotation axis direction vector; when the magnitude is less than a preset minimum threshold, set the rotation matrix as an identity matrix; when the magnitude of the rotation axis is not less than the preset minimum threshold, calculate the cosine and sine values ​​of the principal direction skew angle, and based on the first, second, and third components of the unit rotation axis, the cosine value of the principal direction skew angle, and the sine value of the principal direction skew angle, determine the elements at each position of the three-row, three-column rotation matrix according to the Rodriguez rotation rule.

Images