[Title of Invention]

METHOD FOR SIMULATING THORACIC 4DCT

[Technical Field]

[0001]

The invention relates generally to four-dimensional computed tomography (4DCT), and more particularly to simulating respiration for radiotherapy treatment planning.

[Background Art]

[0002]

The use of four-dimensional computed tomography (4DCT) is a common practice in radiation therapy treatment planning for thoracic regions. The information from 4DCT allows a radiation oncologist to plan accurate treatments for moving tumors, deliver radiation within predetermined certain interval in the breathing cycle, and reduce the risk of treatment-related side effects.

[0003]

Methods for 4DCT include retrospective slice sorting, and prospective sinogram selection. For the retrospective slice sorting method, the projection data are continuously acquired for a time interval longer than a respiratory cycle. Multiple slices corresponding to different times are reconstructed and sorted into respiratory phase bins using various respiratory signals. For the prospective sinogram selection method, the scanner is triggered by the respiratory signal. Then, the projection data within the same phase bin are used to reconstruct CT slices corresponding to that breathing phase.

[0004]

Both methods are time consuming and considerably increase the radiation dose for the patient. For example, the radiation dose of a conventional 4DCT scan is about six times a typical helical CT scan, and 500 times a chest X-ray. Moreover, 4DCT acquisition alone cannot determine the tumor position in-situ. These facts are a major concern in the clinical application of 4DCT, motivating development of advanced 4DCT simulators.

[0005]

Various methods are known for modeling lung inflation and deflation. The first category of methods discretize the soft tissues and bones into masses (nodes), and connect the nodes using springs and dampers (edges) based on a mass-spring-damper system, and CT scan values for spline-based mathematical cardiac-torso (MCAT) phantoms, augmented reality based medical visualization, respiration animation, and tumor motion modeling. A common approach applies affine transformations to control points to simulate respiratory motion. Lungs and body outline are linked to the surrounding ribs for synchronized expansion and contraction. Those simple and fast methods only provide approximate lung deformations when compared with accurate models based on continuum mechanics.

[0006]

The second category of methods use hyperelastic models to describe the non-linear stress-strain behavior of the lung. To simulate lung deformation between two breathing phases (7 ,7 ₊₁ ), the lung shape at phase 7 ₊₁ is used as the contact or constraint surface, and deform the lung at phase 2} based on the predefined mechanical properties of lung. A negative pressure load is applied on the lung surface, and finite element (FE) analysis is used to deform tissues. The lung expands according to the negative pressure, and slide against the contact surface to imitate the pleural fluid mechanism. The pressure can be estimated from the pleural pressure versus lung volume curve, which is measured from a pulmonary compliance test.

[0007]

In addition to lung deformation, the displacements of the rib cage and diaphragm are also very important to design a realistic 4DCT simulator. Rib cage motion can be a rigid transformation. A finite helical axis method can be used to simulate the kinematic behavior of the rib cage. The method can be developed into a chest wall model relating the ribs motion to thorax-outer surface motion for lung simulation.

[0008]

A simple diaphragm model includes central tendon and peripheral muscular fiber. Then, cranio-caudal (CC) forces are applied to each node of the muscular fiber to simulate the contraction of the diaphragm. A Gauchy-Green deformation tensor can model the lung deformation. Organs inside the rib cage can be considered as a convex balloon to estimate an internal deformation field directly by interpolation of skin marker motions.

[Summary of Invention]

[0009]

The embodiments of the invention provide a method for simulating thoracic 4DCT. The method uses a biomechanical motion model that can faithfully simulate deformation of the lung and nearby organs for an entire respiration cycle. Deformations of other organs can also be simulated.

[0010]

As an advantage, the method only requires a single CT scan as input. Loads on the rib cage and the diaphragm constrain the deformation of the lung. This distinguishes the present method from conventional continuum mechanics based methods. The method can also simulate a passive mass-spring model based deformation of abdominal organs due to respiration.

[0011]

The lung stress-strain behavior is modeled using a hyperelastic model. The lung deformation problem is solved by finite element (FE) analysis. Diaphragm control points and spine regions are segmented out to carefully define the boundary conditions and the load definitions, and to improve FE convergence using a mesh optimization procedure.

[0012]

The remaining CT volume is treated as discretized mass points connected by springs and dampers. The motion of liver, bones and other organs can be simulated using finite difference analysis.

[0013]

This heterogeneous design can leverage the advantages of both continuum mechanics and mass-spring-damper system in the way that the lung deformation is modeled with a very high accuracy, while the deformation of the rest of the CT volume is achieved under practical computation constraints.

[Brief Description of the Drawings]

[0014]

[Fig. 1]

Fig. 1 is a flow diagram of a method for simulating thoracic 4DCT according to embodiments of the invention;

[Fig. 2]

Fig. 2 is a flow diagram of a method for simulating thoracic 4DCT according to embodiments of the invention;

[Fig- 3] Fig. 3 is a sample 2D axial view for a maximized region covered a Gaussian function according to embodiments of the invention;

[Fig. 4]

Fig. 4 is an image of a dense 3D point cloud for lung voxels according to embodiments of the invention;

[Fig. 5]

Fig. 5 is an image of the lung after converting the point of Fig. 4 cloud into a weight map according to embodiments of the invention;

[Fig. 6]

Fig. 6 is a sample weight map according to embodiments of the invention; [Fig. 7]

Fig. 7 is a schematic of sample blocks on the lung surface for a weight calculation procedure according to embodiments of the invention;

[Fig. 8]

Fig. 8 is a schematic of masking an inner region of the lung from an outer boundary according to embodiments of the invention;

[Fig. 9]

Fig. 9 is a schematic of a ribbon formed after the masking of Fig. 8 according to embodiments of the invention;

[Fig. 10]

Fig. 10 is a schematic of an irregular surface mesh according to embodiments of the invention; and

[Fig. 11]

Fig. 11 is a schematic of an optimized surface mesh according to embodiments of the invention.

[Description of Embodiments]

[0015] Thoracic 4DCT Simulator Method

Figs. 1 and 2 show a processing pipeline of a thoracic four-dimensional computed tomography (4DCT) simulator method according to embodiments of the invention. The method includes lung deformation and 4DCT scan simulation. The method can be performed in a processor connected to memory and input/output interfaces as known in the art. It should be understood that the invention can also be used to simulate other types of tissue or organ deformations, and utilize other types of medical scans known in the art.

[0016]

Input to the method is a single thoracic CT scan 101. The CT scan is segmented 110 to generate a surface mesh 111 of the lung. Tetrahedralization 120 is applied to the surface mesh to obtain a first volume mesh 121.

[0017]

Finite element(FE) analysis 130, using boundary constraints and load definitions 138, are applied to the first volume mesh to obtain a lung deformation 131. The lung deformation is according to an Ogden model, wherein the lung is modeled using a hyperelastic function based on nonlinear continuum mechanics, such that the tumor motion trajectories inside the lung are accurate.

[0018]

Continuing with in Fig. 2, constrained tetrahedralization 150, using control points 149, is applied to the lung deformation 131 to obtain a second volume mesh 151, which is a refinement of the first volume mesh. The refined mesh coincides with or is partly or entirely outside the first volume mesh. In other words the second mesh corresponds to a refined organ volume. That is, the refined organ volume accurately represent the actual shape of the organ during the simulation.

[0019] Then, mass-spring-damper simulation 160 produces the simulated thoracic 4DCT 161. The mass-spring-damper model is used to deform the CT volume. Thus, segmentation of rig cage, spines, and other organs, as in the prior art, is unnecessary.

[0020]

Our simulator can generate any desired number of lung deformations between a maximal exhale phase and the maximal inhale phase T _in , which are used to synthesize the corresponding 4DCT phases.

[0021]

The method can construct 4DCT phases with or without hysteresis of the tissue deformation, which defines how the lung volume is correlated with the inhale-to-exhale motion path and the exhale-to-inhale motion path. Without hysteresis, the CT scan for the same lung volume states are identical regardless of whether the lung is in the inhale-to-exhale cycle or the inhale-to-exhale cycle.

[0022]

For a particular lung volume, our method assembles the corresponding 4DCT phases to generate a simulated sequence of CT's. In addition to the lung volume, pleural pressure and chest and abdomen height can also be used to synchronize for real-time patient-specific 4DCT, as described below.

[0023]

The method steps described herein can be performed in a processor connected to memory and input/output interfaces as known in the art.

[0024]

Biomechanical Motion Model for Simulation of Lung Deformation

Boundary Constraints

For simplicity of notation, we use x, y, and z to represent lateral, anterioposterior (AP), and superoinferior (SI) direction, respectively. For the lung deformation, we define three types of boundary constraints on the lung surface based on the anatomy and function of the human respiratory system. Because the upper lobes of the lung are constrained by the ribs, the z displacement of the corresponding region are fixed in our simulator to avoid pure translation of the lung in the z direction when simulating the diaphragm contracting on the bottom lobes. Intuitively, the top 5% of the surface vertices in the z direction are selected to represent the upper lobes.

[0025]

During inspiration, the lung sliding against the rib cage mainly occurs in the posterior spine region, while in the anterior region, the lung expands with the increasing of thoracic cavity and the relative sliding is negligible. Thus, we define the boundary conditions for both the front and the back parts of the lung surface to simulate the different sliding actions. Our method fixes the z displacement for all the points, generally 139 in Fig. 1, to simulate the coherent motion of the lung with the thorax expansion on the axial plane. The selection of these vertices is implemented by the heuristics that the points are on or near a convex hull of the lung surface, around the lateral sides of the middle and lower lobes, and have small normal variations, e.g., < 20° .

[0026]

To simulate the pleural sliding in the spine region, our simulator method automatically locates the lung surface vertices near the thoracic vertebrae, and fixes the x and y displacements of these points as the third boundary constraint. Our goal is to find surface vertices close to the spine. Therefore, we use a Gaussian curve fitting process to locate these points instead of adopting a complicated thoracic vertebrae segmentation approach, as in the prior art.

[0027] The idea is to fit a set of Gaussian curves, such that the area cut out by each curve is maximized. This provides a global approximation to the spine shape, and the constraint points can be accurately located.

[0028]

Fig. 3 shows a sample 2D axial view, in which our procedure maximizes the spine region S covered by the Gaussian function

(x-b) ²

f (x) = ae ^2c

[0029]

We formulate the Gaussian function as a constrained multi-variable minimization problem

xmax

min - ∑ /( ), si. f (x) - g(x)≤0, Vx [x _min ,x _max ], a,b,c ^{x =x} _t nin

where the parameters a, b, and c represent a scaling factor, an expected value, and a standard variance of / (x) , x _min and x _max are limits on the lung in the lateral direction, and g(x) is an upper limit for / (x) , and is the minimal value of the lung slice at each x .

[0030]

In our simulator, this constrained minimization problem is solved efficiently by a sequential quadratic programming method, specifically an active-set procedure, which determines a sequential quasi-Newton approximation to the Hessian matrix of the Lagrangian multiplier at each iteration. We extend this 2D procedure to the 3D CT volume by applying the procedure slice by slice.

[0031] Outliers can occur in the top and bottom of the lung where g(x) is only partial constrained for the curve fitting. Our simulator removes the outliers different than a mean Gaussian curve of the set so that correct fittings of the thoracic vertebrae are retained. The missing curves can be estimated by linear interpolation of the remaining curves.

[0032]

Loads Definitions

Because there is no bounding surface during inhalation, we design an extra traction applied on the diaphragm area of the lung besides the negative intra-pleural surface pressure. The pressure force inflates the lung in all directions during inspiration, while the traction allows additional displacement in the z direction to simulate the diaphragm contraction and pleural sliding. Hence, the lung surface mesh is likewise inflated during the simulation.

[0033]

The pressure force is defined from the simulator input. Therefore, we focus on accurately locating the points that are close to the diaphragm for the definition of the traction. We model this as a graph search problem, and solve the problem by a modified shortest closed-path procedure.

[0034]

As shown in Fig. 4, our simulator first determines a dense 3D point cloud by finding the lung voxels at every (x, y ) location with the largest z value.

[0035]

As shown in Fig. 5, we convert the point cloud into a weight map based on the local geometry information, and locate the diaphragm points using our modified shortest closed-path procedure. The left and right lower lobes are treated separately.

[0036] Weight Map

We consider the 3D point cloud as an 2D image with intensity value from the z value of the corresponding point, and run a local line direction discrepancy (LDD) computation on this image to generate a weight map W. Thus, our weight map computation is a special type of image filtering.

[0037]

As shown in Fig. 6 for each line ^d i (x,y) ₀ f _a block centering at ( ^x ^ ), we

1 2 3 2 1 3 construct two sub-lines ^d i ( ^x> y) and ^d i from ( Pi , Pi , Pi ) and ( Pi ,

Pi , Pi ) respectively, - 1 ^> · · ·'4) _{} anc} j determine the LDD as the minimum intersection angle of the four sub-line pairs.

[0038]

Alternatively, we can determine the maximum of the cosine value of these angles to represent the wei ht, which can be determined by a dot product as

where p represents pixel position (x,y ), and the block size is 5 x 5 , for simplicity. Intuitively, regions with high curvature have high positive LDD values, for example, άγ of in Fig. 7, while flat regions have low negative LDD values, for example, ¾ and B ₂ .

[0039]

Diaphragm Point Segmentation

Notice that all outliers are located at the boundary of the weight map. Therefore, we formulate the diaphragm point segmentation as a shortest closed-path (SCP) problem, which finds a optimal cut along the boundary that separates the diaphragm points from the outliers. To construct the graph for SCP, we select four neighborhood connection and set the edge weight as W(q) .

Therefore, E _pq and can have different weights, instead of using the entire weight map to construct the graph.

[0040]

As shown in Fig. 8, we mask out the inner region with morphological operations, and limit an optimal cut between the inner boundary Ω ₂ and the outer boundary Ω . If we use the conventional SCP procedure, then some interior regions can be cut out to favor the lowest cost. Instead, we first sample the outer boundary 6Q _j every ten points, and locate the corresponding points, in terms of Euclidean distance, on the inter boundary <9Ω ₂ ·

[0041]

For the remaining points on the outer boundary 5Ω _{1 ?} we determine their matches on the inner boundary dQ ₂ by linear interpolation of the previous matches, such that there are no crossing matches (lines) and correct ordering is maintained. In this way as shown in Fig. 9, we can "unbend" the ring region between dQ _j and 5Ω ₂ into a flat ribbon by alignment in order, and set the length of the ribbon as the length of the outer boundary θΩ _1? and the width as the shortest distance between outer boundary ΘΩι and the inner boundary <9Ω ₂ . Then, we construct a new adjacency matrix graph from the ribbon for the SCP procedure. This yields an accurate model for the diaphragm for the traction definition.

[0042]

Finite Element Simulation

The final step for biomechanical simulation of lung deformation according to embodiments of the invention is to define the material property of the lung and apply the FE analysis. We assume the lung tissue is homogeneous, isotropic, and use a first-order Ogden model. The well-known Ogden model is a hyperelastic material model used to describe the non-linear stress-strain behavior of complex materials such as biological tissue. The Ogden model, like other hyperelastic material models, assumes that the material behavior can be described by means of a strain energy density function, from which the stress-strain relationships can be derived. These materials can generally be considered to be isotropic, incompressible and strain rate independent. Other models that can be used include a Mooney-Rivlin model, a linear elastic model, a linear visco-elastic model, or combinations thereof.

[0043]

The Ogden model describes a non-linear strain energy density function as

W^ ₁ ,A ₂ , ₃ ,J) = ( ¹ + φ - 3) + (In J) ² ,

«! 2 ₍₃₎ where ½ 2 are the deviatoric principal stretches, μ and <¾ are material constants, J is the Jacobian matrix for the lung deformation, and K is a bulk modulus selected sufficiently high to satisfy near-incompressibility. For the parameters of the Ogden model, we use, e.g., μγ = 0.0329 , and <¾ = 6.82 .

[0044]

We combine all the information (meshes, loads, and boundaries) defined above into a single script file and directly run the FE analysis to simulate the lung deformation using FEBio developed at the University of Utah Musculoskeletal Research Laboratories. FEBio is an open source nonlinear finite element solver specifically designed for biomechanical applications. One important reason that we choose FEBio is that it is specifically designed for biomechanical applications, and offers constitutive models and boundary conditions that are relevant to the modeling of soft tissue deformation. The wide range of boundary interactions supported in FEBio makes the modeling of pleural sliding very straightforward.

[0045]

4DCT Scan Simulation

The deformation of the rest of the CT volume is very complicated if we were to use a prior art continuum mechanics based method, which requires segmenting out all the soft tissues, organs, bones, and other structures, and manually defining different material properties and their corresponding boundary constraints. Instead, we treat the rest of the CT volume as a discrete structure of elements and model the CT volume deformation using a mass-spring simulation system.

[0046]

System Element Definitions

We construct the mass-spring system by sampling control points (masses) in the CT volume excluding the lung region. Then, we connect these points using a constrained Delaunay tetrahedralization procedure to form the links (springs), as a cutting plane. To define mass values for each control points, we first determine their physical densities from the CT image values, i.e., Hounsfield units (HU). The conversion from density to mass is straightforward. The HU versus electron/physical densities curve for human organs and tissues is available for specific CT devices, e.g.,

6.0e ^~ - h(x) + 1.081 otherwise, ^ where p(x) represents the physical density g/cm ³ at position x, and h{x) HU value.

[0047] The guideline to define spring constants for the mass-spring system is to assign large values to relatively dense bones, small values to soft tissues, lung, blood, air, etc. based on their HU values. Considering the HU value for bones ranging from +700 to 3000 (dense bone), we nonlinearly map the HU values [700, 3000] to spring constants [0.8 1] based on the reciprocal of an exponential function, while the spring constants of other materials are mapped to [0.001 0.8]. Other heuristic mappings are also possible and have little impact on the system performance.

[0048]

Dynamics and Forces

We deform the rest of the CT volume based on the lung deformation at each breathing phase. Therefore, our system is an unforced version of a conventional mass-spring-damper system, that is, the external force F _ext = 0. The input is the prescribed displacements of the lung surface vertices, which generates the tensions for the internal springs by Hooke's Law after each lung deformation. We also add damping force to reduce the amplitude of system oscillations. Thus, from

Newton's second law, the total force of the system F _tot is equal to the internal force F _int , and defined as

where F _s is the spring tension defined as F _s = -k , = -R -dy/dt , k is the n x l spring constant vector, y is the displacement vector of the n control points, m is the mass vector corresponding to the control points, and operator denotes the inner product of two vectors. We define the system as a critically damped system, therefore, the damping factor R = 2 m k . Thus, Eqn. 5 can be derived as,

m y + R y + k y = 0 , ^ where the dots above the variables indicate the respective derivatives.

[0049]

We solve this dynamic system using a finite difference method. Between each iteration, we fix the motion of mass points covered by the spine region to simulate the stationary spine movement. After we have the deformation of the control points, we can determine the deformations of all voxels by trilinear interpolation. Other interpolation method, such as, tricubic interpolation, can also be used to obtain smoother and continuous differentiable volumes.

[0050]

Volume Mesh Generation

To generate the lung volume mesh for the FE analysis of the lung deformation, we use 3D tetrahedralization of the surface mesh. We can highlight image regions where user labeling could improve the segmentation results. However, because the tetrahedralization is for general purpose segmentation and not specialized for lung segmentations, additional mesh editing may be required to fix small missing pieces and holes. The number of surface vertices can be reduced to -4000 using quadric edge collapse decimation for simplification. After the surface mesh is given, the lung simulation can run automatically.

[0051]

Fig. 10 shows a portion of a surface mesh. One common problem of the surface mesh generated from segmentation tools is that the mesh is highly irregular, and contains many thin triangles and over-sampled regions. This can lead to badly shaped elements during the 3D Delaunay tetrahedralization process and delay convergence of biomechanical simulations. To optimize the surface mesh, our simulator converts the surface mesh into a solid volume, and then determines a remeshed surface using a marching cube procedure over the volume. This mesh optimization can also fix many other problems with the segmentation results, for instance, self-intersecting faces and non-manifold edges and vertices.

[0052]

Fig. 11 shows that the optimized surface contains well-shaped triangles and uniformly sampled vertices, which constitute a valid basis for 3D Delaunay tetrahedralization.

[0053]

Our biomechanical motion simulation method for lung deformation is very accurate considering a very low spatial resolution of the CT data (0.97 x 0.97 x 2.5 mm), especially in z direction. The method works well in the lower posterior region of the lung where nodal displacement is mostly prominent. This indicates that our lung deformation simulator can provide valuable location-specific tumor motion information to significantly reduce the margins between clinical target volume (CTV) and planning target volume (PTV) to potentially improve survival rates for lung-cancer patients.

[0054]

4DCT Scan Simulation

To generate simulated 4DCT scans, we treat the lung deformation as the equivalent force that drives the deformation of the CT volume because we have already considered all the influencing factors, for instance, the movement of body outline and diaphragm, into the biomechanical simulation of lung deformation. Therefore, or every lung deformation, we can determine its corresponding CT volume deformation through our nonlinear mass-spring-damper system.

[0055] During our simulation, the abdominal organs are squeezed down during rib cage expansion^ while the spines are fixed at the same time. This demonstrates that our simulator is capable of providing realistic deformation and displacement of liver, bones, and other soft tissues besides the accurate respiratory motion of lung.

[0056]

Tumor Embedding

In our simulator, tumors can be embedded into the lung volume mesh through determining barycentric coordinates of their vertices to nearby lung tetrahedra. During X-ray images integration, we can also synthesize the ground truth 3D tumor contours. By putting different types of rumor at different positions of the lung, our simulator is able to provide unbiased ground truth for validating the accuracy of different tumor tracking algorithms. Compared with X-ray videos with metallic markers, our X-ray images remove the side effects of high contrast regions introduced by markers and can provide realistic images for tracking procedures to extract accurate region features. Our simulator also eliminates the need for manually drawing ground truth tumor contours, which can contain large inter- and intra- physician discrepancies.

[0057]

Our simulator can be extended to simulate real-time patient-specific 4DCT scans for particular synchronized pleural pressure and chest and abdomen height using only the single CT scan as input. In this case, the FE analysis of the lung deformation converges very quickly because there is only a small difference of the lung shape between two time steps. The pleural pressure can be estimated from the lung volume readings through the pressure-volume curve, and yields the precise negative pressure for defining the deformation force.

[0058] The synchronized chest and abdomen height yield the patient body outline during respiration, and define the boundary constraints for both lung deformation and 4DCT scan simulation.

[0059]

The invention provides biomechanical motion model based thoracic 4DCT simulation method for patient lung deformation induced by respiration using only one CT scan as input. The method models lung stress-strain behavior using a hyperelastic Ogden model, and treat the rest of the CT volume as discretized mass points connected by springs and dampers.

[0060]

This heterogeneous design leverages the advantages of both continuum mechanics and mass-spring-damper system in the way that the lung deformation is determined with high accuracy, while the deformation of the rest of the CT volume is achieved under practical computational constraints. Extensive analysis and comparisons with the manually labeled dataset demonstrate that our lung deformation results are very accurate.