CN106844994A - The CNV growth prediction method that constitutive model is combined with finite element - Google Patents

Abstract

The invention discloses the CNV growth prediction method that a kind of constitutive model is combined with finite element, image preprocessing；Extracted region and division：Divide the image into is CNV regions, outer retinal layers, 4 regions of interior layer of retina and choroid layer；Gridding：Tetrahedral grid generation is carried out to 4 regions；Modeling：Modeled with super-elasticity biomechanical model and reaction diffusion equation, the quality change after CNV is grown is added in equation as source item, deformation gradient is persistently changed according to the growth of new vessels；Optimized model, calculates optimal accuracy rate, carries out parametric test；One parameter curve of parameter fitting according to each time point prediction, predicts the growth parameter(s) at last time point, is predicted the outcome.The method of the present invention can more flexibly and personalized model organism mechanical model, and the tissue of model hypothesis is orthotropy, also has to non-linear, large deformation region and predicts the outcome well, and the degree of accuracy is high.

Description

The CNV growth prediction method that constitutive model is combined with finite element

Technical field

The present invention relates to computer vision, image processing and analysis technical field, belong to the modeling method of growth prediction, especially It is the modeling side of CNV (CNV) growth prediction for being applied to SD-OCT (domain optical coherence fault imaging) Method.

Background technology

Choroid is located between retina and sclera, is made up of fibr tissue, thin vessels and capillary.Choroidal blood Circulation nutrition ectoretina, crystalline lens and vitreum etc., because CBF is big, itself is soft compared with slow for flow velocity in choroid And thin characteristic, make it influence or impact not tolerating to external world, cause pathogen herein easily be detained, cause choroidal diseases and Various fundus oculi diseases, therefore, the growth of research CNV (CNV) can obtain many important data, be follow-up Medical research provides basic support data.

Because the growth site of new vessels is in the extremely sensitive and fragile ophthalmic retina of environment to external world, enter in real time Being detected to new vessels for invading property is clearly impossible.Existing neovascularization growth Forecasting Methodology is generally only with list One model of growth such as FInite Element is predicted, and this method has certain shortcomings and limitations：1. it is based on isotropism It is assumed that 2. do not account for the different biological characteristic of different tissues structure, 3. can only predict linear change, 4. can only predict The change of small deformation, true to the forecasting inaccuracy of large deformation change, etc..

The content of the invention

The technical problems to be solved by the invention are the defects for overcoming prior art, present invention firstly provides a synthesis Method --- the method that biomechanical model and reaction diffusion equation based on elastic material are combined of property, carries out choroid The growth prediction analysis of new vessels, the degree of accuracy is high.

In order to solve the above technical problems, the present invention provides a kind of CNV three-dimensional generation Forecasting Methodology.

The method mainly includes 6 steps：

Step 1. image preprocessing：Down-sampling is carried out to three-dimensional retina domain optical coherence fault imaging SD-OCT images Treatment；

Step 2. extracted region and division：Divide the image into as CNV CNV regions, outer retinal layers, interior 4 regions of layer of retina and choroid layer；

Step 3. gridding：CNV regions, interior layer of retina, outer retinal layers and choroid layer are carried out multiple dimensioned, adaptive The tetrahedral grid answered is generated；

Step 4. is modeled with super-elasticity biomechanical model with reaction diffusion equation：In super-elasticity biomechanical model knot In structure, integrate the different bulk modulus in each region, modulus of shearing and its response of corresponding volume elasticity, etc. hold elastic response, will Quality change after CNV growth is added in equation as source item, makes deformation gradient according to new vessels Growth persistently change；

Step 5. Optimized model, calculates optimal accuracy rate, carries out parametric test；

Step 6. predicts the life at last time point according to one parameter curve of parameter fitting of each time point prediction Appreciation rate parameter long, is predicted the outcome.

Image preprocessing in step 1 is comprised the following steps：

The treatment of (1a) down-sampling：

(1b) uses CAVASS softwares, carries out OCT image registration.

Using the mesh generator of a three-dimension curved surface and volume to CNV regions, interior layer of retina, outer view in step 3 Film layer and choroid layer carry out multiple dimensioned, self adaptation tetrahedral grid generation.

Step 4 is comprised the following steps：

The invasion and diffusion of (4a) reaction diffusion equation

The reaction diffusion equation that logistic regression type increases is expressed as：

Wherein, t is the time, and D is anisotropic diffusion tensor, comprising three part D_x, D_y, D_z, these three components Respectively x-axis, y-axis, the component of z-axis；ρ represents growing multiplication rate, and N is the cell number of per unit volume, K_maxRepresent maximum cell Capacity,

When the anisotropy concentration of VEGFWhen, there is equation below：

(4b) mass effect and super-elasticity biomechanical model

The strain energy density function of the non-linear stress-strain behavior of material is：

Wherein,Represent Green-Lagrange strain tensor, F^TRepresent the transposition of deformation gradient F Amount, I is unit matrix；J represents the determinant of deformation gradient F, i.e. J=det (F)；κ and μ represent respectively bulk modulus and Modulus of shearing；The mark of Tr representing matrixs；I₁It is the first main invariant of Cauchy-Green's Deformation tensor；

By the gradient force f in formula 4, by cauchy stress tensor σ=J^-1FSF^TWith VEGF it is each to Different in nature concentration c is connected：

Wherein, it is density gradient power that f represents the gradient force produced by c, and ξ represents one needs estimative constant, and this is normal Number depends on the biological characteristics of material；S represents the Piola-Kirchhoff stress tensors of PK2- the 2nd.

In step 5, with single order Markov sequence Optimal Growing appreciation rate ρ parameters：The data at preceding n time point, each Time point is related to a time point before, calculates the value of ρ.

From all of time point before grow the estimation of appreciation rate parameter, then calculated using following formula pre- The accuracy of survey：

It is all time point I based on before_i,c(10) i=1 2 ..., predicts the result that obtains, I_i+1,cWhen being current Between put goldstandard；

Finally, parameter is estimated：ν_i={ D_x,D_y,D_z,ρ；c}.

The coincident pixel point that definition predicts the outcome with goldstandard

Calculate the indices for reflecting the accuracy that predicts the outcome：

(1) recall rate

(2) error rate

(3) similarity DICE coefficients

In step 6, using the growth appreciation rate ρ values at the preceding n time point for drawing, matched curve is drawn, so as to predict (n+1)th growth appreciation rate ρ at time point；The growth appreciation rate ρ at (n+1)th time point that will be predicted substitutes into n-th time Point, the prognostic chart picture for drawing is Chong Die with real (n+1)th image, calculates predictablity rate.

The beneficial effect that the present invention is reached：

Present invention firstly provides a comprehensive growth prediction method --- the biomethanics based on elastic material The method that model and reaction diffusion equation are combined, carries out the growth prediction analysis of CNV.

Method of the present invention combination elastic material and reaction response diffusion equation, being capable of more flexible and personalization mould Formula bio-mechanical model.Material different from the past is the assumption of isotropy, and the tissue of the model hypothesis is orthotropy. Disease location (CNV region) and around tissue (including outer retina and choroid, retinochrome Deng) the different biological characteristic of different tissues structure is simulated, and based on orthotropic it is assumed that to non-linear, big change Shape region also has and predicts the outcome well, and the degree of accuracy is high.

Brief description of the drawings

Fig. 1 is the entire block diagram of CNV prediction；

Fig. 2 a are the results after original image (512 × 1024 × 128) gridding of OCT；

Fig. 2 b are the results after image (64 × 64 × 64) gridding of down-sampling；

Fig. 2 c represent the nodal information and unit information of original image；

Fig. 2 d represent the nodal information and unit information of down-sampled images；

Fig. 3 a are the mesh generations from down-sampled data reconstruct；

Fig. 3 b are the intracellular volumes (ICV) of CNV；

Fig. 3 c are the iso-surface images of ICV；

Fig. 3 d are the displays of stress distribution：In each section, the thick lines in CNV regions are circular and hachure ellipse divides The position of strong stress (being typically distributed across CNV edges) and weak stress (being typically distributed across inside CNV) is not represented；

Fig. 4 is the CNV volume difference schematic diagrames of different patients.

Specific embodiment

Specific embodiment of the invention is described in further detail below：

The fundamental block diagram of this method is as shown in figure 1, mainly include the following steps that：

1) image preprocessing

The treatment of (a) down-sampling：In order to improve treatment effeciency, memory optimization is used, by each pixel 512* in the present embodiment The three-dimensional retina SD-OCT image down samplings of 1024*128 are to pixel 64*64*64.

B () carries out OCT image registration using CAVASS softwares (software that MIPG is developed).

2) extracted region and division

Hand labeled goldstandard, carries out extraction, surrounding tissue (the interior retina in retina area-of-interest (CNV regions) Layer, outer retinal layers, choroid layer) division；

3) gridding, with the mesh generator of a three-dimension curved surface and volume to CNV regions, interior layer of retina, outer view Film layer and choroid layer carry out multiple dimensioned, self adaptation tetrahedral grid generation, see Fig. 2 a, Fig. 2 b, Fig. 2 c, Fig. 2 d；

4) to CNV regions, interior layer of retina, outer retinal layers and choroid layer this 4 region modelings, wherein regional Carry out different super-elasticity biomechanical models carries out parameter training with reaction diffusion equation model, makes it according to different physiology Characteristic is learned to be grown：

The invasion and diffusion of (a) reaction diffusion equation

The reaction diffusion equation that logistic regression type (Logistic) increases can be expressed as：

Wherein, t is the time, and D is anisotropic diffusion tensor, comprising three part D_x, D_y, D_z, these three components Respectively x-axis, y-axis, the component of z-axis.ρ represents growing multiplication rate, and N is the cell number of per unit volume, K_maxRepresent maximum cell Capacity, when the anisotropy concentration of VEGF (VEGF)When, there is equation below：

(b) mass effect and super-elasticity biomechanical model

In model structure, by equation (3) bulk modulus more different than more fully integrating each region, modulus of shearing and Its in constitutive equation corresponding volume elasticity response part, etc. hold elastic response part, the quality after CNV is grown changes to be made For source item is added in the constitutive equation of super-elasticity biomethanics, deformation gradient is set persistently to be become according to the growth of new vessels Change：

Hyperelastic biomechanical model can be used for the prediction of moderate finite deformation, the non-linear stress-strain behavior of material Strain energy density function is：

Wherein,Represent Green-Lagrange strain tensor, F^TRepresent the transposition of deformation gradient F Amount, its component relation is (F^T)_ij=F_ji, I is unit matrix；J represents the determinant of deformation gradient F, i.e. J=det (F)；κ Represent bulk modulus and modulus of shearing respectively with μ；The mark of Tr representing matrixs；I₁It is Cauchy-Green (Cauchy-Green) deformation First main invariant of tensor.

Stress distribution on the mesh generation of data reconstruction and each section, is shown in Fig. 3 a, Fig. 3 b, Fig. 3 c, Fig. 3 d.

Meanwhile, the process is simulated using static balancing equation：

The equation passes through gradient force f, by cauchy stress tensor σ=J^-1FSF^T(wherein S represents the Piola- of PK2- the 2nd Kirchhoff stress tensors), and the anisotropy concentration c of VEGF connects：

Wherein, f represents the gradient force (i.e. density gradient power) produced by c, and ξ represents one needs estimative constant, should Constant depends on the biological characteristics of material.

5) Optimized model, calculates optimal accuracy rate, carries out parametric test.

With single order Markov sequence Optimal Parameters ρ (growth appreciation rate)：Preceding n (n=10 is set in the present embodiment) the individual time The data of point, each time point is related to a time point before, calculates the value of ρ.

The estimation of model parameter (growth appreciation rate) is carried out from all of time point before, is then come using following formula Calculate the accuracy of prediction：

It is all time point I based on before_i,c(i=1 2 ..., n) predicts the result that obtains, I_i+1,cIt is current time The goldstandard of point, w₁、w₂Weight is represented respectively.

According to the physiological property of CNV, the parameter κ of CNV, choroid (choroid) and other parts is set It is respectively with μ：κ_CNV=5kPa, κ_choroid=6kPa, κ_others=3kPa, μ_CNV=μ_choroid=5*10³N/m²,μ_others=1* 10³N/m².Finally, our parameters to be estimated are：ν_i={ D_x,D_y,D_z, ρ；c}.The picture that overlaps that definition predicts the outcome with goldstandard Vegetarian refreshmentsThe indices for reflecting the accuracy that predicts the outcome are calculated, these indexs are just as evaluation and testing model Standard：(1) recall rate(2) error rate(3) similarity DICE coefficients

6) one parameter curve of the parameter fitting according to each time point prediction, predicts the growth ginseng at last time point Number, is predicted the outcome.

Using the growth appreciation rate ρ values at preceding n (n=10 is set in the present embodiment) the individual time point for drawing, matched curve is drawn, So as to predict (n+1)th growth appreciation rate ρ at time point.The growth appreciation rate ρ at (n+1)th time point that will be predicted is substituted into N-th time point, the prognostic chart picture for drawing is Chong Die with real (n+1)th image, calculates predictablity rate.

7) experimental result：

As shown in Figure 4, it is 10 CNV volumes of the six of time point patients, it can be seen that CNV's is special Property substantially, fluctuate it is different.In fact, the CNV of different patient's different time points is no matter in volume size, or in form It is all widely different.

Because retina is to the sensitiveness of light, the display (seeing below accompanying drawing 3d) of stress distribution is also increased in the result, have Help more intuitively show influence of the mechanics to CNV regions.

According to the deformation in first time point, the second time point parameter of prediction；According to the deformation at the first two time point, in advance Survey the 3rd deformation parameter at time point.By that analogy, the song for the growth appreciation rate ρ at preceding n (n=10) individual time point being fitted Line, substitutes into n-th time point, draws and predicts the outcomeIt is compared with actual value (goldstandard at n+1 time points), is drawn pre- Survey result, such as table 1.

Table 1 is predicting the outcome for this method, it is found that except the huge disease of indivedual CNV surrounding times point variations Outside people, remaining predicts the outcome and is obtained for the relatively good degree of accuracy.

So far, a kind of three-dimensional CNV growth prediction models have been carried out and are verified.The present invention proposes a kind of fusion The thought of super-elasticity biomechanical model, reaction diffusion equation and Mass Sources, the Optimal Parameters automatic detection of combined objective function, And the technical scheme of region-wide stress distribution is shown in the result, test result indicate that, growth prediction method energy of the invention Enough relatively good growth predictions for solving the problems, such as CNV (CNV), with prediction accuracy higher.

The above is only the preferred embodiment of the present invention, it is noted that for the ordinary skill people of the art For member, on the premise of the technology of the present invention principle is not departed from, some improvement and deformation can also be made, these improve and deform Also should be regarded as protection scope of the present invention.

Claims (8)

1. the CNV growth prediction method that a kind of constitutive model is combined with finite element, it is characterized in that, including it is following Step：

Step 1. image preprocessing：Three-dimensional retina domain optical coherence fault imaging SD-OCT images are carried out at down-sampling Reason；

Step 2. extracted region and division：Divide the image into is CNV CNV regions, outer retinal layers, interior view 4 regions of film layer and choroid layer；

Step 3. gridding：CNV regions, interior layer of retina, outer retinal layers and choroid layer are carried out multiple dimensioned, self adaptation Tetrahedral grid is generated；

Step 4. is modeled with super-elasticity biomechanical model with reaction diffusion equation：In super-elasticity biomechanical model structure In, integrate the different bulk modulus in each region, modulus of shearing and its response of corresponding volume elasticity, etc. hold elastic response, by arteries and veins Quality change after network film neovascularization growth is added in equation as source item, makes deformation gradient according to new vessels The lasting change of growth；

Step 5. Optimized model, calculates optimal accuracy rate, carries out parametric test；

Step 6. predicts that the growth at last time point increases according to one parameter curve of parameter fitting of each time point prediction Value rate parameter, is predicted the outcome.

2. the CNV growth prediction method that constitutive model according to claim 1 is combined with finite element, its It is characterized in that the image preprocessing in step 1 is comprised the following steps：

The treatment of (1a) down-sampling：

(1b) uses CAVASS softwares, carries out OCT image registration.

3. the CNV growth prediction method that constitutive model according to claim 1 is combined with finite element, its It is characterized in, using the mesh generator of a three-dimension curved surface and volume to CNV regions, interior layer of retina, outer retina in step 3 Layer and choroid layer carry out multiple dimensioned, self adaptation tetrahedral grid generation.

4. the CNV growth prediction method that constitutive model according to claim 1 is combined with finite element, its It is characterized in that step 4 is comprised the following steps：

The invasion and diffusion of (4a) reaction diffusion equation

The reaction diffusion equation that logistic regression type increases is expressed as：

$\frac{\∂N}{\∂t}=div(D\▿N)+\mathrm{\ρ}N(1-\frac{N}{K_{max}})---\left(1\right)$

Wherein, t is the time, and D is anisotropic diffusion tensor, comprising three part D_x, D_y, D_z, these three components difference It is x-axis, y-axis, the component of z-axis；ρ represents growing multiplication rate, and N is the cell number of per unit volume, K_maxRepresent maximum cell appearance Amount,

When the anisotropy concentration of VEGFWhen, there is equation below：

$\frac{\∂c}{\∂t}=\▿\·\left(D\▿c\right)+\mathrm{\ρ}c\left(1-c\right)---\left(2\right)$

(4b) mass effect and super-elasticity biomechanical model

The strain energy density function of the non-linear stress-strain behavior of material is：

$\mathrm{\Ψ}\left(\mathrm{\ϵ}\right)=\frac{1}{2}\mathrm{\κ}{\left(J-1\right)}^2+\frac{1}{2}\mathrm{\μ}Tr(I_1-3)---\left(3\right)$

Wherein,Represent Green-Lagrange strain tensor, F^TRepresent the transposition tensor of deformation gradient F, I It is unit matrix；J represents the determinant of deformation gradient F, i.e. J=det (F)；κ and μ represent bulk modulus and shearing respectively Modulus；The mark of Tr representing matrixs；I₁It is the first main invariant of Cauchy-Green's Deformation tensor；

By the gradient force f in formula 4, by cauchy stress tensor σ=J^-1FSF^TWith the anisotropy of VEGF Concentration c is connected：

$div\left(\mathrm{\σ}\right)+f=0;f=-\mathrm{\ξ}\▿c---\left(4\right)$

Wherein, it is density gradient power that f represents the gradient force produced by c, and ξ represents one needs estimative constant, the constant to take Certainly in the biological characteristics of material；S represents the Piola-Kirchhoff stress tensors of PK2- the 2nd.

5. the CNV growth prediction method that constitutive model according to claim 1 is combined with finite element, its It is characterized in, in step 5, with single order Markov sequence Optimal Growing appreciation rate ρ parameters：The data at preceding n time point, each when Between put it is related to a time point before, calculate ρ value.

6. the CNV growth prediction method that constitutive model according to claim 5 is combined with finite element, its It is characterized in from all of time point before grow the estimation of appreciation rate parameter, is then calculated using following formula pre- The accuracy of survey：

$R\left({\mathrm{\ν}}_i\right)=\underset{i}{\mathrm{\Σ}}{w}_1(1-\frac{{\stackrel{\‾}{I}}_{i,c}\∩I_{i+1}}{I_{i+1,c}})+w_2\left(\frac{{\stackrel{\‾}{I}}_{i,c}-{\stackrel{\‾}{I}}_{i,c}\∩I_{i+1}}{I_{i+1,c}}\right)---\left(5\right)$

It is all time point I based on before_i,c(10) i=1 2 ..., predicts the result that obtains, I_i+1,cIt is current point in time Goldstandard；

Finally, parameter is estimated：ν_i={ D_x,D_y,D_z,ρ；c}.

7. the CNV growth prediction method that constitutive model according to claim 6 is combined with finite element, its It is characterized in, the coincident pixel point that definition predicts the outcome with goldstandard

Calculate the indices for reflecting the accuracy that predicts the outcome：

(1) recall rate

(2) error rate

(3) similarity DICE coefficients

8. the CNV growth prediction side that the constitutive model according to claim 5,6 or 7 is combined with finite element Method, it is characterized in that, in step 6, using the growth appreciation rate ρ values at the preceding n time point for drawing, matched curve is drawn, so that in advance Measure (n+1)th growth appreciation rate ρ at time point；The growth appreciation rate ρ at (n+1)th time point that will be predicted is substituted into n-th Time point, the prognostic chart picture for drawing is Chong Die with real (n+1)th image, calculates predictablity rate.