/*

This file is taken from V^3, the Versatile Volume Viewer.
Copyright (c) 2003-2007 by Stefan Roettger.

The volume viewer is licensed under the terms of the GPL (see
http://www.gnu.org/copyleft/ for more information on the license).
Any commercial use of the code or parts of it requires the explicit
permission of the author! The source code is distributed with
ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY etc.!

The author's contact address is:

   mailto:stefan@stereofx.org
   http://stereofx.org

*/

// (c) by Stefan Roettger

#include "codebase.h" // universal code base
#include "ddsbase.h"

#define DDS_MAXSTR (256)

#define DDS_BLOCKSIZE (1<<20)
#define DDS_INTERLEAVE (1<<24)

#define DDS_RL (7)


//FIXME: endianness conversion seems to be broken, so we perform it in PVMVolumeReader. joerg
#define DDS_ISINTEL 1
    // (*((unsigned char *)(&DDS_INTEL)+1)==0)

#include "voreen/core/io/progressbar.h"
#include "tgt/exception.h"

FILE *DDS_file;

char DDS_ID[] = "DDS v3d\n";
char DDS_ID2[] = "DDS v3e\n";

unsigned int DDS_buffer;
int DDS_bufsize, DDS_bitcnt;

unsigned short int DDS_INTEL = 1;

inline unsigned int DDS_shiftl(const unsigned int value, const int bits) {
    return ((bits >= 32) ? 0 : value << bits);
}

inline unsigned int DDS_shiftr(const unsigned int value, const int bits) {
    return ((bits >= 32) ? 0 : value >> bits);
}

void initbuffer() {
    DDS_buffer = 0;
    DDS_bufsize = 0;
    DDS_bitcnt = 0;
}

void DDS_swapuint(unsigned int *x) {
    unsigned int tmp = *x;

    *x = ((tmp & 0xff) << 24) |
         ((tmp & 0xff00) << 8) |
         ((tmp & 0xff0000) >> 8) |
         ((tmp & 0xff000000) >> 24);
}

void writebits(FILE *file, unsigned int value, int bits) {
    if (bits < 0 || bits > 32)
        ERRORMSG();

    if (bits == 0)
        return ;

    value &= DDS_shiftl(1, bits) - 1;

    if (DDS_bufsize + bits < 32) {
        DDS_buffer = DDS_shiftl(DDS_buffer, bits) | value;
        DDS_bufsize += bits;
    } else {
        DDS_buffer = DDS_shiftl(DDS_buffer, 32 - DDS_bufsize);
        DDS_bufsize += bits - 32;
        DDS_buffer |= DDS_shiftr(value, DDS_bufsize);
        if (DDS_ISINTEL)
            DDS_swapuint(&DDS_buffer);
        if (fwrite(&DDS_buffer, 4, 1, file) != 1)
            ERRORMSG();
        DDS_buffer = value & (DDS_shiftl(1, DDS_bufsize) - 1);
    }

    DDS_bitcnt += bits;
}

void flushbits(FILE *file) {
    if (DDS_bufsize > 0) {
        DDS_buffer = DDS_shiftl(DDS_buffer, 32 - DDS_bufsize);
        if (DDS_ISINTEL)
            DDS_swapuint(&DDS_buffer);
        if (fwrite(&DDS_buffer, (DDS_bufsize + 7) / 8, 1, file) != 1)
            ERRORMSG();
        DDS_bitcnt += (32 - DDS_bufsize) & 7;
    }
}

unsigned int readbits(FILE *file, int bits) {
    unsigned int value;

    if (bits < 0 || bits > 32)
        ERRORMSG();

    if (bits == 0)
        return (0);

    if (bits < DDS_bufsize) {
        DDS_bufsize -= bits;
        value = DDS_shiftr(DDS_buffer, DDS_bufsize);
    } else {
        value = DDS_shiftl(DDS_buffer, bits - DDS_bufsize);
        DDS_buffer = 0;
        size_t read = fread(&DDS_buffer, 1, 4, file);
        if (read == 0)
            return value;//ERRORMSG();  //FIXME: we ignore this error, necessary to load the
                                        //       16-bit bonsai PVM files from The Volume  Library. joerg
        if (DDS_ISINTEL)
            DDS_swapuint(&DDS_buffer);
        DDS_bufsize += 32 - bits;
        value |= DDS_shiftr(DDS_buffer, DDS_bufsize);
    }

    DDS_buffer &= DDS_shiftl(1, DDS_bufsize) - 1;
    DDS_bitcnt += bits;

    return (value);
}

inline int DDS_code(int bits) {
    return (bits > 1 ? bits - 1 : bits);
}

inline int DDS_decode(int bits) {
    return (bits >= 1 ? bits + 1 : bits);
}

// deinterleave a byte stream
void deinterleave(unsigned char *data, unsigned int bytes, unsigned int skip, unsigned int block = 0, BOOLINT restore = FALSE) {
    unsigned int i, j, k;

    unsigned char *data2, *ptr;

    if (skip <= 1)
        return ;

    if (block == 0) {
        if ((data2 = (unsigned char *)malloc(bytes)) == NULL)
            ERRORMSG();

        if (!restore)
            for (ptr = data2, i = 0; i < skip; i++)
                for (j = i; j < bytes; j += skip)
                    * ptr++ = data[j];
        else
            for (ptr = data, i = 0; i < skip; i++)
                for (j = i; j < bytes; j += skip)
                    data2[j] = *ptr++;

        memcpy(data, data2, bytes);
    } else {
        if ((data2 = (unsigned char *)malloc((bytes < skip * block) ? bytes : skip * block)) == NULL)
            ERRORMSG();

        if (!restore) {
            for (k = 0; k < bytes / skip / block; k++) {
                for (ptr = data2, i = 0; i < skip; i++)
                    for (j = i; j < skip*block; j += skip)
                        * ptr++ = data[k * skip * block + j];

                memcpy(data + k*skip*block, data2, skip*block);
            }

            for (ptr = data2, i = 0; i < skip; i++)
                for (j = i; j < bytes - k*skip*block; j += skip)
                    * ptr++ = data[k * skip * block + j];

            memcpy(data + k*skip*block, data2, bytes - k*skip*block);
        } else {
            for (k = 0; k < bytes / skip / block; k++) {
                for (ptr = data + k * skip * block, i = 0; i < skip; i++)
                    for (j = i; j < skip*block; j += skip)
                        data2[j] = *ptr++;

                memcpy(data + k*skip*block, data2, skip*block);
            }

            for (ptr = data + k * skip * block, i = 0; i < skip; i++)
                for (j = i; j < bytes - k*skip*block; j += skip)
                    data2[j] = *ptr++;

            memcpy(data + k*skip*block, data2, bytes - k*skip*block);
        }
    }

    free(data2);
}

// interleave a byte stream
void interleave(unsigned char *data, unsigned int bytes, unsigned int skip, unsigned int block = 0) {
    deinterleave(data, bytes, skip, block, TRUE);
}

// write a Differential Data Stream
void writeDDSfile(char *filename, unsigned char *data, unsigned int bytes, unsigned int skip, unsigned int strip, int nofree) {
    int version = 1;

    unsigned char *ptr1, *ptr2;

    int pre1, pre2,
    act1, act2,
    tmp1, tmp2;

    unsigned int cnt, cnt1, cnt2;
    int bits, bits1, bits2;

    if (bytes < 1)
        ERRORMSG();

    if (bytes > DDS_INTERLEAVE)
        version = 2;

    if (skip < 1 || skip > 4)
        skip = 1;
    if (strip < 1 || strip > 65536)
        strip = 1;

    if ((DDS_file = fopen(filename, "wb")) == NULL)
        ERRORMSG();

    fprintf(DDS_file, (version == 1) ? DDS_ID : DDS_ID2, 0); //supress gcc warning
    //fprintf(DDS_file, (version == 1) ? DDS_ID : DDS_ID2);

    deinterleave(data, bytes, skip, DDS_INTERLEAVE);

    initbuffer();

    writebits(DDS_file, skip - 1, 2);
    writebits(DDS_file, strip++ -1, 16);

    ptr1 = ptr2 = data;
    pre1 = pre2 = 0;

    cnt = cnt1 = cnt2 = 0;
    bits = bits1 = bits2 = 0;

    while (cnt++ < bytes) {
        tmp1 = *ptr1++;
        if (cnt <= strip)
            act1 = tmp1 - pre1;
        else
            act1 = tmp1 - pre1 - *(ptr1 - strip) + *(ptr1 - strip - 1);
        pre1 = tmp1;

        while (act1 < -128)
            act1 += 256;
        while (act1 > 127)
            act1 -= 256;

        if (act1 <= 0)
            for (bits = 0; (1 << bits) / 2 < -act1; bits++)
                ;
        else
            for (bits = 0; (1 << bits) / 2 <= act1; bits++)
                ;

        bits = DDS_decode(DDS_code(bits));

        if (cnt1 == 0) {
            cnt1++;
            bits1 = bits;
            continue;
        }

        if (cnt1 < (1 << DDS_RL) - 1 && bits == bits1) {
            cnt1++;
            continue;
        }

        if (cnt1 + cnt2 < (1 << DDS_RL) && (cnt1 + cnt2)*max(bits1, bits2) < cnt1*bits1 + cnt2*bits2 + DDS_RL + 3) {
            cnt2 += cnt1;
            if (bits1 > bits2)
                bits2 = bits1;
        } else {
            writebits(DDS_file, cnt2, DDS_RL);
            writebits(DDS_file, DDS_code(bits2), 3);

            while (cnt2-- > 0) {
                tmp2 = *ptr2++;
                if (ptr2 - strip <= data)
                    act2 = tmp2 - pre2;
                else
                    act2 = tmp2 - pre2 - *(ptr2 - strip) + *(ptr2 - strip - 1);
                pre2 = tmp2;

                while (act2 < -128)
                    act2 += 256;
                while (act2 > 127)
                    act2 -= 256;

                writebits(DDS_file, act2 + (1 << bits2) / 2, bits2);
            }

            cnt2 = cnt1;
            bits2 = bits1;
        }

        cnt1 = 1;
        bits1 = bits;
    }

    if (cnt1 + cnt2 < (1 << DDS_RL) && (cnt1 + cnt2)*max(bits1, bits2) < cnt1*bits1 + cnt2*bits2 + DDS_RL + 3) {
        cnt2 += cnt1;
        if (bits1 > bits2)
            bits2 = bits1;
    } else {
        writebits(DDS_file, cnt2, DDS_RL);
        writebits(DDS_file, DDS_code(bits2), 3);

        while (cnt2-- > 0) {
            tmp2 = *ptr2++;
            if (ptr2 - strip <= data)
                act2 = tmp2 - pre2;
            else
                act2 = tmp2 - pre2 - *(ptr2 - strip) + *(ptr2 - strip - 1);
            pre2 = tmp2;

            while (act2 < -128)
                act2 += 256;
            while (act2 > 127)
                act2 -= 256;

            writebits(DDS_file, act2 + (1 << bits2) / 2, bits2);
        }

        cnt2 = cnt1;
        bits2 = bits1;
    }

    if (cnt2 != 0) {
        writebits(DDS_file, cnt2, DDS_RL);
        writebits(DDS_file, DDS_code(bits2), 3);

        while (cnt2-- > 0) {
            tmp2 = *ptr2++;
            if (ptr2 - strip <= data)
                act2 = tmp2 - pre2;
            else
                act2 = tmp2 - pre2 - *(ptr2 - strip) + *(ptr2 - strip - 1);
            pre2 = tmp2;

            while (act2 < -128)
                act2 += 256;
            while (act2 > 127)
                act2 -= 256;

            writebits(DDS_file, act2 + (1 << bits2) / 2, bits2);
        }
    }

    flushbits(DDS_file);
    fclose(DDS_file);

    if (nofree == 0)
        free(data);
    else
        interleave(data, bytes, skip, DDS_INTERLEAVE);
}

// read a Differential Data Stream
unsigned char *readDDSfile(char *filename, voreen::ProgressBar* progress, unsigned int *bytes) {
    int version = 1;

    unsigned int skip, strip;

    unsigned char *data = 0, *ptr = 0;

    unsigned int cnt, cnt1, cnt2;
    int bits, act;

    if ((DDS_file = fopen(filename, "rb")) == NULL)
        return (NULL);

    for (cnt = 0; DDS_ID[cnt] != '\0'; cnt++)
        if (fgetc(DDS_file) != DDS_ID[cnt]) {
            fclose(DDS_file);
            version = 0;
            break;
        }

    if (version == 0) {
        if ((DDS_file = fopen(filename, "rb")) == NULL)
            return (NULL);

        for (cnt = 0; DDS_ID2[cnt] != '\0'; cnt++)
            if (fgetc(DDS_file) != DDS_ID2[cnt]) {
                fclose(DDS_file);
                return (NULL);
            }

        version = 2;
    }

    initbuffer();

    skip = readbits(DDS_file, 2) + 1;
    strip = readbits(DDS_file, 16) + 1;

    data = NULL;
    cnt = act = 0;

    while ((cnt1 = readbits(DDS_file, DDS_RL)) != 0) {
        if (progress)
            progress->setProgress(static_cast<float>(cnt) / static_cast<float>(DDS_BLOCKSIZE - 1));

        bits = DDS_decode(readbits(DDS_file, 3));

        for (cnt2 = 0; cnt2 < cnt1; cnt2++) {
            if (cnt <= strip)
                act += readbits(DDS_file, bits) - (1 << bits) / 2;
            else
                act += *(ptr - strip) - *(ptr - strip - 1) + readbits(DDS_file, bits) - (1 << bits) / 2;

            while (act < 0)
                act += 256;
            while (act > 255)
                act -= 256;

            if (cnt % DDS_BLOCKSIZE == 0) {
                if (data == NULL) {
                    if ((data = (unsigned char *)malloc(DDS_BLOCKSIZE)) == NULL)
                        ERRORMSG();
                } else
                    if ((data = (unsigned char *)realloc(data, cnt + DDS_BLOCKSIZE)) == NULL)
                        ERRORMSG();

                ptr = &data[cnt];
            }

            *ptr++ = act;
            cnt++;
        }
    }

    fclose(DDS_file);

    if (cnt == 0)
        return (NULL);

    if ((data = (unsigned char *)realloc(data, cnt)) == NULL)
        ERRORMSG();

    if (version == 1)
        interleave(data, cnt, skip);
    else
        interleave(data, cnt, skip, DDS_INTERLEAVE);

    *bytes = cnt;

    return (data);
}

// write a RAW file
void writeRAWfile(char *filename, unsigned char *data, unsigned int bytes, int nofree) {
    if (bytes < 1)
        ERRORMSG();

    if ((DDS_file = fopen(filename, "wb")) == NULL)
        ERRORMSG();
    if (fwrite(data, 1, bytes, DDS_file) != bytes)
        ERRORMSG();

    fclose(DDS_file);

    if (nofree == 0)
        free(data);
}

// read a RAW file
unsigned char *readRAWfile(char *filename, unsigned int *bytes) {
    unsigned char *data;
    unsigned int cnt, blkcnt;

    if ((DDS_file = fopen(filename, "rb")) == NULL)
        return (NULL);

    data = NULL;
    cnt = 0;

    do {
        if (data == NULL) {
            if ((data = (unsigned char *)malloc(DDS_BLOCKSIZE)) == NULL)
                ERRORMSG();
        } else
            if ((data = (unsigned char *)realloc(data, cnt + DDS_BLOCKSIZE)) == NULL)
                ERRORMSG();

        blkcnt = static_cast<unsigned int>(fread(&data[cnt], 1, DDS_BLOCKSIZE, DDS_file));
        cnt += blkcnt;
    } while (blkcnt == DDS_BLOCKSIZE);

    if (cnt == 0) {
        free(data);
        return (NULL);
    }

    if ((data = (unsigned char *)realloc(data, cnt)) == NULL)
        ERRORMSG();

    fclose(DDS_file);

    *bytes = cnt;

    return (data);
}

void swapshort(unsigned char *ptr, unsigned int size) {
    unsigned int i;
    unsigned char tmp;

    for (i = 0; i < size; i++) {
        tmp = ptr[0];
        ptr[0] = ptr[1];
        ptr[1] = tmp;
        ptr += 2;
    }
}

// write an optionally compressed PNM image
void writePNMimage(char *filename, unsigned char *image, unsigned int width, unsigned int height, unsigned int components, int dds) {
    char str[DDS_MAXSTR];

    unsigned char *data;

    if (width < 1 || height < 1)
        ERRORMSG();

    switch (components) {
    case 1:
        snprintf(str, DDS_MAXSTR, "P5\n%d %d\n255\n", width, height);
        break;
    case 2:
        snprintf(str, DDS_MAXSTR, "P5\n%d %d\n32767\n", width, height);
        break;
    case 3:
        snprintf(str, DDS_MAXSTR, "P6\n%d %d\n255\n", width, height);
        break;
    default:
        ERRORMSG();
    }

    if ((data = (unsigned char *)malloc(strlen(str) + width * height * components)) == NULL)
        ERRORMSG();

    memcpy(data, str, strlen(str));
    memcpy(data + strlen(str), image, width*height*components);

    if (DDS_ISINTEL)
        if (components == 2)
            swapshort(data + strlen(str), width*height);

    if (dds != 0)
        writeDDSfile(filename, data, static_cast<unsigned int>(strlen(str)) + width*height*components, components, width);
    else
        writeRAWfile(filename, data, static_cast<unsigned int>(strlen(str)) + width*height*components);
}

// read a possibly compressed PNM image
unsigned char *readPNMimage(char *filename, voreen::ProgressBar* progress, unsigned int *width, unsigned int *height, unsigned int *components) {
    const int maxstr = 100;

    char str[maxstr];

    unsigned char *data, *ptr1, *ptr2;
    unsigned int bytes;

    int pnmtype, maxval;
    unsigned char *image;

    if ((data = readDDSfile(filename, progress, &bytes)) == NULL)
        if ((data = readRAWfile(filename, &bytes)) == NULL)
            return (NULL);

    if (bytes < 4)
        return (NULL);

    memcpy(str, data, 3);
    str[3] = '\0';

    if (sscanf(str, "P%1d\n", &pnmtype) != 1)
        return (NULL);

    ptr1 = data + 3;
    while (*ptr1 == '\n' || *ptr1 == '#') {
        while (*ptr1 == '\n')
            if (++ptr1 >= data + bytes)
                ERRORMSG();
        while (*ptr1 == '#')
            if (++ptr1 >= data + bytes)
                ERRORMSG();
            else
                while (*ptr1 != '\n')
                    if (++ptr1 >= data + bytes)
                        ERRORMSG();
    }

    ptr2 = ptr1;
    while (*ptr2 != '\n' && *ptr2 != ' ')
        if (++ptr2 >= data + bytes)
            ERRORMSG();
    if (++ptr2 >= data + bytes)
        ERRORMSG();
    while (*ptr2 != '\n' && *ptr2 != ' ')
        if (++ptr2 >= data + bytes)
            ERRORMSG();
    if (++ptr2 >= data + bytes)
        ERRORMSG();
    while (*ptr2 != '\n' && *ptr2 != ' ')
        if (++ptr2 >= data + bytes)
            ERRORMSG();
    if (++ptr2 >= data + bytes)
        ERRORMSG();

    if (ptr2 - ptr1 >= maxstr)
        ERRORMSG();
    memcpy(str, ptr1, ptr2 - ptr1);
    str[ptr2 - ptr1] = '\0';

    if (sscanf(str, "%d %d\n%d\n", width, height, &maxval) != 3)
        ERRORMSG();

    if (*width < 1 || *height < 1)
        ERRORMSG();

    if (pnmtype == 5 && maxval == 255)
        *components = 1;
    else if (pnmtype == 5 && (maxval == 32767 || maxval == 65535))
        *components = 2;
    else if (pnmtype == 6 && maxval == 255)
        *components = 3;
    else
        ERRORMSG();

    if ((image = (unsigned char *)malloc((*width) * (*height) * (*components))) == NULL)
        ERRORMSG();
    if (data + bytes != ptr2 + (*width)*(*height)*(*components))
        ERRORMSG();

    if (DDS_ISINTEL)
        if (*components == 2)
            swapshort(ptr2, (*width)*(*height));

    memcpy(image, ptr2, (*width)*(*height)*(*components));
    free(data);

    return (image);
}

// write a compressed PVM volume
void writePVMvolume(char *filename, unsigned char *volume,
                    unsigned int width, unsigned int height, unsigned int depth, unsigned int components,
                    float scalex, float scaley, float scalez,
                    unsigned char *description,
                    unsigned char *courtesy,
                    unsigned char *parameter,
                    unsigned char *comment) {
    char str[DDS_MAXSTR];

    unsigned char *data;

    unsigned int len1 = 1, len2 = 1, len3 = 1, len4 = 1;

    if (width < 1 || height < 1 || depth < 1 || components < 1)
        ERRORMSG();

    if (description == NULL && courtesy == NULL && parameter == NULL && comment == NULL)
        if (scalex == 1.0f && scaley == 1.0f && scalez == 1.0f)
            snprintf(str, DDS_MAXSTR, "PVM\n%d %d %d\n%d\n", width, height, depth, components);
        else
            snprintf(str, DDS_MAXSTR, "PVM2\n%d %d %d\n%g %g %g\n%d\n", width, height, depth, scalex, scaley, scalez, components);
    else
        snprintf(str, DDS_MAXSTR, "PVM3\n%d %d %d\n%g %g %g\n%d\n", width, height, depth, scalex, scaley, scalez, components);

    if (description == NULL && courtesy == NULL && parameter == NULL && comment == NULL) {
        if ((data = (unsigned char *)malloc(strlen(str) + width * height * depth * components)) == NULL)
            ERRORMSG();

        memcpy(data, str, strlen(str));
        memcpy(data + strlen(str), volume, width*height*depth*components);

        writeDDSfile(filename, data, static_cast<unsigned int>(strlen(str)) + width*height*depth*components, components, width);
    } else {
        if (description != NULL)
            len1 = static_cast<unsigned int>(strlen((char *)description)) + 1;
        if (courtesy != NULL)
            len2 = static_cast<unsigned int>(strlen((char *)courtesy)) + 1;
        if (parameter != NULL)
            len3 = static_cast<unsigned int>(strlen((char *)parameter)) + 1;
        if (comment != NULL)
            len4 = static_cast<unsigned int>(strlen((char *)comment)) + 1;

        if ((data = (unsigned char *)malloc(strlen(str) + width * height * depth * components + len1 + len2 + len3 + len4)) == NULL)
            ERRORMSG();

        memcpy(data, str, strlen(str));
        memcpy(data + strlen(str), volume, width*height*depth*components);

        if (description == NULL)
            *(data + strlen(str) + width*height*depth*components) = '\0';
        else
            memcpy(data + strlen(str) + width*height*depth*components, description, len1);

        if (courtesy == NULL)
            *(data + strlen(str) + width*height*depth*components + len1) = '\0';
        else
            memcpy(data + strlen(str) + width*height*depth*components + len1, courtesy, len2);

        if (parameter == NULL)
            *(data + strlen(str) + width*height*depth*components + len1 + len2) = '\0';
        else
            memcpy(data + strlen(str) + width*height*depth*components + len1 + len2, parameter, len3);

        if (comment == NULL)
            *(data + strlen(str) + width*height*depth*components + len1 + len2 + len3) = '\0';
        else
            memcpy(data + strlen(str) + width*height*depth*components + len1 + len2 + len3, comment, len4);

        writeDDSfile(filename, data, static_cast<unsigned int>(strlen(str)) + width*height*depth*components + len1 + len2 + len3 + len4, components, width);
    }
}

// read a compressed PVM volume
unsigned char *readPVMvolume(char *filename, voreen::ProgressBar* progress,
                             unsigned int *width, unsigned int *height, unsigned int *depth, unsigned int *components,
                             float *scalex, float *scaley, float *scalez,
                             unsigned char **description,
                             unsigned char **courtesy,
                             unsigned char **parameter,
                             unsigned char **comment) {
    unsigned char *data, *ptr;
    unsigned int bytes, numc;

    int version = 1;

    unsigned char *volume;

    float sx = 1.0f, sy = 1.0f, sz = 1.0f;

    unsigned int len1 = 0, len2 = 0, len3 = 0, len4 = 0;

    if ((data = readDDSfile(filename, progress, &bytes)) == NULL)
        throw tgt::FileNotFoundException("Unable to open PVM file for reading", filename);
    if (bytes < 5)
        throw tgt::CorruptedFileException("PVM file corrupted", filename);

    if ((data = (unsigned char *)realloc(data, bytes + 1)) == NULL)
        ERRORMSG();
    data[bytes] = '\0';

    if (strncmp((char *)data, "PVM\n", 4) != 0) {
        if (strncmp((char *)data, "PVM2\n", 5) == 0)
            version = 2;
        else if (strncmp((char *)data, "PVM3\n", 5) == 0)
            version = 3;
        else
            throw tgt::CorruptedFileException("PVM file corrupted", filename);

#ifndef WIN32
        // need to switch to default locale "C" to prevent problems when current locale uses comma
        // instead of decimal point
        const char* oldlocale = setlocale(LC_NUMERIC, "C");
#endif

        if (sscanf((char *)&data[5], "%d %d %d\n%g %g %g\n", width, height, depth, &sx, &sy, &sz) != 6)
            ERRORMSG();

#ifndef WIN32
    setlocale(LC_NUMERIC, oldlocale);
#endif

        if (*width < 1 || *height < 1 || *depth < 1 || sx <= 0.0f || sy <= 0.0f || sz <= 0.0f)
            ERRORMSG();
        ptr = (unsigned char *)strchr((char *) & data[5], '\n') + 1;
    } else {
        if (sscanf((char *)&data[4], "%d %d %d\n", width, height, depth) != 3)
            ERRORMSG();
        if (*width < 1 || *height < 1 || *depth < 1)
            ERRORMSG();
        ptr = &data[4];
    }

    if (scalex != NULL && scaley != NULL && scalez != NULL) {
        *scalex = sx;
        *scaley = sy;
        *scalez = sz;
    }

    ptr = (unsigned char *)strchr((char *)ptr, '\n') + 1;
    if (sscanf((char *)ptr, "%d\n", &numc) != 1)
        ERRORMSG();
    if (numc < 1)
        ERRORMSG();

    if (components != NULL)
        *components = numc;
    else if (numc != 1)
        ERRORMSG();

    ptr = (unsigned char *)strchr((char *)ptr, '\n') + 1;
    if (version == 3)
        len1 = static_cast<unsigned int>(strlen((char *)(ptr + (*width) * (*height) * (*depth) * numc))) + 1;
    if (version == 3)
        len2 = static_cast<unsigned int>(strlen((char *)(ptr + (*width) * (*height) * (*depth) * numc + len1))) + 1;
    if (version == 3)
        len3 = static_cast<unsigned int>(strlen((char *)(ptr + (*width) * (*height) * (*depth) * numc + len1 + len2))) + 1;
    if (version == 3)
        len4 = static_cast<unsigned int>(strlen((char *)(ptr + (*width) * (*height) * (*depth) * numc + len1 + len2 + len3))) + 1;
    if ((volume = (unsigned char *)malloc((*width) * (*height) * (*depth) * numc + len1 + len2 + len3 + len4)) == NULL)
        ERRORMSG();
    if (data + bytes != ptr + (*width)*(*height)*(*depth)*numc + len1 + len2 + len3 + len4)
        ERRORMSG();

    memcpy(volume, ptr, (*width)*(*height)*(*depth)*numc + len1 + len2 + len3 + len4);
    free(data);

    if (description != NULL) {
        if (len1 > 1)
            *description = volume + (*width) * (*height) * (*depth) * numc;
        else
            *description = NULL;
    }

    if (courtesy != NULL) {
        if (len2 > 1)
            *courtesy = volume + (*width) * (*height) * (*depth) * numc + len1;
        else
            *courtesy = NULL;
    }

    if (parameter != NULL) {
        if (len3 > 1)
            *parameter = volume + (*width) * (*height) * (*depth) * numc + len1 + len2;
        else
            *parameter = NULL;
    }

    if (comment != NULL) {
        if (len4 > 1)
            *comment = volume + (*width) * (*height) * (*depth) * numc + len1 + len2 + len3;
        else
            *comment = NULL;
    }

    return (volume);
}

// simple checksum algorithm
unsigned int checksum(unsigned char *data, unsigned int bytes) {
    const unsigned int prime = 271;

    unsigned int i;

    unsigned char *ptr, value;

    unsigned int sum, cipher;

    sum = 0;
    cipher = 1;

    for (ptr = data, i = 0; i < bytes; i++) {
        value = *ptr++;
        cipher = prime * cipher + value;
        sum += cipher * value;
    }

    return (sum);
}

// swap the hi and lo byte of 16 bit data
void swapbytes(unsigned char *data, unsigned int bytes) {
    unsigned int i;
    unsigned char *ptr, tmp;

    for (ptr = data, i = 0; i < bytes / 2; i++, ptr += 2) {
        tmp = *ptr;
        *ptr = *(ptr + 1);
        *(ptr + 1) = tmp;
    }
}

// convert from signed short to unsigned short
void convbytes(unsigned char *data, unsigned int bytes) {
    unsigned int i;
    unsigned char *ptr;
    int v, vmin;

    for (vmin = 32767, ptr = data, i = 0; i < bytes / 2; i++, ptr += 2) {
        v = 256 * (*ptr) + *(ptr + 1);
        if (v > 32767)
            v = v - 65536;
        if (v < vmin)
            vmin = v;
    }

    for (ptr = data, i = 0; i < bytes / 2; i++, ptr += 2) {
        v = 256 * (*ptr) + *(ptr + 1);
        if (v > 32767)
            v = v - 65536;

        *ptr = (v - vmin) / 256;
        *(ptr + 1) = (v - vmin) % 256;
    }
}

// convert from float to unsigned short
void convfloat(unsigned char *data, unsigned int bytes) {
    unsigned int i;
    unsigned char *ptr;
    float v, vmax;

    for (vmax = 1.0f, ptr = data, i = 0; i < bytes / 4; i++, ptr += 4) {
        if (DDS_ISINTEL)
            DDS_swapuint((unsigned int *)ptr);

        v = (float) fabs(*((float *)ptr));
        if (v > vmax)
            vmax = v;
    }

    for (ptr = data, i = 0; i < bytes / 4; i++, ptr += 4) {
        v = (float) fabs(*((float *)ptr)) / vmax;

        data[2*i] = ftrc(65535.0f * v + 0.5f) / 256;
        data[2*i + 1] = ftrc(65535.0f * v + 0.5f) % 256;
    }
}

// helper functions for quantize:

inline int DDS_get(unsigned short int *volume,
                   unsigned int width, unsigned int height, unsigned int /*depth*/,
                   unsigned int i, unsigned int j, unsigned int k) {
    return (volume[i + (j + k*height)*width]);
}

inline double DDS_getgrad(unsigned short int *volume,
                          unsigned int width, unsigned int height, unsigned int depth,
                          unsigned int i, unsigned int j, unsigned int k) {
    double gx, gy, gz;

    if (i > 0)
        if (i < width - 1)
            gx = (DDS_get(volume, width, height, depth, i + 1, j, k) - DDS_get(volume, width, height, depth, i - 1, j, k)) / 2.0;
        else
            gx = DDS_get(volume, width, height, depth, i, j, k) - DDS_get(volume, width, height, depth, i - 1, j, k);
    else
        gx = DDS_get(volume, width, height, depth, i + 1, j, k) - DDS_get(volume, width, height, depth, i, j, k);

    if (j > 0)
        if (j < height - 1)
            gy = (DDS_get(volume, width, height, depth, i, j + 1, k) - DDS_get(volume, width, height, depth, i, j - 1, k)) / 2.0;
        else
            gy = DDS_get(volume, width, height, depth, i, j, k) - DDS_get(volume, width, height, depth, i, j - 1, k);
    else
        gy = DDS_get(volume, width, height, depth, i, j + 1, k) - DDS_get(volume, width, height, depth, i, j, k);

    if (k > 0)
        if (k < depth - 1)
            gz = (DDS_get(volume, width, height, depth, i, j, k + 1) - DDS_get(volume, width, height, depth, i, j, k - 1)) / 2.0;
        else
            gz = DDS_get(volume, width, height, depth, i, j, k) - DDS_get(volume, width, height, depth, i, j, k - 1);
    else
        gz = DDS_get(volume, width, height, depth, i, j, k + 1) - DDS_get(volume, width, height, depth, i, j, k);

    return (sqrt(gx*gx + gy*gy + gz*gz));
}

// quantize 16 bit volume to 8 bit using a non-linear mapping
unsigned char *quantize(unsigned char *volume,
                        unsigned int width, unsigned int height, unsigned int depth,
                        BOOLINT nofree,
                        BOOLINT linear,
                        BOOLINT verbose) {
    unsigned int i, j, k;

    unsigned char *volume2;
    unsigned short int *volume3;

    int v, vmin, vmax;

    double *err, eint;

    err = new double[65536];

    BOOLINT done;

    if ((volume3 = (unsigned short int*)malloc(width * height * depth * sizeof(unsigned short int))) == NULL)
        ERRORMSG();

    vmin = vmax = 256 * volume[0] + volume[1];

    for (k = 0; k < depth; k++)
        for (j = 0; j < height; j++)
            for (i = 0; i < width; i++) {
                v = 256 * volume[2 * (i + (j + k * height) * width)] + volume[2 * (i + (j + k * height) * width) + 1];
                volume3[i + (j + k*height)*width] = v;

                if (v < vmin)
                    vmin = v;
                else if (v > vmax)
                    vmax = v;
            }

    if (!nofree)
        free(volume);

    if (verbose)
        printf("16 bit volume has scalar range=[%d,%d]\n", vmin, vmax);

    if (linear)
        for (i = 0; i < 65536; i++)
            err[i] = 255 * (double)i / vmax;
    else {
        for (i = 0; i < 65536; i++)
            err[i] = 0.0;

        for (k = 0; k < depth; k++)
            for (j = 0; j < height; j++)
                for (i = 0; i < width; i++)
                    err[DDS_get(volume3, width, height, depth, i, j, k)] += sqrt(DDS_getgrad(volume3, width, height, depth, i, j, k));

        for (i = 0; i < 65536; i++)
            err[i] = pow(err[i], 1.0 / 3);

        err[vmin] = err[vmax] = 0.0;

        for (k = 0; k < 256; k++) {
            for (eint = 0.0, i = 0; i < 65536; i++)
                eint += err[i];

            done = TRUE;

            for (i = 0; i < 65536; i++)
                if (err[i] > eint / 256) {
                    err[i] = eint / 256;
                    done = FALSE;
                }

            if (done)
                break;
        }

        for (i = 1; i < 65536; i++)
            err[i] += err[i - 1];

        if (err[65535] > 0.0f)
            for (i = 0; i < 65536; i++)
                err[i] *= 255.0f / err[65535];
    }

    if ((volume2 = (unsigned char *)malloc(width * height * depth)) == NULL)
        ERRORMSG();

    for (k = 0; k < depth; k++)
        for (j = 0; j < height; j++)
            for (i = 0; i < width; i++)
                volume2[i + (j + k*height)*width] = (int)(err[DDS_get(volume3, width, height, depth, i, j, k)] + 0.5);

    delete[] err;
    free(volume3);

    if (verbose)
        printf("quantized volume to 8 bit\n");

    return (volume2);
}
