Update project structure

38 files changed, 2967 insertions, 0 deletions

diff --git a/script.it/layer/canvaPrinter.ml b/script.it/layer/canvaPrinter.ml
new file mode 100755
index 0000000..23cf842
--- /dev/null
+++ b/script.it/layer/canvaPrinter.ml
@@ -0,0 +1,42 @@
+module Path = Brr_canvas.C2d.Path
+module V2 = Gg.V2
+
+type t = Path.t
+
+let create
+  : unit -> t
+  = Path.create
+
+(* Start a new path. *)
+let move_to
+  : Gg.v2 -> t -> t
+  = fun point path ->
+    let x, y = V2.to_tuple point in
+    Path.move_to ~x ~y path;
+    path
+
+let line_to
+  : Gg.v2 -> t -> t
+  = fun  point path ->
+    let x, y = V2.to_tuple point in
+    Path.line_to ~x ~y path;
+    path
+
+let quadratic_to
+  : Gg.v2 -> Gg.v2 -> Gg.v2 -> t -> t
+  = fun  ctrl0 ctrl1 p1 path ->
+    let cx, cy = V2.to_tuple ctrl0
+    and cx', cy' = V2.to_tuple ctrl1
+    and x, y = V2.to_tuple p1 in
+    Path.ccurve_to
+      ~cx ~cy
+      ~cx' ~cy'
+      ~x ~y
+      path;
+    path
+
+let close
+  : t -> t
+  = fun path ->
+    Path.close path;
+    path
diff --git a/script.it/layer/canvaPrinter.mli b/script.it/layer/canvaPrinter.mli
new file mode 100755
index 0000000..0c46448
--- /dev/null
+++ b/script.it/layer/canvaPrinter.mli
@@ -0,0 +1,2 @@
+include Repr.PRINTER 
+  with type t = Brr_canvas.C2d.Path.t
diff --git a/script.it/layer/ductusEngine.ml b/script.it/layer/ductusEngine.ml
new file mode 100755
index 0000000..b943467
--- /dev/null
+++ b/script.it/layer/ductusEngine.ml
@@ -0,0 +1,82 @@
+module Make(Layer: Repr.PRINTER) = struct
+
+  type point = Path.Point.t
+
+  type t =
+    { path: (Layer.t)
+    }
+
+  type repr = Layer.t
+
+  let create_path
+    : 'b -> t
+    = fun _ ->
+      { path = Layer.create ()
+      }
+
+  let start
+    : point -> point -> t -> t
+    = fun p1 p2 { path } ->
+      let path =
+        Layer.move_to (Path.Point.get_coord p1) path
+        |> Layer.line_to (Path.Point.get_coord p2) in
+      { path
+      }
+
+  let line_to
+    : (point * point) -> (point * point) -> t -> t
+    = fun (_, p1) (_, p1') {path} ->
+      let path = Layer.move_to (Path.Point.get_coord p1) path in
+      let path = Layer.line_to (Path.Point.get_coord p1') path in
+      { path
+      }
+
+  let quadratic_to
+    : (point * Gg.v2 * Gg.v2 * point) -> (point * Gg.v2 * Gg.v2 * point) -> t -> t
+    = fun (p0, ctrl0, ctrl1, p1) (p0',  ctrl0', ctrl1', p1') {path} ->
+
+      let bezier =
+        { Shapes.Bezier.p0 = Path.Point.get_coord p0
+        ; ctrl0
+        ; ctrl1
+        ; p1 = Path.Point.get_coord p1
+        }
+
+      and bezier' =
+        { Shapes.Bezier.p0 = Path.Point.get_coord p0'
+        ; ctrl0 = ctrl0'
+        ; ctrl1 = ctrl1'
+        ; p1 = Path.Point.get_coord p1'
+        } in
+
+      (* Mark each point on the bezier curve. The first point is the most
+         recent point  *)
+      let delay =
+        ((Path.Point.get_stamp p0) -. (Path.Point.get_stamp p1))
+        *. 30.
+      in
+      let path = ref path in
+      for i = 0 to ((Int.of_float delay) ) do
+        let ratio = (Float.of_int i) /. delay in
+        let bezier, _ = Shapes.Bezier.slice ratio bezier
+        and bezier', _ = Shapes.Bezier.slice ratio bezier'  in
+
+        let point = Path.Point.mix ratio bezier.Shapes.Bezier.p1 p0 p1
+        and point' = Path.Point.mix ratio bezier'.Shapes.Bezier.p1 p0' p1' in
+
+        path := Layer.move_to (Path.Point.get_coord point) !path;
+        path := Layer.line_to (Path.Point.get_coord point') !path;
+      done;
+
+      { path = !path }
+
+  let stop
+    : t -> t
+    = fun path -> path
+
+
+  let get
+    : t -> Layer.t
+    = fun {path; _} ->
+      path
+end
diff --git a/script.it/layer/ductusEngine.mli b/script.it/layer/ductusEngine.mli
new file mode 100755
index 0000000..e1660f4
--- /dev/null
+++ b/script.it/layer/ductusEngine.mli
@@ -0,0 +1,2 @@
+module Make(R:Repr.PRINTER): 
+  Repr.ENGINE with type repr = R.t
diff --git a/script.it/layer/dune b/script.it/layer/dune
new file mode 100755
index 0000000..3c617ad
--- /dev/null
+++ b/script.it/layer/dune
@@ -0,0 +1,8 @@
+(library
+ (name layer)
+ (libraries 
+   gg
+   brr
+   path
+   )
+ )
diff --git a/script.it/layer/fillEngine.ml b/script.it/layer/fillEngine.ml
new file mode 100755
index 0000000..9a3fe7e
--- /dev/null
+++ b/script.it/layer/fillEngine.ml
@@ -0,0 +1,89 @@
+module Make(Layer: Repr.PRINTER) = struct
+
+  type point = Path.Point.t
+
+  type repr = Layer.t
+
+  type t =
+    { path: Layer.t
+    ; close : Layer.t -> Layer.t
+    }
+
+  let create_path
+    : (Layer.t -> Layer.t) -> t
+    = fun f ->
+      { close = f
+      ; path = Layer.create ()
+      }
+
+  (* Start a new path. *)
+
+  let start
+    : point -> point -> t -> t
+    = fun p1 _ {close ; path } ->
+      let path = Layer.move_to (Path.Point.get_coord p1) path in
+      { close
+      ; path
+      }
+
+  let line_to
+    : (point * point) -> (point * point) -> t -> t
+    = fun (p0, p1) (p0', p1') t ->
+
+      let p0 = Path.Point.get_coord p0
+      and p1 = Path.Point.get_coord p1
+      and p0' = Path.Point.get_coord p0'
+      and p1' = Path.Point.get_coord p1' in
+
+      let path =
+        Layer.move_to p1 t.path
+        |> Layer.line_to p1'
+        |> Layer.line_to p0'
+        |> Layer.line_to p0
+        |> Layer.line_to p1
+        |> Layer.close in
+      let path = t.close path in
+      { t with path}
+
+  let quadratic_to
+    : (point * Gg.v2 * Gg.v2 * point) -> (point * Gg.v2 * Gg.v2 * point) -> t -> t
+    = fun (p0, ctrl0, ctrl1, p1) (p0',  ctrl0', ctrl1', p1') t ->
+
+      let p0 = Path.Point.get_coord p0
+      and p1 = Path.Point.get_coord p1
+      and p0' = Path.Point.get_coord p0'
+      and p1' = Path.Point.get_coord p1'
+      in
+
+      let path =
+        Layer.move_to p1 t.path
+        |> Layer.line_to p1'
+
+        (* Backward *)
+        |> Layer.quadratic_to
+          ctrl1'
+          ctrl0'
+          p0'
+        |> Layer.line_to p0
+
+        (* Forward *)
+        |> Layer.quadratic_to
+          ctrl0
+          ctrl1
+          p1
+        |> Layer.close
+        |> t.close in
+
+
+      { t with path }
+
+  let stop
+    : t -> t
+    = fun t ->
+      t
+
+  let get
+    : t -> Layer.t
+    = fun t ->
+      t.path
+end
diff --git a/script.it/layer/fillEngine.mli b/script.it/layer/fillEngine.mli
new file mode 100755
index 0000000..e1660f4
--- /dev/null
+++ b/script.it/layer/fillEngine.mli
@@ -0,0 +1,2 @@
+module Make(R:Repr.PRINTER): 
+  Repr.ENGINE with type repr = R.t
diff --git a/script.it/layer/lineEngine.ml b/script.it/layer/lineEngine.ml
new file mode 100755
index 0000000..3d15d9c
--- /dev/null
+++ b/script.it/layer/lineEngine.ml
@@ -0,0 +1,68 @@
+module Make(Layer: Repr.PRINTER) = struct
+
+  type point = Path.Point.t
+
+  let mark point path =
+    let open Gg.V2 in
+    let point = Path.Point.get_coord point in
+
+    let dist = 5.
+    and dist' = -5. in
+
+    let path = Layer.move_to (point - (of_tuple (dist, dist))) path
+               |> Layer.line_to ( point + (of_tuple (dist, dist)))
+               |> Layer.move_to (point + (of_tuple (dist', dist)))
+               |> Layer.line_to ( point + (of_tuple (dist, dist')))
+    in
+    path
+
+
+  type t =
+    { path: (Layer.t)
+    }
+
+  type repr = Layer.t
+
+  let create_path
+    : 'b -> t
+    = fun _ ->
+      { path = Layer.create ()
+      }
+
+  let start
+    : point -> point -> t -> t
+    = fun p1 _ { path } ->
+      let path = mark p1 path in
+      { path
+      }
+
+  let line_to
+    : (point * point) -> (point * point) -> t -> t
+    = fun (p0, p1) _ {path} ->
+      let path = Layer.move_to (Path.Point.get_coord p0) path
+                 |> Layer.line_to (Path.Point.get_coord p1)
+                 |> mark p1 in
+      { path
+      }
+
+  let quadratic_to
+    : (point * Gg.v2 * Gg.v2 * point) -> (point * Gg.v2 * Gg.v2 * point) -> t -> t
+    = fun (p0, ctrl0, ctrl1, p1) _ {path} ->
+
+      let path = Layer.move_to (Path.Point.get_coord p0) path
+                 |> Layer.quadratic_to ctrl0 ctrl1 (Path.Point.get_coord p1)
+                 |> mark p1 in
+
+      { path = path }
+
+  let stop
+    : t -> t
+    = fun path -> path
+
+
+  let get
+    : t -> Layer.t
+    = fun {path; _} ->
+      path
+
+end
diff --git a/script.it/layer/lineEngine.mli b/script.it/layer/lineEngine.mli
new file mode 100755
index 0000000..86ef5fb
--- /dev/null
+++ b/script.it/layer/lineEngine.mli
@@ -0,0 +1,2 @@
+module Make(R:Repr.PRINTER):
+  Repr.ENGINE with type repr = R.t
diff --git a/script.it/layer/paths.ml b/script.it/layer/paths.ml
new file mode 100755
index 0000000..d3baf02
--- /dev/null
+++ b/script.it/layer/paths.ml
@@ -0,0 +1,244 @@
+open StdLabels
+(** Common  module for ensuring that the function is evaluated only once *)
+
+(** This represent a single path, which can be transformed throug a [repr]
+    function. *)
+module type PATH = sig
+  type t
+
+  (** Represent the path *)
+  val repr
+    : t -> (module Path.Repr.M
+             with type point = Path.Point.t
+              and type t = 's) -> 's -> 's
+end
+
+type printer =
+  [ `Fill
+  | `Line
+  | `Ductus ]
+
+
+module type P = sig
+  include Path.Repr.M
+
+  type repr
+
+  val create_path
+    : (repr -> repr) -> t
+
+  val get
+    : t -> repr
+end
+
+
+module MakePrinter(M:Repr.ENGINE) : P
+  with type point = M.point
+   and type t = M.t
+   and type repr = M.repr = struct
+
+  type t = M.t
+
+  type point = M.point
+
+  type repr = M.repr
+
+  let get
+    : t -> repr
+    = M.get
+
+  let create_path
+    : (repr -> repr) -> t
+    = M.create_path
+
+  let start
+    : Path.Point.t -> t -> t
+    = fun pt t ->
+      M.start pt pt t
+
+  let line_to
+    : Path.Point.t -> Path.Point.t -> t -> t
+    = fun p0 p1 t ->
+
+      M.line_to
+        ( p0
+        , p1 )
+        ( Path.Point.copy p0 @@ Path.Point.get_coord' p0
+        , Path.Point.copy p1 @@ Path.Point.get_coord' p1 )
+        t
+
+  let quadratic_to
+    : (Path.Point.t * Gg.v2 * Gg.v2 * Path.Point.t) -> t -> t
+    = fun (p0,  ctrl0, ctrl1, p1) t ->
+
+
+      let ctrl0' = Path.Point.get_coord' @@ Path.Point.copy p0 ctrl0
+      and ctrl1' = Path.Point.get_coord' @@ Path.Point.copy p1 ctrl1 in
+      M.quadratic_to
+        (p0, ctrl0, ctrl1, p1)
+        (Path.Point.copy p0 @@ Path.Point.get_coord' p0, ctrl0', ctrl1', Path.Point.copy p1 @@ Path.Point.get_coord' p1)
+
+        t
+
+  let stop = M.stop
+end
+
+(** Transform the two path, into a single one. *)
+module ReprSingle = struct
+
+  type point = Path.Point.t
+
+  type repr =
+    | Move of (point)
+    | Line_to of (point * point)
+    | Quadratic of (point * Gg.v2 * Gg.v2 * point)
+
+  module R = struct
+    type t = repr list
+
+    type point = Path.Point.t
+
+    let start t actions =
+      (Move t)::actions
+
+    let line_to p0 p1 actions =
+      Line_to (p0, p1)::actions
+
+    let quadratic_to
+      : (point * Gg.v2 * Gg.v2 * point) -> t -> t
+      = fun q actions ->
+        (Quadratic q)::actions
+
+    let stop
+      : t -> t
+      = fun v -> v
+
+  end
+
+  let repr
+    : (module PATH with type t = 't) -> 't -> 't -> repr list * repr list
+    = fun (type t) (module P:PATH with type t = t) path back ->
+      let path = P.repr path (module R) []
+      and back = P.repr back (module R) [] in
+      path, back
+end
+
+(* Canva representation *)
+
+module FillCanva = FillEngine.Make(CanvaPrinter)
+module LineCanva = LineEngine.Make(CanvaPrinter)
+module DuctusCanva = DuctusEngine.Make(CanvaPrinter)
+
+(* SVG representation *)
+
+module FillSVG = FillEngine.Make(Svg)
+module DuctusSVG = DuctusEngine.Make(Svg)
+
+
+(** Draw a path to a canva.contents
+
+    The code may seems scary, but is very repetitive. Firt, all points (from the
+    main stroke, and the interior one) are evaluated. Then, they are both rendered
+    using the selected engine.
+*)
+let to_canva
+  : (module PATH with type t = 's) -> 's * 's -> Brr_canvas.C2d.t -> printer -> unit
+  = fun (type s) (module R:PATH with type t = s) (path, back) ctx engine ->
+    let f, b = ReprSingle.repr (module R) path back in
+    match engine with
+    | `Fill ->
+      let t = List.fold_left2 f b
+          ~init:(FillCanva.create_path (fun p -> Brr_canvas.C2d.fill ctx p; p))
+          ~f:(fun ctx f b ->
+              match (f, b) with
+              | ReprSingle.Move p0, ReprSingle.Move p0' -> FillCanva.start p0 p0' ctx
+              | ReprSingle.Line_to l, ReprSingle.Line_to l' -> FillCanva.line_to l l' ctx
+              | ReprSingle.Quadratic q, ReprSingle.Quadratic q' -> FillCanva.quadratic_to q q' ctx
+              | _ -> ctx
+            ) in
+      FillCanva.get t
+      |> Brr_canvas.C2d.stroke ctx
+    | `Line ->
+      let t = List.fold_left2 f b
+          ~init:(LineCanva.create_path (fun p -> Brr_canvas.C2d.fill ctx p; p))
+          ~f:(fun ctx f b ->
+              match (f, b) with
+              | ReprSingle.Move p0, ReprSingle.Move p0' -> LineCanva.start p0 p0' ctx
+              | ReprSingle.Line_to l, ReprSingle.Line_to l' -> LineCanva.line_to l l' ctx
+              | ReprSingle.Quadratic q, ReprSingle.Quadratic q' -> LineCanva.quadratic_to q q' ctx
+              | _ -> ctx
+            ) in
+      LineCanva.get t
+      |> Brr_canvas.C2d.stroke ctx
+    | `Ductus ->
+      let t = List.fold_left2 f b
+          ~init:(DuctusCanva.create_path (fun p -> Brr_canvas.C2d.fill ctx p; p))
+          ~f:(fun ctx f b ->
+              match (f, b) with
+              | ReprSingle.Move p0, ReprSingle.Move p0' -> DuctusCanva.start p0 p0' ctx
+              | ReprSingle.Line_to l, ReprSingle.Line_to l' -> DuctusCanva.line_to l l' ctx
+              | ReprSingle.Quadratic q, ReprSingle.Quadratic q' -> DuctusCanva.quadratic_to q q' ctx
+              | _ -> ctx
+            ) in
+      DuctusCanva.get t
+      |> Brr_canvas.C2d.stroke ctx
+
+
+
+(** Draw a path and represent it as SVG *)
+let to_svg
+  : (module PATH with type t = 's) -> color:Jstr.t -> 's * 's -> printer -> Brr.El.t
+  = fun (type s) (module R:PATH with type t = s) ~color (path, back) engine ->
+    let f, b = ReprSingle.repr (module R) path back in
+    match engine with
+    | `Fill ->
+
+      (* In order to deal with over crossing path, I cut the path in as
+         many segment as there is curve, and fill them all. Then, all of theme
+         are grouped inside a single element *)
+      let paths = ref [] in
+      let init = (FillSVG.create_path
+                    (fun p ->
+                       let repr = Svg.path
+                           ~at:Brr.At.[ v (Jstr.v "d") p ]
+                           [] in
+
+                       paths := repr::!paths;
+                       Jstr.empty)) in
+      let _ = List.fold_left2 f b
+          ~init
+          ~f:(fun ctx f b ->
+              match (f, b) with
+              | ReprSingle.Move p0, ReprSingle.Move p0' -> FillSVG.start p0 p0' ctx
+              | ReprSingle.Line_to l, ReprSingle.Line_to l' -> FillSVG.line_to l l' ctx
+              | ReprSingle.Quadratic q, ReprSingle.Quadratic q' -> FillSVG.quadratic_to q q' ctx
+              | _ -> ctx
+            ) in
+
+      Brr.El.v (Jstr.v "g")
+        ~at:Brr.At.[
+            v (Jstr.v "fill")  color
+          ; v (Jstr.v "stroke") color]
+        !paths
+
+    | `Ductus ->
+      let init = DuctusSVG.create_path (fun _ -> Jstr.empty) in
+      let svg_path = List.fold_left2 f b
+          ~init
+          ~f:(fun ctx f b ->
+              match (f, b) with
+              | ReprSingle.Move p0, ReprSingle.Move p0' -> DuctusSVG.start p0 p0' ctx
+              | ReprSingle.Line_to l, ReprSingle.Line_to l' -> DuctusSVG.line_to l l' ctx
+              | ReprSingle.Quadratic q, ReprSingle.Quadratic q' -> DuctusSVG.quadratic_to q q' ctx
+              | _ -> ctx
+            )
+                     |> DuctusSVG.get in
+
+      Svg.path
+        ~at:Brr.At.[
+            v (Jstr.v "fill")  color
+          ; v (Jstr.v "stroke") color
+          ; v (Jstr.v "d") svg_path ]
+        []
+    | `Line ->
+      raise Not_found
diff --git a/script.it/layer/repr.ml b/script.it/layer/repr.ml
new file mode 100755
index 0000000..552e2b7
--- /dev/null
+++ b/script.it/layer/repr.ml
@@ -0,0 +1,49 @@
+module type PRINTER = sig
+
+  type t
+
+  val create: unit -> t
+
+  (* Start a new path. *)
+  val move_to: Gg.v2 -> t -> t
+
+  val line_to: Gg.v2 -> t -> t
+
+  (** [quadratic_to ctrl0 ctrl1 p1] create a quadratic curve from the current
+      point to [p1], with control points [ctrl0] and [ctrl1] *)
+  val quadratic_to: Gg.v2 -> Gg.v2 -> Gg.v2 -> t -> t
+
+  (** Request for the path to be closed *)
+  val close: t -> t
+
+end
+
+module type ENGINE = sig
+  type t
+
+  type point = Path.Point.t
+
+  type repr
+
+  val get
+    : t -> repr
+
+  val create_path
+    : (repr -> repr) -> t
+
+  val start
+    : point -> point -> t -> t
+
+  val line_to
+    : (point * point) -> (point * point) -> t -> t
+
+  val quadratic_to
+    : (point * Gg.v2 * Gg.v2 * point)
+    -> (point * Gg.v2 * Gg.v2 * point)
+    -> t
+    -> t
+
+  val stop
+    : t -> t
+
+end
diff --git a/script.it/layer/svg.ml b/script.it/layer/svg.ml
new file mode 100755
index 0000000..2394cb8
--- /dev/null
+++ b/script.it/layer/svg.ml
@@ -0,0 +1,64 @@
+(** SVG representation *)
+
+open Brr
+
+module V2 = Gg.V2
+
+let svg : El.cons
+  = fun ?d ?at childs ->
+    El.v ?d ?at (Jstr.v "svg") childs
+
+let path: El.cons
+  = fun ?d ?at childs ->
+    El.v ?d ?at (Jstr.v "path") childs
+
+type t = Jstr.t
+
+let create
+  : unit -> t
+  = fun () -> Jstr.empty
+
+(* Start a new path. *)
+let move_to
+  : Gg.v2 -> t -> t
+  = fun point path ->
+    let x, y = V2.to_tuple point in
+
+    Jstr.concat ~sep:(Jstr.v " ")
+      [ path
+      ; Jstr.v "M"
+      ; Jstr.of_float x
+      ; Jstr.of_float y ]
+
+let line_to
+  : Gg.v2 -> t -> t
+  = fun  point path ->
+    let x, y = V2.to_tuple point in
+    Jstr.concat ~sep:(Jstr.v " ")
+      [ path
+      ; (Jstr.v "L")
+      ; (Jstr.of_float x)
+      ; (Jstr.of_float y) ]
+
+let quadratic_to
+  : Gg.v2 -> Gg.v2 -> Gg.v2 -> t -> t
+  = fun  ctrl0 ctrl1 p1 path ->
+    let cx, cy = V2.to_tuple ctrl0
+    and cx', cy' = V2.to_tuple ctrl1
+    and x, y = V2.to_tuple p1 in
+    Jstr.concat ~sep:(Jstr.v " ")
+      [ path
+      ; (Jstr.v "C")
+      ; (Jstr.of_float cx)
+      ; (Jstr.of_float cy)
+      ; (Jstr.v ",")
+      ; (Jstr.of_float cx')
+      ; (Jstr.of_float cy')
+      ; (Jstr.v ",")
+      ; (Jstr.of_float x)
+      ; (Jstr.of_float y) ]
+
+let close
+  : t -> t
+  = fun path ->
+    Jstr.append path (Jstr.v " Z")
diff --git a/script.it/layer/wireFramePrinter.ml b/script.it/layer/wireFramePrinter.ml
new file mode 100755
index 0000000..81ab271
--- /dev/null
+++ b/script.it/layer/wireFramePrinter.ml
@@ -0,0 +1,80 @@
+module Point = Path.Point
+
+module Make(Repr: Repr.PRINTER) = struct
+  type t = Point.t
+
+  type repr =
+    { back: (Repr.t -> Repr.t)
+    ; path: (Repr.t)
+    ; last_point : Point.t option
+    }
+
+  let create_path
+    : 'b -> repr
+    = fun _ ->
+      { back = Repr.close
+      ; path = Repr.create ()
+      ; last_point = None
+      }
+
+  (* Start a new path. *)
+  let start
+    : Point.t -> repr -> repr
+    = fun t {back; path; _} ->
+      let path = Repr.move_to (Point.get_coord t) path in
+      let line' = Repr.line_to (Point.get_coord' t) in
+      { back = (fun p -> back @@ line' p)
+      ; path
+      ; last_point = Some t
+      }
+
+  let line_to
+    : Point.t -> Point.t -> repr -> repr
+    = fun _ t {back; path; _} ->
+      let line' = Repr.line_to (Point.get_coord' t) in
+      { back = (fun t -> back @@ line' t)
+      ; path = Repr.line_to (Point.get_coord t) path
+      ; last_point = Some t
+      }
+
+  let quadratic_to
+    : Point.t -> Gg.v2 -> Gg.v2 -> Point.t -> repr -> repr
+    = fun p0 ctrl0 ctrl1 p1 t ->
+
+      let ctrl0' = Point.copy p1 ctrl0
+      and ctrl1' = Point.copy p1 ctrl1 in
+
+      let line' path =
+        Repr.quadratic_to
+          (Point.get_coord' @@ ctrl1')
+          (Point.get_coord' ctrl0')
+          (Point.get_coord' p0) path in
+
+      let path = Repr.quadratic_to
+          (Point.get_coord ctrl0')
+          (Point.get_coord ctrl1')
+          (Point.get_coord p1)
+          t.path in
+      { back = (fun p -> t.back @@ line' p)
+      ; path
+      ; last_point = Some p1
+      }
+
+  let stop
+    : repr -> repr
+    = fun {back; path; last_point} ->
+
+      let path =
+        match last_point with
+        | Some point -> Repr.line_to (Point.get_coord' point) path
+        | None -> path in
+
+      { back = (fun x -> x)
+      ; path = back path
+      ; last_point = None }
+
+  let get
+    : repr -> Repr.t
+    = fun {back; path; _} ->
+      back path
+end
diff --git a/script.it/layer/wireFramePrinter.mli b/script.it/layer/wireFramePrinter.mli
new file mode 100755
index 0000000..b198d58
--- /dev/null
+++ b/script.it/layer/wireFramePrinter.mli
@@ -0,0 +1,27 @@
+module Make(Repr:Repr.PRINTER): sig
+
+  type repr 
+
+  type t = Path.Point.t
+
+  val create_path
+    : 'b -> repr
+
+  (* Start a new path. *)
+  val start
+    : Path.Point.t -> repr -> repr
+
+  val line_to
+    : Path.Point.t -> Path.Point.t -> repr -> repr
+
+  val quadratic_to
+    : Path.Point.t -> Gg.v2 -> Gg.v2 -> Path.Point.t -> repr -> repr
+
+  val stop
+    : repr -> repr
+
+
+  val get 
+    : repr -> Repr.t
+
+end
diff --git a/script.it/path/builder.ml b/script.it/path/builder.ml
new file mode 100755
index 0000000..4403599
--- /dev/null
+++ b/script.it/path/builder.ml
@@ -0,0 +1,224 @@
+open StdLabels
+
+(** Signature for points *)
+module type P = sig
+  type t
+
+  val empty : t
+
+  val get_coord : t -> Gg.v2
+
+  val copy : t -> Gg.v2 -> t
+
+end
+
+module Make(Point:P) = struct
+
+  (** Point creation  **)
+
+  type bezier =
+    { p0:Point.t    (* The starting point *)
+    ; p1:Point.t    (* The end point *)
+    ; ctrl0:Gg.v2   (* The control point *)
+    ; ctrl1:Gg.v2 } (* The control point *)
+
+  type t = Point.t list * bezier list
+
+  let get_new_segment connexion0 p5 p4 p3 p2 p1 =
+    let p5' = Point.get_coord p5
+    and p4' = Point.get_coord p4
+    and p3' = Point.get_coord p3
+    and p2' = Point.get_coord p2
+    and p1' = Point.get_coord p1 in
+
+    let points_to_link =
+      [ p1'
+      ; p2'
+      ; p3'
+      ; p4'
+      ; p5' ] in
+    Shapes.Bspline.to_bezier ?connexion0 points_to_link
+
+  let empty = ([], [])
+
+  let add_point
+    : Point.t -> t -> t
+    = fun lastPoint (path, beziers) ->
+      let (let*) v f =
+        match v with
+        | Ok bezier ->
+          if Array.length bezier > 0 then
+            f (Array.get bezier 0)
+          else
+            (lastPoint::path, beziers)
+        | _ ->
+          (lastPoint::path, beziers)
+      in
+
+      let connexion0 = match beziers with
+        | hd::_ -> Some (Point.get_coord hd.p1)
+        | _ -> None in
+
+      match path with
+      | p4::p3::p2::p1::_ ->
+        let* bezier = get_new_segment connexion0
+            lastPoint p4 p3 p2 p1 in
+
+        let bezier_point =
+          { p0 = p2
+          ; p1 = p1
+          ; ctrl0 = bezier.Shapes.Bezier.ctrl1
+          ; ctrl1 = bezier.Shapes.Bezier.ctrl0
+          } in
+
+        (* We remove the last point and add the bezier curve in the list*)
+        let firsts = lastPoint::p4::p3::p2::[] in
+        (firsts, bezier_point::beziers)
+      | _ ->
+        ( lastPoint::path
+        , beziers)
+
+  let replace_last
+    : Point.t -> t -> t
+    = fun lastPoint ((path, beziers) as t) ->
+      match path, beziers with
+      | _::(tl), beziers ->
+
+        ( lastPoint::tl
+        , beziers )
+      | _ ->
+        add_point lastPoint t
+
+  let peek2
+    : t -> (Point.t * Point.t) option
+    = fun (path, _) ->
+      match path with
+      | h1::h2::_ -> Some (h1, h2)
+      | _ -> None
+
+  let peek
+    : t -> Point.t option
+    = fun (path, _) ->
+      match path with
+      | [] -> None
+      | hd::_ -> Some hd
+
+  let repr
+    : t -> (module Repr.M with type point = Point.t and type t = 's) -> 's -> 's
+    = fun (type s) (points, beziers) (module Repr : Repr.M with type point = Point.t and type t = s) path ->
+
+      (* Represent the last points *)
+      let path = match points with
+        | [] ->
+          ( path )
+        | p1::[] ->
+          ( Repr.start p1 path )
+        | p1::p2::[] ->
+          let path =
+            Repr.start p1 path
+            |> Repr.line_to p1 p2 in
+          ( path )
+        | p0::p1::p2::[] ->
+
+          let b0, b1 = Shapes.Bezier.quadratic_to_cubic
+            @@ Shapes.Bezier.three_points_quadratic
+              (Point.get_coord p0)
+              (Point.get_coord p1)
+              (Point.get_coord p2)
+                       |> Shapes.Bezier.slice 0.5
+          in
+          let p0' = Point.copy p0 b0.Shapes.Bezier.p0
+          and p1' = Point.copy p1 b0.Shapes.Bezier.p1
+          and p2' = Point.copy p2 b1.Shapes.Bezier.p1 in
+
+          Repr.start p0 path
+          |> Repr.quadratic_to
+            ( p0'
+            , b0.Shapes.Bezier.ctrl0
+            , b0.Shapes.Bezier.ctrl1
+            , p1' )
+          |> Repr.quadratic_to
+            ( p1'
+            , b1.Shapes.Bezier.ctrl0
+            , b1.Shapes.Bezier.ctrl1
+            , p2' )
+        | (p0::_ as points) ->
+
+          let (let*) v f =
+            match v with
+            | Ok beziers -> f beziers
+            | _ -> path in
+
+          let points' = List.map ~f:Point.get_coord points in
+          let connexion = match beziers with
+            | [] -> None
+            | hd ::_ -> Some (Point.get_coord hd.p1) in
+
+          let* beziers = Shapes.Bspline.to_bezier ?connexion1:connexion points' in
+
+          (* Stdlib does not provide fold_left_i function and we need to map
+             each bezier point with the associated point in the curve.
+
+             So I use references here for keeping each result element
+
+          *)
+          let path = ref path in
+          let point = ref p0 in
+
+          List.iteri
+            points
+            ~f:(fun i pt ->
+
+                (* The first iteration is ignored, as we need both previous and
+                   current point for the two point in the curve.
+
+                   Do not forget that there is always n-1 bezier curve for n
+                   points *)
+                if i > 0 then (
+
+                  let bezier = Array.get beziers (i - 1) in
+
+                  path := Repr.quadratic_to
+                      ( !point
+                      , bezier.Shapes.Bezier.ctrl0
+                      , bezier.Shapes.Bezier.ctrl1
+                      , pt )
+                      (!path);
+                  point := pt;
+                )
+              );
+          ( !path )
+      in
+
+      (* Now represent the already evaluated points. Much easer to do, just
+         iterate on them *)
+      Repr.stop @@ List.fold_left beziers
+        ~init:path
+        ~f:(fun path bezier ->
+            Repr.quadratic_to
+              ( bezier.p0
+              , bezier.ctrl0
+              , bezier.ctrl1
+              , bezier.p1 )
+              path
+          )
+
+  let map
+    : t -> (Point.t -> Point.t) -> t
+    = fun (points, beziers) f ->
+      let points = List.map
+          points
+          ~f
+      and beziers = List.map
+          beziers
+          ~f:(fun bezier ->
+
+              { p0 = f bezier.p0
+              ; p1 = f bezier.p1
+              ; ctrl0 = Point.(get_coord (f ( copy bezier.p0 bezier.ctrl0)))
+              ; ctrl1 = Point.(get_coord (f ( copy bezier.p1 bezier.ctrl1)))
+              }
+            ) in
+      points, beziers
+
+end
diff --git a/script.it/path/builder.mli b/script.it/path/builder.mli
new file mode 100755
index 0000000..2afbd4b
--- /dev/null
+++ b/script.it/path/builder.mli
@@ -0,0 +1,43 @@
+(** Signature for points *)
+module type P = sig
+  type t
+
+  val empty : t
+
+  val get_coord : t -> Gg.v2
+
+  (** Copy a point and all thoses properties to the given location *)
+  val copy : t -> Gg.v2 -> t
+
+end
+
+module Make(Point:P) : sig
+
+  type t
+
+  (** Create an empty path *)
+  val empty: t
+
+  val add_point
+    : Point.t -> t -> t
+
+  (** Replace the last alement in the path by the one given in parameter *)
+  val replace_last
+    : Point.t -> t -> t
+
+  (** Retrieve the last element, if any *)
+  val peek 
+    : t -> Point.t option
+
+  (** Retrieve the last element, if any *)
+  val peek2
+    : t -> (Point.t * Point.t) option
+
+  (** Represent the path *)
+  val repr
+    : t -> (module Repr.M with type point = Point.t and type t = 's) -> 's -> 's
+
+  val map 
+    : t -> (Point.t -> Point.t) -> t
+
+end
diff --git a/script.it/path/dune b/script.it/path/dune
new file mode 100755
index 0000000..863c768
--- /dev/null
+++ b/script.it/path/dune
@@ -0,0 +1,7 @@
+(library
+ (name path)
+ (libraries 
+   gg
+   shapes
+   )
+ )
diff --git a/script.it/path/fixed.ml b/script.it/path/fixed.ml
new file mode 100755
index 0000000..1362ad3
--- /dev/null
+++ b/script.it/path/fixed.ml
@@ -0,0 +1,487 @@
+open StdLabels
+
+(** Signature for points *)
+module type P = sig
+  type t
+
+  val get_coord : t -> Gg.v2
+
+  val id : t -> int
+
+  val copy : t -> Gg.v2 -> t
+
+end
+
+module Make(Point:P) = struct
+
+  type bezier =
+    { ctrl0:Gg.v2   (* The control point *)
+    ; ctrl1:Gg.v2   (* The control point *)
+    ; p1:Point.t    (* The end point *)
+    }
+
+  module type BUILDER = sig
+    type t
+
+    val repr
+      : t -> (module Repr.M with type point = Point.t and type t = 's) -> 's -> 's
+  end
+
+  type path =
+    | Line of Point.t
+    | Curve of bezier
+
+
+  type step =
+    { point : Point.t
+    ; move : path
+    }
+
+  type t = step array
+
+  module ToFixed = struct
+    type point = Point.t
+
+    type t = int * step list
+
+    let create_path () = 0, []
+
+    (* Start a new path. *)
+    let start point t =
+      let _ = point in
+      t
+
+    let line_to
+      : point -> point -> t -> t
+      = fun p1 p2 (i, t) ->
+        ( i + 1
+        , { point = p1
+          ; move = Line p2
+          }:: t )
+
+    let quadratic_to
+      : (point * Gg.v2 * Gg.v2 * point) -> t -> t
+      = fun (p0,  ctrl0, ctrl1, p1) (i, t) ->
+        let curve = Curve
+            { ctrl0
+            ; ctrl1
+            ; p1} in
+        ( i + 1
+        , { point = p0
+          ; move = curve
+          } ::t)
+
+    let stop t = t
+
+    let get
+      : int * step list -> step array
+      = fun (n, t) ->
+
+        (* The array is initialized with a magic number, and just after
+           filled with the values from the list in reverse. All the elements are set.
+        *)
+        let res = Obj.magic (Array.make n 0) in
+        List.iteri t
+          ~f:(fun i elem -> Array.set res (n - i - 1) elem );
+        res
+  end
+
+  let to_fixed
+    : (module BUILDER with type t = 'a) -> 'a -> t
+    = fun (type s) (module Builder: BUILDER with type t = s) t ->
+      Builder.repr t (module ToFixed) (ToFixed.create_path ())
+      |> ToFixed.get
+
+  let repr
+    : t -> (module Repr.M with type point = Point.t and type t = 's) -> 's -> 's
+    = fun (type s) t (module Repr : Repr.M with type point = Point.t and type t = s) repr ->
+      let repr_bezier p p0 bezier =
+        Repr.quadratic_to
+          ( p0
+          , bezier.ctrl0
+          , bezier.ctrl1
+          , bezier.p1 )
+          p in
+
+      let _, repr = Array.fold_left t
+          ~init:(true, repr)
+          ~f:(fun (first, path) element ->
+              let path =  if first then
+                  Repr.start element.point path
+                else path in
+              match element.move with
+              | Line p1 ->
+                ( false
+                , Repr.line_to element.point p1 path )
+              | Curve bezier ->
+                ( false
+                , repr_bezier path element.point bezier )
+            ) in
+      Repr.stop repr
+
+
+  type approx =
+    { distance : float
+    ; closest_point : Gg.v2
+    ; ratio : float
+    ; p0 : Point.t
+    ; p1 : Point.t }
+
+  (** Return the distance between a given point and the curve. May return
+      None if the point is out of the curve *)
+  let distance
+    : Gg.v2 -> t -> approx option
+    = fun point t ->
+
+      Array.fold_left t
+        ~init:None
+        ~f:(fun res step ->
+            match step.move with
+            | Line p1 ->
+              let box = Gg.Box2.of_pts (Point.get_coord step.point) (Point.get_coord p1) in
+              begin match Gg.Box2.mem point box with
+                | false -> res
+                | true ->
+                  (* TODO Evaluate the normal *)
+                  res
+              end
+            | Curve bezier ->
+
+              let bezier' = Shapes.Bezier.(
+
+                  { p0 = Point.get_coord step.point
+                  ; p1 = Point.get_coord bezier.p1
+                  ; ctrl0 = bezier.ctrl0
+                  ; ctrl1 = bezier.ctrl1 }
+                ) in
+              let ratio, point' = Shapes.Bezier.get_closest_point point bezier' in
+              let distance' = Gg.V2.( norm (point - point') ) in
+              match res with
+              | Some {distance; _} when distance < distance' -> res
+              | _ -> Some
+                       { closest_point = point'
+                       ; distance = distance'
+                       ; p0 = step.point
+                       ; p1 = bezier.p1
+                       ; ratio }
+          )
+
+  let map
+    : t -> (Point.t -> Point.t) -> t
+    = fun t f ->
+      Array.map t
+        ~f:(fun step ->
+            match step.move with
+            | Line p2 ->
+              { point = f step.point
+              ; move = Line (f p2)
+              }
+            | Curve bezier ->
+              let point = f step.point in
+              { point
+              ; move = Curve
+                    { p1 = f bezier.p1
+                    ; ctrl0 = Point.get_coord (f (Point.copy step.point bezier.ctrl0))
+                    ; ctrl1 = Point.get_coord (f (Point.copy bezier.p1 bezier.ctrl1))
+                    }
+              }
+          )
+
+  let iter
+    : t -> f:(Point.t -> unit) -> unit
+    = fun t ~f ->
+      Array.iter t
+        ~f:(fun step ->
+            match step.move with
+            | Line p2 -> f step.point; f p2
+            | Curve bezier -> f step.point ; f bezier.p1
+          )
+
+  let get_point'
+    : step -> Point.t
+    = fun { move ; _} ->
+      match move with
+      | Line p1 -> p1
+      | Curve bezier -> bezier.p1
+
+  (** Associate the return from the bezier point to an existing path *)
+  let assoc_point
+    : Shapes.Bezier.t -> step -> step
+    = fun bezier step ->
+      match step.move with
+      | Line p1
+      | Curve {p1; _} ->
+        let p0' = Point.copy step.point bezier.Shapes.Bezier.p0
+        and p1' = Point.copy p1 bezier.Shapes.Bezier.p1 in
+        { point = p0'
+        ; move = Curve
+              { p1 = p1'
+              ; ctrl0 = bezier.Shapes.Bezier.ctrl0
+              ; ctrl1 = bezier.Shapes.Bezier.ctrl1
+              }
+        }
+
+
+  let build_from_three_points p0 p1 p2 =
+    let bezier =
+      Shapes.Bezier.quadratic_to_cubic
+      @@ Shapes.Bezier.three_points_quadratic
+        (Point.get_coord p0)
+        (Point.get_coord p1)
+        (Point.get_coord p2) in
+
+    (* The middle point is not exactly at the middle anymore (it can have been
+       moved), we have the reevaluate it's position *)
+    let ratio, _ = Shapes.Bezier.get_closest_point
+        (Point.get_coord p1)
+        bezier in
+
+    let b0, b1 = Shapes.Bezier.slice ratio bezier in
+    let p0' = Point.copy p0 b0.Shapes.Bezier.p0
+    and p1' = Point.copy p1 b0.Shapes.Bezier.p1
+    and p2' = Point.copy p2 b1.Shapes.Bezier.p1 in
+
+    [| { point = p0'
+       ; move =
+           Curve { ctrl0 = b0.Shapes.Bezier.ctrl0
+                 ; ctrl1 = b0.Shapes.Bezier.ctrl1
+                 ; p1 = p1'
+                 } }
+     ; { point = p1'
+       ; move = Curve { ctrl0 = b1.Shapes.Bezier.ctrl0
+                      ; ctrl1 = b1.Shapes.Bezier.ctrl1
+                      ; p1 = p2' }
+       } |]
+
+  (** Rebuild the whole curve by evaluating all the points *)
+  let rebuild
+    : t -> t option
+    = fun t ->
+
+      match Array.length t with
+      | 0 -> None
+      | 1 ->
+        let step = Array.get t 0 in
+        begin match step.move with
+          | Curve {p1; _}
+          | Line p1 ->
+            Some
+              [|
+                { point = step.point
+                ; move = Line p1 } |]
+        end
+      | 2 ->
+        let p0 = (Array.get t 0).point
+        and p1 = (Array.get t 1).point
+        and p2 = get_point' @@ Array.get t 1 in
+        Some (build_from_three_points p0 p1 p2)
+
+      | _ ->
+
+        (* Convert all the points in list *)
+        let points = List.init
+            ~len:((Array.length t) )
+            ~f:(fun i -> Point.get_coord @@ get_point' (Array.get t i)) in
+        let p0 = Point.get_coord @@ (Array.get t 0).point in
+
+        let points = p0::points in
+
+        (* We process the whole curve in a single block *)
+        begin match Shapes.Bspline.to_bezier points with
+          | Error `InvalidPath -> None
+          | Ok beziers ->
+
+            (* Now for each point, reassociate the same point information,
+               We should have as many points as before *)
+            let rebuilded = Array.map2 beziers t ~f:assoc_point in
+            Some rebuilded
+        end
+
+  let find_pt_index
+    : Point.t -> step array -> int option
+    = fun point path ->
+      (* First search the element to remove. The counter mark the position of
+         the point to remove, not the segment itself. *)
+      let idx = ref None
+      and counter = ref 0  in
+
+      let _ = Array.exists
+          path
+          ~f:(fun element ->
+              let res =
+                if (Point.id element.point) = (Point.id point) then (
+                  idx := Some (!counter) ;
+                  true
+                ) else match element.move with
+                  | Line p1
+                  | Curve {p1;_} when (Point.id p1) = (Point.id point) ->
+                    idx := Some (!counter+1) ;
+                    true
+                  | _ ->
+                    false
+              in
+              incr counter;
+              res) in
+      !idx
+
+  let remove_point
+    : t -> Point.t -> t option
+    = fun t point ->
+
+      match Array.length t with
+      | 0
+      | 1 -> None
+      | 2 ->
+        (* Two segment, we get the points and transform this into a single line *)
+        let p0 = (Array.get t 0).point
+        and p1 = (Array.get t 1).point
+        and p2 = get_point' @@ Array.get t 1 in
+        let elms = List.filter [p0; p1; p2]
+            ~f:(fun pt -> Point.id pt != Point.id point) in
+        begin match elms with
+          | p0::p1::[] ->
+            Some
+              [| { point = p0
+                 ; move = Line p1 }|]
+          | _ -> None
+        end
+      | l ->
+        match find_pt_index point t with
+        | None -> Some t
+        | Some 0 ->
+          (* Remove the first point *)
+          let path = Array.init (l-1)
+              ~f:( fun i -> Array.get t (i+1)) in
+          Some path
+        | Some n when n = (Array.length t) ->
+          (* Remove the last point *)
+          let path = Array.init (l-1)
+              ~f:( fun i -> Array.get t i) in
+          Some path
+        | Some n ->
+          let path' = Array.init (l-1)
+              ~f:(fun i ->
+                  if i < (n-1) then
+                    Array.get t (i)
+                  else if i = (n-1) then
+                    (* We know that the point is not the first nor the last one.
+                       So it is safe to call n-1 or n + 1 point
+
+                       We have to rebuild the point and set that
+                       point_(-1).id = point_(+1).id
+                    *)
+                    let p0 = (Array.get t i).point in
+
+                    match (Array.get t (i+1)).move with
+                    | Line p1 ->
+                      { point = p0
+                      ; move = Line p1 }
+                    | Curve c ->
+                      { point = p0
+                      ; move = Curve c }
+
+                  else
+                    Array.get t (i+1)
+                ) in
+          rebuild path'
+
+  let first_point
+    : step -> Point.t
+    = fun {point; _} -> point
+
+  let replace_point
+    : t -> Point.t -> t option
+    = fun t p ->
+
+      let add_path paths idx f points =
+        if 0 <= idx && idx < Array.length paths then
+          let path = Array.get t idx in
+          Point.get_coord (f path)
+          :: points
+        else points in
+
+      match Array.length t with
+      | 0 -> None
+      | 1 -> (* Only one point, easy ? *)
+        let step = Array.get t 0 in
+        begin match step.move with
+          | Curve {p1; _}
+          | Line p1 ->
+            let p0 = if (Point.id step.point = Point.id p) then p else step.point
+            and p1 = if (Point.id p1 = Point.id p) then p else p1 in
+            Some [|
+              { point = p0
+              ; move = Line p1 }
+            |]
+        end
+
+      | 2 ->
+        let p0 = (Array.get t 0).point
+        and p1 = (Array.get t 1).point
+        and p2 = get_point' @@ Array.get t 1 in
+
+        let p0 = if (Point.id p0 = Point.id p) then p else p0
+        and p1 = if (Point.id p1 = Point.id p) then p else p1
+        and p2 = if (Point.id p2 = Point.id p) then p else p2 in
+        Some (build_from_three_points p0 p1 p2)
+
+      (* More than two segmend, it is ok for a partial reevaluation *)
+      | _ ->
+        match find_pt_index p t with
+        | None -> None
+        | Some n ->
+          let path = Array.copy t in
+
+          let p0, p1 =
+
+            if n < Array.length path then
+              p, get_point' (Array.get path n)
+            else
+              (Array.get path (n -1)).point, p
+          in
+
+          let min_idx = max (n-3) 0 in
+
+          let points =
+            add_path path (n-3) first_point
+            @@ add_path path (n-2) first_point
+            @@ add_path path (n-1) first_point
+            @@ (fun tl -> (Point.get_coord p)::tl)
+            @@ add_path path n get_point'
+            @@ add_path path (n+1) get_point'
+            @@ add_path path (n+2) get_point'
+            @@ [] in
+
+          (* It is impressive how fast it is to evaluate the curve ! Maybe is the
+             worker not required at all…
+          *)
+          let bezier_opt = Shapes.Bspline.to_bezier points in
+          begin match bezier_opt with
+            | Ok paths ->
+              Array.iteri paths
+                ~f:(fun i bezier ->
+                    (* Only take two points before, and two after  *)
+                    let idx = min_idx + i in
+                    if (n-2 < idx) && (idx < n +2) && idx < Array.length path then
+                      Array.set path idx (assoc_point bezier (Array.get path idx))
+                  );
+              Some path
+            | Error _ ->
+              let bezier', _ = Shapes.Bezier.three_points_quadratic
+                  (Point.get_coord p)
+                  (Point.get_coord @@ get_point' (Array.get path 0))
+                  (Point.get_coord @@ get_point' (Array.get path 1))
+                               |> Shapes.Bezier.quadratic_to_cubic
+                               |> Shapes.Bezier.slice 0.5
+              in
+              Array.set path 0
+                { point = p0
+                ; move = (Curve
+                            { ctrl0 = bezier'.Shapes.Bezier.ctrl0
+                            ; ctrl1 = bezier'.Shapes.Bezier.ctrl1
+                            ; p1
+                            })
+                };
+              Some path
+          end
+end
diff --git a/script.it/path/fixed.mli b/script.it/path/fixed.mli
new file mode 100755
index 0000000..111187c
--- /dev/null
+++ b/script.it/path/fixed.mli
@@ -0,0 +1,81 @@
+(** Signature for points *)
+module type P = sig
+  type t
+
+  val get_coord : t -> Gg.v2
+
+  val id : t -> int
+
+  val copy : t -> Gg.v2 -> t
+
+end
+
+module Make(Point:P) : sig
+
+  module type BUILDER = sig
+    type t
+
+    val repr
+      : t -> (module Repr.M with type point = Point.t and type t = 's) -> 's -> 's
+
+  end
+
+  type t 
+
+  (** Create a path from a builder *)
+  val to_fixed
+    : (module BUILDER with type t = 'a) -> 'a -> t
+
+  (** Represent the path *)
+  val repr
+    : t -> (module Repr.M with type point = Point.t and type t = 's) -> 's -> 's
+
+  (** Structure to represent all the required information for evaluating the
+      distance between a point and a path *)
+  type approx =
+    { distance : float
+    ; closest_point : Gg.v2
+    ; ratio : float
+    ; p0 : Point.t
+    ; p1 : Point.t }
+
+  (** Return the distance between a given point and the curve. May return
+      None if the point is out of the curve *)
+  val distance
+    : Gg.v2 -> t -> approx option
+
+  (** Iterate over a path *)
+  val iter
+    : t -> f:(Point.t -> unit) -> unit
+
+  (** Map all the points in the path *)
+  val map
+    : t -> (Point.t -> Point.t) -> t
+
+  (** Reevaluate all the control points on the path in order to get a smooth
+      curve *)
+  val rebuild 
+    : t -> t option
+
+  (** Delete a point in the path. 
+
+      Reconnect the path without the point removed, and reevaluate all the
+      control points from the nodes
+
+      return None if the point is not present in the curve
+  *)
+  val remove_point 
+    : t -> Point.t -> t option
+
+  (** Replace a point by the given one. 
+
+      An existing point with the same id shall be present in the path. 
+
+      The path is not fully evaluated, and rebuild shall be runned in order to
+      get the path completely smooth.
+
+  *)
+  val replace_point 
+    : t -> Point.t -> t option
+
+end
diff --git a/script.it/path/path.ml b/script.it/path/path.ml
new file mode 100755
index 0000000..ea90de4
--- /dev/null
+++ b/script.it/path/path.ml
@@ -0,0 +1,7 @@
+(** Common  module for ensuring that the function is evaluated only once *)
+
+module Point = Point
+module Repr = Repr
+module Path_Builder = Builder.Make(Point)
+module Fixed = Fixed.Make(Point)
+
diff --git a/script.it/path/point.ml b/script.it/path/point.ml
new file mode 100755
index 0000000..ec6f8ad
--- /dev/null
+++ b/script.it/path/point.ml
@@ -0,0 +1,77 @@
+let internal_id = ref 0
+
+type t =
+  { p: Gg.v2
+  ; size : float
+  ; angle: float
+  ; stamp : float
+  ; id : int
+  }
+
+let empty =
+  { p = Gg.V2.of_tuple (0., 0.)
+  ; size = 0.
+  ; angle = 0.
+  ; stamp = 0.
+  ; id = 0
+  }
+
+let create ~angle ~width ~stamp ~x ~y =
+
+  incr internal_id;
+  { p = Gg.V2.v x y
+  ; size = width
+  ; angle = Gg.Float.rad_of_deg (180. -. angle )
+  ; stamp
+  ; id = !internal_id
+  }
+
+let copy point p =
+  { point with
+    p
+  }
+
+let set_angle p angle =
+  { p with angle = Gg.Float.rad_of_deg (180. -. angle) }
+
+let get_angle { angle; _} = 180. -. (Gg.Float.deg_of_rad angle)
+
+let set_width p size =
+  { p with size }
+
+let get_width { size; _} = size
+
+let (+) p1 p2 =
+  { p1 with p = Gg.V2.(+) p1.p p2 }
+
+let get_coord { p; _ } = p
+
+let get_stamp { stamp; _} = stamp
+
+let get_coord'
+  : t -> Gg.v2
+  = fun t ->
+    let open Gg.V2 in
+    let trans = of_polar @@ v t.size t.angle in
+    t.p + trans
+
+let mix
+  : float -> Gg.v2 -> t -> t -> t
+  = fun f point p0 p1 ->
+    let angle0 = p0.angle
+    and angle1 = p1.angle
+    and width0 = get_width p0
+    and width1 = get_width p1
+    and stamp0 = get_stamp p0
+    and stamp1 = get_stamp p1 in
+    let angle = angle0 +. f *. ( angle1 -. angle0 ) in
+    let width = width0 +. f *. ( width1 -. width0 ) in
+    let stamp = stamp0 +. f *. ( stamp1 -. stamp0 ) in
+    { p = point
+    ; size = width
+    ; angle
+    ; stamp
+    ; id = !internal_id
+    }
+
+let id { id; _} = id
diff --git a/script.it/path/point.mli b/script.it/path/point.mli
new file mode 100755
index 0000000..fe4cb45
--- /dev/null
+++ b/script.it/path/point.mli
@@ -0,0 +1,40 @@
+type t
+
+(** Return the point id. Each id is unique *)
+val id 
+  : t -> int
+
+val empty : t
+
+val (+): t -> Gg.v2 -> t
+
+val get_coord : t -> Gg.v2
+
+val get_stamp : t -> float
+
+val create: angle:float -> width:float -> stamp:float -> x:float -> y:float -> t
+
+(** Return a copy of the point at given posistion 
+
+    This is a true copy, the id will be the same for the two points
+    TODO : Should this be renamed set_position ?
+
+*)
+val copy : t -> Gg.v2 -> t
+
+val set_angle : t -> float -> t
+
+val get_angle : t -> float
+
+val set_width: t -> float -> t
+
+val get_width: t -> float
+
+val get_coord'
+  : t -> Gg.v2
+
+(** [mix f point p0 p1] create a new point at the position point, with the
+    characteristics from p0 and p1 *)
+val mix
+  : float -> Gg.v2 -> t -> t -> t
+
diff --git a/script.it/path/repr.ml b/script.it/path/repr.ml
new file mode 100755
index 0000000..17fd914
--- /dev/null
+++ b/script.it/path/repr.ml
@@ -0,0 +1,19 @@
+module type M = sig
+
+  type point
+
+  type t
+
+  (* Start a new path. *)
+  val start
+    : point -> t -> t
+
+  val line_to
+    : point -> point -> t -> t
+
+  val quadratic_to
+    : (point * Gg.v2 * Gg.v2 * point) -> t -> t
+
+  val stop
+    : t -> t
+end
diff --git a/script.it/shapes/bezier.ml b/script.it/shapes/bezier.ml
new file mode 100755
index 0000000..f5f288c
--- /dev/null
+++ b/script.it/shapes/bezier.ml
@@ -0,0 +1,228 @@
+(**
+
+   Bezier curve
+*)
+
+module Utils = Tools.Utils
+
+type quadratic =
+  { p0:Gg.v2     (* The starting point *)
+  ; p1:Gg.v2     (* The end point *)
+  ; ctrl:Gg.v2 } (* The control point *)
+
+
+type t =
+  { p0:Gg.v2      (* The starting point *)
+  ; p1:Gg.v2      (* The end point *)
+  ; ctrl0:Gg.v2   (* The control point *)
+  ; ctrl1:Gg.v2 } (* The control point *)
+
+
+(**
+   Build a control point for a quadratic curve for passing throuht 3 points.
+   taken from https://xuhehuan.com/2608.html
+
+
+   also look to https://pomax.github.io/bezierinfo/#pointcurves
+*)
+let three_points_quadratic
+  : Gg.v2 -> Gg.v2 -> Gg.v2 -> quadratic
+  = fun p0 c1 p1 ->
+
+    let open Gg.V2 in
+
+    let vect_1 = p0 - c1
+    and vect_2 = p1 - c1 in
+    let norm1 = norm vect_1
+    and norm2 = norm vect_2 in
+    let v = (Float.sqrt (norm1 *. norm2)) /. 2. in
+
+    let ctrl = c1 - v * (( vect_1 / norm1) + (vect_2 / norm2)) in
+    {p0; p1; ctrl}
+
+(**
+
+   Convert a cubic bezier curve into a quadratic one
+
+*)
+let quadratic_to_cubic
+  : quadratic -> t
+  = fun {p0; p1; ctrl} ->
+
+    let coef = 2. /. 3. in
+
+    let open Gg.V2 in
+    { p0
+    ; p1
+    ; ctrl0 = mix p0 ctrl coef
+    ; ctrl1 = mix p1 ctrl coef }
+
+
+
+let abc_ratio
+  : int -> float -> float
+  = fun n t ->
+    let n' = Float.of_int n in
+    let bottom = (Float.pow t n') +. (Float.pow  (1. -. t) n') in
+    let top = bottom -. 1. in
+    Float.abs (top /. bottom)
+
+let half_cubic_ratio = abc_ratio 3 0.5
+
+exception Not_found
+
+(**
+
+   https://pomax.github.io/bezierinfo/#pointcurves
+
+*)
+let three_points_cubic
+  : float -> Gg.v2 -> Gg.v2 -> Gg.v2 -> t
+  = fun f p0 p1 p2 ->
+
+    let open Gg.V2 in
+
+    let c = half ( p0 + p2) in
+    let a = p1 + ((p1 - c) / half_cubic_ratio) in
+
+    let vect1_0 = p1 - p0 in
+    let vect2_0 = p2 - p0 in
+
+    let d1 = norm vect1_0
+    and d2 = norm (p2 - p1) in
+    let t = d1 /. (d1 +. d2) in
+
+    let angle_1_0 = angle vect1_0
+    and angle_2_0 = angle vect2_0 in
+
+    (* get our e1-e2 distances *)
+    let angle = mod_float
+        (Gg.Float.two_pi
+         +. angle_2_0
+         -. angle_1_0)
+        Gg.Float.two_pi in
+
+    let distance = (norm vect2_0) *. f in
+
+    let bc =
+      if angle < 0. || angle > Gg.Float.pi then
+        Float.(neg distance)
+      else
+        distance in
+    let de1  = t *. bc
+    and de2 = (1. -. t) *. bc in
+
+    (* get the circle-aligned slope as normalized dx/dy *)
+    let center = Utils.center p0 p1 p2 in
+    match center with
+    | None -> raise Not_found
+    | Some center ->
+      let t' = p1 - center in
+      let tangent0 = v
+          ((x p1) -. (y t'))
+          ((y p1) +. (x t'))
+      and tangent1 = v
+          ((x p1) +. (y t'))
+          ((y p1) -. (x t')) in
+
+      let d = unit (tangent1 - tangent0) in
+
+      (* then set up an e1 and e2 parallel to the baseline *)
+      let e1 = p1 + de1 * d
+      and e2 = p1 - de2 * d in
+
+      (* then use those e1/e2 to derive the new hull coordinates *)
+      let v1 = a + (e1 - a) / (1. -. t)
+      and v2 = a + (e2 - a) / t in
+
+      let ctrl0 = p0 + (v1 - p0) / t
+      and ctrl1 = p2 + (v2 -p2) / (1. -. t) in
+
+      {p0; p1 = p2; ctrl0; ctrl1}
+
+(** Split a bezier curve in two at a given position *)
+let slice
+  : float -> t -> t * t
+  = fun t {p0; p1; ctrl0; ctrl1} ->
+
+    let mix p1 p2 = Gg.V2.mix p1 p2 t in
+
+    let p12 = mix p0 ctrl0
+    and p23 = mix ctrl0 ctrl1
+    and p34 = mix ctrl1 p1 in
+
+    let p123 = mix p12 p23
+    and p234 = mix p23 p34 in
+
+    let p1234 = mix p123 p234 in
+
+    ( { p0
+      ; ctrl0 = p12
+      ; ctrl1 = p123
+      ; p1 = p1234 }
+    , { p0 = p1234
+      ; ctrl0 = p234
+      ; ctrl1 = p34
+      ; p1 } )
+
+
+let get_closest_point
+  : Gg.v2 -> t -> float * Gg.v2
+  = fun point t ->
+
+    let rec f min max t =
+
+      (* First devide the curve in two *)
+      let seq_0, seq_1 = slice 0.5 t in
+      let avg = (min +. max) /. 2. in
+
+      let p0 = t.p0
+      and p1 = t.p1
+      and p01 = seq_0.p1 in (* seq_0.p1 = seq_1.p0 *)
+
+      let open Gg.V2 in
+      let center0 = mix p0 p01 0.5
+      and center1 = mix p01 p1 0.5 in
+
+      if Tools.Utils.equal_point 0.001 p0 p1 then
+        avg, p01
+      else if (norm (point - center0)) < (norm (point - center1)) then
+        f min avg seq_0
+      else
+        f avg max seq_1
+
+    in f 0. 1. t
+
+let reverse
+  : t -> t
+  = fun bezier ->
+    {
+      p0 = bezier.p1
+    ; p1 = bezier.p0
+    ; ctrl0 = bezier.ctrl1
+    ; ctrl1 = bezier.ctrl0 }
+
+(**
+
+   see https://github.com/Pomax/BezierInfo-2/blob/master/docs/js/graphics-element/lib/bezierjs/bezier.js#L504
+
+   let root
+   : t -> 'a
+   = fun bezier ->
+
+    let accept
+      : float -> bool
+      = fun t ->
+        0. <= t && t <= 1. in
+
+    let cuberoot v =
+      if v < 0.  then
+        Float.pow (Float.neg v) ( 1. /. 3.)
+        |> Float.neg
+      else Float.pow v (1. /. 3.) in
+
+
+
+
+    failwith "Non implemented"
+*)
diff --git a/script.it/shapes/bezier.mli b/script.it/shapes/bezier.mli
new file mode 100755
index 0000000..2f5bbcf
--- /dev/null
+++ b/script.it/shapes/bezier.mli
@@ -0,0 +1,40 @@
+type t =
+  { p0:Gg.v2      (* The starting point *)
+  ; p1:Gg.v2      (* The end point *)
+  ; ctrl0:Gg.v2   (* The control point *)
+  ; ctrl1:Gg.v2 } (* The control point *)
+
+type quadratic
+
+(**
+   Build a control point for a quadratic curve for passing throuht 3 points.
+   taken from https://xuhehuan.com/2608.html
+
+
+   also look to https://pomax.github.io/bezierinfo/#pointcurves
+*)
+val three_points_quadratic
+  : Gg.v2 -> Gg.v2 -> Gg.v2 -> quadratic
+
+(**
+   Create a curve from three points.
+
+   This is an implementation for
+   https://pomax.github.io/bezierinfo/#pointcurves
+
+*)
+val three_points_cubic
+  : float -> Gg.v2 -> Gg.v2 -> Gg.v2 -> t
+
+val quadratic_to_cubic
+  : quadratic -> t
+
+(** Split a bezier curve in two at a given position *)
+val slice
+  : float -> t -> t * t
+
+(** Return the closest point to the curve by approximation *)
+val get_closest_point
+  : Gg.v2 -> t -> float * Gg.v2
+
+val reverse: t -> t
diff --git a/script.it/shapes/bspline.ml b/script.it/shapes/bspline.ml
new file mode 100755
index 0000000..bb60227
--- /dev/null
+++ b/script.it/shapes/bspline.ml
@@ -0,0 +1,149 @@
+open StdLabels
+
+type err = [`InvalidPath ]
+
+module M = Matrix.MakeMatrix (struct
+
+    type t = Float.t
+
+    let compare a b =
+
+      let v = Float.compare a b in
+      if v = 0 then Matrix.Order.Equal
+      else if v > 0 then Matrix.Order.Greater
+      else Matrix.Order.Less
+
+    let zero = Float.zero
+    let one = Float.one
+
+    let divide = (/.)
+    let multiply = ( *. )
+    let add = (+.)
+    let subtract = (-.)
+    exception NonElt
+
+
+  end)
+
+type t = Gg.v2 list
+
+let from_points
+  : Gg.v2 array -> (Gg.v2 array, [> `InvalidPath]) Result.t
+  = fun points ->
+
+    let n = (Array.length points - 2) in
+
+    if n <= 1 then
+      Result.error `InvalidPath
+    else
+
+      (* Create the initial matrix.
+
+         The matrix is augmented with two additionals columns, which will be
+         populated with the points from the path.
+      *)
+      let arr = Array.init n ~f:(fun line ->
+          Array.init (n +2) ~f:(fun row ->
+              match row - line with
+              | (-1) -> 1.
+              | 0 -> 4.
+              | 1 -> 1.
+              | _ -> 0.
+            )
+        ) in
+      let matrix = M.from_array arr in
+
+      (* Add the points from the augmented matrix *)
+      let points_array = points in
+      for line = 0 to (n -1) do
+
+        let point =
+          if line = 0 then
+            let p0 = points_array.(0)
+            and p1 = points_array.(1) in
+            Gg.V2.(6. * p1 - p0)
+          else if (line + 1) = n  then
+            let p_n_2 = points_array.(n)
+            and p_n_1 = points_array.(n + 1) in
+            Gg.V2.(6. * p_n_2 - p_n_1)
+          else
+            let n' = line + 1 in
+            Gg.V2.(6. * points_array.(n'))
+        in
+        let x = (Gg.V2.x point)
+        and y = (Gg.V2.y point) in
+
+
+        M.set_elt matrix (line + 1, n + 1) x;
+        M.set_elt matrix (line + 1, n + 2) y;
+      done;
+
+      (* Resolve the matrix *)
+      let res' = M.row_reduce matrix in
+
+      (* Extract the result as points *)
+      let _, col_x = M.get_column res' (n + 1)
+      and _, col_y = M.get_column res' (n + 2) in
+
+      (* Build the result *)
+      let res = Array.make (n + 2) (Array.get points_array (n + 1) ) in
+      for i = 1 to n do
+        let point = Gg.V2.v col_x.(i - 1) col_y.(i - 1) in
+        Array.set res i point;
+      done;
+      Array.set res 0 (Array.get points_array 0);
+      Result.ok res
+
+let (let*) = Result.bind
+
+(** Build a continue curve from path
+
+    see https://www.math.ucla.edu/~baker/149.1.02w/handouts/dd_splines.pdf
+*)
+let to_bezier
+  : ?connexion0:Gg.v2 -> ?connexion1:Gg.v2 -> t -> (Bezier.t array, [> `InvalidPath]) Result.t
+  = fun ?connexion0 ?connexion1 points ->
+
+    let points' = match connexion0 with
+      | None -> points
+      | Some pt -> pt::points in
+
+    let arr_points = match connexion1 with
+      | None -> Array.of_list points'
+      | Some pt ->
+        let arr = Array.make (1 + (List.length points')) pt in
+        List.iteri points'
+          ~f:(fun i value -> Array.set arr i value);
+        arr in
+
+    let* bspline_points = from_points arr_points in
+
+    let start = match connexion0 with
+      | None -> 1
+      | Some _ -> 2
+    and end_ = match connexion1 with
+      | None -> (Array.length bspline_points) - 1
+      | Some _ -> (Array.length bspline_points) - 2 in
+
+    let result = Array.init (end_ - start + 1) ~f:(fun i ->
+
+        let i = i + start in
+
+        let prev_b = Array.get bspline_points (i - 1)
+        and bpoint = Array.get bspline_points i
+        and prev_p = Array.get arr_points (i - 1)
+        and point = Array.get arr_points i in
+        let ctrl0 = Gg.V2.(mix prev_b bpoint (1. /. 3.))
+        and ctrl1 = Gg.V2.(mix prev_b bpoint (2. /. 3.)) in
+
+        let bezier =
+          { Bezier.p0 = prev_p
+          ; Bezier.p1 = point
+          ; Bezier.ctrl0
+          ; Bezier.ctrl1 } in
+
+        bezier
+
+      ) in
+    Result.Ok result
+
diff --git a/script.it/shapes/bspline.mli b/script.it/shapes/bspline.mli
new file mode 100755
index 0000000..a36aa22
--- /dev/null
+++ b/script.it/shapes/bspline.mli
@@ -0,0 +1,24 @@
+type t
+
+type err = 
+  [ `InvalidPath (* Too few points in the path for building the curve *)
+  ]
+
+(** Convert a list of points into a beziers curves. 
+
+    At least 4 points are required for building the path.
+
+    [to_bezier ~connexion points] create a list of beziers segments joining all
+    the points together. 
+
+    [connexion0] add a virtual point in the begining for helping to get the
+    appropriate tangent when connecting path together
+
+    [connexion1] does the same at the end
+
+*)
+val to_bezier
+  :  ?connexion0:Gg.v2
+  -> ?connexion1:Gg.v2
+  -> Gg.v2 list 
+  -> (Bezier.t array, [> err]) Result.t
diff --git a/script.it/shapes/dd_splines.pdf b/script.it/shapes/dd_splines.pdf
new file mode 100755
index 0000000..2618162
--- /dev/null
+++ b/script.it/shapes/dd_splines.pdf
diff --git a/script.it/shapes/dune b/script.it/shapes/dune
new file mode 100755
index 0000000..d03a217
--- /dev/null
+++ b/script.it/shapes/dune
@@ -0,0 +1,7 @@
+(library
+ (name shapes)
+ (libraries 
+   tools
+   matrix
+   )
+ )
diff --git a/script.it/shapes/matrix/EltsI.ml b/script.it/shapes/matrix/EltsI.ml
new file mode 100755
index 0000000..fcfdb50
--- /dev/null
+++ b/script.it/shapes/matrix/EltsI.ml
@@ -0,0 +1,28 @@
+module type ORDERED_AND_OPERATIONAL =
+sig
+
+  (* Exception for from_string. Is raised when from_string is passed something
+   * that is not an elt *)
+  exception NonElt
+
+  type t
+
+  (* The zero element *)
+  val zero : t
+
+  (* The one element *)
+  val one: t
+
+  (* ts must be comparable *)
+  val compare : t -> t -> Order.order
+
+  (* Basic mathematical operations must be possible *)
+  val add: t -> t -> t
+
+  val subtract: t -> t -> t
+
+  val multiply: t -> t -> t
+
+  val divide: t -> t -> t
+
+end
diff --git a/script.it/shapes/matrix/Helpers.ml b/script.it/shapes/matrix/Helpers.ml
new file mode 100755
index 0000000..6980052
--- /dev/null
+++ b/script.it/shapes/matrix/Helpers.ml
@@ -0,0 +1,16 @@
+(* Takes in a string and a separator, and separates the string into a list of
+ * substrings where each substring is between two separators or between a
+ * separator and the beginning/end of the string *)
+let explode (s: string) (space: string) : string list =
+  let rec build (curr: string) (buffer: string) (lst: string list) : string list =
+    let len = String.length curr in
+    if len = 0 then buffer::lst
+    else 
+      let c = String.sub curr (len - 1) 1 in
+      if len = 1 then (c ^ buffer)::lst
+      else 
+        let s' = String.sub curr 0 (len - 1) in
+        if c = space then build s' "" (buffer::lst)
+        else build s' (c ^ buffer) lst in
+  build (String.trim s) "" []
+
diff --git a/script.it/shapes/matrix/Matrix.ml b/script.it/shapes/matrix/Matrix.ml
new file mode 100755
index 0000000..7f1d54b
--- /dev/null
+++ b/script.it/shapes/matrix/Matrix.ml
@@ -0,0 +1,529 @@
+open Order
+
+module Order = Order
+
+(*************** Exceptions ***************)
+
+exception NonSquare
+exception ImproperDimensions
+
+(* Functor so we can Abstract away! *)
+module MakeMatrix (C: EltsI.ORDERED_AND_OPERATIONAL) :
+  (MatrixI.MATRIX with type elt = C.t) =
+struct
+
+
+  (*************** End Exceptions ***************)
+
+  (*************** Types ***************)
+
+  type elt = C.t
+
+  (* A matrix is a pair of dimension (n x p) and a array of arrays
+   * the first array is the row (n) and the second the column (p) *)
+  type matrix = (int * int) * (elt array array)
+
+  (*************** End Types ***************)
+
+  (*************** Base Functions ***************)
+
+  (* catching negative dimensions AND 0 dimensions and too large
+   * of a dimension so we don't have to worry about it later *)
+  let empty (rows: int) (columns: int) : matrix =
+    if rows > 0 && columns > 0 then
+      try
+        let m = Array.make_matrix rows columns C.zero in ((rows,columns),m)
+      with _ ->
+        raise ImproperDimensions
+    else (* dimension is negative or 0 *)
+      raise ImproperDimensions
+
+  (*************** End Base Functions ***************)
+
+  (*************** Helper Functions ***************)
+
+  (* get's the nth row of a matrix and returns (r, row) where r is the length
+   * of the row and row is a COPY of the original row. For example, calling
+   * calling get_row m 1 will return (3, |1 3 4 |)
+   *         ________
+   *    m = | 1 3 4 |
+   *        |*2 5 6 |
+  *)
+  (* aside: we don't check whether n < 1 because of our matrix invariant *)
+  let get_row (((n,p),m): matrix) (row: int) : int * elt array =
+    if row <= n then
+      let row' = Array.map (fun x -> x) m.(row - 1) in
+      (p, row')
+    else
+      raise (Failure "Row out of bounds.")
+
+  (* similar to get_row. For m, get_column m 1 will return (2, |1 2|) *)
+  let get_column (((n,p),m): matrix) (column: int) : int * elt array =
+    if column <= p then
+      begin
+        let column' = Array.make n C.zero in
+        for i = 0 to n - 1 do
+          column'.(i) <- m.(i).(column - 1)
+        done;
+        (n, column')
+      end
+    else
+      raise (Failure "Column out of bounds.")
+
+  (* sets the nth row of the matrix m to the specified array a.
+   * This is done IN-PLACE. Therefore the function returns unit.  You should
+   * nonetheless enfore immutability whenever possible.  For a clarification on
+   * what nth row means, look at comment for get_row above. *)
+  let set_row (((n, p), m): matrix) (row: int) (a: elt array) : unit =
+    if row <= n then
+      begin
+        assert(Array.length a = p);
+        for i = 0 to p - 1 do
+          m.(row - 1).(i) <- a.(i)
+        done;
+      end
+    else
+      raise (Failure "Row out of bounds.")
+
+  (* Similar to set_row but sets the nth column instead *)
+  let set_column (((n,p),m): matrix) (column: int) (a: elt array) : unit =
+    if column <= p then
+      begin
+        assert(Array.length a = n);
+        for i = 0 to n - 1 do
+          m.(i).(column - 1) <- a.(i)
+        done;
+      end
+    else
+      raise (Failure "Column out of bounds.")
+
+  (* returns the ij-th element of a matrix (not-zero indexed) *)
+  let get_elt (((n,p),m): matrix) ((i,j): int*int) : elt =
+    if i <= n && j <= p then
+      m.(i - 1).(j - 1)
+    else
+      raise ImproperDimensions
+
+  (* Changes the i,j-th element of a matrix to e. Is not zero-indexed, and
+   * changes the matrix in place *)
+  let set_elt (((n,p),m): matrix) ((i,j): int*int) (e: elt) : unit =
+    if i <= n && j <= p then
+      m.(i - 1).(j - 1) <- e
+    else
+      raise ImproperDimensions
+
+  (* similar to map, but applies to function to the entire matrix
+   * Returns a new matrix *)
+  let map (f: elt -> elt) (mat: matrix) : matrix =
+    let (dim,m) = mat in
+    (dim, Array.map (Array.map f) m)
+
+  (* Just some wrapping of Array.iter made for Matrices! *)
+  let iter (f: elt -> unit) (mat: matrix) : unit =
+    let _, m = mat in
+    Array.iter (Array.iter f) m
+
+  (* Just some wrapping of Array.iteri. Useful for pretty
+   * printing matrix. The index is (i,j). NOT zero-indexed *)
+  let iteri (f: int -> int -> elt -> unit) (mat: matrix) : unit =
+    let _, m = mat in
+    Array.iteri (fun i row -> Array.iteri (fun j e -> f i j e) row) m
+
+  (* folds over each row using base case u and function f *)
+  (* could be a bit more efficient? *)
+  let reduce (f: 'a -> elt -> 'a) (u: 'a) (((p,q),m): matrix) : 'a =
+    let total = ref u in
+    for i = 0 to p - 1 do
+      for j = 0 to q - 1 do
+        total := f (!total) m.(i).(j)
+      done;
+    done;
+    !total
+
+  let fold_row ~(f: elt array -> 'b) ((_,m): matrix) : 'b list =
+
+    let call_row acc v = (f v)::acc in
+    Array.fold_left call_row [] m
+    |> List.rev
+
+
+
+
+  (* given two arrays, this will calculate their dot product *)
+  (* It seems a little funky, but this is done for efficiency's sake.
+   * In short, it tries to multiply each element by it's respective
+   * element until the one array is indexed out of bounds. If the
+   * other array is also out of bounds, then it returns their value.
+   * Otherwise, the arrays were the wrong size and raises ImproperDimension
+
+     THE ABOVE COMMENT HAS NOT BEEN IMPLEMENTED
+
+     Instead we calculate the length before starting
+  *)
+  let dot (v1: elt array) (v2: elt array) : elt =
+    let rec dotting (i: int) (total: elt) : elt =
+      if i = 0 then total
+      else
+        let curr = C.multiply v1.(i-1) v2.(i-1) in
+        dotting (i - 1) (C.add curr total) in
+    let len1, len2 = Array.length v1, Array.length v2 in
+    if len1 = len2 then dotting len1 C.zero
+    else raise ImproperDimensions
+
+  (* function to expose the dimensions of a matrix *)
+  let get_dimensions (m: matrix) : (int * int) =
+    let ((x,y), _) = m in (x,y)
+
+  (*************** End Helper Functions ***************)
+
+
+  (*************** Primary Matrix Functions ***************)
+
+  (* scales a matrix by the appropriate factor *)
+  let scale (m: matrix) (sc: elt) : matrix = map (C.multiply sc) m
+
+  (* Generates a matrix from a list of lists. The inners lists are the rows *)
+  let from_list (lsts : elt list list) : matrix =
+    let check_length (length: int) (lst: elt list) : int =
+      if List.length lst = length then length
+      else raise ImproperDimensions in
+    let p = List.length lsts in
+    match lsts with
+    | [] -> raise ImproperDimensions
+    | hd::tl ->
+      let len = List.length hd in
+      if List.fold_left check_length len tl = len then
+        ((p,len),Array.map Array.of_list (Array.of_list lsts))
+      else
+        raise ImproperDimensions
+
+  (* Generates a matrix from a list of lists. The inners lists are the rows *)
+  let from_array (arrs : elt array array) : matrix =
+    let check_length (length: int) (arr: elt array) : unit =
+      if Array.length arr = length then ()
+      else raise ImproperDimensions in
+    let p = Array.length arrs in
+    match Array.length arrs with
+    | 0 -> raise ImproperDimensions
+    | _ ->
+      let len = Array.length (Array.get arrs 0) in
+      Array.iter (check_length len) arrs;
+      ((p, len), arrs)
+
+  (* Adds two matrices. They must have the same dimensions *)
+  let add ((dim1,m1): matrix) ((dim2,m2): matrix) : matrix =
+    if dim1 = dim2 then
+      let n, p = dim1 in
+      let (dim', sum_m) = empty n p in
+      for i = 0 to n - 1 do
+        for j = 0 to p - 1 do
+          sum_m.(i).(j) <- C.add m1.(i).(j) m2.(i).(j)
+        done;
+      done;
+      (dim',sum_m)
+    else
+      raise ImproperDimensions
+
+
+  (* Multiplies two matrices. If the matrices have dimensions m x n and p x q, n
+   * and p must be equal, and the resulting matrix will have dimension n x q *)
+  let mult (matrix1: matrix) (matrix2: matrix) : matrix =
+    let ((m,n), _), ((p,q), _) = matrix1, matrix2 in
+    if n = p then
+      let (dim, result) = empty m q in
+      for i = 0 to m - 1 do
+        for j = 0 to q - 1 do
+          let (_,row), (_,column) = get_row matrix1 (i + 1),
+                                    get_column matrix2 (j + 1) in
+          result.(i).(j) <- dot row column
+        done;
+      done;
+      (dim,result)
+    else
+      raise ImproperDimensions
+
+  (*************** Helper Functions for Row Reduce ***************)
+
+  (*
+  (* returns the index of the first non-zero elt in an array*)
+  let zero (arr: elt array) : int option =
+    let index = ref 1 in
+    let empty (i: int option) (e: elt) : int option =
+      match i, C.compare e C.zero with
+      | None, Equal -> (index := !index + 1; None)
+      | None, _ -> Some (!index)
+      | _, _ -> i in
+    Array.fold_left empty None arr
+
+  (* returns the the location of the nth non-zero
+   * element in the matrix. Scans column wise.  So the nth non-zero element is
+   * the FIRST non-zero element in the nth non-zero column *)
+  let nth_nz_location (m: matrix) (_: int): (int*int) option =
+    let ((n,p), _) = m in
+    let rec check_col (to_skip: int) (j: int) =
+      if j <= p then
+        let (_,col) = get_column m j in
+        match zero col with
+        | None -> check_col to_skip (j + 1)
+        | Some i ->
+          if to_skip = 0 then
+            Some (i,j)
+          else (* we want a later column *)
+            check_col (to_skip - 1) (j + 1)
+      else None in
+    check_col (n - 1) 1
+
+  (* returns the the location of the first
+   * non-zero and non-one elt. Scans column wise, from
+   * left to right. Basically, it ignores columns
+   * that are all zero or that *)
+  let fst_nz_no_loc (m: matrix): (int*int) option =
+    let ((_, p), _) = m in
+    let rec check_col (j: int) =
+      if j <= p then
+        let (_,col) = get_column m j in
+        match zero col with
+        | None -> check_col (j + 1)
+        | Some i ->
+          match C.compare col.(i-1) C.one with
+          | Equal -> check_col (j + 1)
+          | _ -> Some (i,j)
+      else None in
+    check_col 1
+    *)
+
+  (* Compares two elements in an elt array and returns the greater and its
+   * index. Is a helper function for find_max_col_index *)
+  let compare_helper (e1: elt) (e2: elt) (ind1: int) (ind2: int) : (elt*int) =
+    match C.compare e1 e2 with
+    | Equal -> (e2, ind2)
+    | Greater -> (e1, ind1)
+    | Less -> (e2, ind2)
+
+  (* Finds the element with the greatest absolute value in a column. Is not
+   * 0-indexed. If two elements are both the maximum value, returns the one with
+   * the lowest index. Returns None if this element is zero (if column is all 0)
+  *)
+  let find_max_col_index (array1: elt array) (start_index: int) : int option =
+    let rec find_index (max_index: int) (curr_max: elt) (curr_index: int)
+        (arr: elt array) =
+      if curr_index = Array.length arr then
+        (if curr_max = C.zero then None
+         else Some (max_index+1)) (* Arrays are 0-indexed but matrices aren't *)
+      else
+        (match C.compare arr.(curr_index) C.zero with
+         | Equal -> find_index max_index curr_max (curr_index+1) arr
+         | Greater ->
+           (let (el, index) = compare_helper (arr.(curr_index))
+                curr_max curr_index max_index in
+            find_index index el (curr_index+1) arr)
+         | Less ->
+           (let abs_curr_elt = C.subtract C.zero arr.(curr_index) in
+            let (el, index) = compare_helper abs_curr_elt curr_max curr_index
+                max_index in
+            find_index index el (curr_index+1) arr))
+    in
+    find_index 0 C.zero (start_index -1) array1
+
+  (* Basic row operations *)
+  (* Scales a row by sc *)
+  let scale_row (m: matrix) (num: int) (sc: elt) : unit =
+    let (_, row) = get_row m num in
+    let new_row = Array.map (C.multiply sc) row in
+    set_row m num new_row
+
+  (* Swaps two rows of a matrix *)
+  let swap_row (m: matrix) (r1: int) (r2: int) : unit =
+    let (len1, row1) = get_row m r1 in
+    let (len2, row2) = get_row m r2 in
+    let _ = assert (len1 = len2) in
+    let _ = set_row m r1 row2 in
+    let _ = set_row m r2 row1 in
+    ()
+
+  (* Subtracts a multiple of r2 from r1 *)
+  let sub_mult (m: matrix) (r1: int) (r2: int) (sc: elt) : unit =
+    let (len1, row1) = get_row m r1 in
+    let (len2, row2) = get_row m r2 in
+    let _ = assert (len1 = len2) in
+    for i = 0 to len1 - 1 do (* Arrays are 0-indexed *)
+      row1.(i) <- C.subtract row1.(i) (C.multiply sc row2.(i))
+    done;
+    set_row m r1 row1
+
+  (*************** End Helper Functions for Row Reduce ***************)
+
+  (* Returns the row reduced form of a matrix. Is not done in place, but creates
+   * a new matrix *)
+  let row_reduce (mat: matrix) : matrix =
+    let[@tailcall] rec row_reduce_h (n_row: int) (n_col: int) (mat2: matrix) : unit =
+      let ((num_row, _), _) = mat2 in
+      if (n_col = num_row + 1) then ()
+      else
+        let (_,col) = get_column mat2 n_col in
+        match find_max_col_index col n_row with
+        | None (* Column all 0s *) -> row_reduce_h n_row (n_col+1) mat2
+        | Some index ->
+          begin
+            swap_row mat2 index n_row;
+            let pivot = get_elt mat2 (n_row, n_col) in
+            scale_row mat2 (n_row) (C.divide C.one pivot);
+            for i = 1 to num_row do
+              if i <> n_row then sub_mult mat2 i n_row (get_elt mat2 (i,n_col))
+            done;
+            row_reduce_h (n_row+1) (n_col+1) mat2
+          end
+    in
+    (* Copies the matrix *)
+    let ((n,p),m) = mat in
+    let (dim,mat_cp) = empty n p in
+    for i = 0 to n - 1 do
+      for j = 0 to p - 1 do
+        mat_cp.(i).(j) <- m.(i).(j)
+      done;
+    done;
+    let _ = row_reduce_h 1 1 (dim,mat_cp) in (dim,mat_cp)
+
+  (*************** End Main Functions ***************)
+
+  (*************** Optional module functions ***************)
+
+  (* calculates the trace of a matrix *)
+  let trace (((n,p),m): matrix) : elt =
+    let rec build (elt: elt) (i: int) =
+      if i > -1 then
+        build (C.add m.(i).(i) elt) (i - 1)
+      else
+        elt in
+    if n = p then build C.zero (n - 1)
+    else raise ImproperDimensions
+
+  (* calculates the transpose of a matrix and retuns a new one *)
+  let transpose (((n,p),m): matrix) =
+    let (dim,m') = empty p n in
+    for i = 0 to n - 1 do
+      for j = 0 to p - 1 do
+        m'.(j).(i) <- m.(i).(j)
+      done;
+    done;
+    assert(dim = (p,n));
+    ((p,n),m')
+
+  (* Returns the inverse of a matrix. Uses a pretty simple algorithm *)
+  let inverse (mat: matrix) : matrix =
+    let ((n, p), _) = mat in
+    if n = p then
+      (* create augmented matrix *)
+      let augmented = empty n (2*n) in
+      for i = 1 to n do
+        let (dim,col) = get_column mat i in
+        let arr = Array.make n C.zero in
+        begin
+          assert(dim = n);
+          arr.(i-1) <- C.one;
+          set_column augmented i col;
+          set_column augmented (n + i) arr
+        end
+      done;
+      let augmented' = row_reduce augmented in
+      (* create the inverted matrix and fill in with appropriate values *)
+      let inverse = empty n n in
+      for i = 1 to n do
+        let (dim, col) = get_column augmented' (n + i) in
+        let _ = assert(dim = n) in
+        let _ = set_column inverse i col in
+        ()
+      done;
+      inverse
+    else
+      raise NonSquare
+
+  (***************** HELPER FUNCTIONS FOR DETERMINANT *****************)
+  (* creates an identity matrix of size n*)
+  let create_identity (n:int) : matrix =
+    let (dim,m) = empty n n in
+    for i = 0 to n - 1 do
+      m.(i).(i) <- C.one
+    done;
+    (dim,m)
+
+  (* Finds the index of the maximum value of an array *)
+  let find_max_index (arr: elt array) (start_index : int) : int =
+    let rec find_index (max_index: int) (curr_index: int) =
+      if curr_index = Array.length arr then max_index+1
+      else
+        match C.compare arr.(curr_index) arr.(max_index) with
+        | Equal | Less -> find_index max_index (curr_index + 1)
+        | Greater -> find_index curr_index (curr_index + 1) in
+    find_index (start_index - 1) start_index
+
+  (* Creates the pivoting matrix for A. Returns swqps. Adapted from
+   * http://rosettacode.org/wiki/LU_decomposition#Common_Lisp *)
+  let pivotize (((n,p),m): matrix) : matrix * int =
+    if n = p then
+      let swaps = ref 0 in
+      let pivot_mat = create_identity n in
+      for j = 1 to n do
+        let (_,col) = get_column ((n,p),m) j in
+        let max_index = find_max_index col j in
+        if max_index <> j then
+          (swaps := !swaps + 1; swap_row pivot_mat max_index j)
+        else ()
+      done;
+      (pivot_mat,!swaps)
+    else raise ImproperDimensions
+
+  (* decomposes a matrix into a lower triangualar, upper triangualar
+   * and a pivot matrix. It returns (L,U,P). Adapted from
+   * http://rosettacode.org/wiki/LU_decomposition#Common_Lisp *)
+  let lu_decomposition (((n,p),m): matrix) : (matrix*matrix*matrix)*int =
+    if n = p then
+      let mat = ((n,p),m) in
+      let lower, upper, (pivot,s) = empty n n, empty n n, pivotize mat in
+      let (_ ,l),(_ ,u), _ = lower,upper,pivot in
+      let ((_, _),mat') = mult pivot mat in
+      for j = 0 to n - 1 do
+        l.(j).(j) <- C.one;
+        for i = 0 to j do
+          let sum = ref C.zero in
+          for k = 0 to i - 1 do
+            sum := C.add (!sum) (C.multiply u.(k).(j) l.(i).(k))
+          done;
+          u.(i).(j) <- C.subtract mat'.(i).(j) (!sum)
+        done;
+        for i = j to n - 1 do
+          let sum = ref C.zero in
+          for k = 0 to j - 1 do
+            sum := C.add (!sum) (C.multiply u.(k).(j) l.(i).(k))
+          done;
+          let sub = C.subtract mat'.(i).(j) (!sum) in
+          l.(i).(j) <- C.divide sub u.(j).(j)
+        done;
+      done;
+      (lower,upper,pivot),s
+    else raise ImproperDimensions
+
+  (* Computes the determinant of a matrix *)
+  let determinant (m: matrix) : elt =
+    try
+      let ((n,p), _) = m in
+      if n = p then
+        let rec triangualar_det (a,mat) curr_index acc =
+          if curr_index < n then
+            let acc' = C.multiply mat.(curr_index).(curr_index) acc in
+            triangualar_det (a,mat) (curr_index + 1) acc'
+          else acc in
+        let ((dim1,l),(dim2,u), _),s = lu_decomposition m in
+        let det1, det2 = triangualar_det (dim1,l) 0 C.one,
+                         triangualar_det (dim2,u) 0 C.one in
+        if s mod 2 = 0 then C.multiply det1 det2
+        else C.subtract C.zero (C.multiply det1 det2)
+      else raise ImproperDimensions
+    with
+    | _ -> C.zero
+
+
+  (*************** Optional module functions ***************)
+
+
+end
diff --git a/script.it/shapes/matrix/MatrixI.ml b/script.it/shapes/matrix/MatrixI.ml
new file mode 100755
index 0000000..fbc4e21
--- /dev/null
+++ b/script.it/shapes/matrix/MatrixI.ml
@@ -0,0 +1,105 @@
+exception NonSquare
+exception ImproperDimensions
+
+module type MATRIX =
+sig
+
+  (******** TYPES ********)
+  type elt
+
+  type matrix
+
+  (* empty matrix of nxp dimensions *)
+  val empty : int -> int -> matrix
+
+  (* Takes a list of lists and converts that to a matrix *)
+  val from_list : (elt list list) -> matrix
+
+  val from_array: elt array array -> matrix
+
+  (******** OPERATIONS ON ONE MATRIX ********)
+  (* Takes in a matrix and returns its dimensions. ie, nxp *)
+  val get_dimensions : matrix -> (int * int)
+
+  (* get's the row of a matrix: Not zero-indexed. *)
+  val get_row : matrix -> int -> (int * elt array)
+
+  (* similar to get_row *)
+  val get_column: matrix -> int -> (int * elt array)
+
+  (* sets the row of a matrix in place! Not zero-index *)
+  val set_row: matrix -> int -> elt array -> unit
+
+  (* similar to set_row, but for a column *)
+  val set_column: matrix -> int -> elt array -> unit
+
+  (* gets the element at the specified index. *)
+  val get_elt: matrix -> (int * int) -> elt
+
+  (* sets the element at the specified index *)
+  val set_elt: matrix -> (int * int) -> elt -> unit
+
+  (* Scales every element in the matrix by another elt *)
+  val scale : matrix -> elt -> matrix
+
+
+  (******** MORE ADVANCED SINGLE MATRIX OPERATIONS ********)
+  (* Returns the row reduced form of a matrix *)
+  val row_reduce: matrix -> matrix
+  (* We will implement the algorithm found in the link above *)
+
+  (* Returns the inverse of a matrix *)
+  val inverse: matrix -> matrix
+
+  (*Transposes a matrix. If the input has dimensions m x n, the output will
+   * have dimensions n x m *)
+  val transpose: matrix -> matrix
+
+  (* Returns the trace of the matrix *)
+  val trace: matrix -> elt
+
+  (******** OPERATIONS ON TWO MATRICES ********)
+  (* Adds two matrices. They must have the same dimensions *)
+  val add : matrix -> matrix -> matrix
+
+  (* Multiplies two matrices. If the matrices have dimensions m x n and p x q, n
+   * and p must be equal, and the resulting matrix will have dimension m x q *)
+  val mult: matrix -> matrix -> matrix
+
+  (**** Other Library Functions ***)
+  (* Function to make over our matrices *)
+  val map : (elt -> elt) -> matrix -> matrix
+
+  (*val iter : (elt -> unit) -> matrix -> unit*)
+
+  (* Returns the LUP decomposition of a matrix *)
+  val lu_decomposition : matrix -> (matrix * matrix * matrix) * int
+
+  (* Returns the determinant of the matrix *)
+  val determinant: matrix -> elt
+
+  (************** Other Library Functions *************)
+  val iter : (elt -> unit) -> matrix -> unit
+
+  val iteri : (int -> int -> elt -> unit) -> matrix -> unit
+
+  (* folds over each row using base case u and function f *)
+  val reduce: ('a -> elt -> 'a) -> 'a -> matrix -> 'a
+
+  val fold_row: f:(elt array -> 'b) -> matrix -> 'b list
+
+  (********** Specific for Simplex Algorithm ***********)
+  (** All of the following functions will raise ImproperDimensions
+   * Exception if the matrix is not the right size for the operation
+   **)
+
+  (* Scales a row *)
+  val scale_row: matrix -> int -> elt -> unit
+
+  (* Swaps two rows *)
+  val swap_row: matrix -> int -> int -> unit
+
+  (* Subtracts a multiple of one row (the 2nd int) from another (the 1st int) *)
+  val sub_mult: matrix -> int -> int -> elt -> unit
+
+end
diff --git a/script.it/shapes/matrix/Order.ml b/script.it/shapes/matrix/Order.ml
new file mode 100755
index 0000000..5f2aa22
--- /dev/null
+++ b/script.it/shapes/matrix/Order.ml
@@ -0,0 +1,2 @@
+(* Defines a general ordering type *)
+type order = Equal | Less | Greater
diff --git a/script.it/shapes/matrix/dune b/script.it/shapes/matrix/dune
new file mode 100755
index 0000000..1c0cab6
--- /dev/null
+++ b/script.it/shapes/matrix/dune
@@ -0,0 +1,3 @@
+(library
+ (name matrix)
+)
diff --git a/script.it/shapes/tools/dune b/script.it/shapes/tools/dune
new file mode 100755
index 0000000..a2c3fdb
--- /dev/null
+++ b/script.it/shapes/tools/dune
@@ -0,0 +1,6 @@
+(library
+ (name tools)
+ (libraries 
+   gg
+   )
+  )
diff --git a/script.it/shapes/tools/utils.ml b/script.it/shapes/tools/utils.ml
new file mode 100755
index 0000000..b8a473f
--- /dev/null
+++ b/script.it/shapes/tools/utils.ml
@@ -0,0 +1,63 @@
+open Gg.V2
+
+let norm_angle vect =
+  mod_float
+    ((angle vect) +. Gg.Float.two_pi)
+    Gg.Float.two_pi
+
+
+let intersection
+  : (Gg.v2 * Gg.v2) -> (Gg.v2 * Gg.v2) -> Gg.v2 option
+  = fun (p0, p1) (p2, p3) ->
+    let i = p1 - p0
+    and j = p3 - p2 in
+
+    let d = (x i *. y j) -. (y i *. x j) in
+    if Float.( (abs d) <= 0.01 ) then
+      None
+    else
+      let m = ((x i) *. (y p0)
+               -. (x i) *. (y p2)
+               -. (y i) *. (x p0)
+               +. (y i) *. (x p2)) /. d in
+      Some (p2 + m * j)
+      (*
+      let k = ((x j) *. (y p0)
+               -. (x j) *. (y p2)
+               -. (y j) *. (x p0)
+               +. (y j) *. (x p2)) /. d in
+      Some (p0 + k * i)
+      *)
+
+
+let center
+  : Gg.v2 -> Gg.v2 -> Gg.v2 -> Gg.v2 option
+  = fun p0 p1 p2 ->
+    (* deltas *)
+    let d1 = p1 - p0
+    and d2 = p2 - p1 in
+
+    (* perpendiculars *)
+    let d1p = ortho d1
+    and d2p = ortho d2 in
+
+    (* Chord midpoints *)
+    let m1 = half (p0 + p1)
+    and m2 = half (p1 + p2) in
+
+    (* midpoint offsets *)
+    let m1n = m1 + d1p
+    and m2n = m2 + d2p in
+
+    intersection (m1, m1n) (m2, m2n)
+
+let rotate
+  : Gg.v2 -> float -> Gg.v2
+  = fun p0 theta ->
+    let r = x (to_polar p0) in
+    of_polar (v r theta)
+
+let equal_point
+  : float -> Gg.v2 -> Gg.v2 -> bool
+  = fun eps p0 p1 ->
+    Gg.V2.equal_f (fun v0 v1 ->  (Float.abs (v1 -. v0)) <= eps ) p0 p1
diff --git a/script.it/shapes/tools/utils.mli b/script.it/shapes/tools/utils.mli
new file mode 100755
index 0000000..4e12906
--- /dev/null
+++ b/script.it/shapes/tools/utils.mli
@@ -0,0 +1,21 @@
+(** Return a normalize angle *)
+val norm_angle 
+  : Gg.v2 -> float
+
+(** return the interesction for two segments *)
+val intersection
+  : (Gg.v2 * Gg.v2) -> (Gg.v2 * Gg.v2) -> Gg.v2 option
+
+(** Return the center of the cercle for three points  
+    None if the point cannot be evaluated
+*)
+val center
+  : Gg.v2 -> Gg.v2 -> Gg.v2 -> Gg.v2 option
+
+(** Rotate the vector by the given angle *)
+val rotate
+  : Gg.v2 -> float -> Gg.v2
+
+(** Test equality between two points *)
+val equal_point
+  : float -> Gg.v2 -> Gg.v2 -> bool