Zilong Chen

Dazhou Zhu

Xiangrong Ren

Hua Cong
email: dconghua@xaas.ac.cn

Cheng Wang
email: wangc@nercita.org.cn

Chunjiang Zhao

Analyzing Thermal Infrared Image Characteristics of Maize Seedling Stage

Keywords: thermal infrared imaging, temperature, corn, drought resistance, image processing

The temperature of crop leaf is closely related to leaf respiration, transpiration and its drought resistance. The rapid measurement of crop temperature is of great significance for breeding and agricultural cultivation management. The traditional infrared thermometer can only obtain the single point temperature. In this study, a high sensitive thermal infrared imager was used for recording the leaf temperature distribution of maize seedling stage. The temperature distribution of 5 different drought resistance varieties of maize in both normal water and drought stress treatment were analyzed. It turned out that there were obvious differences in temperature for different leaf position. The average temperature in the middle part of the leaf is 0.94℃ higher than the leaf apex and 0.81℃ higher than the leaf base. The results also showed that the leaf temperature is higher than that of the stem. And the average temperature of the whole blade is 24.07℃, which is higher than the background soil temperature (19.53℃). For the same maize variety, the leaf temperature of drought stress treatment is higher than that of normal water treatment. It was also found that the temperature of maize with good drought resistance is higher than that of the bad one, furthermore, the better the higher. The results above showed that the characteristics of temperature distribution in crops could be effectively obtained by using thermal infrared images, which also provided potential for the rapid identification of crop drought resistance.

Introduction

Global warming has given rise to increasingly pressing drought. In China, the arid and semi-arid regions account for 52.5% of the land area, seriously affecting the development of agriculture [START_REF] Guohu | Effects on Traits of Corn Hybrids under Drought at the Different Growing Stages[END_REF]. From 1950 to 2007, the drought-led average annual losses of grains in China reached 15.8 billion kg, making up for over 60% of all natural disasters [START_REF] Yuzhong | The engineering technique for agricultural drought and their application[END_REF]. The growth of maize is threatened by the severe drought. For transpiration, leaf blades of plants control the inflow and effusion of moisture by controlling the size of stomata. Changes in transpiration intensity determine the heat loss degree on the blade surface, thus leading to the changes of blade temperature [START_REF] Lourtie | New glyphosate screening technique by infrared thermometry[END_REF].

Under normal circumstances, the blade temperature of plants maintains a relatively stable state through dehydration from transpiration. Under external stress like drought, the variation of blade temperature can be used to detect the stress conditions of plants [START_REF] Laury | Imaging techniques in early detection of plants stress[END_REF].

The traditional temperature measurements are mainly carried out by the infrared thermometer. But this approach can only realize spot measurement rather than the assay on temperature distribution of a specific area. The increasingly sensitive infrared imaging system provides feasibility for the studies on blade stomata changes through the observation of surface temperature variation of blades [START_REF] Leinonen | Combining thermal and visible imagery for estimating canopy temperature and identifying plant stress[END_REF]. The principle of the thermal infrared imaging technology is as follows: different absorption and reflection of various substances exposed under rays with different wavelengths can be reflected on temperature variation. This property can be used in the detection of temperature variation of maize blade. Compared with conventional ways, it presented advantages of speedy, efficient and convenient operation. Stomatal conductance from the thermal imaging measurement update showed satisfying correlation with that from conventional porometer [START_REF] Jones | The use of thermography for quantitative studies of spatial and temporal variation of stomatal conductance over leaf surfaces[END_REF]. By changing the intensity of the irrigation, times of crop temperature variation can be counted. The infrared thermal imaging can be used to assay the stomatal conductance of crops [START_REF] Bajons | Determination of stomatal conductance by means of infrared thermography[END_REF].

In China, scholars like Liu found that when under the threats of drought, temperature variation of blades can effectively reflect the drought resistance of maize in the seedling stage. The blade temperature difference can be used as the preliminary indicator of the screening of drought resistance [START_REF] Liu Ya | Identification of Maize Drought-Tolerance at Seedling Stage Based on Leaf Temperature Using Infrared Thermography[END_REF][START_REF] Xudong | Canopy-air Temperature Difference of Rice in Key Water Requirement Periods[END_REF]; canopy temperature showed great correlation with leaf water potential of rice, setting percentage and biomass.

Higher temperature leads to smaller indicator. Jones conducted remote sensing study on the infrared thermal imaging of crop canopy, so as to quantitatively analyze the reactions of crops on water stress [START_REF] Jones | Thermal infrared imaging of crop canopies for the remote diagnosis and quantification of plant responses to water stress in the field[END_REF]. Sirault applied infrared thermal imaging to quantitatively analyze the permeability of wheat and barley under the salt stress [START_REF] Xavier | A new screening method for osmotic component of salinity tolerance in cereals using infrared thermography[END_REF].

Qiu introduced a reference dry blade (without transpiration), and applied the plant transpiration conversion coefficient to assay the environmental stress around watermelon, tomato and lettuce and other plants, including conditions of hydropenia and high temperature [START_REF] Guoyu | An infrared-based coefficient to screen plant environmental stress: concept, test and applications[END_REF]. Takai discovered that the canopy temperature difference between sunny days and cloudy days can be used to estimate the subtle changes in stomatal conductance of rice [START_REF] Toshiyuki | Canopy temperature on clear and cloudy days can be used to estimate varietal differences in stomatal conductance in rice[END_REF].

Studies

Experiments and Methods

Experimental Material

Five maize species suitable for cultivation in the north were selected, namely, Kuancheng 60(KC 60), Jingdan 68(JD 68), Jingke 528(JK 528), Dantiannuo

No.3(DTN 3) and Yinnuo No.1(YN 1). With the query of species description from breeding units and breeding experts, the features of species were shown in Table 1.

Experimental Design and Plantation Management

The experiment was conducted in the solar greenhouse of Beijing Academy of Agriculture and Forestry. Flower pots with the aperture diameter of 25cm and depth of 15cm were selected, which were filled with 5kg of soils from fields. Quality corn seeds were selected and sowed in one cavity per pot with two seeds in one cavity.

Each species was sowed for 10 pots so it was 50 pots. 2.0g of compound fertilizer carbonate diamine was applied. To the seedling stage, the final singling of seedlings with consistent growth length was made as 1 plant per pot. After sowing, timely ventilation, illumination, watering and fertilization were conducted according to the greenhouse maize management.

A control group was set in the experiment. In the early seedling stage (15days after sowing, with the twin leaves), drought stress was conducted. 5 pots of each species of maize seedlings stopped watering, taken as the experiment group. The other 5 pots were watered regularly as the control group.

Infrared Image Acquisition

The experimental image acquisition device uses a high-resolution portable infrared thermal imaging instrument (VarioCam). The detector type: UFPA; Spectral range: 7.5 to 14µm; thermal sensitivity: at 30℃, the thermal sensitivity<0.065℃;

Frame frequency: 50/60Hz; infrared image: 384x288 pixels.

20d after sowing, we started the image acquisition. Acquisition time: 8:00am and 14:30pm in sunny days were selected until maize plants died under drought stress.

When photographing, the thermal infrared imager was placed on the top of maize plants, keeping the lens perpendicular to the ground (about 0.5 to 0.8m from the top).

The focal length and camera lens position was adjusted so that the canopy image can be clear and recognized. Thermal infrared images of all plants were taken, coded and preserved.

In the experiment, 5 infrared images were collected on October 17, 19, 23, 25

and 31, 2012. The weather during the experiment was sunny with the temperature varying from 28 to 33℃.

Infrared Image Processing

The original image acquired by the infrared camera was .IRB file, and different colors in the image represented different pixel values. Pixels of different colors show different temperature values. At 8: 00 a.m., as the plant leaves went through no illumination at low temperatures overnight, the plant temperature was almost consistent, no showing distinct differences. Hence, the infrared image taken in the afternoon was selected for analysis.

The IRBIS image processing software was used. First, the target image was opened on IRBIS 3. Color distribution of the image was adjusted through the scale calibration plate, so that the maize in the image was obviously separated from objects in the background to facilitate the processing. According to different shapes of blades, oval, round, or self-crossed circle was selected under the measurement label to circle the specific area of the leaves or stems. The maximum, minimum and average temperature in the circled region was read so as to obtain the temperature information of different blades at different leaf positions. The extracted temperature information was then recorded in the table.

Results and Discussion

Temperature Distribution Characteristics of Maize at the Seedling Stage

Seedling stage of maize refers to the stage from germinating to elongation. The seedling stage of maize is a crucial stage for the differentiation and growth of roots, stems, leaves and other vegetative organs and the differentiation of tassels. In this stage, less moisture was needed. Yet, water stress would lead to poor seedling development, influencing the later nutrition absorption and reproductive growth, and eventually reducing the production [START_REF] Zhuang | Discuss on Evaluating Method to Drought-resistance of Maize in Seedling Stage [J][END_REF]. Thus, it can be seen that some characters of the maize at the seedling stage were closely related with the drought resistance of maize. In this study, the temperature distribution characteristics of maize canopy were studied and the relationship between thermal infrared image characteristics at the seedling stage and drought resistance of maize was discussed.

Infrared images acquired at 2:30 p.m. on October 31, 2012 were selected. Fig. 1 shows the thermal infrared image of a maize plant at the seedling stage. Three blades of the leaf canopy were analyzed with the measurement function of the software. In As shown in Fig. 1, the average temperature of three blades at different leaf positions was recorded in Table 1. It should be noted while circling that the circled area should be smaller than leaves so as to ensure that the circled areas were a part of the leaves. The average blade temperature was measured. Each small area was formed by several pixels. The size of circled area will influence the second digit of the temperature while the accuracy has met the requirements. As can be seen from Table 2, for the temperature of each maize leaf, leaf center > leaf base>leaf opex. The average temperature of leaf center was higher than that of leaf opex by 0.94℃, while that of leaf center was higher than that of leaf base by 0.81℃. The plant temperature was much higher than soils in the pot (19.53℃). The stem temperature (23.09℃) was lower than leaves.

Transpiration and respiration affect the airflow and water flow in the leaves of plants. For moisture, with temperature variation, plants will carry out according self-regulation to offset the impact from water stress, and to adapt to the environment change so as to maintain the development and growth. The changes were reflected on the temperature variation of plant leaves. In this experiment, it was found that the temperature at different leaf positions in the maize canopy at the seedling stage showed the following regularities: leaf opex< leaf base< leaf center (as shown in Table 2). In the normal growth, the leaf shape of maize was arched, and the middle of the leaf was facing the sun, receiving most of the sun irradiation and high temperature; two ends of leaves received less sun irradiation and temperature. With the control over the open and closure of stomata, plants can control the moisture diffusion, thus enlarging the biomass of plants. The mechanism stems from the plant adaptation to drought, which is conducive to maintaining the normal growth and metabolism of plants.

The maize temperature difference between the watering group and the drought stress group

Infrared images acquired at 2:30 p.m. on October 31, 2012 were selected, when the leaves under the drought stress showed wilting. In image processing, the self-crossed way was used to circle the canopy leaves of various species in two groups. The average temperature of leaves was recorded in the following table. To better demonstrate the temperature differences, the table was showed in the following curve graph: It can be seen from the above figures that at the seedling stage of the five maize varieties, the average temperature of canopy leaf under drought stress (25.63℃) was higher than that in the watering group (25.06℃).

Under drought stress, the composition and structure of the plant plasma membrane has undergone significant changes. The cell membrane permeability was damaged so that cells cannot be divided and grow. In addition, stomata on the surface were closed, which influenced photosynthesis, the normal growth and development of plants and eventually lower yields [START_REF] Shengmou | Harm by drought and mechanism of drought resistance of plant[END_REF].

In the experiment, two groups of water treatment were set for the simulation of the arid geographical environment in the north of China. 10 pots of each species were planted. Five pots were not watered under drought stress and the rest were watered.

After that, infrared images were acquired, coded and preserved. The images were processed and the average temperature of canopy leaf was extracted. As shown in Fig. 2 that for different maize species at the seedling stage, the blade temperature in the drought group was higher than that of the watering group. Under drought stress, plans will resist drought stress within a certain range through a series of measures, mainly including stomatal regulation, osmotic adjustment, reactive oxygen removal, the protective effect of specific molecules and plant hormones regulation etc. 15. The direct reaction of plants to drought signal is to adjust the size of the stomatal aperture so as to prevent water loss. Under regular watering, maize stomata will be at the on-state and the water loss maintains the status quo. Under the drought stress, the water loss was much less compared with that with regular watering. Meanwhile, water loss will reduce the heat of leaves and lower the temperature. That is, the thermal loss under the drought stress was less than that under the watering. Hence, the blade temperature under the drought stress was lower than that in the watering group.

temperature difference between drought-resistant plants and non-drought resisting plants under the normal condition of watering

Infrared images acquired at 2:30 p.m. on October 31, 2012 were selected, and temperature of plants in the watering group was extracted as the following table 4. To better demonstrate the temperature differences, the table was showed in Fig. 3. From the above chart, in the watering group, plants with good drought resistance had higher blade temperature than those with poor drought resistance. The better the drought resistance, the higher the temperature of blades will be given the normal conditions.

temperature difference between drought-resistant plants and non-drought resisting plants under the drought stress

Infrared images acquired at 2:30 p.m. on October 31, 2012 were selected, and temperature of plants in the drought group was extracted as the following table 5. To better demonstrate the temperature differences, the table was showed in Fig. 6 . From the above chart, in the drought group, plants with good drought resistance had higher blade temperature than those with poor drought resistance. The better the drought resistance, the higher the temperature of blades will be given the normal conditions.

Drought is the main limiting factor in maize growth. Different maize species present different reactions under water stress [START_REF] Bi | Drought Resistance of Different Maize Varieties in Germination Stage[END_REF][START_REF] Bing | Identification of peanut drought-tolerance at seedling stage using infrared thermography[END_REF][START_REF] Ya | Application of Infrared Thermography in the Research of Plant Mechanism Response to Drought[END_REF][START_REF] Zhao Tianxin | The Study of Stomatal Conductance Estimation of Populus deltoids Bartr.x Populus ussuriensis Kom. By Infrared Thermography[END_REF][START_REF] Xuezheng | Application of Near Infrared Reflectance Spectroscopy in Determination of Drought Resistance of Winter Wheat[END_REF]. In this experiment, five maize species suitable for cultivation in the north were selected, and infrared images of plants of different species at the seedling stagewere acquired. As can be seen from This approach can analyze and evaluate the drought resistance of different maize species, thus providing a fast judgment and filtering maize species with properties of water-saving, high yield and high drought resistance. This is of great importance and meanings to alleviate the drought in the current Northern China and improve the current food and water shortage in China.

  at home and abroad achieved fruitful findings in the research of crop temperature. The contrastive analysis in the above studies was conducted based on leaf temperature and the conventional drought resistance parameters. In this study, blade temperature distribution characteristics of maize canopy were analyzed based on thermal infrared image to study on the temperature distribution characteristics at different leaf positions of different blades. In addition, different drought-resistant maize species were processed with water treatment to analyze the temperature variation of plants after the treatment of drought stress and to explore into the relationship between thermal infrared image characteristics of seedling stage and maize drought resistance.

  this study, we measured the temperature of leaf opex, leaf center and leaf base, with the data unit of℃ and the precision of 0.01℃. The right side of the figure shows the scale plate. Different colors show the corresponding temperature range. Before each measurement, the scale plate was adjusted so that the blade color can be easily distinguished compared with the color of background objects so that the measurement region can be found precisely.

Fig. 1 .

 1 Fig. 1. Infrared Image Processing of Maize at the Seedling Stage

Fig. 2 .

 2 Fig. 2. Temperature Comparison between the Watering Group and the Drought Group

Fig. 3 .Fig. 4 .

 34 Fig. 3. The graph of the temperature in the watering group

1 Fig. 5 .Fig. 6 .

 156 Fig. 5. The graph of the temperature in the drought group

Figure 4 ,

 4 Figure 4, in the watering group, plants with good drought resistance had higher blade temperature than those with poor drought resistance, that is, Kuancheng 60>Jingdan 68 > Jingke 528 > Dantiannuo No.3 > Yinnuo No.1, which formed a positive correlation with the drought resistance; as seen from Figure6, the positive correlation with the drought resistance was also formed under drought stress. Thus, given different water treatment, temperature of canopy blade of plants with good drought resistance was higher than that with poor drought resistance. The above analysis presented that temperature of canopy blade of plants at the seedling stage can reflect the drought resistance of plants, which was demonstrated in the different controlling capacities of plants with different drought resistance over the stomatal conductance.Plants with good drought resistance have stronger control over transpiration through stomata opening and closing. As less moisture was taken away by transpiration, the temperature was eventually elevated.

Table 1 .

 1 Information of different maize species

	Species	Drought resistance	Breeding Places
	KC 60	Good	Kuancheng Seed Co., Ltd, Hebei Province
	JD 68	Good	Maize Research Center, Beijing Academy of Agriculture and Forestry Sciences
	JK 528	Good	Maize Research Center, Beijing Academy of Agriculture and Forestry Sciences
	DTN 3	Poor	Horticulture Institute, Dandong Academy of Agricultural Sciences
	YN 1	Poor	Key Laboratory of Beijing Municipal Agricultural Application Technologies

Table 2 .

 2 Average Temperature of Different Blades at Different Blade Positions

	Average Temp. (℃)	leaf opex	leaf center	leaf base
	Blade 1	23.09	24.01	23.88
	Blade 2	23.94	25.12	24.01
	Blade 3	24.11	24.82	23.62
	Average Temp.	23.71	24.65	23.84
	Soil Temp.		19.53
	Stalk Temp.		23.09

Table 3 .

 3 Average temperature in the watering group and the drought group

	Species	Group	1	2	3	4	5	Aver a
	KC 60	the watering the drought	24.25 25.94	25. 81 26. 38	25. 56 26. 40	25. 69 26. 80	25. 15 26. 22	25.29 2 26.34 8
	JD 68	the watering the drought	25.58 26.64	25. 20 26. 67	25. 47 25. 77	24. 17 25. 92	25. 45 26. 25	25.17 4 26.25
	JK 528	the watering the drought	25.19 26.51	25. 34 25. 47	24. 96 25. 41	24. 82 25. 34	24. 68 25. 08	24.99 8 25.56 2
	DTN 3	the watering the drought	25.52 25.80	25. 28 25. 34	24. 82 25. 35	25. 28 25. 03	24. 72 25. 13	25.12 4 25.33
	YN 1	the watering the drought	24.84 24.58	24. 89 25. 62	24. 73 24. 32	24. 96 24. 65	24. 10 24. 17	24.70 4 24.66 8

Table 4 .

 4 Average temperature of canopy leaf of different varieties at seedling stage in the

			watering group
	Species\Temp.	1	2	3	4	5	Aver Temp.
	KC 60	26.33	27.49	27.60	27.62	27.76	27.36
	JD 68	26.73	27.25	27.12	27.40	27.30	27.16
	JK 528	27.02	26.88	26.57	26.98	27.06	26.902
	DTN 3	26.86	26.24	26.37	26.63	26.30	26.48
	YN 1	26.50	26.35	26.32	26.42	25.64	26.246

Table 5 .

 5 Average temperature of canopy leaf of different species at seedling stage in the

			drought group
	Species \ Temp.	1	2	3	4	5	Aver
							Temp.
	KC 60	27.51	27.69	27.01	27.04	26.31	27.112
	JD 68	27.68	26.60	26.19	26.29	25.97	26.546
	JK 528	27.09	26.24	25.91	25.90	25.39	26.106
	DTN 3	26.70	26.36	26.21	26.22	25.59	26.216
	YN 1	25.91	25.88	25.95	25.47	25.70	25.782

Acknowledgment

This research was financially supported by Natural Science Foundation of China (31201125) and Beijing Nova Program (Z111105054511051).


Alice Le Brigant

Marc Arnaudon

Frédéric Barbaresco

Reparameterization invariant metric on the space of curves

 to define a reparameterization invariant metric on the space of immersions M = Imm([0, 1], M ) by pullback of a metric on the tangent bundle TM derived from the Sasaki metric. We observe that such a natural choice of Riemannian metric on TM induces a first-order Sobolev metric on M with an extra term involving the origins, and leads to a distance which takes into account the distance between the origins and the distance between the SRV representations of the curves. The geodesic equations for this metric are given, as well as an idea of how to compute the exponential map for observed trajectories in applications. This provides a generalized theoretical SRV framework for curves lying in a general manifold M .

Introduction

Computing distances between shapes of open or closed curves is of interest in many fields that require shape analysis, from medical imaging to video surveillance, to radar detection. While the shape of an organ or a human contour can be modeled by a closed plane curve, some applications require the manipulation of curves lying in a non flat manifold, such as S 2 -valued curves representing trajectories on the earth or curves in the space of hermitian positive definite matrices, where the values represent covariance matrices of Gaussian processes. The shape space of planar curves has been widely studied ( [START_REF] Michor | Riemannian geometries on spaces of plane curves[END_REF], [START_REF] Michor | An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach[END_REF], [START_REF] Younes | A Metric on Shape Space with Explicit Geodesics[END_REF], [START_REF] Bauer | Constructing reparameterization invariant metrics on spaces of plane curves[END_REF]), and the more general setting of shapes lying in any manifold M has recently met great interest ( [START_REF] Bauer | Sobolev metrics on shape space of surfaces[END_REF], [START_REF] Su | Statistical analysis of trajectories on Riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF], [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF], [START_REF] Zhang | Video-based action recognition using rateinvariant analysis of covariance trajectories[END_REF]). Here we consider open oriented curves in a Riemannian manifold M , more precisely the space of immersions c : [0, 1] → M , M = Imm([0, 1], M ).

Reparameterizations will be represented by increasing diffeomorphisms φ : [0, 1] → [0, 1] (so that they preserve the end points of the curves), and their set is denoted by Diff + ([0, 1]). Then, one way to describe a shape is as the equivalence class of all the curves that are identical modulo reparameterization, and the shape space as the associated quotient space,

S = Imm([0, 1], M )/Diff + ([0, 1]).
The formal principal bundle structure π : M → S induces a decomposition of the tangent bundle T M = V M ⊕ H M into a vertical subspace V M = ker(T π) consisting of all vectors tangent to the fibers of M over S , and a horizontal subspace H M = (V M ) ⊥ G defined as the orthogonal complement of V M according to the metric G that we put on M . We say formal because the manifold structure of the space Imm([0, 1], M ) has not yet been thoroughly studied to our knowledge. We require that G be reparameterization invariant, that is to say that the action of Diff + ([0, 1]) be isometric for

G G c•φ (h • φ, k • φ) = G c (h, k), (1)
for any curve c ∈ M , reparameterization φ ∈ Diff + ([0, 1]), and infinitesimal deformations h, k ∈ T c Mh and k can also be seen as vector fields along the curve c in M . That way, the induced geodesic distance between two curves c 0 and c 1 does not change if we reparameterize them the same way, that is

d (c o • φ, c 1 • φ) = d (c 0 , c 1 ),
for any φ ∈ Diff + ([0, 1]). Also, if this property is satisfied, then G induces a Riemannian metric Ĝ on the shape space, Ĝπ(c

) (T c π(h), T c π(k)) = G c (h H , k H ),
in the sense that the above expression does not depend on the choice of the representatives c, h and k. Here h H , k H denote the horizontal parts of h and k according to the previously mentioned decomposition, as well as the horizontal lifts of T c π(h) and T c π(k), respectively. The geodesic distances d on M and d on S are then simply linked by d

( [c 0 ] , [c 1 ] ) = inf d c 0 , c 1 • φ | φ ∈ Diff + ([0, 1]) ,
where [c 0 ] and [c 1 ] denote the shapes of two given curves c 0 and c 1 , and d verifies the stronger property

d (c 0 • φ, c 1 • ψ) = d (c 0 , c 1 ),
for any reparameterizations φ, ψ ∈ Diff + ([0, 1]). The most natural candidate for a reparameterization invariant metric G on M is the L 2 -metric with integration over arc length, but Michor and Mumford have shown in [START_REF] Michor | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF] that the induced metric Ĝ on the shape space always vanishes. This has motivated the study of Sobolev metrics ( [START_REF] Michor | An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach[END_REF], [START_REF] Bauer | Constructing reparameterization invariant metrics on spaces of plane curves[END_REF], [START_REF] Bauer | Why use Sobolev metrics on the space of curves[END_REF]), and particularly of a first-order Sobolev metric on the space of plane curves,

G c (h, k) = 〈D h ⊥ , D k ⊥ 〉 + 1 4 〈D h , D k 〉 d , (2)
where we integrate according to arc length d = c (t ) dt and 〈•, •〉 denotes the euclidean metric on IR 2 , D h = 1 c h is the derivation of h according to arc length, D s h = 〈D s h, v〉v is the projection of D s h on the unit length tangent vector field v = 1 c c along c, and D s h ⊥ = 〈D s h, n〉n is the projection of D s h on the unit length normal vector field n along c. This particular first-order Sobolev metric is of interest because it can be studied via the square root velocity (SRV) framework, introduced by Srivastava et al. in [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF] and used in several applications ( [START_REF] Laga | Landmark-free statistical analysis of the shape of plant leaves[END_REF], [START_REF] Su | Statistical analysis of trajectories on Riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF]). This framework can be extended to curves in a general manifold by using parallel transport, in a way which allows us to move the computations to the tangent plane to the origin of one of the two curves under comparison, see [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF] and [START_REF] Zhang | Video-based action recognition using rateinvariant analysis of covariance trajectories[END_REF]. In [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF] the transformation used is a generalization of the SRV function introduced by Bauer et al. in [START_REF] Bauer | Constructing reparameterization invariant metrics on spaces of plane curves[END_REF] as a tool to study a more general form of the Sobolev metric [START_REF] Bauer | Why use Sobolev metrics on the space of curves[END_REF]. In [START_REF] Zhang | Video-based action recognition using rateinvariant analysis of covariance trajectories[END_REF] a Riemannian framework is given, including the associated Riemannian metric and the geodesic equations. While our approach in this paper is similar, we feel that the distance we introduce here will be more directly dependent on the "relief" of the manifold, since it is computed in the manifold itself rather than in one tangent plane as in [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF] and [START_REF] Zhang | Video-based action recognition using rateinvariant analysis of covariance trajectories[END_REF]. This enables us to take into account a greater amount of information on the space separating two curves.

New metric on the space of parameterized curves

We consider the square root velocity function (SRVF) introduced in [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF] on the space of curves in M ,

R : M → T M , c → c c ,
where • is the norm associated to the Riemannian metric on M . This function will allow us to define a metric G on M by pullback of a metric G on T M . First, we define the following projections from T T M to T M . Let ξ ∈ T (p,u) T M and (x,U ) be a curve in T M that passes through (p, u) at time 0 at speed ξ. Then we define the vertical and horizontal projections vp (p,u) :

T (p,u) T M → T p M , ξ → ξ V := ∇ x (0) U , hp (p,u) : T (p,u) T M → T p M , ξ → ξ H := x (0).
The horizontal and vertical projections live in the tangent bundle T M and are not to be confused with the horizontal and vertical parts which live in the double tangent bundle T T M and will be denoted by ξ H , ξ V . Furthermore, let us point out that the horizontal projection is simply the differential of the natural projection T M → M , and that according to these definitions, the Sasaki metric ( [START_REF] Sasaki | On the differential geometry of tangent bundles of Riemannian manifolds[END_REF], [START_REF] Sasaki | On the differential geometry of tangent bundles of Riemannian manifolds II[END_REF]) can be written

g S (p,u) (ξ, η) = ξ H , η H + ξ V , η V ,
where 〈•, •〉 denotes the Riemannian metric on M . Now we can define the metric that we put on T M . Let us consider h ∈ T M and ξ, η ∈ T h T M . We define

Gh ξ, η = ξ(0) H , η(0) H + 1 0 ξ(t ) V , η(t ) V dt , (3)
where ξ(t ) H = hp(ξ(t )) and ξ(t ) V = vp(ξ(t )) are the horizontal and vertical projections of ξ(t ) ∈ T T M for all t . Then we have the following result.

Proposition 1. The pullback of the metric G by the square root velocity function R is given by

G c (h, k) = 〈 h(0), k(0) 〉 + ∇ h ⊥ , ∇ k ⊥ + 1 4 ∇ h , ∇ k d , (4)
for any curve c ∈ M and vectors h, k ∈ T c M , where we integrate according to arc length, ∇ h = 1 c ∇ c h is the covariant derivative of h according to arc length, and ∇ h = 〈D h, v〉v and ∇ h ⊥ = ∇ h -∇ h are its tangential and normal components respectively. Remark 1. In the case of curves in a flat space, G is the first-order Sobolev metric (2), studied in [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF], with an added term involving the origins. This extra term guaranties that the induced distance is always greater than the distance between the starting points of the curves in M .

Proof. For any c ∈ M , and h, k ∈ T c M , the metric G is defined by

G c (h, k) = GR(c) (T c R(h), T c R(k)) . For any t ∈ [0, 1], we have T c R(h)(t ) H = h(t ) and T c R(h) V = ∇ h R(c)(t ).
To prove this proposition, we just need to compute the latter. Let s → c(s, •) be a curve in M such that c(0, •) = c and c s (0, •) = h . Here and in all the paper we use the notations c s = ∂c/∂s and c t = ∂c/∂t . Then

∇ h R(c)(t ) = 1 c 1/2 ∇ h c + h c -1/2 c = 1 c t 1/2 ∇ s c t + ∂ s 〈 c t , c t 〉 -1/4 c t = 1 c t 1/2 ∇ t c s - 1 2 〈 c t , c t 〉 -5/4 〈 ∇ s c t , c t 〉 c t = c 1/2 (∇ h) ⊥ + 1 2 〈∇ h , c c 〉 c c ,
where in the last step we use again the inversion ∇ s c t = ∇ t c s .

Fiber bundle structures

Principal bundle over the shape space We already know that we have a formal principal bundle structure over the shape space

π : M = Imm([0, 1], M ) → S = M /Diff + ([0, 1]). which induces a decomposition T M = V M ⊥ ⊕ H M .
Just as in the planar case, the fact that the square root velocity function R verifies the equivariance property

R(c • φ) = φ R(c) • φ for all c ∈ M , h, k ∈ T c M and φ ∈ Diff + ([0, 1]), guaranties that the integral part of G is reparameterization invariant. Remembering that the reparameterizations φ ∈ Diff + ([0, 1]
) preserve the origins of the curves, we notice that G is constant along the fibers, as expressed in equation (1), and so there exists a Riemannian metric Ĝ on the shape space S such that π is (formally) a Riemannian submersion from (M ,G) to (S , Ĝ)

G c (h H , k H ) = Ĝπ(c) (T c π(h), T c π(k)) ,
where h H and k H are the horizontal parts of h and k respectively.

Fiber bundle over the starting points

The special role that plays the starting point in the metric G induces another formal fiber bundle structure, where the base space is the manifold M , seen as the set of starting points of the curves, and the fibers are the set of curves with the same origin. The projection is then

π ( * ) : M → M , c → c(0).
It induces another decomposition of the tangent bundle in vertical and horizontal bundles

V ( * ) c M = ker T π ( * ) = { h ∈ T c M | h(0) = 0 } , H ( * ) c M = V ( * ) c M ⊥ G . Proposition 2. We have the usual decomposition T M = V ( * ) M ⊥ ⊕ H ( * ) M , the hori- zontal bundle H ( * )
c M consists of parallel vector fields along c, and π ( * ) is (formally) a Riemannian submersion for (M ,G) and (M , 〈•, •〉).

Proof. Let h be a tangent vector. Consider h 0 the parallel vector field along c with initial value h 0 (0) = h(0). It is a horizontal vector, since its vanishing covariant derivative along c assures that for any vertical vector l we have G c (h 0 , l ) = 0. The difference h = hh 0 between those two horizontal vectors has initial value 0 and so it is a vertical vector, which gives a decomposition of h into a horizontal vector and a vertical vector. The definition of H ( * ) M as the orthogonal complement of V ( * ) M guaranties that their sum is direct. Now if k is another tangent vector, then the scalar product between their horizontal parts is

G c (h H , k H ) = h H (0) , k H (0) c(0) = 〈 h(0) , k(0) 〉 c(0) = T c π ( * ) (h H ) , T c π ( * ) (k H ) π ( * ) ,
and this completes the proof.

Induced distance on the space of curves

Here we will give an expression for the geodesic distance induced by the metric G. Let us consider two curves c 0 , c 1 ∈ M , and a path of curves s → c(s,

•) linking them in M c(0, t ) = c 0 (t ), c(1, t ) = c 1 (t ),
for all t ∈ [0, 1]. We denote by q(s, •) = R (c(s, •)) the image of this path of curves by the SRVF R. Note that q is a vector field along the surface c in M . Let now q be the raising of q in the tangent plane T c(0,0) M in the following way q(s, t

) = P s,0 c(•,0) • P t ,0 c(s,•) q(s, t ) ,
where we denote by P t 1 ,t 2 γ

: T γ(t 1 ) M → T γ(t 2 ) M the parallel transport along a curve γ from γ(t 1 ) to γ(t 2 ). Notice that q is a surface in a vector space, as illustrated in Figure 1. Lastly, we introduce a vector field (a, τ) → ω s,t (a, τ) in M , which parallel translates q(s, t ) along c(s, •) to its origin, then along c(•, 0) and back down again, as shown in Figure 1. More precisely

ω s,t (a, τ) = P 0,τ c(a,•) • P s,a c(•,0) • P t ,0 c(s,•) q(s, t )
for all b, s. That way the quantity ∇ s ω s,t measures the holonomy along the rectangle of infinitesimal width shown in Figure 1. We can now formulate our result.

Proposition 3.

With the above notations, the geodesic distance induced by the Riemannian metric G between two curves c 0 and c 1 on the space

M = Imm([0, 1], M ) of parameterized curves is given by dist(c 0 , c 1 ) = inf c(0,•)=c 0 ,c(1,•)=c 1 1 0 c s (s, 0) 2 + 1 0 ∇ s q(s, t ) 2 dt ds, ( 5
)
where q = R(c) is the Square Root Velocity representation of the curve c and the norm is the one associated to the Riemannian metric on M . It can also be written

dist(c 0 , c 1 ) = inf c(0,•)=c 0 ,c(1,•)=c 1 1 0 c s (s, 0) 2 + 1 0 qs (s, t ) + Ω(s, t ) 2 dt ds, ( 6
)
where q is the raising of q in the tangent plane T c(0,0) M and the curvature term Ω is given by

Ω(s, t ) = P s,0 c(•,0) • P t ,0 c(s,•) ∇ s ω s,t (s, t ) = P s,0 c(•,0) • P t ,0 c(s,•) t 0 P τ,t c(s,•) R(c τ , c s )P t ,τ c(s,•) q(s, t ) dτ ,
if R denotes the curvature tensor of the manifold M .

Remark 2. Our original motivation for this work was to find a geodesic distance (that is, a distance induced by a Riemannian metric) that resembled the product distance introduced in [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF]. In the first term under the square root of expression (6) we can see the velocity vector of the curve c(•, 0) linking the two origins, and in the second the velocity vector of the curve q linking the TSRVF-images of the curves -Transported Square Root Velocity Function, as introduced by Su et al. in [START_REF] Su | Statistical analysis of trajectories on Riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF]. However there is also a curvature term Ω which, as previously mentionned, measures the holonomy along the rectangle of infinitesimal width shown in Figure 1. If instead we equip the tangent bundle TM with the metric

Gh (ξ, ξ) = ξ h (0) 2 + 1 0 ξ v (t ) - t 0 P τ,t c R(c , ξ h )P t ,τ c h(t ) dτ 2 dt ,
for h ∈ T M and ξ, η ∈ T h T M , then the curvature term Ω vanishes and the geodesic distance on M becomes dist(c 0 ,

c 1 ) = inf c(0,•)=c 0 ,c(1,•)=c 1 1 0 c s (s, 0) 2 + qs (s, •) 2 ds, (7)
where the norm of the second term under the square root is the L 2 -norm, and which corresponds exactly to the geodesic distance associated to the metric on the space C = ∪ p∈M L 2 ([0, 1], T p M ) introduced by Zhang et al. in [START_REF] Zhang | Video-based action recognition using rateinvariant analysis of covariance trajectories[END_REF]. The difference between the two distances ( 5) and ( 7) resides in the curvature term Ω, which translates the fact that in the first one, we compute the distance in the manifold, whereas in the second, it is computed in the tangent space to one of the origins of the curves. Therefore, the first one takes more directly into account the "relief" of the manifold between the two curves under comparison. For example, if there is a "bump" between two curves in an otherwise relatively flat space, the second distance ( 7) might not see it, whereas the first one (5) will thanks to the curvature term.

Remark 3. Let us briefly consider the flat case : if the manifold M is flat, the two distances ( 5) and ( 7) coincide. If two curves c 0 and c 1 in a flat space have the same starting point p, the first summand under the square root vanishes and the distance becomes the L 2 -distance between the two SRV representations q 0 = R(c 0 ) and q 1 = R(c 1 ). If two curves in a flat space differ only by a translation, then the distance is simply the distance between their origins. Proof. Since G is defined by pullback of G by the SRVF R, we know that the lengths of c in M and of q = R(c) in TM are equal and so that

dist(c 0 , c 1 ) = inf c(0,•)=c 0 ,c(1,•)=c 1 1 0 G q s (s, •), q s (s, •) ds, with G q s (s, •), q s (s, •) = c s (s, 0) 2 + 1 0 ∇ s q(s, t ) 2 dt .
Now let us fix t ∈ [0, 1]. Then s → P t ,0 c(s,•) q(s, t ) is a vector field along c(•, 0), and so

∇ s P t ,0 c(s,•) q(s, t ) = P 0,s c(•,0) ∂ ∂s P s,0 c(•,0) • P t ,0 c(s,•) q(s, t ) = P 0,s c(•,0) qs (s, t ) .
We consider the vector field ν along the surface (s, τ) → c(s, τ) that is parallel along all curves c(s, •) and takes value ν(s, t ) = q(s, t ) in τ = t for any s ∈ [0, 1], that is

ν(s, τ) = P t ,τ c(s,•) q(s, t ) ,
for all s ∈ [0, 1] and τ ∈ [0, 1]. That way we know that

∇ s ν(s, t ) = ∇ s q(s, t ), ∇ s ν(s, 0) = P 0,s c(•,0) qs (s, t ) , ∇ τ ν(s, τ) = 0,
for all s, τ ∈ [0, 1]. Then we can express its covariant derivative in the following way

∇ s ν(s, t ) = P 0,t c(s,•) (∇ s ν(s, 0)) + t 0 P τ,t c(s,•) (∇ τ ∇ s ν(s, τ)) d τ = P 0,t c(s,•) • P 0,s c(•,0) qs (s, t ) + t 0 P τ,t c(s,•) R(c τ , c s )P t ,τ c(s,•) q(s, t ) dτ. ( 8
)
Now let us fix s ∈ [0, 1] as well. Notice that the vector field ω s,t defined above verifies ω s,t (s, t ) = q(s, t ),

∇ τ ω s,t (a, τ) = 0, ∇ a ω s,t (a, 0) = 0,
for all a, τ ∈ [0, 1]. Note that unlike ν, we do not have ∇ s ω s,t (s, t ) = ∇ s q(s, t ) because ω s,t (a, t ) = q(a, t ) is only true for a = s. It is easy to verify that the last term of equation ( 8) is precisely the covariant derivative of the vector field ω s,t

∇ s ω s,t (s, t ) = t 0 P τ,t c(s,•) R(c τ , c s )P t ,τ c(s,•) q(s, t ) dτ,
since for any τ ∈ [0, 1], ω s,t (s, τ) = P t ,τ c(s,•) q(s, t ), and finally by composing by P s,0 c(•,0) • P t ,0 c(s,•) , we obtain the second expression [START_REF] Michor | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF], which completes the proof.

Geodesic equation on M

In order to be able to compute the distance given by (5) between two curves, we first need to compute the optimal deformation s → c(s, •) from one to the other. In other words, we need to characterize the geodesics of M for our metric. In order to do so, we use the variational principle. The beginning of the calculations are very similar to those in [START_REF] Zhang | Video-based action recognition using rateinvariant analysis of covariance trajectories[END_REF]. Let us consider two curves c 0 , c 1 ∈ M and a path [0, 1] s → c(s, •) ∈ M going from one to the other. This path c is a geodesic if it minimizes the energy functional

E : M → R + E (c) = 1 0 G ∂c ∂s , ∂c ∂s ds.
Let a → c(a, •, •), a ∈ (-, ), be a proper variation of the path s → c(s, •), meaning that it coincides with c in a = 0, and it preserves its end points

c(0, •, •) = c, c a (a, 0, •) = 0 ∀a ∈ (-, ), c a (a, 1, •) = 0 ∀a ∈ (-, ). Then c is a geodesic of M if d d a a=0 E (c(a, •, •)) = 0 for any proper variation a → c(a, •, •). If we denote by E (a) = E (c(a, •, •))
, for a ∈ (-, ), the energy of a proper variation of c, then we have

E (a) = 〈 c s (a, s, 0), c s (a, s, 0) 〉 ds + ∇ s q(s, t ), ∇ s q(s, t ) dt ds,
where q = c t / c t is the SRV representation of c. Its derivative is given by

1 2 E (a) = 〈 ∇ a c s (a, s, 0), c s (a, s, 0) 〉 ds + ∇ a ∇ s q(a, s, t ), ∇ s q(a, s, t ) dt ds.
Considering that the variation preserves the end points, integration by parts gives 〈 ∇ a c s , c s 〉 ds = -〈∇ s c s , c a 〉 ds ∇ s ∇ a q, ∇ s q ds = -∇ s ∇ s q, ∇ a q ds, and so we obtain

1 2 E (a) = -〈 ∇ s c s , c a 〉| t =0 ds + R(c a , c s )q + ∇ s ∇ a q, ∇ s q dt ds = -〈 ∇ s c s , c a 〉| t =0 ds - R(q, ∇ s q)c s , c a + ∇ s ∇ s q, ∇ a q dt ds.
This quantity has to vanish in a = 0 for all proper variations

∇ s c s | t =0 , c a | a=0,t =0 ds+ R(v, ∇ s q)c s , c a | a=0 + ∇ s ∇ s q, ∇ a q a=0 dt ds = 0.
Unfortunately, we cannot yield any conclusions at this point, because c a (0, s, t ) and ∇ a q(0, s, t ) cannot be chosen independently, since q is not any vector field along c but its image via the Square Root Velocity Function. Using the notations B (s, t ) = R(q, ∇ s q)c s (s, t ),

D(s, t ) = 1 |c t | ∇ s ∇ s q(s, t ) - 1 2 ∇ s ∇ s q , c t |c t | 5/2 c t (s, t ),
and the lighter notation u(t 1 ) t 1 ,t 2 = P t 1 ,t 2 c

(u(t 1 )) to denote the parallel transport along a curve c of a vector field u, the following simple manipulations provide us with a solution B (s, t ), c a (0, s, 0) 0,t + t 0 ∇ t c a (0, s, τ) τ,t dτ dt ds

+ 1 0 1 0 〈D(s, t ), ∇ t c a (0, s, t )〉 dt ds = 1 0 ∇ s c s (s, 0) + 1 0 B (s, t ) t ,0 dt , c a (0, s, 0) ds + 1 0 1 0 t 0 B (s, t ) t ,τ , ∇ t c a (0, s, τ) dτ dt ds + 1 0 1 0 〈D(s, t ), ∇ t c a (0, s, t )〉 dt ds = 1 0 ∇ s c s (s, 0) + 1 0 B (s, t ) t ,0 dt , c a (0, s, 0) ds + 1 0 1 0 D(s, t ) + 1 t B (s, τ) τ,t dτ , ∇ t c a (0, s, t ) dt ds.
Since the variations c a (0, s, 0) and ∇ t c a (0, s, •) can be chosen independently for all s, we obtain the following characterization of the geodesics s → c(s, •) of M ∇ s c s (s, 0) + 1 0 R(q, ∇ s q)c s (s, t ) t ,0 dt = 0, ∀s D(s, t ) + 1 t R(q, ∇ s q)c s (s, τ) τ,t dτ = 0, ∀t , s.

Finally, using the definition of D(s, t ), the geodesic equations are

∇ s c s (s, 0) = r (s, 0), ∀s (9a)
∇ s ∇ s q(s, t ) = q(s, t ) r (s, t ) + r (s, t ) , ∀t , s (9b

)
where q is the SRV representation of c, the vector field r is given by r (s, t ) = -1 t R(q, ∇ s q)c s (s, τ) τ,t dτ, and r = 〈r, v〉 v with v = 1 c t c t , is the tangential component of r .

Exponential map

In this section, we describe an algorithm which allows us to compute the geodesic s → c(s, •) starting from a point c ∈ M at speed u ∈ T c M . This amounts to finding the optimal deformation of the curve c in the direction of the vector field u according to our metric. We initialize this path s → c(s, •) by setting c(0, •) = c and c s (s, 0) = u, and we propagate it using iterations of fixed step > 0. The aim is, given c(s, •) and c s (s, •), to deduce c(s + , •) and c s (s + , •). The first is obtained by following the exponential map on the manifold M c(s

+ , •) = exp c(s,•) c s (s, •),
and the second requires the computation of the variation ∇ s c s (s, •)

c s (s + , •) = [c s (s, •) + ∇ s c s (s, •)] s,s+ ,
where once again, we use the notation w(s) s,s+ = P s,s+ c (w(s)) for the parallel transport of a vector field s → w(s) along a curve s → c(s) in M . If we suppose that at time s we have c(s, •) and c s (s, •), then we also know c t (s, •) and q(s, •) = 1 c t c t (s, •), as well as ∇ t c s (s, •) and

∇ s q(s, •) = ∇ s c t |c t | (s, •) - 1 2 〈 ∇ s c t , c t 〉 |c t | 5/2 c t (s, •), (10)
using the fact that ∇ s c t = ∇ t c s . The variation ∇ s c s (s, •) can then be computed in the following way

∇ s c s (s, t ) = ∇ s c s (s, 0) 0,t + t 0 [∇ s ∇ s c t (s, τ) -R(c t , c s )c s (s, τ)] τ,t dτ for all t ∈ [0, 1],
where ∇ s c s (s, 0) is given by equation (9a), the second order variation ∇ s ∇ s c t (s, •) is given by

∇ s ∇ s c t = |c t | 1/2 ∇ s ∇ s q + 〈∇ t c s , c t 〉 |c t | 2 ∇ t c s + ∇ s ∇ s q, c t |c t | 3/2 - 3 2
〈∇ t c s , c t 〉 2 |c t | 4 + |∇ t c s | 2 |c t | 2 c t , (11)
and ∇ s ∇ s q can be computed via equation (9b). Now, if we no longer have a continuous curve but a series of discrete observations p 0 , p 1 , . . . , p n ∈ M made at discrete times 0 = t 0 < t 1 < . . . < t n = 1, then the optimal deformation of these points in the direction of a series of tangent vectors u 0 , u 1 , . . . , u n , where u k ∈ T p k M for k = 0, . . . , n, can be computed by the following steps.

Algorithm 1 (Discrete Exponential Map)

Initialization : Set c(0, t k ) = p k and c s (0, t k ) = u k for k = 0, . . . , n. Heredity : If c(s, t k ) and c s (s, t k ) are known for all k, then

1. for all 0 ≤ k ≤ n -1, set c t (s, t k ) = 1 t k+1 -t k log c(s,t k ) c(s, t k+1 ), ∇ t c s (s, t k ) = 1 t k+1 -t k c s (s, t k+1 ) t k+1 ,t k -c s (s, t k ) ,
where log denotes the inverse of the exponential map on M , and compute q(s, t k ) = 1 c t c t (s, t k ) and ∇ s q(s, t k ) using equation [START_REF] Sasaki | On the differential geometry of tangent bundles of Riemannian manifolds II[END_REF].

2. Compute r (s, t k ) = n =k R(q, ∇ s q)c s (s, t ) t ,t k for all 0 ≤ k ≤ n and set ∇ s c s (s, t 0 ) = r (s, t 0 ), ∇ s ∇ s q(s, t k ) = q(s, t k ) r (s, t k ) + r (s, t k ) ∥ .

3. For k = 0, . . . , n -1, (a) compute ∇ s ∇ s c t (s, t k ) using equation [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF]

Conclusion

In the same way that the first-order Sobolev metric (2) on the space of plane curves can be obtained as the pullback of the L 2 -metric by the square root velocity function ( [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF]), our metric G can be obtained as the pullback of a natural metric G on the tangent bundle TM by the same SRVF. As such it is reparameterization invariant, and induces a Riemannian metric Ĝ on the shape space S for which the fiber bundle projection is formally a Riemannian submersion. On the other hand, the special role that G gives to the starting points of the curves induces another formal fiber bundle structure, this time over the manifold M seen as the set of starting points of the curves, for which the projection is formally also a Riemannian submersion. The geodesic distance induced by G takes into account the distance between the origins of the curve and the L 2 -distance between the SRV representations, but instead of transporting the computations in a unique tangent plane as in [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF] and [START_REF] Zhang | Video-based action recognition using rateinvariant analysis of covariance trajectories[END_REF], we stay in the manifold. This should allow us to take into account a greater amount of information on its geometry. Finally, explicit equations can be obtained for the geodesics on the space M of parameterized curves for our metric G, and the exponential map can be iteratively computed. Future work will include some numerical simulations to illustrate the work of this paper.

Fig. 1 .

 1 Fig. 1. Illustration of the distance between two curves c 0 and c 1 in the space of curves M

1 0 1 0 1 0 1 0

 1111 〈∇ s c s (s, 0), c a (0, s, 0)〉 ds + 〈B (s, t ), c a (0, s, t )〉 + 〈D(s, t ), ∇ t c a (0, s, t )〉 dt ds = 〈∇ s c s (s, 0), c a (0, s, 0)〉 ds +

  , (b) compute ∇ s c s (s, t k+1 ) using∇ s c s (s, t k+1 ) = [ ∇ s c s (s, t k ) + (t k+1t k ) ( ∇ s ∇ s c t (s, t k ) -R(c s , c t )c s (s, t k ) ) ] t k ,t k+1 .4. Finally, for all0 ≤ k ≤ n, set c(s + , t k ) = exp c(s,t k ) ( c s (s, t k )) c s (s + , t k ) = [ c s (s, t k ) + ∇ s c s (s, t k ) ] t k ,t k+1 .

Acknowledgments

This research was supported by Thales Air Systems and the french MoD DGA.


Weiming Liu

Wang Enguo

Huang Shanyou

He Xianbiao

Keywords: hybrid rice, Yongyou 8, yield trait, path analysis, technical strategy

In order to further clarify the path to high yield cultivation techniques for Yongyou 8, a hybrid late season japonica rice line, we selected 65 sets of data from 11 sites in Zhejiang and Jiangsu Provinces between 2005 and 2011, analyzed profiles of panicle, kernel and weight in Yongyou 8. Through correlation, regression and path analysis, a technical strategy has been proposed for high yield Yongyou 8 cultivation: appropriate control of the number of effective tillers, focus on large panicles, with appropriate attention to seed setting rate and kernel weight.

Introduction

Simple correlation [START_REF] Guo | Correlation analysis and QTL mapping of grain shape and grain weight in rice under upland and lowland environments[J][END_REF][START_REF] Zeng | Canonical correlation analysis between quality and yield characters of hybrid rice[END_REF][START_REF] Zhang | The correlation of rice yield and its components[J][END_REF]or multiple regression [START_REF] Sun | Regressive analysis between leaf characters and yield factors in elongation stage in rice[J][END_REF][START_REF] Zhao | Secondary compounds in rice varieties resistant to Nilaparvata lugen[J]s[END_REF][START_REF] Zhang | Differences of nitrogen absorption and utilization and their influences on grain yield of conventional Indica rice cultivars[J][END_REF][START_REF] Guo | Genetic analysis and utilization of the important agronomic traits on Zhenshan97×Minghui 63 recombinant inbred lines(RIL) in rice (Oryza sativa L.) [J][END_REF]analysis are often employed to study the relationship between plant main characters and crop yield often using. However, simple correlation cannot fully examine the relationship between variables, so that the results might have a certain one-sidedness (Ref??). Multivariate regression analysis, to a certain extent, can eliminate confusion between variables and reveal the real correlation of the independent variables with their dependent variables. The path analysis is able to effectively demonstrate a direct effect of the relevant variables on the outcome and to find out the indirect effect of casual factors on the affected factors, allowing comparison of the relative importance of various causal factors [START_REF] Li | Analysis of main agronomic characters of rice in Ningxia[J][END_REF][START_REF] Cheng | Relevance and path analysis on harvest index and single yield and its composition factors of rice[END_REF][START_REF] Wang | Correlation and path analysis of economic characters of hybrid rice[END_REF]. Such a comprehensive analysis will be a solid ground for improving yield cultivation techniques.

In order to determine cultivation strategies for high yield of intermediate to late season japonica rice, we explored the effect of major economical traits on yield [START_REF] Liu | Correlation, Multiple Regression and path Analysis between Yield Traits and Yield on Intersubspecific Hybrid Rice[J][END_REF][START_REF] Liu | Analysis of high yield and efi ciency technique in intersubspecific hybrid rice[J][END_REF][START_REF] Liu | Study on the High Yield and High Efficiency Cultivation Technique of Indica-japonica Intersubspecific Hybrid Rice as Single-cropping Late Rice[END_REF]. Yongyou 8 , en elite and late season cultivar, is often used as a three-strain hybrid rice line. In 2012, we collated 55 sets of data from 2005 to 2011 in 9 sites of Zhejiang and Jiangsu Provinces, to determinethe relationships between traits and yield, performed regression, correlation and path analysis [START_REF] Liu | Analysis of High Yield and Efficiency Technique on Yongyou 8 Hybrid Rice[END_REF] for Yongyou 8. Later, we have gathered more economical trait data from various sites, and further analyzed the effects of major economical traits on yield of Yongyou 8.

Materials and Methods

Data Collection

Yongyou 8 is a 3-strain hybrid late season japonica rice line with moderate growth period, can be grown late season in areas where rice is harvested once or twice a year. Yongyou 8 has loose tillers, moderate tillering capability, high panicle bearing rate, large spikes, high seed setting rate, moderate lodging resistance. The quality of rice meets the national standard of second excellence grade.

Results from analysis of 65 sets of field data for 7 continuous years from 11 sites in the two Provinces (Table 1) indicate that the yield was negatively correlated with the number of effective tillers per unit area, but positively correlated with the number of total kernels, filled kernels per panicle and 1,000 kernel weight. The seed setting rate was the highest in paddies with a yield at 11250 kg/ha, followed by paddies with a yield at 8250 kg/hm 2 .

Method of Analysis

Sixty five sets of trait data were analyzed for correlation, regression and path to reveal the cultivation technical strategy to high yield.

Results and Analysis

Regression Analysis among Traits

As shown in Table 2, among five traits, yield is negatively correlated with the number of effective tillers not statistically insignificant though, and positively correlated with the remainder 4 traits very significantly. Further analysis of the five traits revealed that the number of effective tillers was negatively correlated at very significant level with total number of kernels per panicle, filled kernels and 1,000 kernel weight, weakly correlated with seed setting rate. The total number of kernels per panicle was very significantly correlated with filled kernels, 1,000 kernel weight, weakly negatively correlated with seed setting rate. Thus, relationships among the five traits are complex, and their relationships with yield vary even negative. Thus, it is the key to understand how to coordinate traits with yield, how to focus on the main factor to fully explore the traits. Note: * denoted difference at α = 0.05. ** denotes deference at α = 0.01.

Gradual correlation analysis

We analyzed 65 sets of data for traits by stepwise regression, maximized coefficient. The regression equation is x = seed setting rate, and 5 x = 1,000 kernel weight.

Through results of stepwise regression analysis, it can be seen that among the five traits, the total kernels per panicle did not contribute to yield. The ranges of the five traits were 1.7700~3.2850 million/hm 2 , 149.60~288.90 kernels/panicle, 114.27~250.74 kernels/panicle, 65.50~89.20% and 25.30~29.32 g for x 1 , x 2 , x 3 , x 4 , x 5 , respectively. The range of yield (y) was 7665.00~13549.50 kg/hm 2 .

Results of variance analysis indicated that the F value for the regression equation was 75.60 with P < 0.01. Analysis of partial correlation revealed that filled kernels per panicle (r (y, x3) = 0.7003) had the largest contribution to yield, followed by effectively panicles (r (y, x 1 ) = 0.6672), weight per thousand kernel (r (y, x 5 ) = 0.5888) and seed setting rate (r (y, x 4 ) = 0.2599). But multiple regression analysis revealed the there are dependences among traits. Since there is colinearity between factors, it is hard to accurately estimate the contribution of each trait on the yield. Path analysis could effectively reveal the direct and indirect effect of these traits on the yield. Thus, path analysis was employed to determine the effect of major economic trait on yield.

Path Analysis

Based on stepwise regression analysis, four main factors were selected for path analysis to further clarify the importance of each trait on yield. Results showed that filled kernels per panicle had the largest effect on yield, followed by effective tillers, 1,000 kernel weight, and seed setting rate (Table 3). The path coefficient due to direct effect of effective tillers was 0.06042. However, effective tillers had a large negative effect through filled kernels per panicle, as much as -0.6141, with small positive effect indirectly through seed setting rate and a small negative effect through the 1,000 kernel weight. The ultimate effect of effective tillers on yield was -0.1498. Thus, the number of effective tillers had a large effect on kernels per panicle. In order to increase the number of filled kernels per panicle, it is necessary to reduce the number of effective tillers.

The number of filled kernels per panicle had the largest direct effect on yield with a path coefficient as high as 0.8581. However, its negative indirect effect through the number of effective tillers was also as high as -0.4324 while its positive indirect effect through the weight per thousand kernel and seed setting rate was small, respectively being 0.2214 and 0.0510. Thus, the number of filled kernels per panicle has a combined effect 0.6982. Therefore, in high yield cultivation, it is extremely important to control the negative effect of the number of panicles through the number of filled kernels, focusing on produce large panicles.

The seed setting rate had the smallest direct coefficient, being 0.1503. Its positive indirect effect through the total number of kernels per panicle was higher, being 0.2909. Its indirect effects through the number of effective tillers and the 1,000 kernel weight were weak, being 0.0143 and 0.0720. Ultimately, the seed setting rate had a combine effect of 0.3835.

The 1,000 kernel weight had a direct path coefficient 0.4147. Its negative effects through the number of effective tillers and seed setting rate were respectively -0.2091 and -0.0261, but the positive indirect effect through the number of filled kernels was 0.4583. Therefore, the 1,000 kernel weight had a combined effect of 0.6378.

Discussion and Conclusion

Our analysis showed that the number of filled kernels had the greatest direct positive effect. The effect of seed setting rate and the 1,000 kernel weight were exerted indirectly through the positive effect of the number of filled kernels. The number of effective tillers was strongly negatively correlated with the number of total kernels, the number of filled kernels, and the weight of thousand kernels, and negatively correlated with yield. The total number of kernels per panicle was very significantly correlated with the number of filled kernels. Therefore, the technical strategy to high yield cultivation of Yongyou 8 is to appropriately control effective tillers, focus on large panicle, and pay attention to seed setting rate and kernel weight.

In the current study, results of analysis of five economic traits were somewhat different from previous report, not as previously reported negative correlation [START_REF] Liu | Study on the High Yield and High Efficiency Cultivation Technique of Indica-japonica Intersubspecific Hybrid Rice as Single-cropping Late Rice[END_REF]. Through stepwise regression analysis, the total number of kernels per panicle fell out to be an important factor. Results of path analysis indicated that yield was affected most significantly by the number of filled kernels per panicle. This is consistent with the previous notion -to appropriate control the number of effective panicle, focus on large panicle and pay attention to kernel weight‖. Data for this analysis is more inclusive, containing those collected from more sites, thus is more representative. Yongyou 8 is an intermediate to late season rice line that can be cultivated in areas where rice is grown once or twice a year. This adaptive line has outstanding performance in Zhejiang and Jiangsu. Zhu has pointed out that Yongyou 8 has displayed a dynamic ultrahigh yield trend of -small at first, stable in the middle and strong at last‖. An ultrahigh yield can be achieved by stabilizing the panicle number, focusing on large panicle, increasing the pool, strengthening the source, preventing lodging [START_REF] Zhu | Formation Rule of Super-high-yielding and Cultural Techniques for Yongyou 8[END_REF]. The current study has found that the large panicle character contributed the largest to ultrahigh yield of Yongyou 8. The number of kernels per panicle had a clear larger contribution to yield than other factors. Because the number of effective tillers was very strongly negatively correlated with the total number of kernels and the number of filled kernels, it is imperative to appropriately control the number of effective tillers.

Table 1 ,

 1 Characteristics of Panicle, Kernel and Weight Traits at Various Yield Level

	Yield Level(kg/hm 2 )	>11250	≤11250->9750	≤9750->8250	≤8250
	No of Paddies	21	19	17	8
	Effective Tillers(10 4 /hm 2 )	226.97	237.33	240.61	241.86
	Total Kernels/Panicle	239.36	229.00	201.58	177.50
	Filled Kernels/Panicle	204.85	173.04	155.68	146.64
	Seed Setting Rate	85.39	75.83	77.89	83.00
	1,000 Kernel Weight(g)	28.28	27.98	27.18	26.56
	Mean Yield(kg/hm 2 )	12633.09	10303.72	8953.82	8050.31

Table 2 ,

 2 Correlation Coefficients for Major Economic Traits

		x	2	x	3	x	4	x	5	Actual Yield
	Effective Tillers	-0.7741 **	-0.7156 **	0.0237	-0.3460 **	-0.1498
	Total Kernels/panicle			0.8546 **	-0.1880	0.6262 **	0.4956 **
	Filled Kernels/Panicle					0.3390 **	0.5340 **	0.6982 **
	Seed Setting Rate							-0.1737	0.3835 **
	1,000 Kernel Weight									0.6378 **

Table 3 ,

 3 Path Coefficient Analysis

	Factor	Direct Effect	Through x1	Through x3 Through x4 Through x5
	Effective Tillers	0.6042		-0.6141	0.0036	-0.1435
	Filled Kernels/Panicle	0.8581	-0.4324		0.0510	0.2214
	Seed Setting Rate	0.1503	0.0143	0.2909		-0.0720
	1,000 Kernel Weight	0.4147	-0.2091	0.4583	-0.0261

Strategies for high yield inferred through path analysis of major economical traits in Yongyou 8, a hybrid late season japonica rice


Weiming Liu

Bao Zuda

Mathematical Modeling and Optimization of Schemes of Major Agronomic Factors for Hybrid Rice 'Yongyou 17'

Keywords: hybrid rice, agronomic scheme, mathematical model, model optimization, nitrogen fertilizer, seedling stage, planting density, yield, net profit

In this study, the computer aided modeling and optimization was used to select agronomic scheme of high yield and high net profit for a single-season late rice hybrid 'Yongyou 17'. A field experiment with three agronomic factors including seedling age, planting density, pure N fertilizer rate were performed and the data were used to construct regression models using yield and net profit as the objective variables. Simulation of those models produced the optimum schemes consisting of seedling age, planting density and pure N fertilizer rate to achieve the highest yield and net profit.

Introduction

Nowadays, the optimal experimentaldesignused has beenwidely used inin agriculturalstudied [START_REF] Chen | Effects of water and fertilizer coupling on seedling biomass of Lespedeza bicolor in sandy land of northwestern Liaoning province[J][END_REF][START_REF] Ruan | Effect of different density and N-application rate on maize quality[END_REF][START_REF] Gao | Application of D-saturation optimum design in culture medium screening for micropropagation of virus-free tube potato[END_REF][START_REF] Liu | Multiplication of subculture for Faba Bean (Vicia faba L.) [J][END_REF].Rotational combination design, the secondary saturatedD-optimal design,and simplexlatticedesignwere all applied in the studies for high-yieldcultivation techniquesin rice, sweet potatoes, corn, soybean,turmeric(Dioscorea)by this research group [START_REF] Liu | Study on optimization combination for intercropping spring corn with spring soybean[J][END_REF][START_REF] Liu | Study on fertiilizer application techniques of annual Dioscoreae zingiberensis[J][END_REF][START_REF] Liu | Research on the transplanting density and nitrogen application for Zhongzheyou 1 planted as single cropping late rice[END_REF][START_REF] Liu | Effect of planting density and rates of nitrogen and potassium fertilizers on yield and economic efficiency of sweet potato (cv. Zheshu 13)[END_REF][START_REF] Liu | Optimized test on comprehensive agronomic measures of Zheshu 13[END_REF][START_REF] Liu | Comprehensive agronomic optimization for Zheqiu soybean 3[J][END_REF]

Materials and methods

The field experiment was set up using the seedling age ( 1x , days), planting density ( 2x , cluster/hm 2 ), pure N fertilizer rate ( 3 x , kg/hm 2 ) as the decision variables in a tertiary quadratic general rotary design [START_REF] Yang | Experimental Design and Statistical Analysis[END_REF]. Table 1 contains the three decision variables, the treatment levels and the corresponding linear code. The yield and net profit were used as the objective variables. 1).

Seeds sowed on May 26th, 30th, June 4th, 10th and 14th produced seedlings of 34, 30, 24, 18 and 14 days, respectively, at the transplanting date on June 28 th . Single plant was planted per hole.

The size of each plot was 16 m 2 (3.2m x 5m). To prevent fertilizer from leaching into adjacent plots, the between plot area was covered with plastic film. In each plot, 450 kg/hm 2 of calcium phosphate was applied as base fertilizer. KCl (300 kg/hm

Results and analysis

According to thequadraticregression orthogonalrotation designprinciples, aregressionmathematical model was obtained as follows:

y ˆ= 0 b +   3 1 i i b i x +   j i ij b i x j x +   3 1 i ii b 2 i x +
In which, ( the Chinese words in the following formula should be changed to English word or letters)

剩 剩 f SS /  
Based on the structureof thematrix and theobtained yield andnet output ofYongyou17, quadraticregression analyseswere carried outfocusing on the yield and profit.

Yield result analysis

Using the tertiary quadratic general rotary regression analysis, models for the three agronomic traits and yield were constructed as follows:

(-1.682≤ i x ≤1.682, i =1, 2,3)
In the model, 2

x , 3 x , 2 2

x and 2 3

x are at the extremely significant level, and 1 3 x are at a significant level. The rest of variables are not significant (α=0.25). To increase accuracy level, all the non-significant variables were excluded from the models. Then the regression model was rewritten as:

1 ŷ =12808.65+490.65 2 x +1247.55 3 x -415.2 2 2 x -688.35 2 3 x +333.6 1 x 3 x (-1.682≤ i x ≤1.682, i =1, 2 ,3)
In this regression models, all the variables are significant, and the lack-of-fit sum of squares is not significant. These models should be very effective for simulation experiment.

Test of significance identified yield is mainly affected by planting density and pure N fertilizer rate, and seedling age has a smaller effect. When compared for interaction effect ( 1 x 3

x ), it was found that the lower nitrogen fertilizer rate should be combined with younger seedlings, and higher fertilizer rate with older seedlings. In general, higher N fertilizer rate was correlated with higher yield. When seedling age and nitrogen fertilizer rate were set around treatment level 1, the yield was at the highest level.

Profit result analysis

In this study, rice, seedling, urea, KCl, and calcium phosphate were priced at 2.80 yuan /㎏, 52 yuan/10,000 seedlings, 2.50 yuan/㎏, 3.60 yuan/㎏ and 0.90 yuan/ ㎏, respectively. When estimating the investment, only the cost of materials and supplies were included. After converting the net profit per plot into profit per hectare, a ternary quadratic general rotary regression model was constructed. After exclusion of variables that are not significant (α=0.25), the regression model for net profit is described as follows:

2 ŷ =32565.60+1178.85 2 x +3085.50 3 x -1162.65 2 2 x -1927.35 2 3 x +934.05 1 x 3 x (-1.682≤ i x ≤1.682, i =1, 2 ,3)
Test of significance showed that, same as in the yield model, 2

x , 3 x , 2 2
x , 2 3

x are extremely significant, and 1

x 3

x is at a significant level. The overall regression model is extremely significant. The relationship between the three agronomic factors and net profit is approximately the same as that in the yield model. Planting density and pure N fertilizer rate are the major factors influencing profit, while the effect from seedling age is relatively small. The interaction effect of 1

x 3

x on net profit is also similar to that on yield. Lower nitrogen fertilizer rate combined with younger seedlings, and higher fertilizer rate with older seedlings should yield a higher net profit. The highest net profit was achieved when the seedling age and nitrogen rate were both set at around treatment level 1.

These results indicate that the three agronomic factors have a consistent effect on yield or net profit. Therefore agronomic schemes aimed at increasing yield also can produce high net profit. Then the optimum agronomic schemes for high yield and high profit were selected using computer simulation of the mathematical models described above.

Model optimization

Optimization of the agronomic schemes for high yield

The high yield agronomic schemes were selected using the yield regression model and selecting ( -1.682≤ i x ≤1.682) as the region for constrain [12 , 13].

Twenty-one schemes were found to produce a yield above 12900kg/hm 2 (Table 2). At the 95% confidence level, the value of 1

x , x 2, 3

x each was 0.456-1.206 (26-31 days), and the yield and net profit of single season hybrid rice 'Yongyou 17'. Model analysis indicates that yield and net profit have similar responses to the three agronomic factors. Among the three agronomic factors, planting density and pure N fertilizer rate have bigger effects whereas that from seedling age is relatively small. Production schemes were constructed using simulation of mathematical models integrating data from a field experiment. Under the production conditions used in this study, for 'Yongyou 17'to achieve high yield and high profit, the optimum scheme consists of the following conditions: seedling age of 26-30 days, planting density of 165000-180000 cluster/hm 2 , pure N fertilizer rate of 210-225kg/hm 2 .

Table 1 ,

 1 Decision variables and the treatment levels

	Linear code	1 x	x	2	x	3
	+1	30	202500	201.15
	-1	18	127500	51.15
	0	24	165000	126.15
	+1.682	34	228075	252.30
	-1.682	14	101925	0
	Seedling age was controlled by sowing seeds at different dates (Table

Model optimization for high net profit

The high profit agronomic schemes were selected by simulating net profit regression model within the same constraint region as for yield. Sixteen schemes were found to produce a net profit exceeding 32700 yuan/ hm 2 (Table 3). The value for 1

x , 2

x , 3

x each (at 95% confidence level) was 0.479-1.368(27-32 days), 0.257-0.999 (174638-202463 clusters/ hm 2 ) and 1.058-1.368 (205.5-228.75 kg/ hm 2 ).

It can be seen that the similar values were selected for the three agronomic factors in the net profit schemes and those aimed at high yield.

Conclusions

In this study, regression models were constructed to simulate the relationship between three agronomic factors (seedling age, planting density, pure N fertilizer rate)


Guoqiang Ma

Naixiang Li
email: blinaixiang@tjau.edu.cn

Xiaojuan Wang
email: cwangxiaojuan6996@163.com

Segmentation of Small Animal Computed Tomography Images using Original CT Values and Converted Grayscale Values

Keywords: Image segmentation, Computed tomography, k-means clustering

Medical image segmentation is the foundation of normal and diseased tissue 3D visualization, operation simulation and visual operation. In this paper, we comprehensively use two values that represent the same object (body tissue), to segment by the same algorithm implementation. The original CT images are downloaded from the web source, the dog rib tumor CT scan image by GE medical systems, all the experiment of dataset of 312 thoracic CT scans. The core of the segmentation is k-means clustering algorithm. The segmentation process consisted of two phases: (1) convert CT value to JPG gray value or not use the original CT value as the data sets for clustering; (2) segmentation bone tissue using the new k-means clustering algorithm program which is implemented with MATLAB 2012a programming language and for two-dimensional data matrix directly. The experiment produced strikingly different results. These results may be indicating that not only the segmentation algorithm of CT image is important, but also the data for segmentation is important too.

Introduction

Image segmentation is a challenging task in image analysis. Image segmentation plays important role in medical image. A large variety of methods of have been proposed in several years. Intracranial image segmentation techniques for CT images include semi automated region-growing [START_REF] Kamdarmr | Intracranial volumes in a large series of healthy children[J][END_REF] , intensity threshold [START_REF] Mardini | Intracranial space,brain, and cerebrospinal fluid volume measurements obtained with the aid of three-dimensional computerized tomography in patients with and without Crouzon syndrome[J][END_REF][START_REF] Strik | Three-dimensional reconstruction and volumetry of intracranial haemorrhage and its mass effect[END_REF] , and automated level set methods [START_REF] Jafarian | Automatic fontanel extraction from newborns' CT-images using a model based level set method[END_REF][START_REF] Kazemi | Automatic fontanel extraction from newborns' CT images using variational level set[END_REF] . The process of grouping a set of physical or abstract objects into classes of similar objects is called clustering. A cluster is a collection of data objects that are similar to one another within the same cluster and are dissimilar to the objects in other clusters. Clustering algorithm is a kind of do not need training, unsupervised (unsupervised) statistical method. In a sense, clustering is a classification of self training. Because of the lack of training sample set, so the iterative clustering algorithm is performed on the image classification and extraction of feature values. Among them, K-means [START_REF] Meena | K Means Segmentation of Alzheimers Disease in PET scan datasets: An implementation[END_REF] , Fuzzy C-Means [START_REF] Valvi | Improved fuzzy c-mean algorithm for Medical Image segmentation[END_REF][START_REF] Funmilola | Fuzzy k-c-means Clustering Algorithm for Medical Image Segmentation[END_REF] , EM (Expectation-Maximization) and hierarchical clustering method [START_REF] Lei | Statistical approach to X-Ray CT imaging and its applications in image analysis -part II: A new stochastic model-based image segmentation technique for X-Ray CT image[END_REF][START_REF] Liang | Parameter estimation and tissue segmentation from multispectral MR images[END_REF] is commonly used in clustering algorithms. K-means clustering algorithm calculate the current average of each class, then according to the new generation of mean pixel to reclassify (mean pixel in recent class), the new generation class iterative execution front steps. In this paper, we propose a new processing framework draws on the existing techniques of k-means clustering algorithm and use the CT value for segmentation in CT images. The main algorithm is still using the K-means algorithm. The best contribution of this paper are: (1) completed a specialized MATLAB language program implementation of k-means clustering algorithm for CT images segmentation;

(2) using the CT value of DICOM format files directly to calculate by k-means clustering algorithm for segmentation. The results show the segmentation quality is better than using the transformational grayscale values.

DICOM format and original CT values

DICOM format

DICOM (Digital Imaging and Communications in Medicine) is a standard for handling, storing, printing, and transmitting information in medical imaging. It includes a file format definition and a network communications protocol. A DICOM data object consists of a number of attributes, including items such as name, ID, etc., and also one special attribute containing the image pixel data (i.e. logically, the main object has no "header" as such: merely a list of attributes, including the pixel data). A single DICOM object can have only one attribute containing pixel data. For many modalities, this corresponds to a single image. But note that the attribute may contain multiple "frames", allowing storage of cine loops or other multi-frame data. The same basic format is used for all applications, including network and file usage, but when written to a file, usually a true "header" (containing copies of a few key attributes and details of the application which wrote it) is added.

CT value

The characteristic of CT is able to distinguish the slight difference of human tissue densities; the standard is based on the body tissues of X -ray absorption coefficient (μ value) to decide. In order to calculate and convenience, the linear attenuation coefficient is divided into 2000 units, known as the CT value, unit of measurement called Hounsfield.

Using water as a value of 0, the upper bound of bone CT value is 1000; the lower air CT value is -1000. At present, most of the CT scanner has a range of 1000 or more than 2000. In fact, CT value is the corresponding X-ray value of the attenuation coefficient with the different human tissue in CT image. In the CT image, bone and air CT values were specified for upper and lower bounds. The calculation formula of CT value is shown in equation 1.

1000

- ) ( _   water water tissue HU Value CT    (1)
The CT values of normal human tissues, organs are shown in Table 1.

_ ) ( _    (2)
CT_Value represents the real CT value of body tissue, pixel_value represents the pixel value of a DICOM formatted CT image file, it can be read out by the dicomread() function.

DICOM values covered to Grayscale values

Jpg/Jpeg is an efficient compression format, the file format is JPEG (Joint Photographic Experts Group) standard products. The depth image are 8bit, 12bit, considering the human eye to gray level of sensory discrimination, the general use of the bit depth is 8bit, the 256 gray level. In the gray image in JPG format, each pixel value is the gray value, the range between 0~255. Conversion principle is very simple, is to use the dicomread() function to read the CT value matrix from a DICOM formatted CT image file and converted to a numerical matrix of 0~255. The conversional relation can be represented in equation 3.

) min( ) max( 255)) min( (

CT

Matlab K-means clustering algorithm program

K-means algorithm is a hard clustering algorithm, it is a typical of the objective function clustering method based on the prototype of representative which uses the distance summation of data points to the prototype (class center) as the optimum objective function, using the method for the extremism of function to get the iterative adjustment rules. The K-means algorithm uses the Euclidean distance as the similarity measure, it calculates an optimal classification of vector The 256 gray levels JPG format image is segmented by the new k-means clustering algorithm program. This process is shown in Figure 1.

T k v v v V ) , , , (

Materials and Method

The original CT images are downloaded from the web source, "DICOM7 medical image database", URL: http://www.dicom7.com/. The compression package is named for the "ANONYMIZE2", the dog rib tumor CT scan image by GE medical systems, all the experiment of dataset of 312 thoracic CT scans. The first step is converting the DICOM format images to JPG format images using the equation 1. The second step is segment the DICOM format images using CT values and the JPG images using gray values by the new k-means function.

The configuration of the computer used for experiments is with Intel Core Duo 3.33 GHz CPU, RAM Memory 2 GB, 32 bits Windows 7 ultimate; the algorithms are used MatLab2012a programming. We analysis the quality of segmentation results through observe segmentation result with the naked eye.

Results

The analysis of computation complexity and comparison of implementation speed for the dog rib tumor segmentation shows that our algorithm takes about 1 second. From the visual effect of view by naked eye, the segmentation effect using of CT value is good to using the gray value by the same segmentation algorithm and the same 2. This could be a complex problem.

Table 2. Three examples of segmentation results

DICOM format files

2 1 

 21  based on an initial clustering center, and keep the evaluation value c J to a minimum. Algorithm often uses the error sum of squares function as the clustering criterion function, the squared error function is defined as  of class i C , the parameter p represents a point of class i C . We do not use the k-means function in the MATLAB statistics toolbox and write a new k-means clustering algorithm program specifically for CT image. The process of k-means clustering algorithm is improved by iterative method. The input parameters are two, respectively is the JPG format image file name I_filename and the number of clusters K; output only corresponding class label matrix IDX_mar. The termination condition function program is that the objective function value of error square sum is no change. The actual amounts of code of this process are only 38 lines using MATLAB programming language. The pseudo -code of the improved k-means function as follows.

Fig. 1 .

 1 Fig.1. Segmentation process of JPG image

  Segmented image program implementation. Three examples of segmentation results are shown in Table

  We have presented a segmentation approach for CT images, which was demonstrated and validated using two dataset of the same CT image scans, consisting of 312 image volumes. The core of the segmentation is k-means clustering algorithm. The segmentation process consisted of two phases: (1) convert CT value to JPG gray value or not use the original CT value as the data sets for clustering; (2) segmentation bone tissue using the new k-means clustering algorithm program which is implemented with MATLAB 2012a programming language and for two-dimensional data matrix directly. In this paper, we comprehensively use two values that represent the same object (body tissue), to segment by the same algorithm implementation. The experiment produced strikingly different results. These results may be indicating that the segmentation algorithm of CT image is important and the data for segmentation is important too. In other words, what kind of data can really represent the real thing?

Table 1 .

 1 CT values of normal human tissues

			Classes		CT values(HU)
			Water			0± 10
			cerebrospinal fluid	3~8
			plasma			3~14
			edema			7~7
			alba			25~32
			ectocinerea		30~40
			blood			13~32
			clot			64~84
			liver			50~70
			spleen			50~65
			pancreas		45~55
			kidney			40~50
			muscle			40~80
			gallbladder		10~30
			fat			-20~-80
			calcification		80~300
			atmosphere		<-200
			bone			>+400
	Matlab R2012a provides dicominfo () and dicomread () function, respectively, used to
	read DICOM format image file header data and metadata. For example, in the
	MATLAB	R2012a		command		window	input:
	info=dicominfo('D:\ANONYMIZE2\DICOM7_000000.dcm'), Running result will
	display detailed information about the file. Two parameters of these detailed
	information illustrate the transform relation of the record value and the real CT value,
	they are "RescaleIntercept" and "RescaleSlope". The conversional relation can be
	represented in equation 2.
		CT	Value	HU	pixel	value	Re	scaleSlope
						Re	scaleInter	cept

Acknowledgment

Funds for this research were provided by the Doctor Foundation of Tianjin Agriculture University, No.2012D09.