Hagedoorn Construction

I’ve been writing GPU code for Kirchhoff Time Migration, a technique applied in seismic imaging. Among the tools used by geophysicists and other workers in oil + gas to perform seismic imaging, the most ubiquitous is probably KTM. This method is directly preceeded by the seismic image construction technique developed primarily by Hagedoorn, although both approaches essentially represent Greens functions solutions to the scalar wave equation. Data surveys are constructed using an array of pulsed sources and corresponding receivers. Data recorded at receiver locations (referred to as ‘traces’) generally correspond to a particular source-receiver pair, and in the simplest case source and receiver are coincident. An image into the earth is constructed by first calculating the traveltime for a given src/rec position and image point, using these values as well as the one way traveltime and velocity at the depth in question. This traveltime is then used as an index into the trace data, to find the amplitude point (energy) which contributes to the image point. Many points from different traces constructively and destructively interfere to create the subsequent image. The initial approach developed by Hagedoorn used little more than a compass and ruler to perform this imaging task, which is to quote the vernacular ‘old school’ and provides great pedagogical insight into KTM. There is a fantastic article by Bleistein here, and below is an example of a Hagedoorn construction performed in octave; the surface is the superposition of the arcs, like figure 2 of the article.

---

black angels

my new favorite band, great psychedelic rock

UNIX redirect,pipes & reconfigurable data processing

In various sciences one is often invested in processing voluminous data which follows a given format, with many examples following (hopefully) the same format. A typical data file might have a line of environment variables, followed by raw data, which repeats in columns and/or rows eg.,

4,23,45,56,78
12,23,45,56,67,78,89,87,…,72
32,34,23,21,23,56,43,23,…,34
34,54,32,45,89,76,54,98,…,58
67,67,88,32,34,21,22,97,…,51

A typical section of C code which would use this might look like:

float config[8];
scanf(“%f,%f,%f,%f,%f”,&config[1],&config[2],&config[3],
&config[4],&config[5]);

printf(“%f,%f,%f,%f,%f\n”,config[1],config[2],config[3],
config[4],config[5]);

int RECORDS = (int) config[1];
int i;
if ( RECORDS < 100){
int SIZE = sizeof(float)*RECORDS;

float *deltaCSMean = malloc(SIZE);
float *deltaCSSigma = malloc(SIZE);
float *CQMean = malloc(SIZE);
float *CQSigma = malloc(SIZE);
float *etaMean = malloc(SIZE);
float *etaSigma = malloc(SIZE);
float *brdF2 = malloc(SIZE);
float *brdF1 = malloc(SIZE);
float *amplitude = malloc(SIZE);

for (i=0; i< RECORDS; i++){
scanf(“%f,%f,%f,%f,%f,%f,%f,%f,%f”,&deltaCSMean[i],
&deltaCSSigma[i],&CQMean[i],&CQSigma[i],
&etaMean[i],&etaSigma[i],
&brdF2[i],&brdF1[i],&amplitude[i]);
printf(“%f,%f,%f,%f,%f,%f,%f,%f,%f\n”,deltaCSMean[i],
deltaCSSigma[i],CQMean[i],CQSigma[i],
etaMean[i],etaSigma[i],brdF2[i],brdF1[i],amplitude[i]);
}
}

Using redirect a simple invocation of executable foo using input foo1.txt and output foo2.txt would of course be:

./foo < foo1.txt > foo2.txt

which is fine until one needs different data from the file. Supposing in this example data repeats along columns, then one can use awk and a pipe to reshape the data input, removing the need to edit and recompile the C source for different data subsets:

gawk ‘BEGIN{FS=”,”; print “4,12,23,44,56″; getline;} {for(i=lowBnd;i<upBnd;i++){printf(“%f,”,$i) }; printf(“%f\n”,$upBnd) }’ foo1.txt | ./foo > foo2.txt

where lowBnd and upBnd correspond to the data column limits in the input file. Note that the environment variables are written during the BEGIN block, allowing the specification of (for instance) a different job for the new data set.