(* I Recherche de l'enveloppe convexe d'un nuage de points *)
(* II.2 Préliminaires *)
(* Question 1. Fonction aire. On calcule la moitié de la composante du produit vectorielle sur u_z *)
let aire (xi, yi) (xj, yj) (xk, yk) = 
  ((yk -. yi) *. (xj -. xi) -. (xk -. xi ) *. (yj -. yi)) /. 2. ;;

(* Question 2. Fonction orient. Il suffit de calculer l'aire signée par la question précédente et d'en déduire la valeur à rendre. *)
let orient pi pj pk = 
  let  a =  aire pi pj pk in
    if a > 0. then  1 else if a = 0. then 0 else -1 ;;

orient (0., 0.) (0., 1.) (1., 1.) ;;
orient (0., 0.) (0., 1.) (1., -1.) ;;
orient (0., 0.) (0., 1.) (-1., 1.) ;;


(* II.3 Approche diviser pour régner : QuickHull *)
(* Question 3 Idée pour gérer la présence de p et q dans cet ordre en tête de liste : on extrait récursivement les points strictement à gauche par une fonction auxiliaire récursive et on rajoute p et q en tête *)
let nuage_gauche nuage p q =
  let rec nuage_gauche_strict nuage =
    match nuage with
    | [] -> []
    | tete :: queue -> if (orient p q tete) = 1 then tete :: (nuage_gauche_strict queue) else nuage_gauche_strict queue 
  in p :: q :: (nuage_gauche_strict nuage) ;;

let test_l = [(0., 0.); (1., 4.); (1., 8.); (4., 1.); (4., 4.); (5., 9.); (5., 6.); (7., -1.); (7., 2.); (8., 5.); (11., 6.); (13., 1.)];;
let test_l2 = [(7., -1.); (11., 6.); (1., 8.); (4., 1.); (4., 4.); (5., 9.); (5., 6.); (0., 0.); (7., 2.); (8., 5.); (1., 4.); (13., 1.)];;
let test_l3 = [(0., 0.); (13., 1.); (1., 8.); (4., 1.); (4., 4.); (5., 9.); (5., 6.); (7., -1.); (7., 2.); (8., 5.); (11., 6.); (1., 4.)];;
let test_l4 = [(1., 4.); (0., 0.); (1., 8.); (4., 1.); (4., 4.); (5., 9.); (5., 6.); (7., -1.); (7., 2.); (8., 5.); (11., 6.); (13., 1.)];;

nuage_gauche test_l (0., 0.) (13., 1.) ;;  
nuage_gauche test_l (13., 1.) (0., 0.) ;;  

(* Question 4. D'après la formule de l'aire d'un triangle 1/2 * base * hauteur, si on considère tous les triangles de base (pq), leurs aires sont proportionnelles à la hauteur qui est précisément la distance de h à la droite (pq). Par ailleurs comme on travaille sur un ensemble avec des points à gauche de pq toutes les aires sont positives. Trouver le point h le plus éloigné est donc bien déterminer celui pour lequel l'aire signée du triangle pqh est maximale ! *)


(* Question 5. *)
let rec plus_loin nuage =
  match nuage with
  | [] | [_] | [_; _] -> failwith "nuage trop petit"
  | p :: q :: m :: [] -> m
  | p :: q :: m :: queue -> let pl = plus_loin (p :: q:: queue) in if aire p q m > aire p q pl then m else pl ;;


plus_loin test_l ;;
plus_loin test_l2 ;;
plus_loin test_l3 ;;
plus_loin test_l4 ;;

(* Question 6. Dans les deux ensembles ignorés les points sont forcément strictement à l'intérieur de l'enveloppe convexe. On n'a donc pas à les traiter *)

(* Question 7. *)
let rec semihull nuage =
  match nuage with
  | [] | [_] -> failwith "nuage trop petit" 
  | [_; _] -> nuage
  | p :: q :: queue -> let h = plus_loin nuage in
                       let nuage1 = nuage_gauche nuage p h in
		       let nuage2 = nuage_gauche nuage h q in
                       let hull1 = semihull nuage1 and hull2 = semihull nuage2 in
                       hull1 @ List.tl hull2;;

semihull test_l3 ;;

(* Question 8. *)
(* Il faut également une fonction d'extraction des points les plus à gauche et à droite *)
(* La fonction rend un couple (P, Q) avec P le plus à gauche, et Q le plus à droite *)
let rec points_extremes nuage = 
  match nuage with
  | [] | [_] -> failwith "nuage trop petit" (* Je laisse ces motifs pour ne faire râler le compilateur *)
  | [tete1; tete2] -> if (fst tete1) < (fst tete2) then (tete1, tete2) else (tete2, tete1) 
  | tete :: queue -> let (gauche, droite) = points_extremes queue in
                     let gauche' = if (fst tete) < (fst gauche) then tete else gauche in
		     let droite' = if (fst tete) > (fst droite) then tete else droite in
                     (gauche', droite') ;;

points_extremes test_l ;;


(* Question 9. QuickHull qui consiste à chercher le premier découpage à l'aide des points extremes, et à appeler semihull sur les deux moitiés et à concaténer les résultats. On suppose sans vérifier qu'il y a au moins deux points dans le nuage*)
let quickhull nuage =
  let (gauche, droite) = points_extremes nuage in
  let nuage1 = nuage_gauche nuage gauche droite in
  let nuage2 = nuage_gauche nuage droite gauche in
  let hull1 = semihull nuage1 and hull2 = semihull nuage2 in
  List.tl hull1 @ List.tl hull2 ;;
  

quickhull test_l ;;

(* Question 10. L'étape de fusion est en O(n), car il faut d'abord trouver le point h (en temps linéaire) et extraire les deux nouveaux ensembles à gauche de PH et HQ ce qui se fait aussi en temps linéaire. L'étape de fusion (concaténation...) se fait en temps linéaire en la taille courante de la liste de points de l'enveloppe, ce qui est inférieur à du O(n). Dès lors les étapes de division/fusion sont en temps linéaire. Avec les hypothèses de l'énoncé on a alors facilement C(n) = 2 * C(n/2) + O(n),  ce qui amène une complexité quasi linéaire, i.e. O(n * log n), comme le tri fusion *)

(* Question 11 *)
(* Ce qui fait baisser la complexité c'est le fait de découper en problème d'à peu près la même taille. On peut facilement imaginer une situation où après détection du point h il n'y a aucun point autres que p et h pour le premier, et tous les autres points dans l'autre ! Dès lors la relation de récurrence est du genre C(n) = C(n-1) + O(n) (le O(n) vient de l'extraction des points à gauche de hq ), ce qui conduit à une complexité en O(n^2) *)

(* Dans ce qui suit vous trouverez un ensemble de fonctions utilitaires pour créer des tableaux ou des listes de points aléatoires, pour les afficher, et pour dessiner également l'enveloppe *)
(* On essaye maintenant avec un nuage aléatoire et une représentation graphique. les points auront des coordonnées flottantes entre 0 et 400 *)
let liste_aleatoire n = 
  let xmin = 50. and xmax = 350. and ymin = 50. and ymax = 350. in
  let rep = Array.make n (xmin, ymin) in
  for i = 0 to n - 1
  do
    let x = xmin +. Random.float (xmax -. xmin) in
    let y = ymin +. Random.float (ymax -. ymin) in
    rep.(i) <- (x, y) 
  done ;
  Array.to_list rep ;;

(* Fonction de visualisation du tableau de points. Pour fermer la fenetre il suffit d'appuyer sur une touche *)
#load "graphics.cma" ;;
open Graphics ;;

(* Fonction de conversion liste vers tableau *)
let to_array liste =
  let n = List.length liste in
  let tab = Array.make n (List.hd liste) in
  let liste_courante = ref (List.tl liste) in
  for i = 1 to n - 1 do
    tab.(i) <- List.hd !liste_courante ;
    liste_courante := List.tl !liste_courante
  done ;
  tab ;;

  
let affiche liste =
  let tableau = to_array liste in
  let n = Array.length tableau in
  open_graph " 400x400" ;
  
  for i = 0 to n - 1
  do
    draw_circle (int_of_float (fst tableau.(i))) (int_of_float (snd tableau.(i))) 3
  done ;
  let _ = read_key () in () ;
  close_graph () ;;

let nuage_test = liste_aleatoire 100 ;;
let env_dpr = quickhull nuage_test ;;

affiche nuage_test ;;
(* Il faut enfin une fonction graphique pour tracer l'enveloppe *)

let trace nuage enveloppe =
  let tableau = to_array nuage in
  let n = Array.length tableau in
  open_graph " 400x400" ;
  
  for i = 0 to n - 1
  do
    draw_circle (int_of_float (fst tableau.(i))) (int_of_float (snd tableau.(i))) 3
  done ;
  let tab_env = to_array enveloppe in
  let n = Array.length tab_env in
  moveto (int_of_float (fst tab_env.(0))) (int_of_float (snd tab_env.(0))) ;
  for i = 1 to n - 1 do
    lineto (int_of_float (fst tab_env.(i))) (int_of_float (snd tab_env.(i)))
  done ;
  lineto (int_of_float (fst tab_env.(0))) (int_of_float (snd tab_env.(0))) ;
  let _ = read_key () in () ;
  close_graph () ;;

trace nuage_test env_dpr;;