Dynamic X-ray ankle complex bony mark point analysis method and device

Abstract

The invention provides a dynamic X-ray ankle complex bony mark point analysis method and device, and belongs to the field of mark point analysis, and the method comprises the steps: obtaining a plurality of X-ray video samples; extracting frame by frame, and marking the bony mark points to obtain a label data set and a label-free data set; constructing an attitude estimation network and carrying out semi-supervised training; extracting a to-be-analyzed dynamic X-ray video frame by frame to obtain a reasoning image sequence, and obtaining a preliminary key point position sequence and motion phase probability distribution of a current frame; performing time sequence smoothing processing on the initial key point position sequence; performing skeleton rigid body characteristic inspection and correction on the filtered key point position sequence to obtain an optimized key point position sequence; kinematics parameter time sequence data are obtained through calculation based on the optimized key point position sequence, and a time history change analysis chart and a cooperative relation analysis chart are generated according to the kinematics parameter time sequence data.

Description

Technical Field

[0001] This invention relates to the field of landmark analysis technology, and in particular to a method and apparatus for dynamic X-ray analysis of bony landmarks in the foot and ankle complex. Background Technology

[0002] The foot and ankle complex, as a crucial weight-bearing and movement hub of the lower limbs, directly impacts an individual's mobility and quality of life. Composed of multiple bones and a complex network of joints, the precise analysis of its coordinated movement function is of paramount importance in clinical diagnosis, treatment planning, orthotic design, postoperative rehabilitation assessment, and biomechanical research. Currently, dynamic motion analysis techniques for the human foot and ankle primarily include: optical motion capture technology based on surface markers, medical imaging techniques based on static computed tomography (CT / MRI), and dynamic imaging analysis techniques based on X-ray fluoroscopy.

[0003] Existing technologies have many limitations when applied to dynamic analysis of the foot and ankle. Optical motion capture technology based on surface markers is affected by "skin motion artifacts" and cannot directly measure the actual movement of deep bones inside the foot and ankle. Although CT and MRI technologies can generate high-resolution bone models, they can only acquire non-weight-bearing static data and cannot capture continuous dynamic responses during real functional activities. Dual-plane dynamic X-ray imaging systems (DFIS) can achieve dynamic capture, but they rely on complex 3D-2D registration processes, requiring additional CT / MRI scans and a large amount of manual intervention, making them technically challenging and time-consuming. Single-plane dynamic DR systems theoretically have the potential to capture skeletal motion trajectories, but current analysis methods rely on time-consuming manual annotation and calculation, making it difficult to achieve efficient synchronous analysis of multiple joints.

[0004] Chinese patent application CN115018977A discloses a joint motion analysis method based on biplane X-rays. This method involves reconstructing a three-dimensional skeletal model from a pre-existing CT scan of the subject, then registering the three-dimensional model with biplane X-ray images, and finally solving for the three-dimensional motion parameters of the joint through phased optimization. However, this method involves manual registration of the initial frame, multi-stage optimization iterations, and manual intervention when automatic matching is ineffective, making the overall process complex and limiting its efficiency and cost in clinical application. Summary of the Invention

[0005] The technical solution of this invention is implemented as follows: On the one hand, the present invention provides a method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex, comprising the following steps: S1. The foot and ankle complex of the subject is continuously scanned in a preset motion mode using a dynamic DR system to obtain multiple X-ray video samples. S2. Frame by frame, extract multiple X-ray video samples to obtain dynamic X-ray image sequences of the foot and ankle complex. Label the bony landmarks on multiple preset bones in the dynamic X-ray image sequences to obtain a labeled dataset. Use the remaining unlabeled dynamic X-ray images as an unlabeled dataset. S3. Construct a pose estimation network by mixing the labeled and unlabeled datasets in proportion to form a training set, inputting it into the pose estimation network for semi-supervised training, and obtaining the trained pose estimation network model. S4. Extract the inference image sequence frame by frame from the dynamic X-ray video to be analyzed, and input the inference image sequence into the trained pose estimation network model to predict frame by frame, so as to obtain the preliminary key point position sequence and the motion phase probability distribution of the current frame. S5. Input the preliminary key point position sequence into the Kalman filter, dynamically adjust the process noise covariance and measurement noise covariance parameters of the Kalman filter according to the motion phase probability distribution, and perform time-series smoothing on the preliminary key point position sequence to obtain the filtered key point position sequence. S6. Perform bone rigidity property verification and correction on the filtered key point position sequence, calculate the deviation between the bone segment length and the preset reference bone length in each frame of dynamic X-ray image, and perform geometric correction when the deviation exceeds the tolerance threshold set according to the motion phase to obtain the optimized key point position sequence. S7. Based on the optimized key point position sequence, calculate the rotation angle of each skeletal segment relative to the preset reference system and the relative distance between the geometric center points of different skeletal segments in each frame of dynamic X-ray image to obtain kinematic parameter time series data, and generate time history change analysis charts and synergistic relationship analysis charts based on the kinematic parameter time series data.

[0006] Based on the above technical solutions, preferably, the pose estimation network includes a feature extraction backbone network, a static analysis branch, a temporal analysis branch, and a lightweight phase classification branch, wherein, The static analysis branch is used to process the initial key point location prediction for a single frame image; The temporal analysis branch is used to process continuous image sequences containing the target frame and its preceding and following frames, and captures the motion correlation of key points in the time dimension through the temporal information processing module and outputs corrected key point position predictions. A lightweight phase classification branch is used to predict gait motion phase.

[0007] Based on the above technical solutions, preferably, the pose estimation network model uses a composite loss function for joint optimization, wherein the composite loss function includes supervised loss, unsupervised temporal consistency loss, unsupervised pose structure loss, skeleton rigid body constraint loss, and motion phase recognition auxiliary loss. The supervised loss is obtained by comparing the mean squared error between the prediction results of the static analysis branch on the labeled dataset and the manually labeled true values; The unsupervised temporal consistency loss calculation is based on the position change of the same key point predicted by the pose estimation network between two temporally adjacent frames. When the position change exceeds the displacement threshold dynamically adjusted according to the motion phase, the part exceeding the threshold is included in the loss term. The unsupervised pose structure loss is based on a pose statistical prior model pre-established from the labeled dataset through principal component analysis. The reconstruction error between the predicted keypoint pose and the pose reconstructed by the pose statistical prior model is calculated and included in the loss term. The rigid body constraint loss of the skeleton is calculated by the deviation between the predicted skeleton segment length and the reference length of the foot and ankle skeleton segment, and the deviation is weighted according to the weight of motion phase modulation. The weighted deviation is included in the loss term. The reference length of the foot and ankle skeleton segment is the median of the lengths of each skeleton segment statistically obtained from the label dataset. The motion phase recognition auxiliary loss lightweight phase classification branch predicts the phase of the current frame in the gait cycle, and calculates the phase by cross-entropy loss with the gait phase pseudo-label automatically generated based on the vertical displacement change of key points. The phase of the gait cycle includes the standing phase, the swinging phase, and the transition phase.

[0008] Based on the above technical solutions, preferably, the unsupervised temporal consistency loss dynamically adjusts the displacement threshold according to the motion phase, specifically including: The phase probability distribution of the current frame is obtained by predicting the phase of the gait cycle in the current frame based on the key point sequence features extracted by the temporal analysis branch using a lightweight phase classification branch. This phase probability distribution includes the standing phase probability. Oscillation probability and transition phase probability ; Set the standing phase weight coefficient Oscillating phase weighting coefficient and transition phase weighting coefficient ,in ; The dynamic displacement threshold is calculated based on the phase probability distribution and weighting coefficients. The position change of the same key point between adjacent frames is calculated. When the position change exceeds the dynamic displacement threshold, the square of the excess is included in the calculation. The temporal consistency loss at this key point in frame number 1; where the dynamic displacement threshold is calculated using the following formula: ; in, Indicates the first The dynamic displacement threshold of the frame. Indicates the reference displacement threshold. Indicates the first The standing phase probability of a frame. Indicates the first The frame swing phase probability, Indicates the first The probability of the transition phase of a frame; The unsupervised temporal consistency loss is obtained by averaging the frames within all time windows of all keypoints and training sets.

[0009] Based on the above technical solutions, preferably, the calculation steps for the rigid body constraint loss of the skeleton specifically include: The lengths of skeletal segments in the labeled frames of the labeled dataset are statistically analyzed, and the median length of each skeletal segment is calculated as the reference length for that bone. , where i represents the skeleton index; During training, the length of each skeletal segment is calculated based on the keypoint locations predicted by the pose estimation network on each frame of dynamic X-ray image. The skeletal segment is formed by connecting two key points on the same bone. Indicates the frame index; Calculate the deviation between the predicted bone segment length and the preset reference bone length for each frame; The phase modulation weights are calculated based on the motion phase probability distribution, and then weighted using these phase modulation weights to account for the skeletal segment length deviations, resulting in the [number of]th [item / section]. The skeletal rigid body constraint loss of each frame is calculated by averaging the losses across all frames in the training set; the calculation formula is as follows: ; ; in, Indicates the first Frame phase modulation weights, This represents the weighting coefficient of the rigid body constraint in the standing phase. This represents the constraint weight coefficient of the oscillating phase rigid body. This represents the weighting coefficient of the rigid body constraint during the transition phase. Indicates the first The frame's skeletal rigid body constraint loss, where N represents the number of bones. .

[0010] Based on the above technical solutions, preferably, the frame-by-frame prediction specifically includes: Determine the position of the target frame to be predicted in the inference image sequence. When the target frame is located near the beginning or end of the inference image sequence and a complete temporal analysis window containing 5 frames (2 frames before and 2 frames after) cannot be constructed, only the static analysis branch is used to predict the target frame, and the prediction result of the static analysis branch is used as the preliminary key point position of the target frame. When the target frame can construct a complete temporal analysis window containing two frames before and after it, the target frame is predicted using both static analysis branch and temporal analysis branch to obtain static prediction results and temporal prediction results respectively. The peak intensity of the heatmap corresponding to the static prediction result and the temporal prediction result is calculated as the confidence score, and the two prediction results are weighted and averaged to obtain the preliminary key point position of the target frame. The lightweight phase classification branch outputs the motion phase probability distribution of the target frame while the timing analysis branch outputs the timing prediction results.

[0011] Based on the above technical solutions, preferably, the step of dynamically adjusting the process noise covariance and measurement noise covariance parameters of the Kalman filter according to the motion phase probability distribution specifically includes: The state vector of the Kalman filter is defined to include the position, velocity, and acceleration of the key points, and a constant acceleration motion model is used for state transition. The standard deviation of process noise and the standard deviation of measurement noise are calculated based on the motion phase probability distribution. The calculation formula is as follows: ; ; in, Indicates the first Standard deviation of frame process noise This represents the standard deviation parameter of noise during the standing phase process. This represents the standard deviation parameter of the noise in the oscillating phase process. This parameter represents the standard deviation of noise during the transition phase process. Indicates the first Frame measurement noise standard, This represents the standard deviation parameter of noise measured during standing phase. This represents the standard deviation parameter of the measured noise in the oscillating phase. This represents the standard deviation parameter of the measurement noise during the transition phase. , ; Using the squares of the process noise standard deviation and the measurement noise standard deviation as the diagonal elements of the process noise covariance matrix and the measurement noise covariance matrix, respectively, we obtain the... Frame process noise covariance matrix and measurement noise covariance matrix ; The process noise covariance matrix and measurement noise covariance matrix The prediction and update steps of the Kalman filter are recursively calculated. The prediction step predicts the current state based on the state transition model and the process noise covariance matrix. The update step calculates the optimal estimated state based on the observed values, the measurement noise covariance matrix, and the predicted state. The position component in the optimal estimated state is extracted as the filtered key point position.

[0012] Based on the above technical solutions, preferably, the specific steps for testing and correcting the rigid body characteristics of the skeleton include: Step a: Calculate the length of each skeletal segment in each frame based on the filtered keypoint position sequence. The skeletal segment is formed by connecting two key points on the same bone; Step b: Calculate the length of the bone segment and the preset reference bone length. deviation ; Step c, according to the first The dominant motion phase is determined by the motion phase probability distribution of the frame. The dominant motion phase is the phase with the highest probability value in the motion phase probability distribution. A tolerance threshold is set based on the dominant motion phase of the current frame. The tolerance threshold for standing is less than the tolerance threshold for swinging. Step d, for the first The first frame Root bones, when deviation Exceeding the tolerance threshold If the skeletal segment violates rigid body properties, the direction vector of the skeletal segment from the starting point to the ending point is calculated, and the normalized direction vector is multiplied by the reference bone length. The corrected skeletal segment vector is obtained. Based on the corrected skeletal segment vector, the coordinates of the key points of the two endpoints of the skeletal segment are updated using a partial correction method. The partial correction method is: keeping the position of the midpoint of the skeletal segment unchanged, and symmetrically adjusting the two endpoints to the position that conforms to the reference bone length along the direction vector. Step e: Repeat step ad for all bones in all frames to obtain the optimized key point position sequence after rigid body property verification and correction.

[0013] Furthermore, step S7 specifically includes: S71. Determine the set of bones and the set of bone pairs for quantitative analysis. The set of bones includes each bone whose rotation angle needs to be calculated. The set of bone pairs includes each bone pair whose center distance needs to be calculated. Establish a two-dimensional Cartesian coordinate system for each frame of dynamic X-ray image. The origin is located at the lower left corner of the image. The X-axis is horizontal to the right and the Y-axis is vertical to the up. S72. In the first frame of dynamic X-ray image, the larger Y-coordinate value of the two endpoints of each bone segment in the skeleton set is determined as the origin of rotation of that bone segment. S73. In each frame of dynamic X-ray image, establish a local coordinate system parallel to the two-dimensional Cartesian coordinate system with the origin of rotation, and calculate the angle between the positive X-axis direction of the local coordinate system and the vector formed by the bone line segment, which is used as the rotation angle of the bone in the current frame. S74. Repeat steps S71-S73 for each bone in each frame of dynamic X-ray image to obtain the set of rotation angles of each bone in each frame of dynamic X-ray image. S75. For each skeletal line segment in each frame of dynamic X-ray image, calculate the arithmetic mean of the coordinates of its two endpoints as the geometric center point of the skeletal line segment. S76. For each bone pair in the set of bone pairs, the Euclidean distance between the geometric center points of the two bone segments in the bone pair is taken as the center distance of the bone pair in the frame. S77. Repeat steps S75-S76 for all selected bone pairs to obtain the set of center distances of each bone pair in each frame of dynamic X-ray image. S78. Organize the set of rotation angles of each bone in each frame of dynamic X-ray image and the set of center distances of each bone pair in each frame of dynamic X-ray image in chronological order to form kinematic parameter time series data. Generate time history change analysis charts and synergy analysis charts based on kinematic parameter time series data.

[0014] On the other hand, the present invention provides a dynamic X-ray foot-ankle complex bony landmark analysis device, employing the dynamic X-ray foot-ankle complex bony landmark analysis method as described above, including: The data acquisition module is used to continuously scan the foot and ankle complex of the subject in a preset movement mode through a dynamic DR system to obtain multiple X-ray video samples. The data preprocessing module is used to extract multiple X-ray video samples frame by frame to obtain a dynamic X-ray image sequence of the foot and ankle complex. It also labels multiple pre-set bony landmarks on bones in the dynamic X-ray image sequence to obtain a labeled dataset. The remaining unlabeled dynamic X-ray images are used as an unlabeled dataset. The model training module is used to build the pose estimation network. It mixes the labeled and unlabeled datasets in a certain proportion to form a training set, which is then input into the pose estimation network for semi-supervised training to obtain the trained pose estimation network model. The key point prediction module is used to extract the inference image sequence frame by frame from the dynamic X-ray video to be analyzed, and input the inference image sequence into the trained pose estimation network model to predict frame by frame, so as to obtain the preliminary key point position sequence and the motion phase probability distribution of the current frame. An adaptive filtering module is used to input the preliminary key point position sequence into a Kalman filter, dynamically adjust the process noise covariance and measurement noise covariance parameters of the Kalman filter according to the motion phase probability distribution, and perform time-series smoothing on the preliminary key point position sequence to obtain the filtered key point position sequence. The posterior verification module is used to check and correct the rigid body characteristics of the skeleton in the filtered key point position sequence. It calculates the deviation between the length of the bone line segment and the preset reference bone length in each frame of dynamic X-ray image. When the deviation exceeds the tolerance threshold set according to the motion phase, geometric correction is performed to obtain the optimized key point position sequence. The analysis module is used to calculate the rotation angle of each skeletal segment relative to a preset reference frame and the relative distance between the geometric center points of different skeletal segments in each frame of dynamic X-ray image based on the optimized key point position sequence, to obtain kinematic parameter time series data, and to generate time history change analysis charts and synergy analysis charts based on the kinematic parameter time series data.

[0015] The dynamic X-ray foot-ankle complex bony landmark analysis method and device of the present invention have the following advantages over the prior art: (1) This invention constructs a pose estimation network based on semi-supervised learning to automatically identify bony landmarks of the foot and ankle complex in dynamic DR video. It combines phase-aware adaptive Kalman filtering with posterior verification of rigid body properties of bones to avoid manual intervention. By calculating kinematic parameters such as bone rotation angle and interosseous center distance, it provides an efficient and accurate quantitative analysis tool for clinical diagnosis, rehabilitation assessment and biomechanical research. (2) This invention learns the temporal motion patterns and anatomical constraint features of key points from a large number of unlabeled dynamic X-ray images by designing unsupervised loss terms such as temporal consistency loss, posture structure loss and skeletal rigid body constraint loss, thereby reducing the number of image annotations and significantly reducing the dependence on annotations; (3) This invention predicts the phase information of gait motion through a lightweight phase classification branch and dynamically adjusts the process noise covariance and measurement noise covariance parameters of the Kalman filter according to the phase probability distribution. It can dynamically balance the weights of the model prediction and the observation value according to the change of motion state, which effectively suppresses the inter-frame jitter phenomenon and avoids the loss of real motion features caused by excessive smoothing, thus improving the stability and accuracy of key point tracking. (4) This invention utilizes the physical property that the length of bones remains constant in a short period of time as an unsupervised constraint. By calculating the deviation between the predicted bone segment length and the reference bone length statistically obtained from the labeled data, and using phase modulation weights for weighting, the network’s dependence on labeled data is effectively reduced, and unreasonable situations such as drastic changes in the length of the same bone between adjacent frames are avoided in the prediction results. This improves the model’s generalization ability and prediction accuracy on unlabeled data. (5) This invention calculates the deviation between the length of the skeletal line segment and the reference bone length in each frame, and sets a tolerance threshold according to the motion phase. It performs geometric correction on the skeletal line segments that exceed the threshold, which can effectively correct the violation of rigid body characteristics caused by accumulated errors or abnormal frames. It ensures that the final output key point position sequence satisfies both temporal smoothness and anatomical geometric invariance constraints, providing a high-quality data foundation for kinematic parameter calculation. Attached Figure Description

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

[0017] Figure 1 This is a flowchart of a dynamic X-ray foot-ankle complex bony landmark analysis method according to the present invention; Figure 2 This is an example of the geometric line segment definition of the target bone in the dynamic X-ray foot-ankle complex bony landmark analysis method of the present invention; Figure 3 This is a flowchart of the posture estimation model training process for a dynamic X-ray foot-ankle complex bony landmark analysis method according to the present invention. Figure 4 This is an example of defining the rotation angle of the target bone in a dynamic X-ray foot-ankle complex bony landmark analysis method of the present invention; Figure 5 This is an example of defining the center distance between target bones in a dynamic X-ray foot-ankle complex bony landmark analysis method of the present invention; Figure 6This is a schematic diagram illustrating the analysis of the rotation angle over time and the synergistic analysis of relative angles in a dynamic X-ray foot and ankle complex bony landmark analysis method of the present invention, taking the tibia and medial cuneiform as examples. Figure 7 This is a schematic diagram illustrating the temporal variation and relative positional synergistic analysis of the interosseous distance between the tibia and calcaneus and between the tibia and talus, as examples of the dynamic X-ray foot and ankle complex bony landmark analysis method of the present invention.

[0018] Figure label: 1- L 胫骨 ;2- L 距骨 ;3- L 跟骨 ;4- L 楔骨 ;5- L 跖骨 6-α tibia ;7-α talus ;8-α calcaneus ;9-α cuneiform ;10-α metatarsal . Detailed Implementation

[0019] The technical solutions of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0020] like Figure 1 As shown, this invention provides a method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex, comprising the following steps: S1. The foot and ankle complex of the subject is continuously scanned in a preset movement mode using a dynamic DR system to obtain multiple X-ray video samples.

[0021] Understandably, the dynamic DR system consists of a fluoroscopy device and an adjustable test platform. The height of the test platform is matched to the X-ray beam position to ensure the stability of the subject's movement on the platform. The fluoroscopy pulse width is 8 milliseconds, the imaging pixel size is 1024×1024, and the image intensifier diameter is 12 inches. Typical exposure parameters for foot-ankle complex scans are 90-110 kVp voltage and 0.5-1.7 mA current, using the system's automatic exposure mode. The X-ray tube distance is 1-1.1 meters, using a grid ratio of 10:1 and a grid focal length of 1 meter. To eliminate artificial motion blur in dynamic images, window averaging and noise reduction functions are disabled. Subjects routinely wear protective equipment during image acquisition.

[0022] Preferably, the preset motion mode may include at least one of the following: Normal pace walking: Subjects walked naturally on the walking path at a comfortable pace of their own choosing.

[0023] Weighted walking: Subjects wore a weighted vest with symmetrical weight distribution in the front and back, with a weight of 50% of the subject's own body weight (including the weight of protective equipment such as lead apron).

[0024] Heel lift-lowering exercise: Subjects performed the complete process of “heel fully lifted - forefoot landing for support - heel landing smoothly”.

[0025] Both feet of the subjects were tested. A fluoroscopic device captured a horizontal view of the foot-ankle complex at a frame rate of 30 frames per second. Subjects were trained before the test. During the walking test, subjects walked three steps continuously at a low speed of approximately 0.5 meters per second to avoid gait abnormalities and imbalances. During the second step, a dynamic DR system captured a horizontal view of the foot-ankle complex.

[0026] S2. Frame by frame, extract multiple X-ray video samples to obtain dynamic X-ray image sequences of the foot and ankle complex. Label the bony landmarks on multiple preset bones in the dynamic X-ray image sequences to obtain a labeled dataset. Use the remaining unlabeled dynamic X-ray images as an unlabeled dataset. Understandably, dynamic X-ray images of the foot and ankle complex are obtained by extracting image sequences frame by frame from a single video. At least two senior orthopedic surgeons or radiology experts then use LabelMe software to annotate multiple frames in the image sequence, jointly and precisely marking two anatomical feature points on several pre-defined key bones to construct an annotated dataset containing image data and the coordinates of the corresponding feature points.

[0027] The anatomical feature points can be defined as follows: Tibia: T1 - the anterior inferior border of the weight-bearing articular surface of the distal tibia; T2 - the posterior inferior border of the weight-bearing articular surface of the distal tibia.

[0028] Talus: A1 - the anterior vertex of the talonavicular joint surface formed by the head of the talus and the navicular bone; A2 - the posterior vertex of the posterior process of the talus.

[0029] Calcaneus: C1 - the origin of the Gissane critical angle; C2 - the anterior vertex of the calcaneal protuberance.

[0030] Medial cuneiform: MC1 - the dorsal (superior) apex of the articular surface formed by it and the navicular bone; MC2 - the dorsal (superior) apex of the articular surface formed by it and the first metatarsal bone.

[0031] Fifth metatarsal: M5-1 - the tuberosity at the base of the fifth metatarsal; M5-2 - the distal end of the head of the fifth metatarsal.

[0032] refer to Figure 2 , Figure 2 Example of defining geometric line segments for the target skeleton. After accurately annotating the key points mentioned above in each frame of the image, the software connects two key points on the same skeleton to form a bone line segment L representing its geometric direction. In this embodiment, five bone line segments are created for five target skeletons: L 胫骨 (L) tibia L 距骨 (L) talus L 跟骨 (L) calcaneus L 楔骨 (L) cuneiform L 跖骨 (L) metatarsal ).

[0033] like Figure 3 As shown in Figure S3, construct a pose estimation network by mixing the labeled and unlabeled datasets in a certain proportion to form a training set, and inputting it into the pose estimation network for semi-supervised training to obtain the trained pose estimation network model. The pose estimation network includes a feature extraction backbone network, a static analysis branch, a temporal analysis branch, and a lightweight phase classification branch. The static analysis branch is used to process the initial key point location prediction for a single frame image; The temporal analysis branch is used to process continuous image sequences containing the target frame and its preceding and following frames, and captures the motion correlation of key points in the time dimension through the temporal information processing module and outputs corrected key point position predictions. A lightweight phase classification branch is used to predict gait motion phase.

[0034] Understandably, the dataset is first loaded using dynamic X-ray images with partially labeled bony landmarks and unlabeled dynamic X-ray images, where the labeled portion is J1. A labeled data loader and an unlabeled data loader are constructed. The labeled data loader randomly samples from the labeled data J1, while the unlabeled data loader extracts frames from the folders after all videos have been extracted, following the order of video names. Specifically, it first randomly selects N frames from the image folder of video 1 as a single sample input for the model, and then performs the same operation on the image folders of videos 2, 3, and 4, repeating the process to ensure that the number of samples for the unlabeled data loader is always 2^N. J1.

[0035] Preferably, when training the pose estimation network, data is extracted from two data loaders simultaneously in each batch, with a batch size of Bn. Each batch includes half labeled data and half unlabeled data. To make full use of the unlabeled data, in each iteration of the model, the labeled portion J1 and part of the unlabeled data J1 are loaded first, and then the original labeled data J1 and the remaining unlabeled data are loaded.

[0036] In one embodiment of the invention, the static branch consists of 3×3 convolutions and 1×1 convolutional blocks, outputting a feature heatmap after the one-dimensional convolution. The feature heatmap has dimensions of (batch / 2, number of keypoints, length, width). The temporal branch first uses a 3x3 convolutional network to extract features from frames T-2, T-1, T, T+1, and T+2 sequentially, resulting in a feature map with dimensions of (batch / 2, number of frames, dimension C, 128, 128). Then, through the feature fusion layer of the Context Head, features are further fused using 3×3 convolutions, while the dimension of the feature map remains unchanged. Next, the temporal information is processed by a bidirectional recurrent neural network (ConvGRU). The bidirectional ConvGRU processes the sequence from the forward direction (from t=0 to t=4) and the backward direction (from t=4 to t=0) respectively, and connects the outputs of the two directions at each time step. Assuming that the number of hidden state channels in each ConvGRU layer is dimension C, the number after connection for bidirectional convolution is 2. The C channel has an output dimension of (batch / 2, frame count, 2). (Dimension C, 128, 128). Finally, a 1×1 convolutional layer converts the feature map into a heatmap prediction, with an output dimension of (batch / 2, number of frames, number of key points, 128, 128). Heatmaps of intermediate frames are then extracted from this Context Head, with dimensions of (batch / 2, number of key points, 128, 128). Both methods, using the Static Head and Context Head, decode and map the final keypoint heatmaps to output the predicted keypoint coordinate sets P1 and P2 for frame T.

[0037] Existing techniques often employ single-frame pose estimation methods (such as fully supervised ResNet), which fail to utilize temporal information and perform poorly in motion-blurred and occluded scenes. Traditional Kalman filter post-processing is merely an independent module and cannot be optimized in conjunction with the model. In dynamic X-ray video scenes, keypoints identified in the skeleton are limited by imaging resolution, and due to the high similarity of keypoints, significant pixel shifts may occur during inter-frame recognition, causing errors. This invention explicitly models temporal dependencies using a CovRNN structure in the Context Head, addressing the occlusion problem in dynamic X-ray videos. Furthermore, by combining Kalman filter post-processing, it further enhances the algorithm's robustness to jitter during pose estimation.

[0038] In one embodiment of the present invention, the pose estimation network model is jointly optimized using a composite loss function, wherein the composite loss function includes supervised loss, unsupervised temporal consistency loss, unsupervised pose structure loss, skeleton rigid body constraint loss, and motion phase recognition auxiliary loss, and the calculation formula is as follows: ; in, Represents the composite loss function. This indicates a loss of oversight. The weighting coefficients represent the unsupervised temporal consistency loss. This represents the unsupervised time series consistency loss. The weights represent the unsupervised pose structure loss. This represents the unsupervised attitude structure loss. The weighting coefficients represent the rigid body constraint loss of the skeleton. This indicates the rigid body constraint loss of the skeleton. The weighting coefficients represent the motion phase recognition auxiliary loss. This indicates motion phase recognition auxiliary loss; The supervised loss is obtained by comparing the mean squared error between the prediction results of the static analysis branch on the labeled dataset and the manually labeled ground truth; the formula for calculating the supervised loss is: ; in, This represents the L2 norm (Euclidean distance). This represents the total number of labeled keypoints during the training process. Predict key group for the current frame. The actual label for the current frame; The unsupervised temporal consistency loss calculation is based on the position change of the same key point predicted by the pose estimation network between two temporally adjacent frames. When the position change exceeds the displacement threshold dynamically adjusted according to the motion phase, the part exceeding the threshold is included in the loss term. The unsupervised pose structure loss is based on a pose statistical prior model pre-established from the labeled dataset through principal component analysis. The reconstruction error between the predicted keypoint pose and the pose reconstructed by the pose statistical prior model is calculated and included in the loss term. The rigid body constraint loss of the skeleton is calculated by the deviation between the predicted skeleton segment length and the reference length of the foot and ankle skeleton segment, and the deviation is weighted according to the weight of motion phase modulation. The weighted deviation is included in the loss term. The reference length of the foot and ankle skeleton segment is the median of the lengths of each skeleton segment statistically obtained from the label dataset. The motion phase recognition auxiliary loss lightweight phase classification branch predicts the phase of the current frame in the gait cycle, and calculates the phase by cross-entropy loss with the gait phase pseudo-label automatically generated based on the vertical displacement change of key points. The phase of the gait cycle includes the standing phase, the swinging phase, and the transition phase.

[0039] This invention learns the temporal motion patterns and anatomical constraint features of key points from a large number of unlabeled dynamic X-ray images by designing unsupervised loss terms such as temporal consistency loss, posture structure loss, and skeletal rigid body constraint loss, thereby reducing the number of image annotations and significantly reducing the dependence on annotations.

[0040] Furthermore, the unsupervised temporal consistency loss dynamically adjusts the displacement threshold based on the motion phase, specifically including: The phase probability distribution of the current frame is obtained by predicting the phase of the gait cycle in the current frame based on the key point sequence features extracted by the temporal analysis branch using a lightweight phase classification branch. This phase probability distribution includes the standing phase probability. Oscillation phase probability and transition phase probability ; Set the standing phase weight coefficient Oscillating phase weighting coefficient and transition phase weighting coefficient ,in ; The dynamic displacement threshold is calculated based on the phase probability distribution and weighting coefficients. The position change of the same key point between adjacent frames is calculated. When the position change exceeds the dynamic displacement threshold, the square of the excess is included in the calculation. The temporal consistency loss at this key point in frame number 1; where the dynamic displacement threshold is calculated using the following formula: ; in, Indicates the first The dynamic displacement threshold of the frame. Indicates the reference displacement threshold. Indicates the first The standing phase probability of a frame. Indicates the first The frame swing phase probability, Indicates the first The probability of the transition phase of a frame; The unsupervised temporal consistency loss is obtained by averaging the values ​​of all keypoints and frames within all time windows in the training set. The calculation formula is as follows: The total unsupervised temporal consistency loss is obtained. : ； ； in, Indicates the frame index. Indicates the key point index. Indicates the first Frame number Temporal consistency loss at key points To determine the total number of frames in the time-series window of the training set, This represents the number of key points.

[0041] In one embodiment of the present invention, the unsupervised attitude structure loss calculation includes the following steps: Step 1: Pre-train the PCA model A matrix is ​​constructed using the keypoint coordinates of all labeled frames. Its dimensions are ,in To indicate the total number of labeled frames, The total dimensions of the keypoint coordinates (each keypoint has...) and Two coordinates, A total of 1 key points (Dimension). For the matrix. Perform PCA dimensionality reduction and retain the previous... Principal components (e.g.) Obtain the projection matrix. (dimension is) ) and mean vector (dimension is) ).

[0042] Step 2: Calculate the reconstructed coordinates During training, predict coordinates for each frame. (dimension is) The deviation between the predicted coordinates and the PCA-reconstructed coordinates is calculated. First, the predicted coordinates are projected onto the PCA space to obtain the coordinates in the PCA space. : ; Then reconstruct it back to the original coordinate space: ; Step 3: Calculate the attitude PCA loss Calculate the deviation between the predicted coordinates and the reconstructed coordinates. When the deviation exceeds a tolerance threshold... The unsupervised pose structure loss for frame t is obtained by incorporating the time into the loss. : ; in, Indicates the first Frame number Predicted coordinates of key points Indicates the first Frame number PCA reconstruction coordinates of key points For tolerance threshold (unit: pixels, e.g.) ).

[0043] The unsupervised pose structure loss is obtained by averaging all frames in the training set. : ; Furthermore, the calculation steps for the rigid body constraint loss of the skeleton specifically include: The lengths of skeletal segments in the labeled frames of the labeled dataset are statistically analyzed, and the median length of each skeletal segment is calculated as the reference length for that bone. , where i represents the skeleton index ( (These correspond to the tibia, talus, calcaneus, cuneiform, and metatarsals, respectively); During training, the length of each skeletal segment is calculated based on the keypoint locations predicted by the pose estimation network on each frame of dynamic X-ray image. The skeletal segment is formed by connecting two key points on the same bone. Indicates the frame index. This represents the coordinates of the second keypoint of the i-th bone in the t-th frame of the dynamic X-ray image. This represents the coordinates of the first keypoint of the i-th bone in the t-th frame of the dynamic X-ray image. Calculate the deviation between the predicted bone segment length and the preset reference bone length for each frame. ; The phase modulation weights are calculated based on the motion phase probability distribution, and then weighted using these phase modulation weights to account for the skeletal segment length deviations, resulting in the [number of]th [item / section]. The skeletal rigid body constraint loss of each frame is calculated by averaging the losses across all frames in the training set; the calculation formula is as follows: ; ; in, Indicates the first Frame phase modulation weights, This represents the weighting coefficient of the rigid body constraint in the standing phase. This represents the constraint weight coefficient of the oscillating phase rigid body. This represents the weighting coefficient of the rigid body constraint during the transition phase. Indicates the first The frame's skeletal rigid body constraint loss, where N represents the number of bones. .

[0044] This invention utilizes the physical property that bones maintain a constant length over a short period of time as an unsupervised constraint. By calculating the deviation between the predicted bone segment length and the reference bone length statistically obtained from labeled data, and using phase modulation weights for weighting, it effectively reduces the network's dependence on labeled data, avoids unreasonable situations where the length of the same bone changes drastically between adjacent frames, and improves the model's generalization ability and prediction accuracy on unlabeled data.

[0045] In one embodiment of the present invention, motion phase recognition-assisted loss calculation includes the following steps: Step 1: Generate pseudo tags Gait phase pseudo-labels are automatically generated based on the vertical displacement changes of key points. Specifically, the vertical coordinates of specific key points (such as C1 point of the calcaneus or A2 point of the talus) are selected. Calculate its vertical displacement relative to the initial frame. .

[0046] Automatically divide the phase based on the variation pattern of vertical displacement: when At that time, it is determined to be a standing phase (false label is 0); when When this occurs, it is determined to be a swing phase (pseudo-label is 2); when When this occurs, it is determined to be a transition phase (pseudo-label is 1); in, and The preset displacement threshold can be determined based on statistics from the dataset (e.g.) , ).

[0047] Step 2: Calculate cross-entropy loss The lightweight phase classification branch outputs the phase probability distribution of the current frame. Cross-entropy loss is calculated using pseudo-labels: ; in, For the first Frame pseudo-label category, For indicator functions, For the predicted first The probability of phase-like phase.

[0048] The total motion phase recognition auxiliary loss is obtained by averaging across all frames in the training set. : ; This invention can automatically and synchronously track the fine movements of multiple key bones (such as the tibia, talus, calcaneus, cuneiform, and metatarsals) in the foot-ankle complex, and provides multi-dimensional quantitative analysis methods. This allows us to reveal complex coordinated movement patterns of multiple joints in the foot that are difficult to capture with previous technologies, providing in-depth insights and precise quantitative tools for clinical diagnosis, rehabilitation assessment, and biomechanical research.

[0049] S4. Extract the inference image sequence frame by frame from the dynamic X-ray video to be analyzed, and input the inference image sequence into the trained pose estimation network model to predict frame by frame, so as to obtain the preliminary key point position sequence and the motion phase probability distribution of the current frame. Specifically, frame-by-frame prediction includes: Determine the position of the target frame to be predicted in the inference image sequence. When the target frame is located near the beginning or end of the inference image sequence and a complete temporal analysis window containing 5 frames (2 frames before and 2 frames after) cannot be constructed, only the static analysis branch is used to predict the target frame, and the prediction result of the static analysis branch is used as the preliminary key point position of the target frame. When the target frame can construct a complete temporal analysis window containing two frames before and after it, the target frame is predicted using both static analysis branch and temporal analysis branch to obtain static prediction results and temporal prediction results respectively. The peak intensity of the heatmap corresponding to the static prediction result and the temporal prediction result is calculated as the confidence score, and the two prediction results are weighted and averaged to obtain the preliminary key point position of the target frame. The lightweight phase classification branch outputs the motion phase probability distribution of the target frame while the timing analysis branch outputs the timing prediction results.

[0050] This invention predicts the phase information of gait motion through a lightweight phase classification branch and dynamically adjusts the process noise covariance and measurement noise covariance parameters of the Kalman filter according to the phase probability distribution. It can dynamically balance the weights of model predictions and observations according to changes in motion state, effectively suppressing inter-frame jitter and avoiding the loss of real motion features caused by excessive smoothing, thereby improving the stability and accuracy of key point tracking.

[0051] Preferably, the dynamic X-ray video to be analyzed is used for prediction using a trained attitude estimation network. The dynamic X-ray video to be analyzed is acquired from a dynamic DR system. The attitude estimation network inference mainly includes: cropping the video segment, inputting the cropped 5 consecutive frames into the attitude estimation network, and when time t is less than 3 or greater than the total number of frames - 3, the attitude estimation network only uses the static branch for prediction; otherwise, at time t, the attitude estimation network uses the static branch to predict the key point group of frame t. The keypoint group results from frame t-2 to t+2 are predicted using time branching, and the keypoint group in frame t is selected. By comparison and Confidence mean adaptive selection of current frame key point group results Based on the estimation results of the current frame and historical information, the noise covariance Q and the observation noise covariance R are adjusted, and then predicted using Kalman filtering. It also updates the current observation matrix and state matrix to facilitate key point prediction in subsequent frames.

[0052] S5. Input the preliminary key point position sequence into the Kalman filter, dynamically adjust the process noise covariance and measurement noise covariance parameters of the Kalman filter according to the motion phase probability distribution, and perform time-series smoothing on the preliminary key point position sequence to obtain the filtered key point position sequence. Specifically, the process noise covariance and measurement noise covariance parameters of the Kalman filter are dynamically adjusted based on the motion phase probability distribution, including: The state vector of the Kalman filter is defined to include the position, velocity, and acceleration of the key points, and a constant acceleration motion model is used for state transition. The standard deviation of process noise and the standard deviation of measurement noise are calculated based on the motion phase probability distribution. The calculation formula is as follows: ; ; in, Indicates the first Standard deviation of frame process noise This represents the standard deviation parameter of noise during the standing phase process. This represents the standard deviation parameter of the noise in the oscillating phase process. This parameter represents the standard deviation of noise during the transition phase process. Indicates the first Frame measurement noise standard, This represents the standard deviation parameter of noise measured during standing phase. This represents the standard deviation parameter of the measured noise in the oscillating phase. This represents the standard deviation parameter of the measurement noise during the transition phase. , ; Using the squares of the process noise standard deviation and the measurement noise standard deviation as the diagonal elements of the process noise covariance matrix and the measurement noise covariance matrix, respectively, we obtain the... Frame process noise covariance matrix and measurement noise covariance matrix ; The process noise covariance matrix and measurement noise covariance matrix The prediction and update steps of the Kalman filter are recursively calculated. The prediction step predicts the current state based on the state transition model and the process noise covariance matrix. The update step calculates the optimal estimated state based on the observed values, the measurement noise covariance matrix, and the predicted state. The position component in the optimal estimated state is extracted as the filtered key point position.

[0053] S6. Perform bone rigidity property verification and correction on the filtered key point position sequence, calculate the deviation between the bone segment length and the preset reference bone length in each frame of dynamic X-ray image, and perform geometric correction when the deviation exceeds the tolerance threshold set according to the motion phase to obtain the optimized key point position sequence. Specifically, the steps for testing and correcting the rigid body properties of bones include: Step a: Calculate the length of each skeletal segment in each frame based on the filtered keypoint position sequence. The skeletal segment is formed by connecting two key points on the same bone; Step b: Calculate the length of the bone segment and the preset reference bone length. deviation ; Step c, according to the first The dominant motion phase is determined by the motion phase probability distribution of the frame. The dominant motion phase is the phase with the highest probability value in the motion phase probability distribution. A tolerance threshold is set based on the dominant motion phase of the current frame. The tolerance threshold for standing is less than the tolerance threshold for swinging. Step d, for the first The first frame Root bones, when deviation Exceeding the tolerance threshold If the skeletal segment violates rigid body properties, the direction vector of the skeletal segment from the starting point to the ending point is calculated, and the normalized direction vector is multiplied by the reference bone length. The corrected skeletal segment vector is obtained. Based on the corrected skeletal segment vector, the coordinates of the key points of the two endpoints of the skeletal segment are updated using a partial correction method. The partial correction method is: keeping the position of the midpoint of the skeletal segment unchanged, and symmetrically adjusting the two endpoints to the position that conforms to the reference bone length along the direction vector. Step e: Repeat step ad for all bones in all frames to obtain the optimized key point position sequence after rigid body property verification and correction.

[0054] This invention calculates the deviation between the length of the skeletal segments and the reference bone length in each frame, and sets a tolerance threshold based on the motion phase. It then performs geometric correction on skeletal segments that exceed the threshold. This effectively corrects the violation of rigid body characteristics caused by accumulated errors or abnormal frames, ensuring that the final output key point position sequence satisfies both temporal smoothness and anatomical geometric invariance constraints, thus providing a high-quality data foundation for kinematic parameter calculation.

[0055] S7. Based on the optimized key point position sequence, calculate the rotation angle of each skeletal segment relative to the preset reference system and the relative distance between the geometric center points of different skeletal segments in each frame of dynamic X-ray image to obtain kinematic parameter time series data, and generate time history change analysis charts and synergistic relationship analysis charts based on the kinematic parameter time series data.

[0056] like Figure 4 , Figure 5 and Figure 6 As shown, where α tibia α is the angle of tibial rotation. talus α is the angle of rotation of the talus. calcaneus α is the angle of rotation of the calcaneus. cuneiform α is the cuneiform rotation angle. metatarsal This represents the angle of metatarsal rotation.

[0057] Specifically, step S7 includes: S71. Determine the set of bones and the set of bone pairs for quantitative analysis. The set of bones includes each bone whose rotation angle needs to be calculated, and the set of bone pairs includes each bone pair whose center distance needs to be calculated. Establish a two-dimensional Cartesian coordinate system for each frame of dynamic X-ray image, with the origin located at the lower left corner of the image, the X-axis pointing horizontally to the right, and the Y-axis pointing vertically upward. The rotation angles include the tibial rotation angle, the talus rotation angle, the calcaneus rotation angle, the cuneiform rotation angle, and the metatarsal rotation angle. S72. In the first frame of dynamic X-ray image, the larger Y-coordinate value of the two endpoints of each bone segment in the skeleton set is determined as the origin of rotation of that bone segment. S73. In each frame of dynamic X-ray image, establish a local coordinate system parallel to the two-dimensional Cartesian coordinate system with the origin of rotation, and calculate the angle α between the positive X-axis direction of the local coordinate system and the vector formed by the bone line segment, as the rotation angle of the bone in the current frame. S74. Repeat steps S71-S73 for each bone in each frame of dynamic X-ray image to obtain the set of rotation angles of each bone in each frame of dynamic X-ray image. S75. For each skeletal line segment in each frame of dynamic X-ray image, calculate the arithmetic mean of the coordinates of its two endpoints as the geometric center point of the skeletal line segment. S76. For each bone pair in the set of bone pairs, the Euclidean distance between the geometric center points of the two bone segments in the bone pair is taken as the center distance of the bone pair in the frame. S77. Repeat steps S75-S76 for all selected bone pairs to obtain the set of center distances of each bone pair in each frame of dynamic X-ray image. S78. Organize the set of rotation angles of each bone in each frame of dynamic X-ray image and the set of center distances of each bone pair in each frame of dynamic X-ray image in chronological order to form kinematic parameter time series data. Generate time history change analysis charts and synergy analysis charts based on kinematic parameter time series data.

[0058] This invention constructs a pose estimation network based on semi-supervised learning to automatically identify bony landmarks of the foot-ankle complex in dynamic DR videos. It combines phase-aware adaptive Kalman filtering with posterior verification of skeletal rigidity properties, avoiding manual intervention. Furthermore, by calculating kinematic parameters such as skeletal rotation angles and interosseous center distances, it provides an efficient and accurate quantitative analysis tool for clinical diagnosis, rehabilitation assessment, and biomechanical research.

[0059] like Figure 6 As shown, in one embodiment of the present invention, taking the tibia and medial cuneiform as examples, a graph analyzing the changes in formation time is generated. Wherein, as... Figure 6 As shown in (a), the rotation angle α of a single bone (such as the tibia) calculated for each frame is plotted as a time-series curve with time or frame number as the horizontal axis to visually demonstrate the independent rotation pattern of that bone; such as Figure 6 As shown in (b), select any two target bones and rotate the angle of one bone by a number of consecutive frames (e.g., α). tibia (x) is the X-axis, and the rotation angle change of another bone in consecutive frames (e.g., α) cuneiform A two-dimensional coordinate system is established with the Y-axis as the axis. Data points throughout the entire motion process are plotted on this system, forming a trajectory curve. The shape, slope, and envelope area of ​​this curve can quantitatively reveal the coupling relationship and cooperative pattern between the two selected bones in rotational motion.

[0060] Understandably, the five key bones identified (calcaneus, tibia, talus, metatarsus, and cuneiform) can be paired to generate C(5,2)=10 different relative angle synergistic analysis charts to achieve a comprehensive assessment of ankle joint complex movement.

[0061] like Figure 7 As shown, in one embodiment of the present invention, a synergistic relationship analysis chart is generated using the tibia-calcaneus center distance and the tibia-talus center distance as examples. Figure 7 As shown in (a), the value of the center distance d between any pair of bones calculated in each frame (e.g., d)tibia-calcaneus and d tibia-talus A time-course curve is plotted with time or frame number as the horizontal axis to visually demonstrate the changing patterns of distance between the bones; for example... Figure 7 As shown in (b), select any two interosseous center distances, and vary the length of a continuous frame along one center distance (e.g., d). tibia-calcaneus (x) is the X-axis, and the other center distance is the change in the length of consecutive frames (e.g., d). tibia-talus Establish a two-dimensional coordinate system with the Y-axis as the Y-axis. Plot the data points throughout the entire motion process on this system to form a trajectory curve. This curve can quantify and reveal more complex positional relationships.

[0062] Understandably, the co-analysis method can be applied to any two of the 10 center distances (C(5,2)=10) formed by five key bones. By combining these 10 center distances in pairs, C(10,2)=45 different co-analysis charts of relative positions can be generated to achieve a comprehensive and in-depth analysis of the relative displacement between multiple bones in the ankle joint.

[0063] This invention also provides a dynamic X-ray foot-ankle complex bony landmark analysis device, employing the dynamic X-ray foot-ankle complex bony landmark analysis method as described above, including: The data acquisition module is used to continuously scan the foot and ankle complex of the subject in a preset movement mode through a dynamic DR system to obtain multiple X-ray video samples. The data preprocessing module is used to extract multiple X-ray video samples frame by frame to obtain a dynamic X-ray image sequence of the foot and ankle complex. It also labels multiple pre-set bony landmarks on bones in the dynamic X-ray image sequence to obtain a labeled dataset. The remaining unlabeled dynamic X-ray images are used as an unlabeled dataset. The model training module is used to build the pose estimation network. It mixes the labeled and unlabeled datasets in a certain proportion to form a training set, which is then input into the pose estimation network for semi-supervised training to obtain the trained pose estimation network model. The key point prediction module is used to extract the inference image sequence frame by frame from the dynamic X-ray video to be analyzed, and input the inference image sequence into the trained pose estimation network model to predict frame by frame, so as to obtain the preliminary key point position sequence and the motion phase probability distribution of the current frame. An adaptive filtering module is used to input the preliminary key point position sequence into a Kalman filter, dynamically adjust the process noise covariance and measurement noise covariance parameters of the Kalman filter according to the motion phase probability distribution, and perform time-series smoothing on the preliminary key point position sequence to obtain the filtered key point position sequence. The posterior verification module is used to check and correct the rigid body characteristics of the skeleton in the filtered key point position sequence. It calculates the deviation between the length of the bone line segment and the preset reference bone length in each frame of dynamic X-ray image. When the deviation exceeds the tolerance threshold set according to the motion phase, geometric correction is performed to obtain the optimized key point position sequence. The analysis module is used to calculate the rotation angle of each skeletal segment relative to a preset reference frame and the relative distance between the geometric center points of different skeletal segments in each frame of dynamic X-ray image based on the optimized key point position sequence, to obtain kinematic parameter time series data, and to generate time history change analysis charts and synergy analysis charts based on the kinematic parameter time series data.

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

Claims

1. A method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex, characterized in that: Includes the following steps: S1. The foot and ankle complex of the subject is continuously scanned in a preset motion mode using a dynamic DR system to obtain multiple X-ray video samples. S2. Frame by frame, extract multiple X-ray video samples to obtain dynamic X-ray image sequences of the foot and ankle complex. Label the bony landmarks on multiple preset bones in the dynamic X-ray image sequences to obtain a labeled dataset. Use the remaining unlabeled dynamic X-ray images as an unlabeled dataset. S3. Construct a pose estimation network by mixing the labeled and unlabeled datasets in proportion to form a training set, inputting it into the pose estimation network for semi-supervised training, and obtaining the trained pose estimation network model. S4. Extract the inference image sequence frame by frame from the dynamic X-ray video to be analyzed, and input the inference image sequence into the trained pose estimation network model to predict frame by frame, so as to obtain the preliminary key point position sequence and the motion phase probability distribution of the current frame. S5. Input the preliminary key point position sequence into the Kalman filter, dynamically adjust the process noise covariance and measurement noise covariance parameters of the Kalman filter according to the motion phase probability distribution, and perform time-series smoothing on the preliminary key point position sequence to obtain the filtered key point position sequence. S6. Perform bone rigidity property verification and correction on the filtered key point position sequence, calculate the deviation between the bone segment length and the preset reference bone length in each frame of dynamic X-ray image, and perform geometric correction when the deviation exceeds the tolerance threshold set according to the motion phase to obtain the optimized key point position sequence. S7. Based on the optimized key point position sequence, calculate the rotation angle of each skeletal segment relative to the preset reference system and the relative distance between the geometric center points of different skeletal segments in each frame of dynamic X-ray image to obtain kinematic parameter time series data, and generate time history change analysis charts and synergistic relationship analysis charts based on the kinematic parameter time series data.

2. The method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex as described in claim 1, characterized in that: The pose estimation network includes a feature extraction backbone network, a static analysis branch, a temporal analysis branch, and a lightweight phase classification branch, wherein... The static analysis branch is used to process the initial key point location prediction for a single frame image; The temporal analysis branch is used to process continuous image sequences containing the target frame and its preceding and following frames, and captures the motion correlation of key points in the time dimension through the temporal information processing module and outputs corrected key point position predictions. A lightweight phase classification branch is used to predict gait motion phase.

3. The method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex as described in claim 1, characterized in that: The pose estimation network model is jointly optimized using a composite loss function, which includes supervised loss, unsupervised temporal consistency loss, unsupervised pose structure loss, skeleton rigid body constraint loss, and motion phase recognition auxiliary loss. The supervised loss is obtained by comparing the mean squared error between the prediction results of the static analysis branch on the labeled dataset and the manually labeled true values; The unsupervised temporal consistency loss calculation is based on the position change of the same key point predicted by the pose estimation network between two temporally adjacent frames. When the position change exceeds the displacement threshold dynamically adjusted according to the motion phase, the part exceeding the threshold is included in the loss term. The unsupervised pose structure loss is based on a pose statistical prior model pre-established from the labeled dataset through principal component analysis. The reconstruction error between the predicted keypoint pose and the pose reconstructed by the pose statistical prior model is calculated and included in the loss term. The rigid body constraint loss of the skeleton is calculated by the deviation between the predicted skeleton segment length and the reference length of the foot and ankle skeleton segment, and the deviation is weighted according to the weight of motion phase modulation. The weighted deviation is included in the loss term. The reference length of the foot and ankle skeleton segment is the median of the lengths of each skeleton segment statistically obtained from the label dataset. The motion phase recognition auxiliary loss lightweight phase classification branch predicts the phase of the current frame in the gait cycle, and calculates the phase by cross-entropy loss with the gait phase pseudo-label automatically generated based on the vertical displacement change of key points. The phase of the gait cycle includes the standing phase, the swinging phase, and the transition phase.

4. The method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex as described in claim 3, characterized in that: The unsupervised temporal consistency loss dynamically adjusts the displacement threshold based on the motion phase, specifically including: The phase probability distribution of the current frame is obtained by predicting the phase of the gait cycle in the current frame based on the key point sequence features extracted by the temporal analysis branch using a lightweight phase classification branch. This phase probability distribution includes the standing phase probability. Oscillation phase probability and transition phase probability ; Set the standing phase weight coefficient Oscillating phase weighting coefficient and transition phase weighting coefficient ,in ; The dynamic displacement threshold is calculated based on the phase probability distribution and weighting coefficients. The position change of the same key point between adjacent frames is calculated. When the position change exceeds the dynamic displacement threshold, the square of the excess is included in the calculation. The temporal consistency loss at this key point in frame number 1; where the dynamic displacement threshold is calculated using the following formula: ； in, Indicates the first The dynamic displacement threshold of the frame. Indicates the reference displacement threshold. Indicates the first The standing phase probability of a frame. Indicates the first The frame swing phase probability, Indicates the first The probability of the transition phase of a frame; The unsupervised temporal consistency loss is obtained by averaging the frames within all time windows of all keypoints and training sets.

5. The method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex as described in claim 1, characterized in that: The calculation steps for the rigid body constraint loss of the skeleton specifically include: The lengths of skeletal segments in the labeled frames of the labeled dataset are statistically analyzed, and the median length of each skeletal segment is calculated as the reference length for that bone. , where i represents the skeleton index; During training, the length of each skeletal segment is calculated based on the keypoint locations predicted by the pose estimation network on each frame of dynamic X-ray image. The skeletal segment is formed by connecting two key points on the same bone. Indicates the frame index; Calculate the deviation between the predicted bone segment length and the preset reference bone length for each frame; The phase modulation weights are calculated based on the motion phase probability distribution, and then weighted using these phase modulation weights to account for the skeletal segment length deviations, resulting in the [number of]th [item / section]. The skeletal rigid body constraint loss of each frame is calculated by averaging the losses across all frames in the training set; the calculation formula is as follows: ； ； in, Indicates the first Frame phase modulation weights, This represents the weighting coefficient of the rigid body constraint in the standing phase. This represents the constraint weight coefficient of the oscillating phase rigid body. This represents the weighting coefficient of the rigid body constraint during the transition phase. Indicates the first The frame's skeletal rigid body constraint loss, where N represents the number of bones. .

6. The method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex as described in claim 5, characterized in that: The frame-by-frame prediction specifically includes: Determine the position of the target frame to be predicted in the inference image sequence. When the target frame is located near the beginning or end of the inference image sequence and a complete temporal analysis window containing 5 frames (2 frames before and 2 frames after) cannot be constructed, only the static analysis branch is used to predict the target frame, and the prediction result of the static analysis branch is used as the preliminary key point position of the target frame. When the target frame can construct a complete temporal analysis window containing two frames before and after it, the target frame is predicted using both static analysis branch and temporal analysis branch to obtain static prediction results and temporal prediction results respectively. The peak intensity of the heatmap corresponding to the static prediction result and the temporal prediction result is calculated as the confidence score, and the two prediction results are weighted and averaged to obtain the preliminary key point position of the target frame. The lightweight phase classification branch outputs the motion phase probability distribution of the target frame while the timing analysis branch outputs the timing prediction results.

7. The method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex as described in claim 6, characterized in that: The process noise covariance and measurement noise covariance parameters of the Kalman filter, which are dynamically adjusted according to the motion phase probability distribution, specifically include: The state vector of the Kalman filter is defined to include the position, velocity, and acceleration of the key points, and a constant acceleration motion model is used for state transition. The standard deviation of process noise and the standard deviation of measurement noise are calculated based on the motion phase probability distribution. The calculation formula is as follows: ； ； in, Indicates the first Standard deviation of frame process noise This represents the standard deviation parameter of noise during the standing phase process. This represents the standard deviation parameter of the noise in the oscillating phase process. This parameter represents the standard deviation of noise in the transition phase process. Indicates the first Frame measurement noise standard, This represents the standard deviation parameter of noise measured in the standing phase. This represents the standard deviation parameter of the measured noise in the oscillating phase. This represents the standard deviation parameter of the measurement noise during the transition phase. , ; Using the squares of the process noise standard deviation and the measurement noise standard deviation as the diagonal elements of the process noise covariance matrix and the measurement noise covariance matrix, respectively, we obtain the... Frame process noise covariance matrix and measurement noise covariance matrix ; The process noise covariance matrix and measurement noise covariance matrix The prediction and update steps of the Kalman filter are recursively calculated. The prediction step predicts the current state based on the state transition model and the process noise covariance matrix. The update step calculates the optimal estimated state based on the observed values, the measurement noise covariance matrix, and the predicted state. The position component in the optimal estimated state is extracted as the filtered key point position.

8. The method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex as described in claim 7, characterized in that: The specific steps for verifying and correcting the rigid body properties of the skeleton include: Step a: Calculate the length of each skeletal segment in each frame based on the filtered keypoint position sequence. The skeletal segment is formed by connecting two key points on the same bone; Step b: Calculate the length of the bone segment and the preset reference bone length. deviation ; Step c, according to the first The dominant motion phase is determined by the motion phase probability distribution of the frame. The dominant motion phase is the phase with the highest probability value in the motion phase probability distribution. A tolerance threshold is set based on the dominant motion phase of the current frame. The tolerance threshold for standing is less than the tolerance threshold for swinging. Step d, for the first The first frame Root bones, when deviation Exceeding the tolerance threshold If the skeletal segment violates rigid body properties, the direction vector of the skeletal segment from the starting point to the ending point is calculated, and the normalized direction vector is multiplied by the reference bone length. The corrected skeletal segment vector is obtained. Based on the corrected skeletal segment vector, the coordinates of the key points of the two endpoints of the skeletal segment are updated using a partial correction method. The partial correction method is: keeping the position of the midpoint of the skeletal segment unchanged, and symmetrically adjusting the two endpoints to the position that conforms to the reference bone length along the direction vector. Step e: Repeat step ad for all bones in all frames to obtain the optimized key point position sequence after rigid body property verification and correction.

9. The method for dynamic X-ray analysis of bony landmarks in the foot-ankle complex as described in claim 1, characterized in that: Step S7 specifically includes: S71. Determine the set of bones and the set of bone pairs for quantitative analysis. The set of bones includes each bone whose rotation angle needs to be calculated. The set of bone pairs includes each bone pair whose center distance needs to be calculated. Establish a two-dimensional Cartesian coordinate system for each frame of dynamic X-ray image. The origin is located at the lower left corner of the image. The X-axis is horizontal to the right and the Y-axis is vertical to the up. S72. In the first frame of dynamic X-ray image, the larger Y-coordinate value of the two endpoints of each bone segment in the skeleton set is determined as the origin of rotation of that bone segment. S73. In each frame of dynamic X-ray image, establish a local coordinate system parallel to the two-dimensional Cartesian coordinate system with the origin of rotation, and calculate the angle between the positive X-axis direction of the local coordinate system and the vector formed by the bone line segment, which is used as the rotation angle of the bone in the current frame. S74. Repeat steps S71-S73 for each bone in each frame of dynamic X-ray image to obtain the set of rotation angles of each bone in each frame of dynamic X-ray image. S75. For each skeletal line segment in each frame of dynamic X-ray image, calculate the arithmetic mean of the coordinates of its two endpoints as the geometric center point of the skeletal line segment. S76. For each bone pair in the set of bone pairs, the Euclidean distance between the geometric center points of the two bone segments in the bone pair is taken as the center distance of the bone pair in the frame. S77. Repeat steps S75-S76 for all selected bone pairs to obtain the set of center distances of each bone pair in each frame of dynamic X-ray image. S78. Organize the set of rotation angles of each bone in each frame of dynamic X-ray image and the set of center distances of each bone pair in each frame of dynamic X-ray image in chronological order to form kinematic parameter time series data. Generate time history change analysis charts and synergy analysis charts based on kinematic parameter time series data.

10. A dynamic X-ray foot-ankle complex bony landmark analysis device, characterized in that: The method for analyzing bony landmarks of the foot-ankle complex using dynamic X-ray as described in any one of claims 1-9 includes: The data acquisition module is used to continuously scan the foot and ankle complex of the subject in a preset movement mode through a dynamic DR system to obtain multiple X-ray video samples. The data preprocessing module is used to extract multiple X-ray video samples frame by frame to obtain a dynamic X-ray image sequence of the foot and ankle complex. It also labels multiple pre-set bony landmarks on bones in the dynamic X-ray image sequence to obtain a labeled dataset. The remaining unlabeled dynamic X-ray images are used as an unlabeled dataset. The model training module is used to build the pose estimation network. It mixes the labeled and unlabeled datasets in a certain proportion to form a training set, which is then input into the pose estimation network for semi-supervised training to obtain the trained pose estimation network model. The key point prediction module is used to extract the inference image sequence frame by frame from the dynamic X-ray video to be analyzed, and input the inference image sequence into the trained pose estimation network model to predict frame by frame, so as to obtain the preliminary key point position sequence and the motion phase probability distribution of the current frame. An adaptive filtering module is used to input the preliminary key point position sequence into a Kalman filter, dynamically adjust the process noise covariance and measurement noise covariance parameters of the Kalman filter according to the motion phase probability distribution, and perform time-series smoothing on the preliminary key point position sequence to obtain the filtered key point position sequence. The posterior verification module is used to check and correct the rigid body characteristics of the skeleton in the filtered key point position sequence. It calculates the deviation between the length of the bone line segment and the preset reference bone length in each frame of dynamic X-ray image. When the deviation exceeds the tolerance threshold set according to the motion phase, geometric correction is performed to obtain the optimized key point position sequence. The analysis module is used to calculate the rotation angle of each skeletal segment relative to a preset reference frame and the relative distance between the geometric center points of different skeletal segments in each frame of dynamic X-ray image based on the optimized key point position sequence, to obtain kinematic parameter time series data, and to generate time history change analysis charts and synergy analysis charts based on the kinematic parameter time series data.

Images