Monotone Cubic Interpolation

(20090728 omission corrected.) The Wikipedia article on monotone cubic interpolation is, regrettably, not of very good quality. The idea is to precondition the derivatives (control points) for the cubic Hermite splines interpolating a monotonic data set so that the resulting piecewise cubic is also monotonic. The algorithm given is extracted from a paper which is not publicly available online. Furthermore, the only example given is in the form of a graphic, without a data set or a resulting array of derivatives, which would have been great to validate an implementation. Even so, I managed to convince myself by visually inspecting a number of trials that the following should do the job.

The first step in the algorithm is to compute the tangents using three-point differences:

let hermite_tangents x y =
  let n = Array.length x in
  if Array.length y != n then invalid_arg "hermite_tangents" else
  let m = Array.make n 0. in
  let d = ref ((y.(1) -. y.(0)) /. (x.(1) -. x.(0))) in
  m.(0) <- !d;
  for i = 1 to n - 2 do
    let d' = (y.(i+1) -. y.(i)) /. (x.(i+1) -. x.(i)) in
    m.(i) <- 0.5 *. (!d +. d');
    d := d'
  done;
  m.(n - 1) <- !d;
  m

Note that I've used a mutable reference to hold the value of the previous divided difference. Note also that this is the standard calculation for Hermite interpolation; this routine can be used as is. Now steps 2 to 5 precondition the tangents to achieve the desired monotonic interpolants:

let monotone_hermite_tangents x y =
  let n = Array.length x in
  let m = hermite_tangents x y in
  for i = 0 to n - 2 do
    let d = (y.(i+1) -. y.(i)) /. (x.(i+1) -. x.(i)) in
    if d = 0. then begin
      m.(i  ) <- 0.;
      m.(i+1) <- 0.
    end else
      let a = m.(i  ) /. d
      and b = m.(i+1) /. d in
      let s = a *. a +. b *. b in
      if s > 9. then begin
        let t = 3. *. d /. sqrt s in
        m.(i  ) <- t *. a;
        m.(i+1) <- t *. b
      end
  done;
  m

In the article's description, the Δ_k appear to be array elements; since no step uses previous values of the variable, it is effectively a loop local. Also, I took the liberty to eliminate some (obvious) common subexpressions. That done, the interpolation code is straightforward, and standard. First, given a point e, I need some way to determine which of the intervals [x_i…x_i+1] contains it. This naturally calls for a binary search:

let bin_search e x =
  let n = Array.length x in
  let l = ref 0
  and h = ref n in
  while !l + 1 != !h do
    let m = (!l + !h) / 2 in
    if e < x.(m) then h := m else l := m
  done;
  !l

For values outside the data range, I return the first or last value, as appropriate. Otherwise, as the page on cubic interpolation explains, the interval containing the point e is normalized to the unit interval [0…1], interpolated there using the Hermite basis, and rescaled to the original interval size via h:

let hermite_interpolation x y m e =
  let n = Array.length x in
  if e <= x.(0)   then y.(0) else
  if e >= x.(n-1) then y.(n-1) else
  let i = bin_search e x in
  let l = x.(i) in
  let h = x.(i+1) -. l in
  let t = (e -. l) /. h in
  let u = t *. t in
  let h00 =  ( 2. *. t -. 3.) *. u +. 1.
  and h10 = ((       t -. 2.) *. t +. 1.) *. t
  and h01 =  (-2. *. t +. 3.) *. u
  and h11 =  (       t -. 1.) *. u in
  y.(i) *. h00 +. y.(i+1) *. h01 +. h *. (m.(i) *. h10 +. m.(i+1) *. h11)