data.c

seismic software,very useful

/*
data object code
data object consists of axes, buffers, and dataset info
four input data formats recognized:
    FORMAT  CONDITIONS      ACTIONS
    vgrid   DataGridHeader()==1 DataLoadFloat() or DataLoadByte()
    segy    else segy=1     DataGetpar(), DataLoadFloat()
    floats  else float=1        DataGetpar(), DataLoadFloat()
    bytes   else            DataGetpar(), DataLoadByte()
*/
#include <stdio.h>
#include "main.h"
#include "axis.h"
#include "data.h"

/* initialize dataset from getpar */
Data DataInit()
{
    Data     data;
    Axis     axis;
    extern int infd;
    int      iaxis, i;
    FILE    *fd;

    NEW(Data, data, 1);
    /* fetch grid parameters from tail record or pars */
    data->vgrid = -1;
    GETPARINT("vgrid", "d", &data->vgrid);
    if (data->vgrid != 0) {
    if (DataGridHeader(data, infd) == 0) {
        DataGetpar(data);
    }
    } else {
    DataGetpar(data);
    }
    lseek(infd, 0, 0);
    strcpy(data->title, "stdin");
    if (GETPARSTRING("title", "s", data->title) == 0)
    GETPARSTRING("in", "s", data->title);
    GETPARSTRING("in", "s", data->file);
    /* get external script file */
    if (GETPARSTRING("script", "s", data->script)) {
    fd = fopen(data->script, "r");
    if (fd == NULL) {
        UIMessage("cant open script file");
    } else {
        for (i = 0; i < AxisSize(data->axis[DATA_AXIS3]) &&
         fgets(AxisScript(data->axis[DATA_AXIS3], i), sizeof(string),
               fd) != NULL; i++) {
        AxisScript(data->axis[DATA_AXIS3],
               i)[strlen(AxisScript(data->axis[DATA_AXIS3], i)) -
                  1] = '\0';
        }
    }
    fclose(fd);
    }
    data->tpow = 0.;
    GETPARFLOAT("tpow", "f", &data->tpow);
    data->gpow = 1.;
    GETPARFLOAT("gpow", "f", &data->gpow);
    data->transp = 0;
    GETPARINT("transp", "d", &data->transp);
    /* perform transpose by swapping axes */
    if (data->transp) {
    NEW(Axis, axis, 1);
    axis[0] = data->axis[DATA_AXIS1][0];
    data->axis[DATA_AXIS1][0] = data->axis[DATA_AXIS2][0];
    data->axis[DATA_AXIS2][0] = axis[0];
    FREE(axis);
    }
    data->value_base = DATA_VALUE_BASE;
    data->value_size = DATA_VALUE_SIZE;
    NEW(Buffer, data->buffer, data->size);
    return (data);
}

/* fetch grid parameters from getpar */
DataGetpar(data)
Data     data;
{
    int      iaxis;
    string   label;

    data->segy = DATA_SEGY;
    GETPARINT("segy", "d", &data->segy);
    data->esize = 0;
    GETPARINT("esize", "d", &data->esize);
    if (data->segy)
    data->esize = 4;
    else if (data->esize == 0)
    data->esize = 1;
    data->size = 1;
    if (data->segy)
    data->hbytes = 3600;
    else
    data->hbytes = 0;
    GETPARINT("hbytes", "d", &data->hbytes);
    for (iaxis = DATA_AXIS1; iaxis <= DATA_AXIS3; iaxis++) {
    /* accumulating size are axis strides */
    data->axis[iaxis] = AxisInit(iaxis, data->size);
    data->size *= AxisSize(data->axis[iaxis]);
    }
    data->axis[DATA_AXIS4] =
    AxisInit2(DATA_AXIS4, data->size, "n4", 1, 0., 1., 1.);
    data->size *= AxisSize(data->axis[DATA_AXIS4]);
    data->axis[DATA_AXIS5] =
    AxisInit2(DATA_AXIS5, data->size, "n5", 1, 0., 1., 1.);
    data->size *= AxisSize(data->axis[DATA_AXIS5]);
    data->max = DATA_HIGH;
    GETPARFLOAT("max", "f", &data->max);
    data->high = data->max;
    GETPARFLOAT("high", "f", &data->high);
    GETPARFLOAT("pclip", "f", &data->high);
    GETPARFLOAT("clip", "f", &data->high);
    if (data->high != DATA_HIGH) {
    data->min = -data->high;
    data->low = data->min;
    } else {
    data->min = DATA_LOW;
    data->low = data->min;
    }
    GETPARFLOAT("low", "f", &data->low);
    GETPARFLOAT("nclip", "f", &data->low);
    GETPARFLOAT("min", "f", &data->min);
    data->cent = DATA_CENT;
    GETPARFLOAT("perc", "f", &data->cent);
    GETPARFLOAT("cent", "f", &data->cent);
    strcpy(label, "samples");
    GETPARSTRING("value", "s", label);
    data->axis[DATA_VALUE] = AxisInit2(DATA_VALUE, 1, label, DATA_VALUE_SIZE,
                       data->low,
                       (data->high -
                    data->low) / (DATA_VALUE_SIZE - 1), 1);
}

/* fetch grid parameters from tail of vgrid format */
DataGridHeader(data, fd)
Data     data;
int      fd;
{
    int      size, n1, n2, n3, n4, n5, size1;
    string   label;

    size = lseek(fd, -sizeof(data->gh), 2);
    if (read(fd, &data->gh, sizeof(data->gh)) < sizeof(data->gh)) {
    lseek(fd, 0, 0);
    data->gh.dtype = 0.;
    return (0);
    }
    if (data->gh.scale == 0.0) {
    data->gh.dtype = 0.;
    return (0);
    }
    size1 = rint(data->gh.n1 / data->gh.scale)
    * rint(data->gh.n2 / data->gh.scale)
    * rint(data->gh.n3 / data->gh.scale)
    * rint(data->gh.n4 / data->gh.scale)
    * rint(data->gh.n5 / data->gh.scale)
    * rint(data->gh.dtype / data->gh.scale);
    if (size != size1) {
    lseek(fd, 0, 0);
    data->gh.dtype = 0.;
    return (0);
    }
    data->vgrid = 1;
    strcpy(label, "time");
    GETPARSTRING("label1", "s", label);
    n1 = rint(data->gh.n1 / data->gh.scale);
    data->axis[DATA_AXIS1] =
    AxisInit2(DATA_AXIS1, 1, label, n1, data->gh.o1 / data->gh.scale,
          data->gh.d1 / data->gh.scale, 1.);
    strcpy(label, "cdp");
    GETPARSTRING("label2", "s", label);
    n2 = rint(data->gh.n2 / data->gh.scale);
    data->axis[DATA_AXIS2] = AxisInit2(DATA_AXIS2, n1, label,
                       n2, (data->gh.o2 / data->gh.scale),
                       (data->gh.d2 / data->gh.scale), 1.);
    strcpy(label, "line");
    GETPARSTRING("label3", "s", label);
    n3 = rint(data->gh.n3 / data->gh.scale);
    GETPARINT("n3", "d", &n3);
    data->axis[DATA_AXIS3] = AxisInit2(DATA_AXIS3, n1 * n2, label,
                       n3, (data->gh.o3 / data->gh.scale),
                       (data->gh.d3 / data->gh.scale), 1.);
    strcpy(label, "n4");
    GETPARSTRING("label4", "s", label);
    n4 = rint(data->gh.n4 / data->gh.scale);
    GETPARINT("n4", "d", &n4);
    data->axis[DATA_AXIS4] = AxisInit2(DATA_AXIS4, n1 * n2 * n3, label,
                       n4, (data->gh.o4 / data->gh.scale),
                       (data->gh.d4 / data->gh.scale), 1.);
    strcpy(label, "n5");
    GETPARSTRING("label5", "s", label);
    n5 = rint(data->gh.n5 / data->gh.scale);
    GETPARINT("n5", "d", &n5);
    data->axis[DATA_AXIS5] = AxisInit2(DATA_AXIS5, n1 * n2 * n3 * n4, label,
                       n5, (data->gh.o5 / data->gh.scale),
                       (data->gh.d5 / data->gh.scale), 1.);
    strcpy(label, "amplitude");
    GETPARSTRING("amplitude", "s", label);
    data->low = data->gh.gmin / data->gh.scale;
    data->high = data->gh.gmax / data->gh.scale;
    GETPARFLOAT("high", "f", &data->high);
    GETPARFLOAT("pclip", "f", &data->high);
    GETPARFLOAT("clip", "f", &data->high);
    GETPARFLOAT("low", "f", &data->low);
    GETPARFLOAT("nclip", "f", &data->low);
    data->axis[DATA_VALUE] = AxisInit2(DATA_VALUE, 1, label,
                       DATA_VALUE_SIZE, data->low,
                       (data->high -
                    data->low) / (DATA_VALUE_SIZE - 1), 1.);
    data->esize = rint(data->gh.dtype / data->gh.scale);
    data->size = n1 * n2 * n3 * n4 * n5;
    data->max = data->high;
    data->min = data->low;
    return (1);
}

/* decide which data reading routine to use */
DataLoad()
{
    extern Data data;

    switch (data->esize) {
       case 4:
    DataLoadFloat(data);
    break;
       case 1:
    DataLoadByte(data);
    break;
    }
    UIMessage("data loaded");
}

/* load byte data; compress to 7 bits- 0-128 */
DataLoadByte(data)
Data     data;
{
    extern int infd;
    byte     table;
    int      i, value;
    register byte tp;
    register Buffer bp, be;

    /* read data */
    UIMessage("loading byte data ...");
    if (data->hbytes)
    read(infd, data->buffer, data->hbytes);
    read(infd, data->buffer, data->size);

    /* compute color compression table table */
    NEW(byte, table, 256);
    /* bounds */
    table[0] = 0;
    /* compute lookup table */
    for (i = 0; i < 256; i++) {
    table[i] = i / 2;
    }
    table[1] = 1;
    for (bp = data->buffer, be = bp + data->size, tp = table; bp < be; bp++) {
    *bp = tp[*bp];
    }
    DataComputeHistogram(data);
    FREE(table);
}


/* load and convert float or segy data */
DataLoadFloat(data)
Data     data;
{
    int      i1, i2, i3, n1, n2, n3, base, size, head;
    float   *buf1, *buf2, DataCent(), *trace, *tgain;
    register float scale, bias, *bp1, *bp2, *fp, *fe, *gp;
    register int datum;
    register byte dp;
    double   pow(), log();
    extern int infd;

    if (data->segy)
    UIMessage("loading segy data ...");
    else
    UIMessage("loading float data ...");
    /* calculate gain parameters: clip, gpow */
    n1 = AxisSize(data->axis[DATA_AXIS1]);
    n2 = AxisSize(data->axis[DATA_AXIS2]);
    n3 = AxisSize(data->axis[DATA_AXIS3])
    * AxisSize(data->axis[DATA_AXIS4])
    * AxisSize(data->axis[DATA_AXIS5]);
    if (data->segy)
    head = 60;
    else
    head = 0;
    NEW(float *, buf1, n1 * n2);
    NEW(float *, buf2, n1 * n2);
    NEW(float *, tgain, n1);

    size = (n1 + head) > DATA_HEAD1 ? (n1 + head) : DATA_HEAD1;
    size = size > data->hbytes ? size : data->hbytes;
    NEW(float *, trace, size);

    fe = trace + n1 + head;
    bp1 = buf1;
    bp2 = buf2;
    /* time gain vector */
    for (i1 = 0; i1 < n1 - 1; i1++) {
    tgain[i1] = AxisValue(data->axis[DATA_AXIS1], i1 + 1);
    tgain[i1] = tgain[i1] != 0.0 ? pow(tgain[i1], data->tpow) : 1.0;
    }
    tgain[n1 - 1] = tgain[n1 - 2];
    /* remember first panel in buf1; copy for percentile */
    if (data->segy) {
    read(infd, trace, DATA_HEAD0 * sizeof(trace[0]));
    read(infd, trace, DATA_HEAD1 * sizeof(trace[0]));
    } else if (data->hbytes) {
    read(infd, trace, data->hbytes);
    }
    for (i2 = 0; i2 < n2; i2++) {
    read(infd, trace, (n1 + head) * sizeof(trace[0]));
    for (fp = trace + head, fe = fp + n1, gp = tgain; fp < fe;
         fp++, bp1++, bp2++, gp++) {
        *bp1 = *fp * *gp;
        *bp2 = *bp1;
    }
    }
    if (data->low == DATA_LOW && data->high == DATA_HIGH) {
    string   label;

    data->min = DataCent(buf2, n1 * n2, DATA_CENT_MIN);
    data->low = DataCent(buf2, n1 * n2, DATA_CENT_LOW);
    data->high = DataCent(buf2, n1 * n2, DATA_CENT_HIGH);
    data->max = DataCent(buf2, n1 * n2, DATA_CENT_MAX);
    strcpy(label, AxisLabel(data->axis[DATA_VALUE]));
    FREE(data->axis[DATA_VALUE]);
    data->axis[DATA_VALUE] =
        AxisInit2(DATA_VALUE, 1, label, DATA_VALUE_SIZE + 1, data->low,
              (data->high - data->low) / DATA_VALUE_SIZE, 1);
    }
    FREE(buf2);
    if (data->low == 0.0 && data->high == 0.) {
    err("first panel appears to be all zeros; set clip");
    }
    if (data->high > data->low) {
    bias = -data->low;
    scale = DATA_VALUE_SIZE / (data->high - data->low);
    } else {
    bias = 0.;
    scale = 0.5 * DATA_VALUE_SIZE / data->high;
    }
    dp = data->buffer;
    for (fp = buf1, fe = buf1 + n1 * n2; fp < fe; fp++, dp++) {
    datum = (*fp + bias) * scale;
    datum = datum > 0 ? datum : 0;
    datum = datum < DATA_VALUE_SIZE ? datum : DATA_VALUE_SIZE - 1;
    datum += DATA_VALUE_BASE;
    *dp = datum;
    }
    FREE(buf1);
    /* convert rest of traces */
    for (i3 = 1; i3 < n3; i3++)
    for (i2 = 0; i2 < n2; i2++) {
        if (read(infd, trace, (n1 + head) * sizeof(trace[0])) < 0)
        break;
        for (fp = trace + head, fe = fp + n1, gp = tgain; fp < fe;
         fp++, dp++, gp++) {
        datum = (*fp * *gp + bias) * scale;
        datum = datum > 0 ? datum : 0;
        datum = datum < DATA_VALUE_SIZE ? datum : DATA_VALUE_SIZE - 1;
        datum += DATA_VALUE_BASE;
        *dp = datum;
        }
    }
    DataComputeHistogram(data);
    FREE(trace);
    FREE(tgain);
}

/*      percentile subroutine
 *      based on Canales, SEP-10
 *      p - percentile <0.,99.999999999999>
 *      x - data
 *      n - vector length
 *  this routine changes data order, so sort a copy
 */
float    DataCent(x, n, p)
float   *x, p;
int      n;
{
    int      q;
    register float *i, *j, ak;
    float   *low, *hi, buf, *k;

    p = p < 99.999 ? p : 99.999;
    p = p > 0.0 ? p : 0.0;
    q = (p * n) / 100.;
    if (q == n)
    q--;
    for (low = x, hi = x + n - 1, k = x + q; low < hi;) {
    ak = *k;
    i = low;
    j = hi;
    do {
        while (*i < ak)
        i++;
        while (*j > ak)
        j--;
        if (i <= j) {
        buf = *i;
        *i++ = *j;
        *j-- = buf;
        }
    } while (i <= j);
    if (j < k)
        low = i;
    if (k < i)
        hi = j;
    }
    return (*k);
}