Bayesian-based VR table tennis motion track backstepping method

Abstract

The invention relates to a Bayesian-based VR table tennis motion track backstepping method. The method comprises the following steps: 1, acquiring original VR motion data; 2, preprocessing the VR motion data; step 3, constructing physiological priori knowledge; 4, constructing angle constraint priori knowledge; 5, solving forward motion; step 6, Bayesian solution is carried out; and step 7, visualizing the motion trail. According to the method, state space modeling is taken as a core, human body biomechanical prior and high-order Markov time sequence constraint are fused, adaptive noise observation and variational Bayesian reasoning are combined, and real-time updating and uncertainty quantification are realized in parallel through a GPU (Graphics Processing Unit).

Description

Technical Field

[0001] This invention relates to a Bayesian-based method for inferring the trajectory of a VR table tennis ball, which is used to achieve high-precision estimation of human kinematic parameters, prediction of table tennis ball trajectory, and assessment of motor skills, and belongs to the field of virtual reality motion analysis technology. Background Technology

[0002] In human-computer interaction scenarios of virtual reality (VR) and its extension to augmented / mixed reality, human posture and motion estimation are core technologies supporting immersive experiences, motion interaction, safety boundaries, and ergonomic assessment. By fusing data from multiple sources of sensors, including the head-mounted display (HMD) and controllers' inertial measurement units (IMUs), wearable motion sensors, inward / outward vision cameras, and laser positioning / depth cameras, the system can reconstruct the joint angles, velocities, and contact states of the human skeleton, thereby driving virtual body shapes, assessing spatial risks, and providing real-time feedback. Compared to traditional offline motion capture, VR scenarios emphasize "real-time, robust, and individualized" capabilities under device constraints, thus placing higher demands on algorithm noise suppression, temporal consistency, and computational efficiency. Traditionally, VR pose estimation has been carried out along three main paths: (1) inverse kinematics (IK) and constrained optimization based on kinematics, inverting upper limb / whole body pose from a small number of controllers and head-mounted display poses; (2) filtering methods based on inertial fusion, such as complementary filtering and extended / unscented Kalman filtering, which take into account both short-term smoothing and long-term drift control; and (3) time-series regression (RNN / TCN / Transformer) or graph networks based on deep learning, estimating three-dimensional joint angles from IMU sequences or videos. Although these methods have made progress under specific conditions, they all face the common challenges of sensitivity to sensor noise and occlusion, insufficient characterization of nonlinear and non-Gaussian errors, poor generalization to individual differences and scene changes, and the difficulty in achieving high accuracy and low latency under limited computing power. In addition, most methods only provide point estimates and lack uncertainty quantification that can be used for system decision-making and safety control.

[0003] Specific challenges in existing technologies include: (1) Observation noise and drift: IMUs suffer from zero-bias drift and magnetic interference. Visual solutions are susceptible to occlusion, field of view and illumination. External positioning may experience momentary loss or jumps. These unstable factors are particularly prominent in rapid movement and large-scale walking, causing posture errors to accumulate over time and making timely self-consistent correction difficult. (2) Multi-source asynchronous and calibration challenges: Different devices have inconsistent sampling rates, timestamp references and communication delays. External parameter calibration (sensor to skeleton, geometry between devices) drifts with changes in wearing and clothing. The lack of online self-calibration and time synchronization mechanisms will introduce systematic biases, affecting fusion accuracy and stability. (3) Insufficient model priors and observability: When relying on only a small number of external pose observations, the underdeterminacy of human joints and posture ambiguity are significant. If biomechanical priors (joint amplitude limiting, limb length consistency, ground contact and dynamic feasibility) are not explicitly introduced, the algorithm is prone to skeletal stretching, cross-intersection and non-physical understanding. (4) Insufficient utilization of temporal and dynamic constraints: Estimating attitude independently frame by frame or using only low-order smoothing makes it difficult to capture non-stationary dynamics such as sudden acceleration, impact, and contact transition, resulting in increased jitter and latency during violent movements, rapid starts and stops, and stamping / grasping events. (5) Computational budget and real-time performance: Complex optimizations or deep networks are difficult to run stably under the constraints of computing power, power consumption, and heat dissipation on the head-mounted display side; cloud inference is also subject to network jitter and privacy restrictions, forming a multi-objective trade-off between accuracy, latency, and stability. (6) Lack of robustness and uncertainty measurement: Most systems do not model for thick-tailed noise, outliers, and missing observations, lacking credible confidence intervals or risk indicators, making it difficult to support safety boundary management, abnormal action alarms, and adaptive control strategies.

[0004] The aforementioned issues collectively limit the upper limits of existing VR pose estimation algorithms in terms of accuracy, efficiency, and adaptability, especially revealing significant shortcomings in fast-paced interactions and long-term usage scenarios. Therefore, there is an urgent need for a Bayesian estimation framework that can integrate biomechanical priors and multi-source observations, is temporally self-consistent, and provides uncertainty quantification, in order to achieve high-precision, low-latency, robust, and interpretable pose and motion estimation under limited computing power. Summary of the Invention

[0005] To overcome the shortcomings of existing pose / motion estimation methods in VR scenarios, such as noise sensitivity, temporal inconsistency, lack of uncertainty measurement, and lack of real-time performance, this invention provides a Bayesian-based method for inferring the trajectory of a VR table tennis ball. This method takes state-space modeling as its core, integrates human biomechanical priors and higher-order Markov temporal constraints, combines adaptive observation noise and variational Bayesian inference, and achieves real-time updates and uncertainty quantification through GPU parallel processing.

[0006] A Bayesian-based method for reverse-engineering the trajectory of a VR table tennis ball includes the following steps:

[0007] Step 1: Acquisition of Raw VR Motion Data: Collect spatial position information, rotational posture information, and event information of the VR headset and controllers to form a data sequence organized by frame; collect IMU data from the VR headset and both controllers to achieve observation coverage of head and hand movements; combine unified timestamp records and device status metadata to provide basic support for subsequent multimodal data fusion and motion reconstruction.

[0008] Step Two: VR Motion Data Preprocessing: The acquired VR data sequence is preprocessed, including time base unification, data framing, coordinate system integration, and extrinsic parameter calibration. Resampling is performed at the target frame rate, and short-term dropped frames are filled in. Quaternion poses are normalized. Anomaly detection, interpolation repair, zero-bias estimation, and drift compensation are applied, outputting frame-by-frame confidence labels. Head position / rotation and eye-tracking channels undergo IQR anomaly detection followed by linear interpolation and Savitzky-Golay smoothing, maintaining synchronization with frame parsing. Self-calibration is then implemented in the Bayesian solver: Limb length priors are provided based on user height / gender; sliding window estimation is used, and noise covariance is fused with the priors, while compensating for time-dependent drift / random walks and environmental interference; finally, the angle sequence is temporally smoothed and biomechanical constraints are applied to obtain robust estimates.

[0009] Step 3: Construction of physiological prior knowledge: Population statistical information and physiological / kinematic knowledge are continuously encoded into "length prior + kinematic prior" using a unified Bayesian expression. Based on the user's height and gender, individualized mean and standard deviation of limb length are generated, proportional and symmetry constraints are established, and the boundaries of the joint physiological mobility domain are set to form computable parameter constraints.

[0010] Step 4: Constructing Prior Knowledge of Angle Constraints: Construct a state vector containing angle, angular velocity, angular acceleration and jerk, introduce first-order, second-order and third-order higher-order Markov constraints, integrate soft physiological mobility domain constraints and temporal smoothing terms, and output window-level scores and penalty intensity through a "stationary-trend-correction" adaptive mode.

[0011] Step 5: Forward Motion Solution: Joint angles and limb lengths are mapped forward geometrically to the position and posture of the hand in the world coordinate system, and further transformed to the head-mounted display coordinate system for analytical predictions compared with sensor readings. Aiming at physiological constraints and observational consistency of human motion, a logarithmic posterior model of head and hand postures is constructed, and joint parameters for each frame are solved using MAP (maximum a posteriori estimation). The model uses the distribution of human limb lengths corresponding to height / gender as the length prior, the continuity and smoothness of angles as the temporal prior, and combines the likelihood function of multi-source IMU observations to form the optimization objective. In numerical implementation, L-BFGS-B optimization is used, and adaptive covariance and biomechanical penalty terms are used to ensure the physical rationality of the solution at joints such as the elbow and shoulder. The algorithm can be solved in parallel across multiple frames under GPU acceleration, and the sliding window history is used to achieve temporal consistency of angular velocity and angular acceleration, ultimately yielding robust angle estimates and trajectory outputs.

[0012] Step Six: Bayesian Solution: Unify the forward motion prediction from Step Five with the priors given in Steps Three and Four into the same probabilistic framework. Construct an observation consistency metric based on multi-channel noise statistics. Achieve maximum a posteriori (MAP) estimation through numerical optimization with physiological and velocimetric boundaries. Combine sliding window iteration and convergence criteria to output parameter estimation results and uncertainty quantification information. The sliding window MAP solution maximizes temporal smoothness and observation consistency within the window and suppresses transient noise and spike interference through an adaptive weighting mechanism. Outliers are repaired using interpolation based on a posteriori consistency. Biomechanical boundary protection and backoff strategies are applied to the smoothing results to ensure that the angle sequence is continuous in time and reasonable in space. Finally, the smoothed angle, angular velocity, and their higher-order derivative sequences, as well as statistical indicators such as continuity score and constraint hit rate, are output, providing high-confidence input for subsequent swing feature analysis and visualization.

[0013] Step Seven: Motion Trajectory Visualization: Based on the joint angles and poses obtained in Step Six, design and implement dynamic visualization of the motion trajectory, presenting the single-round and multi-round movement process in a readable and interactive manner. Provide a frame-level angle reading and aggregation mechanism; generate detailed and simplified angle data, covering shoulder triaxial alignment, racket orientation, swing speed, and related time information; aggregate inverse kinematics statistics (such as convergence, constraint hit rate, anomaly removal ratio, etc.) and key indicators; unify integration during the swing analysis stage, supporting feature extraction, comparison, and report generation, and supporting multi-round batch processing and visualization. This step is not limited to a specific implementation method and can be integrated with different data storage and front-end presentation solutions.

[0014] Step two specifically includes: first, unifying the time base and data framing, synchronizing the timestamps of the head-mounted display and the two-handed controllers to a monotonic clock, resampling according to the target frame rate and filling in short-term dropped frames, and simultaneously completing coordinate system integration and extrinsic parameter calibration, unifying the poses of each device to the world or head-mounted display coordinate system and normalizing the quaternion poses; then, performing quality assessment and anomaly handling, using interquartile range and three times standard deviation rules to detect position and pose anomalies, performing linear or spline interpolation on isolated anomalies and missing segments, performing nearest neighbor interpolation and continuity for pose jumps, estimating zero bias in static segments, compensating for slow drift and random walks within a sliding window, and simultaneously outputting frame-by-frame confidence markers.

[0015] Step three specifically includes:

[0016] Using user height and gender as independent variables, a priori limb length is established, and physiological and velocimetric boundaries are simultaneously determined. Based on the population statistical baseline, the mean lengths of the neck, upper arm, forearm, and shoulder width are adaptively scaled according to height. The mean lengths of each part are linearly adjusted according to the deviation of the user's height from the reference height. Then, a gender correction is applied, with a gender coefficient set for different parts. The scaling factor is 1 for males and 0.9 to 0.95 for females to obtain individualized mean lengths. Then, the standard deviation is adaptively adjusted so that the standard deviation of length is linearly amplified with the degree of height deviation. When the height is close to the reference value, the standard deviation of length is adjusted accordingly. Standard deviation contraction increases as height deviation increases, reflecting individual differences and avoiding over-contraction of uncommon body types. Proportion and symmetry constraints are established, limiting the reasonable proportion range of the upper arm and forearm, setting an upper limit for the symmetry deviation of the left and right upper limbs, applying penalty weights to combinations exceeding limits to ensure geometric consistency, providing physiological mobility domain boundaries, setting upper and lower limits for the shoulder triaxial axis and elbow flexion and extension respectively, and making differentiated configurations according to the racket side and non-racket side to match the characteristics of the project, forming prior outputs, including the length prior distribution and covariance table, physiological and velocity boundary tables, and proportion and symmetry constraint coefficients.

[0017] Step four specifically includes:

[0018] A state vector containing angle, angular velocity, angular acceleration, and jerkiness is constructed. Using historical information from multiple adjacent frames as the temporal basis, higher-order Markov constraints of first-, second-, and third-order historical terms are introduced. Combined with soft physiological mobility domain constraints and temporal smoothing terms, sudden spikes and out-of-bounds behaviors are penalized, and window-level scores and penalty intensity are output. The above temporal priors are adaptively fused using three modes: "stationary," "trend," and "correction," to maintain continuity and interpretability under different states such as high-speed swings and static transitions.

[0019] Step six specifically includes:

[0020] The observations of position and attitude are compared with the forward kinematics predictions, and an observation consistency metric is constructed based on multi-channel noise statistics. The length prior and temporal prior obtained in steps three and four are jointly modeled within the same posterior framework. Numerical optimization with physiological and velocity boundaries is used to obtain the maximum posterior estimate of joint angle and limb length. Preferably, the solution of the previous window is used as the initial value in the form of a sliding window, and the target convergence and gradient threshold are used as stopping criteria. The curvature is approximated at the optimum to obtain the parameter uncertainty, which is used to generate confidence intervals and subsequent visualization confidence bands.

[0021] Step seven specifically includes:

[0022] Based on the joint angle and length parameters obtained through reverse engineering, a rigid body chain model with the shoulder joint as the root node is established. The three-axis rotation of the shoulder, the geometric connection between the upper arm and forearm, and the installation of extrinsic parameters at the end are sequentially connected to obtain the spatial pose of the hand. The "length effect" and "angle effect" are decoupled structurally to facilitate sensitivity analysis, extrinsic parameter correction, and real-time calculation, and it is naturally adapted to the analytical Jacobian to accelerate subsequent solutions. On this basis, a multi-view visualization is generated, rendering the end trajectory and velocity curve, and using color or markers to link key events and stage intervals.

[0023] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0024] The core innovation of this invention lies in constructing a "unified kinematic inference framework integrating physiological priors" and a "real-time estimation method for multi-scale Bayesian parameters," breaking through the limitations of traditional motion capture.

[0025] In terms of physiological prior fusion, based on individualized anatomical parameter modeling, higher-order Markov temporal constraints, and Bayesian posterior probability optimization, synchronous and accurate estimation of joint angles, racket spatial posture, and swing speed is achieved. Combining Gaussian process regression and variational Bayesian inference, intrinsic control patterns of motor skills are extracted, revealing the laws governing hitting strategies and skill formation, and providing a quantitative tool for research on the neuro-motor mechanisms of motor skill acquisition.

[0026] In real-time estimation of multi-scale Bayesian parameters, a dual-sine oscillation model and adaptive sampling strategy are used to achieve visualization of dynamic parameter evolution with millisecond-level time resolution. Real-time display of instantaneous changes in key parameters such as shot cycle, reaction delay, and amplitude enhances the intuitiveness and relevance of data analysis, bridging sports science research and practical training applications.

[0027] This invention overcomes the limitations of traditional motion capture analysis in two aspects: a unified kinematic inference framework integrating physiological priors and a real-time multi-scale Bayesian parameter estimation method, providing a new technical path for human movement science and sports training applications. Firstly, regarding the unified kinematic inference framework integrating physiological priors, this invention achieves simultaneous and accurate estimation of human joint angles, racket spatial posture, and swing speed based on individualized anatomical parameter modeling, higher-order Markov temporal constraints, and Bayesian posterior probability optimization, rather than relying solely on traditional single-point measurements or independent fitting methods. This method not only characterizes the evolution of dynamic features at different stages of movement but also combines Gaussian process regression and variational Bayesian inference to extract intrinsic control patterns of motor skills, thereby revealing deeper-level hitting strategies and skill formation patterns. This analytical framework surpasses traditional static evaluation methods based on experience thresholds, dynamically tracking athletes' technical changes at different training levels, providing a more refined quantitative tool for studying the neuro-motor mechanisms of motor skill acquisition.

[0028] In the area of ​​real-time estimation methods for multi-scale Bayesian parameters, this invention overcomes the limitations of existing methods that can only present post-hoc statistical analysis, proposing a dynamic parameter evolution visualization scheme with millisecond-level temporal resolution. By constructing a dual-sine oscillation model and an adaptive sampling strategy, it can display the instantaneous changes of key parameters such as hitting cycle, reaction delay, and amplitude in real time, thus intuitively reflecting the real-time adjustment of motion control strategies and their evolution over time. This method not only enhances the expressiveness of motion data analysis results, enabling coaches to more intuitively perceive changes in athletes' technical state, but also makes it easier for athletes to understand their own movement characteristics, improving the interpretability and relevance of sports training.

[0029] This invention achieves breakthroughs in improving the accuracy, stability, and real-time visualization capabilities of motion parameter estimation, providing a new methodological foundation for applications in sports science training, skills assessment, injury prevention, and human-computer interaction. This makes motion skill analysis not only more scientifically rigorous but also possesses broader potential for educational applications. By transforming complex biomechanical modeling into an intuitive parameter evolution process, this invention bridges the gap between sports science research and practical training applications. Attached Figure Description

[0030] 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.

[0031] Figure 1This is a flowchart of a Bayesian-based VR table tennis trajectory inversion method according to the present invention.

[0032] Figure 2 This is a flowchart of the Bayesian inference framework for a Bayesian-based VR table tennis trajectory reverse inference method of the present invention.

[0033] Figure 3 This is a comparison chart of the algorithm stability of the Bayesian-based VR table tennis trajectory inversion method of the present invention.

[0034] Figure 4 This is a schematic diagram of the joint angles in a single round of a Bayesian-based VR table tennis trajectory inference method according to the present invention.

[0035] Figure 5 This is a schematic diagram of motion trajectory visualization for a Bayesian-based VR table tennis motion trajectory inversion method according to the present invention. Detailed Implementation

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

[0037] Reference Figure 1 A Bayesian-based method for reverse-engineering the trajectory of a VR table tennis ball includes the following steps:

[0038] Step 1: Acquisition of raw VR motion data: Collect spatial position information, rotation posture information and event information of VR headset and controller, and form a data sequence organized by frame;

[0039] Specifically, the raw VR data received by the system includes multimodal sensor information, covering 6DoF data from the headset and two controllers. Its main features include spatial position information, rotational attitude information, and a timestamp system. Spatial position information includes the 3D coordinates of the head, right hand, and left hand. Rotational attitude information includes head quaternion rotation and hand rotation matrices. The timestamp system uses absolute timestamps to ensure data consistency across sessions, and frame sequence numbers to guarantee data temporal integrity, enabling time / coordinate alignment, self-calibration, and subsequent Bayesian kinematics inference and temporal smoothing during preprocessing.

[0040] Step 2: VR motion data preprocessing: The acquired VR data sequence is preprocessed, including time base unification, data framing, coordinate system integration, and extrinsic parameter calibration. The data sequence is resampled according to the target frame rate and short-term lost frames are filled in. Quaternion poses are normalized. Through anomaly detection, interpolation repair, zero bias estimation, and drift compensation, frame-by-frame confidence labels are output.

[0041] Specifically, outlier detection of the collected VR motion data employs an interquartile range (ICM)-based method. For each spatial dimension, the first and third quartiles are calculated to obtain the ICM. The outlier determination interval is determined based on 1.5 times the ICM. Furthermore, the system employs a physiological constraint-based rationality check, including head movement speed limits and hand reach constraints, to ensure the data conforms to the physiological characteristics of human movement. Data smoothing then utilizes the Savitzky-Golay filtering method. The window size is selected based on one-tenth of the sampling frequency, and the polynomial order is set to 3, effectively suppressing high-frequency noise while preserving signal characteristics. The smoothing effect is evaluated using the variance ratio before and after smoothing and signal fidelity, ensuring that key features of the original signal are preserved while removing noise.

[0042] Outlier range: [Q1 - 1.5 * IQR, Q3 + 1.5 * IQR]

[0043] Where: Q1 is the first quartile, Q3 is the third quartile, and IQR = Q3 - Q1 is the interquartile range.

[0044] Data segmentation is based on match events. The system segments continuous data into meaningful motion segments according to key time points in the mission data, including the start time of the serve, the moment of hitting the ball, and the end time of the rally. Each rally is further divided into six key phases: preparation phase, serve phase, receive preparation phase, hitting phase, return phase, and end phase, providing a structured data foundation for subsequent analysis.

[0045] Step 3: Construction of physiological prior knowledge: Population statistical information and physiological / kinematic knowledge are continuously encoded into "length prior + kinematic prior" using a unified Bayesian expression. Based on the user's height and gender, individualized mean and standard deviation of limb length are generated, proportional and symmetry constraints are established, and the boundaries of the joint physiological mobility domain are set to form computable parameter constraints.

[0046] Specifically,

[0047] Starting with the length domain: For parameters such as neck length, upper and lower arm length, and shoulder width, a linear correction is first applied to the baseline mean using a height-adaptive method, resulting in the following calculation formula:

[0048]

[0049] Where μ base ρ represents the mean baseline length of the population at this location. H Here is the correlation coefficient between body part length and height, where H is the user's height. ref For reference height, the deviation of an individual's height from the population reference is converted into a proportional correction to the mean; the higher the height, the greater the increase in μ.

[0050] Sex correction: Sex correction is used to incorporate male-female body size differences into the individualized estimate of length mean. Using males as the baseline, r is taken as... gender (male) = 1.00; females are given a scaling factor of 0.90–0.95 based on different body parts, and can be recalibrated based on statistical data. The calculation formula is as follows:

[0051] μ final =μ adj ·r gender

[0052] Where μ adj Let r be the mean of the height after adaptive calculation above. gender This is a gender correction factor, configured according to body part, used to map the differences in body shape between men and women onto the mean. It is usually set to 1.00 for men and 0.90 to 0.95 for women.

[0053] The adaptive standard deviation is used to adjust the uncertainty of the prior length according to the degree of deviation from height. The calculation formula is as follows:

[0054] σ obj =σ base ·(0.8+0.4·|HH ref |)

[0055] Where, σ base H represents the baseline standard deviation of this body part within the population, where H is the user's height. ref For reference height. The coefficient 0.8 is defined as when H = H ref Minimum contraction (σ) obj =0.8σ base The coefficient 0.4 determines the degree of deviation from |HH. ref | / H ref The slope is increased linearly. Therefore, the prior is more concentrated when the height is close to the reference value, and more dispersed when the deviation is greater, thus reflecting individual differences and avoiding excessive constraints on uncommon body types.

[0056] After completing height adaptation and gender correction, the individualized mean μ was obtained. final And obtain the individualized standard deviation σ based on the height deviation. obj Then, the length L of any limb segment can be represented as a conditional normal distribution, and the calculation formula is as follows:

[0057]

[0058] Where L is a random variable representing the length of this part, μ final With σ obj These are the individualized mean and standard deviation obtained in the previous two steps, respectively. The significance of this step is to apply macroscopic information such as height and gender to the probability density of the length parameter, which is used to generate confidence intervals and sampling, and serves as the prior input for length in subsequent solutions, thereby maintaining consistency with population statistics while considering individual differences.

[0059] Step 4: Constructing Angle-Constrained Prior Knowledge: Construct a state vector containing angle, angular velocity, angular acceleration, and jerk. Introduce first-, second-, and third-order higher-order Markov constraints, integrate soft physiological mobility domain constraints and temporal smoothing terms, and output window-level scores and penalty intensity through a "stationary-trend-correction" adaptive mode. The physiological mobility domain of the joint and the continuity in time are jointly encoded as angle priors, so that the angle sequence satisfies the anatomical range and is smooth and predictable in time.

[0060] Therefore, the angle sequence is first expressed as a measurable derivative, which serves as the basis for subsequent constraints and priors. The calculation formula is as follows:

[0061]

[0062] Where θ t The joint angle at time t, in degrees or rads, should be consistent in the implementation; Δt is the sampling period, in seconds; v t With a t These are angular velocity and angular acceleration, respectively, with units of deg / s and deg / s^2. This definition converts position quantities into velocity quantities, making it easier to characterize smoothness and the upper limit of change in probabilistic form.

[0063] Next, the physiologically mobile domain ROM is written as a probability factor in the form of a soft penalty, so that the prior automatically decays when the angle exceeds the limit. The calculation formula is as follows:

[0064]

[0065] Where θ min ,θ max These are the lower and upper physiological limits of the joint, in units of θ. t Consistency; λ rom >0 indicates the severity of the penalty for exceeding the boundary; (x) + =max(x,0) is the positive part operator. This term is approximately 1 within the interval, and decays exponentially by the square of the out-of-bounds value when it exceeds the bounds, which is equivalent to encoding the movable field as a priori probability.

[0066] To encourage natural continuity over time, a zero-mean Gaussian prior is given for angular velocity and angular acceleration, controlling the amplitude and frequency of variation. The calculation formula is as follows:

[0067]

[0068] Where σ v With σ a These represent the standard deviations of the prior velocity and acceleration, in units of deg / s and deg / s^2, respectively. Smaller values ​​indicate a stronger requirement for smoothness. They provide soft suppression of rapid abrupt changes while preserving the necessary degrees of freedom of the joints during rapid movements.

[0069] Regarding the temporal evolution of the angle itself, a high-order Markov mixture is used to provide a data-driven prediction for the next frame, taking into account gradual changes, first-order trends, and second-order corrections. The calculation formula is as follows:

[0070] p M (θ t |θ t-1:t-3 )=ω1N(θ t ;μ1,σ1 2 )+ω2N(θ t μ2,σ2 2 )+ω3N(θ t μ3,σ3 2 )

[0071] μ1=ρθ t-1 ,σ1 2 =(1-ρ 2 )σ v 2 Δt 2

[0072] trend = θ t-1 -θ t-2 μ2=θ t-1 +0.5trend t-1 ,σ2 2 =1.5σ v 2 Δt 2

[0073] accel t-1 =trend t-1 -(θ t-2 -θ t-3 ),μ3=θ t-1 +trend t-1 +0.3accel t-1 ,σ3 2 =2.0σ v 2 Δt 2

[0074] Where θ t-1:t-3 ={θ t-1 ,θ t-2 ,θ t-3} represents the angle history of the most recent three frames; ρ∈(-1,1) is the first-order autoregressive correlation coefficient; σ v The meaning of Δt is the same as before; trend t-1 It represents a first-order difference trend, with units consistent with angles; accel t-1 This is a second-order difference term, with units of angle. Weights w1, w2, w3 ≥ 0 and sum to 1, used to adaptively balance the three prediction modes under varying available history lengths. This prior provides robust mean and variance for angle predictions even with short-term historical noise.

[0075] Finally, the physiological boundary and Markov prediction are fused in the probability domain to obtain a comprehensive prior for subsequent inference:

[0076]

[0077] Where θ t The joint angle at time t, in degrees or rads, H t ={θ t-1 ,θ t-2 ,θ t-3} represents the set of historical perspectives used for prediction, Z t As a normalization constant, ensure that the right-hand side is normalized to θ. t The integral is 1, C ROM (θ t ) represents the physiologically mobile domain factor. This fusion expression unifies anatomical feasibility, temporal predictability, and local smoothness into a set of adjustable, interpretable, and implementable angular priors, providing a stable constraint basis for subsequent joint estimation with observations.

[0078] Step 5: Forward motion solution: The joint angles and limb lengths are mapped to the position and posture of the hand in the world coordinate system through forward geometry, and can be further transformed to the head-mounted display coordinate system, which is used as an analytical prediction quantity for comparison with sensor readings;

[0079] The idea is to establish a rigid body chain with the shoulder as the root node, and then sequentially connect the three-axis rotation of the shoulder, the translation of the upper arm, the flexion and extension rotation of the elbow, the translation of the forearm, and the extrinsic parameters of the handle installation to obtain the pose of the hand. This approach can be directly matched with the analytical Jacobian, which is convenient for subsequent numerical optimization or real-time filtering.

[0080] First, we use a simple and stable Euler cascade to represent joint rotation. The shoulder triaxial axis and elbow flexion and extension are written as follows:

[0081] R s =R z(θ s3 )R y (θ s2 )R x (θ s1 ),R e =R x (θ e ).

[0082]

[0083] Where R x (·),R y (·),R z (·) represent the basic rotation matrices about the local x, y, and z axes, respectively; θ s1 ,θ s2 ,θ s3 It is the three-axis joint angle of the shoulder, θ e R represents the elbow flexion-extension angle, measured in degrees or rads and kept consistent throughout the implementation. s With R e Let be the rotation matrices for the shoulder and elbow generated from the aforementioned angles. This expression directly transforms joint angles into rigid body pose, avoiding the overhead of exponential mapping while maintaining sufficient numerical stability.

[0084] Based on a series of rigid body chains, the formulas for calculating the hand position and attitude in the world frame are as follows:

[0085]

[0086] in This refers to the three-dimensional position of the shoulder in the world coordinate system, in meters (m); L upper ,L forearm These are the lengths of the upper arm and forearm, respectively, in meters (m). t represents the unit direction vector of each link in its local coordinate system, usually chosen as its respective axial basis vector; mnt With R mnt It is a fixed external parameter from the end of the forearm to the handle or racket mounting point, used to compensate for assembly geometric offset; and These represent the hand position and posture matrices in the world system. This formula intuitively demonstrates the cumulative effect of two translations after corresponding joint rotations, and also facilitates different assembly and left / right hand scenarios by simply replacing... It can be reused with external parameters.

[0087] To bring the world-system pose into the same coordinate domain as the head-mounted display output, the above pose can be transformed to the head-mounted display coordinate system H. The transformation formula is as follows:

[0088]

[0089] in With R WH These represent the translation of the head-mounted display in the world frame and the rotation from the world frame to the head-mounted display, respectively, derived from the external parameter estimates after time synchronization; This is the predicted pose that matches the sensor readings. The pose can also be converted into a unit quaternion as required by the interface.

[0090] This forward motion modeling method not only explicitly unfolds the skeletal chain geometry with minimal matrix operations, making computational complexity linearly related to the number of joints and meeting real-time application requirements, but also achieves clear decoupling of length and angle effects in its representation: length appears only as linear translation, and angle is only reflected in the rotation matrix, thus facilitating sensitivity analysis, system calibration, and assembly extrinsic parameter correction. Each variable in the model has a clear physical meaning, seamlessly carrying the physiological mobility domain and velocity boundaries given in step four, achieving a unified constraint on anatomical feasibility and motion coherence. Simultaneously, this representation is naturally adapted to the analytical Jacobian, and can be directly used for gradient-based optimization and uncertainty propagation, maintaining good convergence and time consistency in multi-round or sliding window scenarios, providing a stable and efficient computational foundation for subsequent Bayesian estimation, robust tracking, and motion evaluation.

[0091] Step Six: Bayesian Solution: Unify the forward motion prediction from Step Five with the priors given in Steps Three and Four into the same probabilistic framework. Construct an observation consistency metric based on multi-channel noise statistics. Achieve maximum a posteriori estimation through numerical optimization with physiological and velocimetric boundaries. Combine sliding window iteration and convergence criteria to output parameter estimation results and uncertainty quantification information. Incorporate post-processing smoothing into the Bayesian time series estimation model. Use high-order Markov priors to describe the continuity of angle, angular velocity, angular acceleration, and jerk, and introduce physiological boundary constraints to form a joint posterior.

[0092] like Figure 2 As shown, the left side represents the prior modeling. The length prior p(L|H,G) is derived from the user's height H and gender G, and is used to constrain the length vector L of each limb segment; the angle prior p(θ) t Characterizing the joint angle sequence θ tThe framework considers the physiological mobility domain and temporal continuity. The middle section represents the forward kinematics, where joint angles and lengths are substituted into pos = f(θ,L) to obtain the predicted pose pos of the hand (or end effector) (a 6DoF representation of position and attitude). Subsequently, an observation likelihood p(pos|θ,L) is established, and a multivariate Gaussian weighted model of the residual between the predicted pose and sensor observations is used to quantify observation consistency. The rightmost section represents the posterior inference p(θ,L|pos): the observation likelihood is added to two priors in the logarithmic domain and updated using maximum a posteriori or filtering to obtain estimates of joint angles and limb lengths and their uncertainties. This framework achieves a layered transfer of prior constraints, geometric mapping, and data consistency, supporting frame-by-frame or sliding-window recursive updates.

[0093] First, bias correction is performed on the observations and an error vector is defined:

[0094]

[0095] Where y t ∈R m For sensor observations at time t, b t For the calibrated system bias or gradual drift, For the corrected observation, θ t Let θ be the joint angle vector of this frame, L be the limb length vector, f(θ,L) be the prediction of the positive kinematics in the observation space from step five, and r be the joint angle vector of this frame. t This represents the difference between the observed and predicted values. If attitude is included, the rotation error can be expressed as the Lie algebra-minimum rotation vector. And stacked together with the position residual to r t .

[0096] Based on the noise statistics of each channel, a multivariate Gaussian likelihood is constructed, and the calculation formula is as follows:

[0097]

[0098] Where Σ t ∈R m×m To observe the noise covariance matrix, a common approach is to take a block-diagonal approach. σ p With σ R , respectively, represent the standard deviations of position and attitude noise, and I3 is a 3×3 identity matrix. The quadratic form in the exponent is the Mahalanobis distance; the larger the distance, the smaller the likelihood, reflecting the weighting of anisotropy and correlation.

[0099] Within a time window of length T, assuming independent conditions, the formula for calculating the total log-likelihood is as follows:

[0100]

[0101] Where θ 1:T ={θ1,…,θT} represents the angle sequence within the window, log∣Σ t The term | is retained during self-calibration noise to learn the covariance; if the noise is constant, it can be treated as a constant and neglected. The likelihood is added to the prior to obtain the log-posterior objective, which is then solved using the maximum a posteriori (MAP) formula, as follows:

[0102]

[0103] Where p prior (θ t H t H is the angle prior given in step four. t Let L represent the historical set used for prior evaluation, p(L|H,G) be the length prior in step three, and H and G be the user's height and gender, respectively. Maximum a posteriori (MAP) is equivalent to maximizing L. post Alternatively, its negative value can be minimized, and its calculation formula is as follows:

[0104]

[0105] The first term measures the consistency between observations and the model, while the second and third terms inject physiological mobility domain, temporal smoothing, and individualized length constraints into the solution, respectively. In practice, box constraints θ are applied to the variables. min ≤θ t ≤θ max To ensure the safety of the moving domain, iterative updates are performed using a bounded quasi-Newton or Gauss-Newton method; each iteration requires the Jacobian of the positive kinematics. and They approximate the first order of the residual as r. t +J t δθ t +G t δL, thus forming a standard quadratic approximation subproblem and guaranteeing convergence and stability.

[0106] Through the above process, the likelihood characterizes the support of sensor data for the model, the prior provides the inductive bias of physiological boundary and temporal regularity, the two are added in the logarithmic posterior and solved by MAP to obtain the angle and length estimates that are both consistent with observation and satisfy physiological and continuity constraints.

[0107] like Figure 3 As shown, the comparison results of the algorithm stability as a function of noise level are further presented to intuitively illustrate the performance of the above Bayesian solution under actual observation perturbations.

[0108] The experiment was set up as follows: zero-mean noise with different variances was progressively superimposed on the same set of original data while maintaining time synchronization. The stability index was calculated by weighting indicators such as sliding window convergence rate, mean Mahalanobis distance, ROM out-of-bounds rate, and phase drift, and then normalized to the 0–1 interval. The horizontal axis represents the noise level, and the vertical axis represents the stability index. The blue line represents the Bayesian method, and the red line represents the traditional inverse method. It can be seen that the traditional method rapidly destabilizes at extremely low noise levels, with the curve quickly dropping to near zero. The Bayesian method, on the other hand, decreases slowly with increasing noise, consistently showing significantly higher noise levels than the traditional method with less fluctuation in most noise intervals. This advantage stems from three combined factors: First, length and angle priors smoothly incorporate the physiologically mobile domain and temporal sequence into the target, blocking out-of-bounds errors and high-frequency spikes. Second, the covariance-weighted Mahalanobis distance in the likelihood provides directional suppression of multi-channel errors and automatically reduces the weight of low-confidence observations. Third, the sliding window MAP uses historical solutions as initial values, providing convergence traction and suppressing inter-frame error accumulation. Therefore, under the same observation conditions, the Bayesian framework is more robust to jitter, drift and local frame drops, and can stably output angle and length estimates that conform to physiological and continuity constraints.

[0109] Step 7: Motion Trajectory Visualization: Based on the joint angles and poses obtained in Step 6, design and implement dynamic visualization of the motion trajectory to present the single-round and multi-round action process in a readable and interactive way.

[0110] In this step, based on the joint angles and poses obtained in step six, a dynamic visualization of the motion trajectory is designed and implemented to present the single-round and multi-round motion process in a readable and interactive manner. First, by combining the obtained time-series data and event annotations, a visualization aligning joint angles with events is constructed. Figure 4 In this graph, the horizontal axis represents time, and the vertical axis represents angle. The upper panel displays the three-axis angle trajectory of the right shoulder, while the lower panel highlights the periodic changes in the right elbow flexion angle. Key event points and phase intervals are visually marked with dashed lines and color bands, with numerical labels at peaks and troughs to facilitate reading amplitude and rhythm. The physiological range of motion is given as a semi-transparent band, providing real-time indication of any boundary breaches or abnormal spikes. This graph emphasizes the regularity of angles over time and their correspondence with events, facilitating coaches or researchers to quickly identify rhythm deviations, insufficient amplitude, and phase drift. It also supports multi-round overlay comparisons to assess stability and consistency. It constructs a synchronous dynamic visualization of the skeleton and actual movements (e.g., Figure 5(As shown). The left side of the interface displays real video, while the right side presents upper limb movements driven by estimated angles using a 3D skeleton, simultaneously rendering the handle or racket's posture and end effector trajectory. Both sides share a unified timestamp and are aligned frame-by-frame to ensure a one-to-one correspondence between visual actions and the inferred skeleton. To enhance readability, the skeleton color changes with instantaneous velocity, the end effector trajectory displays the movement path as a trail, and the event list and progress bar use different colored segments to indicate each stage. Users can jump to key moments with a single click, loop through any interval, and quickly locate action details through perspective switching, skeleton transparency adjustment, and curve visibility. Frame interpolation is used for smooth transitions between adjacent keyframes to ensure animation continuity and reading comfort. The parameter panel at the bottom of the interface synchronously updates the angle, velocity, and stage indicators of the current frame, providing structured information that can be recorded and exported. Through the dual presentation of joint angle time series and skeleton-video, this step transforms numerical solutions into intuitive evidence, serving both technical evaluation and teaching correction, as well as providing a unified visual foundation for multi-round comparisons, training report generation, and algorithm iteration.

[0111] This invention also has the following advantages:

[0112] By integrating human biomechanical priors with high-order Markov timing constraints, sensor noise and drift are effectively suppressed, non-physical understanding is avoided, and the robustness of attitude estimation is improved.

[0113] Individualized physiological priors are constructed based on user height and gender to adapt to different body type characteristics and enhance the generalization ability of the method.

[0114] A Bayesian framework is used to achieve parameter estimation and uncertainty quantification, providing support for safe decision-making and adaptive control.

[0115] By using GPU parallelism and sliding window iteration, millisecond-level real-time updates are achieved to meet the low-latency requirements of VR scenarios.

[0116] Multi-dimensional visualization solutions intuitively present movement trajectories and key parameters, providing practical tools for sports training, skills assessment, and injury prevention.

[0117] The embodiments of the present invention have been described in detail above with reference to the accompanying drawings, but the present invention is not limited to the described embodiments. For those skilled in the art, various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention, and these variations still fall within the protection scope of the present invention.

Claims

1. A Bayesian-based method for reverse-engineering the trajectory of a VR table tennis ball, characterized in that: Includes the following steps: Step 1: Acquisition of raw VR motion data: Collect spatial position information, rotation posture information and event information of VR headset and controller, and form a data sequence organized by frame; Step 2: VR motion data preprocessing: The acquired VR data sequence is preprocessed, including time base unification, data framing, coordinate system integration, and extrinsic parameter calibration. The data sequence is resampled according to the target frame rate and short-term lost frames are filled in. Quaternion poses are normalized. Through anomaly detection, interpolation repair, zero bias estimation, and drift compensation, frame-by-frame confidence labels are output. Step 3: Construction of physiological prior knowledge: Population statistical information and physiological / kinematic knowledge are continuously encoded into "length prior + kinematic prior" using a unified Bayesian expression. Based on the user's height and gender, individualized mean and standard deviation of limb length are generated, proportional and symmetry constraints are established, and the boundaries of the joint physiological mobility domain are set to form computable parameter constraints. Step 4: Constructing Prior Knowledge of Angle Constraints: Construct a state vector containing angle, angular velocity, angular acceleration and jerk, introduce first-order, second-order and third-order higher-order Markov constraints, integrate soft physiological mobility domain constraints and temporal smoothing terms, and output window-level scores and penalty intensity through a "stationary-trend-correction" adaptive mode. Step 5: Forward motion solution: The joint angles and limb lengths are mapped to the position and posture of the hand in the world coordinate system through forward geometry, and can be further transformed to the head-mounted display coordinate system, which is used as an analytical prediction quantity for comparison with sensor readings; Step Six: Bayesian Solution: Unify the forward motion prediction from Step Five with the priors given in Steps Three and Four into the same probabilistic framework. Construct an observation consistency metric based on multi-channel noise statistics. Achieve maximum a posteriori estimation through numerical optimization with physiological and velocimetric boundaries. Combine sliding window iteration and convergence criteria to output parameter estimation results and uncertainty quantification information. Step 7: Motion Trajectory Visualization: Based on the joint angles and poses obtained in Step 6, design and implement dynamic visualization of the motion trajectory, presenting the movement process in a readable and interactive manner.

2. The Bayesian-based VR table tennis trajectory inversion method according to claim 1, characterized in that: Step two specifically includes: first, unifying the time base and data framing, synchronizing the timestamps of the head-mounted display and the two-handed controllers to a monotonic clock, resampling according to the target frame rate and filling in short-term dropped frames, and simultaneously completing coordinate system integration and extrinsic parameter calibration, unifying the poses of each device to the world or head-mounted display coordinate system and normalizing the quaternion poses; then, performing quality assessment and anomaly handling, using interquartile range and three times standard deviation rules to detect position and pose anomalies, performing linear or spline interpolation on isolated anomalies and missing segments, performing nearest neighbor interpolation and continuity for pose jumps, estimating zero bias in static segments, compensating for slow drift and random walks within a sliding window, and simultaneously outputting frame-by-frame confidence markers.

3. The Bayesian-based VR table tennis trajectory inversion method according to claim 1, characterized in that: Step three specifically includes: Using user height and gender as independent variables, a priori limb length is established, and physiological and velocimetric boundaries are simultaneously determined. Based on the population statistical baseline, the mean lengths of the neck, upper arm, forearm, and shoulder width are adaptively scaled according to height. The mean lengths of each part are linearly adjusted according to the deviation of the user's height from the reference height. Then, a gender correction is applied, with a gender coefficient set for different parts. The scaling factor is 1 for males and 0.9 to 0.95 for females to obtain individualized mean lengths. Then, the standard deviation is adaptively adjusted so that the standard deviation of length is linearly amplified with the degree of height deviation. When the height is close to the reference value, the standard deviation of length is adjusted accordingly. Standard deviation contraction increases as height deviation increases, reflecting individual differences and avoiding over-contraction of uncommon body types. Proportion and symmetry constraints are established, limiting the reasonable proportion range of the upper arm and forearm, setting an upper limit for the symmetry deviation of the left and right upper limbs, applying penalty weights to combinations exceeding limits to ensure geometric consistency, providing physiological mobility domain boundaries, setting upper and lower limits for the shoulder triaxial axis and elbow flexion and extension respectively, and making differentiated configurations according to the racket side and non-racket side to match the characteristics of the project, forming prior outputs, including the length prior distribution and covariance table, physiological and velocity boundary tables, and proportion and symmetry constraint coefficients.

4. The Bayesian-based VR table tennis trajectory inversion method according to claim 1, characterized in that: Step four specifically includes: A state vector containing angle, angular velocity, angular acceleration, and jerkiness is constructed. Using historical information from multiple adjacent frames as the temporal basis, higher-order Markov constraints of first-, second-, and third-order historical terms are introduced. Combined with soft physiological mobility domain constraints and temporal smoothing terms, sudden spikes and out-of-bounds behaviors are penalized, and window-level scores and penalty intensity are output. The above temporal priors are adaptively fused using three modes: "stationary," "trend," and "correction," to maintain continuity and interpretability under different states such as high-speed swings and static transitions.

5. The Bayesian-based VR table tennis trajectory reverse inference method according to claim 1, characterized in that: Step six specifically includes: The observations of position and attitude are compared with the forward kinematics predictions, and an observation consistency metric is constructed based on multi-channel noise statistics. The length prior and temporal prior obtained in steps three and four are jointly modeled within the same posterior framework. Numerical optimization with physiological and velocity boundaries is used to obtain the maximum posterior estimate of joint angle and limb length. Preferably, the solution of the previous window is used as the initial value in the form of a sliding window, and the target convergence and gradient threshold are used as stopping criteria. The curvature is approximated at the optimum to obtain the parameter uncertainty, which is used to generate confidence intervals and subsequent visualization confidence bands.

6. The Bayesian-based VR table tennis trajectory inversion method according to claim 1, characterized in that: Step seven specifically includes: Based on the joint angle and length parameters obtained by reverse engineering, a rigid body chain model with the shoulder joint as the root node is established. The three-axis rotation of the shoulder, the geometric connection between the upper arm and the forearm, and the installation of extrinsic parameters at the end are sequentially connected to obtain the spatial pose of the hand. The "length effect" and "angle effect" are decoupled in structure, which facilitates sensitivity analysis, extrinsic parameter correction and real-time calculation, and is naturally adapted to the analytical Jacobian to accelerate subsequent solutions. On this basis, a multi-view visualization is generated, rendering the end trajectory and velocity curve, and using color or markers to link key events and stage intervals.

Images