CN113887146B

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

Info

Publication number: CN113887146B
Application number: CN202011094195.2A
Authority: CN
Inventors: 朱凤平; 毛颖
Original assignee: Huashan Hospital of Fudan University
Current assignee: Huashan Hospital of Fudan University
Priority date: 2020-07-02
Filing date: 2020-10-14
Publication date: 2025-02-11
Anticipated expiration: 2040-10-14
Also published as: CN113887146A

Abstract

The invention belongs to the field of medical image processing and application, relates to a medical image subsequent processing analysis technology, in particular to an automatic cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, and particularly relates to a three-dimensional hemodynamics automatic calculation flow, wherein the calculation program comprises a model preprocessing module, a trace particle module, a track statistics module, a skeleton extraction module, a gridding sub-module and a fluid solving module; the main program is based on the python language, and the complex functional module is packaged into a dynamic library based on C++ for calling the main program. The flow can automatically execute the flow quantity result of the transfusion flow only by providing a vascular stl model by a user and designating an inlet and outlet surface and corresponding boundary conditions of the blood flow in the model. The invention can complete blood flow simulation calculation by interaction with clinicians as few as possible, judges that the rupture risk of the aneurysm is respectively an extremely high-risk group, a medium-risk group and a low-risk group, is convenient to use a physical simulation technology to assist an operation scheme, and is beneficial to clinical application.

Description

Automatic cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics

Technical Field

The invention belongs to the field of medical image processing and application, relates to a medical image subsequent processing analysis technology, in particular to an automatic cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, and particularly relates to a three-dimensional hemodynamic automatic calculation flow. The blood flow simulation calculation can be completed, the aneurysm rupture risk can be judged, and the physical simulation technology is convenient to assist the operation scheme, thereby being beneficial to clinical application.

Background

The exact cause of the occurrence of the intracranial aneurysm is not completely clear, and the research shows that various factors such as blood flow dynamics (hemodynamics), genetics, vascular biomechanics and the like play an important role in the occurrence of the intracranial aneurysm. Current research generally suggests that hemodynamic factors play an important role in the development, progression, and rupture of aneurysms, making modulation of the hemodynamics of aneurysms one of the targets for intracranial aneurysm treatment.

The prior art hemodynamic experimental study can be divided into two categories, namely animal in-vivo experimental study and in-vitro model experimental study. The research practice shows that the most realistic data can be obtained in vivo experiments, but the in vivo experiments have the limitations of longer period, great influence of individual difference on experimental results, extremely difficult long-time control and monitoring of hemodynamic parameters and the like, in vitro model experiments make up for certain limitations such as a color Doppler detection technology, a magnetic resonance vascular imaging technology and the like to a certain extent, in vitro experiments have the limitations that the model experiments are difficult to approximate to actual flow conditions, and the flow parameters extracted from the model experiments are extremely conventional and limited and the like, so that hemodynamic simulation is necessary, and the defects of the experiments per se can be overcome to a certain extent.

The CFD is a system analysis of related physical phenomena including fluid flow and heat conduction through computer numerical calculation and image display, and the basic idea is to replace the fields of physical quantities originally connected in the time domain and the space domain, such as a speed field and a pressure field, with a set of variable values at a series of limited discrete points, establish an algebraic equation set related to the relationship between the variables at the discrete points through a certain principle and mode, and then solve the algebraic equation set to obtain the approximation of the field variable. CFD can be regarded as a numerical simulation of flow under the control of the flow basis equation, by which the distribution of basic physical quantities (such as velocity, pressure, temperature, concentration, etc.) at various locations in the flow field, and the change over time of these physical quantities, is obtained, which is an extremely complex problem. The current CFD software simulation hemodynamics mainly comprises the following procedures of firstly converting a DICOM format image acquired clinically by using medical image processing software such as Mimics, reconstructing a three-dimensional geometry (STL format) of a patient blood vessel, then guiding the three-dimensional blood vessel geometry into computational fluid meshing software such as ANSYS-ICEM to generate a computational grid for numerical simulation, setting boundary conditions such as blood pressure, blood flow rate, flow and the like for an inlet and an outlet of a model, setting attribute parameters such as density, viscosity, elasticity and the like of blood and a blood vessel wall, finally starting calculation until convergence, extracting visual parameters, and obtaining reliable numerical simulation research results according to whether the steps are completely correct or not, wherein a researcher has to strictly control the quality of each step to obtain the parameters which are the same as the actual conditions of the patient.

In clinical practice, for a previously ruptured intracranial aneurysm with simple morphology, corresponding surgical treatment can be directly selected, but for some previously ruptured complex aneurysms (including multiple aneurysms and malformed aneurysms), simulation analysis can provide great help for clinical diagnosis, help doctors to deeply understand hemodynamic characteristics of the aneurysms, help to judge the maximum stress position and possible rupture position of the wall of the aneurysms, and refer to calculation results of the aneurysms, and by combining the anatomical structures of nearby blood vessels, the doctor can select proper surgical modes (interventional embolism and craniotomy aneurysm occlusion), for example, if numerical simulation predicts that the stress of the saccular aneurysm bottom VonMises is obviously increased, the risk of intracranial hemorrhage caused by puncturing the wall of the aneurysm by a spring bolt is increased, the craniotomy occlusion is preferably selected, and of course, the numerical simulation research results are selected as best as the prior research theory, a great quantity of numerical simulation research on the intracranial aneurysms is required, and the research on the relation between the results of the numerical simulation research and the selection and prognosis of the surgical modes is not formed.

The calculation of the numerical value of the intracranial aneurysm belongs to the research of the interdisciplinary science, the calculation of the blood flow is based on a CFD program, the mechanical analysis of the vascular wall is based on a CAE program, the whole analysis process involves fluid mechanics, mathematics, computer software, computer simulation, a finite element method, materiality, topology, medicine and the like, and the process of analysis by using commercial software is also more superior to the mechanics profession, so that less personnel who fully master the research method are at present, and the method is not promoted in a large amount in clinical work. The most important ring in the simulation process is a region discrete process, namely, a fluid region or a solid region is divided into calculation grids, then calculation solution can be carried out, and the accuracy of a result is directly determined by the grid quality. If the stl model is started, the grid making process and the calculating process are packed into a 'black box' with minimum interaction with a user, the whole simulation process is changed into three stages of image extraction, model making and boundary definition, the settings of the three stages all belong to the professional category of clinicians, and the simulation threshold is greatly reduced.

Based on the current state of the art, the inventor aims to provide an automatic cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, in particular to a three-dimensional hemodynamics automatic calculation flow, and the analysis system can complete blood flow simulation calculation, judge aneurysm rupture risk, is convenient to use a physical simulation technology to assist an operation scheme, and is beneficial to clinical application.

References relevant to the present invention are:

[1]Updegrove,A.,Wilson,N.,Merkow,J.,Lan,H.,Marsden,A.L.and Shadden,S.C.,SimVascular-An open source pipeline for cardiovascular simulation,Annals of Biomedical Engineering(2016).DOI:10.1007/s10439-016-1762-8

[2] Liu Fengrui, huang Jianping, li Zhong, wang Naxiu, etc. a two-stage gravity ball drop algorithm [ J ] core technique, 2017, for single particle size particle packing is generated in a cylindrical vessel.

Disclosure of Invention

The invention aims to provide an automatic cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics based on the current state of the art, in particular to a three-dimensional hemodynamic automatic calculation flow, and the risk analysis system can automatically execute a program of an output blood flow result to complete blood flow simulation calculation, judge the aneurysm rupture risk, is convenient to use a physical simulation technology to assist an operation scheme and is beneficial to clinical application.

Specifically, the invention provides an automatic cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, which comprises an automatic flow for calculating three-dimensional hemodynamics, and is characterized in that a blood flow automatic calculation program in the automatic flow comprises a model preprocessing module, a trace particle module, a track statistics module, a skeleton extraction module, a gridding sub-module and a fluid solving module;

the automatic process only needs a user to provide a vascular stl model, and designates an inlet and outlet surface and corresponding boundary conditions of blood flow in the model, and the process can automatically execute the blood flow output result.

In the invention, the topological information of the vessel model is extracted from the trace particle module based on the particle dynamics method, so that the problem of unstable algorithm of the traditional skeleton extraction algorithm based on topology is avoided.

In the invention, trace statistics module is used for injecting trace particles based on Monte Carlo idea, the initial positions of the particles are randomly distributed on the end face, the initial speeds of the particles are also randomly distributed, the capturing capacity of the particles to the vessel geometry is enhanced by continuously repeating the injection process, and the injection process is ended after the trace category is stable.

In the invention, the initial skeleton obtained by track statistics is subjected to maximum inscribed sphere correction in the skeleton extraction module, and the corrected skeleton provides hemangioma position and volume information.

In the invention, the grid division module uses the division idea to divide the model into the sub-models according to the skeleton, and each sub-model is singly subjected to grid division, and no topological constraint exists between the sub-grids.

In the invention, the fluid solver in the fluid solver module performs coupling solving for different blood flow areas.

According to the invention, the dynamic parameters of blood flow are automatically calculated under the condition of small amount of interaction with a clinician, the clinician is required to provide a blood vessel model in stl format according to a dicom image, the program re-triangulates stl, triangular patches are distributed uniformly as much as possible on the premise of unchanged shape of a model file, subsequent grid division and particle-blood vessel wall collision detection are facilitated, the program is suitable for smoothly approaching the real characteristics of blood vessels in a sharp area in the model, a model interaction interface of the clinician is provided based on a VTK library, and the model interaction interface is used for specifying the inlet and outlet of a blood vessel cluster and related blood flow boundary conditions, which is the only link of program interaction with a user in the invention. Based on the cutting naming operation adopted at present, a clinician is required to qualitatively mark the inlet and outlet positions at an inlet and outlet, for an inlet boundary, the clinician can appoint the inlet surface on a program input card to extend a distance along the external normal direction, in order to ensure the full development of fluid flow, at each inlet and outlet end face of a vascular cluster, the program starts the tracer particle injection operation, particle groups with different batches and different speeds are injected according to the Monte Carlo thought, statistical analysis is carried out on the tracer particle track after the particle track is stable, the vascular branches and the areas with vortex tracks inside are marked, then the maximum inscription sphere correction is carried out on the vascular skeleton obtained by statistics further based on the tracer particle expansion algorithm, the corrected skeleton and the maximum inscription sphere radius distribution can be used for clinic, after skeleton information is obtained, the model is divided based on the vortex area, the divided model set is singly divided into grids with boundary layers, and for the unsuccessfully divided areas, the grid is degenerated into grids without boundary layers, but the grid size setting of the boundary layers is small as a compromise processing mode of grid quality reduction, thereby, the user is required to provide inlet conditions in advance, the calculation is carried out on the special fluid flow condition of the special-purpose fluid which is required to be calculated by a plurality of areas, and the size of the blood flow parameters is required to be calculated and the specific to be calculated;

The whole program execution flow in the invention is shown in figure 3, and the core algorithm is divided into four modules, namely four major parts of a statistic trace particle track/vascular skeleton extraction/model grid division and a fluid solver.

In the invention, for a given stl model, firstly, a program automatically identifies all head and tail end surfaces of a blood vessel cluster, the end surfaces to be processed later are removed so that the model is in an open state at the end surfaces (trace particles are injected later), then the program traverses the head and tail end surfaces of each blood vessel cluster, and by utilizing the thought of Monte Carlo, some spherical particles are randomly thrown out from the end surfaces, so that the particles move forwards along an inlet of the end surface under the action of random initial speed, and an external force field of the whole system is set to drive the particles to move (the action is similar to the pressure difference in a fluid flow equation). At this time, the system only has the acting force between particles and triangular units (namely the acting force between particles and the blood vessel wall), tracks the particle track, can induce the topological characteristic of the particle group by using a statistical method, and provides basic information for the subsequent skeleton extraction and grid division.

The invention provides an automatic cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics, in particular to a three-dimensional hemodynamic automatic calculation flow, which comprises a model preprocessing module, a trace particle module, a track statistics module, a skeleton extraction module, a gridding division module and a fluid solving module. The main program is based on the python language, and the complex functional module is packaged into a dynamic library based on C++ for calling the main program. The flow can automatically execute the flow quantity result of the transfusion flow only by providing a vascular stl model by a user and designating an inlet and outlet surface and corresponding boundary conditions of the blood flow in the model. The invention can complete blood flow simulation calculation by interaction with clinicians as few as possible, judges that the rupture risk of the aneurysm is respectively an extremely high-risk group, a medium-risk group and a low-risk group, is convenient to use a physical simulation technology to assist an operation scheme, and is beneficial to clinical application. The invention greatly reduces the technical threshold of blood flow simulation, can facilitate clinicians to simulate and observe blood flow details before operation, and provides reference information for operation schemes.

Drawings

Fig. 1 is a flow chart of execution of the discrete unit method.

FIG. 2 is a flow chart of fluid computation.

Fig. 3 is a flow of execution of the entire main routine.

Detailed Description

Example 1

Calculating the movement of particles under the action of external force, adopting a particle discrete element method (DEM, a flow is shown in figure 1), and equivalent the stress between the particles and the particle-surface sheet to be a spring-damping model, wherein the spring and the damping are respectively applied in a normal direction n and a tangential direction t, and for two contact particles i and j, the radius is Ri and Rj respectively, and the control equation of the speed and the angular speed of the particles i is as follows:

Where mi, ii represent the mass and moment of inertia, respectively, of particle i, nc is the number of neighboring particles in contact with particle i, and the program calculates the contact force between particles using a simplified Hertz-Mindlin-Deresiewicz contact model:

wherein k and γ represent hardness and damping coefficients; represents the relative speed between particles and μ represents the coefficient of friction. Tangential deformation and normal deformation among particles are respectively as follows:

The stiffness coefficient k and the damping coefficient γ in the tangential and the normal directions are respectively:

γ_n＝γ_t

Wherein R, m and E are each:

e, G and v are Young's modulus, shear modulus and Poisson's ratio, ε is the coefficient of restitution,

The maximum time step is determined by the rayleigh time step:

The position, speed and the like of the particles evolve with time, and are iteratively updated by adopting a leapfrog algorithm (Leap-Frog algorism), the third-order precision of the algorithm is time-reversible, the position of the current time step, the speed of the first half time step and the stress of the particles at the current moment are used for updating the position and the speed,

Wherein the method comprises the steps ofRespectively representing the speed, the position and the stress of the particles at the moment t,

The DEM algorithm searches particle neighbors by adopting a neighbor list mode, strict collision detection is carried out immediately after the neighbor search is finished, particle-particle and particle-vessel wall interaction forces are applied according to the relative positions of the particles and the neighbors, and finally the dynamic information such as displacement, speed, angular speed and the like of each particle is updated by utilizing Newton's law, and the flow is shown in figure 1.

The first characteristic of trace particle track is that for the communication characteristic of different inlets and outlets, similar to the divergence definition of field quantity in differential equation, multiple tracks may be intersected in a certain section area or separated in a certain point, the particle mass flow on a certain section can be used as the mark of linear combination in topology, the place where the mass flow increases represents the merging point of blood flow branches, and vice versa. The tracer particles enter the vascular network from a certain end face, move forward under the action of the particles and the vascular wall until the tracer particles flow out of the vascular network from a certain end face, a complete particle track can represent a vascular access, and a vascular access set can be obtained by traversing all the end faces. Combining the vascular access to form a preliminary vascular model skeleton;

The model skeleton should be able to embody the overall topology of the model in a compact manner while also well preserving the branching structure of the model. In order to improve the quality of model skeleton extraction, it is necessary to ensure that model skeleton points can be well concentrated at the axial position of the model. The method is characterized in that virtual particles with a certain radius are placed on the passages at intervals to ensure that the particle radius is gradually increased, the particle-particle interaction force is cancelled, the particle radius is increased until the corrected vascular passage can be obtained after the movement under the action of the vascular wall and the stabilization are achieved by simply considering the particle-wall interaction force, at the moment, the radius information of the maximum inscribed sphere at the virtual particle position on the vascular passage is obtained, the topological information of the vascular system can be used for clinical application, the algorithm thought is similar to a particle expansion algorithm in a two-stage gravity ball falling algorithm for generating single-particle-diameter particle accumulation in a cylindrical container, and the limitation [2] on the pure topology algorithm is avoided through a particle dynamics algorithm, so that the robustness of the algorithm is increased and the method is more suitable for being embedded into an automatic flow;

The trace particle track is characterized by vortex characteristics in a local range, can be regarded as communication characteristics (such as blood flow vortex in hemangioma) under the condition that the inlet and the outlet of the inside of a blood vessel are overlapped, is similar to the definition of the rotation of field quantity, trace particles sometimes flow out of a certain position of the wall of the blood vessel and then flow in again, the track is a closed curve, the blood flow divergence is 0 at the moment, the rotation is not 0, the mass flow of the trace particles in the vortex is not 0, the mass flow of the trace particles in the upstream and downstream is unchanged, the existence of the vortex cannot be judged, but the corrected blood vessel skeleton is degraded into relatively short branch line segments in the vortex region, so that the mass flow of the trace particles can be detected in the region represented by the branch line segments on the end faces of the non-user marks, and the vortex region is judged.

Example 2

Based on the vortex region which is usually a clinically concerned region and is usually a key position affecting the grid quality, the whole blood vessel model is divided into grids in the prior art, and the grid is difficult to divide properly in all regions, so that a relatively complete whole grid can be obtained by repeated compromise and try for many times, and a large amount of man-machine interaction is needed;

after generating a grid subset for the vessel subset, the fluid solver correspondingly needs to support information transmission of different sub-vessel interfaces (the traditional CFD typical algorithm PISO/SIMPLE only aims at a single flow area), and in the program, the fluid solver adopts the PISO algorithm to respectively carry out iterative calculation on flow in each sub-grid and simultaneously carries out interpolation of field quantities such as speed, pressure and the like on different sub-vessel interfaces:

The above formula for a flow with a constant density is simplified as:

Momentum conservation equation:

The first left term represents the rate of change of the velocity of the micelle, the second left term represents the convection term, the first right term represents the pressure gradient term, the second right term represents the external volumetric force acting on the micelle, the third right term represents the stress of the micelle, the last term represents the gravity, and for an isotropic fluid the fluid-to-fluid relationship is:

the procedure takes a turbulence model based on the reynolds average N-S equation:

Before calculation starts, a clinician is required to specify vascular boundary conditions, which are the only places where an algorithm needs to interact with a user, the user needs to move an interface to move/zoom/rotate a section, the section obtained by cutting the rectangular section and the stl model is used as a blood flow inlet and outlet surface, then flow information such as a type of an input surface (inlet/outlet) of a storage window, an inlet boundary (blood flow/flow velocity distribution model/stretching length), an outlet boundary (outlet pressure) and the like is stored, a fluid solver calculates based on the boundary and physical parameters specified by the user, and for each sub-network area, the momentum equation is solved by the pressure of the given initial pressure or the last iteration step, the pressure is solved by the poisson equation because the solved speed variable does not necessarily meet the continuity equation, and the pressure is obtained by utilizing the solved pressure correction speed, so that the continuity equation can be met. If the velocity does not satisfy the momentum equation, the loop is repeated until convergence, the process flow is as shown in FIG. 2.

The invention can complete blood flow simulation calculation, judges that the rupture risk of the aneurysm is respectively an extremely high-risk group, a medium-risk group and a low-risk group, is convenient to use a physical simulation technology to assist an operation scheme, and is beneficial to clinical application.

Claims

1. An automated cerebral aneurysm rupture risk analysis system based on cerebral hemodynamics comprises an automated process for calculating three-dimensional hemodynamics, and is characterized in that a blood flow automated calculation program in the automated process comprises a model preprocessing module, a trace particle module, a track statistics module, a skeleton extraction module, a meshing module and a fluid solving module;

The automatic process only needs a user to provide a vascular stl model, and designates an inlet and outlet surface and corresponding boundary conditions of blood flow in the model, and the process can automatically execute the blood flow quantity output result;

The topological information of the vessel model is extracted from the trace particle module based on a particle dynamics method;

The trace statistics module is used for injecting trace particles based on the Monte Carlo idea, the initial positions of the particles are randomly distributed on the end face, the initial speeds of the particles are also randomly distributed, the capturing capacity of the particles to the vascular geometry is enhanced in the injection process continuously and repeatedly, and the injection process is finished after the trace category is stable;

In the skeleton extraction module, the initial skeleton obtained by track statistics is subjected to maximum inscribed sphere correction based on a trace particle expansion algorithm, and the corrected skeleton provides hemangioma position and volume information;

the grid division module uses a division idea to divide the model into sub-models according to the skeleton, and each sub-model is divided into grids independently, and no topological constraint exists between the sub-grids;

and the fluid solver in the fluid solver module carries out coupling solving on different blood flow areas, and carries out iterative calculation on the flow in each sub-grid by adopting a PISO algorithm.