Bone drilling force prediction method based on CT and drilling working conditions

Abstract

The invention discloses a bone drilling force prediction method based on CT and drilling working conditions. The bone drilling force prediction method comprises the following steps that 1, a contact surface image sequence of a drill bit and bone tissue on a bone drilling path is extracted from a CT image and preprocessed; step 2, inputting the preprocessed contact surface image and drilling working condition parameters into a multi-physical factor fusion module to obtain an image sequence fused with working condition physical information; 3, inputting the image fused with the physical information of the working condition into an image coding network, and extracting drilling force characteristics in a tangent plane; and step 4, inputting a sequence formed by the drilling force characteristics of all sections on the drilling path into the sequence regression network ATP-UNeXt to obtain prediction of a drill bit stress sequence on the drilling path. The method has the advantages that physical factors related to the drilling force can be comprehensively considered, more accurate preoperative stress prediction is provided, the model reasoning speed is high, the method is suitable for a clinical operation auxiliary system, and the cost of the operation auxiliary system is reduced.

Description

Technical Field

[0001] The present invention relates to the technical field of drilling force prediction in orthopedic surgery, and in particular to a bone drilling force prediction method based on CT (Computed Tomography) and drilling conditions. Background Art

[0002] In orthopedic surgery, precise control and perception of the interaction between surgical tools and bone tissue is the key to achieving precise surgery. Traditional surgical robot systems mainly rely on real-time force sensors to monitor force signals during surgery, but this method has limitations, such as the inability to predict force changes during surgery in advance, which causes the robot to react slowly when faced with complex surgical tasks and is unable to analyze and predict based on preoperative images and real-time perception information like human surgeons.

[0003] Existing physical model-based drilling force prediction methods mainly include analytical methods, voxel methods and finite element model-based methods. Although these methods can simulate and predict drilling forces to a certain extent, their computational efficiency is low and they cannot fully consider the drilling conditions and various physical properties of bone tissue.

[0004] In recent years, the development of artificial intelligence technology has brought new solutions to the field of surgical robots. Neural network models have attracted attention due to their superior approximation ability and high computing speed of the actual drilling physics process at the data level. However, the existing bone drilling force prediction methods based on deep learning only model the relationship between working condition parameters and drilling force, and cannot learn the effects of different bone tissue morphology and material properties on drilling force. In addition, the existing methods are still unable to unify the feature extraction of drilling condition information and bone tissue morphology density information reflected by the image in the neural network, resulting in impaired force prediction accuracy. Summary of the invention

[0005] In order to overcome the shortcomings of the above-mentioned prior art, the present invention provides a bone drilling force prediction method based on CT and drilling conditions. This method incorporates the drilling conditions and bone tissue images into a bone drilling force prediction network based on deep learning, more comprehensively considers the physical factors related to the drilling force, and provides a more accurate preoperative force prediction for bone drilling.

[0006] The present invention adopts the following technical solutions to solve the above problems:

[0007] The present invention provides a bone drilling force prediction method based on CT and drilling conditions, and the specific steps are as follows:

[0008] Step 1: Extract and preprocess the image sequence of the contact surface between the drill bit and the bone tissue on the bone drilling path from the CT image;

[0009] Step 2: Input the preprocessed contact surface image and drilling condition parameters into the multi-physical factor fusion module to obtain an image sequence of fusion condition physical information;

[0010] Step 3: Input the image integrating the physical information of the working condition into the image coding network to extract the drilling force characteristics in the section;

[0011] Step 4: Input the sequence of drilling force features of all sections on the drilling path into the sequence regression network ATP-UNeXt to obtain the prediction of the drill force sequence on the drilling path;

[0012] The training data of the multi-physical factor fusion module and the image encoding network come from the three-dimensional images generated by the simulation program and the forces calculated by the physical model based on the extrusion and cutting effects; the training data of the sequence regression network ATP-UNeXt comes from the CT images scanned in the real bone drilling experiment and the force sensor signals.

[0013] The step 1 comprises:

[0014] Step 1-1, perform MicroCT scanning on the bone tissue before the bone drilling operation and set drilling parameters such as feed rate F, rotation speed S, and drill diameter D; during the training stage, the bone tissue needs to be scanned again after the operation to align the preoperative image with the postoperative image.

[0015] Step 1-2, extracting the contact surface image sequence V between the drill bit and the bone tissue from the preoperative image 1 ,I 2 ,…,I T}. Where T is the number of images extracted from the drilling path, I t ∈R H*W It represents the t-th slice image in the image sequence V, where H and W are the height and width of the image respectively.

[0016] Step 1-3, for each slice image I in the image sequence V t , respectively, perform MinMax normalization processing, as shown in the following formula: Where x is the pixel gray value in the image before normalization, x max With x min are the maximum and minimum values ​​of the pixel grayscale in the image, respectively, and x′ is the pixel grayscale value in the normalized image.

[0017] Steps 1-4: Convert the working condition parameters into the slice image I t Multi-channel 2D matrices of the same size Among them C M is the number of types of operating parameters.

[0018] The step 2 comprises:

[0019] Step 2-1, use the multi-layer gated linear unit GMLP to generate a spatial weight map based on the drilling condition. The process is expressed by the following formula:

[0020] M C =GMLP(M)

[0021] in, is the spatial weight diagram based on drilling conditions, C N is the number of channels of the weight map.

[0022] Step 2-2, weight map M C Each image I in the contact surface image sequence t Perform pixel-by-pixel dot product operations to generate an image sequence that integrates the physical information of the working condition. The process is expressed by the following formula:

[0023]

[0024] in, V is the image after point multiplication and fusion of physical information of the working condition. C For all the drilling paths The spatial weight map is generated by GMLP, and the working condition parameters (such as feed rate and speed) are fused with the CT image pixel by pixel, which solves the limitation of existing methods that only model working conditions or image features independently.

[0025] The image encoding network Enc in step 3 S The cross-sectional image that integrates the physical information of the working condition Converted into drilling force characteristic u in the cutting surface t , the process is expressed by the following formula:

[0026]

[0027] The image coding network in step 3 uses multiple groups of ConvNeXt Blocks for feature extraction. Each ConvNeXt Block consists of a convolution layer, layer normalization, and a GELU activation function. PixelUnShuffle operations are used between each group of ConvNeXt Blocks to upsample the image dimension from C in *H in *W in Transform to C out *H out *W out ,satisfy:

[0028] C out =C in×downsample_factor 2

[0029] H out =H in ÷downsample_factor

[0030] W out =W in ÷downsample_factor

[0031] Among them, downsample_factor is the downsampling factor. Through simulation data pre-training (based on physical models) and experimental data fine-tuning, we can not only use physical prior knowledge to make the model more in line with the actual physical process, but also verify the robustness of the model in real scenarios and improve the generalization ability of the model.

[0032] The overall data processing process of the sequence regression network in step 4 is expressed by the following formula:

[0033]

[0034] Among them, UNX is a sequence regression network, U={u 1 ,u 2 ,…,u T} is the drilling force characteristic u of all cut surfaces on the drilling path t The sequence composed of is the network’s prediction sequence of the drill bit force sequence on the drilling path, where is the force prediction when the drill bit is at the t-th section position, C f is the number of drilling forces to be predicted.

[0035] The sequence regression network in step 4 adopts the UNet architecture suitable for processing multi-channel one-dimensional data, including a multi-layer encoder to extract features of the sequence at different scales. The process is expressed by the following formula:

[0036]

[0037] Among them, Encoder i is the i-th layer encoder, N is the number of layers, is the input of the i-th layer encoder, is the output of the i-th layer encoder.

[0038] Correspondingly, the sequence regression network also contains a multi-layer decoder that receives the jump connection from the encoder at the same level and the decoder output of the previous level as output, and gradually decodes the deep semantic features into a force sequence prediction of the same length as the original input sequence. This process is expressed by the following formula:

[0039]

[0040] Among them, Decoder i is the i-th layer decoder, is the output of the i-th layer decoder, X bottle is the output of the bottleneck layer, and Res is the residual block on the skip connection path.

[0041] The sequence regression network in step 4 introduces a module ATP-Block with a multi-frequency feature extraction and fusion mechanism to semantically fuse the shallow and deep feature sequences in the encoder under different receptive fields. The process is expressed by the following formula:

[0042]

[0043] Among them, W b With b b is the parameter of the fully connected layer, Conv Di is a convolution block with a dilation rate of Di, n is the number of different dilation rates, is the output of the encoder at the Nth layer. The complete force sequence on the drilling path is directly output, avoiding the cumulative error of point-by-point prediction in the traditional single-step prediction method, and improving the overall prediction efficiency by more than 30%.

[0044] Compared with the prior art, the present invention has the following beneficial effects:

[0045] By combining physical models and deep learning technology, the present invention can more accurately predict bone drilling force from CT images and surgical plans. Compared with a single physical model or deep learning model, the prediction accuracy is improved.

[0046] In the network structure designed by the present invention, the multi-physical factor fusion module and the image encoder based on the cutting physical model can learn the physical features within the section that conforms to the physical model, and the sequence regression network can supplement the temporal features between the sections, thereby improving the robustness of the model.

[0047] The model reasoning speed of the present invention is significantly improved compared to the calculation speed of the physical model, and is more suitable for clinical surgery assistance systems.

[0048] The prediction of the present invention is based on preoperative images and bone drilling operation status information, and does not rely on real-time force sensor signals during surgery, thereby reducing the cost of the surgical auxiliary system. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate embodiments consistent with the present application and, together with the description, serve to explain the principles of the present application.

[0050] Figure 1 is a flow chart of a bone drilling force prediction method based on CT and drilling conditions shown in an embodiment of the present invention,

[0051] Figure 2 is a structural diagram of a multi-physical factor fusion module in an embodiment of the present invention,

[0052] Figure 3 is a structural diagram of an image encoder in an embodiment of the present invention,

[0053] Figure 4 is a structural diagram of a sequence regression network in an embodiment of the present invention,

[0054] Figure 5 1 is a set of comparative experimental result diagrams in the examples of the present invention. DETAILED DESCRIPTION

[0055] The present invention is further explained below in conjunction with specific embodiments. It should be understood that these embodiments are only used to illustrate the present invention and are not used to limit the scope of the present invention. After reading the present invention, various equivalent forms of modifications to the present invention by those skilled in the art all fall within the scope defined by the claims attached to this application.

[0056] Embodiment: The bone drilling force prediction method based on CT and drilling conditions described in the embodiment of the present invention is referred to Figure 1 The steps include:

[0057] Step 1: Extract and preprocess the image sequence of the contact surface between the drill bit and the bone tissue on the bone drilling path from the CT image;

[0058] Step 1-1, perform MicroCT scanning on the bone tissue before the bone drilling operation and set drilling parameters such as feed rate F, rotation speed S, and drill diameter D; during the training stage, the bone tissue needs to be scanned again after the operation to align the preoperative image with the postoperative image.

[0059] Step 1-2, extracting the contact surface image sequence V between the drill bit and the bone tissue from the preoperative image 1 ,I 2 ,…,I T}. Where T is the number of images extracted from the drilling path, I t ∈R H*W It represents the t-th slice image in the image sequence V, where H and W are the height and width of the image respectively.

[0060] Step 1-3, for each slice image I in the image sequence V t , respectively, perform MinMax normalization processing, as shown in the following formula: Where x is the pixel gray value in the image before normalization, x maxWith x min are the maximum and minimum values ​​of the pixel grayscale in the image, respectively, and x′ is the pixel grayscale value in the normalized image.

[0061] Steps 1-4: Convert the working condition parameters into the slice image I t Multi-channel 2D matrices of the same size Among them C M is the number of types of operating parameters.

[0062] Step 2: Input the preprocessed contact surface image and drilling condition parameters into the multi-physical factor fusion module to obtain an image sequence of fusion condition physical information;

[0063] Step 2-1, refer to Figure 2 As shown in the figure, a multi-layer gated linear unit GMLP is used to generate a spatial weight map based on the drilling condition to model the influence weight of the physical factors contained in the drilling condition on the cutting micro-element. The process is expressed by the following formula:

[0064] M C =GMLP(M)

[0065] in, is the spatial weight diagram based on drilling conditions, C N is the number of channels of the weight map.

[0066] Step 2-2, weight map M C Each image I in the contact surface image sequence t Perform pixel-by-pixel dot product operations to generate an image sequence that integrates the physical information of the working condition. The process is expressed by the following formula:

[0067]

[0068] in, V is the image after point multiplication and fusion of physical information of the working condition. C For all the drilling paths A sequence of images.

[0069] Step 3: Input the image that integrates the physical information of the working condition into the image coding network to extract the drilling force characteristics within the section. Figure 3 As shown, the image encoding network Enc s The cross-sectional image that integrates the physical information of the working condition Converted into drilling force characteristic u in the cutting surface t , the process is expressed by the following formula:

[0070]

[0071] The image coding network in step 3 uses multiple groups of ConvNeXt Blocks for feature extraction. Each ConvNeXt Block consists of a convolution layer, layer normalization, and a GELU activation function. PixelUnShuffle operations are used between each group of ConvNeXt Blocks to upsample the image dimension from C in *H in *W in Transform to C out *H out *W out ,satisfy:

[0072] C out =C in ×downsample_factor 2

[0073] H out =H in ÷downsample_factor

[0074] W out =W in ÷downsample_factor

[0075] Among them, downsample_factor is the downsampling ratio.

[0076] After being processed by N groups of ConvNeXt Blocks, the feature map is transformed into a feature vector u using the global maximum pooling operation. t As the full path drilling force characteristic sequence U = {u 1 ,u 2 ,…,u T}. In the pre-training stage, the FC (Fully Connected) layer is also used to convert u t The transformation is used to calculate the prediction of drilling force for the physical model, so that the image encoder can fit the cutting physics model.

[0077] Step 4: Input the sequence of drilling force characteristics of all sections on the drilling path into the sequence regression network to obtain the prediction of the drill force sequence on the drilling path. The overall data processing process of this step is expressed by the following formula:

[0078]

[0079] Among them, UNX is a sequence regression network, U={u 1 ,u 2 ,…,u T} is the drilling force characteristic u of all cut surfaces on the drilling path t The sequence composed of is the network’s prediction sequence of the drill bit force sequence on the drilling path, where is the force prediction when the drill bit is at the t-th section position, C f is the number of drilling forces to be predicted.

[0080] Reference Figure 4 As shown in (a), the sequence regression network in step 4 adopts the UNet architecture suitable for processing multi-channel one-dimensional data, including a multi-layer encoder to extract features of the sequence at different scales. The process is expressed by the following formula:

[0081]

[0082] Among them, Encoder i is the i-th layer encoder, N is the number of layers, is the input of the i-th layer encoder, is the output of the i-th layer encoder.

[0083] Correspondingly, the sequence regression network also contains a multi-layer decoder that receives the jump connection from the encoder at the same level and the decoder output of the previous level as output, and gradually decodes the deep semantic features into a force sequence prediction of the same length as the original input sequence. This process is expressed by the following formula:

[0084]

[0085] Among them, Decoder i is the i-th layer decoder, is the output of the i-th layer decoder, X bottle is the output of the bottleneck layer, and Res is the residual block on the skip connection path.

[0086] Reference Figure 4 As shown in (b), the sequence regression network in step 4 introduces a deep separation convolution module DWCBlock, which reduces the amount of calculation and parameters while maintaining a large receptive field. The calculation process is expressed by the following formula:

[0087] X D_out =FC 2 (GELU(FC 1 (LN((DC(X D_in )))))))+X D_in

[0088] Among them, X D_in is the input data of DWC Block, X D_outis the output data of DWC Block, DC is the depth convolution kernel, LN is the layer normalization operation, FC 1 , FC 2 They are the fully connected layers for channel expansion and contraction respectively.

[0089] Reference Figure 4 As shown in (c), the sequence regression network in step 4 introduces a module ATP-Block with a multi-frequency feature extraction and fusion mechanism to semantically fuse the shallow and deep feature sequences in the encoder under different receptive fields. The process is expressed by the following formula:

[0090]

[0091] Among them, W b With b b is the parameter of the fully connected layer, Conv Di is a convolution block with a dilation rate of Di, n is the number of different dilation rates, is the output of the Nth layer encoder.

[0092] The training data of the multi-physical factor fusion module and the image encoding network come from the three-dimensional images generated by the simulation program and the forces calculated by the physical model based on the extrusion and cutting effects; the training data of the sequence regression network ATP-UNeXt comes from the CT images scanned in the real bone drilling experiment and the force sensor signals.

[0093] In the training and testing of the sequence regression network, a data set containing 14 drilling paths from 3 cow bones was constructed, of which 10 paths from 2 bones were used as training and validation sets for 10-fold cross validation, and 4 paths from the remaining bones were used as test sets. During the training process, SmoothL1 was used as the loss function, and the calculation formula was:

[0094]

[0095] Among them, y is the true value, is the model's predicted value.

[0096] Choose to use MAE, MAPE, R 2 As an indicator for evaluating model performance, the calculation formula is:

[0097]

[0098] Where N is the number of samples in the loss function calculation batch.

[0099]

[0100] in, Compute the average of the true values ​​in the batch for the loss function.

[0101] The parameter settings during model training are shown in Table 1.

[0102] Table 1 Experimental hyperparameter settings

[0103]

[0104]

[0105] In this example, the method proposed in the present invention is compared with the existing method through experiments. Table 2 shows the MAE and MAPE indicators of the sequence regression network proposed in the present invention, the existing cutting physical model and other deep learning models on the test set. The results show that the prediction accuracy of ATP-UNeXt is significantly higher than that of the physical model, and it also has the lowest MAE and MAPE error indicators among all the deep learning models in the experiment.

[0106] Table 2 Model performance comparison

[0107]

[0108] Regarding the system in the above embodiment, the specific steps for each module to perform the operation have been described in detail in the embodiment of the method, and will not be elaborated in detail here. The various modules of the above-mentioned multi-stage detection method of CT spinal bone lesions based on deep learning can be implemented in whole or in part by software, hardware and a combination thereof. The above-mentioned modules can be embedded in or independent of the processor in the computer device in the form of hardware, or can be stored in the memory in the computer device in the form of software, so that the processor can call and execute the operations corresponding to the above modules.

[0109] In addition, each functional unit in each embodiment of the present application may be integrated into a processing module, or each unit may exist physically separately, or two or more units may be integrated into one module. The above-mentioned integrated module may be implemented in the form of hardware or in the form of a software functional module. If the integrated module is implemented in the form of a software functional module and sold or used as an independent product, it may also be stored in a computer-readable storage medium.

[0110] The protection content of the present invention is not limited to the above embodiments. Without departing from the spirit and scope of the inventive concept, changes and advantages that can be thought of by those skilled in the art are included in the present invention and are protected by the attached claims.

Claims

1. A bone drilling force prediction method based on CT and drilling conditions, characterized in that: The method comprises the following steps: Step 1: Extract and preprocess the image sequence of the contact surface between the drill bit and the bone tissue on the bone drilling path from the CT image; Step 2: Input the preprocessed contact surface image and drilling condition parameters into the multi-physical factor fusion module to obtain an image sequence of fusion condition physical information; Step 3: Input the image integrating the physical information of the working condition into the image encoder to extract the drilling force characteristics in the section; Step 4: Input the sequence of drilling force characteristics of all sections on the drilling path into the sequence regression network ATP (Atrous Temporal Pyramid)-UNeXt to obtain the prediction of the drill force sequence on the drilling path; The training data of the multi-physical factor fusion module and the image encoding network come from the three-dimensional images generated by the simulation program and the forces calculated by the physical model based on the extrusion and cutting effects; the training data of the sequence regression network ATP-UNeXt comes from the CT images scanned in the real bone drilling experiment and the force sensor signals.

2. The bone drilling force prediction method based on CT and drilling conditions according to claim 1, characterized in that: Step 1 includes: Step 1-1, before the bone drilling operation, perform a MicroCT scan on the bone tissue and set the drilling parameters of feed rate F, rotation speed S, and drill diameter D; during the training phase, the bone tissue needs to be scanned again after the operation to align the preoperative image with the postoperative image. Step 1-2, extracting the contact surface image sequence V between the drill bit and the bone tissue from the preoperative image 1 ,I 2 ,…,I T }, where T is the number of images extracted from the drilling path, I t ∈R H*W represents the t-th slice image in the image sequence V, where H and W are the height and width of the image respectively. Step 1-3, for each slice image I in the image sequence V t , respectively, perform MinMax normalization processing, as shown in the following formula: Where x is the pixel gray value in the image before normalization, x max With x min are the maximum and minimum values ​​of the pixel grayscale in the image, respectively, and x′ is the pixel grayscale value in the normalized image. Steps 1-4: Convert the working condition parameters into the slice image I t Multi-channel 2D matrices of the same size Among them C M is the number of types of operating parameters.

3. The bone drilling force prediction method based on CT and drilling conditions according to claim 1, characterized in that: The multi-physical factor fusion module in step 2 includes a multi-layer gated linear unit GMLP (Gated-MultilayerPerceptron) to generate a spatial weight map based on drilling conditions. C N is the number of channels of the weight map; weight map M C Each image I in the contact surface image sequence t Perform pixel-by-pixel dot product operations to generate an image sequence that integrates the physical information of the working condition. The process is expressed by the following formula: in, V is the image after point multiplication and fusion of physical information of the working condition. C For all the drilling paths A sequence of images.

4. The bone drilling force prediction method based on CT and drilling conditions according to claim 3, characterized in that: The image encoder in step 3 includes a convolutional neural network with a ConvNeXt architecture, which uses the PixelUnShuffle pixel reorganization operation instead of the convolution pooling operation for downsampling, and encodes the image that integrates the physical information of the working condition into the drilling force characteristics within the section; the network uses a staged training strategy: first pre-training with simulation data, and then fine-tuning with experimental measurement data.

5. The bone drilling force prediction method based on CT and drilling conditions according to claim 4, characterized in that: The image encoding network Enc in step 3 S The cross-sectional image that integrates the physical information of the working condition Converted into drilling force characteristic u in the cutting surface t , the process is expressed by the following formula: The image encoding network in step 3 uses multiple groups of ConvNeXt Blocks for feature extraction. Each ConvNeXtBlock consists of a convolution layer, layer normalization, and a GELU activation function. PixelUnShuffle operations are used between each group of ConvNeXtBlocks to upsample the image dimension from C in *H in *W in Transform to C out *H out *W out ,satisfy: C out =C in ×downsample_factor 2 H out =H in ÷downsample_factor W out =W in ÷downsample_factor Among them, downsample_factor is the downsampling ratio.

6. The bone drilling force prediction method based on CT and drilling conditions according to claim 5, characterized in that: In step 4, the sequence regression network ATP-UNeXt takes the image sequence that integrates the physical information of the working condition as input to model the temporal characteristics between sections; the network uses multi-frequency feature extraction modules in the stem layer and the bottleneck layer, and uses depthwise separable convolution kernels in the feature extraction module.

7. The bone drilling force prediction method based on CT and drilling conditions according to claim 6, characterized in that: The overall data processing process of the sequence regression network in step 4 is expressed by the following formula: Among them, UNX is a sequence regression network, U={u 1 ,u 2 ,…,u T } is the drilling force characteristic u of all cut surfaces on the drilling path t The sequence composed of is the network’s prediction sequence of the drill bit force sequence on the drilling path, where is the force prediction when the drill bit is at the t-th section position, C f is the number of drilling forces to be predicted, The sequence regression network in step 4 adopts the UNet architecture suitable for processing multi-channel one-dimensional data, including a multi-layer encoder to extract features of the sequence at different scales. The process is expressed by the following formula: Among them, Encoder i is the i-th layer encoder, N is the number of layers, is the input of the i-th layer encoder, is the output of the i-th layer encoder; The sequence regression network also contains a multi-layer decoder that receives the jump connection from the encoder at the same level and the decoder output of the previous level as output, gradually decoding the deep semantic features into a force sequence prediction of the same length as the original input sequence. The process is expressed by the following formula: Among them, Decoder i is the i-th layer decoder, is the output of the i-th layer decoder, X bottle is the output of the bottleneck layer, and Res is the residual block on the skip connection path.

8. The bone drilling force prediction method based on CT and drilling conditions according to claim 7, characterized in that: The sequence regression network in step 4 introduces a module ATP-Block with a multi-frequency feature extraction and fusion mechanism to semantically fuse the shallow and deep feature sequences in the encoder under different receptive fields. The process is expressed by the following formula: Among them, X bottle is the output of the ATP-Block module, W b With b b is the parameter of the fully connected layer, Conv Di is a convolution block with a dilation rate of Di, n is the number of different dilation rates, is the output of the Nth layer encoder.

Images