[Geométrie] Tracé d'ellipsoïdes

Bonjour,

Je ne parviens pas à (re)trouver comment on trace un ellispoïde rapidement et en étant sûr qu'il n'y a pas de "discontinuité"... J'ai cherché du côté de Bresenham sous Google, mais j'ai rien trouvé !!

D'avance merci !

Mathieu

Je ne comprends pas bien. Il s'agit de tracer une ellipsoïde (3D) sur ton écran (2D)?
Tu utilises une perspective particulière?
S'agit-il de dessiner uniquement le contour ou de donner une idée du relief?

Non non rien de tout ça !

En fait je ne veux pas la visualiser, je veux la tracer dans une image 3D.

Concrètement : je veux tracer des ellipsoides dans une image 3D pour tester des algorithmes de traitement d'image.
J'aimerais bien être sûr d'avoir des surfaces fermées, et il me semblait que Bresenham (ou autre) permettait cela...

J'ai bien penser à utiliser bêtement x²/a² + y²/b² + z²/c² = 1, mais je ne suis pas certain qu'ainsi j'obtienne à coup sûr des surfaces fermées...

C'est plus clair (ou pire) ?

Ben, en fait c'est pire...

Qu'entends-tu par image 3D? Pour moi, une image n'a que 2 dimensions... Surtout que tu as l'air de me dire que ce n'est pas la représentation en 2D d'une réalité en 3D.

En tout cas, l'équation que tu donnes est bien une surface fermée (c'est une ellipsoïde pour toutes valeurs de a,b et c.

Bon je suis pas clair là !

Je réexplique : je ne cherche pas à faire de la visualisation, je bosse en traitement d'images médicales.

En imagerie médicale, les images sont en 3D (ex : les images IRM).
J'ai des algos à faire travailler sur des images 3D où je souhaite modéliser certaines structures par des ellipses.

Le souci : si par exemple je fais de la croissance région, et s'il y a un trou au niveau de la surface (cad si il existe un voxel qui ne touche pas son voxel voisin), l'ago va foirer...

En gros, je dis à l'algo : "fais tes traitements sur telle image". Et cette image, j'aimerais bien qu'elle contienne des ellipsoides...

Donc pas de projection !

A+

Je pense que ce que tu veux faire est bien du Bresenham, mais en trois dimensions au lieu de deux. Il faut donc en chaque point de l'ellipsoïde, calculer le plan tangent en ce point, ce qui se fait très facilement en polarisant la forme quadratique qui constitue l'équation de l'ellipsoïde. Tu obtiens une forme linéaire dont les coefficients sont les coordonnées du vecteur normal à ce plan tangent. Il faut alors déterminer lequel des trois axes de coordonnées à la direction la plus proche de ce vecteur. Il suffit alors de considérer localement ta surface comme le graphe d'une fonction de deux variables qui sont les deux coordonnées restantes, comme on le fait avec Bresenham ordinaire, mais où la fonction en question n'a qu'une variable.

Yep...

Je pensais bien que c'était un truc dans ce goût là !!

Merci !

A+

en fait j'aimerai bien faire la meme chose, mais avec differents objets (ellipse, sinusoide, sphere)
mais je n'ai pas bien compris l'explication DrTopos :oops:
est-ce que vous auriez des ressources (lien, livre)
Merci

Bonjour à tous.

Après avoir un peu plus longuement réfléchi à cette question, je suis en mesure de vous donner des indications supplémentaires.

Rappelons pour commencer que Bresenham à l'origine s'applique au tracé de droites dans le plan. L'exigence fondamentale (comme l'a suggéré Mathieu) étant de ne pas avoir de trou dans le tracé. Dans le cas d'une droite (dans le plan) d'équation:

Code:

---

1 2	ax + by = c

---

qui peut s'écrire (en supposant b != 0)

Code:

---

1 2	y = (-a/b)x + (c/b)

---

on voit que la pente de la droite est (-a/b). Son module est donc |a/b|, et il est plus petit que 1 si et seulement si le module de b est plus grand que le module de a.

Si tel est le cas, on voit que y varie moins vite que x, et donc qu'en considérant cette droite comme le graphe de la fonction x |-> (-a/b)x + (c/b), et en traçant tout simplement ce graphe point par point, on est sûr de ne pas avoir de trou. Dans le cas contraire, c'est à dire |a| > |b|, on échange les rôles de x et y tout simplement.

Evidemment, dans le cas des cercles, ellipses, quadriques (ellipsoïdes, etc...) c'est un peu plus compliqué. Par exemple, pour le cercle de rayon 1 dans le plan (ensemble des complexes de module 1), on a l'équation:

Code:

---

1 2	x^2 + y^2 = 1

---

Pour dessiner ce cercle, on peut le décomposer en plusieurs arcs, et toujours pour la même raison, on n'aura pas de trou dans le tracé, si chacun de ces arcs est le graphe d'une fonction dont le module de la dérivée ne dépasse pas 1.

On voit facilement que pour le cercle, il faut faire un découpage en 4 arcs, les points de découpage correspondant aux angles suivants (en radians):

Code:

---

1 2	pi/4 3pi/4 5pi/4 et 7pi/4

---

Il y a deux arcs 'plutôt horizontaux', qu'on va tracer comme graphes des fonctions (où +- veut dire 'plus ou moins'):

Code:

---

1 2	y = +- sqrt(1 - x^2)

---

et deux arcs 'plutôt verticaux', qu'on va tracer comme graphes des fonctions:

Code:

---

1 2	x = +- sqrt(1 - y^2)

---

Les domaines de définitions des fonctions étant dans chaque cas l'intervalle [-cos(pi/4),cos(pi/4)].

Dans le cas d'une figure à trois dimensions, le principe est le même, c'est à dire décomposer la surface à tracer en morceaux 'plutôt horizontaux' quand on choisit bien la composante 'verticale' (c'est-à-dire l'un des trois axes de coordonnées). Etre 'plutôt horizontal' pour une fonction de deux variables z = f(x,y), signifie que quand x ou y varie de 1 (un pixel dans le plan des (x,y), z varie de moins que 1. C'est bien sûr ce qui assure l'absence de trous dans le graphe. En mathématique cette propriété porte un nom. On dit que la fonction est 1-lipschitzienne.

Le cas d'une surface quelconque est désespéré, mais le cas d'un ellipsoîde est tout à fait maîtrisable. C'est essentiellement dû au fait que l'ellipsoïde est convexe et aux bonnes propriétés des formes quadratiques. Voyons à titre d'exemple le cas de la sphère d'équation:

Code:

---

1 2	x^2 + y^2 + z^2 = 1

---

La décomposition va consister en 6 morceaux qu'on peut facilement se représenter. En effet, plaçons cette sphère dans le cube de coté 2 et dont les arêtes sont parallèles aux axes de coordonnées. En fait ce cube (intérieur compris) a pour équations:

Code:

---

1 2	|x| =< 1 et |y| =< 1 et |z| =< 1

---

Projetons maintenant les 6 faces du cube sur la sphère par une projection centrale (autrement-dit, chaque point du cube, qui est un vecteur de R^3, est envoyé sur le vecteur de même direction mais de module 1). On obtient les 6 morceaux de surface qui ont les bonnes propriétés pour le tracé.

Bien entendu, ce cas est trop particulier et il faut une méthode pour le cas général, et en particulier une méthode pour déterminer les domaines de définition de nos 6 fonctions. On va voir qu'il s'agit à chaque fois de l'intersection des intérieurs de deux ellipses. Le parcours point par point d'un tel domaine est plutôt facile à réaliser.

Suite dans le prochain post...

Dans le cas général, on a donc un ellipsoïde dans R^3. Son équation est de la forme:

Code:

---

1 2 3 4 5

ax^2 + by^2 + cz^2 + 2dxy + 2eyz + 2fzx + 2gx + 2hy + 2iz + j = 0

---

J'ai intentionnellement mis les termes du même degré sur la même ligne. Par ailleurs, les facteurs 2 sont inessentiels, mais c'est l'habitude de les mettre, car ils simplifient l'écriture de la forme polaire.

En fait, rien ne dit que l'équation ci-dessus est celle d'un ellipsoïde. Cette surface est certainement une quadrique, mais elle peut aussi bien être un hyperboloïde ou un paraboloïde. Dans le problème qui est posé par Mathieu, il se peut qu'on sache a priori que la surface est un ellipsoïde. C'est ce que je vais supposer. Toutefois, il existe une méthode pour le savoir (que je ne vais pas exposer dans ce post).

On aura bien compris que c'est la direction du plan tangent à la surface en un point qui détermine le morceau auquel ce point doit appartenir. Dans le cas d'une quadrique, l'équation du plan tangent s'obtient par polarisation. Toutefois, le membre de gauche de l'équation en tête de ce post n'est pas une forme quadratique, car cette formule n'est pas homogène de degré 2. Pour en faire une forme quadratique, il faut la rendre homogène par l'introduction d'une 4ème variable:

Code:

---

1 2 3 4 5

ax^2 + by^2 + cz^2 + 2dxy + 2eyz + 2fzx + 2gxt + 2hyt + 2izt + jt^2 = 0

---

Cette nouvelle équation est celle de la même quadrique, mais vue dans l'espace projectif P^3 au lieu de l'espace affine (euclidien en l'occurrence) R^3. Elle est de la forme:

Code:

---

1 2	q(x,y,z,t) = 0

---

q est bien une forme quadratique sur R^4. Sa forme polaire, c'est à dire l'unique forme bilinéaire symétrique f dont cette forme quadratique est la condensation (q(X) = f(X,X), avec X = (x,y,z,t)) est la suivante:

Code:

---

1 2 3 4 5

f((x_1,y_1,z_1,t_1),(x_2,y_2,z_2,t_2)) = ax_1x_2 + by_1y_2 + cz_1z_2 + jt_1t_2 + d(x_1y_2 + x_2y_1) + e(y_1z_2 + y_2z_1) + f(z_1x_2 + z_2x_1) + g(x_1t_2 + x_2t_1) + h(y_1t_2 + y_2t_1) + i(z_1t_2 + z_2t_1)

---

X_0 étant fixé (dans P^3), l'équation (d'inconnue X):

Code:

---

1 2	f(X_0,X) = 0

---

est par définition celle du 'plan polaire' du point X_0 par rapport à la quadrique. Dans le cas où X_0 est sur la quadrique (c'est-à-dire si q(X_0) = 0), alors ce plan est le plan tangent à la quadrique en X_0.

En posant X_0 = (x_0,y_0,z_0,t_0), on voit donc que l'équation (homogène) du plan (projectif) tangent en ce point est:

Code:

---

1 2 3 4 5

ax_0x + by_0y + cz_0z + jt_0t + d(x_0y + xy_0) + e(y_0z + yz_0) + f(z_0x + zx_0) + g(x_0t + xt_0) + h(y_0t + yt_0) + i(z_0t + zt_0) = 0

---

On peut revenir maintenant dans R^3 en faisant t = t_0 = 1. On voit alors que le vecteur normal à ce plan tangent (en X_0) a pour coordonnées:

Code:

---

1 2 3 4

n_x = ax_0 + dy_0 + fz_0 + g n_y = by_0 + dx_0 + ez_0 + h n_z = cz_0 + ey_0 + fx_0 + i

---

Les trois coordonnées x, y et z jouant des rôles semblables, je vais supposer que c'est l'axe des z qui est 'vertical', ce qui va nous permettre de tracer 2 des 6 morceaux de la surface.

L'équation de la quadrique dans R^3 est un trinôme du second degré en z. On peut donc résoudre cette équation par rapport à z, ce qui va donner une formule de la forme:

Code:

---

1 2	z = u(x,y) +- sqrt(v(x,y))

---

(que je laisse à Mathieu le soin de calculer explicitement) c'est-à-dire en fait deux fonctions, une pour chacun des deux morceaux 'plutôt horizontaux'. Je vais noter F(x,y) l'une de ces deux fonctions. On remarquera (en se remémorant la formule (-b +- sqrt(b^2 - 4ac))/2a), que u(x,y) est de degré 1 alors que v(x,y) est de degré 2.

Reste maintenant le plus important, à savoir le domaine de définition de F. Il y a deux conditions à satisfaire pour que le point de coordonnées (x_0,y_0,z_0) soit dans l'un de nos 2 morceaux. Ce sont:

Code:

---

1 2	|n_z| >= |n_x| et |n_z| >= |n_y|

---

La première de ces conditions s'écrit encore comme ceci:

Code:

---

1 2	|cF(x_0,y_0) + ey_0 + fx_0 + i| >= |ax_0 + dy_0 + fF(x_0,y_0) + g|

---

Pour éviter les valeurs absolues, elle peut être remplacée par les deux conditions:

Code:

---

1 2 3

-cF(x_0,y_0) - ey_0 - fx_0 - i <= ax_0 + dy_0 + fF(x_0,y_0) + g ax_0 + dy_0 + fF(x_0,y_0) + g <= cF(x_0,y_0) + ey_0 + fx_0 + i

---

La première se réécrit:

Code:

---

1 2	(f+c)F(x_0,y_0) + (e+d)y_0 + (f+a)x_0 +(g+i) >= 0

---

En arrangeant tout cela et en éliminant la racine carrée, on obtient l'équation d'une conique (qui ne peut être qu'une ellipse compte tenu de nos hypothèses).

Notre domaine de définition apparait donc comme une intersection d'intérieurs d'ellipses (a priori 4, mais peut-être moins). Comme les intérieurs d'ellipse sont convexes, il en est de même de cette intersection.

Pour être sûr de parcourir exhaustivement ce domaine convexe du plan des (x,y), on peut démarrer en un point du domaine. Appelons-le A. La droite parallèle à l'axe des x et passant par A coupe le domaine selon un segment (intervalle) de la forme [B,C], avec A entre B et C. On traite successivement les points de ce segment qui sont à gauche de A puis ceux qui sont à droite de A.

A partir de chacun de ces points, on recommence mais selon la direction de l'axe des y. On parcourt ainsi tout le domaine de définition de F. Il ne reste plus qu'à dessiner (dans R^3) les points du graphe.

Il y a encore une question qui est comment trouver le point de départ A ? On peut prendre un point (x,y), tel que le plan tangent à l'ellipsoïde soit vraiment horizontal. Les conditions sont donc:

Code:

---

1 2	n_x = n_y = 0

---

c'est-à-dire:

Code:

---

1 2 3

ax_0 + dy_0 + fF(x_0,y_0) + g = 0 by_0 + dx_0 + eF(x_0,y_0) + h = 0

---

Après élimination des radicaux, il s'agit de deux équations du second degré. Un tel système d'équations a a priori 4 solutions. Il faudra y regarder de plus près, pour savoir laquelle des 4 peut convenir (pas forcément la même pour les deux morceaux 'plutôt horizontaux').

Note: j'ai écrit un cours sur les formes quadratiques, qu'on trouvera sur ma page enseignement.

bonjour,

Je me demande en fait si une solution assez basique ne pourrait pas faire l'affaire :
1 - generer un ellipsoide plein dans l'image, en utilisant l'inegalite x^2/a^2 + y^2/b^2 + z^2/c^2 <= 1
2 - selectionner tous les pixels de l'ellipsoide qui sont connexes a un pixel du fond -> ca donne une surface discrete sans trou.

Ca a l'avantage de rester assez simple, et de laisser le choix de la connexite entre les voxels, ce qui est utile si on veut faire du remplissage par la suite.

A moins que je soie hors-sujet ?

A+

Kangourou >> ce que tu propose est sans doute faisable, mais pour des ellipsoïde de grand volume, cela risque d'être très lent par rapport à ce que j'ai proposé.

Citation:

---

Envoyé par DrTopos

mais pour des ellipsoïde de grand volume, cela risque d'être très lent par rapport à ce que j'ai proposé.

---

mmoui....
Pour des images 3D, le temps d'execution est pas toujours facile a prevoir.
J'ai quand meme tendance a privilegier les methodes simples, ca permet d'eviter les heures de debug.
Et si c'est pour generer une image de test, le temps d'execution est pas forcement critique.

faut voir en fonction de l'appli.

a+

Merci DrTopos
je vais m'y lancer