INTERNAL: Explanation of Cubic Keyword to INTERPOLATE

Topic:
Subject: Explanation of Cubic

Cubic convolution is an interpolation method that closely approximates the theoretically optimum sinc interpolation function using cubic polynomials. According to sampling theory, details of which are beyond the scope of this document: if the original signal, f, is a band-limited signal, with no frequency component larger than Omega_0, and f is sampled with spacing less than or equal to 1/(2*Omega_0), then f can be reconstructed by convolving with a sinc function. Sinc(x) = sin(pi * x) / (pi * x).

In the one dimensional case, four neighboring points are used, while in the two dimensional case 16 points are used. Cubic convolution interpolation is significantly slower than bilinear interpolation.

For further details see: Rifman, S.S. and D.M. McKinnon, "Evaluation
of Digital Correction Techniques for ERTS Images --- Final Report",
Report 20634-6003-TU-00, TRW Systems, Redondo Beach, Calif, July 1974.Discussion:
To Ray Sterner From David Stern

I saw your posting mentioning that INTERPOLATE with /CUBIC has a bug. Well, it doesn't as far as I know. The people that complain that it sometimes returns negative results given only positive inputs, or that it has wiggles, don't understand the desired result.

When you sample a continuous signal with a rectangular window you smear the high frequencies. "Cubic convolution" is an approximation to the ideal sampling filter, sinc(x) = sin(pi*x)/(pi*x), and is meant to restore these high frequencies. Any filter that enhances high frequencies, must have negative components, hence the possibility of a negative output given only positive inputs.

My references are:

An Introduction to Digital Image Processing, Niblak, Prentice/Hall
International, 1985, Pages 142-149.

Digital Image Warping, Wolberg, IEEE Computer Society Press, 1990.
Page 129.

It should not be used on artificial data, digital graphics, improperly sampled data, data with Poison statistics, or contrived examples. Your DEM datasets are probably not properly sampled, that's why they seem to look best with bilinear interpolation. For example, if the smallest variation that you are concerned with can occur in X meters, you must sample the elevation every X/2 meters.

I hope that this clarifies things.

From Ray Sterner to the Newsgroup
In an earlier post I gave a warning about using the CUBIC keyword in INTERPOLATE, CONGRID, ROT:

Warning: INTERPOLATE has a keyword to do cubic interpolation. It's easy to use, just add /cubic to the call. However, don't assume this is better. Check it very carefully before you rely on it. A number of people have reported problems with /cubic in other routines (CONGRID, ROT). If the result of the cubic interp. is carefully examined it will be found to have high frequency wiggles in the values (with a period of roughly 20 some pixels for the case I just checked). As of V3.6 beta this problem still exists (for INTERPOLATE anyway). The default, bilinear, works very well.>/UL>I have been advised that this option is actually working correctly. I now believe that to be true. However I suspect a number of people misunderstood how this option actually works. Here is some code to demonstrate:Solution:


   ; Some random data.
   x=indgen(4) 
   y=[1.,2,5,10]
 
   ; Interpolate between points.
   x2=findgen(100)/100.*3. 
 
   ; Linear (L).
   yl=interpolate(y,x2) 
 
   ; Nearest Neighbor (NN).
   yn=interpolate(y,round(x2)) 
 
   ; Cubic convolution (CC).
   yc=interpolate(y,x2,/cubic) 
 
   ; Original data (4 points).
   plot,x,y 
 
   ; NN
   oplot,x2,yn 
   ; L
   oplot,x2,yl,psym=-4 
   ; CC
   oplot,x2,yc 
 
;The four original points are connected by 3 straight lines.
;Nearest Neighbor interpolation, yn, is the worst case.
;Linear interpolation, yl, is better.
;Cubic interpolation, yc, might be expected to give a nice 
; smooth curve through the points, 
;   but that is a misunderstanding. The cubic
;interpolation is working correctly but it is not a spline fit.
 
;The following code plots the expected curve.
 
   ; Set up spline.
   r=nr_spline(x,y) 
   ; Do spline interpolation.
   ycs=nr_splint(x,y,r,x2) 
   oplot,x2,ycs
 
;Unfortunately this doesn't help for 2-d interpolation.
 
;Besides making processed images look better, the cubic convolution
;option does have other useful applications. Here is some code that
;shows how high spatial frequencies are restored better using the
;/cubic keyword:
 
   ; Generate a curve.
   x = findgen(1000) 
   y = sin(x/(1+x/200.))
   ; Plot entire curve.
   plot,x,y 
 
   ;Just look at the high frequency part of the data:
   plot,x,y,xran=[0,50],yran=[-1.2,1.2]
 
;Note that it is somewhat ragged due to undersampling. 
;Now interpolate up to a larger number of points using /cubic:
 
   x2 = findgen(20000)/20.
   y2 = interpolate(y,x2,/cubic)
   oplot,x2,y2
 
;Note that the interpolated curve appears better and might also be 
;more useful for locating peaks. The peaks are not exactly where 
;they should be but are closer than the uninterpolated peaks. Also 
;note that the interpolated curve overshoots the true values as can 
;be better seen by adding horizontal lines:
 
   plots,[0,1000],[1,1]
   plots,[0,1000],-[1,1]