CN102314710A - Medical tissue dynamic simulation method based on force asynchronous diffusion model - Google Patents

Medical tissue dynamic simulation method based on force asynchronous diffusion model Download PDF

Info

Publication number
CN102314710A
CN102314710A CN201110287127A CN201110287127A CN102314710A CN 102314710 A CN102314710 A CN 102314710A CN 201110287127 A CN201110287127 A CN 201110287127A CN 201110287127 A CN201110287127 A CN 201110287127A CN 102314710 A CN102314710 A CN 102314710A
Authority
CN
China
Prior art keywords
voxel
model
particle
physical model
zone
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201110287127A
Other languages
Chinese (zh)
Inventor
袁志勇
司伟鑫
段昭亮
廖祥云
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN201110287127A priority Critical patent/CN102314710A/en
Publication of CN102314710A publication Critical patent/CN102314710A/en
Pending legal-status Critical Current

Links

Images

Abstract

The invention relates to a three-dimensional medical soft tissue model deformation simulation method based on a force asynchronous diffusion model, which includes a preprocessing part for providing a convenient model for a subsequent deformation process, a two-stage collision detection part based on model preprocessing and a high-efficiency model solution part. In order to eliminate the affection of invalid links on the deformation effect, the invention puts forward a regional diffusion algorithm to differentiate internal and external voxels, removes the redundant external voxels of a physical model, which are obtained by division, and achieves the purpose of simplifying the original calculation process under the premise of guaranteeing the surface precision of the physical model. The three-dimensional tissue deformation simulation method presented by the invention can be perfectly applied to the complex organ tissues of the human body.

Description

Medical tissue dynamics simulation method based on the asynchronous diffusion model of power
Technical field
The invention belongs to medical tissue dynamics simulation and virtual reality technology field, particularly relate to a kind of medical tissue dynamics simulation method based on the asynchronous diffusion model of power.
Background technology
Human body complex organ Soft Tissue Deformation is one of gordian technique of medical science dynamics simulation, for example virtual operation emulation, the operative training based on image, operation tracking and location etc.
Operation emulation based on virtual reality is just more and more universal as substituting traditional medical training means, and preceding having obtained widely with the surgery planning field of medical modeling and simulation technology orthopaedic surgical operations operation used.Research contents comprises that the simulations mutual and visual, the various variations of virtual human body organ under the virtual operation instrument effect of medical data reach the simulation of operating personnel's various sensory feedback (like visual feedback and force feedback etc.) etc., is the important application of virtual reality technology in operative training.Pertinent literature: Florian Schulze; Katja B ü hler, Andr é Neubauer, Armin Kanitsar; Leslie Holton; Stefan Wolfsberger. Intra-operative virtual endoscopy for image guided endonasal transsphenoidal pituitary surgery. International Journal of Computer Assisted Radiology and Surgery, 2010,5 (2): 143-154.
In virtual operation emulation, the real-time sense of reality deformation simulation of organ-tissue is an important research contents, is one of gordian technique of virtual operation emulation.1996, Mr.Stava proposed the notion about three generations's medical science analogue system framework in the 4th medical science virtual reality meeting.First generation medical science analogue system focuses on the geometrical property of performance human body, with the roaming in the virtual reality technology with immerse notion and be applied to the human anatomy data set, limited user interactions is provided, and in medical personnel's education and training, is applied.Second generation medical science analogue system is the medical science emulation on the physical significance; The anatomical structure that not only visually presents tissue, but also physical features, the biomechanics characteristic of performance organ-tissue are like the viscoelasticity of organ-tissue; Do the time spent as external force, deformation will take place in organ-tissue.Third generation physical medicine analogue system makes simulated effect more true to nature, has added the physical characteristics of human body as biosome, has caused the Study of calculation model based on physical model and principle.Pertinent literature: Changmok Choi; Jungsik Kim, Hyonyung Han, Bummo Ahn; Jung Kim. Graphic and Haptic Modeling of the Oesophagus for VR-based Medical Simulation; Int J Med Robotics Comput Assist Surg, 5:257-266,2009.
Finite Element Method (FEM) is the most a kind of method of fluid emulation field utilization, also extensively applies to the deformation simulation field.The FEM method can access higher precision, but solving equation larger often, has reduced the real-time of emulation.Halic proposes a kind of Soft Tissue Deformation method based on particle-spring model (MSM); Contrast MSM and FEM method are for the pros and cons of Soft Tissue Deformation; Traditional M SM model can not show the soft tissue physical characteristics fully; When system restriction is not enough, can cause system unstable, and deformation that can not the accurate description soft tissue when distortion is big.In recent years, a plurality of elastic deformation model based on complex geometric shapes and complicated boundary condition are suggested and have obtained good deformation effect.Thomas has proposed a kind of efficiently based on the distorted pattern of data-driven to soft tissue model, reflected the physical characteristics and the physiology characteristic of tissue truly.Faure has proposed a kind of new deformation method that is directed to the complex geometry curved surface and has the isomery material behavior, and geometric properties and material behavior through automatic computing node distribution situation dividing tissue solve the nonuniform stiffness problem effectively.Zhu is implemented in the high resolving power elastic deformation model on 8 cores and the 16 core SMP, can the processing complexity geometric configuration and complicated border condition, and make its concrete needs that adapt to figure and animation application, and the extensibility of multiprocessor shared drive has been discussed.Elastic model under this multi grid framework can efficient emulation different materials parameter with random geometry arbitrarily under distortion, but its collision detection algorithm becomes and finds the solution bottleneck of performance.Gilles will combine to propose a kind of new distorted pattern-Frame-based method with the frame-based method that Muller proposed in 2005 based on the continuum mechanics model of physics; The distance that is about to sampled point current location and target location replaces interaction force, geometrical-restriction relation is replaced the geometry deformation model of energy.Gilles adopts the Lagrangian mechanical equation that reflects model internal physical characteristic to each sampled point on the research basis of Muller; From the sense of reality distortion in complicacy that is deformed to of simple unimodal curve description, the Frame-based elastic deformation model can compare favourably with Finite Element Method.Pertinent literature: Lenka Jerabkova, Torsten Kuhlen, Timm P. Wolter; Norbert Pallua; A Voxel Based Multiresolution Technique for Soft Tissue Deformation. VRST ' 04, November 10-12,2004; Hong Kong, Pages:158-161.
But existing emulation mode all lacks the suitable treatments means to organ-tissue peripheral recesses position, causes precision not enough, and this area is needed new emulation mode badly and occurred.
Summary of the invention
In order to remedy the deficiency of existing deformation simulation method, the present invention proposes a kind of three-dimensional tissue's deformation simulation method based on the asynchronous diffusion model of power.
Technical scheme of the present invention is a kind of medical tissue dynamics simulation method based on the asynchronous diffusion model of power, may further comprise the steps:
Step 1 is set up geometric model, physical model and transition model for three-dimensional soft tissue, and the mode of specifically setting up comprises the geometric model of setting up three-dimensional soft tissue surfaces; Set up physical model through the subdivision geometric model, physical model is made up of voxel, and voxel is made up of the spring between particle and the particle; Set up transition model physical model is mapped to geometric model;
Step 2 by the inside and outside voxel of regional broadcast algorithm differentiation physical model, is removed the outside voxel of physical model, keeps the boundary voxel and the inner voxel of physical model;
The concrete mode of said regional broadcast algorithm does, at first in physical model, gets an inner voxel arbitrarily, is called seed voxels, and the zone that directly links to each other of definition voxel is the zone of the six direction of front and back up and down of this voxel; For seed voxels, extend toward its zone that directly links to each other, if its some voxel that directly link to each other is not a boundary voxel, then must be inner voxel; For newfound inner voxel, extend toward its zone that directly links to each other, if its some voxel that directly link to each other is not a boundary voxel, then must be inner voxel, by that analogy, up to extending to all boundary voxel of physical model;
Step 3 adopts two-stage collision mode to carry out crash tests, confirms collision detection point position according to the excellent position of operation; Said two-stage collision mode comprises finds the solution the regional area that bumps earlier, confirms the position of collision detection point at regional area again;
Step 4 is found the solution the Gaussian curvature of every bit on the geometric model;
Step 5 according to step 4 gained Gaussian curvature, is radiation radius of every bit definition on the geometric model; After soft tissue applied power; The asynchronous diffusion process of physical model stress and deformation is to be that basic point begins toward external diffusion with step 3 gained collision detection point; Diffuse to for the first time the point in the radiation areas of collision detection point, then the point in the radiation areas is all added the zone of action; Zone of action after the diffusion is as first asynchronous area for the first time; For the point of new adding zone of action, in diffusion for the second time, with them as basic point toward external diffusion; If a point is inner in their radiation areas; Then this point is joined the zone of action, the newly-increased point of all basic points has constituted second asynchronous area, by that analogy; All particles up to physical model all add the zone of action, have only the particle in the zone of action to move; The radiation areas of basic point obtain according to the radiation radius of respective point on the geometric model on the said physical model;
Step 6 according to the asynchronous area of step 5 gained physical model, is found the solution the position of each particle in asynchronous diffusion process of physical model, and the deformation effect on the physical model is embodied on geometric model through transition model and draws; The mode that the particle postition is found the solution does; On the physical model spring between the particle is being classified; According to spring classification analysis particle stressing conditions, analyze draw that particle is suffered and make a concerted effort after, utilize particle current time and last a position constantly to find the solution next position constantly of particle.
And, in the step 1,
The geometric model of setting up comprises triangular apex set
Figure 861922DEST_PATH_IMAGE001
and triangle surface set
Figure 2011102871272100002DEST_PATH_IMAGE002
; Wherein
Figure 41230DEST_PATH_IMAGE003
is the set of triangular apex index; is the triangular apex coordinate set;
Figure 503304DEST_PATH_IMAGE005
is the set of triangle surface index, the index value on
Figure 2011102871272100002DEST_PATH_IMAGE006
expression Atria summit;
The physical model of setting up comprises particle set
Figure 84459DEST_PATH_IMAGE007
and set of voxels
Figure 2011102871272100002DEST_PATH_IMAGE008
; Wherein
Figure 606576DEST_PATH_IMAGE009
is the three-dimensional index set of particle;
Figure 335497DEST_PATH_IMAGE004
is the triangular apex coordinate set; is the three-dimensional index set of voxel,
Figure 351995DEST_PATH_IMAGE011
be the three-dimensional index set on eight summits of voxel;
The transition model of setting up comprises voxel-vertex set
Figure 2011102871272100002DEST_PATH_IMAGE012
, summit-set of voxels
Figure 36923DEST_PATH_IMAGE013
and voxel-triangle set ; Wherein
Figure 746253DEST_PATH_IMAGE010
is the three-dimensional index set of voxel;
Figure 696891DEST_PATH_IMAGE015
is the relative position of surface vertices in voxel; is the set of triangular apex index; is the set of triangle surface index; is triangular apex number in the voxel, and
Figure 1992DEST_PATH_IMAGE017
is the triangle number relevant with voxel.
The present invention proposes a kind of based on the asynchronous diffusion model (FADM of power; Force Asynchronous Diffusion Model) three-dimensional Simulation of Soft Tissue Deformation method; It is divided into three parts: for follow-up deformation process provides model preprocessing part easily, based on the pretreated two-stage collision detection of model part and model solution part efficiently., use the inside and outside voxel of regional broadcast algorithm differentiation and remove redundant voxel in the application aspect the distortion of human body complicated tissue organ to FADM, reached the effect of under the situation that guarantees physical model surface accuracy, simplifying original computation process at preprocessing part.Be to embody the physiological characteristic of organ-tissue, the present invention partly introduces surface model Gaussian curvature parameter at model solution, realizes the asynchronous diffusion of organ-tissue, makes grid density degree on the corresponding curved surface of organ-tissue along with model structure characteristic and geometric properties variation.Experimental result shows that the deformation simulation method that the present invention proposes has been removed the invalid link in the concaver, visually can present organ-tissue realistically and receive external force to make the deformation effect of time spent, realizes high-precision medical science dynamics deformation simulation.
Description of drawings
Fig. 1 is the model synoptic diagram of the embodiment of the invention.
Fig. 2 is the two-stage collision detection synoptic diagram of the embodiment of the invention.
Fig. 3 asks for synoptic diagram for the tri patch summit Gaussian curvature of the embodiment of the invention.
Fig. 4 is that synoptic diagram is spread in the zone based on Gaussian curvature in the embodiment of the invention.
Fig. 5 organizes spring and is connected the spring synoptic diagram for the embodiment of the invention.
Fig. 6 is the virtual spring synoptic diagram of the embodiment of the invention.
Fig. 7 is the process flow diagram of the embodiment of the invention.
Embodiment
Simulation of Soft Tissue Deformation relates to the complex mathematical physical process, needs lots of data and sets up the metaplasia model that satisfies physical property and real-time.The present invention proposes a kind of deformation simulation model of novelty-based on the asynchronous diffusion model FADM of power.Referring to Fig. 1, FADM is divided into three parts: model preprocessing part, collision detection part and realistic model are found the solution part.To the soft tissue inner structure, the spring of different effects is classified and is obtained the reference mark through collision detection between the FADM particle.Yet FADM inevitably can be in the face of complicated human organ model, for example organize models such as heart, liver in operation emulation.When FADM is used for complicated human body organ model and carries out deformation simulation; Its geometric model subdivision is being set up in the process of physical model, had invalid connection, original unallied zone passage spring is coupled together at organ-tissue peripheral recesses position; This influences in the deformation simulation of convex body surface model not quite; But to concaver, the existence meeting of invalid connection directly influences the accurate of deformation effect, and the present invention proposes regional broadcast algorithm and distinguishes inside and outside voxel; The outside redundant voxel of physical model that removal obtains because of subdivision, and reached the original computation process of simplification under the situation that guarantees physical model surface accuracy.
The Soft Tissue Deformation physical motion shows as the stressed position change of particle and two processes are successively spread in the moving region.Power begins successively to spread from the reference mark, and the zone of power diffusion process becomes the zone of action, and it is stressed to have only particle in the zone of action just to begin, and the position changes.Every frame traversal all particles calculating location in the classic method, FADM only need calculate the particle of current time zone of action.Division for asynchronous area is mainly judged with the variation severe according to the geometric properties of model surface.Confirm the variation severe of surface model through introducing the surface vertices Gaussian curvature, divide the asynchronous area scope.Change comparatively violent and the complicated place of model area at model surface, it is more careful that asynchronous area is divided, and the zone is less, i.e. asynchronous area division is more careful; Change mild and the uncomplicated place of model area for model surface, the asynchronous area scope division is more greatly that asynchronous area is divided more coarse.
Specify technical scheme of the present invention below in conjunction with accompanying drawing and embodiment.Referring to Fig. 7, the embodiment of the invention may further comprise the steps:
Step 1 is set up geometric model, physical model and transition model for three-dimensional soft tissue, and the mode of specifically setting up comprises the geometric model of setting up three-dimensional soft tissue surfaces; Set up physical model through the subdivision geometric model, physical model is made up of voxel, and voxel is made up of the spring between particle and the particle; Set up transition model physical model is mapped to geometric model.This step belongs to the model preprocessing part.
Embodiment reads model information from stl file, and the geometric jacquard patterning unit surface model of setting up three-dimensional soft tissue is described the topological structure of soft tissue; Complete subdivision geometric model is set up the physical model that has fixed connection between voxel in three-dimensional voxel; The power that at first will act on the soft tissue is converted into the power that acts on the physical model, after the physical model generation deformation, sets up transition model the deformation on the physical model is mapped on the geometric model, to embody the deformation of soft tissue.Stl file is a kind of 3D file layout, gathers with triangle and representes the object external outline shape, is made up of the definition of a plurality of tri patchs, and the definition of each tri patch comprises the three-dimensional coordinate on each summit of triangle and the method vector of tri patch.
Geometric model is the surface model of three-dimensional soft tissue, reflection soft tissue space topological structure.Physical model is that a plurality of square voxels are formed, and voxel is made up of the particle at place, square summit and the spring between particle.Distortion on physical model is found the solution very simple, and can reflect the distortion situation of organization internal, makes distortion have more physical characteristics.For the deformation on the physical model is mapped on the geometric model, the present invention introduces transition model, in order to the deformation of coordinate model and geometric model.Implementation model recovers in step 6, utilizes transition model exactly, and the deformation effect on the physical model is embodied the deformation effect through transition model on geometric model.
The geometric model that the present invention sets up three-dimensional soft tissue to data is described the topological structure of soft tissue; Be divided into two parts: triangular apex set
Figure 2011102871272100002DEST_PATH_IMAGE018
and triangle surface set
Figure 626878DEST_PATH_IMAGE019
, as shown in table 1.
Table 1
Wherein, triangular apex set
Figure 985178DEST_PATH_IMAGE001
reflected triangular apex spatial information and with the leg-of-mutton topological relation in its place.Where is the apex of the triangle index set
Figure 758148DEST_PATH_IMAGE021
; as a triangle vertex coordinates the collection,
Figure 2011102871272100002DEST_PATH_IMAGE022
;
Figure 845369DEST_PATH_IMAGE023
for the triangular facets of the index set.
Figure 2011102871272100002DEST_PATH_IMAGE024
expression triangular apex call number; triangular apex coordinate, expression triangle call number.Triangle surface set has reflected the topological relation on triangle and its three summits, the index value on
Figure 2011102871272100002DEST_PATH_IMAGE028
,
Figure 912399DEST_PATH_IMAGE029
,
Figure 2011102871272100002DEST_PATH_IMAGE030
expression Atria summit.
Physical model is made up of a plurality of three-dimensional voxels, fully the subdivision geometric model.Physical model can be regarded the particIe system that fixed connection is arranged as, and each particle can be used as particle.And different, each voxel contains 8 summits as particle to interparticle acting force with interparticle distance difference, and 12 limits are as spring.The residing cube in voxel as shown in Figure 2 summit
Figure 878081DEST_PATH_IMAGE031
; Cubical 8 summits promptly are particles; Article 12, the limit is as spring, and summit
Figure 59663DEST_PATH_IMAGE031
is one of 8 summits.
As shown in table 2, the present invention is divided into two parts with physical model: particle set
Figure 2011102871272100002DEST_PATH_IMAGE032
and set of voxels
Figure 436287DEST_PATH_IMAGE033
.
Table 2
Figure 2011102871272100002DEST_PATH_IMAGE034
Wherein, The three-dimensional index of
Figure 273793DEST_PATH_IMAGE035
expression particle;
Figure 839903DEST_PATH_IMAGE022
is the triangular apex coordinate set; The three-dimensional index of
Figure 2011102871272100002DEST_PATH_IMAGE036
expression voxel, the three-dimensional index set on eight summits of expression voxel.The three-dimensional call number of
Figure 2011102871272100002DEST_PATH_IMAGE038
expression voxel,
Figure 626780DEST_PATH_IMAGE039
representes the call number on eight summits of voxel respectively.
As shown in table 3; The present invention is divided into three parts with transition model: voxel-vertex set
Figure 2011102871272100002DEST_PATH_IMAGE040
, summit-set of voxels and voxel-triangle set
Figure 2011102871272100002DEST_PATH_IMAGE042
.
Table 3
Figure 508334DEST_PATH_IMAGE043
Wherein, Voxel-vertex set
Figure 398930DEST_PATH_IMAGE012
expression voxel with its related triangle; Relative position
Figure 994700DEST_PATH_IMAGE015
= of summit-set of voxels
Figure 385865DEST_PATH_IMAGE013
presentation surface summit in voxel, triangle number and its index that voxel-triangle set
Figure 105876DEST_PATH_IMAGE014
expression one voxel is comprised.Wherein,
Figure 303508DEST_PATH_IMAGE010
is the three-dimensional index set of voxel;
Figure 271464DEST_PATH_IMAGE003
is the set of triangular apex index;
Figure 305279DEST_PATH_IMAGE005
is the set of triangle surface index;
Figure 16883DEST_PATH_IMAGE016
is triangular apex number in the voxel, and
Figure 69022DEST_PATH_IMAGE017
is the triangle number relevant with voxel.
Figure 207879DEST_PATH_IMAGE003
=
Figure 791307DEST_PATH_IMAGE045
The expression voxel the set of related triangular apex index,
Figure 2011102871272100002DEST_PATH_IMAGE046
The expression summit is in voxel X, y, zThe relative coordinate value.
If the coordinate on eight summits of summit place voxel
Figure 2011102871272100002DEST_PATH_IMAGE048
is
Figure 150930DEST_PATH_IMAGE049
;
Figure 2011102871272100002DEST_PATH_IMAGE050
, then relative position
Figure 398372DEST_PATH_IMAGE044
does
Figure 653117DEST_PATH_IMAGE051
If the coordinate on eight summits of voxel after the deformation
Figure 706524DEST_PATH_IMAGE048
is
Figure 483987DEST_PATH_IMAGE049
;
Figure 964647DEST_PATH_IMAGE050
; Then the reposition of summit
Figure 444038DEST_PATH_IMAGE047
is
Figure 2011102871272100002DEST_PATH_IMAGE052
Figure 301136DEST_PATH_IMAGE053
Step 2 by the inside and outside voxel of regional broadcast algorithm differentiation physical model, is removed the outside voxel of physical model, keeps the boundary voxel and the inner voxel of physical model.This step belongs to the model preprocessing part.
When handling concaver to FADM, the existence meeting of invalid connection directly influences the accuracy of deformation effect, the big limitations of this circumscribed existence its in medical Application in Modeling.The present invention proposes regional broadcast algorithm and distinguishes inside and outside voxel, removes the outside redundant voxel of the physical model that produces because of subdivision, under the situation that guarantees physical model surface accuracy, simplifies original computation process.The process of removing invalid link mainly is to distinguish inside and outside voxel process, keeps the inner voxel of boundary voxel and model, removes outside voxel, can remove redundant invalid link like this, strengthens deformation effect.During practical implementation, the inside and outside voxel of distinguishing the sealing three-dimensional model can also adopt rays method, prior aries such as anglec of rotation method and area symbol diagnostic method.For example basic rays method is a voxel injection line in model, by voxel inside and outside the number judgment of object and ray intersection.But this method exists a large amount of asking to hand over to be calculated, and when handling critical condition, can produce mistake.Therefore the regional broadcast algorithm of the present invention's proposition is an optimum implementation.
Distinguish inside and outside voxel process and can use following method: initial physical model is made as square, chooses the suitable physical model voxel length of side through experiment and covers the geometric model surface fully, makes boundary voxel that inside and outside voxel is separated.If the voxel length of side is too short, boundary voxel just contains the voxel of geometric model surface vertices, not link mutually of boundary voxel then possibly occur, causes outside voxel and inner voxel to communicate, and can influence the differentiation of inside and outside voxel like this.Cutting apart the inside and outside best prerequisite of voxel is that the border separates inside and outside voxel.In case after having distinguished inside and outside voxel, again boundary voxel and inner voxel are kept, remove outside voxel.The specific algorithm of voxel inside and outside distinguishing is at first sought an inner voxel in model center voxel zone, is called seed voxels, and the zone that directly links to each other of definition voxel is the zone of the six direction of front and back up and down of some voxels.Then, extend toward its zone that directly links to each other, if its some voxel that directly link to each other is not a boundary voxel for seed voxels; Then must be inner voxel,, six voxels that directly link to each other judged for newfound inner voxel; If what link to each other is not boundary voxel, then must be inner voxel, adopt same method to judge for new inside voxel; Up to extending to all boundary voxel, can obtain a recursive algorithm of finding the solution inner voxel thus.
Step 3 adopts two-stage collision mode to carry out crash tests, confirms collision detection point position according to the excellent position of operation; Said two-stage collision mode comprises finds the solution the regional area that bumps earlier, confirms the position of collision detection point at regional area again.This step belongs to two-stage collision detection part.
The two-stage collision detection in the model pre-service, through transition model, has been set up the mapping relations of geometric model and physical model based on the model pre-service.In the first step; At first detect the cube voxel (being not limited to the cube voxel of physical model) at operating theater instruments place; Can know through inquiry voxel-vertex set
Figure 933106DEST_PATH_IMAGE040
whether this voxel is the boundary voxel that comprises the geometric model summit; If this voxel is not a boundary voxel; Then do not bump; Then continue to move operating theater instruments, till finding boundary voxel.
As shown in Figure 2, to find after the boundary voxel, second step was sought the triangle set that is associated with this boundary voxel through inquiry voxel-triangle set, and the operating theater instruments and the point of impingement are located on one of them triangle in the triangle set.Regard operating theater instruments as straight line; Adopting M ller-Trumbore method to come calculated line and leg-of-mutton intersection point, promptly is intersection point
Figure 771618DEST_PATH_IMAGE055
.And intersection
Figure 816934DEST_PATH_IMAGE055
Recent voxel vertex
Figure 149826DEST_PATH_IMAGE031
it as a collision detection point.Pertinent literature: Tomas M ller, Ben Trumbore. Fast, minimum storage ray-triangle intersection. Doktorsavhandlingar vid Chalmers Tekniska Hogskola, 1425:109-115,1998.
Traditional collision detection need travel through all summits, sees whether collision detection has taken place, and the present invention adopts the classification collision checking method, can reduce the time complexity of collision detection greatly.
Step 4 is found the solution the Gaussian curvature of every bit on the geometric model.
After soft tissue applied power; The deformation effect of soft tissue zones of different is different; The diffusion of power zones of different on soft tissue also is different, and soft tissue model is divided into different zones, and power is also different because of diffusion radius difference rate of propagation in these different zones; Each zone is called asynchronous area on the soft tissue, and the diffusion process of power on asynchronous area is called the asynchronous diffusion of power.
Division for asynchronous area is mainly judged according to the variation severe of model surface; Asynchronous area is according to the difference of the variation severe of surface model; Its division scope is also different, and power velocity of propagation organizationally is also different, and the scope of zone of action increases also different.Change comparatively violent and the complicated place of model area at model surface, it is more careful that asynchronous area is divided, and the zone is less, i.e. asynchronous area division is more careful; Change mild and the uncomplicated place of model area for model surface, the asynchronous area scope division is more greatly that asynchronous area is divided more coarse.
Gaussian curvature is described the degree of crook of curved surface, and the present invention is through introducing the standard that Gaussian curvature is divided as asynchronous area.The evaluation method of discrete Gaussian curvature has multiple.In these methods, people such as Schneider utilize four PDE to find the solution Gaussian curvature; People such as Mayer utilize the Laplace-Beltrami operator to find the solution Gaussian curvature; People such as Meyer M utilize Voronoi unit and Finite Element Method to find the solution Gaussian curvature.The present invention adopts the Voronoi method of people's propositions such as Meyer M to find the solution Gaussian curvature.Relevant paper: Meyer M; Desbrun M, Schroder P, Barr A H. Discrete differential-geometry operators for triangulated 2-manifolds. In Visualization and Mathematics; Berlin; Germany, 2002,52~58.
During the enforcement reference, the explanation that provides embodiment to find the solution Gaussian curvature is following:
To selected boundary voxel; Get its inner triangle table vertex of surface
Figure 2011102871272100002DEST_PATH_IMAGE056
; If there are
Figure 885570DEST_PATH_IMAGE057
,
Figure 2011102871272100002DEST_PATH_IMAGE058
,
Figure 973612DEST_PATH_IMAGE059
in the summit that is adjacent ... As shown in Figure 3; Its Gaussian curvature
Figure 2011102871272100002DEST_PATH_IMAGE060
does
Wherein
Figure 2011102871272100002DEST_PATH_IMAGE062
is the angle of limit
Figure 567109DEST_PATH_IMAGE063
and ;
Figure 908091DEST_PATH_IMAGE065
is the projected area of point
Figure 2011102871272100002DEST_PATH_IMAGE066
, and.
Figure 353985DEST_PATH_IMAGE067
The definition of
Figure 2011102871272100002DEST_PATH_IMAGE068
is following: if
Figure 45997DEST_PATH_IMAGE069
is oxygon; Then
Figure 783009DEST_PATH_IMAGE068
is the area of polygon
Figure 2011102871272100002DEST_PATH_IMAGE070
, and wherein
Figure 493345DEST_PATH_IMAGE071
and
Figure 2011102871272100002DEST_PATH_IMAGE072
is the point that hangs down.If
Figure 860873DEST_PATH_IMAGE069
is that the angle that obtuse triangle and
Figure 23870DEST_PATH_IMAGE056
locate is an acute angle; Then
Figure 564572DEST_PATH_IMAGE068
is the half the of triangle
Figure 145726DEST_PATH_IMAGE069
area; If the angle that
Figure 480893DEST_PATH_IMAGE056
locates is the obtuse angle, then
Figure 417273DEST_PATH_IMAGE068
is 1/4th of triangle
Figure 496087DEST_PATH_IMAGE069
area.
Step 5 according to step 4 gained Gaussian curvature, is radiation radius of every bit definition on the geometric model; After soft tissue applied power; The asynchronous diffusion process of physical model stress and deformation is to be that basic point begins toward external diffusion with step 3 gained collision detection point; Diffuse to for the first time the point in the radiation areas of collision detection point, then the point in the radiation areas is all added the zone of action; Zone of action after the diffusion is as first asynchronous area for the first time; For the point of new adding zone of action, in diffusion for the second time, with them as basic point toward external diffusion; If a point is inner in their radiation areas; Then this point is joined the zone of action, the newly-increased point of all basic points has constituted second asynchronous area, by that analogy; All particles up to physical model all add the zone of action, have only the particle in the zone of action to move; The radiation areas of basic point obtain according to the radiation radius of respective point on the geometric model on the said physical model.
Among the embodiment; According to step 4 result; The Gaussian curvature of every bit on the known geometric model; And be that it defines a radiation radius
Figure 931748DEST_PATH_IMAGE073
; Then its radiation areas are
Figure 2011102871272100002DEST_PATH_IMAGE074
; By the Gaussian curvature decision of certain point , this diffusion radius is used for the generation of auxiliary asynchronous area.
Through the Gaussian curvature
Figure 778667DEST_PATH_IMAGE077
of other voxel summits in the computational activity zone
Figure 2011102871272100002DEST_PATH_IMAGE076
, calculate the diffusion radius on each summit.The diffusion function of radius is:
Wherein
Figure 2011102871272100002DEST_PATH_IMAGE080
is the initial propagations radius; Set its numerical value according to experiment;
Figure 200607DEST_PATH_IMAGE081
is the radius of influence of voxel summit
Figure 80838DEST_PATH_IMAGE076
;
Figure 2011102871272100002DEST_PATH_IMAGE082
is the mean-gaussian curvature of voxel summit
Figure 705724DEST_PATH_IMAGE076
; is respectively geometric model maximum gaussian curvature and minimal gaussian curvature; is coefficient of diffusion, obtains through experiment.
The radiation areas of setting up an office
Figure 536593DEST_PATH_IMAGE056
are ; The asynchronous diffusion of distortion begins diffusion from the collision detection point; Diffuse to for the first time the point in the radiation areas of the point of impingement, then these points of putting in all radiation areas are all added the zone of action.If collision detection point is
Figure 499575DEST_PATH_IMAGE085
; Then its radiation radius is ; If the zone of action is
Figure 661566DEST_PATH_IMAGE087
; When then spreading for the first time, the point in
Figure 988642DEST_PATH_IMAGE085
radiation areas is joined in the zone of action:
Figure DEST_PATH_IMAGE088
Zone of action after the diffusion is as asynchronous area
Figure 459943DEST_PATH_IMAGE089
for the first time.Point for new adding zone of action; In diffusion for the second time; With them as basic point toward external diffusion; If a point
Figure DEST_PATH_IMAGE090
is inner in their radiation areas, then this point is joined the zone of action:
Figure 810153DEST_PATH_IMAGE091
The part that zone of action after spreading for the second time increases newly is as asynchronous area
Figure DEST_PATH_IMAGE092
.Similarly, we can see first Subdiffusion after the active area for the new part of the asynchronous region
Figure DEST_PATH_IMAGE094
.Initiate like this point further adds the zone of action with the point in its radiation areas again, so just constantly shows diffusion effect, finally can arrive The model.As shown in Figure 4; The asynchronous diffusion of distortion begins diffusion from collision detection point
Figure 878789DEST_PATH_IMAGE095
; Then the particle
Figure DEST_PATH_IMAGE096
in point
Figure 520992DEST_PATH_IMAGE095
radiation areas is joined the zone of action; In the diffusion particle in the radiation areas of particle is being joined the zone of action afterwards next time; With the method, up to all particles are joined the zone of action.
Step 6 according to the asynchronous area of step 5 gained physical model, is found the solution the position of each particle in asynchronous diffusion process of physical model, and the deformation effect on the physical model is embodied on geometric model through transition model and draws; The mode that the particle postition is found the solution does; On the physical model spring between the particle is being classified; According to spring classification analysis particle stressing conditions, analyze draw that particle is suffered and make a concerted effort after, utilize particle current time and last a position constantly to find the solution next position constantly of particle.
On physical model, carry out finding the solution of realistic model, can be divided into three parts: the stressed of particle found the solution, and asynchronous area is divided, the asynchronous diffusion of power.Like Fig. 5 and shown in Figure 6, there are 26 in abutting connection with particle around each particle, 27 particles altogether, 9 every layer.For true embodiment soft tissue inner structure more, the embodiment of the invention is classified to the spring of different effects between the particle:
Spring between the inner particle of organized layer is referred to as to organize spring
Figure 924609DEST_PATH_IMAGE097
; Embody the inner action effect of three-dimensional organ different tissues layer, shown in dot-and-dash line among Fig. 5; The spring of line organization's layer is referred to as to connect spring
Figure DEST_PATH_IMAGE098
between between organized layer; Embody the action effect of three-dimensional organ different tissues interlayer, shown in dotted line among Fig. 5; Through virtual spring
Figure 144719DEST_PATH_IMAGE099
is set; Make the interior transition of three-dimensional tissue's different tissues interlayer and layer more level and smooth, shown in double dot dash line among Fig. 6.Because the viscoelastic property of biological tissue, soft tissue suffered power in deformation process has operation external force
Figure DEST_PATH_IMAGE100
, damping force
Figure 708555DEST_PATH_IMAGE101
and elastic force
Figure DEST_PATH_IMAGE102
.Elastic force has three kinds; They are tissue elasticity power, connect elastic force and virtual elastic force that correspondence is organized spring , connected spring
Figure 386847DEST_PATH_IMAGE098
and virtual spring
Figure 543022DEST_PATH_IMAGE099
respectively respectively.
Stressed being expressed as for particle
Figure 277760DEST_PATH_IMAGE103
Figure DEST_PATH_IMAGE104
If the quality of all particles all is
Figure 73546DEST_PATH_IMAGE105
;
Figure DEST_PATH_IMAGE106
is the acceleration of particle
Figure 184722DEST_PATH_IMAGE103
; If particle
Figure 195403DEST_PATH_IMAGE103
is positioned at
Figure 350310DEST_PATH_IMAGE107
; 26 sets of particles around it are combined into
Figure DEST_PATH_IMAGE108
; The elongation of elastic force and spring is proportional, can be expressed as following formula
Figure 384125DEST_PATH_IMAGE109
Wherein, The position of
Figure DEST_PATH_IMAGE110
expression particle
Figure 20030DEST_PATH_IMAGE103
; The position of
Figure 150797DEST_PATH_IMAGE111
expression spring other end particle; For organizing spring, connecting spring and virtual spring;
Figure DEST_PATH_IMAGE112
is the former length of respective springs; According to the relative position situation of particle, can judge
Figure 227338DEST_PATH_IMAGE113
for organizing spring
Figure 810766DEST_PATH_IMAGE097
, connecting a kind of spring in spring
Figure 247432DEST_PATH_IMAGE098
and the virtual spring
Figure 232706DEST_PATH_IMAGE099
.The movement velocity v of damping force and particle is ratio, and available following formula is represented:
Figure DEST_PATH_IMAGE114
is ratio of damping in the formula;
Figure 480147DEST_PATH_IMAGE115
is time step,
Figure DEST_PATH_IMAGE116
speed of approximate expression particle.Time step is corresponding with each time interval in step of the asynchronous diffusion process of step 5.During practical implementation, can be provided with by those skilled in the art, the time interval is more little, and simulation accuracy is high more.
Through analyzing the particle stressing conditions, after solving that particle is suffered and making a concerted effort, adopt the Verlet integral method to utilize particle current time and last a position constantly to find the solution next position constantly of particle.The Verlet integral method is a kind of existing method that solves equation of motion numerical integration, and the speed of having avoided direct explicit each particle of calculating so that reduce the computation complexity in the deformation process, has obtained widespread use in molecular dynamics simulation emulation.According to the Verlet integral method:
Figure 737822DEST_PATH_IMAGE117
Just can obtain next position constantly of particle thus.The pass of deformation effect on the physical model through transition model tied up to embody the deformation effect on the geometric model, and adopt existing OpenGL shape library software to draw out.
Step 4,5,6 belongs to the model solution part.Owing to set up the kernel function relevant with the surface model geological information; Be used to control the asynchronous area diffusion; Generate the node on the parameter field; Make grid density degree on the corresponding curved surface of organ-tissue along with model structure characteristic and geometric properties variation, visually can present organ-tissue more realistically and receive external force to make the deformation effect of time spent.
Specific embodiment described herein only is that the present invention's spirit is illustrated.Person of ordinary skill in the field of the present invention can make various modifications or replenishes or adopt similar mode to substitute described specific embodiment, but can't depart from spirit of the present invention or surmount the defined scope of appended claims.

Claims (2)

1. the medical tissue dynamics simulation method based on the asynchronous diffusion model of power is characterized in that, may further comprise the steps:
Step 1 is set up geometric model, physical model and transition model for three-dimensional soft tissue, and the mode of specifically setting up comprises the geometric model of setting up three-dimensional soft tissue surfaces; Set up physical model through the subdivision geometric model, physical model is made up of voxel, and voxel is made up of the spring between particle and the particle; Set up transition model physical model is mapped to geometric model;
Step 2 by the inside and outside voxel of regional broadcast algorithm differentiation physical model, is removed the outside voxel of physical model, keeps the boundary voxel and the inner voxel of physical model;
The concrete mode of said regional broadcast algorithm does, at first in physical model, gets an inner voxel arbitrarily, is called seed voxels, and the zone that directly links to each other of definition voxel is the zone of the six direction of front and back up and down of this voxel; For seed voxels, extend toward its zone that directly links to each other, if its some voxel that directly link to each other is not a boundary voxel, then must be inner voxel; For newfound inner voxel, extend toward its zone that directly links to each other, if its some voxel that directly link to each other is not a boundary voxel, then must be inner voxel, by that analogy, up to extending to all boundary voxel of physical model;
Step 3 adopts two-stage collision mode to carry out crash tests, confirms collision detection point position according to the excellent position of operation; Said two-stage collision mode comprises finds the solution the regional area that bumps earlier, confirms the position of collision detection point at regional area again;
Step 4 is found the solution the Gaussian curvature of every bit on the geometric model;
Step 5 according to step 4 gained Gaussian curvature, is radiation radius of every bit definition on the geometric model; After soft tissue applied power; The asynchronous diffusion process of physical model stress and deformation is to be that basic point begins toward external diffusion with step 3 gained collision detection point; Diffuse to for the first time the point in the radiation areas of collision detection point, then the point in the radiation areas is all added the zone of action; Zone of action after the diffusion is as first asynchronous area for the first time; For the point of new adding zone of action, in diffusion for the second time, with them as basic point toward external diffusion; If a point is inner in their radiation areas; Then this point is joined the zone of action, the newly-increased point of all basic points has constituted second asynchronous area, by that analogy; All particles up to physical model all add the zone of action, have only the particle in the zone of action to move; The radiation areas of basic point obtain according to the radiation radius of respective point on the geometric model on the said physical model;
Step 6 according to the asynchronous area of step 5 gained physical model, is found the solution the position of each particle in asynchronous diffusion process of physical model, and the deformation effect on the physical model is embodied on geometric model through transition model and draws; The mode that the particle postition is found the solution does; On the physical model spring between the particle is being classified; According to spring classification analysis particle stressing conditions, analyze draw that particle is suffered and make a concerted effort after, utilize particle current time and last a position constantly to find the solution next position constantly of particle.
2. according to claim 1 based on the medical tissue dynamics simulation method of the asynchronous diffusion model of power, it is characterized in that: in the step 1,
The geometric model of setting up comprises triangular apex set
Figure 575451DEST_PATH_IMAGE001
and triangle surface set
Figure 2011102871272100001DEST_PATH_IMAGE002
; Wherein
Figure 676612DEST_PATH_IMAGE003
is the set of triangular apex index;
Figure 2011102871272100001DEST_PATH_IMAGE004
is the triangular apex coordinate set;
Figure 394033DEST_PATH_IMAGE005
is the set of triangle surface index, the index value on
Figure 2011102871272100001DEST_PATH_IMAGE006
expression Atria summit;
The physical model of setting up comprises particle set
Figure 976193DEST_PATH_IMAGE007
and set of voxels
Figure 2011102871272100001DEST_PATH_IMAGE008
; Wherein
Figure 790565DEST_PATH_IMAGE009
is the three-dimensional index set of particle;
Figure 550711DEST_PATH_IMAGE004
is the triangular apex coordinate set;
Figure 2011102871272100001DEST_PATH_IMAGE010
is the three-dimensional index set of voxel,
Figure 4695DEST_PATH_IMAGE011
be the three-dimensional index set on eight summits of voxel;
The transition model of setting up comprises voxel-vertex set
Figure 2011102871272100001DEST_PATH_IMAGE012
, summit-set of voxels
Figure 141278DEST_PATH_IMAGE013
and voxel-triangle set ; Wherein is the three-dimensional index set of voxel;
Figure 990471DEST_PATH_IMAGE015
is the relative position of surface vertices in voxel;
Figure 744800DEST_PATH_IMAGE003
is the set of triangular apex index;
Figure 419495DEST_PATH_IMAGE005
is the set of triangle surface index;
Figure 2011102871272100001DEST_PATH_IMAGE016
is triangular apex number in the voxel, and
Figure 132761DEST_PATH_IMAGE017
is the triangle number relevant with voxel.
CN201110287127A 2011-09-26 2011-09-26 Medical tissue dynamic simulation method based on force asynchronous diffusion model Pending CN102314710A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201110287127A CN102314710A (en) 2011-09-26 2011-09-26 Medical tissue dynamic simulation method based on force asynchronous diffusion model

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201110287127A CN102314710A (en) 2011-09-26 2011-09-26 Medical tissue dynamic simulation method based on force asynchronous diffusion model

Publications (1)

Publication Number Publication Date
CN102314710A true CN102314710A (en) 2012-01-11

Family

ID=45427846

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201110287127A Pending CN102314710A (en) 2011-09-26 2011-09-26 Medical tissue dynamic simulation method based on force asynchronous diffusion model

Country Status (1)

Country Link
CN (1) CN102314710A (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106202247A (en) * 2016-06-30 2016-12-07 哈尔滨理工大学 A kind of collision checking method based on longitude and latitude
CN106781941A (en) * 2016-11-24 2017-05-31 北京理工大学 A kind of method and its system for simulating microtrauma puncture operation
CN108961907A (en) * 2018-08-17 2018-12-07 深圳先进技术研究院 Virtual micro- ophthalmologic operation training method and system

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050018885A1 (en) * 2001-05-31 2005-01-27 Xuesong Chen System and method of anatomical modeling
CN1975784A (en) * 2006-12-28 2007-06-06 上海交通大学 Point particle spring deformation simulating method based on skeleton linear net
CN102044086A (en) * 2010-11-30 2011-05-04 华北水利水电学院 Soft tissue deformation simulation method

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050018885A1 (en) * 2001-05-31 2005-01-27 Xuesong Chen System and method of anatomical modeling
CN1975784A (en) * 2006-12-28 2007-06-06 上海交通大学 Point particle spring deformation simulating method based on skeleton linear net
CN102044086A (en) * 2010-11-30 2011-05-04 华北水利水电学院 Soft tissue deformation simulation method

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
WEIXIN SI ET AL.: "3D soft tissue warping dynamics simulation based on force asynchronous diffusion model", 《COMPUTER ANIMATION AND VIRTUAL WORLDS》 *
袁志勇等: "基于改进质点-弹簧模型的图像变形仿真方法", 《华中科技大学学报(自然科学版)》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106202247A (en) * 2016-06-30 2016-12-07 哈尔滨理工大学 A kind of collision checking method based on longitude and latitude
CN106781941A (en) * 2016-11-24 2017-05-31 北京理工大学 A kind of method and its system for simulating microtrauma puncture operation
CN108961907A (en) * 2018-08-17 2018-12-07 深圳先进技术研究院 Virtual micro- ophthalmologic operation training method and system

Similar Documents

Publication Publication Date Title
Delingette et al. Craniofacial surgery simulation testbed
Je et al. PolyDepth: Real-time penetration depth computation using iterative contact-space projection
Qian et al. Virtual reality based laparoscopic surgery simulation
TW201013595A (en) A haptic response simulation method and system for real-time haptic and imaging responses
Saha et al. Fuzzy object skeletonization: theory, algorithms, and applications
Bornik et al. A hybrid user interface for manipulation of volumetric medical data
Mirhosseini et al. Benefits of 3D immersion for virtual colonoscopy
Gonçalves et al. A benchmark study on accuracy-controlled distance calculation between superellipsoid and superovoid contact geometries
CN103345774B (en) A kind of modeling method of three-dimensional multi-scale vector quantization
Wei et al. Centerline extraction of vasculature mesh
Fukuhara et al. Proposition and evaluation of a collision detection method for real time surgery simulation of opening a brain fissure
CN102314710A (en) Medical tissue dynamic simulation method based on force asynchronous diffusion model
Yi et al. The implementation of haptic interaction in virtual surgery
Rossignac et al. HelSweeper: Screw-sweeps of canal surfaces
Vlasov et al. Haptic rendering of volume data with collision detection guarantee using path finding
Selmi et al. 3D interactive ultrasound image deformation for realistic prostate biopsy simulation
Noborio et al. Comparison of GPU-based and CPU-based Algorithms for Determining the Minimum Distance between a CUSA Scalper and Blood Vessels.
Su et al. Haptic‐based virtual reality simulator for lateral ventricle puncture operation
Chen et al. Computer-aided liver surgical planning system using CT volumes
Bogoni et al. Haptic technique for simulating multiple density materials and material removal
Tsai et al. Accurate visual and haptic burring surgery simulation based on a volumetric model
Yureidini et al. Local implicit modeling of blood vessels for interactive simulation
Wang et al. A prediction method for assembly surface contact considering form error
Audette et al. A review of mesh generation for medical simulators
Yan et al. Real-time bone sawing interaction in orthopedic surgical simulation based on the volumetric object

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C12 Rejection of a patent application after its publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20120111