Homework 3 |
In this exercise, we will use a slice out of a 3-D CT-scan of a carbonate rock sample, shown in Figure 1a. Notice microfracture channels.
circle,rotate
Figure 1. 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 1b 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 (alternatively, Python script rotate.py.) 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 2 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)
nearest,linear
Figure 2. 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 method %s",interp); break; } } } /* write result */ sf_floatwrite(rotd[0],n1*n2,out); exit(0); } |
#!/usr/bin/env python import sys import math import numpy as np import m8r # initialize par = m8r.Par() inp = m8r.Input() out = m8r.Output() # get dimensions from input n1 = inp.int('n1') n2 = inp.int('n2') # get parameters from command line angle = par.float('angle',90.) # rotation angle interp = par.string('interp','nearest') # [n,l,c] interpolation type # convert degrees to radians angle = angle*math.pi/180. cosa = math.cos(angle) sina = math.sin(angle) orig = np.zeros([n1,n2],'f') rotd = np.zeros([n1,n2],'f') # read data inp.read(orig) # central point c1 = (n1-1)*0.5 c2 = (n2-1)*0.5 for i2 in range(n2): for i1 in range(n1): # rotated coordinates x1 = c1+(i1-c1)*cosa-(i2-c2)*sina x2 = c2+(i1-c1)*sina+(i2-c2)*cosa # nearest neighbor k1 = int(math.floor(x1)) k2 = int(math.floor(x2)) x1 -= k1 x2 -= k2 if interp[0] == 'n': # nearest neighbor if x1 > 0.5: k1 += 1 if x2 > 0.5: k2 += 1 if k1 >=0 and k1 < n1 and k2 >=0 and k2 < n2: rotd[i2,i1] = orig[k2,k1] else: rotd[i2,i1] = 0. elif interp[0] == 'l': # bilinear if k1 >=0 and k1 < n1-1 and k2 >=0 and 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. elif interp[0] == 'c': # cubic convolution # !!! ADD CODE !!! break else: sys.stderr.write('Unknown method "%s"' % interp) sys.exit(1) # write result */ out.write(rotd) sys.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 # exe = Program('rotate.c') # UNCOMMENT ABOVE AND COMMENT BELOW IF YOU WANT TO USE C exe = Command('rotate.exe','rotate.py','cp $SOURCE $TARGET') AddPostAction(exe,Chmod(exe,0o755)) rotate = str(exe[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() |
Your task:
scons viewto reproduce the figures on your screen.
scons nearest.vpland
scons linear.vplto see movies of incremental slice rotation with different methods.
Homework 3 |