10921 Rate this article:

Upsampling images using a Lagrange polynomial method

A few weeks ago I posted about using the Lanczos kernel for resampling images to a higher resolution. This week I am continuing with the same example, but adding in the Lagrange resampling method. Both Lagrange and Lanczos have some similar characteristics in that they show better detail than a purely linear interpolation. Both methods can also be adapted to an irregularly gridded dataset instead of the raster images used in my examples here. The code produces 4 upsampled images using different methods, and the results are shown below.
function lanczos, data
  xval = [-3:3:.25]
  lanc3 = 3*sin(!pi*xval)*(sin(!pi*xval/3d)/!pi/!pi/xval/xval)
  lanc3[where(xval eq 0)] = 1
  l2d = lanc3 # lanc3
  ; high resolution version
  msk = fltarr(data.dim*4)
  msk[0:*:4,0:*:4] = data
  hi = convol(msk, l2d, /edge_trunc)
  hi = byte(round(hi>0<255))
  return, hi
function lagrange, a, x, y
  compile_opt idl2, logical_predicate
  xf = floor(x)
  yf = floor(y)
  x1 = x - xf
  y1 = y - yf
  off = [-1,0,1,2]
  retval = replicate(0., size(x, /DIM))
  weightx = replicate(1., [size(x1, /DIM),4])
  weighty = replicate(1., [size(x1, /DIM),4])
  for i=0,3 do begin
    for j=0,3 do begin
      if i ne j then begin
        weightx[*,*,i] *= (x1 - off[j]) / (off[i] - off[j])
        weighty[*,*,i] *= (y1 - off[j]) / (off[i] - off[j])
  for i=0,3 do begin
    for j=0,3 do begin
      retval += weightx[*,*,j] * weighty[*,*,i] * a[xf+off[j], yf+off[i]]
  return, retval
pro upsample_example
  compile_opt idl2,logical_predicate
  ; Read the original image data
  f = filepath('moon_landing.png', subdir=['examples','data'])
  data = read_png(f)
  dim = data.dim
  window, xsize=dim[0], ysize=dim[1], 0, title='Original full size'
  tv, data
  ; Define a zoomed in are on the flag.
  xs = 120
  ys = 105
  dx = 60
  dy = 100
  ; display upsampled 4x with pixel replication
  window, xsize=4*dx, ysize=4*dy, 1, title='CONGRID pixel-replication'
  tv, congrid(data[xs:xs+dx-1,ys:ys+dy-1],4*dx,4*dy)
  ; display upsampled 4x with bilinear interpretation
  window, xsize=4*dx, ysize=4*dy, 2, title='CONGRID linear'
  tv, congrid(data[xs:xs+dx-1,ys:ys+dy-1],4*dx,4*dy,/interp)
  ; display upsampled 4x with Lanczos convolution
  window, xsize=4*dx, ysize=4*dy, 3, title='Lanczos'
  tv, (lanczos(data))[xs*4:xs*4+dx*4-1,ys*4:ys*4+dy*4-1]
  ; Lagrange
  window, xsize=4*dx, ysize=4*dy, 4, title='Lagrange'
  xcoord = [float(xs):xs+dx:0.25]
  ycoord = [float(ys):ys+dy:0.25]
  tv, byte(lagrange(float(data), $
    xcoord # replicate(1,1,ycoord.length), $
    replicate(1,xcoord.length) # ycoord)>0<255)
 Pixel replication


Bi-linear interpolation

Lanczos resampling

Lagrange resampling