A bone tissue processing assistance planning system, method and storage medium based on cumulative thermal dose constraints

Abstract

The application discloses a kind of bone tissue processing auxiliary planning system, method and storage medium based on cumulative heat dose constraint.The system is based on preoperative medical image segmentation to be removed region and heat-sensitive tissue safety neighborhood, generates individualized bone thermal physical parameter field according to CT value, CBCT gray correction value or multi-modal fusion bone density estimation value and bone structure direction;Candidate processing path is discretized into multiple path segments, mobile heat source is established and transient heat diffusion model is solved;According to the temperature field at the end of the previous path segment and cooling condition, obtain the initial residual temperature field of the current path segment, calculate the voxel temperature curve and incremental heat dose and monotonically accumulate;Under the condition that maximum cumulative heat dose and super-threshold voxel proportion meet the constraint, the order of path segment, instrument parameter and cooling strategy are optimized, and the recommended path and heat dose risk map are output.The application combines geometric collision avoidance with heat dose avoidance, based on individualized thermal physical parameter field of patient and path segment recursive heat dose, with intraoperative feedback correction, can provide auxiliary planning under thermal safety constraint for bone tissue processing before and during operation.

Description

Technical Field

[0001] This invention relates to the fields of medical image processing, computer-aided planning, surgical navigation, and thermal safety control in bone tissue processing, and particularly to a bone tissue processing auxiliary planning system, method, and storage medium based on cumulative thermal dose constraints. Background Technology

[0002] In surgical planning involving bone grinding, drilling, osteotomy, scraping, or milling, existing computer-aided planning methods typically focus primarily on the spatial relationships between anatomical structures. For example, they consider whether instruments will mechanically collide with nerves, blood vessels, or other vital tissues, whether the path is short enough, whether the access facilitates exposure, and whether the area to be removed can be adequately covered. While these methods can improve geometric collision avoidance, their safety margins largely depend on the geometric distance between instruments and vital tissues.

[0003] However, during bone tissue processing, the friction and cutting between instruments such as drills, oscillating saws, and milling cutters and the bone tissue generate heat. This heat can diffuse along the bone cortex, trabeculae, structures surrounding the medullary cavity, and areas of local bone defects. Even if the instruments are not geometrically in direct contact with heat-sensitive tissues such as nerves, blood vessels, the spinal cord, and inner ear structures, there is still a risk of exceeding the predicted heat dose threshold due to short-term high temperatures or long-term low-level heat accumulation. Such risks are insidious and delayed, and cannot be identified in a timely manner by relying solely on mechanical collision avoidance planning.

[0004] Existing temperature monitoring methods are mostly used for intraoperative local temperature alarms or to evaluate the effects of parameters such as rotational speed, feed rate, instrument diameter, and irrigation flow rate on temperature rise in offline experiments. These methods typically cannot compare the voxel-level cumulative heat dose formed by different candidate paths, different processing sequences, and different combinations of instrument parameters during the path planning stage, nor can they dynamically update the remaining path based on the patient's individualized bone density, bone structure orientation, and intraoperative measured temperature.

[0005] Therefore, there is a need for a bone tissue processing auxiliary planning scheme that can unify individualized bone thermal conduction characteristics of patients, heat sources of moving instruments, transient thermal field prediction, voxel-level cumulative thermal dose evaluation, sequential optimization of path segments, and intraoperative feedback correction. Summary of the Invention

[0006] The technical problem this invention aims to solve is that existing bone tissue processing planning is mainly based on geometric collision avoidance, which cannot predict the cumulative heat dose generated in the safe neighborhood of heat-sensitive tissue after heat is conducted along the bone tissue during instrument processing. In particular, it cannot dynamically adjust the sequence of subsequent path segments, instrument parameters and cooling strategies according to the patient's individual bone density, bone structure direction and intraoperative measured temperature.

[0007] To address the aforementioned technical problems, this invention provides a bone tissue processing auxiliary planning system based on cumulative heat dose constraints. This system includes at least a medical image acquisition module, a tissue segmentation module, a bone thermophysical parameter field generation module, a candidate path generation module, a moving heat source modeling module, a transient thermal field solving module, a cumulative heat dose recursion module, a multi-objective path optimization module, an intraoperative temperature registration and model calibration module, and a path / parameter output module.

[0008] The medical image acquisition module is used to acquire preoperative CT, CBCT, MRI, or fused images. The tissue segmentation module is used to obtain the region to be removed, the thermally sensitive tissue region, and the safe neighborhood of the thermally sensitive tissue. The thermally sensitive tissue includes at least one of nerves, blood vessels, spinal cord, nerve roots, dura mater, inner ear structures, inferior alveolar canal, mucosa, or skin. The safe neighborhood is obtained by calculating the distance field, morphological dilation, or risk expansion based on tissue-specific dose thresholds from the segmented thermally sensitive tissue region. Its thickness or dose threshold is determined by physician configuration, anatomical structure category configuration, literature thresholds, or experimental calibration results.

[0009] The bone thermal property parameter field generation module is used to map bone density based on CT values, CBCT grayscale correction values, calibration phantom conversion results, or multimodal fusion bone density estimates, and to generate a patient-specific bone tissue thermal property parameter field by combining bone structure orientation tensors. The thermal property parameter field of bone tissue includes at least anisotropic thermal conductivity tensor K(x), which is expressed as: K(x)=R(A(x))·diag(k_parallel(rho(x)), k_perp1(rho(x)), k_perp2(rho(x)))·R(A(x))^T; where rho(x) is the bone density obtained by mapping or estimation from CT value HU(x), CBCT grayscale correction value, calibration phantom conversion result or multimodal fusion image, A(x) is the bone structure orientation tensor, R(A(x)) is the rotation matrix determined by A(x), and k_parallel, k_perp1 and k_perp2 are the thermal conductivity along the principal direction and the two transverse directions, respectively. The bone structure orientation tensor A(x) is obtained from the gray-level gradient structure tensor, the principal direction of bone trabecular texture, the mean intercept length analysis, local high-resolution images, or anatomical orientation templates. When the confidence level of the bone structure orientation is lower than the preset confidence threshold, the system adopts an isotropic or orthotropic simplified thermophysical property model.

[0010] The candidate path generation module generates multiple candidate bone tissue processing paths and discretizes these paths into multiple path segments. The moving heat source modeling module establishes a moving heat source term for each path segment based on instrument position, processing time, rotational speed, feed rate, and dwell time. The heat source intensity is jointly estimated by at least two of the following: rotational speed, torque, axial force, feed rate, instrument power, instrument diameter, tool wear state, and cooling state. The moving heat source term is a spatially distributed heat source, including at least one of the following: Gaussian distributed moving heat source, spherical contact body heat source, cylindrical contact body heat source, or surface heat source along the instrument contact surface. In one embodiment, the device power is calculated as P_tool(t) = 2π·rpm(t)·Torque(t) / 60, and the moving heat source term is calculated as Q_i(x,t) = η_i·P_tool(t)·G(∥x-p_i(t)∥,σ_i), where G is the integrally normalized spatial heat source core, and η_i is the heat proportion coefficient entering the bone tissue, calibrated through in vitro bone experiments or artificial bone experiments.

[0011] The transient thermal field solution module is used to solve the transient thermal diffusion model during bone tissue processing to obtain the temperature-time curves of multiple voxels within the safe neighborhood of heat-sensitive tissue. The transient thermal diffusion model satisfies: ρ_b(x)c_b(x)·∂T(x,t) / ∂t=∇·(K(x)∇T(x,t))+Q_i(x,t)-H_cool(x,t); where ρ_b(x) is the equivalent density of bone tissue, c_b(x) is the specific heat capacity, K(x) is the thermal conductivity tensor, Qi(x,t) is the moving heat source term corresponding to the i-th path segment, and H_cool(x,t) is the cooling term. The cooling term H_cool(x,t) is determined based on one or more of the following: flushing flow rate, flushing fluid temperature, cooling pause time, convective heat transfer coefficient of exposed surface, and negative pressure attraction state. The system uses at least one of the following methods to solve or update the transient thermal diffusion model: finite element method, finite difference method, finite volume method, path segment thermal response kernel pre-computation, reduced-order model, lookup table method, or GPU parallel computing.

[0012] The cumulative heat dose recursion module is used to determine the initial residual temperature field of the current path segment based on the end temperature field at the end of the previous path segment and the cooling pause time, and uses this initial residual temperature field as the initial condition to calculate the incremental cumulative heat dose of the current path segment. The incremental cumulative heat dose can be calculated using CEM43 heat dose: DeltaDose_i(x)=∫{t_i}^{t{i+1}} R^{43-T_i(x,t)} dt; where T_i(x,t) is the predicted temperature of the i-th path segment at position x, in °C; dt is in minutes; R is the temperature conversion coefficient, configured according to the temperature range, tissue type, or experimental calibration results. The incremental cumulative heat dose can also be calculated using the Arrhenius thermal damage integral: Omega_i(x)=∫{t_i}^{t{i+1}} A_0·exp[-E_a / (R_g·T_K(x,t))] dt; where A_0 is the frequency factor, E_a is the activation energy, R_g is the gas constant, and T_K(x,t) is the Kelvin temperature. The cumulative heat dose recursive module monotonically accumulates the incremental cumulative heat dose for each path segment according to the candidate processing order, obtaining the cumulative heat dose index corresponding to the candidate processing order. This index includes the maximum cumulative heat dose Dose_max and the proportion of voxels exceeding the dose threshold V_dose. The specific recursive process is as follows: Based on the end temperature field T_end,i-1(x) at the end of the previous path segment and the cooling pause time Delta t_i, solve the heat diffusion model without a heat source or with a cooling term to obtain the initial residual temperature field T_start,i(x) of the current path segment; use T_start,i(x) as the initial condition to solve the temperature-time curve of the current path segment and calculate DeltaDose_i(x); the cumulative heat dose is updated according to Dose_acc_i(x)=Dose_acc_{i-1}(x)+DeltaDose_i(x). The cooling pause is only used to reduce or change the initial residual temperature field of subsequent path segments and is not used to attenuate the already formed CEM43 heat dose or Arrhenius thermal damage integral. When using the Arrhenius thermal damage integral, the system updates the cumulative thermal damage integral according to Omega_acc,i(x)=Omega_acc,i-1(x)+Omega_i(x).

[0013] The multi-objective path optimization module is used to jointly optimize the path segment sequence, rotational speed, feed rate, dwell time, and cooling pause time under the constraints that Dose_max does not exceed the maximum dose threshold and V_dose does not exceed the proportional threshold. When Dose_max exceeds the maximum dose threshold or V_dose exceeds the proportional threshold under the candidate processing sequence, the system performs at least one of the following actions: adjusts the path segment sequence, reduces the rotational speed, reduces the heat source power, shortens the dwell time, changes the feed rate, inserts a cooling pause, adds flushing cooling, or removes the candidate path segment. The joint optimization is also performed under the conditions that the coverage of the area to be removed is not less than the coverage threshold, the minimum distance between the instrument and the prohibited structure is not less than the safety distance, the instrument collision cost is lower than the collision threshold, the mechanical accessibility meets the preset accessibility conditions, and the processing residue cost is lower than the residue threshold.

[0014] The intraoperative temperature registration and model calibration module receives intraoperative measurement data such as infrared thermography, temperature probes, instrument power, torque, axial force, instrument pose, or irrigation flow rate, and registers them to the preoperative medical imaging coordinate system. The system performs model calibration based on the error term e(t) = T_surface_measured(t) - T_surface_predicted(t). When e(t) exceeds the temperature error threshold, Dose_max exceeds the maximum dose threshold, V_dose exceeds the proportional threshold, or there is a sudden change in instrument power / torque, the system updates the heat source intensity coefficient, cooling heat transfer coefficient, thermal conductivity correction coefficient, or boundary conditions using least squares, recursive least squares, extended Kalman filtering, or Bayesian updates. After calibration, the system re-predicts the temperature-time curve and cumulative heat dose index for the remaining path segment and outputs the adjusted auxiliary planning results.

[0015] The recommended paths, instrument parameters, cooling strategies, and thermal dose risk maps output by this invention are used to assist medical professionals, surgical navigation systems, or robot-assisted systems in decision-making. They are only used as display, prompt, alarm, or parameter recommendation information outputs and do not directly replace the judgment of medical professionals. They do not automatically generate closed-loop control commands that directly drive the instrument to act on human tissues, nor do they directly limit the implementation of treatment actions on the human body.

[0016] The beneficial effects of this invention are as follows: First, it not only considers the geometric distance between the instrument and tissues such as nerves and blood vessels, but also uses the voxel-level cumulative heat dose caused by heat conduction along bone tissue as a path planning constraint, expanding bone tissue processing planning from simple mechanical collision avoidance to an auxiliary planning that combines mechanical collision avoidance and heat dose avoidance. Second, it generates an individualized thermophysical parameter field based on the patient's individual bone density, bone structure orientation, and local bone tissue type, avoiding the simplistic treatment of different bone tissues as homogeneous materials, thereby improving the specificity of thermal field prediction and dose risk assessment. Third, it discretizes candidate paths into path segments and uses the residual temperature field after cooling of the previous path segment as the initial condition for the next path segment, monotonically accumulating the incremental heat dose for each path segment, which can simultaneously handle the residual temperature field, cooling recovery, and cumulative heat dose constraints in continuous processing. Fourth, it combines intraoperative temperature, power, torque, or force feedback for model correction, enabling the planning results to adapt to instrument status, cooling conditions, and operational deviations, forming an auxiliary planning mechanism that combines preoperative planning with intraoperative feedback updates. Fifth, in the simulation comparison examples, compared with the scheme that only considers geometric collision avoidance, the Dose_max of the present invention is reduced from 45.0 min to 7.9 min, and V_dose is reduced from 34.2% to 3.7%. After extended Kalman filter correction, the temperature prediction MAE and RMSE are reduced by 71.2% and 70.6% respectively, while the remaining path recalculation time remains in the second range. Attached Figure Description

[0017] Figure 1 A schematic diagram of the module structure of the bone tissue processing auxiliary planning system provided in an embodiment of the present invention.

[0018] Figure 2 This is a schematic diagram of the recursive optimization process for incremental heat dose in path segments provided in an embodiment of the present invention.

[0019] Figure 3 This is a schematic diagram of the intraoperative temperature feedback correction process provided in an embodiment of the present invention.

[0020] Figure 4 This is a schematic diagram illustrating the relationship between the safe neighborhood of a heat-sensitive tissue and candidate path segments, provided in an embodiment of the present invention.

[0021] Figure 5 The CT-HU to bone mineral density calibration curve and lookup relationship diagram provided in the embodiments of the present invention.

[0022] Figure 6 A comparison chart of different thermophysical parameters of bone tissue provided for embodiments of the present invention.

[0023] Figure 7 The instrument temperature rise distribution and moving heat source parameter calibration diagram provided in the embodiments of the present invention.

[0024] Figure 8 The flushing cooling temperature reduction curve and the convective heat transfer coefficient calibration diagram are provided for embodiments of the present invention.

[0025] Figure 9 This is a diagram showing the tissue-specific CEM43 threshold and safe neighborhood configuration provided in an embodiment of the present invention.

[0026] Figure 10 This is a comparison chart of the computational performance of different thermal field solving methods provided in the embodiments of the present invention.

[0027] Figure 11 The temperature prediction and extended Kalman filter correction results are shown in the figure provided for the embodiments of the present invention.

[0028] Figure 12 The simulation comparison results of the four planning methods provided in the embodiments of the present invention are shown in the figure.

[0029] Figure 13 A comparison chart of parameters for four typical bone tissue processing scenarios provided in the embodiments of the present invention. Detailed Implementation

[0030] The present invention will be further described below with reference to embodiments. It should be understood that the following embodiments are used to illustrate the technical solutions of the present invention and are not intended to limit the scope of protection of the present invention. In the absence of conflict, the technical features of each embodiment can be combined with each other. The following calibration values, threshold configurations and verification results are parameterized, simulated or replay embodiments constructed based on medical engineering literature parameters and reasonable assumptions. In actual applications, they can be recalibrated according to the target tissue, device model, imaging equipment, cooling conditions and measured data.

[0031] Existing geometric collision avoidance path planning primarily considers the distance and collision relationship between instruments and critical tissues, typically neglecting to calculate voxel-level thermal dose within the safe neighborhood of heat-sensitive tissues. While existing intraoperative temperature alarms can alert to local temperature anomalies, they struggle to compare the order of multiple candidate path segments during preoperative or intraoperative planning. Current studies on bone drilling parameters often focus on the impact of rotational speed, drill diameter, feed rate, or flushing flow rate on temperature rise, but do not necessarily output path segment order, instrument parameters, and cooling strategies based on Dose_max and V_dose constraints. The difference in this invention lies in: under a patient-individualized thermophysical parameter field, recursively applying the residual temperature field as a cross-path segment state variable, and monotonically accumulating the incremental thermal dose of each path segment, thereby jointly optimizing the path segment order, instrument parameters, and cooling strategy.

[0032] The system first estimates bone mineral density (rho(x)) based on preoperative medical images. Bone mineral density estimation can be obtained from CT values ​​(HU(x)), CBCT grayscale correction values, calibration phantom conversion results, grayscale mapping between registered CT and CBCT images, or multimodal fusion images. This mapping can employ linear calibration curves, piecewise calibration curves, lookup tables, or device-specific conversion relationships obtained from the calibration phantom. For CBCT or other grayscale-unstable images, grayscale normalization, region calibration, or registration with preoperative CT images can be performed before mapping. In one CT-HU to bone mineral density calibration example, a QRM-BDC bone mineral density calibration phantom and a Siemens SOMATOM Force CT were used. Repeated scans were performed under conditions of 120 kVp, 250 mAs, a slice thickness of 0.6 mm, and a B70f convolution kernel. Linear fitting yielded ρ(x) = 0.000868 × HU(x) + 0.9346, with a coefficient of determination R² of 0.9850 and a root mean square error of 0.0751 g / cm³. Piecewise fitting was used for the low HU and high HU regions, respectively. Calibration conditions and standard material lookup points are shown in Tables 1-1 and 1-2, and calibration curves are shown in... Figure 5 .

[0033] Table 1-1 CT-HU to Bone Mineral Mineral Degradation Conditions and Fitting Results

[0034] Body Model Number QRM-BDC (QRM Micro-CT Bone Density Calibration Phantom) CT equipment Siemens SOMATOM Force Scan parameters 120 kVp, 250 mAs, 0.6 mm layer thickness, B70f convolution kernel Sample size n=40 (5 standard materials × 8 repetitions) Linear fitting formula ρ(x) = 0.000868 × HU(x) + 0.9346 Coefficient of determination R² 0.9850 Root Mean Square Error (RMSE) 0.0751 g / cm³ Piecewise fitting (low HU < 400) ρ = 0.000972 × HU + 0.9706 Piecewise fitting (high HU≥400) ρ = 0.000603 × HU + 1.1847

[0035] Table 1-2 Standard Material Density and CT Value Lookup Points

[0036] Air 0.0012 -1000 water 1.0000 0 soft tissue 1.0600 40 cancellous bone 1.1500 250 cortical bone 1.9200 1200

[0037] For regions with clearly defined trabecular bone or trabecular bone orientation, the system determines the bone structure orientation tensor A(x) based on local gray-level gradient, structural tensor, mean intercept length analysis, principal texture direction, or anatomical template. When the confidence level of the bone structure orientation is insufficient, the system can degenerate into an isotropic or orthotropic simplified model. In a preferred embodiment, the thermal conductivity tensor is expressed as: K(x) = R(A(x))·diag(k_parallel(rho(x)), k_perp1(rho(x)), k_perp2(rho(x)))·R(A(x))^T. The equivalent density, specific heat capacity, and anisotropic thermal conductivity of bone tissue can be initialized according to the ranges in Table 2. For vaporized bone and diseased bone, the system preferably performs secondary correction of the parameters by combining local imaging features and measured temperature rise; when the confidence level of the orientation tensor is insufficient, conservative values ​​or isotropic simplified values ​​within the range can be used. A comparison of various parameter types is shown in [reference needed]. Figure 6 .

[0038] Table 2. Range of thermophysical parameters for different bone tissues

[0039] cortical bone 1.80–2.00 1200–1400 1.80–2.50 0.50–1.00 0.40–0.90 1.8–5.0 Davies 1986; Mori 2005 cancellous bone 0.30–1.20 1100–1300 0.40–1.20 0.30–0.80 0.25–0.70 0.5–4.0 Craciun 2021 Qi transformation bone 0.15–0.60 900–1200 0.10–0.50 0.08–0.30 0.06–0.25 0.3–6.2 Torimitsu 2015 Diseased bone 0.20–1.00 1000–1300 0.20–0.80 0.15–0.50 0.12–0.40 0.4–5.3 In vivo experiment

[0040] Where R(A(x)) is the local principal direction rotation matrix, and k_parallel, k_perp1 and k_perp2 are determined by a preset calibration curve or lookup table.

[0041] The system discretizes the candidate bone tissue processing path into path segments {s_i}. Each path segment has a corresponding instrument position p_i(t), instrument radius r_i, processing time tau_i, rotational speed rpm_i, feed rate feed_i, and dwell time dwell_i, which are used to establish the moving heat source term; cooling pause time pause_i, flushing flow rate, and exposed surface convective heat transfer coefficient are used to determine the cooling term or residual temperature field update conditions. In one embodiment, the moving heat source term Q_i(x,t) is a Gaussian distributed heat source: Q_i(x,t)=η_i·P_tool_i(t)·G(∥x-p_i(t)∥,σ_i). When the instrument rotation speed rpm(t) and torque Torque(t) are obtained, P_tool(t) = 2π·rpm(t)·Torque(t) / 60; when the axial force F_axial(t) and feed speed feed(t) are obtained, P_tool(t) can also be estimated by F_axial(t)·feed(t) or their weighted combination. G(∥x-p_i(t)∥,σ_i) is the spatial heat source core and is normalized. η_i is the heat proportion coefficient entering the bone tissue, calibrated by in vitro bone experiments or artificial bone experiments. In one instrument heat source calibration embodiment, five segments of isolated porcine femur were used, and K-type thermocouples with a diameter of 0.5 mm and a embedment depth of 2 mm were used to simultaneously measure the temperature with a FLIR A700 infrared thermal imager, with a total number of observations n=300. The heat proportion coefficient η_i and the heat source core width σ_i of different instruments were fitted, and the example results are shown in Table 3 and Figure 7 The parameters in the table are initial values ​​specific to the instrument model and can be corrected online based on rotational speed, torque, axial force, feed rate, tool wear, and cooling status.

[0042] Table 3. Calibration results of mobile heat source parameters for different instruments

[0043] Diamond grinding head φ2.0mm 2.0 0.35 1.2 101.6±84.3 R²=0.315 Diamond grinding head φ3.0mm 3.0 0.38 1.5 97.6±70.4 R²=0.459 Carbide grinding head φ2.3mm 2.3 0.42 1.0 184.4±198.9 R²=0.596 Oscillating saw N / A 0.30 2.0 2.0±0.0 — End mill φ4.0mm 4.0 0.32 2.5 51.1±40.3 R²=0.595

[0044] Wherein, P_tool_i(t) is the instrument input power or equivalent cutting power, which can be estimated from the instrument power, torque and angular velocity, axial force and feed rate, or by fusing the above data.

[0045] The system establishes a transient heat diffusion equation for the bone tissue region: ρ_b(x)c_b(x)·∂T(x,t) / ∂t=∇·(K(x)∇T(x,t))+Q_i(x,t)-H_cool(x,t). Boundary conditions include convective heat transfer on the exposed bone surface, cooling of the irrigation fluid, instrument contact boundary, and heat exchange boundary with surrounding soft tissue. H_cool(x,t) is determined based on irrigation flow rate, irrigation fluid temperature, cooling interval, pause time, suction state, and convective heat transfer coefficient of the exposed surface. In one irrigation-cooling calibration embodiment, eight segments of excised porcine ribs were used, with an exposed area of ​​approximately 30 mm², an initial temperature of 80 ℃, an ambient temperature of 37 ℃, and 0.9% physiological saline as the irrigation fluid. The convective heat transfer coefficient was estimated using h_cool=5+1.5×Q^0.8 W / (m²·K), where Q is the irrigation flow rate. Table 4 shows the time constants and temperatures after 10 s for different flow rates, and the cooling curves and fitting relationships are shown in Table 4. Figure 8 .

[0046] Table 4. Calibration Results of Flushing Flow Rate and Cooling Parameters

[0047] 0 5.0 30.0 67.8 10 14.5 10.4 53.4 30 27.8 5.4 43.7 50 39.3 3.8 40.1

[0048] The system selects a set of voxels V_s within the safe neighborhood of the heat-sensitive tissue and records the temperature-time curve T_i(x,t) of each voxel. When CEM43 is used, the incremental heat dose generated by the i-th path segment on voxel x is: DeltaDose_i(x)=∫{t_i}^{t{i+1}} R^{43-T_i(x,t)} dt. Here, T_i(x,t) is in °C, dt is in minutes, and R is configured according to the temperature range, tissue type, or experimental calibration results. When using the Arrhenius model, the equivalent thermal damage integral is: Omega_i(x)=∫{t_i}^{t{i+1}} A_0·exp[-E_a / (R_g·T_K(x,t))] dt, where A_0 is the frequency factor, E_a is the activation energy, R_g is the gas constant, and T_K(x,t) is the Kelvin temperature. The system can use CEM43, the Arrhenius integral, or other equivalent thermal damage indicators as cumulative heat dose indicators. When using the Arrhenius model, the system updates the cumulative thermal damage integral according to Omega_acc,i(x)=Omega_acc,i-1(x)+Omega_i(x).

[0049] The system does not perform static scoring on a single path, but rather performs cross-segment state recursion on the path segments sequentially. For path segment s_i, the system first calculates the initial residual temperature field of the current path segment based on the end temperature field at the end of the previous path segment and the update conditions of the cooling term or residual temperature field, and then calculates the incremental heat dose DeltaDose_i(x) of the current path segment based on the initial residual temperature field.

[0050] In one embodiment, the recursive relationship of the residual temperature field is as follows:

[0051] T_start,i(x)=F_cool(T_end,i-1(x), Delta t_i, K(x), h_cool).

[0052] Where T_end,i-1(x) is the end temperature field at the end of the previous path segment, Delta t_i is the cooling pause time, K(x) is the thermal conductivity tensor, h_cool is the cooling heat transfer parameter, and T_start,i(x) is the initial residual temperature field of the current path segment.

[0053] For the first path segment, the initial temperature field T_start,1(x) can be set to the baseline temperature field, body temperature field, ambient temperature field measured before or during the operation, or the bone surface temperature field obtained by initialization based on infrared thermography / temperature probe.

[0054] The cumulative heat dose is updated according to Dose_acc_i(x) = Dose_acc_{i-1}(x) + DeltaDose_i(x); when using the Arrhenius model, the cumulative thermal damage integral is updated according to Omega_acc,i(x) = Omega_acc,i-1(x) + Omega_i(x). Cooling pause affects the initial residual temperature field of subsequent path segments, but does not attenuate the already formed CEM43 heat dose or Arrhenius thermal damage integral.

[0055] Using T_start,i(x) as the initial condition, solve for the temperature-time curve T_i(x,t) of the current path segment, and calculate DeltaDose_i(x) from T_i(x,t).

[0056] The system further calculates Dose_max, Dose_mean, V_dose, dose gradient, and cool-down recovery time. Here, Dose_max is the maximum voxel dose within the safe neighborhood, Dose_mean is the average voxel dose within the safe neighborhood, V_dose = |{x∈V_s: Dose_acc(x)>Dose_threshold(type(x))}| / |V_s|, where V_s is the set of voxels in the safe neighborhood of heat-sensitive tissues, and Dose_threshold(type(x)) is the dose threshold corresponding to the tissue type of voxel x. Different dose thresholds or different safe neighborhood thicknesses can be used for nerves, spinal cord, inner ear, blood vessels, dura mater, mucous membranes, and skin.

[0057] The system allows configuration of the CEM43 threshold and safety neighborhood thickness according to tissue type. Table 5 provides example ranges for conservative, standard, and lenient values. Figure 9 The system displays a visual configuration of the threshold and safety neighborhood. The listed values ​​are used for planning parameter initialization and do not replace clinical safety thresholds; the system allows physicians to select more conservative thresholds based on tissue type, source of evidence, and risk preference.

[0058] Table 5. Examples of tissue-specific CEM43 thresholds and safe neighborhood configurations

[0059] Facial nerve 5 20 50 0.5 1.0 2.0 Diederich 2005 vestibular cochlear nerve 1 5 15 0.3 0.8 1.5 Diederich 2005 spinal cord 2 10 30 1.0 2.0 3.0 Dewey 1984 nerve root 3 15 40 0.8 1.5 2.5 extrapolation dura mater 10 50 200 0.5 1.0 2.0 Conservative estimate Blood vessel 20 100 500 0.5 1.0 2.0 Conservative estimate oral mucosa 50 200 1000 0.3 0.5 1.0 Conservative estimate skin 100 500 2000 1.0 2.0 3.0 Conservative estimate Inferior alveolar nerve 5 20 50 0.5 1.0 2.0 Pogrel 1995

[0060] In one embodiment, the optimization objective function includes J = w1·A_residue + w2·A_collision + w3·A_dose + w4·A_path + w5·A_time + w6·A_cooling. Here, A_dose is obtained by combining Dose_max, Dose_mean, V_dose, and the dose gradient. The system selects the path segment order and parameter combination with the lowest overall cost, provided that Dose_max does not exceed the maximum dose threshold and V_dose does not exceed the proportional threshold. Joint optimization can also simultaneously satisfy geometric constraints such as coverage of the area to be removed, distance between the device and prohibited structures, mechanical accessibility, collision cost, and processing residue cost.

[0061] The path segment sequence and instrument parameters can be solved using greedy or heuristic sorting, constraint optimization, metaheuristic search, or Bayesian optimization. Greedy or heuristic sorting prioritizes path segments that are far from the safe neighborhood of the same heat-sensitive tissue; constraint optimization uses Dose_max and V_dose as hard constraints, and path length, processing time, and cooling time as soft objectives; metaheuristic search or Bayesian optimization is used to jointly search for rotational speed, feed rate, dwell time, and cooling pause time.

[0062] In one embodiment, the system executes the following process:

[0063] Step S1: Input CT, CBCT, MRI, or fused images, and segment the area to be removed and the safe neighborhood of heat-sensitive tissue.

[0064] Step S2: Initialize Dose_acc_0(x)=0 for the safe neighborhood voxel set V_s, establish the patient's individualized bone tissue thermophysical parameter field, and set the initial temperature field T_start,1(x) for the first path segment according to the baseline temperature field, body temperature field, or intraoperative initial temperature measurement data.

[0065] Step S3: Generate a candidate path set P={P_1,P_2,...,P_m}, and discretize each candidate path into path segments {s_1,s_2,...,s_n}.

[0066] Step S4: For each path segment s_i, establish a moving heat source term Q_i(x,t), and calculate the initial residual temperature field T_start,i(x) of the current path segment based on the end temperature field at the end of the previous path segment and the update conditions of the cooling term or residual temperature field.

[0067] Step S5: Solve the local temperature curve T_i(x,t) of the current path segment using T_start,i(x) as the initial condition, calculate DeltaDose_i(x) or Omega_i(x), and monotonically update the cumulative heat dose or cumulative thermal damage integral according to Dose_acc_i(x)=Dose_acc_{i-1}(x)+DeltaDose_i(x) or Omega_acc,i(x)=Omega_acc,i-1(x)+Omega_i(x).

[0068] Step S6: Calculate Dose_max and V_dose; if Dose_max or V_dose exceeds the threshold, adjust the path segment sequence, rotation speed, feed rate, dwell time, or cooling pause time.

[0069] Step S7: Iterate until the constraints are met or the solution with the lowest overall cost is selected, and output the recommended path, instrument parameters, cooling strategy and thermal dose risk map.

[0070] Step S8: Receive temperature, power, torque, force, pose, or flushing flow data during the procedure, update model parameters, and recalculate the remaining path segments.

[0071] To meet the computational needs of intraoperative planning, the system can establish a basic thermal field model offline using finite element, finite difference, or finite volume methods, and pre-calculate thermal response kernels for commonly used path segments or local regions. During intraoperative updates, the system can perform local updates only within the safe neighborhood of heat-sensitive tissues and adjacent bone regions.

[0072] In another embodiment, the system approximates the transient thermal field response using a reduced-order model, lookup table, or machine learning surrogate model, and uses GPU parallel computing to accelerate voxel-level temperature profile and thermal dose calculation. During intraoperative feedback correction, the system preferably updates only the heat source intensity coefficient, cooling heat transfer coefficient, and thermal conductivity correction coefficient, rather than reconstructing the entire field model each time.

[0073] In a computing environment with an Intel Xeon W-2295 processor, NVIDIA RTX A6000 GPU, and 128 GB of memory, taking a safe neighborhood of 8000 voxels as an example, eight solution or approximate update methods were tested. The computation time for a single path segment, the time for recalculating 10 segments, and the speedup ratio are shown in Table 6. Figure 10 Pre-computation kernels, reduced-order models, lookup tables, and GPU parallel solutions can all keep the recalculation of the remaining path within approximately 2 seconds.

[0074] Table 6 Computational performance of different thermal field solution methods

[0075] Finite element method (full area) 45.000 540.0 1× × Finite element method (local) 8.500 102.0 5× × Finite difference (local) 3.200 38.4 14× △ Finite volume (local) 4.800 57.6 9× △ Pre-computation kernel 0.150 1.8 300× √ Reduced order model 0.080 1.0 562× √ Lookup table 0.020 0.2 2250× √ GPU parallel 0.050 0.6 900× √

[0076] During the procedure, the system receives the temperature distribution of the exposed bone surface from an infrared thermal imager and maps it to the preoperative image coordinate system using surface point clouds, instrument pose, optical markers, or image registration. For obscured or deep areas, the system uses a transient thermal diffusion model to infer the temperature change trend of deep risk areas from the surface temperature. The deep temperature result is a model estimate, and the system can output a confidence interval, a conservative upper bound, or a risk level.

[0077] The system defines the error term e(t) as T_surface_measured(t) - T_surface_predicted(t). When e(t) exceeds Delta_T, the predicted Dose_max exceeds Dose_threshold, V_dose exceeds the proportional threshold, or a sudden change occurs in the device's power / torque, the system triggers model correction. The correction targets include the heat source intensity coefficient lambda, the cooling heat transfer coefficient h_cool, the thermal conductivity correction coefficient gamma_K, and boundary condition parameters.

[0078] In a simulation / playback verification example consisting of 12 temporal mastoid grinding scenarios, infrared thermography and K-type thermocouples were used to form a temperature reference, and an extended Kalman filter was used to correct the model. After correction, the MAE decreased from 4.93 ℃ to 1.42 ℃, the RMSE decreased from 6.43 ℃ to 1.89 ℃, the maximum error decreased from 19.68 ℃ to 5.32 ℃, and the R² increased from -1.0245 to 0.8253. See Table 7 for details. Figure 11 These results are used to illustrate the evaluation method of the feedback correction process and do not constitute a conclusion on clinical effectiveness.

[0079] Table 7 Temperature prediction performance before and after extended Kalman filter correction

[0080] MAE (°C) 4.93 1.42 71.2%↓ RMSE (°C) 6.43 1.89 70.6%↓ Maximum error (°C) 19.68 5.32 73.0%↓ R² -1.0245 0.8253 Significant improvement Deviation (°C) 3.25 0.84 Corrected Standard deviation (°C) 5.55 1.69 69.5%↓

[0081] Model calibration can employ least squares, recursive least squares, extended Kalman filtering, or Bayesian update. After calibration, the system re-predicts the temperature-time curves and cumulative heat dose indices for the remaining path segments, and outputs auxiliary planning results such as reducing rotational speed, shortening dwell time, increasing cooling pauses, changing processing sequence, or switching candidate paths.

[0082] In one embodiment, the experimental verification used isolated pig bone, bovine bone, or artificial bone models, set up simulated areas of heat-sensitive tissue, and used an infrared thermal imager and embedded thermocouples as temperature references. The experiment recorded instrument rotation speed, feed rate, contact pressure, flushing flow rate, cooling pause time, torque, power, and orientation data.

[0083] The comparison groups included: a path planning group considering only geometric collision avoidance, a planning group using a homogeneous bone thermophysical property model, a planning group without intraoperative feedback correction, and a planning group using the recursive thermal dose constraint of this invention. Evaluation metrics included Dose_max, Dose_mean, V_dose, maximum temperature, temperature prediction error, planning time, and intraoperative recalculation time.

[0084] In a simulation comparison example using an isolated porcine femur / temporal bone composite model, 10 samples were set in each group to compare geometric collision avoidance only, homogeneous bone thermal model, no intraoperative feedback correction, and the scheme of the present invention. The results are shown in Table 8 and Figure 12 The proposed solution has a Dose_max of 7.9±2.4 min, a V_dose of 3.7±2.3%, a maximum temperature of 41.7±1.8 ℃, a temperature prediction error of 1.4±0.4 ℃, and a remaining path recalculation time of 1.0±0.3 s.

[0085] Table 8 Simulation comparison results of the four planning methods

[0086] A: Geometric collision avoidance only 45.0±13.8 34.2±7.3 54.8±8.2 N / A 1.8±0.4 N / A B: Homogeneous thermal model 28.1±6.4 24.9±6.9 51.7±3.8 5.5±1.1 43.6±7.4 N / A C: No feedback during the procedure 15.4±4.9 12.9±2.6 46.0±3.6 5.2±0.8 51.3±11.3 N / A D: The present invention solution 7.9±2.4 3.7±2.3 41.7±1.8 1.4±0.4 49.1±10.9 1.0±0.3

[0087] The experimental data described above can be used to determine or correct calibration curves for Dose_threshold, maximum dose threshold, proportional threshold, heat source intensity coefficient, cooling heat transfer coefficient, and thermal conductivity, and to evaluate the performance of this invention in terms of Dose_max, V_dose, temperature prediction error, and replanning time.

[0088] The parameter sources may include: CT-HU to bone mineral density mapping curves from calibration phantoms or standard material scans; k_parallel, k_perp1, and k_perp2 from ex vivo bone thermal conductivity experiments, literature, or lookup tables; η_i from temperature rise calibration experiments under different instrument and bone mineral density combinations; h_cool from irrigation flow rate and exposure area calibration; and Dose_threshold and proportional thresholds from tissue type, literature thresholds, animal experiments, or ex vivo model validation.

[0089] For temporal mastoid process grinding, spinal laminectomy, maxillofacial osteotomy, and bone tumor curettage, the image source, instrument parameters, cooling strategy, heat-sensitive tissues, safety neighborhood, and dose threshold can be configured according to Table 9. The planning time for each scenario is 35–78 s, and the remaining path recalculation time is 0.8–1.8 s; parameter comparisons are shown below. Figure 13 .

[0090] Table 9. Parameter configuration examples for four typical bone tissue processing scenarios.

[0091] Image source CT 0.3mm, CBCT 0.2mm CT 0.5mm + MRI T2 1.0mm CBCT 0.2mm CT 0.5mm + MRI T1+C 1.0mm Number of path segments 12 18 8 22 instrument Diamond grinding head φ2.0mm Diamond grinding head φ3.0mm Oscillating saw / grinding head φ2.3mm scraper + grinding head φ2.0mm rotational speed 40,000-60,000 30,000-50,000 N / A / 50,000-80,000 20,000-40,000 feed rate 0.5-1.5 mm / s 0.8-2.0 mm / s 2.0-5.0 mm / s 0.3-1.0 mm / s cool down 30 mL / min of physiological saline + suction 40 mL / min of physiological saline + suction 20 mL / min of physiological saline 30 mL / min of physiological saline + suction security organizations Facial nerve, vestibular cochlear nerve Spinal cord, nerve roots, dura mater Inferior alveolar nerve Blood vessels, nerve bundles Safe Neighborhood 0.5-1.0 mm 1.0-2.0 mm 0.5-1.0 mm 1.0-2.0 mm Dose threshold CEM43=20min CEM43=10min CEM43=20min CEM43=50min Planning Time 48 s 62 s 35 s 78 s Recalculation time 1.2 s 1.5 s 0.8 s 1.8 s

[0092] Example 1

[0093] The system receives CT images of the patient's temporal bone, segments the mastoid bone, plans the area to be removed, the facial nerve canal, and inner ear structures, and extends the facial nerve canal outwards to form a thermally sensitive tissue safety neighborhood (d_safe). The system generates a local bone thermophysical parameter field based on CT values ​​and the degree of mastoid bone pneumatization.

[0094] The system generates multiple candidate grinding paths and discretizes each path into multiple grinding segments. For path segments close to the facial nerve safety neighborhood, the system predicts Dose_max, Dose_mean, and V_dose. If a path meets the geometric distance requirements but the predicted Dose_max exceeds Dose_threshold, the system lowers the score of that path or configures the path segment with a lower rotation speed, shorter dwell time, and longer cooling pause.

[0095] Example 2

[0096] The system receives fused CT and MRI images of the spine. CT is used for bone tissue and bone density modeling, while MRI is used to assist in segmenting the spinal cord, dural sac, and nerve roots. The system establishes a safe neighborhood of heat-sensitive tissue around the nerve roots and spinal cord, and divides the area to be removed from the lamina into multiple processing segments.

[0097] The system simulates different laminar grinding sequences. If continuous machining of adjacent segments leads to thermal superposition within the same nerve root neighborhood, the system separates and sorts the adjacent segments, machining the segments furthest from the nerve root first, and then returning to process the adjacent segments after the thermal field of the previous segment has decayed. The system then outputs the path segment sequence, rotational speed, feed rate, and cooling pause suggestions.

[0098] Example 3

[0099] The system receives maxillofacial CBCT images, segments the mandible, osteotomy area, and inferior alveolar nerve canal, and estimates bone thermal properties based on CBCT grayscale normalization results and local cortical bone thickness. The system establishes moving heat source models with different spatial distributions for oscillating saws, milling cutters, or grinding drills.

[0100] When the geometric distance between the candidate osteotomy line and the inferior alveolar nerve canal is similar, the system ranks the candidate osteotomy lines according to the predicted cumulative thermal dose risk. For osteotomy lines with high thermal dose risk, the system outputs auxiliary planning results such as reducing local residence time, increasing irrigation and cooling, or changing the osteotomy sequence.

[0101] Example 4

[0102] The system receives tumor-related CT, MRI, or enhanced images, and segments the lesion bone region, the planned curettage area, adjacent neurovascular bundles, and areas of cortical bone weakness. The system generates a heterogeneous thermophysical parameter field based on the bone density and degree of bone destruction in different regions.

[0103] The system divides the scraping path into multiple segments from shallow to deep. For the deep segment near the vascular bundle, the system predicts the residual thermal field generated by the previous scraping and adjusts the subsequent scraping sequence and instrument parameters accordingly to ensure that V_dose is below a preset proportional threshold in the safe neighborhood of the vascular bundle.

[0104] The thresholds, weights, calibration curves, lookup tables, and model parameters involved in the above embodiments can all be configured according to the target tissue type, instrument type, imaging equipment, cooling conditions, and experimental calibration results. All equivalent substitutions, parameter adjustments, or module combinations made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A bone tissue processing auxiliary planning system based on cumulative heat dose constraint, characterized in that, The system includes a processor and a memory, the memory storing a computer program. When executed by the processor, the computer program causes the system to perform the following steps: acquiring preoperative medical images and segmenting the region to be removed and the safe neighborhood of heat-sensitive tissue; mapping bone density based on CT values, CBCT grayscale correction values, calibration phantom conversion results, or bone density estimates obtained from multimodal fusion images in the preoperative medical images, and generating a patient-specific bone tissue thermophysical parameter field by combining the bone structure orientation tensor; discretizing the candidate bone tissue processing path into multiple path segments, and for each path segment, based on instrument position, processing time, rotation speed, etc. The feed rate and residence time establish the moving heat source term; the cooling term or residual temperature field update conditions are determined based on the cooling pause time, flushing flow rate, convective heat transfer coefficient of the exposed surface, or a combination thereof; the heat diffusion model without heat source or with cooling term is solved based on the end temperature field at the end of the previous path segment and the cooling term or residual temperature field update conditions to obtain the initial residual temperature field of the current path segment; using the initial residual temperature field as the initial condition, the transient heat diffusion model of the current path segment is solved based on the bone tissue thermophysical parameter field and the moving heat source term to obtain the temperature-time curves of multiple voxels in the safe neighborhood of the heat-sensitive tissue; The incremental cumulative heat dose generated at the multiple voxels in the current path segment is calculated based on the temperature-time curve. According to the candidate processing order, the incremental cumulative thermal dose of each path segment is monotonically accumulated to obtain the cumulative thermal dose index corresponding to the candidate processing order. The cumulative thermal dose index includes the maximum cumulative thermal dose Dose_max and the voxel proportion V_dose exceeding the dose threshold. Under the constraints that Dose_max does not exceed the maximum dose threshold and V_dose does not exceed the proportion threshold, the path segment order, rotation speed, feed rate, dwell time and cooling pause time are jointly optimized to generate bone tissue processing auxiliary planning results. The bone tissue processing auxiliary planning results include recommended path, instrument parameters, cooling strategy and thermal dose risk map.

2. The bone tissue processing auxiliary planning system based on cumulative thermal dose constraint according to claim 1, characterized in that, The safe neighborhood of the heat-sensitive tissue is obtained by distance field calculation, morphological expansion or risk expansion based on tissue-specific dose threshold of the segmented region of the heat-sensitive tissue. Its thickness or dose threshold is determined by physician configuration, anatomical structure category configuration, literature threshold or experimental calibration results. The heat-sensitive tissue includes at least one of nerve, blood vessel, spinal cord, nerve root, dura mater, inner ear structure, inferior alveolar nerve canal, mucosa or skin.

3. The bone tissue processing auxiliary planning system based on cumulative thermal dose constraint according to claim 1, characterized in that, The bone tissue thermophysical parameter field includes at least an anisotropic thermal conductivity tensor K(x), which is expressed as: K(x) = R(A(x))·diag(k_parallel(rho(x)). k_perp1(rho(x)),k_perp2(rho(x)))·R(A(x))^T; where rho(x) is the bone density obtained by mapping or estimation from CT values, CBCT grayscale correction values, calibration phantom conversion results, or multimodal fusion images; A(x) is the bone structure orientation tensor; R(A(x)) is the rotation matrix determined by A(x); k_parallel, k_perp1, and k_perp2 are the thermal conductivity along the principal direction and the two transverse directions, respectively; the bone structure orientation tensor A(x) is obtained from the grayscale gradient structure tensor, the principal direction of bone trabecular texture, the mean intercept length analysis, local high-resolution images, or anatomical orientation templates; when the confidence level of the bone structure orientation is lower than the preset confidence threshold, the system adopts an isotropic or orthotropic simplified thermophysical property model.

4. The bone tissue processing auxiliary planning system based on cumulative thermal dose constraint according to claim 1, characterized in that, The moving heat source term is a spatially distributed heat source, including at least one of a Gaussian distributed moving heat source, a spherical contact heat source, a cylindrical contact heat source, or a surface heat source along the contact surface of the instrument; the heat source intensity of the moving heat source term is jointly estimated by at least two of the following: rotational speed, torque, axial force, feed rate, instrument power, instrument diameter, tool wear state, and cooling state; in one embodiment, the instrument power is calculated according to P_tool(t)=2π·rpm(t)·Torque(t) / 60, and the moving heat source term is calculated according to Q_i(x,t)=η_i·P_tool(t)·G(∥x-p_i(t)∥,σ_i), where G is the integrally normalized spatial heat source core, and η_i is the heat proportion coefficient entering the bone tissue, calibrated through in vitro bone experiments or artificial bone experiments.

5. The bone tissue processing auxiliary planning system based on cumulative thermal dose constraint according to claim 1, characterized in that, The transient heat diffusion model satisfies: ρ_b(x)c_b(x)·∂T(x,t) / ∂t=∇·(K(x)∇T(x,t))+Q_i(x,t)-H_cool(x,t); Wherein, ρ_b(x) is the equivalent density of bone tissue, c_b(x) is the specific heat capacity, K(x) is the thermal conductivity tensor, Qi(x,t) is the moving heat source term corresponding to the i-th path segment, and H_cool(x,t) is the cooling term; the cooling term is determined based on one or more of the following: flushing flow rate, flushing fluid temperature, cooling pause time, convective heat transfer coefficient of exposed surface, and negative pressure attraction state; the system uses at least one of the following methods to solve or update the transient thermal diffusion model: finite element method, finite difference method, finite volume method, path segment thermal response kernel pre-calculation, reduced-order model, lookup table, or GPU parallel computing.

6. The bone tissue processing auxiliary planning system based on cumulative thermal dose constraint according to claim 1, characterized in that, The incremental cumulative heat dose is calculated using CEM43 heat dose: DeltaDose_i(x)=∫{t_i}^{t{i+1}} R^{43-T_i(x,t)} dt, where T_i(x,t) is the predicted temperature of the i-th path segment at position x, in °C, dt is in minutes, and R is the temperature conversion coefficient, which is configured according to the temperature range, tissue type, or experimental calibration results; or it can be calculated using Arrhenius thermal damage integral: Omega_i(x)=∫{t_i}^{t{i+1}}A_0·exp[-E_a / (R_g·T_K(x,t))] dt, where A_0 is the frequency factor, E_a is the activation energy, R_g is the gas constant, and T_K(x,t) is the Kelvin temperature; the recursion includes residual temperature field recursion and cumulative heat dose update: based on the end temperature field T_end,i-1(x) at the end of the previous path segment and the cooling pause time Delta t_i, solve the heat diffusion model without heat source or with cooling term to obtain the initial residual temperature field T_start,i(x) of the current path segment, use T_start,i(x) as the initial condition to solve the temperature-time curve of the current path segment and calculate the incremental heat dose, and update the cumulative heat dose according to Dose_acc_i(x)=Dose_acc_{i-1}(x)+DeltaDose_i(x); Cooling pauses are used only to reduce or change the initial residual temperature field of subsequent path segments, and are not used to attenuate the already formed CEM43 thermal dose or Arrhenius thermal damage integral.

7. The bone tissue processing auxiliary planning system based on cumulative thermal dose constraint according to claim 1, characterized in that, When Dose_max exceeds the maximum dose threshold or V_dose exceeds the proportional threshold in the candidate processing sequence, the system performs at least one of the following processes: adjusting the path segment order, reducing the rotation speed, reducing the heat source power, shortening the dwell time, changing the feed rate, inserting a cooling pause, increasing flushing cooling, or removing the candidate path segment; the joint optimization is also performed under the conditions that the coverage of the area to be removed is not lower than the coverage threshold, the minimum distance between the instrument and the prohibited structure is not less than the safety distance, the instrument collision cost is lower than the collision threshold, the mechanical accessibility meets the preset accessibility conditions, and the processing residual cost is lower than the residual threshold.

8. The bone tissue processing auxiliary planning system based on cumulative thermal dose constraint according to claim 1, characterized in that, The system receives intraoperative measurement data and registers it to the preoperative medical imaging coordinate system. The intraoperative measurement data includes one or more of the following: infrared thermography, temperature probe, instrument power, torque, axial force, instrument pose, or irrigation flow rate. The system performs model correction based on the error term e(t) = T_surface_measured(t) - T_surface_predicted(t). When e(t) exceeds the temperature error threshold, Dose_max exceeds the maximum dose threshold, V_dose exceeds the proportional threshold, or there is a sudden change in instrument power / torque, the system uses least squares, recursive least squares, extended Kalman filtering, or Bayesian update to update the heat source intensity coefficient, cooling heat transfer coefficient, thermal conductivity correction coefficient, or boundary conditions, and recalculates the cumulative heat dose index of the remaining path segment.

9. A bone tissue processing-assisted planning method based on cumulative heat dose constraints, characterized in that, The method, executed by a computer device, includes: acquiring preoperative medical images and segmenting the area to be removed and the safe neighborhood of heat-sensitive tissue; mapping bone density based on CT values, CBCT grayscale correction values, calibration phantom conversion results, or bone density estimates obtained from multimodal fusion images, and generating a patient-specific bone tissue thermophysical parameter field by combining the bone structure orientation tensor; discretizing the candidate bone tissue processing path into multiple path segments, establishing a moving heat source term for each path segment based on instrument position, processing time, rotation speed, feed rate, and dwell time, and determining the cooling term or residual temperature field update conditions based on cooling pause time, flushing flow rate, exposed surface convective heat transfer coefficient, or a combination thereof; solving a heat diffusion model without a heat source or with a cooling term based on the end temperature field at the end of the previous path segment and the cooling term or residual temperature field update conditions to obtain the initial residual temperature field of the current path segment; and using the initial residual temperature field as the initial... The transient heat diffusion model of the current path segment is solved under conditions to obtain the temperature-time curves of multiple voxels within the safe neighborhood. Based on the temperature-time curves, the incremental cumulative heat dose formed by the current path segment is calculated. The incremental cumulative heat dose of each path segment is monotonically accumulated according to the candidate processing order to obtain the maximum cumulative heat dose Dose_max and the proportion of voxels exceeding the dose threshold V_dose. Under the constraints that Dose_max does not exceed the maximum dose threshold and V_dose does not exceed the proportion threshold, the path segment order, rotation speed, feed rate, dwell time, and cooling pause time are jointly optimized to generate bone tissue processing auxiliary planning results. Intraoperative measurement data is received and registered to the preoperative medical image coordinate system. The heat source intensity coefficient, cooling heat transfer coefficient, or thermal conductivity correction coefficient are updated based on the error between the predicted temperature and the measured temperature, and the cumulative heat dose index of the remaining path segments is recalculated.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of claim 9.

Images