data.c

seismic software,very useful

            AxisInit2(DATA_VALUE, 1, label, DATA_VALUE_SIZE + 1, data->low,
                      (data->high - data->low) / DATA_VALUE_SIZE, 1.0);
    }
    if( buf2 ){
        free(buf2);
        buf2 = 0;
        if( memwatch ){
            printf("free %s\n", "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;
    }
    if( buf1 ){
        free(buf1);
        buf1 = 0;
        if( memwatch ){
            printf("free %s\n", "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);

    if( trace ){
        free(trace);
        trace = 0;
        if( memwatch ){
            printf("free %s\n", "trace");
        }
    };
    if( tgain ){
        free(tgain);
        tgain = 0;
        if( memwatch ){
            printf("free %s\n", "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(float *x, int n, float p)
{
    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);
}

/* return long data name */
char    *DataLabel(Data data)
{
    static Message message;

    if( !data ){
        return (0);
    }
    sprintf(message, "%s: %dx%dx%dx%dx%d=%d samples ", data->title,
            AxisSize(data->axis[DATA_AXIS1]),
            AxisSize(data->axis[DATA_AXIS2]),
            AxisSize(data->axis[DATA_AXIS3]),
            AxisSize(data->axis[DATA_AXIS4]),
            AxisSize(data->axis[DATA_AXIS5]), data->size);
    return (message);
}

/* return short data name */
char    *DataTitle(Data data)
{
    if( !data ){
        return (0);
    }
    return (data->title);
}

/* return data file */
char    *DataFile(Data data)
{
    if( !data ){
        return (0);
    }
    return (data->file);
}

/* return short name of dataset from filename */
char    *DataShortName(Data data)
{
    static string name;
    char    *c, *strrchr(const char *, int), *strchr(const char *, int);

    /* truncate leading pathname */
    if( (c = strrchr(data->file, '/')) != NULL ){
        strcpy(name, c + 1);
    } else {
        strcpy(name, data->file);
    }
    /* truncate trailing tail */
    if( (c = strchr(name, '.')) != NULL ){
        *c = '\0';
    }
    return (name);
}

/* return data buffer */
Buffer   DataBuffer(Data data)
{
    if( !data ){
        return (0);
    }
    return (data->buffer);
}

/* return axis reference */
Axis     DataAxis(Data data, int iaxis)
{
    if( !data || iaxis > DATA_NAXIS ){
        return (0);
    }
    return (data->axis[iaxis]);
}

/* return data size */
int DataSize(Data data)
{
    return (data->size);
}

/* return longest dimension */
int DataMaxDim(Data data)
{
    int      max, iaxis;

    if( !data ){
        return (0);
    }
    for( iaxis = 0, max = 1; iaxis < DATA_NAXIS - 1; iaxis++ ){
        max =
            AxisSize(data->axis[iaxis + 1]) >
            max ? AxisSize(data->axis[iaxis + 1]) : max;
    }
    return (max);
}

/* return data amplitude */
float    DataValue(Data data, int value)
{
    return (data->low +
            ((data->high - data->low) * ((value & 0x7F) - data->value_base)) /
            DATA_VALUE_SIZE);
}

/* return data value size */
int DataValueSize(Data data)
{
    if( !data ){
        return (0);
    } else {
        return (data->value_size);
    }
}

/* return base data value */
int DataValueBase(Data data)
{
    if( !data ){
        return (0);
    } else {
        return (data->value_base);
    }
}

/* return data high clip */
float    DataHigh(void)
{
    extern Data data;

    if( !data ){
        return (0.0);
    }
    return (data->high);
}

/* return data low clip */
float    DataLow(void)
{
    extern Data data;

    if( !data ){
        return (0.0);
    }
    return (data->low);
}

/* return integer percentage of histogram */
int DataHistogram(int i)
{
    extern Data data;

    if( !data || i < 0 || i > 255 ){
        return (0);
    }
    return ((int) (100 * data->histogram[i]));
}

void DataComputeHistogram(Data data)
{
    register Buffer bp, be;
    int      i;
    float    max;

    /*----------------------------*/
    /* clear the histogram buffer */

    /*----------------------------*/

    for( i = 0; i < 256; i++ ){
        data->histogram[i] = 0;
    }

    /*------------------------------*/
    /* calculate total for each bin */

    /*------------------------------*/

    for( bp = data->buffer, be = bp + data->size; bp < be; bp++ ){

        data->histogram[*bp] += 1.0;
    }

    /*------------------------------*/
    /* rescale totals to percentage */

    /*------------------------------*/

    max = 0;

    for( i = 0; i < ColorSize(); i++ ){
        max = data->histogram[i] > max ? data->histogram[i] : max;
    }

    for( i = 0; i < ColorSize(); i++ ){
        data->histogram[i] /= max;
    }
}

/* print data information */
void DataInfo(Data data)
{
    Message  message;

    sprintf(message,
            "Data: %s: in=%s %dx%dx%dx%dx%d=%d esize=%d segy=%d hbytes=%d script=%s",
            data->title, data->file, AxisSize(data->axis[DATA_AXIS1]),
            AxisSize(data->axis[DATA_AXIS2]), AxisSize(data->axis[DATA_AXIS3]),
            AxisSize(data->axis[DATA_AXIS4]), AxisSize(data->axis[DATA_AXIS5]),
            data->size, data->esize, data->segy, data->hbytes, data->script);
    UIMessage(message);
}

/* data value parameters */
void DataValueInfo(Data data)
{
    Message  message;

    sprintf(message,
            "Value: min=%g low=%g high=%g max=%g tpow=%g gpow=%g base=%d size=%d",
            data->min, data->low, data->high, data->max, data->tpow, data->gpow,
            data->value_base, data->value_size);
    UIMessage(message);
}

/* save data parameters */
void DataSavePar(Data data)
{
    int      iaxis;
    Message  message;

    if( !data ){
        return;
    }
    sprintf(message,
            "Data: %s: in=%s %dx%dx%dx%dx%d=%d esize=%d segy=%d hbytes=%d script=%s",
            data->title, data->file, AxisSize(data->axis[DATA_AXIS1]),
            AxisSize(data->axis[DATA_AXIS2]), AxisSize(data->axis[DATA_AXIS3]),
            AxisSize(data->axis[DATA_AXIS4]), AxisSize(data->axis[DATA_AXIS5]),
            data->size, data->esize, data->segy, data->hbytes, data->script);
    UISaveMessage(message);
    sprintf(message,
            "Value: min=%g low=%g high=%g max=%g tpow=%g gpow=%g base=%d size=%d",
            data->min, data->low, data->high, data->max, data->tpow, data->gpow,
            data->value_base, data->value_size);
    UISaveMessage(message);
    for( iaxis = DATA_AXIS1; iaxis <= DATA_AXIS5; iaxis++ ){
        AxisSavePar(data->axis[iaxis]);
    }
}

/* dump data bytes */
void DataDumpBytes(Data data, char *file, int fd)
{
    register byte bp, be;

    DrawWatch(1);
    for( bp = data->buffer, be = bp + data->size; bp < be; ){
        *bp++ <<= 1;
    }
    write(fd, data->buffer, data->size);
    for( bp = data->buffer, be = bp + data->size; bp < be; ){
        *bp++ >>= 1;
    }
    DataDumpHeader(data, file, fd, 1);
    DrawWatch(0);
}

/* dump data vgrid header */
void DataDumpHeader(Data data, char *file, int datafd, int esize)
{
    Message  message;
    char     parfile[80];
    extern FILE *savefd;

    if( data->gh.dtype == 0. ){ /* build a header */
        data->gh.scale = 1.;
        data->gh.n1 = AxisSize(data->axis[DATA_AXIS1]) * data->gh.scale;
        data->gh.n2 = AxisSize(data->axis[DATA_AXIS2]) * data->gh.scale;
        data->gh.n3 = AxisSize(data->axis[DATA_AXIS3]) * data->gh.scale;
        data->gh.n4 = AxisSize(data->axis[DATA_AXIS4]) * data->gh.scale;
        data->gh.n5 = AxisSize(data->axis[DATA_AXIS5]) * data->gh.scale;
        data->gh.d1 = AxisDelta(data->axis[DATA_AXIS1]) * data->gh.scale;
        data->gh.d2 = AxisDelta(data->axis[DATA_AXIS2]) * data->gh.scale;
        data->gh.d3 = AxisDelta(data->axis[DATA_AXIS3]) * data->gh.scale;
        data->gh.d4 = AxisDelta(data->axis[DATA_AXIS4]) * data->gh.scale;
        data->gh.d5 = AxisDelta(data->axis[DATA_AXIS5]) * data->gh.scale;
        data->gh.o1 = AxisFirst(data->axis[DATA_AXIS1]) * data->gh.scale;
        data->gh.o2 = AxisFirst(data->axis[DATA_AXIS2]) * data->gh.scale;
        data->gh.o3 = AxisFirst(data->axis[DATA_AXIS3]) * data->gh.scale;
        data->gh.o4 = AxisFirst(data->axis[DATA_AXIS4]) * data->gh.scale;
        data->gh.o5 = AxisFirst(data->axis[DATA_AXIS5]) * data->gh.scale;
        data->gh.dcdp2 = 1.;
        data->gh.dline3 = 1.;
        data->gh.ocdp2 = 1.;
        data->gh.oline3 = 1.;
        data->gh.gmin = data->low * data->gh.scale;
        data->gh.gmax = data->high * data->gh.scale;
    }
    data->gh.dtype = esize * data->gh.scale;
    write(datafd, &data->gh, sizeof(data->gh));
    close(datafd);
    sprintf(parfile, "%s.H", file);
    if( (savefd = fopen(parfile, "w")) == NULL ){
       UIMessage("cant create vgrid parfile");
    }
    DataSavePar(data);
    sprintf(message, "Byte data dump into %s", parfile);
    UIMessage(message);
    fclose(savefd);
}

/* dump data floats */
void DataDumpFloats(Data data, char *file, int fd)
{
    int      nn, n3, i;
    Buffer   dbuf;
    float   *fbuf, min, max;

    DrawWatch(1);
    nn = AxisSize(data->axis[DATA_AXIS1]) * AxisSize(data->axis[DATA_AXIS2]);
    n3 = data->size / nn;

    {
        extern int _alloc;

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

    /* convert and dump a plane at a time */
    min = data->high;
    max = data->low;
    for( i = 0, dbuf = data->buffer; i < n3; i++, dbuf += nn ){
        Data2Float(dbuf, fbuf, nn, &min, &max);
        write(fd, fbuf, nn * sizeof(fbuf[0]));
    }
    data->gh.gmin = min * data->gh.scale;
    data->gh.gmax = max * data->gh.scale;
    DataDumpHeader(data, file, fd, 4);

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

void Data2Float(Buffer dbuf, float *fbuf, int n, float *min, float *max)
{
    register Buffer bp, be;
    register float *fp, *mp, Min, Max;
    extern Data data;

    Min = *min;
    Max = *max;
    /* convert */
    for( mp = AxisValues(data->axis[DATA_VALUE]) - DATA_VALUE_BASE, bp =
         dbuf, fp = fbuf, be = dbuf + n; bp < be; bp++, fp++ ){
        *fp = mp[*bp];
        Max = Max > *fp ? Max : *fp;
        Min = Min < *fp ? Min : *fp;
    }
    *min = Min;
    *max = Max;
}

int      irint(float x)
{
    int      i;

    i = x;
    return ((x - i) < (1 - (x - i)) ? i : i + 1);
}