#open "graphics";;
include "int2hex.ml";;
(*****************************************************************)
(* tri fusion :
cmp_inf_stt : 'a -> 'a -> bool -- operateur de compraison (ex <)
v0 : 'a vect -- les elements a trier*)
(*****************************************************************)
let tri_fusion cmp_inf_stt v0 =
let permute cmp_inf_strict v = let n = vect_length v in
for i=0 to (n/2 - 1) do
let (a,b) = (v.(2*i),v.(2*i+1)) in
if cmp_inf_strict b a then (v.(2*i) <- b ;v.(2*i+1) <- a)
done;
v in
let n = vect_length v0
in let v = ref (permute cmp_inf_stt (copy_vect v0))
and w = ref (make_vect n v0.(0))
and t= ref v0 and p =ref 2 and c = ref 0
in
while !p < n do
c:=0;
while !c<n do
let j= ref (min (!c + !p) n) and i = ref !c
and k1 = min (!c + !p ) n and k2 = min (!c + 2 * !p ) n
and stop = ref false
in while not !stop do
(match (!i=k1,!j=k2) with
| (false,false) -> (if cmp_inf_stt (!v).(!i) (!v).(!j)
then ((!w).(!c) <- (!v).(!i);incr i)
else ((!w).(!c) <- (!v).(!j);incr j))
| (true,false) -> ((!w).(!c) <- (!v).(!j);incr j)
| (false,true) -> ((!w).(!c) <- (!v).(!i);incr i)
| (true,true) -> (decr c;stop:=true));
incr c
done;
done;
p:=!p * 2;
t:=!v;v:=!w;w:=!t;
done;
!v;;
(*****************************************************)
(*****************************************************)
(*****************************************************)
(*****************************************************)
(* Attendre la prochaine action *)
(*****************************************************)
let wait ()= wait_next_event [Key_pressed;Button_up];;
(*****************************************************)
(* Ca c'est Pi *)
(*****************************************************)
let Pi = 4. *. (atan 1.);;
(*****************************************************)
(*****************************************************)
(*****************************************************)
(*****************************************************
Outils pour la sauvegarde de l'ecran
et la lecture aussi
*****************************************************)
(************ Conversions ****************************)
(* coordonnee RGB d'un entier *)
let rgb2mem col =
let v=make_vect 3 0 in
v.(0) <- col land 255;
v.(1) <- (col lsr 8) land 255;
v.(2) <- (col lsr 16) land 255;
v;;
let r_g_b c = let v=rgb2mem c in (v.(2),v.(1),v.(0));;
(* utilitaires *)
let mem2rgb b g r = rgb r g b;;
let matsize m =((vect_length m),(vect_length m.(0)));;
let mem_of_int t = ((t land 0xFF) lsl 24) lor
((t land 0xFF00) lsl 8) lor
((t land 0xFF0000) lsr 8) lor
((t land 0xFF000000) lsr 24);;
let read_byte chan = int_of_char (input_char chan);;
(***************** Ecrire un point ***************)
let writeFrgb chan col =
let v = rgb2mem col in
for i=0 to 2 do output_byte chan v.(i) done;;
(***************** Sauvegarder une portion d'ecran ********)
(*
filename : string : nom du fichier destination (format TGA)
x,y : coordonnées du point en bas a gauche du rectangle
a sauvegarder.
width,height : largeur hauteur du sus-dit rectangle *)
let write_image filename x y width height =
let chan=open_out_bin filename
and n=height - 1 and p=width - 1 in
output_binary_int chan (1 lsl 9);output_binary_int chan 0;
output_binary_int chan 0;
output_byte chan width;output_byte chan (width lsr 8);
output_byte chan height;output_byte chan (height lsr 8);
output_byte chan 24;output_byte chan 0;
for i=0 to n do
for j=0 to p do
writeFrgb chan (point_color (x+j) (y+i));
done;
done;
close_out chan
;;
(********** Lire une image ********************************)
(*
filename : string : nom du fichier destination (format TGA)
x,y : coordonnees du point en bas a gauche dans lequel sera
afficher l'image
*)
let read_image filename x y=
let inchan = open_in_bin filename in
seek_in inchan 12;
let v = [| (read_byte inchan) ;(read_byte inchan) ;
(read_byte inchan) ;(read_byte inchan) |] in
let (p,n) = (v.(3) lor (v.(2) lsl 8),v.(1) lor (v.(0) lsl 8))
in seek_in inchan 18;
for i=0 to n-1 do
for j=0 to p-1 do
set_color (rgb (read_byte inchan)
(read_byte inchan)
(read_byte inchan));
plot (x+j) (y+i);
done
done;
close_in inchan;
;;
(*****************************************************)
(*****************************************************)
(*****************************************************)
(********* Operation sur les vecteurs/points ***********)
type point == float * float * float ;;
let sqr a = a *. a;;
let vadd (a,b,z1) (c,d,z2) = (a +. c , b +. d , z1 +. z2);;
let vdiff (a,b,z1) (c,d,z2) = (a -. c , b -. d , z1 -. z2);;
let opp (x,y,z) = (-.x,-.y,-.z);;
let vscal l (x,y,z) = (l *. x,l *. y,l *. z);;
let dist (x1,y1,z1) (x2,y2,z2) = sqrt(sqr(x2 -. x1) +. sqr (y2 -. y1) +. sqr (z2 -. z1));;
let norme (x,y,z) = sqrt ((sqr x) +. (sqr y) +. (sqr z));;
let pscal (x1,y1,z1) (x2,y2,z2) = x1 *. x2 +. y1 *. y2 +. z1 *. z2;;
let pvect (x1,y1,z1) (x2,y2,z2) =
( y1 *. z2 -. y2 *. z1 ,
z1 *. x2 -. z2 *. x1 ,
x1 *. y2 -. x2 *. y1 );;
let normalize v = if v = (0.,0.,0.) then
v
else
vscal (1. /. (norme v)) v ;;
(***************** Position de l'observateur *****************)
let zobs = ref 800.
and xobs = ref 0.
and yobs = ref 0.;;
let observateur = ref (!xobs, !yobs, !zobs);;
(***************** Position de l'eclairage *******************)
let spot = ref (100. , 900. , 900.);;
let pdfz = ref (-.1000.);;
let pdf = ref (0.,500.,!pdfz);;
(***************** Vecteurs unitaires ************************)
let x0 = (1.,0.,0.);;
let y0 = (0.,1.,0.);;
let z0 = (0.,0.,1.);;
(***************** Projection d'un point 3D -> 2D ************)
(*projeter : float * float * float -> float * float *)
let projeter ((x,y,z) as v) =
let k = (!zobs /. (!zobs -. z)) in
let (a,b,_) = (vscal k (vdiff v !observateur)) in
(320 + (int_of_float (a -. !xobs)),280 + (int_of_float (b -. !yobs)));;
(*let projeter ((x,y,z) as v) =
let vf (x,y,_) = normalize (vdiff v !pdf) in
let k = (!pdfz -. z) /. !pdfz in
let (a,b,_) = vadd !pdf (vscal k (vdiff v !pdf)) in
(320 + (int_of_float (a -. !xobs)),280 + (int_of_float (b -. !yobs)));;
*)
(*****************************************************)
(*****************************************************)
(*****************************************************)
(***********Tracer la ligne [a,b] ********************)
(*
a,b : float * float * float
*)
let LTracer (a,b) =
let (x1,y1) = projeter a
and (x2,y2) = projeter b in
moveto x1 y1 ; lineto x2 y2
;;
(*******Representation des polygones*******************)
(* Vis : vecteur orthogonal a la face*)
type face = {T : point vect ; Vis : point};;
(**** Barycentre d'un tableau de points -> point*******)
let barycentre T =
let n = vect_length T
and r = ref (0.,0.,0.) in
for i = 0 to n-1 do
r := vadd !r T.(i)
done;
vscal (1. /. (float_of_int n)) !r
;;
(**** Calcul du vecteur d'eclairage d'un polygone ******)
let eclairage f =
let v = (vdiff !spot (barycentre f.T)) in
normalize v;;
(**** Calcul du vecteur de vision d'un polygone ********)
let vision f =
let v = (vdiff !observateur (barycentre f.T)) in
normalize v;;
(**** Le polygone f est il visible ?********************)
let visible f =
(pscal (vision f)
f.Vis) >= 0. ;;
(**********Afficher un point ***************************)
let print_point (x,y,z) =
print_string "(";
print_float x;
print_string ";";
print_float y;
print_string ";";
print_float z;
print_string ")";;
(**********Tracer les contours d'un polygone************)
let FTracer f =
let n = vect_length f.T in
for i=0 to n-2 do
LTracer (f.T.(i),f.T.(i+1))
done;
LTracer (f.T.(n-1),f.T.(0))
;;
(*********Determiner la couleur d'une face**************)
(*** proportionnel a l'angle (Vis,eclairage) ***********)
let FColor f =
let vtemp = if (visible f) then
f.Vis
else opp f.Vis in
let arc = acos(pscal vtemp (eclairage f)) in
(* let c = if arc <= (Pi /. 2.) then
max (255 - int_of_float (255. *. arc /. (Pi /. 2.))) 0x40
else 0x40*)
let c = max (255 - int_of_float (191. *. arc /. (Pi /. 2.))) 0x40
in
(* print_point vtemp;print_point (eclairage f);
print_string " | ";print_endline (i2h c);*)
rgb c c c;;
(*********Remplir une face avec fond et un cadre de *****)
(********* couleur cadre ********************************)
let FcRemplir f fond cadre=
(****** Si c'est une ligne c'est particulier*************)
if vect_length f.T = 2 then
(set_color fond;
LTracer (f.T.(0),f.T.(1)))
else
let T2 = map_vect projeter f.T in
(set_color fond;
fill_poly T2;
set_color cadre;
FTracer f);
;;
(*** Produit vectoriel normalisee pour creer le *********)
(*** vecteur Vis d'une face *****************************)
let TrigoVis a b c =
let v = pvect (vdiff b a) (vdiff c a) in
vscal (1. /. (norme v)) v;;
(*** Afficher une pyramide de 4 faces *******************)
(**** obsolete*******************************************)
let PTracer (f1,f2,f3,f4) =
FTracer f1;
FTracer f2;
FTracer f3;
FTracer f4;;
(********* Extension du type face pour trier les polygones***)
type face2 = {F:face ; Fill_c : int ; Bord_c : int};;
(*****************************************************)
(*****************************************************)
(*****************************************************)
(******* Matrice de rotation / (1,0,0) ***************)
(******* t est l'angle en degre de la rotation *******)
let m_rot_i t =
let phi = (float_of_int t) *. Pi /. 180. in
[|[| 1. ; 0. ; 0. |];
[| 0. ; cos phi ; -. (sin phi) |];
[| 0. ; sin phi ; cos phi |]|];;
(******* Matrice de rotation / (0,1,0) ***************)
(******* t est l'angle en degre de la rotation *******)
let m_rot_j t =
let phi = (float_of_int t) *. Pi /. 180. in
[|[| cos phi ; 0. ; -. (sin phi) |];
[| 0. ; 1. ; 0. |];
[| sin phi ; 0. ; cos phi |]|];;
(******* Matrice de rotation / (0,0,1) ***************)
(******* t est l'angle en degre de la rotation *******)
let m_rot_k t =
let phi = (float_of_int t) *. Pi /. 180. in
[|[| cos phi ; -. (sin phi) ; 0.|];
[| sin phi ; cos phi ; 0.|];
[| 0. ; 0. ; 1.|]|];;
(******* Produit de 2 matrices (multipliables) *******)
let pmat m1 m2 =
let n1 = vect_length m1 and p1 = vect_length m1.(0)
and n2 = vect_length m2 and p2 = vect_length m2.(0)
in
(if (p1<>n2) then failwith "pas de multiplication");
let m = make_matrix n1 p2 0. in
for i=0 to n1-1 do
for j=0 to p2-1 do
for k=0 to n2-1 do
m.(i).(j) <- m.(i).(j) +. m1.(i).(k) *. m2.(k).(j)
done;
done;
done;
m;;
(******* Multiplication d'une matrice par un scalaire*)
let scalmat l m =
let n = vect_length m and p = vect_length m.(0) in
let m2 =make_matrix n p 0. in
for i = 0 to n-1 do
for j = 0 to p-1 do
m2.(i).(j) <- l *. m.(i).(j)
done
done;
m2;;
(******* Somme de 2 matrices *************************)
let summat m1 m2 =
let n1 = vect_length m1 and p1 = vect_length m1.(0)
and n2 = vect_length m2 and p2 = vect_length m2.(0)
in
(if (p1<>p2) or (n1<>n2) then failwith "addition impossible");
let m2 =make_matrix n1 p1 0. in
for i = 0 to n1-1 do
for j = 0 to p1-1 do
m2.(i).(j) <- m1.(i).(j) +. m2.(i).(j)
done
done;
m2;;
(******* Appliquer une matrice sur un point ***********)
let apply m (x,y,z) =
(m.(0).(0) *. x +. m.(0).(1) *. y +. m.(0).(2) *. z,
m.(1).(0) *. x +. m.(1).(1) *. y +. m.(1).(2) *. z,
m.(2).(0) *. x +. m.(2).(1) *. y +. m.(2).(2) *. z);;
(******* appliquer une rotation sur un polygone *******)
let app_rot r {T = T1 ; Vis = v} =
let T2 = map_vect (apply r) T1
and v2 = apply r v in
{T = T2 ; Vis = v2};;
(******* matrice de rotation tout axes ******************)
let rot_mat i j k = pmat (pmat (m_rot_i i) (m_rot_j j)) (m_rot_k k);;
(*****************************************************)
(*****************************************************)
(*****************************************************)
(****** ajouter les 4 faces d'une pyramides a la *****)
(****** liste des faces a trier en appliquant la *****)
(****** rotation r et en determinant les couleurs ****)
let ajouter_pyram (f1,f2,f3,f4) l r =
let (g1,g2,g3,g4) = (app_rot r f1 , app_rot r f2 ,
app_rot r f3 , app_rot r f4 ) in
let c = FColor g1 in
l:={F = g1 ; Fill_c = c ; Bord_c = black}::(!l);
let c = FColor g2 in
l:={F = g2 ; Fill_c = c ; Bord_c = black}::(!l);
let c = FColor g3 in
l:={F = g3 ; Fill_c = c ; Bord_c = black}::(!l);
let c = FColor g4 in
l:={F = g4 ; Fill_c = c ; Bord_c = black}::(!l);
;;
(****** Trier les faces de la liste l*****************)
let trier_faces l =
let v = vect_of_list l
and cmp f1 f2 = (let (_,_,z) = (barycentre f1.F.T) in z) <.
(let (_,_,z) = (barycentre f2.F.T) in z) in
tri_fusion cmp v;;
let trier_faces l =
let milieu (x1,y1,z1) (x2,y2,z2) = ((x1+.x2)/.2.,(y1+.y2)/.2.,(z1+.z2)/.2.) in
let normal F =
let M = milieu F.T.(0) F.T.(1) in
let T = [| F.T.(0) ; F.T.(1) ; vadd M (normalize (vdiff F.T.(2) M))|] in
let (_,_,z) = barycentre T in z
in let v = vect_of_list l
and cmp f1 f2 = normal f1.F <. normal f2.F in
tri_fusion cmp v;;
let trier_faces l =
let milieu (x1,y1,z1) (x2,y2,z2) = ((x1+.x2)/.2.,(y1+.y2)/.2.,(z1+.z2)/.2.) in
let normalize u = vscal (1. /. norme(u)) u in
let normal F =
let M = milieu F.T.(0) F.T.(1) in
let T = [| F.T.(0) ; F.T.(1) |] in
let (_,_,z) = barycentre T in z
in let v = vect_of_list l
and cmp f1 f2 = normal f1.F <. normal f2.F in
tri_fusion cmp v;;
(****** afficher un tableau de liste *****************)
let afficher_faces v =
let n = vect_length v in
for i = 0 to n-1 do
FcRemplir v.(i).F v.(i).Fill_c v.(i).Bord_c
done;;
(******** "polygones" representant les axes **********)
let axe_i l=
let A = ((l /. 2.),0.,0.)
and B = ((-. l/. 2.),0.,0.) in
{F = { T = [|A ; B|] ; Vis = (0.,0.,0.)} ; Fill_c = black ; Bord_c = black};;
let axe_j l=
let A = (0.,(l /. 2.),0.)
and B = (0.,(-. l/. 2.),0.) in
{F = { T = [|A ; B|] ; Vis = (0.,0.,0.)} ; Fill_c = black ; Bord_c = black};;
let axe_k l=
let A = (0.,0.,(l /. 2.))
and B = (0.,0.,(-. l/. 2.)) in
{F = { T = [|A ; B|] ; Vis = (0.,0.,0.)} ; Fill_c = black ; Bord_c = black};;
(*****************************************************)
(*****************************************************)
(*****************************************************)
(********** matrice symetrique aleatoire d'ordre n****)
(* max ext l'amplitude des valeurs aleatoires ********)
(* M(i,j) E [-max ; max ] ****************************)
let rand_smat n max=
let m = make_matrix n n 0. in
for i = 0 to n-1 do
for j=0 to i do
m.(i).(j) <- (random__float (2.*.max)) -. max;
m.(j).(i) <- m.(i).(j)
done
done;
m;;
(********** matrice aleatoire d'ordre n *************)
(* max ext l'amplitude des valeurs aleatoires ********)
(* M(i,j) E [-max ; max ] ****************************)
let rand_mat n max =
let m = make_matrix n n 0. in
for i = 0 to n-1 do
for j = 0 to n-1 do
m.(i).(j) <- (random__float (2.*.max)) -. max
done
done;
m;;
(************** signe d'un flottant (|a|/a) **********)
let signe a =
if a =. 0. then
1.
else a /. (abs_float a)
;;
(*****************************************************)
(*****************************************************)
(*****************************************************)
(************** dessiner une matrice m ***************)
(*
lx,ly,lz : largeur de la boite dans laquelle est
contenue la matrice.
rot : rotation a lui appliquer
*)
let plotmat m lx ly lz rot=
let n = vect_length m in
let stepx = lx /. (float_of_int n)
and stepz = lz /. (float_of_int n)
and l_faces = ref [] (*ref [axe_i lx;axe_j ly;axe_k lz]*)
and xstart = -.lx /. 2. and zstart = -.lz /. 2. in
let Sxstart = (-.lx +. stepx) /. 2.
and Szstart = (-.lz +. stepz) /. 2. in
for i = 0 to n-1 do
for j = 0 to n-1 do
let (i1,j1) = (float_of_int i , float_of_int j) in
let (i2,j2) = (i1 +. 1. , j1 +. 1.) in
let M = ( xstart +. j1 *. stepx , 0. , zstart +. i1 *. stepz)
and P = ( xstart +. j1 *. stepx , 0. , zstart +. i2 *. stepz)
and Q = ( xstart +. j2 *. stepx , 0. , zstart +. i2 *. stepz)
and R = ( xstart +. j2 *. stepx , 0. , zstart +. i1 *. stepz)
and S = ( Sxstart +. j1 *. stepx, m.(i).(j) , Szstart +. i1 *. stepz)
and coef = signe m.(i).(j) in
(*let B = { T = [|M ; P ; Q ; R|] ; Vis = vscal (-. coef) (TrigoVis M P Q) }*)
let F1 = { T = [|M ; P ; S|] ; Vis = vscal coef (TrigoVis M P S) }
and F2 = { T = [|P ; Q ; S|] ; Vis = vscal coef (TrigoVis P Q S) }
and F3 = { T = [|Q ; R ; S|] ; Vis = vscal coef (TrigoVis Q R S) }
and F4 = { T = [|R ; M ; S|] ; Vis = vscal coef (TrigoVis R M S) } in
ajouter_pyram (F1,F2,F3,F4) l_faces rot;
done
done;
afficher_faces (trier_faces !l_faces)
;;
(*****************************************************)
(*****************************************************)
(*****************************************************)
random__init 100000;;
spot := (-900. , 100. , 900.);;
zobs := 500.;;
xobs := 0.;;
yobs := 0.;;
observateur := (!xobs, !yobs, !zobs);;
pdfz := -.1000.;;
pdf := (0.,500.,!pdfz);;
let tmat = rand_smat 15 100.;;
open_graph "";;
set_color (rgb 0 0 0xFF);;
(*fill_rect 0 0 640 480;;*)
plotmat tmat 400. 100. 400. (rot_mat 30 40 0);;
(*write_image "test1.tga" 0 0 640 480;;*)
wait();;