Interpolation after coordinate transformation

In this exercise, we will use a slice out of a 3-D CT-scan of a carbonate rock sample, shown in Figure 6a . Notice microfracture channels.


circle,rotate Figure 6. Slice of a CT-scan of a carbonate rock sample. (a) Original. (b) After clockwise rotation by $90^{\circ }$ .

The goal of the exercise is to apply a coordinate transformation to the original data. A particular transformation that we will study is coordinate rotation. Figure 6b shows the original slice rotated by 90 degrees. A 90-degree rotation in this case amounts to simple transpose. However, rotation by a different angle requires interpolation from the original grid to the modified grid.

The task of coordinate rotation is accomplished by the C program rotate.c. Two different methods are implemented: nearest-neighbor interpolation and bilinear interpolation.

To test the accuracy of different methods, we can rotate the original data in small increments and then compare the result of rotating to $360^{\circ}$ with the original data. Figure 7 compares the error of the nearest-neighbor and bilinear interpolations after rotating the original slice in increments of $20^{\circ}$ . The accuracy is comparatively low for small discontinuous features like microfracture channels.

To improve the accuracy further, we need to employ a longer filter. One popular choice is cubic convolution interpolation, invented by Robert Keys (a geophysicist, currently at ConocoPhillips). The cubic convolution filter can be expressed as the filter (Keys, 1981)


nearest,linear Figure 7. Error of different interpolation methods computed after full circle rotation in increments of 20 degrees. (a) Nearest-neighbor interpolation. (b) Bi-linear interpolation.

/* Rotate around. */
#include <rsf.h>

int main(int argc, char* argv[])
{
    int n1, n2, i1, i2, k1, k2;
    float x1, x2, c1, c2, cosa, sina, angle;
    float **orig, **rotd;
    char *interp;
    sf_file inp, out;

    /* initialize */
    sf_init(argc,argv);
    inp = sf_input("in");
    out = sf_output("out");

    /* get dimensions from input */
    if (!sf_histint(inp,"n1",&n1)) sf_error("No n1= in inp");
    if (!sf_histint(inp,"n2",&n2)) sf_error("No n2= in inp");

    /* get parameters from command line */
    if (!sf_getfloat("angle",&angle)) angle=90.;
    /* rotation angle */

    if (NULL == (interp = sf_getstring("interp"))) 
	interp="nearest";
    /* [n,l,c] interpolation type */

    /* convert degrees to radians */
    angle *= SF_PI/180.;
    cosa = cosf(angle);
    sina = sinf(angle);

    orig = sf_floatalloc2(n1,n2);
    rotd = sf_floatalloc2(n1,n2);

    /* read data */
    sf_floatread(orig[0],n1*n2,inp);

    /* central point */
    c1 = (n1-1)*0.5;
    c2 = (n2-1)*0.5;

    for (i2=0; i2 < n2; i2++) {
	for (i1=0; i1 < n1; i1++) {

	    /* rotated coordinates */
	    x1 = c1+(i1-c1)*cosa-(i2-c2)*sina;
	    x2 = c2+(i1-c1)*sina+(i2-c2)*cosa;

	    /* nearest neighbor */
	    k1 = floorf(x1); x1 -= k1;
	    k2 = floorf(x2); x2 -= k2;

	    switch(interp[0]) {
		case 'n': /* nearest neighbor */
		    if (x1 > 0.5) k1++;
		    if (x2 > 0.5) k2++;
		    if (k1 >=0 && k1 < n1 &&
			k2 >=0 && k2 < n2) {
			rotd[i2][i1] = orig[k2][k1];
		    } else {
			rotd[i2][i1] = 0.;
		    }
		    break;
		case 'l': /* bilinear */
		    if (k1 >=0 && k1 < n1-1 &&
			k2 >=0 && k2 < n2-1) {
			rotd[i2][i1] = 
			    (1.-x1)*(1.-x2)*orig[k2][k1]   +
			    x1     *(1.-x2)*orig[k2][k1+1] +
			    (1.-x1)*x2     *orig[k2+1][k1] +
			    x1     *x2     *orig[k2+1][k1+1];
		    } else {
			rotd[i2][i1] = 0.;
		    }
		    break;
		case 'c': /* cubic convolution */
		    /* !!! ADD CODE !!! */
		    break;
		default:
		    sf_error("Unknown interpolation %s",
			     interp);
		    break;
	    }
	}
    }

    /* write result */
    sf_floatwrite(rotd[0],n1*n2,out);

    exit(0);
}

from rsf.proj import *

# Download data
Fetch('slice.rsf','ctscan')
Flow('circle','slice','dd type=float')

grey = 'grey wanttitle=n screenratio=1 bias=128 clip=105'

Result('circle',grey)

# Rotate program
program = Program('rotate.c')
rotate = str(program[0])

# Rotate by 90 degrees
Flow('rotate',['circle',rotate],
     './${SOURCES[1]} angle=90 interp=nearest')

Result('rotate',grey)

# Mask for the circle
Flow('mask','circle',
     '''
     put d1=1 o1=-255.5 d2=1 o2=-255.5 |
     math output="sqrt(x1*x1+x2*x2)" |
     mask min=255.5 | dd type=float | 
     smooth rect1=3 rect2=3 |
     mask max=0 | dd type=float |
     put d1=0.1914 o1=0 d2=0.1914 o2=0
     ''')

for case in ('nearest','linear'): # !!! MODIFY ME !!!
    new = 'circle'
    rotates = []
    for r in range(18):
        old = new
        new = '%s-circle%d' % (case,r)
        Flow(new,[old,rotate],
             './${SOURCES[1]} angle=20 interp=%s' % case)
        Plot(new,grey)
        rotates.append(new)

    # Movie of rotating circle
    Plot(case,rotates,'Movie',view=1)

    # Plot error
    Result(case,[new,'circle','mask'],
           '''
           add scale=1,-1 ${SOURCES[1]} |
           add mode=p ${SOURCES[2]} |
           %s bias=0 scalebar=y clip=12
           wanttitle=y title="Error of %s Interpolation"
           ''' % (grey,case.capitalize()))

End()

$\displaystyle Z^{\sigma} \approx C(Z)$	$\textstyle =$	$\displaystyle -\frac{\sigma\,(1-\sigma)^2}{2}\,Z^{-1} + \frac{(1-\sigma)\,(2 + 2\,\sigma - 3 \sigma^2)}{2} +$
		$\displaystyle \frac{\sigma\,(1 + 4\,\sigma - 3\,\sigma^2)}{2}\,Z - \frac{(1-\sigma)\,\sigma^2}{2}\,Z^2\;.$	(8)