Computer Graphics & Geometry
S. Czanner, R.Durikovic and H. Inoue
Software Department, The University of Aizu
AizuWakamatsu City, Japan
email:
The growth of the brain of human embryo is changing through a long time in the body of mother. So it is very difficult to observe and to understand that process. Therefore, the embryologists found the realistic human organ models and animations necessary for their studies. But to create the realistic human embryo brain models and to do the animations it requires an appropriate methodology. The aim of this paper is to present a developing methodology based on the functional representation and the convolution of the surfaces [4,5,6]. We employed this to create a growth simulation of the human embryo brain. The idea of this technique is the following. As a first step the 2D central skeleton is created from an artistic drawing and then the 3D skeleton is modeled by adding the thickness. As the next step, the skeletons representing keyframe models are used to create an animation. At the end the gap between the keyframe models is filled by the suitable interpolation techniques and finally the animation is composed.
Keywords: human embryo, brain, convolution surfaces, Frep, geometric modeling, blending, skeleton, morphing, animation.
To understand the development of the human organs, the realistic models are necessary. One of the significant organs during a human embryo development is the human brain. The shape of the brain is organic and has many folds that are hard to model with conventional techniques; it is therefore necessary to create a new methodology.
In our research we have modified the techniques developed by Durikovic, Kaneda and Yamashita [1]. There was developed a system in which multiple organs and the environment are separate processes, information from each of which is transmitted using communication modules. The methodology, which is based on algebraic Lsystems with cycles, was demonstrated on a growth model of human embryo stomach. The proposed method is not well designed for use in the growth simulation of other organs, therefore we decided to use for our growth simulation theory of implicit surfaces.
There is other semiautomatic reconstruction method based on implicit isosurfaces generated by skeleton [11] that can be used on noisy scattered points of medical organs. The basis of the method is minimization of the energy that represents the distance between the data and the surface, during a series of skeleton refinements. In our approach, the skeleton consist several triangulated tubular shapes. To generate the models for growth simulation we used a convolution over the triangles.
In the process of the drawing of special objects, computer graphics algorithms must perform a variety of geometrical calculations. The mathematical representation and manipulation of shapes is, therefore, crucial to the geometric modelling process. Various forms are used as geometrical primitives ranging from parametric surfaces through implicit surfaces to solids.
A parametric surface is given by a special position function p(u,v) = [x_{1}(u,v), x_{2}(u,v), x_{3}(u,v)]. In practice, the functions are splines defined by pieces of polynomials, or ratios of polynomials.
An implicit surface defined by function f is the set of points {p � R^{3}:f(p) = 0}, p = p[x_{1}, x_{2},x_{3}] The function f is called an implicit function [4,5]. The implicit surfaces are often used on convolution to round a solid model.
Let us have the following modelling equation in an implicit form:
N S i = 1 
F_{i}(p)T = 0. (1) 
The convolution surface defined by Eq. 1 is called an isosurface, where F_{i} are the source potentials and T is the isopotential threshold value.
To create a model of human embryo brain we used a convolution with piecewise planar skeleton. The skeleton, a standard geometric modeling representation, has become a popular construction for implicit design. It consists of an arbitrary number of elements and each element generates an associated implicit primitive (point, line, triangle, ).
Let have function f be representing the geometry of an implicit primitive. Let h be a potential function that describes the field generated by a single point of that primitive.
The total amount of field at point r, generated by the whole primitive is
F(r) =  � � 
R^{3} 
f(p)h(rp)dp. 
which is a convolution of two functions f and h.
The convolution integral can be conveniently written as an integral of the potential function h(r) over the volume of the primitive:
F(r) =  � � 
V 
h(rv)dv. 
McCormack and Sherstyuk [2,3] define a new kernel that allows analytical calculation of field equation simplified as:
F(r) =  � � 
V 
dv
(1+s^{2}r^{2}(v))^{2} 
where r = rp is the distance from the point p and coefficient s controls the width of the kernel, see Figure .
Figure 1: The convolution kernel.
The shape of a convolution surface can be varied in several ways: by changing the skeleton, the changing the convolution kernel, and by spatial deformation.
Let us consider closed subsets of ndimensional Euclidian space E^{n} with the definition:
f(x_{1}, x_{2},..., x_{n}) � 0, (2)
where f is a real continuous function defined on E^{n}. The above inequality is called a function representation (Frep) of a geometric object and function f is called the defining function. In threedimensional case the boundary of such a geometric object is called implicit surface. The major requirement on the function is to have at least C^{0} continuity. In our model the defining function is the function F(r) defined over a 3D skeleton, which has the C^{1} continuity.
The set of points X_{i}(x_{1}, x_{2},..., x_{n}), i = 0,...,N in E^{n} associated with Eq. 2 can be classified as follow:
f(Xi) >0 if X_{i} is inside the object,
f(Xi) =0 if X_{i} is on the boundary of the object,
f(Xi) <0 if X_{i} is outside the object.
Let binary operations on geometric objects represented by function be defined by
F(f_{1}(X),f_{2}(X)) � 0, (3)
where F is a continuous real function of two variables [5]. Such operations are closed on the set of function representations. The resulting object will have the descriptive function as follows:
For object union
f_{3} = f_{1}f_{2} � 
1
1+a 
(f_{1}+f_{2}+  � 
f_{1}^{2}+f_{2}^{2}2af_{1}f_{2} 
). 
For object intersection
f_{3} = f_{1}&f_{2} � 
1
1+a 
(f_{1}+f_{2}  � 
f_{1}^{2}+f_{2}^{2}2af_{1}f_{2} 
). 


There is an example of using the union operation on Figure . We connected three line segments with following parameters: horizontal line, T = 0.7, S = 0.5, vertical left line: T = 0.7, S = 0.7, vertical right line: T = 0.7, S = 0.5. The connection of the line segments is not smooth enough and there is not a possibility to correct or control the shape of connection, therefore we decided to use blending union operation.
a)  b) 
Figure 2: An example of using the settheoretic operation union. There are three line segments connected by operation union. a) In the middle of each line segment is rendered its skeleton. b) The line segments were rendered as solid objects.
a)  c) 
b)  d) 
Figure 3: Operation Blending union with several blend parameters. We connected three line segments with the same parameters as was showed on Figure 2. For blending we used the same parameters a_{0}, a_{1} and a_{2} a) a_{0} = a_{1} = a_{2} = 0.01 b) a_{i} = 0.05 c) a_{i} = 0.25 and d) a_{i} = 0.5.
The mathematical description of the blending union operation is as follow:
_______ �f_{1}^{2}+f_{2}^{2} 
+ 
a_{0}

where f_{1} and f_{2} are defining functions of the objects which are blended. The absolute value a_{0} defines the total displacement of the blending surface from two initial surfaces. The values a_{1} > 0 and a_{2} > 0 are proportional to the distance between blending surface and the original surface defined by f_{1} and f_{2} respectively, see Figure 3.
The previous chapter discussed the theory of Frep and convolution surfaces. As a next, we will show a method for modeling the organics shapes by Frep.
First step in the model creation process is to obtain the size measurements of brain stages from embryological atlas. Embryological atlas contains handdrawing pictures ordered by age. In the presented growth simulation we used the models from day 28  56 days old brain. There are several possibilities, how to obtain the brain measurements. In our case the brain on the pictures has been measured by ruler. Then upper and lower lines were created and the correspondent points from upper and lower lines were connected. At the end triangular patches were drawn over the pictures and photographs, see Figure .
Figure 4: The conversion of drawing object to central skeleton.
The result of the measurements and the conversion made in previous section is 2D model containing the triangular patches. We call it the central skeleton. In order to be more detailed in the description of some brain models and also because of unwanted blending problems the central skeleton was divided to several parts. An example is shown on the Figure . Namely the part I described rhombencephalon, part II mesencephalon and part III prosencephalon.
Figure 5: Dividing the central skeleton to 3 parts. The line in the middle of the central skeleton is called central line.
The next important step for generating a spatial brain model is to calculate a central line. Central line is the line through the center of each part of central skeleton. It is a base for defining the thickness of the model. The control points of central line were calculated as the middles of corresponding points from vertical segments.
By adding the thickness to 2D model we obtained the 3D skeleton. Then we slightly scaled a multiple number of copies of central skeleton a shifted them to the left and right sides from the central skeleton, see Figure .
Figure 6: Adding the thickness by scaling and shifting the central skeleton.
Each of side skeletons is scaled to fit the ellipses whose center is on the central line. Ellipse has two radiuses a and b. Radius a is a distance to the central line from the border of the central skeleton. Radius b follows the equation: b = aa , where a is a ratio between a and b which is used to create enough thickness. As next step the side skeletons are translated by t, which has been obtained from t = ccosq. The description of the variables c and q is shown on Figure . Finally, side skeletons are connected with a central skeleton or with other side skeletons by a triangular mesh. After erasing all interior triangular patches, which was done by triangular patches indexing we obtain multiple tubular shapes forming together the entire skeleton of the brain.
Figure 7: A 3D skeleton for 36 days old human embryo brain.
The 3D skeleton model, shown on the Figure 7 is a 3D object created from union of three tubular shapes used for quick preview.
The last step in the model creation process is to create a smooth convolution surface defined over the skeleton. Final convolution surface is calculated by convolution operation between the kernel and the skeleton triangles for each tubular skeleton part. Geometric primitive corresponding to respective skeleton parts are given together by set theoretic or by arithmetic operations between defining functions.
In order to create brain model with convolution surfaces, we use HyperFun [6,7] as modelling software and POVRay [8] as rendering software.
To generate a convolution surface over the triangles we used the HyperFun command hfConvTriangle. For the entire brain model, which is shown on Figure right, the convolution parameter S has been set to 0.5 and T to 0.6. We have set the brain thickness parameters to 1.0 at parts I and II ant to 1.2 at part III. To create the entire model all three parts are connected with blendunion operation using the HyperFun command hfBlendUni. In our example the blending parameters a_{1},a_{2} and a_{3} are equal to 0.2.
Figure 8: The human embryo brain, age: 36 days. Left: 3D skeleton, right: entire brain model, which was defined by function representation.
Figure 9: Composition of entire brain model (age: 28 days) with RTG image of human embryo.
In order to make animation, the keyframe models and morphing interpolation between them are important.
The embryological atlas consists of 8 artistic drawings of the human embryo brain development. We have used first 7 images as keyframe models for the animation process. The keyframe models were named Stage1,..., Stage7. The base measurements and statistical information is collected in Table .
Stages  Age  Size (mm)  Number of  
[days]  [mm]  skeleton triangles  parts  
Stage1  28  3.5  296  1 
Stage2  32  5  476  1 
Stage3  36  9  588  3 
Stage4  42  11  548  3 
Stage5  49  15  724  4 
Stage6  56  27  704  4 
Stage7  72  56  710  4 
Table 1: All stages of brain model data.
The shape transformations can be applied to the various types of graphical objects, such as 2D drawings, images, surfaces and volumes. Given a certain specification, there are various choices in implementing the transformation, not only in selecting the type of modeling, but also in deciding how the object data will be transformed. There does not exist a general morphing technique we could use in our case because the available techniques use strictly the shape and topology information but they neglect the known growth processes, movements and knowledge of embryologists. A work described in [10] proposed a skeleton feature vectors but they still use a blending technique to hide the topology errors that occur during their 3D morphing step.
We decided to propose the featuredbased 3D morphing technique with a simple user interaction to control the complicated growth process. Our technique uses both the global deformations to roughly match the global movements of the brain during the growth and the local morphing technique to correct the shape details.
To generate the models between keyframe brain stages we used various types of interpolation. Between Stage1 and Stage2, Stage2, Stage4 and Stage4, Stage6 and Stage6, Stage7 the tricubic interpolation based on CatmullRom interpolating curves [9] was used. We used also this interpolation technique to produce a new temporary skeleton points inbetween frame skeletons.
The equation for CatmullRom splines for onedimensional case, can be expressed by the following matrix formula:
C(u) = UMP^{T}, (4)
where U = [u^{3}, u^{2}, u, 1], P = [p_{i1}, p_{i}, p_{i+1}, p_{i+2}] , and
M=  � � � � � � 
0.5  1.5  1.5  0.5  � � � � � � 
1.0  2.5  2.0  0.5  
0.5  0  0.5  0  
0  1  0  0 
C(u) is the interpolated value, p_{i1}, p_{i}, p_{i+1}, p_{i+2} are four consecutive data points and u � [0,1] is a parameter that defines the fractional position between p_{i} and p_{i+1}.
CatmullRom interpolation is used to interpolate position between skeleton points, values of convolution filter parameters S and T and thickness parameters a, q and last blending parameters a_{0}, a_{1} and a_{2}. This technique generates good inbetween models. Some of them are shown on Figure .
Figure 10: A growing sequence of human embryo brain. The age of embryo is 56 to 72 days.
The visual result of modeling the human embryo brain using function representation is an animation.
The number of inbetween keyframe models depended on the age of keyframe models. We calculated 10 frames for 1 day of the growth. The speed of final animation corresponds to 12 frames per development day. A single frame from the final animation is on the Figure 10.
Figure 11: Composition of entire brain model (age: 28 days) with RTG image of human embryo.
We succeeded to model human embryo brain using convolution surfaces and to create growth animation between 28 days old human embryo brain and 56 days old human embryo brain. The same approaches can be used for other human organs. We have used and improved an idea to model the organic shapes by function representation and convolution surfaces. Some of the advantages of this approach are to avoid the topology artifacts during the animation of inbetween shapes and easier key framing by a simple animation of the skeleton of convolution surface.
In the future work we implement branches and details for brain models and create grownup brain models and growth animation between older human embryos.
Authors wish to thank Mineo Yasuda from Medical School of Hiroshima University for sharing the knowledge as an embryologist. The images were rendered by the POVRay raytracing program programmed by POVRay Team. This research was sponsored by grants from the Fukushima Prefectural Foundation in Japan for the Advancement of Science and Education.
Computer Graphics & Geometry