Previous Up Next

Chapitre 26  Tracer de courbes algébriques avec les bases de Bernstein

Ce TP porte sur le tracer des courbes algébriques implicites définies par une équation algébrique f(x,y)=0. Cette courbe est considérée comme la courbe d’intersection entre le plan z=0 et la surface z=f(x,y), avec 0≤ x, y ≤ 1; de plus le polynôme f est exprimé dans la base tensorielle de Bernstein: (B0(n)(x), B1(n)(x),… Bn(n)(x)) × (B0(m)(y), B1(m)(y),… Bm(m)(y)), si bien que la surface est une surface de Bézier; elle est donc incluse dans l’enveloppe convexe de ses points de contrôle.

Quelques mots d’explication:

La fonction est considérée comme:

z=
i=n
i=0
j=m
j=0
 zi,j Bi(n)(xBj(m)(y)

draw_curve_be prend en arguments un tableau à 2 dimensions des zi,j (les coefficients dans la base de Bernstein) et un niveau de récursion maximum pour la subdivision.

casteljau : applique la méthode de de Casteljau sur un tableau à 1 dimension de coefficients dans la base de Bernstein. Elle rend une paire: le tableau des coefficients de la moitié gauche, le tableau des coefficients de la moitié droite. Elle utilise la fonction "pyramide" qui pour un tableau t rend la liste: t, vecteur des milieux de t, etc.

casteljau_mat applique la méthode de de Castejau sur une matrice (2 dimensions) de coefficients dans la base de Bernstein; elle coupe en 2 chaque ligne de la matrice, et rend la matrice des moitiés gauches et des moitiés droites.

casteljau_mat’ : idem, mais la division se fait selon l’autre axe. Par paresse, j’utilise la transposition de matrices.

casteljau_surface coupe en 4 la matrice des coefficients; chaque matrice fille est la matrice des coefficients d’un quartier: nord-est, nord-ouest, sud-est, sud-ouest.

Il y a (au moins) deux façons de convertir vers la base de Bernstein.

26.1  Conversion vers la base de Bernstein

26.1.1  Conversion par interpolation

Les valeurs de la fonction f(x,y) sont calculées aux sommets d’une grille régulière (la même grille que celle des valeurs de contrôle de la surface de Bézier). Les valeurs de contrôle de la surface de Bézier, i.e., les coefficients dans la base de Bernstein, en sont déduits.

draw_implicit_curve (x0,y0) cote xdegre ydegre equation dessine une courbe implicite. equation est une fonction qui attend 2 arguments (flottants), x et y et rend la valeur (flottante, le z) de la fonction en (x,y). Cette fonction est supposée polynomiale... Les autres arguments décrivent le domaine x, y dans lequel afficher la courbe.

grille_des_z (x0, x1) (y0,y1) xdeg ydeg fonc : étant donnée une fonction: fonc, cette fonction calcule une grille des valeurs de la fonction, aux sommets de la grille. Les bornes de la grille sont (x0, x1) en x, et (y0,y1) en y. xdeg et ydeg sont les degrés du polynôme. Avec cette grille

Gijfonc(j/xdegi/ydeg)

il est possible de calculer la matrice des coefficients de la fonction dans la base de Bernstein.

coeff_bernstein_curve attend un tableau des valeurs d’une fonction f (f est un polynôme de degré n) en i/n, pour i=0, … n. Elle retourne les valeurs des coefficients bi de f dans la base de Bernstein, tels que

f(x)=
n
i=0
 bi Bi(n)(x)

coeff_bernstein_surface attend une matrice des valeurs d’une fonction f (f est un polynôme de degré n en x, et m en y). La matrice a n+1 lignes, où n est le degré en x. Elle a m+1 colonnes, où m est le degré en y. Chaque ligne de la matrice décrit un tableau en y de m éléments. Cette matrice est Gl,c=f(c/n,l/m), avec . Elle retourne une matrice des coefficients bij de f dans la base tensorielle de Bernstein, tels que

f(x,y)=
n
i=0
m
j=0
 bij Bi(n)(xBj(m)(y)

Avec par exemple n=3 et m=2 :

f(x,y)=
B0(n)(x), B1(n)(x), B2(n)(x),  B3(n)(x)




b00b01b02
b10b11b12
b20b21b22
b30b31b32







B0(m)(y)
B1(m)(y)
B2(m)(y)



Et donc:

  




f(0,0)f(0,1/2)f(0,1) 
f(1/3,0)f(1/3,1/2)f(1/3,1) 
f(2/3,0)f(2/3,1/2)f(2/3,1) 
f(1,0)f(1,1/2)f(1,1) 
 



  




B0(3)(0)B1(3)(0)B2(3)(0)B3(3)(0) 
B0(3)(1/3)B1(3)(1/3)B2(3)(1/3)B3(3)(1/3) 
B0(3)(2/3)B1(3)(2/3)B2(3)(2/3)B3(3)(2/3) 
B0(3)(1)B1(3)(1)B2(3)(1)B3(3)(1) 
 







b00b01b02
b10b11b12
b20b21b22
b30b31b32







B0(2)(0)B0(2)(1/2)B0(2)(1)  
B1(2)(0)B1(2)(1/2)B1(2)(1) 
B2(2)(0)B2(2)(1/2)B2(2)(1) 



ou plus concisément: G=Bx b By. Donc la matrice des coefficients dans la base de Bernstein est: b=Bx−1 G By−1. Remarque: il est possible de calculer Bx−1 et By−1 sans recours à l’inverse, mais l’auteur était débordé... Dans le source, Bx est appelée mat_x , By est appelée mat_y.

Post scriptum: ces matrices Bx−1 et By−1 sont remarquables; bien sûr elles sont rationnelles, et les multiplier par des entiers appropriés donne des matrices entières (à l’imprécision près des calculs). Il existe vraisemblablement une formule qui donne la valeur de chaque case, sans avoir à calculer explicitement l’inverse. Ces matrices entretiennent vrisemblablement des liens avec des matrices spéciales de Vandermonde, les matrices de Pascal, et leurs inverses.

Post scriptum 2: une erreur fréquente dessine les nappes z=f(x,y) au mauvais domaine (x, y). Pour éviter cette erreur, le plus simple et le plus radical est d’encoder dans la grille des points de contrôle non seulement les valeurs zij, mais aussi les coordonnées xij et yij. Pour éliminer tout risque d’erreur, on peut encoder les points x, y, par un enregistrement plutôt que par un tableau. Il est facile de modifier la fonction recherchant le min-max en z pour qu’elle rende les coordonnées x, y, z des sommets min et max de la boite englobante. Il est certes un peu ridicule d’appliquer la méthode de de Casteljau sur les coordonnées x et y (elles varient linéairement, alors que la fonction f(x,y) peu être de plus haut degré), mais cela ne modifie pas la complexité théorique. De plus cette solution simple donne rapidement un programme correct, que vous pouvez toujours améliorer par la suite. Cf variante brute plus bas.

Post scriptum 3: Un calcul d’inverse par la méthode LUP, et de résolution de système linéaire, est fourni, dans le fichier lu_decompose.ml.

26.1.2  Conversion de la base canonique à la base tensorielle de Bernstein

Il est aussi possible de convertir de la base canonique à la base de Bernstein. Dans ce paragraphe, C(k,n) est le nombre de sous ensembles de k éléments dans un ensemble de n éléments.

Pour x entre 0 et 1, la contribution de xk à Bin(x) est nulle si k>n, et sinon c’est C(k,i) / C(k,n).

Si l’intervalle étudié n’est pas [0, 1], mais [x0, x1], la contribution de xk à Bin(x) est

k
j=0
 C(j,k)× (x1x0)j × x0kj × C(j,i) / C(j,n)

La contribution d’un monôme xayb à Bim(x)Bjn(y) est le produit des contributions de xa à Bim(x) et de yb à Bjn(y).

Le fichier to_bernstein.ml fournit une conversion qui ne nécéssite pas de résolution de système linéaire.

26.2  Remarques finales

verbose: trace-t-on les rectangles étudiés (sans les remplir).

Ce programme est très largement améliorable :

- on peut ne pas couper toujours en 4; par exemple, quand la fonction est constante ou presque selon x, il est inutile de couper selon x;

- la subdivision peut être stoppée dès que le réseau de contrôle est monotone en x et en y; variante plus ou moins équivalente : calculer les dérivées en x et en y; quand le vecteur gradient (fx, fy) ne s’annule pas dans le domaine, il est inutile de subdiviser : la courbe est "simple" dans le domaine.

- il est possible de réduire en considérant l’enveloppe convexe 2D des 2 projections des points de contrôle, sur Oxz d’une part, et sur Oyz d’autre part.

- plutôt que l’algorithme de de Casteljau, il est possible de multiplier par les matrices de de Casteljau.

- etc...

Le principe de la méthode (recours au bases de Bernstein, subdivision, enveloppe convexe) est aussi utilisable pour résoudre f(x,y)=g(x,y)=0. Une optimisation importante est de préconditionnement : dans un domaine, f(x,y)=g(x,y)=0 a les mêmes racines que a f(x,y)+ b g(x,y)=c f(x,y)+ d g(x,y)=0; a, b, c, d sont choisies pour que les deux courbes : a f + b g=0 et c f + d g=0 soient "parallèles" aux axes, autant que possible; les coefficients sont donnés par l’inverse du jacobien, au centre du domaine étudié : si la fonction est linéaire, alors les deux courbes obtenues seront des droites parallèles aux axes. Avec la réduction par l’enveloppe convexe, la convergence devient super quadratique, (preuve: cf Bernard Mourrain).

Ce problème f(x,y)=g(x,y)=0 se pose lors du calcul de l’intersection entre une droite et une surface algébrique paramétrée, en lancer de rayons par exemple. Cette méthode est aussi utilisable pour le calcul des points d’intersection entre 3 surfaces paramétrées (bicubic patches), ou des points particuliers (ex: normale verticale, points singuliers...) sur une courbe d’intersection entre 2 surfaces paramétrées.

Cette méthode est utilisable en 3D, pour trianguler des surfaces implicites. Elle est très fiable (elle n’oublie pas de composante connexe) mais un peu lente; il vaut mieux utiliser cette méthode pour trouver des points germes sur la surface (les points singuliers, les points avec une normale parallèle à un axe), et ensuite recourir à une triangulation par suivi : triangles marcheurs, cubes marcheurs, tétraèdres marcheurs (marching triangles, marching cubes, marching tetrahedra); les points particuliers sont choisis de façon à garantir qu’aucune composante connexe n’est oubliée.

26.3  Les sources

26.3.1  Conversion de la base canonique vers la base de Bernstein: to_bernstein.ml

(* THIS IS THE FILE to_bernstein.ml *) module A=Array;; module L=List;; let foi= float_of_int;; let rec puiss a n = if n=0 then 1. else if n=1 then a else if 0= n mod 2 then puiss (a *. a) (n/2) else a *. (puiss (a *. a) (n/2));; let matrice nblig nbcol fonc_lc = A.init nblig (function lig-> A.init nbcol (function col -> fonc_lc lig col));; let sum_float i0 i1 fonc = let rec add i accu = if i>i1 then accu else add (i+1) (accu +. (fonc i)) in add i0 0.;; let product_float i0 i1 fonc = let rec multiply i accu = if i>i1 then accu else multiply (i+1) (accu *. (fonc i)) in multiply i0 1.;; module Choices= struct let rec int_choices i n = if i=0 then 1 else if i>n then 0 else if i=1 then n else ((int_choices (i-1) n) * (n - i + 1)) / i ;; let float_choices i n = foi (int_choices i n);; let max_degree=21;; let table_float = matrice (max_degree+1) (max_degree+1) (fun i n -> float_choices i n );; let table_int = matrice (max_degree+1) (max_degree+1) (fun i n -> int_choices i n );; let fchoices i n = if n<= max_degree then table_float.(i).(n) else float_choices i n;; let choices i n = if n<= max_degree then table_int.(i).(n) else int_choices i n;; end;; (* returns the contribution of x^k to B_i^(n)(x) *) let contribution_01 k i n = let module C = Choices in if k> n then 0. else (C.fchoices k i) /. (C.fchoices k n);; let contribution_interval x0 x1 k i n= (* cas 1D : ici x0 et x1 sont des nombres *) let module C = Choices in let a=x1 -. x0 and b=x0 in sum_float 0 k (function j-> (C.fchoices j k) *. (puiss a j) *. (puiss b (k-j)) *. (C.fchoices j i) /. (C.fchoices j n));; (* returns the multi dimensional contribution of a degree x_0^d_0 x_1^d_1... to the B_i0^(n0)(x0) B_i1^(n1)(x1)... *) let md_contribution ks is ns = product_float 0 (A.length ks - 1) (function i -> contribution_01 ks.(i) is.(i) ns.(i)) (* It is : (contribution ks.(0) is.(0) ns.(0)) *.(contribution ks.(1) is.(1) ns.(1)) *. etc *) let md_contribution_box x0s x1s ks is ns= (* cas n dimensions: x0, x1, etc sont des vecteurs *) product_float 0 (A.length ks - 1) (function i -> contribution_interval x0s.(i) x1s.(i) ks.(i) is.(i) ns.(i))

26.3.2  Le dessin de courbes algébriques: draw_implicite_curve.ml

(* FILE draw_implicit_curve.ml *) module A=Array;; module L=List;; module Lup=Lu_decompose;; open Lup;; let milieu a b = (a +. b) /. 2. ;; let vmilieux tab = vecteur (A.length tab - 1) (function i-> milieu tab.(i) tab.(i+1));; let rec choix i n = if i=0 then 1 else if i>n then 0 else if n-i < i then choix (n-i) n else ((choix (i-1) n) * (n-i+1))/ i ;; let foi= float_of_int ;; let rec puiss x k = if (k<0) then puiss (1. /. x) (0 - k) else if k=0 then 1. else if k mod 2 = 0 then puiss (x *. x) (k/2) else x *. (puiss x (k-1));; (* rend B_k(x) pour le degre n *) let bernstein n k x = (foi (choix k n )) *. (puiss x k) *. (puiss (1. -. x) (n - k));; let minmax_tab tab = let mini = ref tab.(0).(0) in let maxi = ref tab.(0).(0) in for i=0 to A.length tab - 1 do for j=0 to A.length tab.(i) - 1 do mini := min !mini tab.(i).(j); maxi := max !maxi tab.(i).(j) done done; (!mini, !maxi);; (* version longue let casteljau polygoneControle= let n= Array.length polygoneControle in let couches= Array.make n polygoneControle in for niveau=1 to n-1 do couches.(niveau) <- Array.make (n-niveau) couches.(niveau-1).(0); for i=0 to n-niveau-1 do couches.(niveau).(i) <- milieu couches.(niveau-1).(i) couches.(niveau-1).(i+1) done done; let pol1=Array.make n polygoneControle.(0) in let pol2=Array.make n polygoneControle.(n-1) in for niveau=0 to n-1 do pol1.(niveau) <- couches.(niveau).(0) done; for j=0 to n-1 do pol2.(j) <- couches.(n-j-1).(j) done; pol1, pol2 ;; *) (* pyramide t rend la liste: [ t; les milieux de t; les milieux des milieux.... ] *) let rec pyramide t = if A.length t = 1 then [t] else t::(pyramide (vmilieux t));; let casteljau tableau = let pyra = A.of_list (pyramide tableau) in let n=A.length tableau in (vecteur n (function i-> pyra.(i).(0)), vecteur n (function i-> pyra.(n-1-i).(i)) );; let casteljaumat mat = let vect_of_pairs = Array.map casteljau mat in let mat1 = Array.map fst vect_of_pairs in let mat2 = Array.map snd vect_of_pairs in mat1, mat2 ;; let casteljaumat' mat = let m1,m2 = casteljaumat (transpose mat) in transpose m1, transpose m2;; let casteljau_surface mat = let m1,m2=casteljaumat mat in let m11,m12=casteljaumat' m1 in let m21,m22=casteljaumat' m2 in m11,m12,m21,m22 ;; module G=Graphics;; let palegoldenrod=G.rgb 238 232 170;; let trace_rectangle (x0,y0) cote =G.draw_rect x0 y0 cote cote;; let remplir_rectangle (x0,y0) cote =G.fill_rect x0 y0 cote cote;; let verbose= ref true;; (* les z_control sont dans la base de Bernstein *) let draw_curve_be z_control nivmax = let rec study (x0,y0) cote coeffs niveau = if !verbose then trace_rectangle (x0,y0)cote; let (mi,ma)= minmax_tab coeffs in if mi > 0. || ma < 0. then () else if niveau=nivmax then remplir_rectangle (x0,y0) cote else begin let so,no,se,ne= casteljau_surface coeffs in let demicote = cote / 2 in let xm=(x0+x0+cote)/2 and ym=(y0+y0+cote)/2 in study (x0,y0) demicote so (niveau+1); study (x0,ym) demicote no (niveau+1); study (xm,y0) demicote se (niveau+1); study (xm,ym) demicote ne (niveau+1); end in G.clear_graph(); G.set_color palegoldenrod; G.fill_rect 0 0 (G.size_x()) (G.size_y()) ; G.set_color G.black; study (0,0) (min (G.size_x()) (G.size_y())) z_control 0;; (*====================================================================================*) (* equation est une (fun x y -> ... ); il faut donner le degre, car comment le deviner? cette fonction rend le tableau des points de controle dans la base de Bernstein elle calcule une grille reguliere des valeurs de la fonction et en deduit les points de controle. Voir le polycopie. Remarque 1: les matrices inversees sont toujours les memes, elles ne dependent que du degre. Elles ont une tete sympathique, il doit etre possible de trouver une formule et de la prouver par recurence. Remarque 2: Il y a d'autres facons de proceder, sans doute davantage stables numeriquement. Une autre nethode est donnee plus loin, dans bernstein_of_polynomial *) let to_bernstein (x0,x1) (y0,y1) xdeg ydeg fonc_xy = let coeff_bernstein_surface z_ij = let nl=A.length z_ij in let xdeg=nl-1 in let xdegf=foi xdeg in let nc=A.length z_ij.(0) in let ydeg=nc-1 in let ydegf=foi ydeg in let mat_x=matrice nl nl (fun lig col -> bernstein xdeg col (foi lig /. xdegf)) in let mat_y= matrice nc nc (fun lig col -> bernstein ydeg lig (foi col /. ydegf)) in let inv_mat_x = Lup.inverse_mat mat_x in let inv_mat_y = Lup.inverse_mat mat_y in mat_mat inv_mat_x (mat_mat z_ij inv_mat_y) in let xdegf= foi xdeg and ydegf= foi ydeg in let zij= matrice (xdeg+1) (ydeg+1) (fun i j -> let xx=(foi i)/. xdegf in let x= xx *. (x1 -. x0) +. x0 in let yy=(foi j)/. ydegf in let y= yy *. (y1 -. y0) +. y0 in fonc_xy x y) in (* bon, c'est dans le mauvais sens, donc je prend la transposee *) transpose (coeff_bernstein_surface zij);; let draw_implicit_curve (x0,y0) cote xdeg ydeg equation = draw_curve_be (to_bernstein (x0,x0 +. cote) (y0, y0 +. cote) xdeg ydeg equation) 10;; (*====================================================================================*) (* rend la matrice de Bernstein: M.(l).(c) = bernstein_c( l/degre) pour le degre specifie *) let mat_bernstein degre = matrice (degre+1) (degre+1) (fun lig col -> bernstein degre col (foi lig /. (foi degre)));; (*====================================================================================*) G.open_graph " 513x513";; let demo()= verbose:=true; draw_implicit_curve (-1.2,-1.2) 2.4 2 2 (fun x y -> x *. x +. y *. y -. 1.); draw_implicit_curve ( -3., -3.) 6. 2 2 (fun x y -> x *. x +. y *. y /. 3. -. 1.); draw_implicit_curve ( -3., -3.) 6. 1 1 (fun x y -> x *. y -. 1. ); draw_implicit_curve ( -3., -3.) 6. 2 1 (fun x y -> y -. x *. x /. 3. ); draw_implicit_curve ( -3., -3.) 6. 1 2 (fun x y -> y *. y /. 3. -. x ); draw_implicit_curve (-3.,-3.) 6. 3 3 (fun x y -> x *. x *. x +. y *. y *. y -. 1.); draw_implicit_curve (-3.,-3.) 6. 3 3 (fun x y -> x *. x *. x +. y *. y -. 1.); verbose := false; draw_implicit_curve (-3.,-3.) 6. 3 3 (fun x y -> x *. x *. x -. y *. y -. 1.); draw_implicit_curve (-3.,-3.) 6. 3 3 (fun x y -> x *. x *. x -. y *. y *. x -. 1.); draw_implicit_curve (-3.,-3.) 6. 3 3 (fun x y -> x *. x *. y -. y *. y *. x -. 1.); draw_implicit_curve (-3.,-3.) 6. 3 3 (fun x y -> x *. x *. y -. y *. y *. x +. x *. y -. 1. ); draw_implicit_curve (-3.,-3.) 6. 2 4 (fun x y -> x *. x *. y *. y -. y *. y *. x +. x *. y -. 1. ); draw_implicit_curve (-3.,-3.) 6. 3 3 (fun x y -> x *. x *. x *. y *. y -. y *. y *. x +. x *. y -. 1. ); draw_implicit_curve (-3.,-3.) 6. 10 10 (fun x y -> (puiss x 10) +. (puiss y 10) -. 10. *. (puiss x 9) *. (puiss y 8) +. 3. *. (puiss x 9) -. 1. );; if not !Sys.interactive then demo();; (*====================================================================================*) (* ici , on utilise des polynomes= des listes de monomes, et le module To_bernstein pour les convertir dans la base de Bernstein *) module Be=To_bernstein;; type polynom = (float * int * int) list;; (* le monome est: (a coeff, degre en x, degre en y *) let degre_x polynome = L.fold_left (fun deg (_,d,_)-> max deg d) 0 polynome;; let degre_y polynome = L.fold_left (fun deg (_,_,d)-> max deg d) 0 polynome;; let bernstein_of_polynomial (x0,x1) (y0,y1) polynom = let xdeg= degre_x polynom in let ydeg= degre_y polynom in let control_points = matrice (xdeg+1) (ydeg+1) (fun a b -> L.fold_left (fun somme (coeff, kx,ky) (* coeff * x^kx y^ky *) -> let contribution_monome = Be.md_contribution_box [|x0;y0|][|x1;y1|] [|kx;ky|] [|a;b|] [|xdeg;ydeg|] in somme +. coeff *. contribution_monome ) 0. polynom ) in transpose control_points ;; let draw_algebraic_curve (x0,x1) (y0,y1) polynom= draw_curve_be (bernstein_of_polynomial (x0,x1) (y0,y1) polynom) 9;; let demo2 ()= let cercle=[ (1., 2,0); (1., 0,2); (-1., 0,0) ] in draw_algebraic_curve (-1.1, 1.1) (-1.1,1.1) cercle; let parabole= [ 1.,0,1; -1.,2,0] in draw_algebraic_curve (-1.,1.)(-0.5,1.5) parabole; ;; if not !Sys.interactive then demo2();;

26.3.3  Variante brute: draw_implicite_curve_brute.ml

Dans la version brute, les valeurs de contrôle (les éléments de la grille) sont des points (x, y, z). Ce fichier est complet.

(* FILE draw_implicit_curve_brute.ml *) (* pour s'eviter de dessiner les champs d'altitude au dessus des mauvais domaines..., chaque point de controle comporte son x et son y *) module A=Array;; module L=List;; let vecteur=A.init;; let matrice nl nc fonc = vecteur nl (function lig -> vecteur nc (function col -> fonc lig col));; let transpose m = matrice (A.length m.(0)) (A.length m) (fun l c -> m.(c).(l));; type point=float*float*float;; let milieu_point (ax,ay,az) (bx,by,bz) = ((ax +. bx) /. 2., (ay +. by)/. 2., (az +. bz) /. 2.) ;; let vmilieux tab = vecteur (A.length tab - 1) (function i-> milieu_point tab.(i) tab.(i+1));; let rec choix i n = if i=0 then 1 else if i>n then 0 else if n-i < i then choix (n-i) n else ((choix (i-1) n) * (n-i+1))/ i ;; let foi= float_of_int ;; let iof=int_of_float;; let choices i n = foi( choix i n);; let rec puiss x k = if (k<0) then puiss (1. /. x) (0 - k) else if k=0 then 1. else if k mod 2 = 0 then puiss (x *. x) (k/2) else x *. (puiss x (k-1));; (* rend B_k(x) pour le degre n *) let bernstein n k x = (foi (choix k n )) *. (puiss x k) *. (puiss (1. -. x) (n - k));; let xof (x,_,_) = x;; let yof (_,y,_) = y;; let zof (_,_,z) = z;; let minimum_l fonc liste = L.fold_left (fun mini x-> min mini (fonc x)) (fonc (L.hd liste)) (L.tl liste);; let maximum_l fonc liste = L.fold_left (fun maxi x-> max maxi (fonc x)) (fonc (L.hd liste)) (L.tl liste);; (* rend les coordonnees des cois de la boite englobante : (x0,y0,z0), (x1,y1,z1) *) let xyz_minmax_mat tab = let zmini = ref (zof tab.(0).(0)) in let zmaxi = ref (zof tab.(0).(0)) in for i=0 to A.length tab - 1 do for j=0 to A.length tab.(i) - 1 do zmini := min !zmini (zof tab.(i).(j)); zmaxi := max !zmaxi (zof tab.(i).(j)) done done; let n0 = A.length tab in let n1=A.length tab.(0) in let sommets = [ tab.(0).(0); tab.(0).(n1-1); tab.(n0-1).(0); tab.(n0-1).(n1-1)] in let xmini = minimum_l xof sommets in let xmaxi = maximum_l xof sommets in let ymini = minimum_l yof sommets in let ymaxi = maximum_l yof sommets in ((xmini, ymini, !zmini), (xmaxi, ymaxi, !zmaxi));; (* pyramide t rend la liste: [ t; les milieux de t; les milieux des milieux.... ] *) let rec pyramide t = if A.length t = 1 then [t] else t::(pyramide (vmilieux t));; let casteljau tableau = let pyra = A.of_list (pyramide tableau) in let n=A.length tableau in (vecteur n (function i-> pyra.(i).(0)), vecteur n (function i-> pyra.(n-1-i).(i)) );; let casteljaumat mat = let vect_of_pairs = Array.map casteljau mat in let mat1 = Array.map fst vect_of_pairs in let mat2 = Array.map snd vect_of_pairs in mat1, mat2 ;; let casteljaumat' mat = let m1,m2 = casteljaumat (transpose mat) in transpose m1, transpose m2;; let casteljau_surface mat = let m1,m2=casteljaumat mat in let m11,m12=casteljaumat' m1 in let m21,m22=casteljaumat' m2 in m11,m12,m21,m22 ;; module G=Graphics;; G.open_graph " 513x513";; let palegoldenrod=G.rgb 238 232 170;; type affichage= { mutable tx : float; mutable ty: float; mutable ech : float };; let aff= { tx =0.; ty=0.; ech=1. };; let set_window (x0,y0)(x1,y1)= let sx= float_of_int (G.size_x()) and sy= float_of_int (G.size_y()) in let dx = x1 -. x0 and dy = y1 -. y0 in aff.ech <- min (sx/. dx) (sy /. dy ); aff.tx <- 0.5 *. ( sx -. aff.ech *. (x0 +. x1)); aff.ty <- 0.5 *. ( sy -. aff.ech *. (y0 +. y1));; let to_screen (x,y) = (iof (aff.ech *. x +. aff.tx), iof (aff.ech *. y +. aff.ty ));; let dessiner_rectangle fonc (x0,y0) (x1,y1) = let (x0,x1)= (min x0 x1, max x0 x1) in let (y0,y1)= (min y0 y1, max y0 y1) in let xx0,yy0= to_screen (x0,y0) in let xx1,yy1= to_screen (x1,y1) in fonc xx0 yy0 (xx1-xx0) (yy1-yy0);; let trace_rectangle= dessiner_rectangle G.draw_rect;; let remplir_rectangle = dessiner_rectangle G.fill_rect;; let verbose= ref true;; (* les pts_control sont dans la base de Bernstein *) let draw_curve_be pts_control nivrec = let rec study coeffs niveau = let (x0,y0,z0),(x1,y1,z1)= xyz_minmax_mat coeffs in if !verbose then trace_rectangle (x0,y0) (x1,y1); if z0 > 0. || z1 < 0. then () else if niveau=0 then remplir_rectangle (x0,y0) (x1, y1) else begin let so,no,se,ne= casteljau_surface coeffs in study so (niveau-1); study no (niveau-1); study se (niveau-1); study ne (niveau-1); end in G.clear_graph(); G.set_color palegoldenrod; G.fill_rect 0 0 (G.size_x()) (G.size_y()) ; G.set_color G.black; let (x0,y0,z0),(x1,y1,z1)= xyz_minmax_mat pts_control in set_window (x0,y0) (x1,y1); study pts_control nivrec;; let cercle= [| [| 0.,0.,-2.; 200.,0., 0.; 400.,0., -2. |]; [| 0.,200., 0.; 200.,200., 3.; 400.,200., 0.|]; [| 0.,400., -2.; 200.,400., 0.; 400.,400., -2. |]|];; draw_curve_be cercle 8;; (*====================================================================================*) (* ici , on utilise des polynomes= des listes de monomes, et le module To_bernstein pour les convertir dans la base de Bernstein *) let sum_float i0 i1 fonc = let rec add i accu = if i>i1 then accu else add (i+1) (accu +. (fonc i)) in add i0 0.;; let product_float i0 i1 fonc = let rec multiply i accu = if i>i1 then accu else multiply (i+1) (accu *. (fonc i)) in multiply i0 1.;; (* returns the contribution of x^k to B_i^(n)(x) *) let contribution_01 k i n = if k> n then 0. else (choices k i) /. (choices k n);; let contribution_interval x0 x1 k i n= (* cas 1D : ici x0 et x1 sont des nombres *) let a=x1 -. x0 and b=x0 in sum_float 0 k (function j-> (choices j k) *. (puiss a j) *. (puiss b (k-j)) *. (choices j i) /. (choices j n));; (* returns the multi dimensional contribution of a degree x_0^d_0 x_1^d_1... to the B_i0^(n0)(x0) B_i1^(n1)(x1)... *) let md_contribution ks is ns = product_float 0 (A.length ks - 1) (function i -> contribution_01 ks.(i) is.(i) ns.(i)) ;; (* It is : (contribution ks.(0) is.(0) ns.(0)) *.(contribution ks.(1) is.(1) ns.(1)) *. etc *) let md_contribution_box x0s x1s ks is ns= (* cas n dimensions: x0, x1, etc sont des vecteurs *) product_float 0 (A.length ks - 1) (function i -> contribution_interval x0s.(i) x1s.(i) ks.(i) is.(i) ns.(i)) (*====================================================================================*) type polynom = (float * int * int) list;; (* le monome est: (a coeff, degre en x, degre en y *) let degre_x polynome = L.fold_left (fun deg (_,d,_)-> max deg d) 0 polynome;; let degre_y polynome = L.fold_left (fun deg (_,_,d)-> max deg d) 0 polynome;; let bernstein_of_polynomial (x0,x1) (y0,y1) polynom = let xdeg= degre_x polynom in let ydeg= degre_y polynom in let control_points = matrice (xdeg+1) (ydeg+1) (fun a b -> let z = L.fold_left (fun somme (coeff, kx,ky) (* coeff * x^kx y^ky *) -> let contribution_monome = md_contribution_box [|x0;y0|][|x1;y1|] [|kx;ky|] [|a;b|] [|xdeg;ydeg|] in somme +. coeff *. contribution_monome ) 0. polynom in let x= x0 +. (x1 -. x0) *. (foi a) /. (foi xdeg) in let y= y0 +. (y1 -. y0) *. (foi b) /. (foi ydeg) in (x,y,z) ) in (* transpose *) control_points ;; let draw_algebraic_curve (x0,x1) (y0,y1) polynom= draw_curve_be (bernstein_of_polynomial (x0,x1) (y0,y1) polynom) 9;; let demo2 ()= let cercle=[ (1., 2,0); (1., 0,2); (-1., 0,0) ] in draw_algebraic_curve (-1.1, 1.1) (-1.1,1.1) cercle; let parabole= [ 1.,0,1; -1.,2,0] in draw_algebraic_curve (-1.,1.)(-0.5,1.5) parabole; ;; if not !Sys.interactive then demo2();;

Previous Up Next