Evolution of Brain Tumor and Stability of Geometric Invariants

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.

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 crosssectional 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 T2weighted 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.

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 I(x, y) 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 m and a direction φ in the image m = (I(x + 1, y) − I(x, y)) 2 + (I(x, y + 1) − I(x, y)) 2 , 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 "s" 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 "s" for the norm of the gradient, that is, all points on contours such that m > s (see Figure 2).
International Journal of Telemedicine and Applications Figure 2: Segmentation of the brain tumor and the contour points. Segmentation validated by François Cotton (professor of radiology).

Surfaces
A surface is a topological space in which each point has a neighborhood that is homeomorph to the unit disk {(x, y) ∈ The boundary of a surface S, denoted by ∂S, is the set of points of this surface that do not have a neighborhood homeomorph to the unit disk. The interior of S is complementary to its boundary. A surface is therefore abstract; the only thing we know is the neighborhoods of each point.

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 e 0 , . . . , e n be a set of n + 1 linearly independent vectors in the Euclidean space R m , m ≥ n ≥ 0. One calls n-simplex of n vertices e 0 , . . . , e n the convex hull σ of these points. n is the dimension of σ.
Example 1. In R 3 , a 0-simplex is a vertex, a 1-simplex is an edge, a 2-simplex is a triangle, and a 3-simplex is a tetrahedron (see Figure 4).

Definition 3.
Let S be a set of linearly independent vectors and σ its convex hull. Then, the convex hull τ of all subset T of S 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
then all its faces are in K; (ii) let σ i , σ j ∈ K, then σ i ∩ σ j = ∅ or σ i ∩ σ j ∈ K (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 K be a simplicial complex in R m . The union |K | of all the simplices in K with the topology of the subsets of R m is called the polyhedron of K. Definition 6. The triangulation of a topologic space X is a simplicial complex K such that its polyhedron |K | is homeomorph to X. If such simplicial complex exists, one said that X is triangulated.

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

Voronoi Diagram and Delaunay Triangulation
Definition 7. Given a set S of N sites, S = {p i ∈ R 2 /i = (1, . . . , N)} such that no four sites lie on a common circle. One defines V (i) as the set of points closer to site p i than to any other site in S: International Journal of Telemedicine and Applications V (i) is called the Voronoi cell associated to p i . The union over all the V (i) is called the Voronoi diagram of S (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 S 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  (ii) The Delaunay triangulation is unique.
(iii) The number of triangles in the Delaunay triangulation is at most 2N − 5, where N 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).

Triangulation of Two Parallel Sections.
Let us consider two parallel sections (contours) C 1 , C 2 , with the same number of points on each one. On C 1 , we have the sequence of points A 0 , A 1 , . . . , A m , and on C 2 the points B m+1 , B m+2 , . . . , B n . In other words, if m = 49, then we have 50 points on C 1 , as a result n = 99, and we have 50 points on C 2 . Note that the points A 0 , A 1 , . . . , A m , B m+1 , B m+2 , . . . , B n on each contour are distributed in the clockwise direction (see Figure 9).
(i) First step. Let A 0 be the starting point in our triangulation, the idea of this method is based on finding the nearest point to A 0 in the adjacent section.
In other words, we calculate the distance where d is the Euclidean distance between two points.
(iii) Third step. We connect by a segment the point A 0 to the nearest point in the adjacent section C 2 which is the point B k . The first segment of the triangulation will be [A 0 B k ], then on the section C 1 we take the neighborhood point of A 0 in the clockwise direction, let A 1 be this point, by a segment, we connect A 1 to the point B k in the adjacent section to form with A 0 the first triangle in the triangulation which has the vertices A 0 A 1 B k (see Figure 10).
(iv) Fourth step. We take the point A 1 , and we determine the nearest point to A 1 in the adjacent section without regarding the point B k . As the sections are parallel, the nearest point will be B k+1 . By repeating the same procedure in the third step, we find the second triangle which has the vertices A 1 A 2 B k+1 . Automatically, the triangle which has the vertices B k B k+1 A 1 will form between the two triangles: Figures 11,12).
By repeating these steps m times, we obtain the triangulation T 0 and the volume V 0 delimited between the two sections C 1 , C 2 .

Changing the Starting Point.
Let us initialize the triangulation and change the starting point, it means instead of A 0 as a starting point, we take the point A 1 and apply steps 1, 2, 3, 4, . . . to get the triangulation T 1 and the volume V 1 between C 1 , C 2 .

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

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 149.084, with 100 points, it is equal to

Volume of a Domain Limited by a Surface
Let D be a domain in R 3 limited by the surface S.
where D z is the intersection of D with the plane z = const.
Finally, the Green-Ostrogradsky theorem reduces the calculus of the volume to an integral surface: where ∂D is the boundary of D, and n is the unit normal vector at dS oriented to the exterior of D.

Discrete Gauss Curvature
We will give the definition of the Gauss curvature at a vertex in a triangulation (see [9]). Let T be a triangulation, p is a vertex of T, T T (p) is the set of triangles of T that having p as a vertex. We denote by 30 where α σ (p) is an angle at p of the triangle σ (see Figure 17). The characteristic Euler of a triangulation T is given by where N S is the number of vertices of T, N a is the number of edges, and N f 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 p ∈ S • T is (ii) The total Gauss curvature of T is

Discrete Mean Curvature
Let T be an oriented triangulation, a is an interior edge of T, n 1 and n 2 are the unit oriented normal vectors of the two triangles adjacent to a (see Figure 18).
where l(a) is the length of a.

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    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.

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 [10][11][12]).
(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.

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

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 4π, 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.