fpcompress.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <unistd.h>
#include "raster3d_intern.h"

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

#define XDR_DOUBLE_LENGTH 8
#define XDR_DOUBLE_NOF_EXP_BYTES 2
#define XDR_FLOAT_LENGTH 4
#define XDR_FLOAT_NOF_EXP_BYTES 1

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

void Rast3d_fpcompress_print_binary(char *c, int numBits)
{
    unsigned char bit;
    int i;

    bit = 1 << (numBits - 1);

    for (i = 0; i < numBits; i++) {
        printf("%d", (*((unsigned char *)c) & bit) != 0);
        bit >>= 1;
    }
}

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

void Rast3d_fpcompress_dissect_xdr_double(unsigned char *numPointer)
{
    char sign, exponent;

    sign = *numPointer >> 7;
    exponent = (*numPointer << 1) | (*(numPointer + 1) >> 7);

    printf("%f: sign = ", *((float *)numPointer));
    Rast3d_fpcompress_print_binary(&sign, 1);
    printf("   exp = ");
    Rast3d_fpcompress_print_binary(&exponent, 8);
    printf("   mantissa = ");
    Rast3d_fpcompress_print_binary((char *)(numPointer + 1), 7);
    Rast3d_fpcompress_print_binary((char *)(numPointer + 2), 8);
    Rast3d_fpcompress_print_binary((char *)(numPointer + 3), 8);
    printf("\n");
}

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

static unsigned char clearMask[9] =
    { 255, 128, 192, 224, 240, 248, 252, 254, 255 };

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

#define ALL_NULL_CODE 2
#define ZERO_NULL_CODE 1
#define SOME_NULL_CODE 0

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

static void
G_fpcompress_rearrangeEncodeFloats(unsigned char *src, int size,
                                   int precision, unsigned char *dst,
                                   int *length, int *offsetMantissa)
{
    unsigned int nNullBits, nBits;
    register unsigned char *srcStop;
    register unsigned char *cp0, *cp1, *cp2, *cp3, *nullBits;
    unsigned char mask, isNull;
    int gt8, gt16, srcIncrement, nofNull;
    float *f;

    srcStop = src + size * XDR_FLOAT_LENGTH;

    if ((precision >= 23) || (precision == -1)) {
        cp3 = dst;
        cp2 = cp3 + size;
        cp1 = cp2 + size;
        cp0 = cp1 + size;
        while (srcStop != src) {
            *cp3++ = *src++;    /* sign + 7 exponent bits */
            *cp2++ = *src++;    /* 1 exponent bit + 7 ms mantissa bits */
            *cp1++ = *src++;    /* 8 mantissa bits */
            *cp0++ = *src++;    /* 8 ls mantissa bits */
        }

        *length = size * XDR_FLOAT_LENGTH;
        *offsetMantissa = size;

        return;
    }

    f = (float *)src;
    nofNull = 0;
    while (srcStop != (unsigned char *)f)
        nofNull += Rast3d_is_xdr_null_float(f++);

    if (nofNull == size) {
        *dst = (unsigned char)ALL_NULL_CODE;

        *length = 1;
        *offsetMantissa = 1;

        return;
    }

    precision += 1;             /* treat the ls exponent bit like an addl mantissa bit */

    *dst = (unsigned char)(nofNull == 0 ? ZERO_NULL_CODE : SOME_NULL_CODE);

    gt16 = precision > 16;
    gt8 = precision > 8;
    srcIncrement = 1 + (!gt8) + (!gt16);

    precision %= 8;

    nullBits = dst + 1;
    if (nofNull)
        cp0 = nullBits + size / 8 + ((size % 8) != 0);
    else
        cp0 = nullBits;
    cp3 = cp0 + size - nofNull;
    cp2 = cp3 + size - nofNull;
    cp1 = cp3 + (gt8 + gt16) * (size - nofNull);

    mask = clearMask[precision];
    nBits = nNullBits = 0;

    while (srcStop != src) {
        if (nofNull) {
            isNull = Rast3d_is_xdr_null_float((float *)src);

            if (nNullBits) {
                *nullBits |= ((unsigned char)isNull << nNullBits++);
                if (nNullBits == 8) {
                    nullBits++;
                    nNullBits = 0;
                }
            }
            else {
                *nullBits = (unsigned char)isNull;
                nNullBits++;
            }

            if (isNull) {
                src += XDR_FLOAT_LENGTH;
                continue;
            }
        }

        /* printf ("write src cp0 %d %d (%d %d) %d\n", *src, *cp0, src, cp0, nullBits); */

        *cp0++ = *src++;
        if (gt8)
            *cp3++ = *src++;
        if (gt16)
            *cp2++ = *src++;

        if (nBits && precision) {
            *cp1 |= (unsigned char)((unsigned char)(*src & mask) >> nBits);

            /*printf ("%d\n", ((*src & mask) >> nBits) << nBits); */

            if (8 - nBits < precision) {
                cp1++;

                /*printf ("%d %d\n", *cp1, (*src & mask) << (8 - nBits)); */

                *cp1 =
                    (unsigned char)((unsigned char)((*src & mask)) <<
                                    (8 - nBits));
                nBits += precision - 8;
            }
            else {
                nBits = (nBits + precision) % 8;
                if (!nBits)
                    cp1++;
            }
        }
        else {
            *cp1 = (unsigned char)(*src & mask);
            /* printf ("%d %d %d\n", *cp1, *src, nBits); */
            nBits = (nBits + precision) % 8;
            if (!nBits)
                cp1++;
        }

        src += srcIncrement;
    }

    *length = 1;                /* null-bit-vector indicator-byte */

    if (nofNull)                /* length of null-bit-vector */
        *length += size / 8 + ((size % 8) != 0);

    /* length of data */
    *length += (gt8 + gt16 + (precision == 0) + 1) * (size - nofNull) +
        ((precision * (size - nofNull)) / 8) +
        (((precision * (size - nofNull)) % 8) != 0);

    *offsetMantissa = size - nofNull;
}

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

static void
G_fpcompress_rearrangeEncodeDoubles(unsigned char *src, int size,
                                    int precision, unsigned char *dst,
                                    int *length, int *offsetMantissa)
{
    unsigned int nNullBits, nBits;
    unsigned char isNull;
    register unsigned char *srcStop;
    register unsigned char *cp0, *cp1, *cp2, *cp3;
    register unsigned char *cp4, *cp5, *cp6, *cp7, *nullBits;
    unsigned char mask;
    int gt8, gt16, gt24, gt32, gt40, gt48, srcIncrement, nofNull;
    double *d;

    srcStop = src + size * XDR_DOUBLE_LENGTH;

    if ((precision >= 52) || (precision == -1)) {
        cp7 = dst;
        cp6 = cp7 + size;
        cp5 = cp6 + size;
        cp4 = cp5 + size;
        cp3 = cp4 + size;
        cp2 = cp3 + size;
        cp1 = cp2 + size;
        cp0 = cp1 + size;

        while (srcStop != src) {
            *cp7++ = *src++;    /* sign + 7 ms exponent bits */
            *cp6++ = *src++;    /* 4 exponent bits + 4 ms mantissa bits */
            *cp5++ = *src++;    /* 8 mantissa bits */
            *cp4++ = *src++;
            *cp3++ = *src++;
            *cp2++ = *src++;
            *cp1++ = *src++;
            *cp0++ = *src++;    /* 8 ls mantissa bits */
        }

        *length = size * XDR_DOUBLE_LENGTH;
        *offsetMantissa = size;

        return;
    }

    precision += 4;             /* treat the 4 ls exponent bits like addl mantissa bits */

    d = (double *)src;
    nofNull = 0;
    while (srcStop != (unsigned char *)d)
        nofNull += Rast3d_is_xdr_null_double(d++);

    if (nofNull == size) {
        *dst = (unsigned char)ALL_NULL_CODE;

        *length = 1;
        *offsetMantissa = 1;

        return;
    }

    *dst = (unsigned char)(nofNull == 0 ? ZERO_NULL_CODE : SOME_NULL_CODE);

    gt48 = precision > 48;
    gt40 = precision > 40;
    gt32 = precision > 32;
    gt24 = precision > 24;
    gt16 = precision > 16;
    gt8 = precision > 8;
    srcIncrement =
        1 + (!gt8) + (!gt16) + (!gt24) + (!gt32) + (!gt40) + (!gt48);

    precision %= 8;

    nullBits = dst + 1;
    if (nofNull)
        cp0 = nullBits + size / 8 + ((size % 8) != 0);
    else
        cp0 = nullBits;
    cp7 = cp0 + size - nofNull;
    cp6 = cp7 + size - nofNull;
    cp5 = cp6 + size - nofNull;
    cp4 = cp5 + size - nofNull;
    cp3 = cp4 + size - nofNull;
    cp2 = cp3 + size - nofNull;
    cp1 = cp7 + (gt8 + gt16 + gt24 + gt32 + gt40 + gt48) * (size - nofNull);

    mask = clearMask[precision];
    nBits = nNullBits = 0;

    while (srcStop != src) {
        if (nofNull) {
            isNull = Rast3d_is_xdr_null_double((double *)src);

            if (nNullBits) {
                *nullBits |= ((unsigned char)isNull << nNullBits++);
                if (nNullBits == 8) {
                    nullBits++;
                    nNullBits = 0;
                }
            }
            else {
                *nullBits = (unsigned char)isNull;
                nNullBits++;
            }

            if (isNull) {
                src += XDR_DOUBLE_LENGTH;
                continue;
            }
        }

        *cp0++ = *src++;
        if (gt32) {
            *cp7++ = *src++;
            *cp6++ = *src++;
            *cp5++ = *src++;

            if (gt32)
                *cp4++ = *src++;
            if (gt40)
                *cp3++ = *src++;
            if (gt48)
                *cp2++ = *src++;
        }
        else {
            if (gt8)
                *cp7++ = *src++;
            if (gt16)
                *cp6++ = *src++;
            if (gt24)
                *cp5++ = *src++;
        }

        if (nBits && precision) {
            *cp1 |= (unsigned char)((unsigned char)(*src & mask) >> nBits);
            if (8 - nBits < precision) {
                cp1++;
                *cp1 =
                    (unsigned char)(((unsigned char)(*src & mask)) <<
                                    (8 - nBits));
                nBits += precision - 8;
            }
            else {
                nBits = (nBits + precision) % 8;
                if (!nBits)
                    cp1++;
            }
        }
        else {
            *cp1 = (unsigned char)(*src & mask);
            nBits = (nBits + precision) % 8;
            if (!nBits)
                cp1++;
        }

        src += srcIncrement;
    }

    *length = 1;

    if (nofNull)
        *length += size / 8 + ((size % 8) != 0);

    *length +=
        (1 + gt8 + gt16 + gt24 + gt32 + gt40 + gt48 +
         (precision ==
          0)) * (size - nofNull) + ((precision * (size - nofNull)) / 8) +
        (((precision * (size - nofNull)) % 8) != 0);

    if (gt8)
        *offsetMantissa = 2 * (size - nofNull);
    else
        *offsetMantissa = *length;
}

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

static void
G_fpcompress_rearrangeDecodeFloats(unsigned char *src, int size,
                                   int precision, unsigned char *dst)
{
    unsigned int nNullBits, nBits;
    register unsigned char *dstStop;
    register unsigned char *cp0, *cp1, *cp2, *cp3, *nullBits;
    unsigned char mask, isNull;
    int gt8, gt16, dstIncrement, nofNull;
    float *f, *fStop;

    if ((precision != -1) && (precision <= 15)) {       /* 23 - 8 */
        cp3 = dst + 3;
        dstStop = dst + XDR_FLOAT_LENGTH * size + 3;
        while (dstStop != cp3) {
            *cp3 = 0;
            cp3 += XDR_FLOAT_LENGTH;
        }

        if (precision <= 7) {
            cp3 = dst + 2;
            dstStop = dst + XDR_FLOAT_LENGTH * size + 2;
            while (dstStop != cp3) {
                *cp3 = 0;
                cp3 += XDR_FLOAT_LENGTH;
            }
        }
    }

    dstStop = dst + size * XDR_FLOAT_LENGTH;

    if ((precision >= 23) || (precision == -1)) {
        cp3 = src;
        cp2 = cp3 + size;
        cp1 = cp2 + size;
        cp0 = cp1 + size;
        while (dstStop != dst) {
            *dst++ = *cp3++;
            *dst++ = *cp2++;
            *dst++ = *cp1++;
            *dst++ = *cp0++;
        }

        return;
    }

    if (*src == (unsigned char)ALL_NULL_CODE) {
        f = (float *)dst;
        while (dstStop != (unsigned char *)f)
            Rast3d_set_xdr_null_float(f++);

        return;
    }

    precision += 1;             /* treat the ls exponent bit like an addl mantissa bit */

    gt16 = precision > 16;
    gt8 = precision > 8;
    dstIncrement = 1 + (!gt8) + (!gt16);

    precision %= 8;

    nofNull = 0;
    nullBits = src + 1;
    nNullBits = 0;
    if (*src == (unsigned char)SOME_NULL_CODE) {
        f = (float *)src;
        fStop = (float *)(src + size * XDR_FLOAT_LENGTH);
        while (fStop != f++) {
            nofNull += ((*nullBits & ((unsigned char)1 << nNullBits++)) != 0);
            if (nNullBits == 8) {
                nullBits++;
                nNullBits = 0;
            }
        }
    }

    nullBits = src + 1;
    if (nofNull)
        cp0 = nullBits + size / 8 + ((size % 8) != 0);
    else
        cp0 = nullBits;
    cp3 = cp0 + size - nofNull;
    cp2 = cp3 + size - nofNull;
    cp1 = cp3 + (gt8 + gt16) * (size - nofNull);

    mask = clearMask[precision];
    nBits = nNullBits = 0;

    while (dstStop != dst) {
        if (nofNull) {
            isNull = *nullBits & ((unsigned char)1 << nNullBits++);

            if (nNullBits == 8) {
                nullBits++;
                nNullBits = 0;
            }

            if (isNull) {
                Rast3d_set_xdr_null_float((float *)dst);
                dst += XDR_FLOAT_LENGTH;
                continue;
            }
        }

        *dst++ = *cp0++;
        if (gt8)
            *dst++ = *cp3++;
        if (gt16)
            *dst++ = *cp2++;

        if (nBits && precision) {
            *dst = (unsigned char)((*cp1 << nBits) & mask);

            if (8 - nBits < precision) {
                cp1++;
                *dst |= (unsigned char)((*cp1 >> (8 - nBits)) & mask);
                nBits += precision - 8;
            }
            else {
                nBits = (nBits + precision) % 8;
                if (!nBits)
                    cp1++;
            }
        }
        else {
            *dst = (unsigned char)(*cp1 & mask);
            nBits = (nBits + precision) % 8;
            if (!nBits)
                cp1++;
        }

        dst += dstIncrement;
    }
}

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

static void
G_fpcompress_rearrangeDecodeDoubles(unsigned char *src, int size,
                                    int precision, unsigned char *dst)
{
    unsigned int nNullBits, nBits;
    register unsigned char *dstStop;
    register unsigned char *cp0, *cp1, *cp2, *cp3;
    register unsigned char *cp4, *cp5, *cp6, *cp7, *nullBits;
    unsigned char mask, isNull;
    int gt8, gt16, gt24, gt32, gt40, gt48, dstIncrement, offs, nofNull;
    double *d, *dStop;

    if ((precision != -1) && (precision <= 44)) {
        for (offs = 7; offs >= (precision + 19) / 8; offs--) {
            cp7 = dst + offs;
            dstStop = dst + XDR_DOUBLE_LENGTH * size + offs;
            while (dstStop != cp7) {
                *cp7 = 0;
                cp7 += XDR_DOUBLE_LENGTH;
            }
        }
    }

    dstStop = dst + size * XDR_DOUBLE_LENGTH;

    if ((precision >= 52) || (precision == -1)) {
        cp7 = src;
        cp6 = cp7 + size;
        cp5 = cp6 + size;
        cp4 = cp5 + size;
        cp3 = cp4 + size;
        cp2 = cp3 + size;
        cp1 = cp2 + size;
        cp0 = cp1 + size;

        while (dstStop != dst) {
            *dst++ = *cp7++;
            *dst++ = *cp6++;
            *dst++ = *cp5++;
            *dst++ = *cp4++;
            *dst++ = *cp3++;
            *dst++ = *cp2++;
            *dst++ = *cp1++;
            *dst++ = *cp0++;
        }

        return;
    }

    if (*src == (unsigned char)ALL_NULL_CODE) {
        /*printf ("all null\n"); */
        d = (double *)dst;
        while (dstStop != (unsigned char *)d)
            Rast3d_set_xdr_null_double(d++);

        return;
    }

    precision += 4;             /* treat the 4 ls exponent bits like addl mantissa bits */

    gt48 = precision > 48;
    gt40 = precision > 40;
    gt32 = precision > 32;
    gt24 = precision > 24;
    gt16 = precision > 16;
    gt8 = precision > 8;

    dstIncrement =
        1 + (!gt8) + (!gt16) + (!gt24) + (!gt32) + (!gt40) + (!gt48);

    precision %= 8;

    nofNull = 0;
    nullBits = src + 1;
    nNullBits = 0;
    if (*src == (unsigned char)SOME_NULL_CODE) {
        d = (double *)src;
        dStop = (double *)(src + size * XDR_DOUBLE_LENGTH);
        while (dStop != d++) {
            nofNull += ((*nullBits & ((unsigned char)1 << nNullBits++)) != 0);
            if (nNullBits == 8) {
                nullBits++;
                nNullBits = 0;
            }
        }
    }

    nullBits = src + 1;
    if (nofNull)
        cp0 = nullBits + size / 8 + ((size % 8) != 0);
    else
        cp0 = nullBits;
    cp7 = cp0 + size - nofNull;
    cp6 = cp7 + size - nofNull;
    cp5 = cp6 + size - nofNull;
    cp4 = cp5 + size - nofNull;
    cp3 = cp4 + size - nofNull;
    cp2 = cp3 + size - nofNull;
    cp1 = cp7 + (gt8 + gt16 + gt24 + gt32 + gt40 + gt48) * (size - nofNull);

    mask = clearMask[precision];
    nBits = nNullBits = 0;

    while (dstStop != dst) {
        if (nofNull) {
            isNull = *nullBits & ((unsigned char)1 << nNullBits++);

            if (nNullBits == 8) {
                nullBits++;
                nNullBits = 0;
            }

            if (isNull) {
                Rast3d_set_xdr_null_double((double *)dst);
                dst += XDR_DOUBLE_LENGTH;
                continue;
            }
        }

        *dst++ = *cp0++;
        if (gt32) {
            *dst++ = *cp7++;
            *dst++ = *cp6++;
            *dst++ = *cp5++;

            if (gt32)
                *dst++ = *cp4++;
            if (gt40)
                *dst++ = *cp3++;
            if (gt48)
                *dst++ = *cp2++;
        }
        else {
            if (gt8)
                *dst++ = *cp7++;
            if (gt16)
                *dst++ = *cp6++;
            if (gt24)
                *dst++ = *cp5++;
        }

        if (nBits && precision) {
            *dst = (unsigned char)((*cp1 << nBits) & mask);

            if (8 - nBits < precision) {
                cp1++;
                *dst |= (unsigned char)((*cp1 >> (8 - nBits)) & mask);
                nBits += precision - 8;
            }
            else {
                nBits = (nBits + precision) % 8;
                if (!nBits)
                    cp1++;
            }
        }
        else {
            *dst = (unsigned char)(*cp1 & mask);
            nBits = (nBits + precision) % 8;
            if (!nBits)
                cp1++;
        }

        dst += dstIncrement;
    }
}

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

int
Rast3d_fpcompress_write_xdr_nums(int fd, char *src, int nofNum, int precision,
                          char *compressBuf, int isFloat)
{
    int status;
    int nBytes;
    int offsetMantissa;

    if (isFloat)
        G_fpcompress_rearrangeEncodeFloats((unsigned char *)src, nofNum, precision,
                        (unsigned char *)(compressBuf + 1),
                                           &nBytes, &offsetMantissa);
    else
        G_fpcompress_rearrangeEncodeDoubles((unsigned char *)src, nofNum, precision,
                                (unsigned char *)(compressBuf + 1),
                                            &nBytes, &offsetMantissa);

        *compressBuf = 0;
        status = G_write_compressed(fd, (unsigned char *)compressBuf, nBytes + 1, 2);

    if (status < 0) {
        Rast3d_error("Rast3d_fpcompress_write_xdr_nums: write error");
        return 0;
    }

    return 1;
}

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

int
Rast3d_fpcompress_read_xdr_nums(int fd, char *dst, int nofNum, int fileBytes,
                         int precision, char *compressBuf, int isFloat)
{
    int status;
    int lengthEncode, lengthDecode;
    int nBytes;
    char *src, *dest, *srcStop;
    nBytes = (isFloat ? XDR_FLOAT_LENGTH : XDR_DOUBLE_LENGTH);

    status = G_read_compressed(fd, fileBytes, (unsigned char *)compressBuf,
                         nofNum * nBytes + 1, 2);

    if (status < 0) {
        Rast3d_error("Rast3d_fpcompress_read_xdr_nums: read error");
        return 0;
    }

    /* This code is kept for backward compatibility */
    if (*compressBuf++ == 1) {
        status--;
        Rast3d_rle_decode(compressBuf, dst, nofNum * nBytes, 1,
                     &lengthEncode, &lengthDecode);
        if (*dst == ALL_NULL_CODE)
            Rast3d_fatal_error("Rast3d_fpcompress_read_xdr_nums: wrong code");

        if (status == nofNum * nBytes)
            status -= lengthDecode - lengthEncode;

        src = compressBuf + status - 1;
        srcStop = compressBuf + lengthEncode - 1;
        dest = compressBuf + (status - lengthEncode) + lengthDecode - 1;
        while (src != srcStop)
            *dest-- = *src--;

        src = dst;
        srcStop = src + lengthDecode;
        dest = compressBuf;
        while (src != srcStop)
            *dest++ = *src++;
    }

    if (isFloat)
        G_fpcompress_rearrangeDecodeFloats((unsigned char *)compressBuf, nofNum, precision,
                        (unsigned char *)dst);
    else
        G_fpcompress_rearrangeDecodeDoubles((unsigned char *)compressBuf, nofNum, precision,
                        (unsigned char *)dst);

    return 1;
}