Mod\`ele g\'eom\'etrique pour une Caract\'erisation morphom\'etrique d'individus dans des conditions \textit{in situ}

TL;DR

This paper introduces a non-invasive geometric model that reconstructs and characterizes individual morphometric features using dual orthogonal images and cylindrical coordinate interpolation.

Contribution

It presents a novel method for morphometric analysis based on modeling individuals as polyarticulated chains with image-based reconstruction.

Findings

Accurate morphometric data can be obtained non-invasively.

The method effectively reconstructs segment borders using cylindrical coordinates.

The approach enables detailed morphometric characterization in situ.

Abstract

We develop a non invasive method for determine morphometric caracteristic of individuals. We modelize a individual by a polyarticulated chain of rigid segments and each of these segments corresponds to a volume. For this, each segment is photographed under two orthogonal plans and each image is digitalized. By assuming that the border of each segment can be described under cylindric coordinates, the border is reconstructed by interpolation. Finally, the principal morphometric data are determined.

Tables

Table 1. Table 1. Pourcentages du volume approché par rapport au volume réel.

$N$
$33$	$108.26$	$100.03$	$99.87$	$99.91$
$17$	$108.02$	$99.81$	$99.64$	$99.69$
$9$	$107.27$	$99.11$	$98.95$	$99.00$
$5$	$105.35$	$97.33$	$97.18$	$97.22$
$3$	$101.58$	$93.86$	$93.71$	$93.75$
$Q$	$4$	$8$	$16$	$32$

Table 2. Table 2. Erreurs commises sur la position du centre de gravité.

$N$
$33$	${\mathrm{1.514\hspace{0.17em}10}}^{-16}$	${\mathrm{4.139\hspace{0.17em}10}}^{-17}$	${\mathrm{6.565\hspace{0.17em}10}}^{-17}$	${\mathrm{1.209\hspace{0.17em}10}}^{-17}$
$17$	${\mathrm{7.751\hspace{0.17em}10}}^{-17}$	${\mathrm{1.123\hspace{0.17em}10}}^{-17}$	${\mathrm{3.395\hspace{0.17em}10}}^{-17}$	${\mathrm{2.76\hspace{0.17em}10}}^{-17}$
$9$	${\mathrm{6.206\hspace{0.17em}10}}^{-17}$	${\mathrm{2.786\hspace{0.17em}10}}^{-17}$	${\mathrm{2.790\hspace{0.17em}10}}^{-17}$	${\mathrm{1.839\hspace{0.17em}10}}^{-17}$
$5$	${\mathrm{8.746\hspace{0.17em}10}}^{-17}$	${\mathrm{4.637\hspace{0.17em}10}}^{-17}$	${\mathrm{2.841\hspace{0.17em}10}}^{-17}$	${\mathrm{2.839\hspace{0.17em}10}}^{-17}$
$3$	${\mathrm{9.256\hspace{0.17em}10}}^{-17}$	${\mathrm{2.762\hspace{0.17em}10}}^{-17}$	${\mathrm{1.508\hspace{0.17em}10}}^{-17}$	${\mathrm{1.499\hspace{0.17em}10}}^{-17}$
$Q$	$4$	$8$	$16$	$32$

Table 3. Table 3. Erreurs commises sur les inerties.

$N$
$33$	$0.0964$	$0.0064$	$0.0050$	$0.0032$
$17$	$0.0881$	$0.0140$	$0.0127$	$0.0109$
$9$	$0.0682$	$0.0378$	$0.0364$	$0.0346$
$5$	$0.0679$	$0.0977$	$0.0963$	$0.0945$
$3$	$0.1390$	$0.2103$	$0.2087$	$0.2069$
$Q$	$4$	$8$	$16$	$32$

Table 4. Table 4. Valeurs de N 𝑁 N , ( c i ) 0 ≤ i ≤ N subscript subscript 𝑐 𝑖 0 𝑖 𝑁 {(c_{i})}_{0\leq i\leq N} , ( α i ) 0 ≤ i ≤ N subscript subscript 𝛼 𝑖 0 𝑖 𝑁 {(\alpha_{i})}_{0\leq i\leq N} et ( β i ) 0 ≤ i ≤ N subscript subscript 𝛽 𝑖 0 𝑖 𝑁 {(\beta_{i})}_{0\leq i\leq N} en fonction de A 𝐴 A , B 𝐵 B , a 𝑎 a et b 𝑏 b .

$A$	$B$	$N$	${\left(c_i\right)}_{0\le i\le N}\in$	${\left({\alpha }_i\right)}_{0\le i\le N}$	${\left({\beta }_i\right)}_{0\le i\le N}$
0	0	0	$\left(1\right)$	$\left(0\right)$	$\left(0\right)$
1	0	1	$\left(a\mathrm{,1}\right)$	$\left(\mathrm{0,1}\right)$	$\left(\mathrm{0,0}\right)$
0	1	1	$\left(b\mathrm{,1}\right)$	$\left(\mathrm{0,0}\right)$	$\left(\mathrm{0,1}\right)$
2	0	2	$\left(a^2\mathrm{,2}a\mathrm{,1}\right)$	$\left(\mathrm{0,1,2}\right)$	$\left(\mathrm{0,0,0}\right)$
0	2	2	$\left(b^2\mathrm{,2}b\mathrm{,1}\right)$	$\left(\mathrm{0,0,0}\right)$	$\left(\mathrm{0,1,2}\right)$
1	1	3	$(ab,b,a\mathrm{,1})$	$\left(\mathrm{0,1,0,1}\right)$	$\left(\mathrm{0,0,1,1}\right)$

Equations

$$(\theta ,z)\mapsto r(\theta ,z).$$

$$\partial\mathcal{D}(z)=\Bigl{\{}M(x,y,z),\quad x=r(\theta,z)\cos\theta,\quad y=r(\theta,z)\sin\theta,\quad\theta\in]-\pi,\pi]\Bigr{\}}.$$

$$\partial\mathcal{D}(z)=\Bigl{\{}M(x,y,z),\quad x=r(\theta,z)\cos\theta+x^{0}(z),\quad y=r(\theta,z)\sin\theta+y^{0}(z),\quad\theta\in]-\pi,\pi]\Bigr{\}}.$$

$$\forall k\in \{1,\dots ,P\},\exists {\theta }_k\in ]-\pi ,\pi ],$$

$$y_k=r({\theta }_k,z_k)\sin {\theta }_k+y^0(z_k).$$

$$\forall i\in \{1,\dots ,N\},\forall j\in \{1,\dots ,n_i\},\exists {\theta }_{i,j}\in ]-\pi ,\pi ],$$

$$y_j=r({\theta }_{i,j},z_i)\sin {\theta }_{i,j}+y^0(z_i).$$

$$\forall j\in \{1,\dots ,n_i\},\exists {\theta }_j\in ]-\pi ,\pi ],$$

$$y_j=r_i({\theta }_j)\sin {\theta }_j+y_i^0.$$

$$\forall \theta \in ]-\pi ,\pi ],$$

$$x^0(z_i)=x_i^0,$$

$$y^0(z_i)=y_i^0.$$

$$x_j=r_{i,j}\cos {\theta }_j+x_i^0,$$

$$y_j=r_{i,j}\sin {\theta }_j+y_i^0.$$

$$\forall l\in \{1,\dots ,Q\},\rho ({\theta }_l)=R_l.$$

$$R_0=R_Q,$$

$$\forall l\in \{0,\dots ,Q\},\rho ({\theta }_l)=R_l.$$

$$\forall l\in \{0,\dots ,Q\},{\rho }^{\prime }({\theta }_l)=0.$$

$$S={\left((x_k-X_k{)}^2+(y_k-Y_k{)}^2-R^2\right)}^2.$$

$$S^{\prime }=\sqrt{(x_k-X_k{)}^2+(y_k-Y_k{)}^2}$$

$$x_k=R_k\cos {\theta }_k+x^0,$$

$$y_k=R_k\sin {\theta }_k+y^0.$$

$$\forall k\in \{1,\dots ,P\},\rho ({\theta }_k,z_k)=R_k.$$

$$\hat{I}={\iiint }_{\mathcal{D}}f(x,y,z)dm,$$

$$I={\iiint }_{\mathcal{D}}f(x,y,z)dV,$$

$$dm=\mu dV.$$

$$\hat{I}=\mu I.$$

$$\{\begin{array}{ll}& x=r\cos \theta +x^0(z),\\ & y=r\sin \theta +y^0(z).\end{array}$$

$$J=\det \left(\begin{array}{ccc}\frac{\partial x}{\partial r}& \frac{\partial x}{\partial \theta }& \frac{\partial x}{\partial z}\\ \frac{\partial y}{\partial r}& \frac{\partial y}{\partial \theta }& \frac{\partial y}{\partial z}\\ \frac{\partial z}{\partial r}& \frac{\partial z}{\partial \theta }& \frac{\partial z}{\partial z}\end{array}\right).$$

$$J=\det \left(\begin{array}{ccc}\cos \theta & -r\sin \theta & {x^0}^{\prime }\\ \sin \theta & r\cos \theta & {y^0}^{\prime }\\ 0& 0& 1\end{array}\right)=r{\cos }^2\theta +r{\sin }^2\theta =r,$$

$$dV=dxdydz=rdrd\theta dz.$$

$$I$$

$$={\int }_{z=z_{\min }}^{z_{\max }}{\iint }_{\begin{array}{c}\theta \in [-\pi ,\pi ]\\ r\in [0,r(\theta ,z)]\end{array}}f(r\cos \theta +x^0(z),r\sin \theta +y^0(z),z)rdrd\theta dz,$$

$$I={\int }_{z_{\min }}^{z_{\max }}F(z)dz,$$

$$F(z)={\iint }_{\begin{array}{c}\theta \in [-\pi ,\pi ]\\ r\in [0,r(\theta ,z)]\end{array}}f(r\cos \theta +x^0(z),r\sin \theta +y^0(z),z)rdrd\theta .$$

$$f(x,y,z)={\left(x-\hat{x}\right)}^A{\left(y-\hat{y}\right)}^B{\left(z-\hat{z}\right)}^C$$

Full text

1. Introduction

L’une des problématiques essentielles en biomécanique humaine ou animale est de déterminer les caractéristiques morphométriques des individus en mouvement. Ces caractéristiques sont la position des centres de masse des segments de l’individu modélisé sous la forme d’un système poly-articulé, les masses de ces segments et leurs matrices d’inertie. La connaissance de ces données est essentielle dans le sens où elle nous permettent de calculer à chaque instant la position du centre de gravité de l’individu ainsi que par dynamique inverse ou dynamique directe les forces et moments interarticulaires au cours du mouvement. Il existe dans la littérature des tables anthropométriques [ZS83, Win09] ou morphométriques [WD87, LHB12]. Cependant, ces tables concernent des individus moyens et ne peuvent donc rendre compte des variations morphométriques des individus en fonction de leur environnement ou encore au cours de leur croissance.

L’idée de ce papier est de reconstruire la frontière de chacun des segments de l’individu modélisé en utilisant les coordonnés cylindriques. Fondamentalement, ce choix est guidé par le fait que chaque membre est proche d’un volume de révolution autour d’un axe et donc, qu’à altitude $z=z_{0}$ constante, la frontière de ce membre est proche d’un cercle, donc définie en coordonnées polaires par la fonction $r$ : $\theta\mapsto r(\theta)$, où $r$ ne varie que << peu >>.

Une autre idée consisterait à utiliser les fonctions $B$-splines comme dans [HNA07] (on pourra aussi consulter les ouvrages [HM13, Boo01]), où la frontière recherchée est définie par une surface $B$-spline. Cette méthode n’a pas été envisagée essentiellement pour deux raisons :

• —

Comme écrit dans [HNA07] : <<Once the initial bounding volume is constructed, a set of spatial points is generated uniformly over the volume [NF97]. For each of these spatial points, the closest data point is identified and the $B$-spline solid is deformed to make the data point and corresponding spatial point coincident. In general, this method produces good initial shapes when the number of control points are chosen to be fewer than the number of data points. Subsequent manual adjustments of the control points can be made for shape refinement.>> L’absence de gestion automatique des points de contrôle qui définissent les splines impose une intervention humaine et fastidieuse.

• —

Une spline est définie par les points de contrôle et non par les points de passages, qui sont les seuls expérimentalement connus. Un système linéaire permettrait de passer de l’un à l’autre et donc de définir la spline grâce à la connaissance des points de passage. L’inversion de ce système linéaire risque cependant d’être source de longs calculs et ce d’autant plus que les points expérimentaux de passage risquent d’être en nombre important.

Dans [Hen99], Henderson utilise une technique dont la nôtre s’inspire, mais qui n’est pas satisfaisante en l’état. Elle a été adaptée et améliorée. À partir de deux vues, l’une latérale et l’autre frontale (voir figure 1, extraite de [Hen99]), il détermine 4 points et leurs coordonnées, par mesure de l’intersection de segments avec le profil étudié. Il en déduit donc ensuite une ellipse passant par ces quatre points et, grâce à l’ensemble des ces ellipses, il reconstitue le profil 3D, en en déduisant le calcul du centre de masse et du volume, en supposant la masse volumique constante. Cette méthode semble être apparemment utile pour reconstruire notre profil et en déduire les données morphométriques et servira de point de départ ; cependant, on ne la conservera pas telle quelle pour plusieurs raisons :

• —

Cette méthode ne fonctionne que si l’animal étudié présente un plan sagittal de symétrie, ce qui ne sera pas toujours le cas, notamment si des ondulations sont présentes, comme par exemple la queue de la figure 9 ;

• —

Seuls quatre points sont considérés et cette méthode ne permet pas d’en rajouter, ce qui pourra être fait pour celle présentée, pour prendre en compte certains détails ;

• —

Enfin, dans cette méthode, seuls la position du centre et les volumes sont déterminés ; il manque les inerties, ce que l’on propose de déterminer.

Nous reviendrons néanmoins sur cette méthode que l’on comparera à celle proposée.

Enfin, une dernière idée aurait pu être utilisée, à l’instar des photos aériennes prises par l’IGN, la stéréoscopie (voir http://fr.wikipedia.org/wiki/Stéréoscopie). Cette technique nécessite de faire plusieurs paires de photos, avec de légers décalages, afin de couvrir l’ensemble de la surface du membre étudié et utilise ensuite des techniques de calculs plus sophistiquées. Elle n’est pas pertinente ici, puisque l’on essaye de proposer une méthodologie recquérant peu de prises de vues et de mesures et des moyens de calcul relativement sommaires. De plus, cette méthode ne permet pas d’exploiter le fait que, à $z$ fixé, la frontière est proche d’un cercle.

Nous présenterons en section 2 l’acquisition des données (relevé géométriques des points de la frontière), puis, en section 3, comment cette frontière est reconstruite, et enfin, en section 4, le calcul des principales données morphométriques (volume, position du centre de masse et matrices d’inerties). Puisque la masse volumique est supposée constantes, ce sont en fait des données géométriques qui sont calculées. Enfin, quelques simulations seront présentées en section 5.

2. Acquisition des données

On se donne deux photos de dos et de profil du membre considéré (ici la tête). Voir figure 2.

Chaque segment est digitalisé selon un nombre de coupes $n$ définies par l’utilisateur. Il en résulte $N=n+1$ plans de coupe. Le plan de coupe pour la vue de profil correspond au plan $(y,z)$. Le plan de coupe pour la vue dorsale correspond au plan $(x,z)$. L’axe $z$ correspond à l’axe longitudinal du segment (voir l’image 3). Les valeurs de $z$ sont donc celles des distances entre l’origine (fixée au premier plan de coupe à l’extrémité distale ou caudale du segment) et le plan de coupe considéré. Après digitalisation, toutes les coordonnées sont normalisées par la distance entre les plans de coupe des extrémités proximales et distales du segment.

Dans les cas présentés par la suite, seuls quatre points symétriques ont été considéres comme dans [Hen99], de telle sorte que l’on aurait pu utiliser les ellipses pour reconstituer notre volume. En fait, nous procéderons autrement, ce qui permettra un ajout possible de points. Nous montrerons néanmoins que nos résultats sont très proches de ceux obtenus par [Hen99] en section 6.

Pour chaque coupe d’altitude $(z_{i})_{1\leq i\leq N}$, on connaît donc les coordonnées $(x_{j},y_{j},z_{i})_{(1\leq j\leq n_{i})}$ de points de la frontière. On connaît donc un nombre $P=\sum_{i}n_{i}$ de points de la frontière du membre considéré, notés aussi $(x_{k},y_{k},z_{k})_{1\leq k\leq P}$.

Le tronc de l’image 3 sera utilisé tout au long du papier ou du transparents comme fil conducteur qui illustre la méthode.

Remarque 2.1*.*

Cette façon de procéder permet aussi d’enrichir les données en rajoutant des points caractéristiques, visibles sur les deux clichés.

Remarque 2.2*.*

Chaque section sera mesurée soit par rapport à un centre correspondant à une droite fixée, soit par rapport au barycentre des points expérimentaux mesurés.

Remarque 2.3*.*

Les différentes longueurs obtenues sont telles que l’altitude mimimale est nulle et l’altitude maximale est égale à 1. Il faudra donc multiplier toutes les longeurs obtenues par un facteur d’échelle $\zeta$, égal par exemple (pour le tronc de la figure 3) à $0.063075$ pour obtenir les longueurs en mètres. Ainsi, les différentes données morphométriques présentées en section 4 doivent être elles aussi multipliées respectivement par $\zeta^{3}$ pour obtenir le volume réel, $\zeta$ pour les coordonnées du centre de masse, et par $\zeta^{5}$ pour la matrice d’inertie.

3. Reconstruction de la frontière

On se donne un volume $\mathcal{D}$ de l’espace, dont on suppose la frontière $\partial\mathcal{D}$ connue sous la forme de la surface paramétrée donnée en coordonnées cylindrique par la fonction $r$ de $]-\pi,\pi]\times[z_{\min},z_{\max}]$ dans $\mathbb{R}$ (voir figure 4) :

[TABLE]

On note $P_{z}$ le plan d’altitude $z$. Ainsi, pour $z\in[z_{\min},z_{\max}]$ fixé, la frontière $\partial\mathcal{D}(z)=\partial\mathcal{D}\cap P_{z}$ est décrite en coordonnées polaires par

[TABLE]

Remarque 3.1*.*

Conformément à la remarque 2.2, on peut autoriser un centre variable de coordonnées $x^{0}(z)$ et $y^{0}(z)$ (dépendant de $z$) pour le repérage des coordonnées par section (voir figure 5) : on remplacera donc le paramétrage polaire (3.2) par

[TABLE]

Les fonctions $x^{0}$ et $y^{0}$ seront définies (constantes) par l’utilisateur ou mesurées expérimentalement.

On cherche donc une fonction $r$ et deux fonctions $x^{0}$ et $y^{0}$ vérifiant

[TABLE]

Ce problème est équivalent à

[TABLE]

Nous allons proposer plusieurs méthodes et conserver celle qui semblera la mieux adaptée au problème donné.

3.1. Reconstruction de la frontière par coupe

Le problème (3.5) peut se découpler ainsi : pour chaque $i\in\{1,\ldots,N\}$, on cherche $x^{0}_{i}$, $y^{0}_{i}$ et une fonction $r_{i}$ de $]-\pi,\pi]$ dans $\mathbb{R}$ vérifiant

[TABLE]

Les fonctions $r$, $x^{0}$ et $y^{0}$ seront alors connues uniquement aux points d’altitude $z_{i}$ avec

[TABLE]

On détermine alors pour chaque $i\in\{1,\ldots,N\}$ et pour chaque $j\in\{1,\ldots,n_{i}\}$, $r_{i,j}$ et $\theta_{j}$ tels que

[TABLE]

Bref, on cherche donc à trouver une fonction $\rho$ de $]-\pi,\pi]$ dans $\mathbb{R}$ telle que, pour $Q\in\mathbb{N}^{*}$ donné, ${(\theta_{l})}_{1\leq l\leq Q}\in]-\pi,\pi]^{Q}$ et ${(R_{l})}_{1\leq l\leq Q}\in\mathbb{R_{\,+}}^{Q}$ donnés, on ait

[TABLE]

Sans perte de généralité, on supposera les $\theta_{l}$ strictement croissants. La fonction $\rho$ est naturellement $2\pi$-périodique : on peut donc rajouter l’angle $\theta_{0}=\theta_{Q}-2\pi$ et la valeur $R_{0}=R_{Q}$ et remplacer (3.8) par

[TABLE]

Pour éviter la présence de points anguleux, on cherchera une fonction $\rho$ dérivable. On est donc exactement dans le cas de l’interpolation périodique. Voir l’annexe A et [Bas14]. On aura donc une fonction polynomiale cubique par morceaux.

Deux solutions d’offrent à nous :

• —

Soit, on cherche une fonction $\rho$ de classe $C^{1}$ uniquement, mais qui respecte la monotonie des données, ce que permet la fonction perpchip. Cette méthode sera notée la méthode A. La fonction $\rho(\theta)$ ne variant que << peu >>, cette méthode sera à privilégier puisqu’elle limite les variations d’ampltitude de la fonction $\rho$.

• —

Soit, on cherche une fonction $\rho$ de classe $C^{2}$, ce que permet la fonction perspline. Cette méthode sera notée la méthode B.

Une troisième possibilité est aussi envisagée : on considère que les points mesurés correspondent à des extrémaux locaux de $\rho$, ce qui est vérifié quand les points extrémaux sont considérés pour chaque profil : voir figure 6. On a donc la condition supplémentaire

[TABLE]

On cherchera donc l’unique fonction $\rho$ polynomiale cubique par morceaux sur chaque intervalle $[\theta_{l},\theta_{l+1}]$, de classe $C^{1}$ et vérifiant (3.9) et (3.10). Cette méthode sera notée la méthode C.

Remarque 3.2*.*

Une autre idée serait d’utiliser le fait qu’à $z$ fixé, la définition en polaire de la frontière $\mathcal{D}(z)$ étant définie par la fonction $r$ : $\theta\mapsto r(\theta)$, où $r$ ne varie que << peu >>. On pourrait alors approcher cette frontière par un cercle, grâce à une méthode de moindre carré, comme on a fait dans [Cre+12] :

<< The characteristics of the circle can be determined using a least squares method. The center $\Omega$ and radius $R$ of the circle minimizing the distance sum of squares between experimental data $(x_{k},y_{k})$ and theoretical data $(X_{k},Y_{k})$ were researched. Firstly, a direct method was used to minimize the sum of

[TABLE]

Secondly, the sum

[TABLE]

was minimized by using an iterative method. For this method, the results of the direct optimization were used as initial values for $\Omega$ and $R$. This final optimization was performed with the library Matlab least squares geometric element software, available at http://www.eurometros.org/gen_report.php?category=distributions&pkey=14. For more details, see [Raz97, Raz98, ARW01, EMM07, MK91, Fan90, CJ89]. >>

Dans ce cas, $x^{0}_{i}$ et $y^{0}_{i}$ seront déterminés, ainsi que la fonction constante $\rho=\rho_{i}$. Cette méthode n’a pas été implémentée car elle ne permet pas que la frontière passe exactement par chacun des points expérimentaux.

3.2. Reconstruction de la frontière par interpolation à deux dimensions

On peut aussi considérer de nouveau le problème (3.4). Dans le cas où $x^{0}$ et $y^{0}$ sont constants, on détermine, pour chaque $k\in\{1,\ldots,P\}$, $R_{k}$ et $\theta_{k}$ tels que

[TABLE]

On doit donc déterminer cette fois une fonction $\rho$ de $]-\pi,\pi]\times[z_{\min},z_{\max}]$ dans $\mathbb{R}$ telle que,

[TABLE]

On utilisera alors les deux fonctions interp2 et griddata de matlab qui permettent de déterminer une telle fonction $\rho$, sur une grille aussi fine que l’on veut et donc de calculer les valeurs de la frontière de $\mathcal{D}$. Cette méthode sera notée la méthode D.

Cette méthode permet de déterminer des points de la frontière à toute altitude $z$, contrairement aux méthodes A, B ou C, qui permet de déterminer des points de la frontière uniquement aux altitudes $z_{i}$. Mais, elle est plus longue à calculer. De plus, la périodicité en $\theta$ n’est plus assurée, ce qui fait apparaître des points anguleux sur la frontière reconstituée.

4. Calculs des données morphométriques

Nous déterminons dans cette section les principales données morphométriques (volume, position du centre de masse et matrices d’inerties), grâce à la détermination de $\mathcal{D}$ et en supposant la masse volumique $\mu$ de chacun des segments connue et homogène.

4.1. Généralités et calcul en coordonnées cylindriques

Soit une fonction $f$ est définie (et intégrable) de $\mathcal{D}$ dans $\mathbb{R}$. On considère les intégrales

[TABLE]

où $dm$ est l’élément de masse élémentaire autour du point de coordonnées $(x,y,z)$ et $dV$ le volume élémentaire autour de ce point, liés par

[TABLE]

Si $\mu$ est connue et homogène, on a donc

[TABLE]

Déterminons maintenant l’intégrale $I$. On considère le changement de variable correspondant à (3.3), c’est-à-dire, l’application qui à $(x,y,z)$ associe $(r,\theta,z)$ tels que, pour $z\in[z_{\min},z_{\max}]$ fixé,

[TABLE]

Les nouvelles variables sont $(r,\theta,z)$ et le jacobien $J=\frac{D(x,y,z)}{D(r,\theta,z)}$ est donné par

[TABLE]

On a alors

[TABLE]

et on retrouve donc le fait que, comme en coordonnées cylindriques,

[TABLE]

On intègre donc, finalement, par tranches :

[TABLE]

Bref,

[TABLE]

Enfin, on suppose que l’on se place dans le cas particulier où $A$, $B$ et $C$ sont trois entiers naturels, que $\widehat{x}$, $\widehat{y}$ et $\widehat{z}$ sont trois réels et que et que

[TABLE]

On note

[TABLE]

et

[TABLE]

4.1.1. Calcul approché de l’intégrale fournissant $F$

À $z$ fixé, posons

[TABLE]

On considère le polynôme à deux variables $P$ défini par

[TABLE]

En le développant, on constate qu’il existe $N\in\mathbb{N}^{*}$, ${(c_{i})}_{0\leq i\leq N}\in\mathbb{R}^{N+1}$ (qui dépendent de $z$), ${(\alpha_{i})}_{0\leq i\leq N}\in\mathbb{R}^{N+1}$ et ${(\beta_{i})}_{0\leq i\leq N}\in\mathbb{R}^{N+1}$ (qui ne dépendent que de $A$ et de $B$), tels que

[TABLE]

Voir l’annexe B où les calculs de $N\in\mathbb{N}^{*}$, ${(c_{i})}_{0\leq i\leq N}$, ${(\alpha_{i})}_{0\leq i\leq N}$ et ${(\beta_{i})}_{0\leq i\leq N}$ en fonction des valeurs particulières de $A$ et de $B$ utiles pour la suite sont donnés. Ainsi, pour $z$ fixé, on a donc

[TABLE]

et donc

[TABLE]

Ainsi, il vient

[TABLE]

Finalement, on a

[TABLE]

D’après la reconstruction de la frontière faite en section 3, on peut donc considérer un entier $M\in\mathbb{N}^{*}$ et ${(z_{k})}_{1\leq k\leq M}$ (qui ne sont pas nécessairement les points expérimentaux) vérifiant

[TABLE]

On considère les pas ${(h_{k})}_{1\leq k\leq M}$ définis par

[TABLE]

On approche alors l’intégrale $I(A,B,C)$ de la façon suivante

[TABLE]

et donc

[TABLE]

et $c_{i}(z_{k})$ est défini par le tableau B (où $a$ et $b$ sont définis en fonction de $z_{k}$ grâce à (4.9)). Cette intégrale sera approchée par exemple par la méthode des trapèzes (voir par exemple [BM03]).

4.1.2. Calcul exact de l’intégrale fournissant $F$

La formule (4.7)-(4.8) peut se réécrire, compte tenu de (4.9), pour chaque valeur de $z$ :

[TABLE]

La détermination de la fonction $r(\theta,z)$ (à $z$ fixé, pour les méthodes A à C) se fait par l’intermédiaire de fonctions qui fournissent un polynôme de degré 3 par morceaux. On peut donc supposer que, sur chacun des intervalles considérés du type $[\theta_{0},\theta_{1}]$, on a

[TABLE]

Il faut donc déterminer, pour chacun de ces intervalles,

[TABLE]

Dans l’intégrale en $\theta$, on fait le changement de variable $\theta^{\prime}=\theta-\theta_{0}$. Il vient donc

[TABLE]

On est capable de déterminer explicitement cette fonction (en utilisant par exemple matlab symbolique). On a donc

[TABLE]

On approchera finalement l’intégrale $I(A,B,C)$ comme précédemment :

[TABLE]

4.1.3. Calcul exact de l’intégrale fournissant $F$ dans le cadre des moindres carrés

Si on utilise la méthode de la remarque 3.2, le calcul est plus simple. On utilise de nouveau les formules (4.15). Pour le calcul de $I_{i,k}$, on utilise le fait que $r(.,z_{k})$ est constant et égal à $\rho_{k}$. On a donc

[TABLE]

cette intégrale pouvant se calculer explicitement (en utilisant matlab symbolique par exemple).

4.2. Applications

On choisit tout d’abord

[TABLE]

Le volume $V$ de $\mathcal{V}$ est donné par

[TABLE]

et donc donc finalement

[TABLE]

Par exemple, l’abcisse $x_{G}$ du centre de masse de $\mathcal{V}$, de masse $\mathcal{M}$, est donnée par

[TABLE]

et donc, de façon plus générale, pour la position du centre de masse, on a

[TABLE]

On choisit maintenant

[TABLE]

La matrice d’inertie est donnée par

[TABLE]

et il vient finalement

[TABLE]

Souvent, $\mu$ est de l’ordre de 1000. On normalisera en choisissant $\mu=1$ dans les simulations présentées plus loin.

Remarque 4.1*.*

Les inerties obtenues devront donc être multipliée par la masse volumique réelle. Voir aussi la remarque 2.3.

5. Simulations

5.1. Validation : comparaison avec un ellipsoïde

Comme dans [Hen99, p. 96], nous avons comparé nos résultats avec les résultats exacts d’un ellipsoïde, d’équation

[TABLE]

où les trois demi-axes sont donnés. En plus de la comparaison avec le volume analogue à celle de [Hen99], nous donnons aussi la comparaison avec le centre de masse et les inerties. On peut montrer que le volume de cet ellipsoïde est égal à

[TABLE]

tandis que le centre de gravité est nul et que

[TABLE]

Comme dans [Hen99], nous avons fait varier $N-1\in\{2,4,8,16,32\}$ le nombre total de coupes et $Q\in\{4,8,16,32\}$, le nombre total de valeurs d’angles choisis. Nous avons aussi choisi $a=1.0$, $b=1.6$ et $c=2.4$.

Voir la figure 7 où quelques volumes ont été tracés. On voir sur cette figure que les volumes semblent dépendre surtout de $N$.

Voir aussi le tableau 1 où ont été donnés les pourcentages du volume approché par rapport au volume réel. On constate sur ce tableau que les pourcentages se rapprochent de 100, comme dans le tableau 1 de [Hen99]. La différence est que les résulats sont plus précis et que, contrairement à Henderson, les pourcentages sont parfois supérieurs à 100.

On a aussi affiché dans les tableaux 2 et 3, les erreurs commises sur les évaluations du centre de gravité, très faibles, par raison de symétrie mais aussi sur les inerties qui se rapprochent de zéro quand $N$ et $Q$ augmentent.

Enfin, pour $N=1000$ et $Q=1000$, on obtient les erreurs suivantes sur le volume, le centre de gravité et l’inertie :

[TABLE]

ce qui est très faible.

5.2. Exemple du Pogona

En choisissant un centre de section fixe et nul (voir remarque 2.2), on a pu utiliser les 4 méthodes. Les calculs des données morphologiques (volume, position du centre de masse et matrices d’inerties données par (4.22),(4.23) et (4.25), notés avec des indices X pour X dans $\{\text{A},\text{B},\text{C},\text{D}\}$ en fonction des quatre méthodes présentées) sont présentées ci-dessous, pour les différentes méthodes utilisées. On obtient pour le volume

[TABLE]

pour la position du centre de masse

[TABLE]

et pour les matrices d’inerties

[TABLE]

On observe donc des résultats très proches pour les quatre méthodes.

On présente en figure 8 les résultats des quatre méthodes appliquées au Pogona de la figure 2.

En choisissant maintenant, un centre de section égal au centre de masse des point expérimentaux (voir remarque 2.2), on a pu utiliser les 3 méthodes.

Les calculs des données morphologiques sont présentées ci-dessous, pour les trois méthodes utilisées. On obtient pour le volume

[TABLE]

pour la position du centre de masse

[TABLE]

et pour les matrices d’inerties

[TABLE]

On observe donc des résultats très proches pour les trois méthodes.

5.3. Comparaison de quelques méthodes

Conformément à ce que l’on a présenté en remarque 2.2 et en section 4.1, on peut faire le calcul avec plusieurs possibilités :

• —

en utilisant comme centre de section, le barycentre des points mesurés ou un point fixe (choisi de coordonnées nulles) ;

• —

en calculant l’intégrale en $\theta$ de façon exacte ou approchée.

En choisissant un calcul approché de l’intégrale en $\theta$, on a comparé le calcul avec calcul avec barycentre (résultats obtenus dans (5.7), (5.8) et (5.9)) et sans barycentre. L’écart obtenu entre les deux méthodes vaut

• $\bullet$

pour le volume $5.527\,10^{-4}$ ;

• $\bullet$

pour la position du centre de gravité $5.985\,10^{-3}$ ;

• $\bullet$

pour la matrice d’inertie $1.349\,10^{-4}$.

Ces écarts étant faibles, on utilisera donc par la suite la méthode avec barycentre, ce qui autorise aussi un éventule décalage des sections par rapport à un axe de référence.

En choisissant le calcul avec calcul avec barycentre, on a aussi comparé le calcul approché de l’intégrale en $\theta$ avec $1001$ points (résultats obtenus dans (5.7), (5.8) et (5.9)) avec le calcul exact de cette intégrale. On obtient un écart total (pour le volume, centre de gravité et inerties) égal à $7.197\,10^{-11}$ et un rapport de temps d’exécution (exact/approché), égal à $1.652867$. De même, si on choisit le calcul avec origine nulle de section, on obtient un écart total égal à $1.286\,10^{-12}$ et un rapport de temps d’exécution égal à $5.230257$. On conservera donc pour la suite le calcul approché de l’intégrale en $\theta$ avec $1001$ points.

5.4. Exemple du lézard

On présente en figure 10 les résultats de la méthode conservée (c’est-à-dire la A) pour le lézard (Anolis ferreus : queue) dont le membre étudié est représenté en figure 9.

On obtient pour le volume

[TABLE]

On présente en figure 12 les résultats de la méthode conservée (c’est-à-dire la A) pour le lézard (Anolis ferreus : tronc) dont le membre étudié est représenté en figure 11.

On obtient pour le volume

[TABLE]

5.5. Comparaison avec le cylindre

Remarque 5.1*.*

Pour le tronc du lézard, proche d’un cylindre, on a aussi calculé les données morphométriques en considérant un cylindre moyen. Ce cylindre est de hauteur $h$ (sur l’axe $z$) et de rayon $R$. On vérifie que

[TABLE]

et donc

[TABLE]

De même,

[TABLE]

Enfin,

[TABLE]

et donc

[TABLE]

Par symétrie (ou imparité) les termes non diagonaux de la matrice d’inertie sont nuls et on a

[TABLE]

et donc

[TABLE]

Les résultats sont du même ordre que ceux obtenus précédemment (voir (5.11)). On obtient pour le volume

[TABLE]

pour la position du centre de masse

[TABLE]

et pour les matrices d’inerties

[TABLE]

Sur la figure 13, a été tracée la frontière $\partial\mathcal{D}(z_{0})$ correspondant à $z_{0}=0.2667$. On constate sur cette figure que les trois méthodes donnent des résultats très proches.

Pour comparer aussi avec la méthode D, on a tracé en 14, a été tracée la frontière $\partial\mathcal{D}(z_{0})$ correspondant à $z_{0}=0.2667$ avec centre de section nul. On constate sur cette figure que les quatre méthodes donnent des résultats très proches. Comme annoncé, la frontière correspondant à la méthode D présente un point anguleux. Cette image conforte l’idée que chaque section est proche d’un cercle et, qu’en coordonnées polaires, $r(\theta)$, ne varie que peu en $\theta$.

Sur la figure 15, a été tracée la surface $(\theta,z)\mapsto\rho(\theta,z)$ de la méthode D (avec centre de section nul) et les points expérimentaux.

6. Comparaison avec les ellipses de Henderson

Puisque notre méthode correspond à la méthode de Henderson quand $Q=4$, montrons que les résultats que l’on obtient sont proches de ceux qu’auraient obtenus Henderson.

Pour cela, il suffit de remarquer que l’utilisation de la fonction perpchip avec 4 points symétriques donne une fonction $r(\theta)$ dont les dérivées sont nulles aux valeurs de $\theta$ multiples de $\pi/2$. Voir la figure 16(b).

Pour comparer les deux solutions obtenus, il suffit de calculer les deux fonctions $r_{1}$ et $r_{2}$ correspondant aux coordonnées polaires respectivement donnés pour la solution avec perpchip et avec l’ellispe.

On peut par symétrie, supposer que $\theta$ décrit $[0,\pi/2]$. La fonction $r_{1}$ donnée par perpchip et donc l’unique fonction $f$ de degré 3 sur $[0,\pi/2]$ telle que

[TABLE]

Grâce à la théorie de l’interpolation (voir par exemple [BM03]), on montre que nécessairement, on a

[TABLE]

L’ellipse de demi-axes $a$ et $b$ est d’équation

[TABLE]

Si on pose $x=r_{2}\cos\theta$ et $y=r_{2}\sin\theta$, on a donc

[TABLE]

Pour caluler l’écart relatif entre les deux courbes, il suffit de calculer l’écart maximun entre $r_{1}$ et $r_{2}$. puisque $r_{1}(0)=r_{2}(0)$ et $r_{1}(\pi/2)=r_{2}(\pi/2)$, Il suffit de donc de résoudre $r^{\prime}_{1}(\theta)-r^{\prime}_{2}(\theta)=0$ sur $]0,\pi/2[$, soit encore l’équation

[TABLE]

que l’on résoud par un solveur et on en déduit l’écart maximum.

On peut aussi comparer l’aire entre les deux courbes. Pour l’ellipse on a, en se ramenant à un quart de l’aire totale

[TABLE]

Pour l’ellipse donnée par (6.1), on a, en passant en coordonnées polaires

[TABLE]

On obtient donc

[TABLE]

On détermine cette intégrale en symbolique par exemple et on obtient

[TABLE]

On n’a plus qu’à calculer

[TABLE]

On vérifie que si $a=b$, les deux erreurs sont nulles.

Pour $a=1.2$ et $b=1.0$, on a obtenu un écart relatif en $r$, $\varepsilon_{r}=1.27\%$ et en surface égal à $\Delta S=1.24\%$, soit donc petits et il y aura donc peu d’écart entre les résultats obtenus avc notre méthode et les ellipses de Henderson. Voir aussi la figure 17.

On a aussi fait varier $a$ dans $[a_{\min},a_{\max}]$ et $b$ dans $[b_{\min},b_{\max}]$ avec respectivement $n_{a}$ et $n_{b}$ valeurs. Pour les paramètres suivants,

[TABLE]

l’écart maximum obtenu en $r$ vaut $6.17\%$ et en surface égal à $6.19\%$, soit donc petits.

Voir aussi la figure 18.

Profitons-en aussi pour tracer des ellipses en faisant varier $b$ : on choisit $a=1.0$ et $b$ prenant ses valeurs dans $\{1.3,1.5,2.0,4.0\}$. Voir la figure 19.

Dans certains cas, les ellipses de Henderson paraissent donner des meilleurs résultats, notamment quand les ellipses sont trop dissymétrique (voir par exemple figures 19(c) et 19(d)), où les profils ne semblent pas convexes. Cela n’est pas contraire à la monotonie de $r$. On peut établir par un critère si des profils ne sont pas convexes. En effet (on renvoie à [GAJ94, p. 96] ou [Bas07, p. 24]) : une fonction définie en polaire par $r(\theta)$ est convexe et l’origine appartient à la convavité en au point défini par l’angle $\theta$ si et seulement si

[TABLE]

ce que l’ont sait établir si $r$ est un polynôme de degré 3.

Par exemple, aucun des profils des figures 7 n’est convexe, ce qui est surtout visible pour les petites valeurs de $N$.

Pour la figure 10, seul le profil d’indice $5$ n’est pas convexe.

Pour la figure 12, seul le profil d’indice $16$ n’est pas convexe.

En l’absence de convexité, il suffira de rajouter des points de mesure comme le permet notre méthode.

7. Conclusion

Comme observé dans la section 5, les 4 méthodes donnent des résultats très semblables. La méthode A nous semble donc à privilégier dans la mesure où, découplée par plan de coupe, elle permet de réduire les calculs aux petits nombres de données expérimentalement par coupe.

Annexe A Interpolation polynomiale d’hermite cubique périodique par morceaux : les fonctions perpchip et perspline

For more details, the reader is referred to [Bas14]111\citefield46580percubint_Jbastien_2014note.

Let $N\in\mathbb{N}^{*}$ with $n\geq 2$ and ${(x_{i})}_{1\leq i\leq N}$, ${(y_{i})}_{1\leq i\leq N}$ be $2N$ reals numbers. We assume without loss of generality that

[TABLE]

We try to find a function $f$ such that

[TABLE]

This interpolation problem is very classical and there are a lot of solutions. For example, cubic piecewise interpolation are provided by the following matlab functions :

• $\bullet$

spline, written by Carl de Boor ; see [Boo01] or [Mol04] ;

• $\bullet$

pchip ( so called for << Piecewise Cubic Hermite Interpolating Polynomial >>) ; see [FC80, chap. 3], [KMN88, p. 104] or [Mol04].

These both functions ensure that $f$ is of class $C^{1}$ on the interval $[x_{1},x_{N}]$. Moreover, (see the help of function pchip) the functions supplied by spline and pchip are constructed in exactly the same way, except that the slopes at the $x_{i}$ are chosen differently. This has the following effects: spline is smoother, i.e., $f^{\prime\prime}$ is continuous. spline is more accurate if the data are values of a smooth function. pchip has no overshoots and less oscillation if the data are not smooth. pchip is less expensive to set up. Moreover, we assume that now

[TABLE]

Thus, we try to determine $f$ which is $T$-periodical (with $T=x_{N}-x_{1}$). This property implies that

[TABLE]

The functions pchip and spline of matlab are adapted to the periodical case: perpchip and perspline.

See for example figures 21 and 22.

Annexe B Développement du polynôme $P$ en fonction des valeurs particulière de $A$ et de $B$

Le polynôme $P$ est donné par (4.10) et (4.11). On peut utiliser la formule du binôme de newton pour donner la forme générale des coefficients $c_{i}$ et des exposants $\alpha_{i}$ et $\beta_{i}$, mais seules un nombre finis de valeurs de $A$ et de $B$ seront utilisés. On vérifie donc à la main que

[TABLE]

On en déduit donc les valeurs de $N$, ${(c_{i})}_{0\leq i\leq N}$, ${(\alpha_{i})}_{0\leq i\leq N}$ et ${(\beta_{i})}_{0\leq i\leq N}$ dans le tableau 4.

1

• Paris 1

• Dunod

1

• Springer-Verlag, New York

1

• Bréal

1

• Ellipse

1

• Prentice Hall

2

• Society for IndustrialApplied Mathematics, Philadelphia, PA

1

• New York 1

• John Wiley and Sons

Références

• [ARW01] Sung Joon Ahn, Wolfgang Rauh and Hans-Jürgen Warnecke

“Least-squares orthogonal distances fitting of circle, sphere, ellipse, hyperbola, and parabola”

In Pattern Recognition 34.12, 2001, pp. 2283–2303

• [Bas07] Jérôme Bastien

“Applications de l’algèbre et de l’analyse à la géométrie.” Notes de cours de l’UV MT25 de l’UTBM, disponible sur le web : http://utbmjb.chez-alice.fr/UTBM/index.html, rubrique MT25, 2007

URL: http://cel.archives-ouvertes.fr/docs/00/55/54/85/PDF/coursMT25_P07.pdf

• [Bas14] Jérôme Bastien

“Periodical Piecewise Cubic Hermite Interpolating Polynomial: the functions perpchip and perspline.” Les fonctions perpchip et perspline, ainsi que l’aide periodical_cubic_interpolation.pdf, sont disponibles sur http://www.mathworks.com/matlabcentral/fileexchange/46580-periodical-cubic-interpolation, 2014

• [BM03] Jérôme Bastien and Jean-Noël Martin

“Introduction à l’analyse numérique” Ouvrage disponible à la bibliothèque Sciences de Lyon 1 (cote : 519.4 BAS, 4 ième étage), 2003

• [Boo01] Carl Boor

“A practical guide to splines” 27, Applied Mathematical Sciences, 2001, pp. xviii+346

• [CJ89] Maurice G. Cox and Helen M. Jones

“An algorithm for least-squares circle fitting to data with specified uncertainty ellipses”

In IMA J. Numer. Anal. 9.3, 1989, pp. 285–298

DOI: 10.1093/imanum/9.3.285

• [Cre+12] Thomas Creveaux, Jérôme Bastien, Clément Villars and Pierre Legreneur

“Model of joint displacement using sigmoid function. Experimental approach for planar pointing task and squat jump”, 2012

arXiv: http://arxiv.org/abs/1207.2627

• [EMM07] H. Endo, T. Murahashi and E. Marui

“Accuracy estimation of drilles holes with small diameter and influence of drill parameter on the machining accuracy when drillin in mild steel shett”

In Machine Tools and manufacture 47, 2007, pp. 175–181

• [Fan90] De Liang Fan

“On formulas for calculating parameters of least square circles”

In J. Southeast Univ. 20.6, 1990, pp. 96–101

• [FC80] F.. Fritsch and R.. Carlson

“Monotone piecewise cubic interpolation”

In SIAM J. Numer. Anal. 17.2, 1980, pp. 238–246

DOI: 10.1137/0717021

• [GAJ94] Daniel Guinin, François Aubonnet and Bernard Joppin

“Précis de mathématiques, Tome 5 : Géométrie”, 1994

• [Hen99] Donald M. Henderson

“Estimating the masses and centers of mass of extinct animals by 3-D mathematical slicing”

In Paleobiology 25.1, 1999, pp. 88–106

• [HM13] Frédéric Holweck and Jean-Noël Martin

“Géométries pour l’ingénieur”, 2013

• [HNA07] John R Hutchinson, Victor Ng-Thow-Hing and Frank C Anderson

“A 3D interactive method for estimating body segmental parameters in animals: application to the turning and running performance of Tyrannosaurus rex.”

In J. Theor. Biol. 246.4, 2007, pp. 660–680

DOI: 10.1016/j.jtbi.2007.01.023

• [KMN88] David Kahaner, Cleve Moler and Stephen Nash

“Numerical Methods and Software”, 1988

• [LHB12] P. Legreneur, D.. Homberger and V. Bels

“Assessment of the mass, length, center of mass, and principal moment of inertia of body segments in adult males of the Brown Anole (Anolis sagrei) and Green, or Carolina, Anole (Anolis carolinensis)”

In Journal of Morphology 273, 2012, pp. 765–775

• [MK91] L. Moura and R. Kitney

“A direct method for least-squares circle fitting”

In Comput. Phys. Comm. 64.1, 1991, pp. 57–63

DOI: 10.1016/0010-4655(91)90049-Q

• [Mol04] Cleve B. Moler

“Numerical computing with MATLAB”, 2004, pp. xii+336

• [NF97] V. Ng-Thow-Hing and E. Fiume

“Interactive display and animation of $B$-spline solids as muscle shape primitives”

In Proc. Comp. Anim. Simul., 1997, pp. 81–97

• [Raz97] A. Razet

“Résolution analytique d’un cercle de moindres carrés pour une utilisation en interferométrie” Bureau National de Métrologie, Conservatoire National des Arts et Métier, 292 rue Saint-Martin, 75141 Paris Cedex 03, France

In Bulletin du Bureau National de Métrologie 108, 1997, pp. 39–48

• [Raz98] A. Razet

“Analytical resolution of least-square applications for the circle in interferometry and radiometry”

In Metrologia 35, 1998, pp. 143–149

• [WD87] J.. Wells and D.. DeMenthon

“Measurement of body segment mass, center of gravity, and determination of moments of inertia by double pendulum in Lemur fulvus”

In American Journal of Primatology 12, 1987, pp. 299–308

• [Win09] David Winter

“Biomechanics and Motor Control of Human Movement”, 2009

• [ZS83] V.. Zatsiorsky and V. Seluyanov

“The mass and inertia characteristics of the main segments of the human body”

In Biomechanics 8B, 1983, pp. 1152–1159