Craniofacial skeleton digital twin growth prediction system and method

Abstract

The invention discloses a craniofacial skeleton digital twin growth prediction system and method, and relates to the technical field of digital medical treatment and orthodontic intelligent auxiliary decision making. The system comprises a data acquisition module for acquiring multi-modal craniomaxillofacial data and dynamic biomechanical time sequence data; the data processing module is used for preprocessing and registering the acquired data; the digital twinborn body construction module constructs an individualized three-dimensional digital twinborn body model; the mixed growth prediction module adopts a mode of coupling data driving and mechanism simulation to simulate an accumulative growth effect under the action of dynamic load circulation and generate a growth deformation prediction result; and the visualization module visually displays the deformation process of the individualized three-dimensional digital twin model on a preset time sequence based on the growth deformation prediction result. According to the method, precise individualized dynamic prediction of the craniomaxillofacial growth trajectory and quantitative prediction of the treatment intervention effect are achieved, and the scientificity of clinical diagnosis and treatment scheme design, the reasonability of treatment opportunity selection and the reliability of prognosis evaluation are improved.

Description

Technical Field

[0001] This invention relates to the field of digital medical and orthodontic intelligent auxiliary decision-making technology, and in particular to a craniofacial skeleton digital twin growth prediction system and method. Background Technology

[0002] Accurate prediction of craniofacial growth and development is a core aspect of clinical diagnosis and treatment in orthodontics, orthognathic surgery, and other related fields. The stage of growth and development directly determines the choice of treatment plan, the timing of treatment, the assessment of prognosis, and the evaluation of relapse risk. Therefore, constructing an accurate predictive model that can characterize an individual's growth trajectory over time and reflect the coordinated evolution of bone and soft tissue is of significant theoretical and clinical value for improving clinical decision-making and optimizing treatment outcomes.

[0003] Currently, existing craniofacial growth prediction technologies still have the following limitations. First, most existing technologies rely on static cross-sectional CBCT data for malocclusion classification studies, lacking the ability to model individual growth trajectories over time. While existing technologies achieve high accuracy in predicting mandibular growth trends using lateral cephalometric radiographs, these two-dimensional images cannot retain three-dimensional spatial structural information and facial soft tissue morphology data. Therefore, single-modal or static data cannot fully depict the dynamic evolution of the craniofacial three-dimensional structure over time, resulting in incomplete prediction information and limited accuracy. Second, existing technologies primarily focus on regression prediction of morphological features, failing to incorporate core factors such as the biomechanical environment (e.g., masticatory force, muscle force distribution) and individual genetic predisposition that determine growth direction and rate into the modeling system. This renders the prediction model a black box, lacking physical rationality and clinical interpretability, and making it difficult to simulate the coordinated evolution of bone and soft tissue under growth and treatment forces. Third, existing technical solutions are mostly based on the statistical laws of the population to build predictive models, ignoring the significant differences between individuals in terms of anatomical morphology, growth potential, and functional environment. Although they can reflect the macro growth trend, they are difficult to accurately depict the unique growth trajectory of an individual and cannot provide reliable quantitative basis for key clinical decisions such as orthodontic treatment plan design and surgical timing. Summary of the Invention

[0004] To address the shortcomings of existing technologies, such as limited data dimensions, lack of biomechanical mechanisms, and insufficient individualization, this invention proposes a craniofacial skeletal digital twin growth prediction system and method. Through multimodal data fusion, individualized 3D digital twin construction, hybrid prediction combining data-driven and mechanism simulation, dynamic functional simulation, and closed-loop interactive optimization, it achieves accurate individualized dynamic prediction of craniofacial growth trajectory and quantitative prediction of treatment intervention effects, thereby improving the scientific nature of clinical treatment plan design, the rationality of treatment timing selection, and the reliability of prognostic assessment.

[0005] The present invention achieves the above objectives through the following technical solutions:

[0006] A craniofacial skeletal digital twin growth prediction system includes:

[0007] The data acquisition module is used to acquire multimodal craniofacial data and dynamic biomechanical time-series data under physiological activity states; the multimodal craniofacial data includes at least cone-beam CT data, three-dimensional facial surface scan data, clinical phenotype data and genetic information data;

[0008] The data processing module is used to preprocess multimodal craniofacial data and dynamic biomechanical time-series data, and to register the preprocessed cone-beam CT data with the preprocessed three-dimensional facial surface scan data.

[0009] The digital twin construction module is used to construct an individualized three-dimensional digital twin model containing bones, jawbone, and facial soft tissue based on processed cone-beam CT data and processed three-dimensional facial surface scan data.

[0010] The hybrid growth prediction module is used to simulate the cumulative growth effect under dynamic load cycles based on an individualized three-dimensional digital twin model, combined with processed clinical phenotypic data, processed dynamic biomechanical time series data and processed genetic information data. It adopts a data-driven and mechanism simulation coupled approach to generate growth deformation prediction results of the craniofacial region.

[0011] The visualization module is used to visualize the deformation process of an individualized 3D digital twin model over a preset time series based on the growth deformation prediction results.

[0012] As a preferred embodiment of the present invention, the data acquisition module is equipped with an oral cavity flexible sensor and / or a mandibular motion capture device;

[0013] The clinical phenotypic data includes tooth arrangement characteristics, occlusal parameters, and craniofacial growth type classification data.

[0014] The genetic information data includes genome-wide association analysis data or multigene risk score data;

[0015] The dynamic biomechanical time series data includes chewing force time series and mandibular movement trajectory time series.

[0016] As a preferred embodiment of the present invention, the data processing module includes:

[0017] The data preprocessing unit is used to preprocess multimodal craniofacial data and dynamic biomechanical time-series data; the preprocessing includes outlier detection, noise reduction and format conversion;

[0018] The cross-modal registration unit is used to unify the preprocessed cone-beam CT data and the preprocessed 3D facial surface scan data into the same spatial coordinate system based on the deep learning feature point matching algorithm, thus completing the registration.

[0019] As a preferred embodiment of the present invention, the digital twin construction module is configured to perform the following steps:

[0020] Based on the processed cone-beam CT data, surface mesh models of the skull and jawbone are reconstructed as initial static geometric models of the skeleton.

[0021] The morphological contour, elastic modulus, and Poisson's ratio parameters of facial soft tissue were obtained from the processed three-dimensional facial surface scanning data. A static geometric model of facial soft tissue was constructed based on the attachment relationship between soft and hard tissues.

[0022] The mechanical coupling relationship between the initial static geometric model of the skeleton and the static geometric model of the facial soft tissue was established through finite element analysis.

[0023] Based on the mechanical coupling relationship, the initial static geometric model of the skeleton is fused with the static geometric model of the facial soft tissue to construct an individualized three-dimensional digital twin model containing the skeleton, jawbone and facial soft tissue.

[0024] As a preferred embodiment of the present invention, the mixed growth prediction module includes:

[0025] The data-driven prediction unit is used to input the individualized three-dimensional digital twin model, processed clinical phenotypic data, and processed genetic information data into the deep learning model and output the initial growth deformation field.

[0026] The mechanism-driven simulation unit is used to construct a time-series load spectrum from the processed dynamic biomechanical time-series data. The time-series load spectrum is applied as a boundary condition to the individualized three-dimensional digital twin model. The model is solved by a finite element solver to obtain the time-series stress field distribution and instantaneous deformation field. Based on the solution results and the biomechanical theory of bone metabolism, the bone micro-remodeling rate field is obtained.

[0027] The coupled computational unit is used to regularize and correct the initial growth deformation field by taking the temporal stress field distribution, instantaneous deformation field and bone micro-remodeling rate field as physical constraints, and generate the final growth deformation prediction result.

[0028] As a preferred embodiment of the present invention, the instantaneous deformation field includes at least the three-dimensional displacement values ​​of the global mesh nodes of the individualized three-dimensional digital twin model; the global mesh nodes include bone hard tissue mesh nodes and facial soft tissue mesh nodes;

[0029] The time-series stress field distribution includes at least the equivalent stress values ​​of the bone hard tissue mesh nodes.

[0030] As a preferred embodiment of the present invention, the mechanism-driven simulation unit obtains the bone micro-remodeling rate field and is configured to perform the following steps:

[0031] Based on the three-dimensional displacement values ​​in the instantaneous deformation field, the equivalent stress values ​​of bone hard tissue mesh nodes whose deformation exceeds the physiological limit index are eliminated; the physiological limit index is set according to the physiological range of craniofacial bones.

[0032] Extract all remaining equivalent stress values ​​of each bone hard tissue mesh node throughout the entire time series, and weight the time points according to the chewing cycle to calculate the average equivalent stress value of each bone hard tissue mesh node;

[0033] Based on the average equivalent stress value and the preset bone remodeling stress threshold, the stress partition type of each bone hard tissue mesh node is determined; the preset bone remodeling stress threshold is set based on the physiological remodeling characteristics of craniofacial bones.

[0034] Based on the stress partitioning type and average equivalent stress value of each bone hard tissue mesh node, a preset nonlinear bone remodeling rate model is used to calculate the bone micro remodeling rate.

[0035] The number, three-dimensional displacement coordinates, average equivalent stress value, stress partition type, and bone micro-remodeling rate of each bone hard tissue mesh node are associated and mapped to the global mesh nodes of the individualized three-dimensional digital twin model to obtain the bone micro-remodeling rate field.

[0036] In a preferred embodiment of the present invention, the coupling calculation unit is configured to perform the following steps:

[0037] The initial growth deformation field, temporal stress field distribution, instantaneous deformation field, and bone micro-remodeling rate field are normalized.

[0038] Based on the normalized temporal stress field distribution, instantaneous deformation field, initial growth deformation field, and bone micro-remodeling rate field, a regularized objective function incorporating spatially adaptive physical constraints is constructed:

[0039]

[0040] in, The value of the regularization objective function; The normalized and corrected growth deformation field to be solved; This represents the initial growth deformation field after normalization. To correspond with the normalized equivalent stress distribution of the growth deformation field; The normalized temporal stress field distribution; The spatially adaptive weight matrix is ​​set based on clinical anatomical data and biomechanical experimental results; To normalize and correct the Laplace operator of the growth deformation field; This is the stress constraint coefficient; The coefficient is the spatial smoothing constraint coefficient. The deformation gradient constraint coefficient; The coefficient constraining the rate of bone remodeling; To normalize and correct the displacement gradient values ​​of the growth deformation field; This represents the displacement gradient value of the instantaneous deformation field after normalization. To normalize the first derivative of the growth deformation field with respect to time; The normalized skeletal microremodeling rate field; This is element-wise multiplication; It is the Euclidean norm;

[0041] The regularization objective function is iteratively optimized using the conjugate gradient method. When the change in deformation field between adjacent iterations is less than the preset deformation threshold, the normalized corrected deformation field is obtained.

[0042] The normalized modified deformation field is denormalized to output the final growth deformation prediction result.

[0043] As a preferred embodiment of the present invention, it further includes:

[0044] The interaction and feedback module is used to trigger the hybrid growth prediction module to regenerate the growth deformation prediction results based on the user's input treatment intervention parameters when the user inputs treatment intervention parameters, and to update the regenerated growth deformation prediction results to the visualization module in real time for display.

[0045] The training data augmentation module is used to use generative adversarial networks to complete the missing image data between adjacent acquisition time points in the collected historical patient clinical follow-up data, and to construct a continuous longitudinal time sequence for training the deep learning model in the hybrid growth prediction module.

[0046] A method for predicting the growth of a digital twin of the craniofacial skeleton includes:

[0047] Collect multimodal craniofacial data and dynamic biomechanical time-series data under physiological activity states; the multimodal craniofacial data includes at least cone-beam CT data, three-dimensional facial surface scan data, clinical phenotype data and genetic information data;

[0048] Multimodal craniofacial data and dynamic biomechanical time-series data were preprocessed, and the preprocessed cone-beam CT data were registered with the preprocessed three-dimensional facial surface scan data.

[0049] Based on the processed cone-beam CT data and the processed three-dimensional facial surface scan data, an individualized three-dimensional digital twin model containing bones, jawbone and facial soft tissues was constructed.

[0050] Based on a personalized three-dimensional digital twin model, and combined with processed clinical phenotypic data, processed dynamic biomechanical time series data and processed genetic information data, a data-driven and mechanism simulation coupled approach is used to simulate the cumulative growth effect under dynamic load cycle and generate growth deformation prediction results for the craniofacial region.

[0051] Based on the growth deformation prediction results, the deformation process of the individualized 3D digital twin model over a preset time series is visualized.

[0052] The beneficial effects of this invention are as follows: Through the synergistic effect of the data acquisition module and the data processing module, comprehensive acquisition of multimodal static anatomical data, dynamic biomechanical functional data, genetic information data, and clinical phenotypic data of the craniofacial region is achieved. After preprocessing and registration techniques, a standardized dataset with a unified format and logical correlation is output, providing a complete and accurate data foundation for subsequent modeling and prediction. This effectively solves the problem of limited prediction accuracy caused by the single information dimension and spatiotemporal misalignment of data in existing technologies. Based on this, the digital twin construction module, by constructing individualized three-dimensional digital twin models, concretizes abstract data into digital entities with precise anatomical structures and biomechanical properties, providing a simulable and driveable physical carrier for subsequent growth prediction. Furthermore, the hybrid growth prediction module adopts a data-driven and mechanism simulation coupled approach to simulate the cumulative growth effect under dynamic load cycles. This fully leverages the powerful ability of deep learning to uncover individualized patterns while ensuring that the prediction results conform to the laws of biomechanical conservation, achieving a balance between accuracy and reliability. It successfully transforms the prediction model from a black box into a white box with physical interpretability. Meanwhile, by adopting a scientific approach of modeling first and then predicting, the prediction results are directly reflected in the vertex displacement field of the individualized 3D digital twin model. This effectively avoids the shortcomings of existing technologies that involve prediction followed by reconstruction, such as the disconnect between the results and the model, and the information loss introduced by secondary reconstruction. It achieves a precise depiction from the average pattern of the population to the unique growth trajectory of each individual. Finally, the visualization module intuitively displays the continuous deformation of the individualized 3D digital twin model over a preset time series, realizing a dynamic and quantitative presentation of the growth process. This provides efficient, accurate, and interactive technical support for clinical diagnosis, treatment design, and doctor-patient communication. Attached Figure Description

[0053] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the description of the embodiments 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. Wherein:

[0054] Figure 1 A schematic diagram of the modular structure of a craniofacial skeleton digital twin growth prediction system provided by the present invention;

[0055] Figure 2 This is a schematic diagram of the internal structure of the hybrid growth prediction module in an embodiment of the present invention;

[0056] Figure 3 The flowchart illustrates a method for predicting the growth of a digital twin of the craniofacial skeleton provided by this invention. Detailed Implementation

[0057] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the described embodiments of the present invention are within the scope of protection of the present invention.

[0058] like Figure 1 As shown, this is an embodiment of the present invention, which provides a craniofacial skeleton digital twin growth prediction system. The system includes a data acquisition module, a data processing module, a digital twin construction module, a hybrid growth prediction module, a visualization module, an interaction and feedback module, and a training data augmentation module. The specific structure and function of each module are as follows:

[0059] The data acquisition module is used to collect multimodal craniofacial data and dynamic biomechanical time-series data under physiological activity states.

[0060] The multimodal craniofacial data includes at least cone-beam computed tomography (CBCT) data, three-dimensional facial surface scan data, clinical phenotype data, and genetic information data; the clinical phenotype data includes tooth arrangement characteristics, occlusal parameters, and craniofacial growth type classification data; and the genetic information data includes genome-wide association analysis data or multigene risk score data.

[0061] Specifically, cone-beam CT data is acquired using a dedicated dental cone-beam CT scanner, containing three-dimensional anatomical data of hard tissues such as the skull, jawbone, dentition, and alveolar bone, and stored in DICOM format after acquisition. Three-dimensional facial surface scan data is acquired using a structured light 3D scanner, containing the morphological contours, texture information, and surface deformation data under static or dynamic facial expressions of facial soft tissues. Geometric data in STL format and texture data in JPG format are output simultaneously after acquisition. Clinical phenotypic data is extracted using a combination of oral endoscopy and cephalometric radiographs. This data can be structured or unstructured text data. The oral endoscope is used to record tooth alignment characteristics such as crowding and rotation, as well as occlusal parameters such as overbite height and overjet distance. The cephalometric radiographs are used to calculate indices such as SNA angle, SNB angle, and ANB angle to determine craniofacial growth type classifications such as vertical growth type, horizontal growth type, and average growth type. If the data is unstructured text, the data preprocessing unit will extract structured data through natural language processing. Genetic information data were collected using an Illumina GSA chip and a matching detection platform to characterize an individual's craniofacial growth genetic potential.

[0062] The data acquisition module is equipped with a flexible oral sensor and / or a mandibular motion capture device; the dynamic biomechanical time series data includes masticatory force time series and mandibular motion trajectory time series.

[0063] Specifically, the masticatory force time series consists of masticatory force vectors acquired sequentially over time, each vector containing information on the direction, magnitude, and frequency of force application. When using a flexible oral sensor for acquisition, a flexible thin-film pressure sensor array is attached to the occlusal surface of the teeth and the corresponding periodontal mucosa area, allowing for the capture of force signals during occlusion and mastication without delay. The mandibular movement trajectory time series consists of mandibular spatial position data acquired sequentially over time, including the trajectory coordinates of mandibular opening and closing, protrusion, and lateral movements. When using a mandibular movement capture device for acquisition, the spatial position and movement trajectory of the mandible are obtained in real time through non-invasive methods such as extraoral optical capture, ultrasonic positioning, or electromagnetic tracking, without affecting the patient's normal physiological activities.

[0064] In one specific embodiment, the oral flexible sensor uses a flexible pressure sensor array made of polyimide material with a resolution of 0.1 N. It is attached to the occlusal surface of the upper and lower jaws and the sampling frequency is set to 100 Hz. It simultaneously collects the direction, magnitude, and temporal distribution of masticatory force. The 100 Hz sampling frequency is set based on the physiological frequency characteristics of craniofacial masticatory movements, which can completely capture the dynamic changes of masticatory force. The mandibular motion capture device uses the OptiTrack Flex 13 optical motion capture system, equipped with 6 infrared cameras. The sampling frequency is set to 120 Hz, and the spatial tracking accuracy is set to 0.1 mm. By attaching 3 reflective markers to the patient's mandibular angle and chin, the three-dimensional trajectory temporal data of mandibular opening and closing, protrusion, and lateral movements are recorded without interference. The 120 Hz sampling frequency is higher than the masticatory movement frequency to avoid loss of motion trajectory data, and the 0.1 mm is in line with the measurement requirements of craniofacial anatomy.

[0065] All dynamic data is sampled continuously at equal time intervals to form a standardized time-series data stream with timestamps. After data collection is complete, all types of data are uniformly transmitted to a buffer with data backup capabilities to prevent data loss.

[0066] The data processing module is used to preprocess multimodal craniofacial data and dynamic biomechanical time-series data, and to register the preprocessed cone-beam CT data with the preprocessed three-dimensional facial surface scan data.

[0067] The data preprocessing unit is used to preprocess multimodal craniofacial data and dynamic biomechanical time-series data under physiological activity states. Preprocessing includes outlier detection, noise reduction, and format conversion.

[0068] Specifically, for outlier detection, the operation is as follows: For cone-beam CT data, an outlier detection algorithm based on statistical distribution is used to detect outliers with gray values ​​exceeding the specified threshold. The range of voxels is considered an outlier, where, The average grayscale value. The standard deviation of the grayscale value is used. For 3D facial surface scan data, an out-of-radius point removal algorithm is used, with a neighborhood search radius of 0.5 mm. Points with fewer than 30 neighboring points are identified as isolated outliers. For clinical phenotype data, if the data is structured, a box plot method is used to mark data exceeding 1.5 times the interquartile range as abnormal; if the data is unstructured text, a natural language processing algorithm is used to identify invalid information such as garbled characters and duplicate records. For genetic information data, PLINK software is used to perform quality control on whole-genome SNP data acquired from Illumina GSA chips, removing data with a genotyping success rate of less than 95%, a minimum allele frequency of less than 1%, and a Hardy-Weinberg equilibrium P-value that is too small. SNP loci were identified, and individuals with a sample missing rate greater than 5% were labeled. Anomaly detection was performed using a sliding window method. Specifically, a sliding window of fixed time length was set, and local means were calculated based on the data within the current sliding window. and local standard deviation It will exceed Data within a certain range is considered abnormal. The sliding window size is set based on the sampling frequency; a value of 100 milliseconds is recommended.

[0069] For noise reduction, the specific operations are as follows: For cone-beam CT data, 3×3×3 adaptive median filtering is used in abnormal areas, and Gaussian filtering with a standard deviation of 0.8 is used in non-abnormal areas to preserve bone edge details; for 3D facial surface scan data, moving least squares is used for smoothing to preserve the natural contours and texture details of facial soft tissues; for clinical phenotype data, if it is structured data, nearest neighbor interpolation with a K value of 5 is used to correct outliers; if it is unstructured text data, special characters and whitespace characters are removed using regular expressions, and jieba word segmentation is used to segment the text and filter stop words; for genetic information data, IMPUTE2 software is used, based on 1000 Genomes Project Phases. The III database fills in missing SNP sites, retaining only sites with an info value greater than 0.4 after filling. For dynamic biomechanical time-series data, linear interpolation is used to correct outliers, and a Butterworth low-pass filter is used to filter high-frequency noise from the overall data. The cutoff frequency is set to 10Hz, which is determined based on the highest frequency of craniofacial physiological movements, effectively filtering out high-frequency interference generated during equipment acquisition.

[0070] The specific operations for format conversion are as follows: Cone-beam CT data is converted to NIfTI format, and grayscale values ​​are normalized to the range [0, 255]. For 3D facial surface scan data, geometric data is converted to PLY format, and texture data is converted to PNG format. The two are linked through vertex indexing to ensure consistency in texture mapping. For clinical phenotype data, if it is structured data, it is converted to CSV format, with fields including patient ID, tooth alignment features, occlusal parameters, craniofacial growth type classification, and acquisition timestamp; if it is unstructured text data, structured information is extracted and converted to JSON format, and the data is linked to the CSV format data through the patient ID. For genetic information data, SNP data is converted to VCF format, and polygenic risk score data is converted to TXT format, with fields including patient ID, SNP locus number, genotype, and polygenic risk score value. Dynamic biomechanical time-series data is converted to HDF5 format, indexed by patient ID, data type, and timestamp, and the sampling frequency is uniformly set to 100Hz.

[0071] The cross-modal registration unit is used to unify the preprocessed cone-beam CT data and the preprocessed 3D facial surface scan data into the same spatial coordinate system based on the deep learning feature point matching algorithm, thus completing the registration.

[0072] Specifically, the preprocessed cone-beam CT data is segmented using a 3D U-Net neural network to extract feature points such as the nasal root, bilateral infraorbital foramina, bilateral tragus, bilateral mandibular angles, and chin apex. The three-dimensional coordinates of each feature point in the cone-beam CT coordinate system are then recorded. All coordinate units are millimeters. Simultaneously, the PointNet++ point cloud feature extraction network is used to automatically extract matching feature points in the corresponding anatomical regions of the preprocessed 3D facial surface scan data, and the 3D coordinates of each feature point in the facial scan coordinate system are recorded. All coordinate units are millimeters. Based on the correlation matrix between feature point pairs, the rigid body transformation matrix, which includes the rotation matrix, is solved using the least squares method. Translation vector The translation vector is in millimeters, and the calculation formula is:

[0073]

[0074] Substituting the coordinates of all vertices of the 3D facial surface scan data into the aforementioned rigid transformation matrix, spatial registration with the cone-beam CT data is completed. Using a unified timestamp as a reference, the preprocessed clinical phenotypic data and processed genetic information data are used as attribute information and associated with the corresponding patient ID and time node, ultimately outputting a standardized dataset with spatiotemporal alignment and uniform format.

[0075] The digital twin construction module is used to build an individualized 3D digital twin model containing bones, jawbone, and facial soft tissues based on processed cone-beam CT data and processed 3D facial surface scan data.

[0076] The digital twin building block is configured to perform the following steps:

[0077] Based on the processed cone-beam CT data, surface mesh models of the skull and jawbone are reconstructed as initial static geometric models of the skeleton.

[0078] Specifically, threshold segmentation is performed on cone-beam CT data, and a threshold range for bone tissue is set according to Henle units (e.g., voxels with Henle unit values ​​greater than 300 are initially labeled as bone tissue) to distinguish bone tissue, tooth structure, alveolar bone, and surrounding cavity structures. The moving cube algorithm is used for three-dimensional surface reconstruction to generate initial triangular mesh models of the skull and jawbone. The initial triangular mesh models are smoothed (e.g., Laplacian smoothing), deburred, simplified (e.g., quadratic edge folding simplification algorithm), and repaired to eliminate non-manifold edges, voids, and overlapping patches, ensuring that the mesh model quality meets the requirements of finite element simulation. The reconstruction results include the complete anatomical morphology of the skull base, maxilla, mandible, temporomandibular joint region, dentition, and alveolar bone, and finally output an initial static geometric model of the skeleton in STL format.

[0079] The morphological contour, elastic modulus, and Poisson's ratio parameters of facial soft tissue are obtained from the processed 3D facial surface scan data. A static geometric model of facial soft tissue is constructed based on the attachment relationship between soft and hard tissues.

[0080] Specifically, facial soft tissue contour data is extracted from three-dimensional facial surface scan data, and biomechanical material properties such as elastic modulus and Poisson's ratio are obtained by combining clinical common knowledge and standard parameter libraries in this field. Based on the facial soft tissue contour data, a closed curved surface mesh is generated to ensure that the resolution of the mesh matches that of the hard tissue. Finally, an STL format static geometric model of facial soft tissue is constructed, which can intuitively present the morphological features of areas such as lips, nose, chin, and cheeks.

[0081] The mechanical coupling relationship between the initial static geometric model of the skeleton and the static geometric model of the facial soft tissue was established through finite element analysis.

[0082] Specifically, the nearest neighbor search algorithm is used to search for the node with the closest Euclidean distance on the inner surface of the facial soft tissue static geometric model, with the surface node of the initial skeletal static geometric model as the target point, and to establish a spatial correspondence table between the skeletal hard tissue mesh node and the facial soft tissue mesh node.

[0083] Based on this spatial correspondence table, in anatomical attachment areas such as alveolar ridge, nasal base, infraorbital margin, and mandibular margin, binding contact pairs are defined, with the bone surface as the master surface and the soft tissue inner surface as the slave surface, and the contact attribute is set as binding constraint; in non-attachment areas, separable contact pairs are defined, with the normal as hard contact and the tangential as penalty friction, allowing limited slippage of soft tissue relative to bone.

[0084] A load node set is created on the occlusal surface of the teeth and the periodontal ligament region. This node set is associated with the load boundary of the initial skeletal static geometry model and is called by the mechanism-driven simulation unit of the hybrid growth prediction module when applying the time-series load spectrum.

[0085] Based on the mechanical coupling relationship, the initial static geometric model of the skeleton is fused with the static geometric model of the facial soft tissue to construct an individualized three-dimensional digital twin model containing the skeleton, jawbone and facial soft tissue.

[0086] Specifically, the initial static geometric model of the skeleton and the static geometric model of the facial soft tissue are combined in a unified cone-beam CT coordinate system to construct a complete craniofacial anatomical model that includes both hard and soft tissues. At the same time, the spatial correspondence table established in the mechanical coupling relationship, the defined contact attributes, and the created load node set are stored as additional attributes of the anatomical model, ultimately forming an individualized three-dimensional digital twin model. This model has accurate anatomical structure, biomechanical material properties, and mechanical coupling relationship between hard and soft tissues, and can be directly used for subsequent dynamic load simulation and bone micro-remodeling simulation.

[0087] The hybrid growth prediction module is used to simulate the cumulative growth effect under dynamic load cycles based on an individualized three-dimensional digital twin model, combined with processed clinical phenotypic data, processed dynamic biomechanical time series data, and processed genetic information data. It adopts a data-driven and mechanism simulation coupled approach to generate growth deformation prediction results for the craniofacial region.

[0088] like Figure 2 As shown, the hybrid growth prediction module includes a data-driven prediction unit, a mechanism-driven simulation unit, and a coupled computation unit. The specific structure and workflow of each unit are as follows:

[0089] The data-driven prediction unit is used to input the individualized 3D digital twin model, processed clinical phenotypic data, and processed genetic information data into the deep learning model and output the initial growth deformation field.

[0090] Specifically, the deep learning model employs a hybrid architecture based on 3DCNN and Transformer to extract local geometric features and global correlation features from multimodal data, thereby improving the prediction accuracy of the initial growth deformation field. The internal structure of the deep learning model is as follows:

[0091] The input layer contains three input branches, corresponding to three types of input data: the individualized 3D digital twin model is converted into a 1×256×256×128 4D floating-point tensor, the processed clinical phenotype data is converted into a 1×21 2D floating-point tensor, and the processed genetic information data is converted into a 1×87 2D floating-point tensor.

[0092] The feature extraction layer comprises three independent feature extraction sub-networks, each corresponding to one of the three branches of the input layer, and all output a fixed-dimensional one-dimensional feature vector. Branch 1 is a 3DCNN feature extraction sub-network that extracts high-dimensional growth-related features from low-dimensional geometric features through convolutional layers, pooling layers, dilated convolutional layers, and global pooling layers, outputting a 1024-dimensional morphological feature vector, including jawbone growth potential, facial asymmetry, and quantitative indicators of craniofacial growth patterns. Branch 2 is a fully connected feature extraction sub-network that extracts core features directly related to craniofacial growth from clinical phenotypic data through max-min normalization and feature mapping, outputting a 256-dimensional phenotypic feature vector, including tooth crowding, occlusal plane inclination, and quantitative values ​​of craniofacial growth patterns. Branch 3 is an attention-based selection and fully connected feature extraction sub-network that selects 50 core SNP loci through an attention mechanism and extracts core genetic features, outputting a 512-dimensional genetic feature vector, including information on bone metabolism capacity, craniofacial morphological development potential, and genetic predisposition to growth rate.

[0093] The multimodal feature fusion layer adopts a single-module cascaded structure, comprising a feature concatenation layer and a deep fusion layer. The feature concatenation layer concatenates 1024-dimensional, 256-dimensional, and 512-dimensional feature vectors into a 1792-dimensional fused feature tensor. The deep fusion layer is a single fully connected layer with 1024 neurons, followed by a ReLU activation layer and a Dropout layer, completing the deep mapping and dimensionality reduction of the features, outputting a 1×1024 multimodal deep feature vector. The dropout rate of the Dropout layer can be set to 0.2.

[0094] The Transformer global association learning layer employs a lightweight Transformer encoder structure, comprising a positional encoding layer and two Transformer encoder layers. Each Transformer encoder layer includes an 8-head multi-head self-attention mechanism and a 2048-neuron feedforward neural network to capture global associations among morphological, phenotypic, and genetic features. The positional encoding layer adds virtual positional information through sine and cosine encoding to enhance the association representation. The final output is a 1×1024 feature vector that integrates global association patterns.

[0095] The prediction output layer adopts a structure of fully connected layers and a dimension matching module connected in series. The fully connected layer consists of three layers: the first layer has 1024 neurons, the second layer has 512 neurons, and the third layer has 3N neurons. Each layer is connected to a ReLU activation layer, where N is the total number of grid nodes in the individualized 3D digital twin model. The dimension matching module uses a grid node index mapping table to decompose the 1×3N tensor into an N×3 displacement matrix, where each row of the matrix corresponds to the growth 3D displacement value of a grid node. Finally, the initial growth deformation field in VTK format is output.

[0096] The deep learning model uses mean squared error as the loss function and takes the error between the predicted three-dimensional displacement value of the initial growth deformation field and the measured displacement value during clinical follow-up as the optimization objective. The Adam optimizer is used with a learning rate of 0.0001, a batch size of 10, and 1000 training rounds. An early stopping strategy is introduced, triggered by the validation set loss not decreasing for 10 consecutive rounds. Training is completed based on a continuous longitudinal time series constructed by the training data augmentation module to ensure the physiological rationality and individualized adaptability of the prediction results.

[0097] The mechanism-driven simulation unit is used to construct a time-series load spectrum from the processed dynamic biomechanical time-series data. The time-series load spectrum is applied as a boundary condition to the individualized three-dimensional digital twin model. The model is solved by a finite element solver to obtain the time-series stress field distribution and instantaneous deformation field. Based on the solution results and the biomechanical theory of bone metabolism, the bone micro-remodeling rate field is obtained.

[0098] Specifically, a cubic spline interpolation algorithm is used to synchronize the chewing force time series and the mandibular movement trajectory time series in the dynamic biomechanical time series data to the same time axis, resulting in a discrete data point sequence, which is then converted into a discrete time-load table that can be recognized by the finite element solver, i.e., a time series load spectrum. The fields include timestamp, chewing force magnitude, direction, and coordinates of mandibular feature points.

[0099] Mechanical boundary conditions, displacement boundary conditions, and ENCASTRE fixed constraints are all applied to the mesh nodes of the individualized 3D digital twin model, and the mechanical boundary conditions are precisely correlated with the previously constructed load node set. Based on the parameters of each time node in the time-series load spectrum, the mechanical boundary conditions apply concentrated forces to the occlusal contact point mesh nodes corresponding to the maxillary and mandibular dentition in the individualized 3D digital twin model. The direction and magnitude of the force strictly match the time-series load spectrum data. The occlusal contact point mesh nodes are determined as follows: based on the patient's clinical occlusal record and the anatomical morphology of the model's dentition, a spatial distance detection algorithm is used to extract mesh nodes with a spatial distance less than or equal to 0.1 mm on the occlusal surface of the maxillary and mandibular dentition as occlusal contact points, and these are matched vertically according to the dentition alignment to form a set of masticatory force application points. Based on the coordinates of mandibular feature points at each time node of the time-series load spectrum, displacement boundary conditions are applied to the corresponding feature point mesh nodes of the mandible in the individualized 3D digital twin model. This drives the mandible to complete opening and closing, protrusion, and lateral movements synchronized with masticatory forces, while simultaneously constraining the relative positions of the internal mesh nodes of the mandible to remain unchanged, ensuring its rigid body motion characteristics. The feature point mesh nodes are the nearest neighbor mesh nodes of key anatomical landmarks of the mandible, such as the bilateral mandibular angles, chin apex, and bilateral condylar apexes. The spatial positions of anatomical landmarks are located using the standard craniofacial anatomical coordinate library already labeled in the individualized 3D digital twin model. The mesh node with the smallest Euclidean distance to that position in the individualized 3D digital twin model is then extracted as the node to which the displacement boundary conditions are applied.

[0100] An ENCASTRE fixed constraint is applied to the occipital region of the skull base to completely lock all six degrees of freedom (translation in the X, Y, and Z directions and rotation around the X, Y, and Z axes) of all mesh nodes in this region. This simulates the physiological rigid connection between the skull and the cervical spine, prevents overall rigid body displacement, and ensures the stability and physiological rationality of the finite element solution.

[0101] An explicit dynamic finite element method (FEB) solver, using the central difference method as the core algorithm, was employed to numerically solve the dynamic mechanical response of the craniofacial region, yielding the temporal stress field distribution and instantaneous deformation field. The temporal stress field distribution includes fields such as timestamp, hard tissue mesh node number, equivalent stress value, and maximum principal stress direction vector, reflecting the temporal variation of stress under dynamic loads and focusing on capturing the characteristics of stress concentration areas such as the mandibular angle, condyle, and alveolar ridge. The instantaneous deformation field includes fields such as timestamp, global mesh node number, three-dimensional displacement value, and displacement gradient value, reflecting the instantaneous deformation characteristics under the combined action of mandibular movement and masticatory forces. Global mesh nodes refer to all uniquely numbered hard and soft tissue mesh nodes of the craniofacial region in the individualized three-dimensional digital twin model at each standard timestamp in the discrete-time-load table; these are the smallest units for the finite element solver to calculate the mechanical response. In one specific embodiment, the duration is calculated with time steps of 0.01 seconds and from 0 seconds to 180 seconds to fully capture at least one mandibular opening and closing and occlusal contact cycle.

[0102] After the solution is completed, the global mesh nodes of the individualized 3D digital twin model are classified and labeled for craniofacial skeletal hard tissue and facial soft tissue. Only the skeletal hard tissue mesh nodes are extracted for subsequent derivation of the bone remodeling rate field.

[0103] The bone micro-remodeling rate field is derived only for hard bone tissue; facial soft tissue is not included in the calculation. Based on temporal stress field distribution and instantaneous deformation field, and combined with bone metabolism biomechanical theory, the derivation is performed step-by-step. The specific steps are as follows:

[0104] Based on the three-dimensional displacement values ​​in the instantaneous deformation field, the equivalent stress values ​​of bone hard tissue mesh nodes whose deformation exceeds the physiological limit index are eliminated. The physiological limit index is set according to the physiological range of motion of the craniofacial bones. In one specific embodiment, the physiological limit index is set as follows: the three-dimensional displacement value of the mandibular mesh node exceeds 10 mm, or the abrupt change in displacement gradient of adjacent alveolar bone mesh nodes is greater than 0.5 mm per step.

[0105] Grouped by grid node number, all valid equivalent stress values ​​for each bone hard tissue grid node over the entire time series are extracted. The time points are then weighted according to the chewing cycle period to calculate the average equivalent stress value for each bone hard tissue grid node. The calculation formula is as follows:

[0106]

[0107] in, For the first The average equivalent stress value of each bone hard tissue mesh node, in megapascals; For time points The corresponding weighting coefficient is set based on the load percentage at that point in the chewing cycle. For the first Each bone hard tissue mesh node at time point The equivalent stress value, in megapascals; This represents the total number of time steps.

[0108] Based on the average equivalent stress value and the preset bone remodeling stress threshold, the stress zoning type of each bone hard tissue mesh node is determined. The preset bone remodeling stress threshold is set based on the physiological remodeling characteristics of the craniofacial skeleton. In a specific embodiment, the bone remodeling stress zoning thresholds for each anatomical location are shown in Table 1.

[0109] Table 1 Bone remodeling stress zoning thresholds for each anatomical site

[0110]

[0111] Based on the stress partition type and average equivalent stress value of each bone hard tissue mesh node, a pre-defined nonlinear bone remodeling rate model is used to calculate the bone micro remodeling rate node by node.

[0112] In one specific embodiment, the preset nonlinear bone remodeling rate model adopts a nonlinear logarithmic function, which fits the nonlinear saturation relationship between bone remodeling rate and stress, and its specific form is as follows:

[0113]

[0114] in, The rate of bone microremodeling is expressed in millimeters per year. The remodeling rate coefficient for the tension deposition zone can be set at a base value of 0.05 mm / year to 0.1 mm / year, and can be finely adjusted based on individual age and craniofacial growth pattern. The remodeling rate coefficient for the pressure absorption zone can be set at a base value of 0.06 mm / year to 0.12 mm / year, and can be finely adjusted based on individual age and craniofacial growth pattern. This represents the lower limit threshold of the equivalent stress in the tension deposition zone, in megapascals (MPa). The lower limit threshold of equivalent stress in the pressure absorption zone is given in megapascals. The equivalent stress value is the physiological limit, and the unit is megapascal. It is set based on the mechanical tolerance limit of the craniofacial skeleton. It is the natural logarithm function.

[0115] The data is correlated with the node number, 3D displacement coordinates, average equivalent stress value, stress partition type, and bone microremodeling rate of each bone hard tissue mesh node to form a complete dataset. Based on the spatial coordinate system of the individualized 3D digital twin model, this dataset is mapped to the global mesh nodes to ensure a one-to-one correspondence between bone microremodeling rates and anatomical locations. Finally, it is converted into a bone microremodeling rate field in VTK format, with fields including mesh node coordinates, mesh node number, remodeling rate scalar field, stress partition type, and average equivalent stress value. This format is consistent with the temporal stress field distribution and instantaneous deformation field formats and can be directly used for coupled computational units.

[0116] The coupled computational unit is used to regularize and correct the initial growth deformation field by taking the temporal stress field distribution, instantaneous deformation field and bone micro-remodeling rate field as physical constraints, and generate the final growth deformation prediction result.

[0117] The coupled computation unit is configured as follows:

[0118] The initial growth deformation field, temporal stress field distribution, instantaneous deformation field, and bone micro-remodeling rate field were normalized.

[0119] In one specific embodiment, for the initial growth deformation field, the three-dimensional displacement value of each grid node in the initial growth deformation field is divided by the displacement normalization reference value to obtain the normalized initial growth deformation field; for the instantaneous deformation field, the three-dimensional displacement value is divided by the displacement normalization reference value to obtain the normalized instantaneous deformation field; for the temporal stress field distribution, the equivalent stress value of each bone hard tissue grid node is divided by the stress normalization reference value to obtain the normalized temporal stress field distribution; for the bone microremodeling rate field, the bone microremodeling rate value of each grid node is divided by the bone remodeling rate normalization reference value to obtain the normalized bone microremodeling rate field. Simultaneously, the spatial coordinates of all grid nodes in the individualized 3D digital twin model are also normalized by dividing by the displacement normalization reference value, and are only used for intermediate calculations; the actual geometric mesh of the individualized 3D digital twin model remains unchanged. The time axis of the growth prediction is normalized by dividing by the time normalization reference value.

[0120] Among them, the displacement normalization reference value is the global maximum amplitude of the three-dimensional displacement value of growth, which is recommended to be set to 15.0 mm; the stress normalization reference value is the peak equivalent stress, which is recommended to be set to 50.0 MPa; the bone remodeling rate normalization reference value is the maximum absolute value of the scalar value of the bone remodeling rate, which is recommended to be set to 2.0 mm / year; and the time normalization reference value is recommended to be set to 1 year.

[0121] Based on the normalized temporal stress field distribution, instantaneous deformation field, skeletal micro-remodeling rate field, and initial growth deformation field, a regularized objective function incorporating spatially adaptive physical constraints is constructed:

[0122]

[0123] in, The value of the regularization objective function; The normalized and corrected growth deformation field to be solved; This represents the initial growth deformation field after normalization. To correspond with the normalized equivalent stress distribution of the growth deformation field; The normalized temporal stress field distribution; The spatially adaptive weight matrix is ​​set based on clinical anatomical data and biomechanical experimental results; To normalize and correct the Laplace operator of the growth deformation field; This is the stress constraint coefficient; The coefficient is the spatial smoothing constraint coefficient. The deformation gradient constraint coefficient; The coefficient constraining the rate of bone remodeling; To normalize and correct the displacement gradient values ​​of the growth deformation field; The spatial gradient of the instantaneous deformation field after normalization; To normalize the first derivative of the growth deformation field with respect to time; The normalized skeletal microremodeling rate field; This is element-wise multiplication; It is the Euclidean norm.

[0124] Specifically, the spatial adaptive weight matrix This is a one-dimensional row matrix with the same number of elements as the number of skeletal hard tissue mesh nodes. Each element is bound to a single skeletal hard tissue mesh node, and the values ​​are assigned based on clinical anatomical data and biomechanical experimental results, ranging from 0.2 to 1.0. In a specific embodiment, the core areas for masticatory force transmission, such as the mandibular angle, condyle, alveolar ridge, and mandibular ramus, as well as key areas for craniofacial growth, are set as high-weight areas, with values ​​ranging from 0.8 to 1.0; the secondary core areas for stress or growth, such as the maxillary body, zygomatic bone, and mandibular body, are set as medium-weight areas, with values ​​ranging from 0.5 to 0.7; and the areas of resting stress, such as the cranial parietal bone and frontal bone, as well as areas with slow growth, are set as low-weight areas, with values ​​ranging from 0.2 to 0.4.

[0125] Stress constraint coefficient Spatial smoothing constraint coefficient Deformation gradient constraint coefficient and bone remodeling rate constraint coefficient The value ranges from 0 to 1. The basic preset value is set based on the clinical needs of craniofacial growth prediction and can be fine-tuned according to individual circumstances. In one specific embodiment, the stress constraint coefficient... Set to 0.8, spatial smoothing constraint coefficient Set to 0.5, deformation gradient constraint coefficient Set to 0.6, the bone remodeling rate constraint coefficient. Set it to 1.0.

[0126] As a basic term, the deviation between the corrected deformation field and the initial deformation field is quantified to avoid over-correction and loss of individualized growth trend; As a stress constraint term, the stress constraint in the key area is strengthened so that the stress distribution of the modified deformation field conforms to the dynamic load mechanics law; This is a spatial smoothing constraint term that constrains the spatial continuity of the deformation field, preventing non-physiological facial deformities caused by abrupt displacement changes between mesh nodes. As a deformation gradient constraint term, the modified deformation gradient matches the instantaneous deformation characteristics of the dynamic load, ensuring the coordination between growth deformation and physiological movement; The bone remodeling rate constraint term ensures that the modified growth deformation rate conforms to the biomechanical laws of bone metabolism, and is the core link connecting mechanical simulation and growth prediction.

[0127] The regularized objective function is iteratively optimized using the conjugate gradient method. When the change in deformation field between adjacent iterations is less than the preset deformation threshold, the normalized corrected deformation field is obtained.

[0128] Specifically, the classic conjugate gradient method is used to iteratively solve the problem by minimizing the objective function value. In each iteration, based on the normalized correction growth deformation field of the current iteration step, the derived quantities are calculated as follows:

[0129] The displacement gradient value of the normalized modified growth deformation field is obtained by performing spatial differentiation on the normalized modified growth deformation field; the Laplace operator of the normalized modified growth deformation field is obtained by performing second-order spatial differentiation on the normalized modified growth deformation field; the calculation method of the normalized equivalent stress distribution is consistent with the calculation method of the equivalent stress value in the above-mentioned time-series stress field distribution, that is, based on the normalized modified growth deformation field of the current iteration step, the normalized stress tensor is calculated through the strain-displacement relationship and linear elastic constitutive relationship of the individualized three-dimensional digital twin model, and then converted by the equivalent stress formula.

[0130] The iteration termination condition is: the maximum change in the 3D displacement value of the global mesh nodes of the individualized 3D digital twin model in adjacent iteration steps is less than a preset deformation threshold. This threshold is set based on the physiological accuracy of the craniofacial bone micro-remodeling and the mesh resolution of the individualized 3D digital twin model. It is a general basic value and does not need to be fine-tuned according to individual circumstances. The recommended value is 0.001 mm. After the iteration terminates, the normalized corrected deformation field is obtained.

[0131] The normalized modified deformation field is denormalized to output the final growth deformation prediction result.

[0132] Specifically, the final growth deformation prediction results include a deformation field file in VTK format and a quantitative analysis report. The VTK deformation field file contains corrected 3D displacement values ​​for all mesh nodes, which can be directly loaded into an individualized 3D digital twin model for visualization, intuitively presenting the growth deformation trends of various regions of the craniofacial region. The quantitative analysis report includes quantitative data on the growth amount and rate of key anatomical regions such as the mandibular angle, chin apex, nasolabial angle, and alveolar ridge, as well as predicted growth deformation values ​​at different time scales based on a preset time series, providing quantitative evidence for clinicians to develop orthodontic and orthognathic intervention plans.

[0133] In one specific embodiment, after the regularized iterative solution is completed, the coupled computation unit multiplies the corrected three-dimensional displacement value in the obtained normalized deformation field by the displacement normalization reference value and performs inverse normalization to restore its true physical scale, in millimeters. The displacement normalization reference value uses the same global maximum amplitude of the grown three-dimensional displacement value as the normalization value, and is recommended to be set to 15.0 millimeters.

[0134] The visualization module is used to visualize the deformation process of an individualized 3D digital twin model over a preset time series based on the growth deformation prediction results.

[0135] Specifically, the preset time series offers two settings: equal time interval and clinical growth stage. The zero point of time is uniformly set to the time of patient clinical data collection. The equal time interval can be set to a step size of 1 month, 6 months, or 1 year. The clinical growth stage matches key clinical growth nodes such as the primary dentition period, mixed dentition period, adolescent growth peak period, and adult stable period.

[0136] The visualization module is implemented based on 3D graphics rendering engines such as OpenGL, DirectX, or Unreal Engine. The specific technical implementation process is as follows:

[0137] Load the static mesh data of the individualized 3D digital twin model, including vertex coordinates, triangle patch indices and normal vectors, and construct the basic rendering environment such as 3D scene, camera, and light source.

[0138] For each preset time point, the vertex displacement field is loaded into the graphics processing unit cache; in the vertex shader, the static vertex coordinates are superimposed with the displacement field to generate the deformed vertex coordinates, thereby realizing real-time deformation rendering of the mesh.

[0139] The stress field distribution or growth remodeling intensity data is normalized to the 0 to 1 range; a predefined color lookup table is used, for example, using a red-yellow-green-blue gradient, where red corresponds to high stress or high remodeling intensity areas and blue corresponds to low stress or low remodeling intensity areas; in the fragment shader, the color lookup table is indexed according to the stress value of each vertex to generate surface color, thereby realizing the visualization of stress distribution and growth remodeling intensity.

[0140] A time slider is set in the user interface, with the slider position corresponding to a preset time point index. When the user drags the slider, the vertex shader is triggered to load the displacement field of the corresponding time point, updating the mesh deformation in real time. It also supports a continuous animation playback mode, loading data of each time point in sequence according to a set frame rate to realize dynamic simulation of the growth process.

[0141] Based on a virtual camera model, the camera's position, target, and up vector can be adjusted via mouse or touch events to perform rotation, scaling, and translation operations. By defining a clipping plane, fragments located on one side of the clipping plane can be discarded in the fragment shader to achieve cross-sectional cutting at any position. In addition, individualized 3D digital twin model data from two time points can be loaded simultaneously, and before-and-after comparisons can be achieved through semi-transparent rendering or side-by-side display. Predictions before and after treatment can be loaded, and differences can be displayed through color differentiation or semi-transparent overlay.

[0142] It supports split-screen display mode and simultaneously presents multiple visualization formats such as 3D deformation view, stress heat map, and time series curve, providing users with multi-dimensional growth analysis tools.

[0143] The interaction and feedback module is used to trigger the hybrid growth prediction module to regenerate the growth deformation prediction results based on the user's input treatment intervention parameters when the user inputs them. The regenerated growth deformation prediction results are then updated in real time to the visualization module for display.

[0144] Specifically, users can input or select treatment intervention parameters such as appliance type, traction direction, treatment duration, and functional intervention method through the interactive interface of the visualization module, either by drop-down selection or manual input. The interaction and feedback module encapsulates these parameters into intervention loads and boundary conditions according to the craniofacial clinical biomechanics conversion rules. The converted format is consistent with the temporal load spectrum of the mechanism-driven simulation unit. The clinical biomechanics conversion rules specifically include: traction direction corresponding to load direction, traction force corresponding to load magnitude, and appliance wearing area corresponding to load application node. The encapsulated intervention loads and boundary conditions are then sent to the hybrid growth prediction module, triggering the system to regenerate growth deformation prediction results based on the treatment intervention parameters, and simultaneously pushing them to the visualization module. The visualization module then loads the new deformation field file, realizing real-time dynamic updates of the prediction results, allowing users to intuitively observe the craniofacial growth deformation trends under different treatment intervention parameters.

[0145] The training data augmentation module is used to use generative adversarial networks to fill in missing image data between adjacent acquisition time points in the collected historical patient clinical follow-up data, and to construct a continuous longitudinal time sequence for training deep learning models.

[0146] Specifically, historical patient clinical follow-up data is retrieved from the clinical follow-up database. This data includes a time series of data collected at different growth stages and time points, including at least cone-beam CT data, three-dimensional facial surface scan data, clinical phenotype data, and genetic information data. All data are classified and labeled according to patient ID, age, gender, and craniofacial growth type.

[0147] A Wasserstein generative adversarial network with gradient penalty is used as the core network architecture to learn the nonlinear temporal evolution of craniofacial growth. The Adam optimizer is employed with a learning rate of 0.0002, a batch size of 8, and 500 training epochs. The loss function is WGAN-GP loss to avoid mode collapse.

[0148] Early and late-stage multimodal data from each patient in a historical patient population are input into a generative adversarial network (GAN). Through adversarial training, this network learns craniofacial growth patterns across different ages, sexes, growth types, and genetic backgrounds. Based on these patterns, it fills in missing craniofacial morphology, bone structure, and soft tissue contour data between adjacent data collection points. The completed data is then fused with the original clinical follow-up data to construct a continuous, smooth, and physiologically consistent longitudinal time-series sequence. This sequence maintains the same format and spatial coordinates as the actual clinical follow-up data and is stored by patient ID and age. This constructed continuous longitudinal time-series sequence is used as a training dataset to train the deep learning model in the hybrid growth prediction module, expanding the training sample size and improving the model's prediction accuracy, robustness, and generalization ability.

[0149] like Figure 3 As shown, another embodiment of the present invention provides a method for predicting the growth of a craniofacial digital twin. This method is based on the aforementioned craniofacial digital twin growth prediction system and includes:

[0150] S1. Acquire multimodal craniofacial data and dynamic biomechanical time-series data under physiological activity states. Among them, the multimodal craniofacial data includes at least cone-beam computed tomography (CBCT) data, three-dimensional facial surface scan data, clinical phenotype data, and genetic information data.

[0151] S2 preprocesses multimodal craniofacial data and dynamic biomechanical time-series data, and registers the preprocessed cone-beam CT data with the preprocessed three-dimensional facial surface scan data.

[0152] S3 constructs an individualized three-dimensional digital twin model containing bones, jawbone, and facial soft tissues based on processed cone-beam CT data and processed three-dimensional facial surface scan data.

[0153] S4, based on an individualized three-dimensional digital twin model, and combined with processed clinical phenotypic data, processed dynamic biomechanical time-series data, and processed genetic information data, adopts a data-driven and mechanism simulation coupled approach to simulate the cumulative growth effect under dynamic load cycles and generate growth deformation prediction results for the craniofacial region.

[0154] S5, based on the growth deformation prediction results, visualizes the deformation process of the individualized 3D digital twin model over a preset time series.

[0155] In summary, this invention proposes a craniofacial skeletal digital twin growth prediction system and method. Through the full-domain fusion of multimodal data, the construction of individualized 3D digital twins, hybrid prediction combining data-driven and mechanistic simulation, dynamic functional simulation, and closed-loop interactive optimization, it achieves precise individualized dynamic prediction of craniofacial growth trajectory and quantitative prediction of treatment intervention effects. This improves the scientific rigor of clinical treatment plan design, the rationality of treatment timing selection, and the reliability of prognostic assessment. Specifically, the data acquisition module and data processing module work collaboratively to comprehensively capture static anatomical data, dynamic biomechanical time-series data, genetic information data, and clinical phenotypic data of the craniofacial region, achieving full-dimensional information coverage of morphology, function, and genetics. Through cross-modal registration techniques such as timestamp association, differentiated preprocessing, and deep learning feature point matching, scattered data is transformed into a standardized dataset with unified format, spatiotemporal alignment, and logical association. This fundamentally solves the core pain points of existing technologies, such as single data dimension, spatiotemporal misalignment, modeling difficulties, and limited prediction accuracy. The digital twin construction module adopts a model-through-the-prediction design, constructing an individualized 3D digital twin model that includes bones, jawbone, and facial soft tissues, and possesses a mechanical coupling relationship between hard and soft tissues. This design completely eliminates the shortcomings of existing technologies that rely on prediction followed by reconstruction, avoiding the disconnect between prediction results and the model, and the information loss introduced by secondary reconstruction. It allows prediction results to be directly presented in the form of vertex displacement fields of the individualized 3D digital twin model, achieving a precise depiction from the average growth pattern of the population to the unique growth trajectory of each individual, ensuring lossless information transmission and higher prediction accuracy. Building upon this, the hybrid growth prediction module employs an innovative architecture driven by dynamic load coupling and both physics and data. First, dynamic biomechanical time-series data is transformed into a time-series load spectrum, which serves as the load node set applied as dynamic boundary conditions to an individualized 3D digital twin model. This realistically recreates the long-term cumulative impact of daily functional activities such as chewing and speech on craniofacial growth and development. Then, an initial growth deformation field is output through a deep learning model using a hybrid architecture of 3DCNN and Transformer. This is combined with the time-series stress field and instantaneous deformation field obtained from finite element simulation, as well as the bone micro-remodeling rate field derived from bone metabolism biomechanics theory. A regularized objective function with spatially adaptive physical constraints is then constructed to correct the initial growth deformation field through multi-dimensional physical constraints. This architecture ensures that the prediction results simultaneously satisfy statistical data laws and biomechanical conservation laws, successfully transforming the prediction model from a black box into a physically interpretable white box. This guarantees both prediction speed and individualized accuracy while ensuring that the results conform to physiological laws, significantly improving the accuracy and reliability of the predictions.The visualization module dynamically presents the craniofacial growth and remodeling process based on a preset time series, providing an efficient and intuitive support tool for clinical diagnosis, treatment design, and doctor-patient communication. The interaction and feedback module allows doctors to input treatment intervention parameters, which are encapsulated into intervention loads and boundary conditions through clinical biomechanical transformation rules. This triggers the hybrid growth prediction module to recalculate in real time, achieving a closed loop from parameter adjustment to effect simulation and treatment optimization, enabling the system to continuously learn and evolve in clinical applications. Furthermore, the training data augmentation module constructs continuous longitudinal time series through generative adversarial networks, effectively addressing the industry challenges of scarce longitudinal medical data and sparse time points, providing ample high-quality samples for training deep learning models.

[0156] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various variations or substitutions within the technical scope disclosed in this application, and these should all be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A craniofacial skeleton digital twin growth prediction system, characterized in that, include: The data acquisition module is used to collect multimodal craniofacial data and dynamic biomechanical time-series data under physiological activity states; The multimodal craniofacial data includes at least cone-beam CT data, three-dimensional facial surface scan data, clinical phenotype data, and genetic information data; The data processing module is used to preprocess multimodal craniofacial data and dynamic biomechanical time-series data, and to register the preprocessed cone-beam CT data with the preprocessed three-dimensional facial surface scan data. The digital twin construction module is used to construct an individualized three-dimensional digital twin model containing bones, jawbone, and facial soft tissue based on processed cone-beam CT data and processed three-dimensional facial surface scan data. The hybrid growth prediction module is used to simulate the cumulative growth effect under dynamic load cycles based on an individualized three-dimensional digital twin model, combined with processed clinical phenotypic data, processed dynamic biomechanical time series data and processed genetic information data. It adopts a data-driven and mechanism simulation coupled approach to generate growth deformation prediction results of the craniofacial region. The visualization module is used to visualize the deformation process of an individualized 3D digital twin model over a preset time series based on the growth deformation prediction results.

2. The craniofacial skeleton digital twin growth prediction system according to claim 1, characterized in that, The data acquisition module is equipped with an oral cavity flexible sensor and / or a mandibular motion capture device; The clinical phenotypic data includes tooth arrangement characteristics, occlusal parameters, and craniofacial growth type classification data. The genetic information data includes genome-wide association analysis data or multigene risk score data; The dynamic biomechanical time series data includes chewing force time series and mandibular movement trajectory time series.

3. The craniofacial skeleton digital twin growth prediction system according to claim 1, characterized in that, The data processing module includes: The data preprocessing unit is used to preprocess multimodal craniofacial data and dynamic biomechanical time-series data; the preprocessing includes outlier detection, noise reduction and format conversion; The cross-modal registration unit is used to unify the preprocessed cone-beam CT data and the preprocessed 3D facial surface scan data into the same spatial coordinate system based on the deep learning feature point matching algorithm, thus completing the registration.

4. The craniofacial skeleton digital twin growth prediction system according to claim 1, characterized in that, The digital twin building module is configured to perform the following steps: Based on the processed cone-beam CT data, surface mesh models of the skull and jawbone are reconstructed as initial static geometric models of the skeleton. The morphological contour, elastic modulus, and Poisson's ratio parameters of facial soft tissue were obtained from the processed three-dimensional facial surface scanning data. A static geometric model of facial soft tissue was constructed based on the attachment relationship between soft and hard tissues. The mechanical coupling relationship between the initial static geometric model of the skeleton and the static geometric model of the facial soft tissue was established through finite element analysis. Based on the mechanical coupling relationship, the initial static geometric model of the skeleton is fused with the static geometric model of the facial soft tissue to construct an individualized three-dimensional digital twin model containing the skeleton, jawbone and facial soft tissue.

5. The craniofacial skeleton digital twin growth prediction system according to claim 1, characterized in that, The hybrid growth prediction module includes: The data-driven prediction unit is used to input the individualized three-dimensional digital twin model, processed clinical phenotypic data, and processed genetic information data into the deep learning model and output the initial growth deformation field. The mechanism-driven simulation unit is used to construct a time-series load spectrum from the processed dynamic biomechanical time-series data. The time-series load spectrum is applied as a boundary condition to the individualized three-dimensional digital twin model. The model is solved by a finite element solver to obtain the time-series stress field distribution and instantaneous deformation field. Based on the solution results and the biomechanical theory of bone metabolism, the bone micro-remodeling rate field is obtained. The coupled computational unit is used to regularize and correct the initial growth deformation field by taking the temporal stress field distribution, instantaneous deformation field and bone micro-remodeling rate field as physical constraints, and generate the final growth deformation prediction result.

6. The craniofacial skeleton digital twin growth prediction system according to claim 5, characterized in that, The instantaneous deformation field includes at least the three-dimensional displacement values ​​of the global mesh nodes of the individualized three-dimensional digital twin model; the global mesh nodes include skeletal hard tissue mesh nodes and facial soft tissue mesh nodes; The time-series stress field distribution includes at least the equivalent stress values ​​of the bone hard tissue mesh nodes.

7. A craniofacial skeleton digital twin growth prediction system according to claim 6, characterized in that, The mechanism-driven simulation unit obtains the bone micro-remodeling rate field and is configured to perform the following steps: Based on the three-dimensional displacement values ​​in the instantaneous deformation field, the equivalent stress values ​​of bone hard tissue mesh nodes whose deformation exceeds the physiological limit index are eliminated; the physiological limit index is set according to the physiological range of craniofacial bones. Extract all remaining equivalent stress values ​​of each bone hard tissue mesh node throughout the entire time series, and weight the time points according to the chewing cycle to calculate the average equivalent stress value of each bone hard tissue mesh node; Based on the average equivalent stress value and the preset bone remodeling stress threshold, the stress partition type of each bone hard tissue mesh node is determined. The preset bone remodeling stress threshold is set based on the physiological remodeling characteristics of the craniofacial skeleton. Based on the stress partitioning type and average equivalent stress value of each bone hard tissue mesh node, a preset nonlinear bone remodeling rate model is used to calculate the bone micro remodeling rate. The number, three-dimensional displacement coordinates, average equivalent stress value, stress partition type, and bone micro-remodeling rate of each bone hard tissue mesh node are associated and mapped to the global mesh nodes of the individualized three-dimensional digital twin model to obtain the bone micro-remodeling rate field.

8. The craniofacial skeleton digital twin growth prediction system according to claim 5, characterized in that, The coupled computation unit is configured to perform the following steps: The initial growth deformation field, temporal stress field distribution, instantaneous deformation field, and bone micro-remodeling rate field are normalized. Based on the normalized temporal stress field distribution, instantaneous deformation field, initial growth deformation field, and bone micro-remodeling rate field, a regularized objective function incorporating spatially adaptive physical constraints is constructed: ； in, The value of the regularization objective function; The normalized and corrected growth deformation field to be solved; This represents the initial growth deformation field after normalization. To correspond with the normalized equivalent stress distribution of the growth deformation field; The normalized temporal stress field distribution; The spatially adaptive weight matrix is ​​set based on clinical anatomical data and biomechanical experimental results; To normalize and correct the Laplace operator of the growth deformation field; This is the stress constraint coefficient; The coefficient is the spatial smoothing constraint coefficient. The deformation gradient constraint coefficient; The coefficient constraining the rate of bone remodeling; To normalize and correct the displacement gradient values ​​of the growth deformation field; This represents the displacement gradient value of the instantaneous deformation field after normalization. To normalize the first derivative of the growth deformation field with respect to time; The normalized skeletal microremodeling rate field; This is element-wise multiplication; It is the Euclidean norm; The regularization objective function is iteratively optimized using the conjugate gradient method. When the change in deformation field between adjacent iterations is less than the preset deformation threshold, the normalized corrected deformation field is obtained. The normalized modified deformation field is denormalized to output the final growth deformation prediction result.

9. A craniofacial skeleton digital twin growth prediction system according to claim 5, characterized in that, Also includes: The interaction and feedback module is used to trigger the hybrid growth prediction module to regenerate the growth deformation prediction results based on the user's input treatment intervention parameters when the user inputs treatment intervention parameters, and to update the regenerated growth deformation prediction results to the visualization module in real time for display. The training data augmentation module is used to use generative adversarial networks to complete the missing image data between adjacent acquisition time points in the collected historical patient clinical follow-up data, and to construct a continuous longitudinal time sequence for training the deep learning model in the hybrid growth prediction module.

10. A method for predicting the growth of a craniofacial digital twin, based on the craniofacial digital twin growth prediction system according to any one of claims 1 to 9, characterized in that, include: Collect multimodal craniofacial data and dynamic biomechanical time-series data under physiological activity states; The multimodal craniofacial data includes at least cone-beam CT data, three-dimensional facial surface scan data, clinical phenotype data, and genetic information data; Multimodal craniofacial data and dynamic biomechanical time-series data were preprocessed, and the preprocessed cone-beam CT data were registered with the preprocessed three-dimensional facial surface scan data. Based on the processed cone-beam CT data and the processed three-dimensional facial surface scan data, an individualized three-dimensional digital twin model containing bones, jawbone and facial soft tissues was constructed. Based on a personalized three-dimensional digital twin model, and combined with processed clinical phenotypic data, processed dynamic biomechanical time series data and processed genetic information data, a data-driven and mechanism simulation coupled approach is used to simulate the cumulative growth effect under dynamic load cycle and generate growth deformation prediction results for the craniofacial region. Based on the growth deformation prediction results, the deformation process of the individualized 3D digital twin model over a preset time series is visualized.

Images