Homework 3

Next: Your own data Up: Homework 3 Previous: Fourier compression

# 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 .

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 with the original data. Figure 7 compares the error of the nearest-neighbor and bilinear interpolations after rotating the original slice in increments of . 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)

 (8)

and is designed to approximate the ideal sinc-function interpolator.

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 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() 

1. Change directory to hw4/rotate
2. Run
scons view

to reproduce the figures on your screen.
3. Additionally, you can run
scons nearest.vpl

and
scons linear.vpl

to see movies of incremental slice rotation with different methods.
4. Modify the rotate.c program and the SConstruct file to implement the cubic convolution interpolation and to compare its results with the two other methods.
5. EXTRA CREDIT for implementing an interpolation algorithm, which is more accurate than cubic convolution.

 Homework 3

Next: Your own data Up: Homework 3 Previous: Fourier compression

2014-10-02