Directory|Libraries|E-learning|Directions|Jobs|Calendar|U tv|Contact
En
UMONS>Faculty of Science>Department of Mathematics>Numerical Analysis Team>Software>OCaml>Circle Packing

Packing circles into a rectangle

Given a rectangle and a number \(n\), the problem of circle packing is to choose the center and radius of \(n\) circles so that

Mathematically, the problem translates as follows. Let \(w\) be the rectangle width, \(h\) be its height, and let \((x_i,y_i)\) and \(r_i\), \(i=1,\dots,n\), denote the center and the radius of the circles. The goal is to maximize the function \[ (x_1,\dots,x_n, y_1,\dots,y_n, r_1,\dots,r_n) \mapsto \frac{\pi}{wh} \, \sum_{i=1}^n r_i^2 \] under the constraints that all the circles are inside the rectangle: \[ x_i + r_i \le w,\qquad x_i - r_i \ge 0,\qquad y_i + r_i \le h,\qquad y_i - r_i \ge 0,\qquad (i=1,\dots,n) \] and that they do not overlap: \[ \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2} \ge r_i + r_j,\qquad \text{for all } i, j = 1,\dotsc,n,\text{ such that } i \lt j. \]

Following this blog post, we will use a penalization formulation. We will thus maximize the following function \(f\) w.r.t. the parameters \(z = (x_1,\dots,x_n, y_1,\dots,y_n, r_1,\dots,r_n)\): \[ \begin{split} f(z) := \tfrac{1}{2} \sum_{i=1}^n r_i^2 - \tfrac{1}{2} \eta \Biggl( &\sum_{i=1}^n \bigl((x_i + r_i - w)^+\bigr)^2 + \sum_{i=1}^n \bigl((r_i - x_i)^+\bigr)^2 \\ &+ \sum_{i=1}^n \bigl((y_i + r_i - h)^+\bigr)^2 + \sum_{i=1}^n \bigl((r_i - y_i)^+\bigr)^2 \\ &+ \sum_{i=1}^n \sum_{j \lt i} \bigl( (r_i + r_j - \sqrt{(x_i - x_j)^2 + (y_i - y_j)^2})^+ \bigr)^2 \Biggr) \end{split} \] where \(\eta \gt 0\) is a penalization parameter and \(\xi^+ := \max\{\xi, 0\}\). Let us remark that, if all points \((x_i, y_i)\) are inside the rectangle, the vector \(z\) where all \(r_i = 0\) is a critical point of \(f\) — corresponding to its minimum — and thus cannot be used as a starting point (no descent direction).

Implementation

The implementation uses Lbfgs.F.max from the OCaml binding of the L-BFGS-B routines by J. Nocedal & al. The API is very simple thanks to OCaml optional arguments. The graphics are done thanks to Archimedes.

The main function to write is the penalized objective function f_df. Lbfgs.F.max requires that f_df z g returns the value of the function \(f\) at point z and stores its gradient in g. The quantities ofsy and ofsr are the offsets of the \((y_1,\dots,y_n)\) and \((r_1,\dots,r_n)\) components which allows to define sub-arrays for easier access (the content of the sub-arrays x, y and r is shared with the main array z). The code is a straightforward implementation of the formula for \(f\) — and for \(\nabla f\). To understand it, it suffices to know that

  let f_df (z: vec) (grad: vec) =
    let x = Array1.sub z 1 n
    and y = Array1.sub z ofsy n
    and r = Array1.sub z ofsr n in     (* z = (x,y,r), content shared *)
    let gx = Array1.sub grad 1 n       (* ∇_x *)
    and gy = Array1.sub grad ofsy n    (* ∇_y *)
    and gr = Array1.sub grad ofsr n in (* ∇_r *)
    let main_obj = Vec.sqr_nrm2 r in   (* ∑ r²_i *)
    Array1.fill grad 0.;
    ignore(copy ~y:gr r);
    (* x_i + r_i ≤ w *)
    let tmp = copy x in
    axpy tmp ~x:r;  sub_const tmp w;  pos tmp; (* tmp.{i} = (x_i + r_i - w)⁺ *)
    let c1 = Vec.sqr_nrm2 tmp in
    axpy gx ~alpha:(-. eta) ~x:tmp;
    axpy gr ~alpha:(-. eta) ~x:tmp;
    (* x_i - r_i ≥ 0 *)
    ignore(copy ~y:tmp r);
    axpy tmp ~alpha:(-1.) ~x;  pos tmp; (* tmp.{i} = (r_i - x_i)⁺ *)
    let c2 = Vec.sqr_nrm2 tmp in
    axpy gr ~alpha:(-. eta) ~x:tmp;
    axpy gx ~alpha:eta      ~x:tmp;
    (* y_i + r_i ≤ h *)
    ignore(copy ~y:tmp y);
    axpy tmp ~x:r;  sub_const tmp h;  pos tmp; (* tmp.{i} = (y_i + r_i - h)⁺ *)
    let c3 = Vec.sqr_nrm2 tmp in
    axpy gy ~alpha:(-. eta) ~x:tmp;
    axpy gr ~alpha:(-. eta) ~x:tmp;
    (* y_i - r_i ≥ 0 *)
    ignore(copy ~y:tmp r);
    axpy tmp ~alpha:(-1.) ~x:y;  pos tmp; (* tmp.{i} = (r_i - y_i)⁺ *)
    let c4 = Vec.sqr_nrm2 tmp in
    axpy gr ~alpha:(-. eta) ~x:tmp;
    axpy gy ~alpha:eta      ~x:tmp;
    (* √((x_i - x_j)² + (y_i - y_j)²) ≥ r_i + r_j *)
    let c5 = ref 0. in
    for i = 1 to n do
      for j = 1 to i-1 do
        let dx = x.{i} -. x.{j}
        and dy = y.{i} -. y.{j} in
        let dij = sqrt(dx *. dx +. dy *. dy) in
        let c = r.{i} +. r.{j} -. dij in
        if c > 0. then ( (* constraint not satisfied *)
          c5 := !c5 +. c *. c;
          let m = eta *. c /. dij in
          gx.{i} <- gx.{i} +. m *. dx;
          gx.{j} <- gx.{j} -. m *. dx;
          gy.{i} <- gy.{i} +. m *. dy;
          gy.{j} <- gy.{j} -. m *. dy;
          gr.{i} <- gr.{i} -. eta *. c;
          gr.{j} <- gr.{j} -. eta *. c;
        )
      done
    done;
    0.5 *. (main_obj -. eta *. (c1 +. c2 +. c3 +. c4 +. !c5))
For the initial vector z, the points are distributed randomly in the rectangle and a small radius is selected (the initial circles may overlap).

The drawing prodedure is simple. One requests Archimedes to open a viewport drawing on a backend bk, then one creates a path p consisting of all the circles, and finally one draws the path p.

let draw w h (z: vec) bk =
  let vp = A.init bk ~w:800. ~h:800. in
  A.xrange vp 0. w;
  A.yrange vp 0. h;
  A.Axes.box vp;
  A.set_color vp A.Color.red;
  let n = Vec.dim z / 3 in
  let p = A.Path.make() in
  for i = 1 to n do
    let xi = z.{i} and yi = z.{n + i} and ri = z.{2 * n + i} in
    A.Path.move_to p (xi +. ri) yi;
    A.Path.arc p ri 0. (2. *. pi);
  done;
  A.Viewport.stroke vp `Data ~path:p;
  A.close vp

Unlike the implementation in R, this code handles 500 circles in a matter of seconds. Unfortunately, some of the final circles have radius 0. The following figures show, left, the initial vector z and, right, the same vector after optimization (click them to enlarge).

Full Code

open Printf
open Bigarray
open Lacaml.D
module A = Archimedes

let pi = 4. *. atan 1.

(** [pos v] transform [v] into the vector [max v.{i} 0]. *)
let pos v =
  for i = 1 to Vec.dim v do if v.{i} < 0. then v.{i} <- 0. done

let add_const (v: vec) c =
  for i = 1 to Vec.dim v do v.{i} <- v.{i} +. c done

let sub_const (v: vec) c =
  for i = 1 to Vec.dim v do v.{i} <- v.{i} -. c done

let packing_density ~w ~h (z: vec) =
  let n = Vec.dim z / 3 in
  pi /. (w *. h) *. Vec.sqr_nrm2 z ~n ~ofsx:(2 * n + 1)

(** [pack_circles ~w ~h z] modifies the initial guess [z] until an
    optimal sphere packing is reached.

   [z] packs the x_i, y_i and r_i in the following way z =
   (x₁,...,x_n, y₁,....,y_n, r₁,...,r_n). *)
let pack_circles ?(eta=256.) ~w ~h (z: vec) =
  let n = Vec.dim z / 3 in
  let ofsy = n + 1 and ofsr = 2 * n + 1 in
  let f_df (z: vec) (grad: vec) =
    let x = Array1.sub z 1 n
    and y = Array1.sub z ofsy n
    and r = Array1.sub z ofsr n in     (* z = (x,y,r), content shared *)
    let gx = Array1.sub grad 1 n       (* ∇_x *)
    and gy = Array1.sub grad ofsy n    (* ∇_y *)
    and gr = Array1.sub grad ofsr n in (* ∇_r *)
    let main_obj = Vec.sqr_nrm2 r in   (* ∑ r²_i *)
    Array1.fill grad 0.;
    ignore(copy ~y:gr r);
    (* x_i + r_i ≤ w *)
    let tmp = copy x in
    axpy tmp ~x:r;  sub_const tmp w;  pos tmp; (* tmp.{i} = (x_i + r_i - w)⁺ *)
    let c1 = Vec.sqr_nrm2 tmp in
    axpy gx ~alpha:(-. eta) ~x:tmp;
    axpy gr ~alpha:(-. eta) ~x:tmp;
    (* x_i - r_i ≥ 0 *)
    ignore(copy ~y:tmp r);
    axpy tmp ~alpha:(-1.) ~x;  pos tmp; (* tmp.{i} = (r_i - x_i)⁺ *)
    let c2 = Vec.sqr_nrm2 tmp in
    axpy gr ~alpha:(-. eta) ~x:tmp;
    axpy gx ~alpha:eta      ~x:tmp;
    (* y_i + r_i ≤ h *)
    ignore(copy ~y:tmp y);
    axpy tmp ~x:r;  sub_const tmp h;  pos tmp; (* tmp.{i} = (y_i + r_i - h)⁺ *)
    let c3 = Vec.sqr_nrm2 tmp in
    axpy gy ~alpha:(-. eta) ~x:tmp;
    axpy gr ~alpha:(-. eta) ~x:tmp;
    (* y_i - r_i ≥ 0 *)
    ignore(copy ~y:tmp r);
    axpy tmp ~alpha:(-1.) ~x:y;  pos tmp; (* tmp.{i} = (r_i - y_i)⁺ *)
    let c4 = Vec.sqr_nrm2 tmp in
    axpy gr ~alpha:(-. eta) ~x:tmp;
    axpy gy ~alpha:eta      ~x:tmp;
    (* √((x_i - x_j)² + (y_i - y_j)²) ≥ r_i + r_j *)
    let c5 = ref 0. in
    for i = 1 to n do
      for j = 1 to i-1 do
        let dx = x.{i} -. x.{j}
        and dy = y.{i} -. y.{j} in
        let dij = sqrt(dx *. dx +. dy *. dy) in
        let c = r.{i} +. r.{j} -. dij in
        if c > 0. then ( (* constraint not satisfied *)
          c5 := !c5 +. c *. c;
          let m = eta *. c /. dij in
          gx.{i} <- gx.{i} +. m *. dx;
          gx.{j} <- gx.{j} -. m *. dx;
          gy.{i} <- gy.{i} +. m *. dy;
          gy.{j} <- gy.{j} -. m *. dy;
          gr.{i} <- gr.{i} -. eta *. c;
          gr.{j} <- gr.{j} -. eta *. c;
        )
      done
    done;
    0.5 *. (main_obj -. eta *. (c1 +. c2 +. c3 +. c4 +. !c5))
  in
  (* Bounds *)
  let l = Vec.make0 (Vec.dim z) in
  ignore(Lbfgs.F.max f_df z ~l ~factr:1e5)


let draw w h (z: vec) bk =
  let vp = A.init bk ~w:800. ~h:800. in
  A.xrange vp 0. w;
  A.yrange vp 0. h;
  A.Axes.box vp;
  A.set_color vp A.Color.red;
  let n = Vec.dim z / 3 in
  let p = A.Path.make() in
  for i = 1 to n do
    let xi = z.{i} and yi = z.{n + i} and ri = z.{2 * n + i} in
    A.Path.move_to p (xi +. ri) yi;
    A.Path.arc p ri 0. (2. *. pi);
  done;
  A.Viewport.stroke vp `Data ~path:p;
  A.close vp

let pack_and_draw_circles ?(w=1.) ?(h=1.) z ?(ext=".png") fbase =
  let fbase = Filename.concat Filename.temp_dir_name fbase in
  draw w h z (A.backend_of_filename(fbase ^ "0" ^ ext));
  pack_circles ~w ~h z;
  let n = Vec.dim z / 3 in
  let nr = Vec.fold (fun n r -> if r > 0. then n+1 else n) 0 z ~ofsx:(2*n+1) in
  printf "n = %i (not empty = %i): packing density = %g%% (\"%s*%s\")\n%!"
         n nr (100. *. packing_density ~w ~h z) fbase ext;
  draw w h z (A.backend_of_filename(fbase ^ "1" ^ ext))

let random_guess ?(w=1.) ?(h=1.) n =
  let r = 0.01 *. min w h in
  let z = Vec.make (3 * n) r in
  for i = 1 to n do
    z.{i} <- r +. Random.float (w -. 2. *. r);
    z.{n + i} <- r +. Random.float (h -. 2. *. r);
  done;
  z

let () =
  Random.self_init();
  pack_and_draw_circles (random_guess 500) "circle"

If this code is put in a file circle_pack.ml, it can be compiled with: ocamlfind ocamlopt -annot -o circle_pack -package lacaml,lbfgs,archimedes -linkpkg circle_pack.ml

Powered by MathJax