Interpolation Methods for Color Values

[matt7]

I recently posted to the forum a gradient drawing library (here) and a gradient editing program (here). The gradient library interpolates colors in the RGB space, which can be undesirable at times because interpolating R, G, and B values independently can often cause unexpected or unintended intermediate colors in the gradient between defined gradient stops (such as red to cyan creating gray in the middle). There are alternatives (see this article: The Secrets of Colour Interpolation), but the color space that is being interpolated through is not what this post is going to address.

Instead, I plan to keep interpolation in the RGB space (at least for now), but I want to improve and provide more options than the 5th order polynomial easing function that is currently implemented. I finally have a working cubic interpolation function (it was simpler than I thought!) which I wanted to showcase below. The real trick for using cubic interpolation for colors though is that you want to prevent overshoot outside of the range [0-1] and you want to preserve monotonicity where possible. With the help of Wikipedia (Cubic Hermite splines and Monotone cubic interpolation), I finally got that working as well.

I have not yet updated the gradient library, but wanted to share the quick and dirty script I put together that compares a variety of interpolation methods. Compared in this script are the following interpolation curves:

• Linear Interpolation
• Cosine Interpolation
• Parameterized 5th Order Polynomial Interpolation (select acceleration values from 0 to 30)
• Cubic Hermite Spline (Monotonic)
• Cubic Hermite Spline

Here is a comparison between linear, cosine, and the 5th order polynomial:



Again, the 5th order polynomial is currently what the gradient library uses. You can see where it isn't as smooth when values are not oscillating, since the slope at every gradient stop is 0, causing the curve to have to level out to 0 when it could just pass through the point. This is the main problem that cubic interpolation solves.

Here is interpolation over the same points using cubic hermite splines (one with tangents estimated using finite difference and the other with tangents adjusted for preserving monotonicity):



The "Cubic Monotone Splines" are what I'm most excited about implementing, because I think they will look the smoothest and most natural. Hopefully this method won't have too big of an impact on the gradient drawing speed.

Code: Select all

DEF GET_TANGENTS_FINITE_DIFFERENCE (x(), p())

  GET DIM x XSIZE n
  GET DIM p XSIZE np

  IF n <> np THEN RETURN -1
  DIM m(n)

  FOR k = 0 TO n-1
    IF k = 0 THEN
      m(k) = (p(k+1) - p(k)) / (x(k+1) - x(k))
    ELSE ! IF k = n-1 THEN
      m(k) = (p(k) - p(k-1)) / (x(k) - x(k-1))
    ELSE
      mL = (p(k) - p(k-1)) / (x(k) - x(k-1))
      mR = (p(k+1) - p(k)) / (x(k+1) - x(k))
      m(k) = (mL + mR)/2
    END IF ! END IF
  NEXT k

END DEF

DEF GET_TANGENTS_MONOTONIC (x(), p())

  GET DIM x XSIZE n
  GET DIM p XSIZE np

  IF n <> np THEN RETURN -1
  DIM m(n)

  ' Start with average of secant line slopes
  GET_TANGENTS_FINITE_DIFFERENCE(x, p)
  COPY1(GET_TANGENTS_FINITE_DIFFERENCE.m, m)

  ' Adjust tangents to preserve curve mononocity
  FOR k = 1 TO n-2
    mL = (p(k) - p(k-1)) / (x(k) - x(k-1))
    mR = (p(k+1) - p(k)) / (x(k+1) - x(k))
    IF SIGN(mL) <> SIGN(mR) THEN m(k) = 0
  NEXT k

  FOR k = 0 TO n-2
    mR = (p(k+1) - p(k)) / (x(k+1) - x(k))
    IF mR = 0 THEN
      m(k) = 0
      m(k+1) = 0
    ELSE
      ak = m(k)/mR
      bk = m(k+1)/mR
      IF ak < 0 OR bk < 0 THEN
        m(k) = 0
      ELSE ! IF (ak*ak + bk*bk) > 9 THEN
        ' Prevent overshoot
        tk = 3/(SQRT(ak*ak + bk*bk))
        m(k) = tk*ak*mR
        m(k+1) = tk*bk*mR
      END IF ! END IF
    END IF
  NEXT k

END DEF

DEF CUBIC_HERMITE_SPLINE (t, x1, x2, p1, p2, m1, m2)

  dx = x2 -x1
  t2 = t*t
  t3 = t2*t

  ' Hermite basis functions
  h00 = 2*t3 - 3*t2 + 1
  h10 = t3 - 2*t2 + t
  h01 = -2*t3 + 3*t2
  h11 = t3 - t2

  RETURN h00*p1 + h10*dx*m1 + h01*p2 + h11*dx*m2

END DEF


'==============================================================

.pi = 2*ACOS(0)

DEF LERP (t, a, b) = a + t*(b-a)

DEF COSERP (t, a, b) = LERP((1 - COS(.pi*t))/2, a, b)

DEF EASE_POLY5_A (t, a) = t*t*(t*(t*(t*(6-a)+(5*a-30)/2)+2*(5-a))+a/2)

DEF COPY1 (from(), to())

  GET DIM from XSIZE n
  GET DIM to XSIZE m

  IF n <> m THEN RETURN -1

  FOR i = 0 TO n-1
    to(i) = from(i)
  NEXT i

END DEF

'==============================================================

GRAPHICS CLEAR 1,1,1

SET ORIENTATION PORTRAIT
GET SCREEN SIZE scrW, scrH

xMin = scrW*0.1
xMax = scrW*0.9
yMin = xMin
yMax = xMax

sz = xMax - xMin

n = 9
DIM xk(n), pk(n), mk(n)

xk(0) = sz*0.05
pk(0) = sz*0.2

xk(1) = sz*0.19
pk(1) = sz*0.5

xk(2) = sz*0.31
pk(2) = sz*0.5

xk(3) = sz*0.5
pk(3) = sz*1

xk(4) = sz*0.6
pk(4) = sz*0.9

xk(5) = sz*0.65
pk(5) = sz*0.4

xk(6) = sz*0.7
pk(6) = sz*0.1

xk(7) = sz*0.8
pk(7) = sz*0.5

xk(8) = sz*0.95
pk(8) = sz*0.3

GET_TANGENTS_FINITE_DIFFERENCE(xk, pk)
COPY1(GET_TANGENTS_FINITE_DIFFERENCE.m, mk)
mk(0) = 0
mk(n-1) = 0

DRAW COLOR 0, 0, 0
fs = FONT_SIZE()

SPRITE "LERP" BEGIN scrW,scrH
FILL COLOR .5, .5, 1
FOR k = 0 TO n-1
  x = xMin + xk(k)
  y = yMin + sz - pk(k)
  FILL CIRCLE x,y SIZE 3
NEXT k
FILL COLOR 0,0,0
FOR x = 0 TO sz STEP 0.25
  IF x < xk(0) THEN
    y = pk(0)
  ELSE ! IF x > xk(n-1) THEN
    y = pk(n-1)
  ELSE
    FOR k = 0 TO n-2
      IF x < xk(k+1) THEN BREAK k
    NEXT k
    t = (x - xk(k)) / (xk(k+1) - xk(k))
    y = LERP(t, pk(k), pk(k+1))
  END IF ! END IF
  y = sz - y
  FILL CIRCLE xMin+x, yMin+y SIZE 0.5
NEXT x
DRAW TEXT "Linear Interpolation" AT xMin, yMax
SPRITE END

SPRITE "COSERP" BEGIN scrW,scrH
FILL COLOR .5, .5, 1
FOR k = 0 TO n-1
  x = xMin + xk(k)
  y = yMin + sz - pk(k)
  FILL CIRCLE x,y SIZE 3
NEXT k
FILL COLOR 0,0,0
FOR x = 0 TO sz STEP 0.25
  IF x < xk(0) THEN
    y = pk(0)
  ELSE ! IF x > xk(n-1) THEN
    y = pk(n-1)
  ELSE
    FOR k = 0 TO n-2
      IF x < xk(k+1) THEN BREAK k
    NEXT k
    t = (x - xk(k)) / (xk(k+1) - xk(k))
    y = COSERP(t, pk(k), pk(k+1))
  END IF ! END IF
  y = sz - y
  FILL CIRCLE xMin+x, yMin+y SIZE 0.5
NEXT x
DRAW TEXT "Cosine Interpolation" AT xMin, yMax
SPRITE END

ea = 0
SPRITE "POLY5_0" BEGIN scrW,scrH
FILL COLOR .5, .5, 1
FOR k = 0 TO n-1
  x = xMin + xk(k)
  y = yMin + sz - pk(k)
  FILL CIRCLE x,y SIZE 3
NEXT k
FILL COLOR 0,0,0
FOR x = 0 TO sz STEP 0.25
  IF x < xk(0) THEN
    y = pk(0)
  ELSE ! IF x > xk(n-1) THEN
    y = pk(n-1)
  ELSE
    FOR k = 0 TO n-2
      IF x < xk(k+1) THEN BREAK k
    NEXT k
    t = (x - xk(k)) / (xk(k+1) - xk(k))
    t = EASE_POLY5_A(t, ea)
    y = LERP(t, pk(k), pk(k+1))
  END IF ! END IF
  y = sz - y
  FILL CIRCLE xMin+x, yMin+y SIZE 0.5
NEXT x
DRAW TEXT "5th Order Polynomial" AT xMin, yMax
DRAW TEXT "(ea="&ea&")" AT xMin, yMax+ fs*1.2
SPRITE END

ea = 6
SPRITE "POLY5_6" BEGIN scrW,scrH
FILL COLOR .5, .5, 1
FOR k = 0 TO n-1
  x = xMin + xk(k)
  y = yMin + sz - pk(k)
  FILL CIRCLE x,y SIZE 3
NEXT k
FILL COLOR 0,0,0
FOR x = 0 TO sz STEP 0.25
  IF x < xk(0) THEN
    y = pk(0)
  ELSE ! IF x > xk(n-1) THEN
    y = pk(n-1)
  ELSE
    FOR k = 0 TO n-2
      IF x < xk(k+1) THEN BREAK k
    NEXT k
    t = (x - xk(k)) / (xk(k+1) - xk(k))
    t = EASE_POLY5_A(t, ea)
    y = LERP(t, pk(k), pk(k+1))
  END IF ! END IF
  y = sz - y
  FILL CIRCLE xMin+x, yMin+y SIZE 0.5
NEXT x
DRAW TEXT "5th Order Polynomial" AT xMin, yMax
DRAW TEXT "(ea="&ea&")" AT xMin, yMax+ fs*1.2
SPRITE END

ea = 14
SPRITE "POLY5_14" BEGIN scrW,scrH
FILL COLOR .5, .5, 1
FOR k = 0 TO n-1
  x = xMin + xk(k)
  y = yMin + sz - pk(k)
  FILL CIRCLE x,y SIZE 3
NEXT k
FILL COLOR 0,0,0
FOR x = 0 TO sz STEP 0.25
  IF x < xk(0) THEN
    y = pk(0)
  ELSE ! IF x > xk(n-1) THEN
    y = pk(n-1)
  ELSE
    FOR k = 0 TO n-2
      IF x < xk(k+1) THEN BREAK k
    NEXT k
    t = (x - xk(k)) / (xk(k+1) - xk(k))
    t = EASE_POLY5_A(t, ea)
    y = LERP(t, pk(k), pk(k+1))
  END IF ! END IF
  y = sz - y
  FILL CIRCLE xMin+x, yMin+y SIZE 0.5
NEXT x
DRAW TEXT "5th Order Polynomial" AT xMin, yMax
DRAW TEXT "(ea="&ea&")" AT xMin, yMax+ fs*1.2
SPRITE END

ea = 24
SPRITE "POLY5_24" BEGIN scrW,scrH
FILL COLOR .5, .5, 1
FOR k = 0 TO n-1
  x = xMin + xk(k)
  y = yMin + sz - pk(k)
  FILL CIRCLE x,y SIZE 3
NEXT k
FILL COLOR 0,0,0
FOR x = 0 TO sz STEP 0.25
  IF x < xk(0) THEN
    y = pk(0)
  ELSE ! IF x > xk(n-1) THEN
    y = pk(n-1)
  ELSE
    FOR k = 0 TO n-2
      IF x < xk(k+1) THEN BREAK k
    NEXT k
    t = (x - xk(k)) / (xk(k+1) - xk(k))
    t = EASE_POLY5_A(t, ea)
    y = LERP(t, pk(k), pk(k+1))
  END IF ! END IF
  y = sz - y
  FILL CIRCLE xMin+x, yMin+y SIZE 0.5
NEXT x
DRAW TEXT "5th Order Polynomial" AT xMin, yMax
DRAW TEXT "(ea="&ea&")" AT xMin, yMax+ fs*1.2
SPRITE END

ea = 30
SPRITE "POLY5_30" BEGIN scrW,scrH
FILL COLOR .5, .5, 1
FOR k = 0 TO n-1
  x = xMin + xk(k)
  y = yMin + sz - pk(k)
  FILL CIRCLE x,y SIZE 3
NEXT k
FILL COLOR 0,0,0
FOR x = 0 TO sz STEP 0.25
  IF x < xk(0) THEN
    y = pk(0)
  ELSE ! IF x > xk(n-1) THEN
    y = pk(n-1)
  ELSE
    FOR k = 0 TO n-2
      IF x < xk(k+1) THEN BREAK k
    NEXT k
    t = (x - xk(k)) / (xk(k+1) - xk(k))
    t = EASE_POLY5_A(t, ea)
    y = LERP(t, pk(k), pk(k+1))
  END IF ! END IF
  y = sz - y
  FILL CIRCLE xMin+x, yMin+y SIZE 0.5
NEXT x
DRAW TEXT "5th Order Polynomial" AT xMin, yMax
DRAW TEXT "(ea="&ea&")" AT xMin, yMax+ fs*1.2
SPRITE END

SPRITE "CUBIC_HERMITE" BEGIN scrW,scrH
FILL COLOR .5, .5, 1
FOR k = 0 TO n-1
  x = xMin + xk(k)
  y = yMin + sz - pk(k)
  FILL CIRCLE x,y SIZE 3
NEXT k
FILL COLOR 0,0,0
FOR x = 0 TO sz STEP 0.25
  IF x < xk(0) THEN
    y = pk(0)
  ELSE ! IF x > xk(n-1) THEN
    y = pk(n-1)
  ELSE
    FOR k = 0 TO n-2
      IF x < xk(k+1) THEN BREAK k
    NEXT k
    t = (x - xk(k)) / (xk(k+1) - xk(k))
    y = CUBIC_HERMITE_SPLINE(t, xk(k), xk(k+1), pk(k), pk(k+1), mk(k), mk(k+1))
  END IF ! END IF
  y = sz - y
  FILL CIRCLE xMin+x, yMin+y SIZE 0.5
NEXT x
DRAW TEXT "Cubic Hermite Splines" AT xMin, yMax
SPRITE END

SPRITE "CUBIC_MONOTONIC" BEGIN scrW,scrH
FILL COLOR .5, .5, 1
FOR k = 0 TO n-1
  x = xMin + xk(k)
  y = yMin + sz - pk(k)
  FILL CIRCLE x,y SIZE 3
NEXT k

GET_TANGENTS_MONOTONIC(xk, pk)
COPY1(GET_TANGENTS_MONOTONIC.m, mk)
mk(0) = 0
mk(n-1) = 0

FILL COLOR 0,0,0
FOR x = 0 TO sz STEP 0.25
  IF x < xk(0) THEN
    y = pk(0)
  ELSE ! IF x > xk(n-1) THEN
    y = pk(n-1)
  ELSE
    FOR k = 0 TO n-2
      IF x < xk(k+1) THEN BREAK k
    NEXT k
    t = (x - xk(k)) / (xk(k+1) - xk(k))
    y = CUBIC_HERMITE_SPLINE(t, xk(k), xk(k+1), pk(k), pk(k+1), mk(k), mk(k+1))
  END IF ! END IF
  y = sz - y
  FILL CIRCLE xMin+x, yMin+y SIZE 0.5
NEXT x
DRAW TEXT "Cubic Monotone Splines" AT xMin, yMax
SPRITE END

num_methods = 9
DIM n$(num_methods)
n$(0) = "LERP"
n$(1) = "COSERP"
n$(2) = "POLY5_0"
n$(3) = "POLY5_6"
n$(4) = "POLY5_14"
n$(5) = "POLY5_24"
n$(6) = "POLY5_30"
n$(7) = "CUBIC_MONOTONIC"
n$(8) = "CUBIC_HERMITE"

GRAPHICS

DRAW COLOR 0.5, 0.5, 0.5
DRAW RECT xMin,yMin TO xMax,yMax
DRAW FONT SIZE fs*0.6
DRAW TEXT "Tap to cycle interpolation methods" AT xMin, yMax*1.5
DRAW TEXT "Use two touches to exit" AT xMin, yMax*1.5 + fs*1.2

index = 0
SPRITE n$(index) SHOW

WHILE TOUCH_X(1) = -1
  IF NOT sameTouch AND TOUCH_X(0) > -1 THEN
    sameTouch = 1
    prevIndex = index
    index = (index + 1) % (num_methods)

    SPRITE n$(index) SHOW
    SPRITE n$(prevIndex) HIDE
  ELSE ! IF TOUCH_X(0) = -1 THEN
    sameTouch = 0
  END IF ! END IF

  SLOWDOWN
END WHILE
  

[Henko]

Nice post!
I like the monotone cubic spline as well.
No worries about the speed / needed time

IMG_1540.PNG (367.49 KiB) Viewed 2637 times

Code: Select all

option base 1

n=9 ! ni=20 ! nret=(n-1)*ni+1
dim x(n),y(n),xret(nret),yret(nret)

for i=1 to n ! read x(i),y(i) ! next i
data -7,-5, -5,9, -1,8, .5,6, 2,6, 3,10, 4,5.5, 6,5, 8,2
time reset
m_c_spline(n,x,y,ni,xret,yret)
tt=time()
graph_it(n,x,y,nret,xret,yret)
draw text "input  : "&n&" points" at 270,30
draw text "output : "&nret&" interpolation points" at 200,60
draw text "time = "&tt&" sec." at 240,100 ! refresh
end

' Montone cubic spline
' n   - # of x,y points
' x() - x coordinates of points
' y() - y coordinates of points
' ni  - number of intervals between two successive x coordinates
' xret() - x coordinates of resulting fit
' yret() - y coordinates of resulting fit
'
def m_c_spline(n,x(),y(),ni,xret(),yret())
dim sec(n-1),m(n)
for i=1 to n-1
  sec(i)=(y(i+1)-y(i))/(x(i+1)-x(i))
  if i=1 then ! m(1)=sec(1)
    else
    if sec(i-1)*sec(i)<0 then ! m(i)=0
      else ! m(i)=(sec(i-1)+sec(i))/2
      end if
    end if
  next i
for i=1 to n-1
  if sec(i)=0 then ! m(i)=0 ! m(i+1)=0 ! end if
  next i
m(n)=sec(n-1)
dt=1/ni
k=1 ! xret(1)=x(1) ! yret(1)=y(1) 
for i=1 to n-1
  a=y(i) ! b=m(i) ! xx=x(i) ! dx=x(i+1)-xx ! dy=y(i+1)-y(i) 
  c=3*dy-2*m(i)-m(i+1) ! d=m(i)+m(i+1)-2*dy
  for t=dt to 1 step dt
    k+=1 ! xx+=dx*dt ! xret(k)=xx
    yy=a+t*(b+t*(c+d*t)) ! yret(k)=yy
    next t
  next i
end def

def graph_it(n,x(),y(),ni,xi(),yi())
graphics ! graphics clear .7,.7,.7 ! refresh off
get screen size sw,sh ! xc=sw/2 ! yc=sh/2 ! gw=320
draw color 0,0,0 ! draw size 3
draw line xc-gw,yc to xc+gw,yc
draw line xc,yc-gw to xc,yc+gw
draw size 1 ! draw alpha .3 ! ds=gw/10
for i=-10 to 10
draw line xc-gw,yc-i*ds to xc+gw,yc-i*ds
draw line xc-i*ds,yc-gw to xc-i*ds,yc+gw
next i
fill color 1,0,0 ! draw alpha 1 ! gw=320 ! sc=gw/10
for i=1 to n
xx=xc+x(i)*sc ! yy=yc-y(i)*sc
fill circle xx,yy size 3
next i
draw color 0,0,1 ! draw size 1
draw to xc+xi(1)*sc,yc-yi(1)*sc
for i=2 to ni ! draw line to xc+xi(i)*sc,yc-yi(i)*sc ! next i
refresh
end def

[matt7]

Wow! This is very cool! I will definitely be checking this out, especially since your algorithm seems to be must faster than mine. I am calling my function every time I want to get an interpolated value along the curve, but I like the idea of just generating all the points over a specified range all at once. Also, the curve that is drawn looks great since there are no gaps where the slope is steep.

[Henko]

I think it is the best method for an application where the entire value range has to be processed sequenrially. I understand that this is the case in your gradient application.

If an application requires frequent processing of random values from a given range, i would again set up the full table first, then for each random value of x, use a binary search to find the segment where the sought x-value resides, and then do a simple lineair interpolation in that segment. I quess that this last step is hardly necessary if the table has a fine enough subdivision of the segments (the value of the variable "ni"), but then again, a linear interpolation in one segment is near to nothing.

By the way, when putting ni=1 in the code, a clear linear interpolation spline results (as should be expected).

The needed time grows linear with the value of "ni". On my iPad air2: T = 0.4 + 0.0284*ni msec.

[Henko]

I played somewhat with steep slopes and noticed a severe under-/overshoot. See the picture below.

IMG_1543.PNG (369.61 KiB) Viewed 2613 times

Reason enough to add the overshoot protection as given on the internet an used by Matt. The adapted code:

Code: Select all

option base 1
n=9 ! ni=20 ! nret=(n-1)*ni+1
dim x(n),y(n),xret(nret),yret(nret)

for i=1 to n ! read x(i),y(i) ! next i
data -7,-5, -5,9, -1,8, .5,6, 2,6, 3,10, 3.2,5.5, 6,5, 8,2
time reset
m_c_spline(n,x,y,ni,xret,yret)
tt=time()
graph_it(n,x,y,nret,xret,yret)
draw text "input  : "&n&" points" at 270,30
draw text "output : "&nret&" interpolation points" at 200,60
draw text "time = "&tt&" sec." at 240,100 ! refresh
end

' Monotone cubic spline
' n   - # of x,y input points
' x() - x coordinates of points
' y() - y coordinates of points
' ni  - # of intervals between two successive points
' xret() - x coordinates of resulting fit
' yret() - y coordinates of resulting fit
'
def m_c_spline(n,x(),y(),ni,xret(),yret())
dim sec(n-1),m(n)
for i=1 to n-1
  sec(i)=(y(i+1)-y(i))/(x(i+1)-x(i))
  if i=1 then ! m(1)=sec(1)
    else
    if sec(i-1)*sec(i)<0 then ! m(i)=0
      else ! m(i)=(sec(i-1)+sec(i))/2
      end if
    end if
  next i
for i=1 to n-1
  if sec(i)=0 then ! m(i)=0 ! m(i+1)=0
    else
    r=(m(i)/sec(i))^2+(m(i+1)/sec(i))^2
    if r>9 then ! t=3/sqrt(r) ! m(i)*=t ! m(i+1)*=t ! end if
  end if
  next i
m(n)=sec(n-1)
dt=1/ni
k=1 ! xret(1)=x(1) ! yret(1)=y(1) 
for i=1 to n-1
  xx=x(i) ! dx=x(i+1)-x(i) ! dxx=dx*dt ! dy=y(i+1)-y(i)
  a=y(i) ! b=m(i) ! c=3*dy-2*m(i)-m(i+1) ! d=m(i)+m(i+1)-2*dy
  for t=dt to 1 step dt
    k+=1 ! xx+=dxx ! xret(k)=xx ! yret(k)=a+t*(b+t*(c+d*t))
    next t
  next i
end def

def graph_it(n,x(),y(),ni,xi(),yi())
graphics ! graphics clear .7,.7,.7 ! refresh off
get screen size sw,sh ! xc=sw/2 ! yc=sh/2 ! gw=320 ! sc=gw/10
draw color 0,0,0 ! draw size 3
draw line xc-gw,yc to xc+gw,yc
draw line xc,yc-gw to xc,yc+gw
draw size 1 ! draw alpha .3 ! ds=gw/10
for i=-10 to 10
draw line xc-gw,yc-i*ds to xc+gw,yc-i*ds
draw line xc-i*ds,yc-gw to xc-i*ds,yc+gw
next i
fill color 1,0,0 ! draw alpha 1
for i=1 to n
xx=xc+x(i)*sc ! yy=yc-y(i)*sc
fill circle xx,yy size 3
next i
draw color 0,0,1 ! draw size 1
draw to xc+xi(1)*sc,yc-yi(1)*sc
for i=2 to ni ! draw line to xc+xi(i)*sc,yc-yi(i)*sc ! next i
refresh
end def

And the output, the undershoot beiing eliminated.

IMG_1544.PNG (367.1 KiB) Viewed 2613 times

[matt7]

EDIT (2/11/2019): Added a missing "dx*" to "b=" in the simplified line of code.

Henko,

Thanks again for sharing this routine with me for generating and returning all the interpolated values of a spline with one function call. It is way faster than the method I was using!

I was working through your code and convincing myself of the math you were doing because I wanted to understand it before making some tweaks for my own application, and I think I found two mistakes in the function m_c_spline.

1) For the coefficients a, b, c, and d calculated in m_c_spline, there are some terms that need to be multiplied by dx that are not. This affects the shape of your splines.

Here is what I believe the code should be (my corrections in red):

a=y(i) ! b=dx*m(i) ! c=3*dy-2*dx*m(i)-dx*m(i+1) ! d=dx*m(i)+dx*m(i+1)-2*dy

Which simplifies to:

Code: Select all

a=y(i) ! b=dx*m(i) ! c=3*dy-dx*(2*m(i)+m(i+1)) ! d=dx*(m(i)+m(i+1))-2*dy

I don't know if you were the one to work these out or what reference material was used, but it looks like this mistake comes from accidentally dropping the "h" in the below equation (from the Wikipedia page on monotonic cubic interpolation that I linked to in the OP) when solving for the coefficients:

where h = x_upper - x_lower (which is "dx" in your code)
(Note: I'm not referring to h00, h10, h01, or h11, but just the "h" that appears in front of the second and fourth terms.)

Compare before and after these corrections:


2) Some values of ni will result in the last point in xret and yret being 0,0. (The smallest value where this occurs is 52, so this is not really a big deal since ni will almost always be less than 52, but I still thought you should be aware.) This is because of the way you loop through t using STEP dt. In my code I am fixing this by replacing the line "for t=dt to 1 step dt" with:

Code: Select all

t=0
for j = 1 to ni
t+=dt

This increases processing time ever so slightly, but for me that's negligible. However, if you wanted the absolute fastest version of the function then just be aware of the possibility that the last value in xret and yret could be 0 for higher values of ni.

[Henko]

Hi Matt,

I am in the middle of developing a "neural network" program, which is a bit complicated. So, i will look at the possible error in my interpolation code after i finish the AI (study) tool. Anyhow thanks for bringing it to attention.
Henk