/*
 * $Id: readZimage.cxx,v 4.1 2004/05/25 06:40:20 hut66au Exp $
 *
 * Imview, the portable image analysis application
 * http://www.cmis.csiro.au/Hugues.Talbot/imview
 * ----------------------------------------------------------
 *
 *  Imview is an attempt to provide an image display application
 *  suitable for professional image analysis. It was started in
 *  1997 and is mostly the result of the efforts of Hugues Talbot,
 *  Image Analysis Project, CSIRO Mathematical and Information
 *  Sciences, with help from others (see the CREDITS files for
 *  more information)
 *
 *  Imview is Copyrighted (C) 1997-2001 by Hugues Talbot and was
 *  supported in parts by the Australian Commonwealth Science and 
 *  Industry Research Organisation. Please see the COPYRIGHT file 
 *  for full details. Imview also includes the contributions of 
 *  many others. Please see the CREDITS file for full details.
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *  
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *  
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA.
 * */

/*------------------------------------------------------------------------
 *
 * Z-image reader.
 *
 * Based on the Z-image reader for XV by the same author
 *
 * Hugues Talbot	14 Jan 1998
 *      
 *-----------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <errno.h>
#include "imview.hxx"
#include "imSystem.hxx"
#include "../imageIO.hxx"
#include "readZimage.hxx"
#include "../machine.hxx"

extern imageIO *IOBlackBox;



/* prototypes for static functions */
static void *readCHARZimage(FILE *fp, int nx, int ny, int nz);
static void *readINTZimage(FILE *fp, int nx, int ny, int nz, endianness e);
static void *readDOUBLEZimage(FILE *fp, int nx, int ny, int nz, endianness e);
static int needswapint(int *buf, long buflen);
static int swapint(int i);
static int needswapdouble(double *buf, long buflen);
static double swapdouble(double i);

// the `arithmetic exception' signal handler
static void arithExceptHandler(int sig);
// portable longjump buffer
im_jmp_buf AEenv;

int readZImage(const char *fname,
		int imageindex)
{
    int     res, start[3], end[3];
    int     thepixtype, theimgtype, spp, nbcomp;
    void   *buffp;
    
    // call the low-level Z-IMAGE reader
    res = load_zimage(fname, imageindex, start, end,
		      &thepixtype, &theimgtype, &spp, &nbcomp,
		      &buffp);

    if (res == 0) { // all went well
	// the new buffer needs to be set first before
	// the dimensions are changed because the previous buffer
	// might be freed, and we will need its dimensions.
	IOBlackBox->setCurrBuffp((void **)buffp);
	IOBlackBox->setCurrImgWidth(end[0]-start[0]+1);
	IOBlackBox->setCurrImgHeight(end[1]-start[1]+1);
	IOBlackBox->setCurrImgThickness(end[2]-start[2]+1);
	IOBlackBox->setXOffset(start[0]);
	IOBlackBox->setYOffset(start[1]);
	IOBlackBox->setZOffset(start[2]);
	IOBlackBox->setCurrZPos(0); // first plane
	IOBlackBox->setCurrImgType((imgtype)theimgtype);
	IOBlackBox->setCurrPixType((pixtype)thepixtype);
	IOBlackBox->setCurrImgNbComps(nbcomp); // Z-IMAGE supports multiple independent images in a single file
	IOBlackBox->setCurrImgNbSamples(spp);
	IOBlackBox->setImgDesc("Z-IMAGE");
    }

    return res;
}

int readZOverlay(const char *fname)
{
    void *p;
    int   res, start[3], end[3];
    int   thepixtype, theimgtype, imspp, nbcomp;

    // call the low-level Z-IMAGE reader
    res = load_zimage(fname, 0, start, end,
		      &thepixtype, &theimgtype, &imspp, &nbcomp,
		      &p);
    if (res == 0) {
	res = IOBlackBox->renderOverlay(fname, p, imspp, start, end, (pixtype)thepixtype);
	IOBlackBox->freeRawBuffers(p, imspp);
    }

    return res;
}


// This function returns the number of frames
// in a Z-IMAGE file. See description of load_zimage
// to understand what a `frame' is in the context
// of Z-IMAGE. This is in general NOT returning the number
// or components in the image, but rather the
// number of image `groups', or number of multicomponent sets.
// returns a negative number in case of error.
int zimagenbsubfiles(const char *fname)
{
    FILE *fp;
    int   nbread, Zversion, firstcomp;
    int   *firstR, *lastR, *firstC, *lastC, *firstP, *lastP, *imKind;
    int    i, n, nbim;
    long   filesize;
    char   alias[12], kind[10];

    
    fp = im_fopen(fname, "rb");

    if (!fp) {
	dbgprintf("zimagenbsubfiles: cannot read file %s\n", fname);
	return -1;
    }
    
    /* compute file length */
    fseek(fp, 0L, SEEK_END);
    filesize = ftell(fp);
    fseek(fp, 0L, SEEK_SET);

    /* read the Header now */
    nbread = fscanf(fp, "CSIROZ version %d\n", &Zversion);
    if (nbread != 1) {
	errprintf("zimagenbsubfiles: faulty Z-IMAGE format %s\n", fname);
	return(-2);
    }
    /* read the number of images in the file */
    if ((nbread = fscanf(fp, "%d\n", &nbim)) != 1) {
	errprintf("zimagenbsubfiles: can't read nb of components in %s\n", fname);
	return(3);
    }

    /* now read all the data in the headers */
    firstR = new int[nbim];
    lastR =  new int[nbim];
    firstC = new int[nbim];
    lastC =  new int[nbim];
    firstP = new int[nbim];
    lastP =  new int[nbim];
    
    imKind =  new int[nbim];

    /* read component images */
    for (i = 0 ; i < nbim ; i++) {
	if ((nbread = fscanf(fp, "%d %s %d %d %d %d ",
			     &n, alias, firstR+i, lastR+i, firstC+i, lastC+i
	    )) != 6) {
	    delete[] firstR;
	    delete[] lastR;
	    delete[] firstC;
	    delete[] lastC;
	    delete[] firstP;
	    delete[] lastP;
	    delete[] imKind;
	    errprintf("zimagenbsubfiles: error in first part of image header for component %d, file %s\n", i,fname);
	    return(-4);
	} else {
	    /* now try to decide whether this is a 3D image */
	    if ((nbread = fscanf(fp, "%d %d %s\n", firstP+i, lastP+i, kind)) != 3)  {
		/* this is not a 3D image, try again */
		firstP[i] = 0;
		lastP[i] = 0;
		if ((nbread = fscanf(fp, "%s\n", kind)) != 1) {
		    /* now this is an error */
		    delete[] firstR;
		    delete[] lastR;
		    delete[] firstC;
		    delete[] lastC;
		    delete[] firstP;
		    delete[] lastP;
		    delete[] imKind;
		    errprintf("zimagenbsubfiles: error in second part of image header for component %d, file %s\n", i, fname);
		    return(-5);
		}
	    }	
	    if (strcmp(kind, "BINARY") == 0)
		imKind[i] = IM_BINARY;
	    else if ((strcmp(kind, "CHAR") == 0) || (strcmp(kind, "GREY") == 0))
		imKind[i] = IM_UINT1; // UNSIGNED char is the default!
	    else if (strcmp(kind, "INT") == 0)
		imKind[i] = IM_INT;
	    else if (strcmp(kind, "DOUBLE") == 0)
		imKind[i] = IM_DOUBLE;
	    else {
		delete[] firstR;
		delete[] lastR;
		delete[] firstC;
		delete[] lastC;
		delete[] firstP;
		delete[] lastP;
		delete[] imKind;
		errprintf("zimagenbsubfiles: broken header for component %d, file %s\n", i, fname);
		return(-6);
	    }
	}
    }


    int tmpindex = 0, imageindex = 0;
    i = 1;
    do {
	firstcomp = tmpindex;
	while (i < nbim) { 
	    if (imKind[tmpindex] != imKind[i])
		break;
	    if ((lastR[i] != lastR[tmpindex])
		|| (lastC[i] != lastC[tmpindex])
		|| (lastP[i] != lastP[tmpindex]))
		break;
	    if ((firstR[i] != firstR[tmpindex])
		|| (firstC[i] != firstC[tmpindex])
		|| (firstP[i] != firstP[tmpindex]))
		break;
	    i++;
	}
	tmpindex = i;
	imageindex++;
    } while (i < nbim);

    // cleanup
    fclose(fp);
    
    delete[] firstR;
    delete[] lastR;
    delete[] firstC;
    delete[] lastC;
    delete[] firstP;
    delete[] lastP;
    delete[] imKind;
    
    return imageindex;
}

int load_zimage(const char *fname,	 /* input file name  */
		int         imageindex,	 /* index in the image file */
		int         start[3],	 /* coordinates of the starting point */
		int         end[3],	 /* coordinates of the end point */
		int        *thepixtype,	 /* pixel type (BINARY, CHAR, etc) */
		int        *theimgtype,	 /* image type (MULTI, RGB, etc) */
		int        *spp,         /* samples per pixels (read number of components) */
		int        *nbcomp,      /* total number of components in file, including unread ones */
		void      **inbuffp)	 /* pixel data (will be allocated) */
{
/** Reads a multispectral Z-IMAGE file format

    RETURN VALUE:	int 

    DESCRIPTION:

    This function reads a Z-IMAGE file and returns a multispectral
    image in a series of buffers.  The difficulties with Z-IMAGE is
    that the format does not indicate either the bit ordering or the
    image interpretation (RGB, etc). What this function does is try to
    work out both of these with the following strategy (heuristics):
    1) guess the bit ordering for INT and DOUBLE by finding out which
    is yielding the smallest range
    2) guess the photometric interpretation by reading the header
    information and collating together all the sub-images which have
    the same x,y and z dimensions, and the same offsets. If there are
    3 such images, the image is deemed RGB, if there are more than one
    (but not 3), it is deemed a MULTISPECTRAL image, if there is only
    one such image, it is deemed a SINGLE component image.

    If there are more components after the last one that has been
    interpreted as being part of a multispectral image, they are just
    ignored.

    It is possible to read all the sub-images in a file by specifying
    an image index `n' greater than 0. The caller will get the nth
    multispectral set.

    HISTORY:
    Written by Hugues Talbot	19 Jan 1998

    TESTS:

    REFERENCES:

    KEYWORDS:

**/
    FILE *fp;
    char  alias[12], kind[10], marker;
    int   res = 0;
    int   i, nbim, n, firstcomp, lastcomp, nbread, Zversion;
    long  filesize;
    int	 *firstR, *lastR, *firstC, *lastC, *firstP, *lastP;
    int  *imKind;
    void **buffp = 0;
    int   nx, ny, nz;
    long offs = 0;
    endianness e;

    
    fp = im_fopen(fname, "rb");
    if (!fp) {
	dbgprintf("load_zimage: cannot read file %s\n", fname);
	return 1;
    }
    
    /* compute file length */
    fseek(fp, 0L, SEEK_END);
    filesize = ftell(fp);
    fseek(fp, 0L, SEEK_SET);

    /* read the Header now */
    nbread = fscanf(fp, "CSIROZ version %d\n", &Zversion);

    dbgprintf("CSIRO Z version %d\n", Zversion);

    if (nbread != 1) {
	errprintf("load_zimage: faulty Z-IMAGE format %s\n", fname);
	return(2);
    }
    /* read the number of images in the file */
    if ((nbread = fscanf(fp, "%d\n", &nbim)) != 1) {
	errprintf("load_zimage: can't read nb of components in %s\n", fname);
	return(3);
    }

    /* now read all the data in the headers */
    firstR = new int[nbim];
    lastR =  new int[nbim];
    firstC = new int[nbim];
    lastC =  new int[nbim];
    firstP = new int[nbim];
    lastP =  new int[nbim];
    
    imKind =  new int[nbim];

    /* read component images */
    for (i = 0 ; i < nbim ; i++) {
	if ((nbread = fscanf(fp, "%d %s %d %d %d %d ",
			     &n, alias, firstR+i, lastR+i, firstC+i, lastC+i
	    )) != 6) {
	    delete[] firstR;
	    delete[] lastR;
	    delete[] firstC;
	    delete[] lastC;
	    delete[] firstP;
	    delete[] lastP;
	    delete[] imKind;
	    errprintf("load_zimage: error in first part of image header for component %d, file %s\n", i,fname);
	    return(4);
	} else {
	    /* now try to decide whether this is a 3D image */
	    if ((nbread = fscanf(fp, "%d %d %s\n", firstP+i, lastP+i, kind)) != 3)  {
		/* this is not a 3D image, try again */
		firstP[i] = 0;
		lastP[i] = 0;
		if ((nbread = fscanf(fp, "%s\n", kind)) != 1) {
		    /* now this is an error */
		    delete[] firstR;
		    delete[] lastR;
		    delete[] firstC;
		    delete[] lastC;
		    delete[] firstP;
		    delete[] lastP;
		    delete[] imKind;
		    errprintf("load_zimage: error in second part of image header for component %d, file %s\n", i, fname);
		    return(5);
		}
	    }	
	    if (strcmp(kind, "BINARY") == 0)
		imKind[i] = IM_BINARY;
	    else if ((strcmp(kind, "CHAR") == 0) || (strcmp(kind, "GREY") == 0))
		imKind[i] = IM_UINT1; // Unsigned char is the default!
	    else if (strcmp(kind, "INT") == 0)
		imKind[i] = IM_INT;
	    else if (strcmp(kind, "DOUBLE") == 0)
		imKind[i] = IM_DOUBLE;
	    else {
		delete[] firstR;
		delete[] lastR;
		delete[] firstC;
		delete[] lastC;
		delete[] firstP;
		delete[] lastP;
		delete[] imKind;
		errprintf("load_zimage: broken header for component %d, file %s\n", i, fname);
		return(6);
	    }
	}
    }

    // endianness issue
    if (Zversion > 5) {
        if (strncmp(alias, "\"(null)\"", 8) == 0) {
            dbgprintf("File advertises little-endian data\n");
            e = ENDIAN_LITTLE;
        } else {
            dbgprintf("File advertises big-endian data\n");
            e = ENDIAN_BIG;
        }
    } else {
        dbgprintf("Z version too old, unknown data endianness\n");
        e = ENDIAN_UNKNOWN;
    }


    /* now work out what to do */
    if (imageindex >= nbim) {
	dbgprintf("load_zimage: required image index greater than number of components: %d >= %d\n", imageindex, nbim);
	imageindex = 0;
    }
    
    int tmpindex = 0;
    i = 1;
    *nbcomp = 0;
    do {
	firstcomp = tmpindex;
	while (i < nbim) { 
	    if (imKind[tmpindex] != imKind[i])
		break;
	    if ((lastR[i] != lastR[tmpindex])
		|| (lastC[i] != lastC[tmpindex])
		|| (lastP[i] != lastP[tmpindex]))
		break;
	    if ((firstR[i] != firstR[tmpindex])
		|| (firstC[i] != firstC[tmpindex])
		|| (firstP[i] != firstP[tmpindex]))
		break;
	    i++;
	}
	(*nbcomp)++;
	tmpindex = i;
    } while ((imageindex-- > 0) && (i < nbim));

    if ((imageindex != -1) && (i==nbim)) {
	// caller asked for a component beyond the end of the file
	// this is legitimate! just an error code reported
	// cleanup
	fclose(fp);

	delete[] firstR;
	delete[] lastR;
	delete[] firstC;
	delete[] lastC;
	delete[] firstP;
	delete[] lastP;
	delete[] imKind;

	*inbuffp = 0;
	return(7);
    }
    
    lastcomp = i-1;
    *spp = lastcomp-firstcomp+1; // this is the real number of samples per pixel
    *thepixtype = imKind[firstcomp];

    if (*spp == 3)
	*theimgtype = IM_RGB;       // could be IM_HLS, but there is no way of knowing
    else 
	*theimgtype = IM_SPECTRUM;  // this includes single-component images

    start[1] = firstR[firstcomp];
    end[1] = lastR[firstcomp];
    start[0] = firstC[firstcomp];
    end[0] = lastC[firstcomp];
    start[2] = firstP[firstcomp];
    end[2] = lastP[firstcomp];

    nx = end[0] - start[0] +1;
    ny = end[1] - start[1] +1;
    nz = end[2] - start[2] +1;

    // start reading the data
    
    // we must skip over the nth first components,
    // which are of no interest to us...

    dbgprintf("Position in Z file: %ld\n", ftell(fp));
    // apparently Win32 eats the 0 by itself. This is weird
    // but first we have to jump over this nasty '\0'
    if ((marker = skipzero(fp)) != '\0') {
	errprintf("load_zimage: wrong first data byte in file: got %d\n", marker);
	return 8;
    }
    dbgprintf("Position in Z file after eating '\\0': %ld\n", ftell(fp));
    
    for (i = 0 ; i < firstcomp ; i++) {
	int nnx = lastR[i] - firstR[i] + 1;
	int nny = lastC[i] - firstC[i] + 1;
	int nnz = lastP[i] - firstP[i] + 1;
	switch(imKind[i]) {
	  case IM_BINARY:
	  case IM_CHAR:
	  case IM_UINT1:
	    offs = nnx*nny*nnz*sizeof(char);
	    break;
	    
	  case IM_UINT4:
	  case IM_INT4:
	    offs = nnx*nny*nnz*sizeof(int);
	    break;
	    
	  case IM_DOUBLE:
	    offs = nnx*nny*nnz*sizeof(double);
	    break;

	  default:
	    offs = 0;
	    break;
	}
	if ((res = fseek(fp, offs, SEEK_CUR)) != 0) {
	    errprintf("Can't jump over unwanted components\n%s",
		      strerror(errno));
	    // cleanup
	    fclose(fp);
	    delete[] firstR;
	    delete[] lastR;
	    delete[] firstC;
	    delete[] lastC;
	    delete[] firstP;
	    delete[] lastP;
	    delete[] imKind;
	    *inbuffp =  0;
	    return(9);
	}
    }

    dbgprintf("Position in Z file after skipping unwanted components: %ld\n",
	      ftell(fp));
    dbgprintf("Should have skipped %ld bytes\n", offs);
    
    buffp = (void **) malloc(*spp * sizeof(void **));
    
    for (i = firstcomp ; i <= lastcomp ; i++) {
	switch (*thepixtype) {
	  case IM_BINARY:
	  case IM_UINT1:
	  case IM_INT1:
	    buffp[i-firstcomp] = readCHARZimage(fp, nx, ny, nz);
	    break;

	  case IM_UINT4:
	  case IM_INT4:
	    buffp[i-firstcomp] = readINTZimage(fp, nx, ny, nz, e);
	    break;

	  case IM_DOUBLE:
	    buffp[i-firstcomp] = readDOUBLEZimage(fp, nx, ny, nz, e);
	    break;
	    
	  default:
	    errprintf("load_zimage: type %d reader not implemented yet\n", *thepixtype);
	    break;
	}
	if (buffp[i-firstcomp] == 0) {
	    // an error occured
	    for (int j = firstcomp ; j < i ; j++)
		free(buffp[i-firstcomp]);
	    free(buffp);
	    buffp = 0;
	    res = 10;
	}
    }

    // cleanup
    fclose(fp);

    delete[] firstR;
    delete[] lastR;
    delete[] firstC;
    delete[] lastC;
    delete[] firstP;
    delete[] lastP;
    delete[] imKind;

    *inbuffp = buffp;
    return res;
}

endianness checkendian(void)
{
    /* Are we little or big endian?  From Harbison&Steele.  */
    union
    {
        long l;
        char c[sizeof (long)];
    } u;
    u.l = 1;
    if (u.c[sizeof (long) - 1] == 1)
        return(ENDIAN_BIG);
    else
        return(ENDIAN_LITTLE);
}

static void *readCHARZimage(FILE *fp, int nx, int ny, int nz)
{
    uchar *retbuf;
    long imlength, nbread;

    imlength = nx*ny*nz;

    retbuf = (uchar *)malloc(nx*ny*nz * sizeof(uchar)); // stupid malloc, to be compatible with the TIFF reader

    if (retbuf == 0) {
	errprintf("readCHARZimage: not enough memory to load file\n");
	return 0;
    }

    nbread = fread(retbuf, sizeof(uchar), imlength, fp);
    if (nbread != imlength) {
	errprintf("load_zimage: incomplete image\n");
	// this is only a warning, complete with zeroes
	memset(retbuf+nbread, 0, (imlength-nbread)*sizeof(uchar));
    }

    return retbuf;
}



// this is not as easy because of the big-endian vs. little-endian stuff
static void *readINTZimage(FILE *fp, int nx, int ny, int nz, endianness e)
{
    int  *retbuf;
    long  imlength, nbread;

    imlength = nx*ny*nz;

    retbuf = (int *)malloc(nx*ny*nz * sizeof(int)); // stupid malloc, to be compatible with the TIFF reader

    if (retbuf == 0) {
	errprintf("readINTZimage: not enough memory to load file\n");
	return 0;
    }

    nbread = fread(retbuf, sizeof(int), imlength, fp);
    if (nbread != imlength) {
	errprintf("load_zimage: incomplete image\n");
	// this is only a warning, complete with zeroes
	memset(retbuf+nbread, 0, (imlength-nbread)*sizeof(int)); 
    }

    if (((e != ENDIAN_UNKNOWN) && (e != checkendian()))
        || ((e == ENDIAN_UNKNOWN) && (needswapint(retbuf, imlength)))) {
	dbgprintf("Image buffer will be swapped\n");
	swapintbuf(retbuf, imlength); // in-place processing 
    } else
	dbgprintf("Image buffer will not be swapped\n");
    
    return retbuf;
}

// this is not as easy because of the big-endian vs. little-endian stuff
static void *readDOUBLEZimage(FILE *fp,int nx, int ny, int nz, endianness e)
{
    double *retbuf;
    int  ns;
    long imlength, nbread;

    imlength = nx*ny*nz;

    retbuf = (double *)malloc(nx*ny*nz * sizeof(double)); // stupid malloc, to be compatible with the TIFF reader

    if (retbuf == 0) {
	errprintf("readDOUBLEZimage: not enough memory to load file\n");
	return 0;
    }

    nbread = fread(retbuf, sizeof(double), imlength, fp);
    if (nbread != imlength) {
	errprintf("readDOUBLEZimage: incomplete image\n");
	// this is only a warning, complete with zeroes
	memset(retbuf+nbread, 0, (imlength-nbread)*sizeof(double));
    }

    ns = 2;
    if (((e != ENDIAN_UNKNOWN) && (e != checkendian()))
        || ((e == ENDIAN_UNKNOWN) && ((ns = needswapdouble(retbuf, imlength)) == 1))) {
	dbgprintf("Image buffer will be swapped\n");
	swapdoublebuf(retbuf, imlength); // in-place processing 
    } else if (ns == -1) {
        /* if needswapdouble was called but caused SIGFPE on both 
           swapped and unswapped versions */
        errprintf("readDoubleZImage: cannot make sense of image\n because of invalid double precision values present\n");
	free(retbuf);
	retbuf = 0;
    } else {
	dbgprintf("Image buffer will not be swapped\n");
    }

    return retbuf;
}



static int needswapint(int *buf, long buflen)
{
    int *p, swp, *end, minval, maxval;
    int swminval, swmaxval;
    unsigned int range, swrange;
    int sample = trivmin((buflen/503 + 1), 997); // largish prime number

    
    p = buf;
    end = buf+buflen;
    minval = maxval = *p;

    // try the unswapped version
    while (p < end) {
	if (*p > maxval)
	    maxval = *p;
	if (*p < minval)
	    minval = *p;
	p += sample;    // no need to sample all the pixels...
    }
    range = maxval - minval;
    dbgprintf("needswapint: unswapped range = %u\n", range);
	    
    // now do the same thing with swap on
    p = buf;
    swmaxval = swminval = swapint(*p);
    while (p < end) {
	swp = swapint(*p); 
	if (swp > swmaxval)
	    swmaxval = swp;
	if (swp < swminval)
	    swminval = swp;
	p += sample;    // no need to sample all the pixels...
    }
    swrange = swmaxval - swminval;
    dbgprintf("needswapint: swapped range = %u\n", swrange);

    if ((swrange == 0) && (range == 0)) 
	return 0; // we don't know so we don't swap
    else if (swrange < range)
	return 1; // yes, need swapping!
    else if ((swrange == range) && (swmaxval < maxval))
	return 1 ; // yes: still need swapping
    else return 0; // no need to swap
}

static void arithExceptHandler(int sig)
{
    // reset the signal handler
    signal(sig, arithExceptHandler);
    im_longjmp(AEenv, 1);
}

// log are used in there otherwise FP exception
// might be raised when comparing extreme numbers
// that would only be possible with loss of precision
// I.e: comparing 10^-300 and 10^+300

typedef void (signalHandler)(int);

static int needswapdouble(double *buf, long buflen)
{
    volatile void *vp;
    double *p, lp, swp, lswp, *end, minval, maxval; vp = &p; vp = &lp; vp = &swp; vp = &lswp; vp= &end; vp = &minval; vp = &maxval;
    double swminval, swmaxval; vp = &swminval; vp = &swmaxval;
    double range, swrange; vp = &range ; vp = &swrange;
    int sample = trivmin((buflen/503 + 1),997); vp = &sample;// largish prime number 
    bool raisedWhenSwapped, raisedWhenNotSwapped; vp = &raisedWhenSwapped; vp = &raisedWhenNotSwapped;
    int  retval = 0; vp = &retval;
    
    p = buf;
    end = buf+buflen;

    signalHandler *oldHandler = signal(SIGFPE, arithExceptHandler);
    if (oldHandler == SIG_ERR) {
	errprintf("needswapdouble: Could not set up FP exception handler.\n");
	// continue anyway!
    }

    if (im_setjmp(AEenv, 1) == 0) {
	// try the unswapped version
	minval = maxval = 0;
	// look for first non-null pixel, if any 
	while (p < end) {
	    if (*p != 0.0) {
		lp = log(fabs(*p));
		if ((lp > 1000) || (lp < -1000)) {
		    p += sample;
		    continue;
		}
		/* first value found */
		minval = maxval = *p;
		break;
	    }
            ++p;
	}
	while (p < end) {
	    //skip the zeroes
	    if (*p != 0.0) {
		lp = log(fabs(*p));
		if ((lp > 1000) || (lp < -1000)) {
		    p += sample;
		    continue;
		}
	        
		if (lp > maxval)
		    maxval = lp;
		if (lp < minval)
		    minval = lp;
	    }
	    p += sample;    // no need to sample all the pixels...
	}
	range = maxval - minval;
	dbgprintf("needswapdouble: unswapped log range = %g\n", range);
	raisedWhenNotSwapped = false;
    } else {
	// we got there because of an exception:
	raisedWhenNotSwapped = true;
    }
    
    // now do the same thing with swap on
    p = buf;
    swmaxval = swminval = 0;
    if (im_setjmp(AEenv, 1) == 0) {
	// look for first non-null pixel, if any 
	while (p < end) {
	    swp = swapdouble(*p);
	    if (swp != 0.0) {
		lswp = log(fabs(swp));
		if ((lswp > 1000) || (lswp < -1000)) {
		    p += sample;
		    continue;
		}
		/* first value found */
		minval = maxval = lswp;
		break;
	    }
            ++p;
	}
	while (p < end) {
	    swp = swapdouble(*p);
	    // skip over the zeroes
	    if (swp != 0.0) {
		lswp = log(fabs(swp));
		if ((lswp > 1000) || (lswp < -1000)) {
		    p += sample;
		    continue;
		}
		if (lswp > swmaxval)
		    swmaxval = lswp;
		if (lswp < swminval)
		    swminval = lswp;
	    }
	    p += sample;    // no need to sample all the pixels...
	}
	swrange = swmaxval - swminval;
	dbgprintf("needswapdouble: swapped log range = %g\n", swrange);
	raisedWhenSwapped = false;
    } else {
	raisedWhenSwapped = true;
    }

    // reset the FPE handler to its previous state
    if (oldHandler != SIG_ERR) {
	if (signal(SIGFPE, oldHandler) == SIG_ERR) {
	    errprintf("needswapdouble: Could not restore original FP exception handler\n");
	    //continue anyway!
	}
	
    }

    if (raisedWhenSwapped && raisedWhenNotSwapped)
	retval = -1; // this image is basically unreadable
    else if (raisedWhenSwapped)
	retval = 0; // do not swap
    else if (raisedWhenNotSwapped)
	retval = 1; // do swap, then
    else  if (swrange < range)
	retval = 1; // yes, need swapping!
    else if ((swrange == 0) && (range == 0))
	retval = 0; // we have a constant value in the image, we don't swap because it could be anything
    else if ((swrange == range) && (swmaxval < maxval))
	retval = 1 ; // yes: still need swapping
    else retval = 0; // no need to swap

    return retval;
}

static int swapint(int i)
{
    uchar *from, *to;
    int    j;
    
    from = (uchar *)&i;
    to = (uchar *)&j;

    to[0] = from[3];
    to[1] = from[2];
    to[2] = from[1];
    to[3] = from[0];
    
//  j = ((int)(p[0]) << 24) +
//	((int)(p[1]) << 16) +
//	((int)(p[2]) << 8)  +
//	((int)(p[3]));

    return j;
}

static double swapdouble(double d)
{
    uchar *from, *to;
    double b;

    from = (uchar *)&d;
    to   = (uchar *)&b;
    to[0] = from[7];
    to[1] = from[6];
    to[2] = from[5];
    to[3] = from[4];
    to[4] = from[3];
    to[5] = from[2];
    to[6] = from[1];
    to[7] = from[0];
    
    return b;
}


// no questions asked...
// in-place processing
void swapintbuf(int *buf, long buflen)
{
    int *p, *end;

    dbgprintf("Swapping int buffer\n");

    p = buf;
    end = p + buflen;

    while (p < end) {
	*p = swapint(*p);
	p++;
    }

    return;
}

void swapdoublebuf(double *buf, long buflen)
{
    double *p, *end;

    dbgprintf("Swapping double buffer\n");

    p = buf;
    end = p + buflen;

    while (p < end) {
	*p = swapdouble(*p);
	p++;
    }

    return;
}