2008-06-06

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 [xixi+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)

3 comments:

Jeremie said...

Hi Matias,

I hope you'll read this comment on a old post. My name is Jeremie, I was looking for some information on the Monotone cubic interpolation and the best I could find was your website.

I'm trying to apply this algorythm to a 3D spline.
I've been able to reproduce the Catmull-Rom and the Kochanek–Bartels spline but this one seems way over my mathematical skills. I'll be very glad if you can contact me so I can ask you a few question to see how I could get such curve in a 3D software.

my email is geerem-at-hotmail.com

Cheers,
Jeremie

Håvard Berland said...

Thanks for reviewing the Wikipedia article I submitted some years ago, out of frustration of only having a non-publicly available paper as a source. If you feel that the text there could be improved to your liking, please do so.

There is now a link to a GPL-implementation in C++, which might serve your purpose, that particular code was used to produce the data in the plot.

Regards

Matías Giovannini said...

Håvard,
It was ungracious of me to critizice the article, especially since anything is always better than nothing. I apologize. I also realize belatedly that supplemental material is not really appropriate for a Wikipedia article.