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 <math.h>

#include "par.h"

#include "main.h"
#include "axis.h"
#include "data.h"
#include "color.h"
#include "draw.h"
#include "ui_menu.h"
#include "ui_window.h"

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

    cwp_String parString;

    {
        extern int _alloc;

        data = (Data) malloc((size_t) sizeof(data[0]));
        _alloc += sizeof(data[0]);
        if( data == 0 ){
            err("cant allocate %d bytes for  data; %d already allocated",
                sizeof(data[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " data", sizeof(data[0]));
        }
    };
    /* fetch grid parameters from tail record or pars */
    data->vgrid = -1;
    data->overlay_mode  = 0;

    getparint("vgrid", &data->vgrid);
    if( data->vgrid != 0 ){
        if( DataGridHeader(data, infd) == 0 ){
            DataGetpar(data);
        }
    } else {
        DataGetpar(data);
    }
    lseek64(infd, 0, 0);
    strcpy(data->title, "stdin");
    if( getparstring("title", &parString) == 0 ){
        if( getparstring("in", &parString ) ){
            strcpy( data->title ,parString );
        }
    }else{
        strcpy( data->title ,parString );
    }

    if( getparstring("in", &parString ) ){
        strcpy( data->file ,parString );
    }
 
    /* get external script file */

    if( getparstring("script", &parString) ){
        strcpy( data->script ,parString );
        fd = fopen(data->script, "r");
        if( fd == NULL ){
            UIMessage("cant open script file");
        } else {
            int n = AxisSize(data->axis[DATA_AXIS3]);
            for( i = 0; i < n &&
                 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", &data->tpow);
    data->gpow = 1.;
    getparfloat("gpow", &data->gpow);
    data->transp = 0;
    getparint("transp", &data->transp);
    /* perform transpose by swapping axes */
    if( data->transp ){

        {
            extern int _alloc;

            axis = (Axis) malloc((size_t)sizeof(axis[0]));
            _alloc += sizeof(axis[0]);
            if( axis == 0 ){
                err("cant allocate %d bytes for  axis; %d already allocated",
                    sizeof(axis[0]), _alloc);
            }
            if( memwatch ){
                (void) printf("malloc %s=%d\n", " axis", sizeof(axis[0]));
            }
        };
        axis[0] = data->axis[DATA_AXIS1][0];
        data->axis[DATA_AXIS1][0] = data->axis[DATA_AXIS2][0];
        data->axis[DATA_AXIS2][0] = axis[0];
        if( axis ){
            free(axis);
            axis = 0;
            if( memwatch ){
                printf("free %s\n", "axis");
            }
        };
    }
    data->value_base = DATA_VALUE_BASE;
    data->value_size = DATA_VALUE_SIZE;
    {
        extern int _alloc;

        data->buffer = (Buffer) malloc((size_t) (data->size) * sizeof(data->buffer[0]));
        _alloc += (data->size) * sizeof(data->buffer[0]);
        if( data->buffer == 0 ){
            err
                ("cant allocate %d bytes for  data->buffer; %d already allocated",
                 (data->size) * sizeof(data->buffer[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " data->buffer",
                          (data->size) * sizeof(data->buffer[0]));
        }
    };
    return (data);
}

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

    data->segy = DATA_SEGY;
    getparint("segy", &data->segy);
    data->esize = 0;
    getparint("esize", &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", &data->hbytes);
    for( iaxis = DATA_AXIS1; iaxis <= DATA_AXIS3; iaxis++ ){

        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", &data->max);
    data->high = data->max;
    getparfloat("high", &data->high);
    getparfloat("pclip", &data->high);
    getparfloat("clip", &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", &data->low);
    getparfloat("nclip", &data->low);
    getparfloat("min", &data->min);
    data->cent = DATA_CENT;
    getparfloat("perc", &data->cent);
    getparfloat("cent", &data->cent);
    strcpy(label, "samples");
    if( getparstring("value", &parString) ){
       strcpy(label ,parString );
    }
    data->axis[DATA_VALUE] = AxisInit2(DATA_VALUE ,1 ,label 
                            ,DATA_VALUE_SIZE ,data->low 
                            ,(data->high - data->low) / (DATA_VALUE_SIZE-1) 
                            ,1.0);
}

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

    size = lseek64(fd, -(long long)sizeof(data->gh), 2);
    if( read(fd, &data->gh, sizeof(data->gh)) < sizeof(data->gh) ){
        lseek64(fd, 0, 0);
        data->gh.dtype = 0.;
        return (0);
    }
    if( data->gh.scale == 0.0 ){
        data->gh.dtype = 0.;
        return (0);
    }
    size1 = irint(data->gh.n1 / data->gh.scale)
        * irint(data->gh.n2 / data->gh.scale)
        * irint(data->gh.n3 / data->gh.scale)
        * irint(data->gh.n4 / data->gh.scale)
        * irint(data->gh.n5 / data->gh.scale)
        * irint(data->gh.dtype / data->gh.scale);
    if( size != size1 ){
        lseek64(fd, 0, 0);
        data->gh.dtype = 0.;
        return (0);
    }
    data->vgrid = 1;
    if( irint( data->gh.gtype /data->gh.scale ) == 5 ){
       data->overlay_mode = 1;
    }
    strcpy(label, "time");

    if( getparstring("label1", &parString ) ){
        strcpy( label ,parString );
    }
    n1 = irint(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");
    if( getparstring("label2", &parString ) ){
        strcpy( label ,parString );
    }
    n2 = irint(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");
    if( getparstring("label3", &parString ) ){
        strcpy( label ,parString );
    }
    n3 = irint(data->gh.n3 / data->gh.scale);
    getparint("n3", &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");
    if( getparstring("label4", &parString ) ){
        strcpy( label ,parString );
    }
    n4 = irint(data->gh.n4 / data->gh.scale);
    getparint("n4", &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");
    if( getparstring("label5", &parString ) ){
        strcpy( label ,parString );
    }
    n5 = irint(data->gh.n5 / data->gh.scale);
    getparint("n5", &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");
    if( getparstring("amplitude", &parString ) ){
        strcpy( label ,parString );
    }
    data->low = data->gh.gmin / data->gh.scale;
    data->high = data->gh.gmax / data->gh.scale;
    getparfloat("high", &data->high);
    getparfloat("pclip", &data->high);
    getparfloat("clip", &data->high);
    getparfloat("low", &data->low);
    getparfloat("nclip", &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 = irint(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 */
void DataLoad(void)
{
    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 */
void DataLoadByte(Data data)
{
    extern int infd;
    byte     table;
    int      i;
    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 */

    {
        extern int _alloc;

        table = (byte) malloc((size_t) (256) * sizeof(table[0]));
        _alloc += (256) * sizeof(table[0]);
        if( table == 0 ){
            err("cant allocate %d bytes for  table; %d already allocated",
                (256) * sizeof(table[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " table", (256) * sizeof(table[0]));
        }
    };

    table[0] = 0;

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

    if( table ){
        free(table);
        table = 0;
        if( memwatch ){
            printf("free %s\n", "table");
        }
    };
}

void DataLoadFloat(Data data)
{
    int      i1, i2, i3, n1, n2, n3, size, head;
    float   *buf1, *buf2, DataCent(float *x, int n, float p), *trace, *tgain;
    float scale, bias, *bp1, *bp2, *fp, *fe, *gp;
    int datum;
    byte dp;
    extern int infd;

    if( data->segy ){
        UIMessage("loading segy data ...");
    } else
        UIMessage("loading float data ...");

    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;
    }

    {
        extern int _alloc;

        buf1 = (float *) malloc((size_t) (n1 * n2) * sizeof(buf1[0]));
        _alloc += (n1 * n2) * sizeof(buf1[0]);
        if( buf1 == 0 ){
            err("cant allocate %d bytes for  buf1; %d already allocated",
                (n1 * n2) * sizeof(buf1[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " buf1",
                          (n1 * n2) * sizeof(buf1[0]));};
        }
    {
        extern int _alloc;

        buf2 = (float *) malloc((size_t) (n1 * n2) * sizeof(buf2[0]));
        _alloc += (n1 * n2) * sizeof(buf2[0]);
        if( buf2 == 0 ){
            err("cant allocate %d bytes for  buf2; %d already allocated",
                (n1 * n2) * sizeof(buf2[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " buf2",
                          (n1 * n2) * sizeof(buf2[0]));};
        }
    {
        extern int _alloc;

        tgain = (float *) malloc((size_t) (n1) * sizeof(tgain[0]));
        _alloc += (n1) * sizeof(tgain[0]);
        if( tgain == 0 ){
            err("cant allocate %d bytes for  tgain; %d already allocated",
                (n1) * sizeof(tgain[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " tgain", (n1) * sizeof(tgain[0]));
        }
    };

    size = (n1 + head) > DATA_HEAD1 ? (n1 + head) : DATA_HEAD1;
    size = size > data->hbytes ? size : data->hbytes;
    {
        extern int _alloc;

        trace = (float *) malloc((size_t) (size) * sizeof(trace[0]));
        _alloc += (size) * sizeof(trace[0]);
        if( trace == 0 ){
            err("cant allocate %d bytes for  trace; %d already allocated",
                (size) * sizeof(trace[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " trace",
                          (size) * sizeof(trace[0]));};
        }

    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]));

        if( data->axis[DATA_VALUE] ){
            free(data->axis[DATA_VALUE]);
            data->axis[DATA_VALUE] = 0;
            if( memwatch ){
                printf("free %s\n", "data->axis[DATA_VALUE]");
            }
        };
        data->axis[DATA_VALUE] =