Abstract

This paper presents a method to reconstruct and to calculate geometric invariants on brain tumors. The geometric invariants considered in the paper are the volume, the area, the discrete Gauss curvature, and the discrete mean curvature. The volume of a tumor is an important aspect that helps doctors to make a medical diagnosis. And as doctors seek a stable calculation, we propose to prove the stability of some invariants. Finally, we study the evolution of brain tumor as a function of time in two or three years depending on patients with MR images every three or six months.

1. Introduction

Cancer is a disease that starts in our cells. Our bodies are made up of millions of cells grouped together to form tissues and organs such as muscles and bones, the lungs, and the liver. Genes inside each cell order it to grow, work, reproduce, and die. Normally, our cells obey these orders, and we remain healthy. But sometimes the instructions get mixed up, causing the cells to form lumps or tumors, or spread through the bloodstream and lymphatic system to other parts of the body. Tumors can be either benign (noncancerous) or malignant (cancerous). Benign tumor cells stay in one place in the body and are not usually life-threatening. Malignant tumor cells are able to invade nearby tissues and spread to other parts of the body. Cancer cells that spread to other parts of the body are called metastases.

Brain cancer starts in the cells of the brain. The brain is a soft mass of nerves (neurons) and supportive tissue (glial cells), surrounded by membranes (meninges) and protected by the skull. The brain has 3 main areas (see Figure 1).

(i)The cerebrum is the largest part of the brain and is made up of the right and left cerebral hemispheres. It allows us to see, feel, think, speak, and move. The right side of our brain controls the left side of our body and vice versa.(ii) The cerebellum is located in the back of the brain and controls balance and coordination.(iii) The brain stem controls our vital bodily functions, like heartbeat, breathing, and reflexes. It connects the brain to the spinal cord. The skull is hard and cannot expand, so as a tumor grows the pressure can damage or destroy delicate brain cells. Brain cancer can involve either the neurons or the glial cells. Most adult cancers start in the glial cells and are called astrocytomas or gliomas.

Neurological specialists perform several diagnostic tests for brain tumors. (http://www.medecinenet.com). These tests include the following.

(i)Electroencephalogram (EEG).(ii)Lumbar puncture (spinal tap).(iii)RN (radionuclide).(iv)Computerized axial tomographer (CT or CAT).(v)MRI.(vi)PetScan.(vii)Biopsy. Since the development of tomography, computers have been used extensively in medical diagnosis. Unlike classical radiology, tomography requires complex mathematical calculations in order to obtain a two-dimensional image. Computers are used for image treatment, visualization, and archiving, and also for 3D reconstruction. In [1], it is proposed a method to reconstruct a 3D model of certain organs from a number of 2D cross-sectional images. This method enables a better understanding of spatial structures, and also open the way to new applications like radiation therapy planing and surgical planing.

In [2], it is proposed a new model to simulate the growth of glioblastomas multiforma (GBM), the most aggressive glial tumors. This simulation has a different medical applications, including an optimized dosimetry in radiotherapy or a better neurosurgical planing in case of tumor resection.

In [3], it is proposed a high-resolution three-dimensional (3D) connectivity, surface construction, and display algorithms that detect, extract, and display the surface of a brain from contiguous magnetic resonance (MR) images. The algorithms identify the external brain surface and create a 3D image, showing the fissures and surface convolutions of the cerebral hemispheres, cerebellum, and brain stem. For the purposes of the 3D reconstruction, it is shown that T1-weighted images give better contrast between the surface of the brain and the cerebral spinal fluid than T2-weighted images. 3D reconstruction of MR data provides a noninvasive procedure for examination of the brain surface and other anatomical features.

In this paper, we reconstruct the tumor from 2D sections, coming from MRI sequences, for patients having a brain tumor. The distance between two parallel sections is 5 mm. The tumor has been segmented semiautomatically by a software provided to us by (INRIA, Cedex, France). and based on the gradient method. After the identification of points on the contours, 50 points were spot on each contour. The reconstruction method is based on these points. Once the tumor is reconstructed, we calculate the geometric invariants (volume, area, and discrete Gauss curvature). The doctors seek a stable calculation that does not depend on the number of points or their distributions on the contours. This is why the distribution and the number of points are changed, and the calculation of geometric invariants is redone to show the stability. The discrete mean curvature is calculated on the edges of the triangulation surface. However, this invariant is not stable according to the distribution and the number of points on the contours.

2. Detection and Segmentation of Brain Tumors

Detection in MR image with brain tumor is an important image processing technique applied in radiology for 3D reconstruction. Indeed, contours are rich indexes for any subsequent interpretation of the image. Contours in image are due to

(i)discontinuities of the reflectance function,(ii)discontinuities of depth (boundaries of the object). The contours are characterized by the discontinuities of the intensity function. Therefore, the principle of contours detection is based on the study of the derivatives of the intensity function in the image. The contours characterize the boundaries of the objects, and generally they are defined as a transition zones between two regions of different characteristics presented simultaneously to within a single digital image.

Definition 1. Let be the intensity function of an original image; the gradient of the image is defined by the vector This gradient is characterized by a module and a direction in the image

The direction of the gradient maximizes the directional derivative. The basic steps of contours detection are thus to calculate the derivatives of the intensity function, then to specify the contour points.

The collaboration with Professor François Cotton was a very important step, especially to locate and to precise the position of the brain tumor in each 2D sequence for all patients. In our applications, the boundary of the tumor was detected by François Cotton, and so the final step before reconstructing the tumor was to segment it by using a software provided to us by INRIA [4]. The idea of this algorithm is to surround the tumor by a circle, to regard the minimal and maximal intensities on the boundary of the tumor, and to fix thresholds for them, then a threshold “” for the gradient is fixed to determine the tumor's boundary. Once these parameters are fixed, we apply an ultra classical procedure of minimization to obtain the boundary of the brain tumor by moving the circle.

Once the gradient is evaluated, the points on the contours which are characterized by local extrema are identified. The idea is to select the pixels by using the threshold “” for the norm of the gradient, that is, all points on contours such that (see Figure 2).

3. Surfaces

A surface is a topological space in which each point has a neighborhood that is homeomorph to the unit disk or to the half unit disk . The boundary of a surface , denoted by , is the set of points of this surface that do not have a neighborhood homeomorph to the unit disk. The interior of is complementary to its boundary. A surface is therefore abstract; the only thing we know is the neighborhoods of each point.

3.1. Triangulations

We will now give some definitions of discrete surfaces, that is, surfaces defined by a finite number of points. Indeed, in computer science we work often on such surfaces, which in generally are approximations of real or ideal surfaces, for example, the rabbit in Figure 3—the famous bunny of stanford—is a mesh surface, created by scanning a real model in a clay (the scanner detects the geometric position of some number of points of the model, and then these points are related three by three to form triangles, using an adequate algorithm).

Definition 2. Let be a set of linearly independent vectors in the Euclidean space , . One calls -simplex of vertices the convex hull of these points. is the dimension of .

Example 1. In , a -simplex is a vertex, a -simplex is an edge, a -simplex is a triangle, and a -simplex is a tetrahedron (see Figure 4).

Definition 3. Let be a set of linearly independent vectors and its convex hull. Then, the convex hull of all subset of is a simplex subset of . One says that is a face of , and one obtains .

Definition 4. A simplicial complex is a finite set of simplices such that
(i)if , then all its faces are in ;(ii) let , then or (the two simplices have a common face) (see Figure 5).
The dimension of a simplicial complex is the dimension of its greatest simplex. We will now establish a link between the topology of a set of points and combinatorial topology, more precisely between the notions of topological space and the simplicial complex.

Definition 5. Let be a simplicial complex in . The union of all the simplices in with the topology of the subsets of is called the polyhedron of .

Definition 6. The triangulation of a topologic space is a simplicial complex such that its polyhedron is homeomorph to . If such simplicial complex exists, one said that is triangulated.

4. Reconstruction Method

The method presented in this paragraph is a combination of the 2D Delaunay triangulation of contours and the maximizing volume method.

4.1. Voronoi Diagram and Delaunay Triangulation

Definition 7. Given a set of sites, such that no four sites lie on a common circle. One defines as the set of points closer to site than to any other site in : is called the Voronoi cell associated to . The union over all the is called the Voronoi diagram of (see Figure 6).

The boundaries of the regions are referred to as Voronoi edges, the joints of three Voronoi edges are called Voronoi vertices. Each Voronoi edge is associated with two adjacent Voronoi cells, a Voronoi vertex is equidistant to three sites (see [5]).

If we draw a line segment between each pair of sites whose Voronoi cells share an edge, we obtain a triangulation of the points in called the Delaunay triangulation (see Figure 7).

A Voronoi vertex represents a Delaunay triangle; more precisely, it is the center of its circumcircle. Each Voronoi edge corresponds to an edge in the Delaunay triangulation despite the fact that they may not even intersect. This geometric difference between Voronoi diagram and Delaunay triangulation becomes important in the reconstruction issue. Some further properties of the Delaunay triangulation are the following.

(i)The boundary of the Delaunay triangulation is the convex hull of its sites.(ii)The Delaunay triangulation is unique.(iii)The number of triangles in the Delaunay triangulation is at most , where is the number of vertices in the triangulation.(iv)A Delaunay triangle does not contain any other site in its circumcircle (empty circle property).

In our application, the object contours are given as a set of straight line segments, forming one simple closed polygon. However, just triangulating the polygons vertices—or calculating the Voronoi diagram of point sites—may result in a triangulation where contour segments are not guaranteed to be edges of the triangulation. Since our goal is to get a 3D polyhedron whose intersection with the given cross-sections yields the original contours, our triangulation has to satisfy the following requirement:

All contour segments have to appear as Delaunay edges in the Delaunay triangulation (contour containment condition).

To obtain a triangulation that satisfies this condition, we simply calculate the Delaunay triangulation of the polygons vertices. Then, we check each contour segment to be contained in the Delaunay triangulation. Segments that do not appear as Delaunay edges are split into two parts by adding a new vertex in the middle. We add the new vertex to the Delaunay triangulation and verify the contour containment condition one more.

It can be shown that such a procedure terminates and yields a Delaunay triangulation that satisfies the containment condition. The contour shape is not changed, since we add vertices onto contour segments (see [6]).

Once all contour segments are in the triangulation, we eliminate the edges which are outside the contour, the vertices added to the contour, and the edges related to these vertices without touching the original segments of the contour (see Figure 8).

4.2. Triangulation of Two Parallel Sections

Let us consider two parallel sections (contours) , with the same number of points on each one. On , we have the sequence of points , and on the points . In other words, if , then we have points on , as a result , and we have points on . Note that the points , on each contour are distributed in the clockwise direction (see Figure 9).

(i)First step. Let be the starting point in our triangulation, the idea of this method is based on finding the nearest point to in the adjacent section. In other words, we calculate the distance(ii)Second step. We choose the point such thatwhere is the Euclidean distance between two points.(iii)Third step. We connect by a segment the point to the nearest point in the adjacent section which is the point . The first segment of the triangulation will be , then on the section we take the neighborhood point of in the clockwise direction, let be this point, by a segment, we connect to the point in the adjacent section to form with the first triangle in the triangulation which has the vertices (see Figure 10).(iv)Fourth step. We take the point , and we determine the nearest point to in the adjacent section without regarding the point . As the sections are parallel, the nearest point will be . By repeating the same procedure in the third step, we find the second triangle which has the vertices . Automatically, the triangle which has the vertices will form between the two triangles: and (see Figures 11, 12).

By repeating these steps times, we obtain the triangulation and the volume delimited between the two sections .

4.3. Changing the Starting Point

Let us initialize the triangulation and change the starting point, it means instead of as a starting point, we take the point and apply steps to get the triangulation and the volume between .

The change of the starting point is done times to reconstruct triangulations (), and volumes () delimited between (see Figure 13).

4.4. Selection of the Best Triangulation

The triangulation which maximizes the volume of the polyhedron , gives the optimal approximation of the surface provided by a pair of closed contour segments (see [7]).

Remark 1. This algorithm is true for convex contours and concave contours (see Figure 14).

5. Surface Area

In our applications, we calculate an area approximation of a smooth surface by the sum of triangles area. To compute the area of a piecewise linear 2-dimensional space, one divides it in a partition of triangles, computes the area of each triangle with the familiar formula—half the product of the basis by the heigh—and then adds all these areas (see Figure 16). The only point to check is that the result is independent of the triangulation. In Figure 15, with 50 points on the contours, the area is equal to , with points, it is equal to and with points, it is equal to (the unit is ); see [8].

6. Volume of a Domain Limited by a Surface

Let be a domain in limited by the surface .

Definition 8. The volume of is given by where is the intersection of with the plane .

Finally, the Green-Ostrogradsky theorem reduces the calculus of the volume to an integral surface:where is the boundary of , and is the unit normal vector at oriented to the exterior of .

7. Discrete Gauss Curvature

We will give the definition of the Gauss curvature at a vertex in a triangulation (see [9]).

Let be a triangulation, is a vertex of , is the set of triangles of that having as a vertex. We denote by the set of all interior vertices of and by the set of all boundary vertices of . The angle at a vertex is the real number defined bywhere is an angle at of the triangle (see Figure 17).

The characteristic Euler of a triangulation is given bywhere is the number of vertices of , is the number of edges, and is the number of faces. A simple calculus gives the following formula of Gauss-Bonnet: This formula will motivate the following definitions.

(i)The discrete Gauss curvature at a vertex is(ii)The total Gauss curvature of is

8. Discrete Mean Curvature

Let be an oriented triangulation, is an interior edge of , and are the unit oriented normal vectors of the two triangles adjacent to (see Figure 18).

(i)We call dihedral angle of at the angle(ii)We define the index convexity of the integer by (iii) The mean curvature of the edge is(iv)The total mean curvature of iswhere is the length of .

9. Patients and Images Modality

The following description precises the different parameters of the images modality, and they are given by the second author.

All patients had partial seizures and a first MRI compatible with brain tumors. Low-grade gliomas were confirmed with neurosurgical biopsy, and patients were then followed with serial MRI every 6 months. MR imaging was performed on Philips Intera 1.5 Tesla MRI system (Philips Medical Systems, Erlangen, The Netherlands) with a standard head coil. Localizing sagital T1-weighted images of the brain was obtained initially, followed by axial FLAIR and T2-weighted images. Thereafter, diffusion- and perfusion-weighted MR images were obtained (diffusion-weighted image comprised an echo-planar spin-echo sequence, with the following parameters, TR = 4247, TE = 95, EPI factor = 77, field of view (FOV) = 230 mm, slice thickness = 5 mm, slice gap = 1 mm, number of excitations = 1, matrix = 77 256, number of slices = 22, acquisition time = 30 seconds). The images in diffusion-weighted image (DWI) were acquired for values of b (diffusion gradient factor) equal to 0 and 1000 . An isotropic image was constructed in real time, pixel by pixel, as an average of the signal intensities of three orthogonal directions. Apparent diffusion coefficient (ADC) maps were then calculated from the diffusion-weighted images. Dynamic susceptibility contrast enhanced MRI (DSC-MRI) using the PRESTO sequence with the following parameters was employed: TR = 17, TE = 8, flip angle = 7, EPI factor = 17, FOV) = 220 mm, slice thickness = 3.5 mm, slice gap = 0 mm, number of excitations = 1, matrix = 64 64, number of slices = 60 series of 22 slices, acquisition time = 77 seconds. The PRESTO technique is highly sensitive to T changes due to its very long TE. Image acquisition began simultaneously with contrast agent injection. For all patients, a standard dose (0.2 mL/kg = 0.1 mmol/kg) of gadobenate dimeglumine (multihance; Bracco Imaging, Milan, Italy) was injected using an automated power injector, at a flow rate 6 mL/sec followed by a 40 mL saline flush at the same rate. Contrast agent administrations were in all cases performed using a 20-gauge intravenous catheter. Finally, postcontrast 3D T1-weighted images, with an isotropic voxel of 0.9 mm, were obtained.

10. Results

The collaboration with professor François Cotton (IRM, Hospices Civils de Lyon, Lyon, France) consists of the three following steps.

(i)First step. Series of 3 typical biopsy-proved low-grade gliomas followed with serial MRI every 6 months were first selected for the purpose of the study. FLAIR sequence was used for the segmentation because it has high sensitivity for the tumor detection, which typically appears on hypersignal compared to the normal brain parenchyma. DICOM images including FLAIR sequence were then burned on a CD.(ii)Second step. The validation of the segmentation was done firstly by comparing the medical data with the slope of the tumor growth obtained by linear measurement—secondly according to neuroexperience of the tumor growth. In medical experience and literature, the mean tumor diameters of untreated low-grade gliomas (linear measurements) showed a linear and constant growth before anaplastic transformation with an average slope of 4 mm per year (2 to 8 mm/year) (see [1012]).(iii)Third step. The exchange between mathematicians and specialists in medical imaging was based on knowledge of MRI sequences, tumor growth, macroscopic and microscopic structure, microvascularity and angiogenesis, density and contour of brain tumors, and mathematical model that could be used to analyze these parameters.

Our goal is to provide some tools to help the doctor to make his diagnosis, then to show that some geometric invariants (like the volume, the area, and the Gauss curvature) are relatively stable by adding points at the segmentation. Others, like the mean curvature, are sensitive to the number of used points. Finally, we would like to sensitize the community of the medical imaging in the importance of calculating geometric invariants in order to verify that the segmentation and the 3D reconstruction of the tumor are good from this point of view.

We present the results established on some patients in Figures 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, and 33.

10.1. Evolution Curves

We present the evolution curves of the brain tumor as a function of time in Figures 34, 35, and 36.

10.2. Interpretation and Conclusion

In this section, we give an interpretation of the results previously shown. Patients 1 and 3 have curves tumor growth which agree with the expectations of doctors and the standard models of growth tumor. The curve of the second patient is characterized by the volume evolution of the tumor after a surgical operation. After a phase of growth, the two last dates show that the tumor has been partially removed and that a remission phase appears.

Results are consistent with the clinical history and natural course of low-grade gliomas. 3D segmentation is more sensitive and consistent than linear measurement. A fully automatic method of volume and contour segmentation will be of great interest for the clinician and patient management. It is now recognizing that all low-grade gliomas (except type 1 such as pilocytic astrocytoma) will present an anaplastic transformation which highly correlates with tumor growth.

After this study, we remark that some geometric invariants are stable if we increase the number of points, or if we change the distribution of points. The Gauss curvature is always equal to , which is normal, because the triangulated surface is a deformable sphere without holes. However, the mean curvature depends heavily on the number of points. Once we increase the number of points, the surface becomes increasingly disturbed.

Our study helps to precise the importance of the geometric invariants during the reconstruction of brain tumors in medical imaging. It should be noted that another algorithm was tested to reconstruct the tumor but topological problems were found, especially for the Gauss curvature for some 3D objects. It seems essential to verify the stability and the values of geometric invariants in order to give valid estimations to the doctors.

Acknowledgments

This work is supported by the ANR JCLAMA no. 05-BLAN-0140-01. The authors would like to thank Mustapha Khamiss, INRIA, and the Hospital of Lyon for providing the materials of this paper.