An automated cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics

Abstract

The present invention belongs to the field of medical image processing and application, and relates to a medical image subsequent processing and analysis technology, and specifically relates to an automated cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, especially a three-dimensional hemodynamics automated calculation process, whose calculation program includes a model preprocessing module, a tracer particle module, a track statistics module, a skeleton extraction module, a gridding module, and a fluid solution module; its main program is based on the python language, and the complex function module is encapsulated as a dynamic library based on C++ for the main program to call. The process only requires the user to provide a blood vessel stl model, and specify the inlet and outlet surfaces of the blood flow in the model and the corresponding boundary conditions, and the blood flow field quantity results can be automatically executed and output. The present invention can complete the blood flow simulation calculation with as little interaction as possible with the clinician, and judge the aneurysm rupture risk as extremely high risk, high risk, medium risk, and low risk groups, which is convenient for using physical simulation technology to assist the surgical plan and is helpful for clinical application.

Description

Technical Field

[0001] The present invention belongs to the field of medical image processing and application, and relates to the subsequent processing and analysis technology of medical images, and specifically relates to an automated cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, especially a three-dimensional hemodynamics automated calculation process, whose calculation program includes a model preprocessing module, a tracer particle module, a track statistics module, a skeleton extraction module, a gridding module, and a fluid solution module, which can automatically execute the program that outputs the blood flow field quantity results. It can complete the blood flow simulation calculation, judge the risk of aneurysm rupture, facilitate the use of physical simulation technology to assist the surgical plan, and contribute to clinical application. Background Art

[0002] Research reports that the exact cause of intracranial aneurysms is still unclear, but studies have shown that multiple factors, including hemodynamics, genetics, and vascular biomechanics, play an important role in the occurrence of intracranial aneurysms. Current research generally believes that hemodynamic factors play an important role in the occurrence, development, and rupture of aneurysms, making the adjustment of aneurysm hemodynamics one of the goals of intracranial aneurysm treatment.

[0003] The existing hemodynamic experimental research can be divided into two categories: in vivo animal experimental research and in vitro model experimental research. Research practice shows that in vivo experiments can obtain the most realistic data, but in vivo experiments have limitations such as long cycles, experimental results are greatly affected by individual differences, and long-term control and monitoring of hemodynamic parameters are extremely difficult; in vitro model experiments have made up for some limitations to some extent, such as the use of color Doppler detection technology and magnetic resonance vascular imaging technology, but in vitro experiments still have limitations such as it is difficult to approximate model experiments with actual flow conditions, and the flow parameters extracted from model experiments are extremely conventional and limited; therefore, hemodynamic simulation is very necessary, which can make up for the shortcomings of the experiment itself to a certain extent.

[0004] CFD is a systematic analysis of related physical phenomena including fluid flow and heat conduction through computer numerical calculation and image display. Its basic idea is to replace the fields of physical quantities originally connected in the time domain and space domain, such as velocity field and pressure field, with a set of variable values ​​at a finite number of discrete points, and establish a set of algebraic equations about the relationship between the variables at these discrete points through certain principles and methods, and then solve the algebraic equations to obtain the approximate values ​​of the field variables. CFD can be regarded as a numerical simulation of flow under the control of the basic flow equation. Through this numerical simulation, the distribution of basic physical quantities (such as velocity, pressure, temperature, concentration, etc.) at various positions in the flow field of extremely complex problems and the changes of these physical quantities over time are obtained. At present, the simulation of hemodynamics using CFD software is mainly divided into the following processes: first, use medical image processing software such as Mimics to convert the DICOM format images collected from the clinic, reconstruct the patient's vascular three-dimensional geometry (STL format) and save it; then import the three-dimensional vascular geometry into computational fluid meshing software such as ANSYS-ICEM to generate a computational grid for numerical simulation; then set boundary conditions for the model inlet and outlet, such as blood pressure, blood flow velocity, flow rate, etc.; set the property parameters of blood and blood vessel walls, such as density, viscosity, elasticity, etc.; finally, start the calculation and extract the visualization parameters after convergence; the reliability of the numerical simulation study depends on whether the above steps are carried out completely correctly. In order to obtain reliable numerical simulation research results, researchers must strictly control the quality of each step and strive to obtain the same parameters as the actual situation of the patient at each step.

[0005] In clinical practice, for simple intracranial aneurysms that have ruptured before, corresponding surgical treatment can be directly selected. However, for some complex aneurysms that have ruptured before (including multiple aneurysms and malformed aneurysms), simulation analysis can provide great help for clinical diagnosis and treatment, helping doctors to have a deeper understanding of the hemodynamic characteristics of aneurysms, helping to determine the maximum stress position of the aneurysm wall and the possible rupture position. Doctors can refer to the calculation results and combine them with the anatomical structure of nearby blood vessels to select the appropriate surgical method (interventional embolization, craniotomy and aneurysm clipping); for example: if numerical simulation predicts that the VonMises stress at the base of a saccular aneurysm is significantly increased, the choice of interventional embolization may increase the risk of spring plugs piercing the aneurysm wall and causing intracranial hemorrhage. In this case, craniotomy and clipping are the best choice. Of course, previous studies have not yet reached a conclusion on which surgical method is the best based on the results of numerical simulation studies. Further numerical simulation studies of a large number of intracranial aneurysms and studies on the relationship between the results of numerical simulation studies and surgical method selection and prognosis are needed.

[0006] Numerical calculation of intracranial aneurysm is an interdisciplinary research. Blood flow calculation is based on CFD program, and mechanical analysis of blood vessel wall is based on CAE program. The whole analysis process involves fluid mechanics, mathematics, computer software, computer simulation, finite element method, materials science, topology, medicine, etc. The process of analysis using commercial software is also more biased towards mechanics. At present, there are few people who fully master the above research methods, and the methods have not been widely promoted in clinical work. The most important part of the simulation process is the regional discretization process, that is, the fluid area or solid area needs to be divided into computational grids before the calculation and solution can be performed. The quality of the grid directly determines the accuracy of the result. If starting from the stl model, the grid production process and the calculation process are packaged into a "black box" with minimal interaction with the user, then the whole simulation process becomes three stages: image extraction, model production, and boundary definition. The settings of these three stages are all within the professional scope of clinicians, and the threshold of simulation will be greatly reduced.

[0007] Based on the current status of the prior art, the inventors of this application intend to provide an automated cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, especially a three-dimensional hemodynamic automated calculation process. The analysis system can complete blood flow simulation calculations, determine the risk of aneurysm rupture, facilitate the use of physical simulation technology to assist surgical plans, and facilitate clinical applications.

[0008] References related to the present invention include:

[0009] [1]Updegrove,A.,Wilson,N.,Merkow,J.,Lan,H.,Marsden,ALand Shadden,SC,SimVascular-An open source pipeline for cardiovascular simulation,Annalsof Biomedical Engineering(2016).DOI:10.1007 / s10439-016-1762-8

[0010] [2] Liu Fengrui, Huang Jianping, Li Zhong, Wang Naxiu, et al. Two-stage gravity ball algorithm for generating single-size particle accumulation in cylindrical container[J]. Nuclear Technology, 2017. Summary of the invention

[0011] The purpose of the present invention is to provide an automated cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics based on the current state of the prior art, especially a three-dimensional hemodynamic automated calculation process. The risk analysis system of the present invention can automatically execute the program for outputting blood flow field quantity results, complete blood flow simulation calculations, determine the risk of aneurysm rupture, facilitate the use of physical simulation technology to assist surgical plans, and contribute to clinical applications.

[0012] Specifically, the present invention provides an automated cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, which includes an automated process for calculating three-dimensional hemodynamics, characterized in that the blood flow automated calculation program in the automated process includes a model preprocessing module, a tracer particle module, a track statistics module, a skeleton extraction module, a gridding module, and a fluid solution module;

[0013] The automated process only requires the user to provide a blood vessel stl model and specify the inlet and outlet surfaces of the blood flow in the model and the corresponding boundary conditions, and the program can automatically execute and output the blood flow field quantity results.

[0014] In the present invention, the tracer particle module extracts the topological information of the blood vessel model based on the particle dynamics method, thereby avoiding the problem of algorithm instability of the traditional topology-based skeleton extraction algorithm.

[0015] In the present invention, the track statistics module injects tracer particles based on the Monte Carlo concept, the initial positions of the particles are randomly distributed on the end surface, the initial velocities of the particles are also randomly distributed, and the injection process is repeated continuously to enhance the particles' ability to capture vascular geometry. The injection process ends when the track category is stable.

[0016] In the present invention, the skeleton extraction module performs maximum inscribed sphere "correction" on the initial skeleton obtained by track statistics, and the corrected skeleton provides the location and volume information of the hemangioma.

[0017] In the present invention, the mesh division module uses the segmentation concept to divide the model into sub-models according to the skeleton, and meshes each sub-model separately, and there is no topological constraint between the sub-grids.

[0018] In the present invention, the fluid solver in the fluid solving module performs coupled solving for different blood flow areas.

[0019] In the present invention, the dynamic parameters of blood flow are automatically calculated with little interaction with the clinician. The clinician is required to provide a vascular model in STL format based on the DICOM image. The program re-triangulates the STL and makes the distribution of triangular faces as uniform as possible under the premise that the shape of the model file remains unchanged, which is convenient for subsequent mesh division and particle-vascular wall collision detection; for sharp areas in the model, the program appropriately smoothes the model to approximate the real characteristics of the blood vessel; based on the VTK library, a model interaction interface is provided to the clinician for specifying the inlet and outlet of the vascular cluster and the related blood flow boundary conditions. This is also the only link in the present invention where the program needs to interact with the user. Based on the current cutting and naming operation, clinicians are required to qualitatively cut at the inlet and outlet to mark the inlet and outlet positions. For the inlet boundary, clinicians can specify on the program input card that the inlet surface is extended by a certain distance along the external normal direction. In order to ensure the full development of fluid flow, at each inlet and outlet end face of the vascular cluster, the program starts the tracer particle injection operation, and injects different batches of particles at different speeds according to the Monte Carlo idea. After the particle trajectory is stable, a statistical analysis is performed on the tracer particle trajectory to mark the vascular branches and the areas with vortex trajectories inside. Then, the maximum inscribed sphere correction is performed on the statistically obtained vascular skeleton based on the tracer particle expansion algorithm. The corrected bone The distribution of the skeleton and the maximum inscribed sphere radius can be used clinically; after the skeleton information is obtained, the model is segmented based on the vortex area, and the segmented model set is meshed with a boundary layer separately. For the area where the segmentation is unsuccessful, it degenerates into a grid without a boundary layer, but as a compromise for the degradation of the grid quality, the grid size without the boundary layer is set to a small size; the meshing requires the user to provide the inlet blood flow conditions in advance, so as to determine the grid size parameters such as the thickness of the grid boundary layer based on the blood flow Reynolds number; finally, the boundary conditions specified by the clinician and the generated partitioned grid are combined for fluid calculation, and its specific execution requires calling a dedicated fluid solver for multi-zone grids;

[0020] The entire program execution process in the present invention is as follows Figure 3 As shown; the core algorithm is divided into four modules: statistical tracer particle trajectory / vascular skeleton extraction / model meshing and fluid solver.

[0021] In the present invention, for a given stl model, the program automatically identifies all the head and end faces of the vascular cluster, removes the end faces to be processed later, and makes the model open at the end faces (subsequently injects tracer particles), and then the program traverses the head and end faces of each vascular cluster, and uses the Monte Carlo idea to randomly throw some spherical particles from the end faces, allowing these particles to move forward along the end face entrance under the action of random initial velocity, and sets the external force field of the entire system to drive the particle movement (the effect is similar to the pressure difference in the fluid flow equation). At this time, the system only has particle-particle forces and particle-triangular unit forces (i.e., particle-vascular wall effects), and the particle trajectory is tracked. The topological characteristics of the particle group can be summarized by statistical methods, providing basic information for subsequent skeleton extraction and meshing.

[0022] The present invention provides an automated cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, in particular, a three-dimensional hemodynamics automated calculation process, wherein the calculation program includes a model preprocessing module, a tracer particle module, a track statistics module, a skeleton extraction module, a gridding module, and a fluid solution module, which are designed to automatically execute a program that outputs blood flow field quantity results. The main program is based on the Python language, and the complex function module is encapsulated as a dynamic library based on C++ for the main program to call. The process only requires the user to provide a vascular STL model, and specify the inlet and outlet surfaces of the blood flow in the model and the corresponding boundary conditions, and then the blood flow field quantity results can be automatically executed. The present invention can complete the blood flow simulation calculation with as little interaction as possible with the clinician, and judge the aneurysm rupture risk as extremely high risk, high risk, medium risk, and low risk groups, which is convenient for using physical simulation technology to assist the surgical plan, and is helpful for clinical application. The present invention greatly reduces the technical threshold of blood flow simulation, and can facilitate clinicians to simulate and observe blood flow details before surgery, and provide reference information for the surgical plan. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] Figure 1 This is the execution flow of the discrete element method.

[0024] Figure 2 It is a flow chart of fluid calculation.

[0025] Figure 3 It is the overall execution flow of the main program. DETAILED DESCRIPTION

[0026] Example 1

[0027] The particle discrete element method (DEM) is used to calculate the motion of particles under external forces. Figure 1As shown in the figure, the forces between particles and between particles and patches are equivalent to a spring-damper model, where the spring and damping are applied in the normal direction n and the tangential direction t respectively; for two contact particles i and j, with radii Ri and Rj respectively, the control equation for the angular velocity of particle i is:

[0028]

[0029] Where mi, Ii represent the mass and moment of inertia of particle i, respectively, and Nc is the number of neighboring particles in contact with particle i. The program uses a simplified Hertz–Mindlin–Deresiewicz contact model to calculate the contact force between particles:

[0030]

[0031] Where k and γ represent the stiffness and damping coefficient; represents the relative velocity between particles; μ represents the friction coefficient. The tangential deformation and normal deformation between particles are:

[0032]

[0033] The stiffness coefficient k and damping coefficient γ in the tangential and normal directions are:

[0034]

[0035]

[0036] γ n =γ t

[0037] Where R*, m* and E* are:

[0038]

[0039] E, G and υ are Young's modulus, shear modulus and Poisson's ratio respectively, ε is the coefficient of restitution,

[0040] The maximum time step size is determined by the Rayleigh time step:

[0041]

[0042] The position and velocity of the particles evolve over time and are iteratively updated using the Leap-Frog algorithm, which has third-order accuracy and is time-reversible. The position and velocity are updated using the position of the current time step, the velocity of the first half of the time step, and the force on the particle at the current moment.

[0043]

[0044] in They represent the velocity, position and force of the particle at time t,

[0045] The DEM algorithm uses a neighbor list to search for particle neighbors. After completing the neighbor search, it performs strict collision detection and applies particle-particle and particle-vascular wall interaction forces based on the relative positions between the particles and their neighbors. Finally, Newton's laws are used to update the displacement, velocity, angular velocity and other dynamic information of each particle. The process is as follows: Figure 1 shown.

[0046] The first characteristic of the tracer particle trajectory is the connectivity characteristics of different inlets and outlets, which is similar to the definition of the divergence of field quantities in differential equations. Multiple trajectories may intersect in a certain area or separate at a certain point. The particle mass flow rate on a certain section can be used as a sign of line division and merging in topology. The place where the mass flow rate increases represents the merging point of the blood flow branch, and vice versa. The tracer particle enters the vascular network from a certain end face, and moves forward under the action between particles and vascular walls until it flows out of the vascular network from a certain end face. During this period, a complete particle trajectory can represent a vascular pathway, and a set of vascular pathways can be obtained by traversing all end faces. Merging vascular pathways can form a preliminary vascular model skeleton;

[0047] The model skeleton should be able to reflect the overall topological structure of the model in a concise way, while also maintaining the branching structure of the model. In order to improve the quality of model skeleton extraction, it is necessary to ensure that the model skeleton points can be well concentrated on the model's central axis. The set of pathways obtained by statistically tracing particle trajectories depends on the tracing method in the Monte Carlo process, and is not necessarily close to the center of the blood vessel. Therefore, it is necessary to "correct" it based on the three-dimensional maximum inscribed sphere principle. The specific method is to place virtual particles of a certain radius on these pathways at intervals, so that the particle radius gradually increases, the particle-particle interaction is cancelled, and the particle-wall interaction force is simply considered. The particle radius is increased until it moves under the action of the blood vessel wall and reaches stability. The corrected vascular pathway can be obtained. At this time, the radius information of the maximum inscribed sphere at the position of the virtual particle on the vascular pathway has been obtained, and the topological information of the vascular system can be used for clinical application. The algorithm is similar to the particle expansion algorithm in the "two-stage gravity ball algorithm for generating single-size particle accumulation in a cylindrical container". The particle dynamics algorithm is used to avoid some limitations of the pure topological algorithm [2], which increases the robustness of the algorithm and is more suitable for embedding into the automation process.

[0048] Another characteristic of the tracer particle trajectory is the vortex characteristic in the local area, which can be regarded as the connectivity characteristic when the inlet and outlet inside the blood vessel overlap (such as the blood flow vortex inside the hemangioma). Similar to the definition of the curl of the field quantity, the tracer particles sometimes flow out from a certain position on the blood vessel wall and then flow in again. Its trajectory is a closed curve. At this time, the blood flow divergence is 0 but the curl is not 0. The mass flow rate of the tracer particles inside the vortex is not 0, and the mass flow rate of the tracer particles in the upstream and downstream is unchanged. The existence of the vortex cannot be determined. However, at this time, the corrected vascular skeleton will degenerate into a relatively short branch line segment in the vortex area. Therefore, for these non-user-marked end faces, the tracer particle mass flow rate can be detected in the area represented by the branch line segment to determine the vortex area.

[0049] Example 2

[0050] Since the vortex area is often the area of ​​clinical concern and is usually the key position that affects the quality of the grid, it is difficult to properly divide all areas by traditionally dividing the entire vascular model into grids. It is often necessary to make compromises and try many times to obtain a relatively complete overall grid, which requires a lot of human-computer interaction. In the present invention, the program divides the stl model of the blood vessel according to the vortex area, and each subset after the division is independently meshed. The grid consistency is not required between different vascular units. Therefore, the topological shape and blood flow characteristics of each grid unit are relatively simple, and a fluid grid that conforms to its flow characteristics can be generated for each vascular subset, which is the key to the automation of the entire simulation process.

[0051] After the grid subset is generated for the blood vessel subset, the fluid solver needs to support the information transmission of different sub-vessel interfaces (the traditional CFD typical algorithm PISO / SIMPLE is only for a single flow area). The fluid solver in the program of the present invention adopts the PISO algorithm to iteratively calculate the flow in each sub-grid, and interpolate the velocity, pressure and other field quantities at the interface of different sub-vessels; based on computational fluid dynamics (CFD) to calculate fluid dynamics information, the control equations followed mainly include the continuity equation:

[0052]

[0053] For a flow with constant density, the above equation simplifies to:

[0054]

[0055] Momentum conservation equation:

[0056]

[0057] The first term on the left represents the rate of change of the micro-group velocity; the second term on the left represents the convection term; the first term on the right represents the pressure gradient term; the second term on the right represents the external volume force acting on the micro-group; the third term on the right represents the micro-group stress; the last term represents gravity; for isotropic fluids, the fluid relationship is:

[0058]

[0059] The program adopts the turbulence model based on the Reynolds-averaged NS equations:

[0060]

[0061] Before the calculation begins, the clinician needs to specify the vascular boundary conditions. This is the only place where the algorithm needs to interact with the user. The user needs to use the mobile interface to move / scale / rotate the section, and use the section obtained by cutting this rectangular section with the stl model as the blood flow inlet and outlet surface. Then, in the save window, enter the surface type (inlet / outlet), inlet boundary (blood flow / velocity distribution model / stretching length), outlet boundary (outlet pressure) and other flow information; the fluid solver performs calculations based on the user-specified boundaries and physical parameters. For each subnet area, the momentum equation is first solved with the given initial pressure or the pressure of the previous iteration step. Because the obtained velocity variable does not necessarily satisfy the continuity equation, the pressure Poisson equation needs to be solved to obtain the pressure, and the obtained pressure is used to correct the velocity so that it can satisfy the continuity equation. If the velocity does not satisfy the momentum equation, repeat the cycle until convergence. The program flow is as follows: Figure 2 shown.

[0062] The present invention can complete blood flow simulation calculations and determine the risk of aneurysm rupture into extremely high risk, high risk, medium risk, and low risk groups, making it convenient to use physical simulation technology to assist surgical plans and facilitating clinical applications.

Claims

1. An automated cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, comprising an automated process for calculating three-dimensional hemodynamics, characterized in that: The blood flow automation calculation program in the automation process includes a model preprocessing module, a tracer particle module, a track statistics module, a skeleton extraction module, a meshing module and a fluid solution module; The automated process only requires the user to provide a blood vessel stl model and specify the inlet and outlet surfaces of the blood flow in the model and the corresponding boundary conditions, and the program can automatically execute and output the blood flow field quantity results; The tracer particle module extracts topological information of the vascular model based on a particle dynamics method; In the track statistics module, tracer particles are injected based on the Monte Carlo concept. The initial positions of the particles are randomly distributed on the end surface, and the initial velocities of the particles are also randomly distributed. The injection process is repeated continuously to enhance the particle's ability to capture the vascular geometry. The injection process is terminated when the track category is stable. The DEM particle discrete element method is used to calculate the movement of the particles. In the skeleton extraction module, the initial skeleton obtained by track statistics is corrected by maximum inscribed sphere based on the tracer particle expansion algorithm, and the corrected skeleton provides the location and volume information of the hemangioma; The meshing module uses the idea of ​​segmentation to segment the model into sub-models based on the skeleton, and meshes each sub-model separately, without topological constraints between sub-grids. The fluid solver in the fluid solving module performs coupled solving for different blood flow areas; the fluid solver adopts the PISO algorithm to perform iterative calculation on the flow in each sub-grid respectively.

Images