/*
* zi2spi.c : convert ZI scanner tif files to SPIDER format
* Copyright (C) 2006-2008  Health Research Inc.
* 
* HEALTH RESEARCH INCORPORATED (HRI),
* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455
* 
* Email:  spider@wadsworth.org
* 
* 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.
*/
/*
 * compile on SGI:
 * cc  -w -Xcpluscomm -o zi2spi.irix zi2spi.c -lm
 *
 * Linux:
 * gcc -w -DINTEL -o zi2spi.linux.XXX zi2spi.c -lm
 * where XXX is result of 'uname -m', e.g., x86_64
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/*  #include <conio.h> */

typedef unsigned char uint8;
typedef short int16;
typedef unsigned short uint16;
typedef unsigned long uint32;
typedef long int32;

#define BUFSIZE 20000

/* this macro is for Intergraph data, which is always little-endian,
even if the rest of the TIF header is bigendian   */
// #define ig2bytes(a,b) ( ((a) << 8) + (b) )


uint8 bigendian;

/* get: read n bytes from position pos */

int get(FILE *fd, int32 pos, uint8 *buf, int32 n) {
	int i;
	fseek(fd, pos, SEEK_SET);
	i = fread(buf, sizeof(uint8), n, fd);
	return(i);
}

uint16 twobytes(uint8 a, uint8 b) {
	if (bigendian) 	   return((a << 8) + b);
	else	   return((b << 8) + a);
}


uint32 fourbytes(uint8 a, uint8 b, uint8 c, uint8 d) {
	if (bigendian) return( (a << 24) + (b << 16) + (c << 8) + d );
	else return( (d << 24) + (c << 16) + (b << 8) + a );
}

uint16 ig2bytes(uint8 a, uint8 b) {
	return ((b << 8) + a);
}

void rawifd(uint8 *b) {
	printf("%x%x %x%x %x%x%x%x %x%x%x%x\n",
		b[0],b[1],b[2],b[3],b[4],b[5],b[6],b[7],b[8],b[9],b[10],b[11]);
}

void printifd(uint8 *b) {
	uint16 i;
	printf("%4d", twobytes(b[0],b[1]));
	printf(" %2d", twobytes(b[2],b[3]));
	printf(" %2x%2x%2x%2x", b[4],b[5],b[6],b[7]);
	printf(" %2x%2x%2x%2x", b[8],b[9],b[10],b[11]);
	printf("\n");
}

/*  TAGVALUE
	Returns value if count == 1;
	Returns offsets to first value if count > 1
*/
uint32 tagvalue(uint8 *b) {
  uint16  tag, type;
  uint32  count, value;

	type = twobytes(b[2],b[3]);
	count = fourbytes(b[4],b[5],b[6],b[7]);

	value = 0;
	if (count == 1) {
		if (type == 1) { value = b[8]; }
		else if (type == 3) { value = twobytes(b[8],b[9]); }
		else if (type == 4) { value = fourbytes(b[8],b[9],b[10],b[11]); }
	}
	else {
		value = fourbytes(b[8],b[9],b[10],b[11]);
	}

	return(value);
}

void transmiss(uint8 *b, FILE *fp, double *dtmin, double *dtmax) {
	uint16 tag;
	uint32 offset;
	int i;
	double *dbuf;
    uint8 btmp;

	union {
  	 double d;
  	 uint8 b[sizeof (double)];
	} un;


	tag = twobytes(b[0],b[1]);

	if (tag == 33918) {
		offset = fourbytes(b[8],b[9],b[10],b[11]);
		get(fp, offset, b, 8L);

		offset = offset + 8;

/*   get transmissivity min & max. They're doubles, so use union "un"  */
		dbuf = (double *) malloc(10 * sizeof(double));
		fseek(fp, offset+176, SEEK_SET);
		fread(dbuf, sizeof(double), 2, fp);
        un.d = dbuf[0];
/*
 the following works on Unix but not on Windows. Intergraph packet data does not necessarily follow the same bigendian/littlendian conventions as the rest of the tif header.  
If using on Linux, compile with  gcc -DINTEL
*/
#ifndef INTEL
	        if (!bigendian) {
              btmp = un.b[0]; un.b[0] = un.b[7]; un.b[7] = btmp;
              btmp = un.b[1]; un.b[1] = un.b[6]; un.b[6] = btmp;
              btmp = un.b[2]; un.b[2] = un.b[5]; un.b[5] = btmp;
              btmp = un.b[3]; un.b[3] = un.b[4]; un.b[4] = btmp;
            }
#endif
			
            *dtmin = un.d;
            un.d = dbuf[1];
			
#ifndef INTEL
	        if (!bigendian) {
              btmp = un.b[0]; un.b[0] = un.b[7]; un.b[7] = btmp;
              btmp = un.b[1]; un.b[1] = un.b[6]; un.b[6] = btmp;
              btmp = un.b[2]; un.b[2] = un.b[5]; un.b[5] = btmp;
              btmp = un.b[3]; un.b[3] = un.b[4]; un.b[4] = btmp;
            }
#endif
			
            *dtmax = un.d;
		
		free(dbuf);

	}   /* end if (tag == 33918) */

}
	
/* ***********************************************************************
*
* main
*
************************************************************************* */


void main (argc, argv)
int argc;
char *argv[];
{
 FILE *fp, *ofp;
 char *filename, *outfile, c;
 uint8  *buf, *b, *line, *tline, mybyt8a, mybyt8b;
 int entries, done, overview, fromget, indx;
 int labrec, lenbyt, labbyt, nrec, inb;
 uint16 tag, nbits, type, mybyt16, *line16;
 uint32 imsize, offset, n, firstoffset, xsize, ysize, linoff, fromget32;
 uint32 RowsPerStrip, tileByteOff, *TileByteCounts;
 uint32 TileWidth, TileLength, *TileOffsets, tileOff, tileCount, tilenum;
 uint32 nxtiles, nytiles, nxrem, nyrem, i, j, k, ti, x, nvalue, xrem, yrem;
 float *sphdr, *fline;
 double dtmax, dtmin, dm, db, dt;
 
 
   buf = (unsigned char *) malloc(BUFSIZE * sizeof(unsigned char));
   b = (unsigned char *) malloc(BUFSIZE * sizeof(unsigned char));
   filename = (char *) malloc(60 * sizeof(char));
   outfile = (char *) malloc(60 * sizeof(char));

   if (argc < 3) 
	 { printf("Usage: zi2spi infile outfile [overview#]\n"); exit(0); }
   else if (argc == 3) {
	 filename = argv[1];  
     outfile = argv[2];
     overview = 1;
   }
   else {
	 filename = argv[1];  
     outfile = argv[2];
     overview = atoi(argv[3]);
   }

   if ((fp = fopen(filename, "rb")) == NULL) {
	   printf("unable to open %s\n", filename); return; }
   if ((ofp = fopen(outfile, "wb")) == NULL) {
	   printf("unable to open %s\n", outfile); return; }


/* get the 1st 8 bytes: byte order, TIFF number, and first IFD */

	offset = 0;
	get(fp, offset, buf, 8L);

	if (buf[0] == 'M') { bigendian = 1; }
	else if (buf[0] == 'I') { bigendian = 0; }

/*	printf("byte order: %c %c\n", buf[0], buf[1]); */
	n = twobytes(buf[2], buf[3]);
/*    printf("TIFF number: %d\n", n); */

	offset = fourbytes(buf[4], buf[5], buf[6], buf[7]);
	firstoffset = offset;

	dtmin = dtmax = -1;

/* get transmissivity from the first Image File Header */

	get(fp, offset, buf, 2L);
	entries = twobytes(buf[0], buf[1]);

	offset += 2;
	for (i = 0; i < entries; i++) {
		get(fp, offset, buf, 12L);
	/*	 printifd(buf); */
		transmiss(buf, fp, &dtmin, &dtmax);
		offset += 12;
	}

/* follow the offset pointers in each IFD to get to the requested overview  */ 

	offset = firstoffset;

	for (k = 1; k < overview; k++) {

		get(fp, offset, buf, 2L);
		entries = twobytes(buf[0], buf[1]);
	/*	printf("%d entries:\n", entries);  */

		offset += 2;
		for (i = 0; i < entries; i++) {
			get(fp, offset, buf, 12L);
	/*		 printifd(buf); */
			transmiss(buf, fp, &dtmin, &dtmax);
			offset += 12;
		}

        get(fp, offset, buf, 4L);
	    offset = fourbytes(buf[0], buf[1], buf[2], buf[3]);
		if (offset == 0) {
			printf("data does not contain overview # %d\n", overview);
			exit(0);
		}

	}  /* end for k : offset points to overview */


/*	printf("transmissivity: lo %lf, hi %lf\n", dtmin, dtmax);     */
	if (dtmax == -1)
	 { printf("Error - no transmissivity values\n"); exit(1); }

	if (overview == 1) offset = firstoffset; 

    get(fp, offset, buf, 2L);
    entries = twobytes(buf[0], buf[1]);

	offset += 2;

/* get Xsize, Ysize, nbits */
	done = 0; 
	tileOff = 0;

	for (i = 0; i < entries; i++) {
		get(fp, offset, buf, 12L);
		tag = twobytes(buf[0],buf[1]);

		if (tag == 256) { xsize = tagvalue(buf); }
		if (tag == 257) { ysize = tagvalue(buf); }
		if (tag == 258) { nbits = tagvalue(buf); }
		if (tag == 322) { TileWidth = tagvalue(buf); }
		if (tag == 323) { TileLength = tagvalue(buf); }
		if (tag == 324) { 
			type = twobytes(buf[2],buf[3]); 
			tileCount = fourbytes(buf[4],buf[5],buf[6],buf[7]); 
			tileOff = fourbytes(buf[8],buf[9],buf[10],buf[11]); 
		}
		/*
		if (tag == 325) { 
			type = twobytes(buf[2],buf[3]); 
			tileCount = fourbytes(buf[4],buf[5],buf[6],buf[7]); 
			tileByteOff = fourbytes(buf[8],buf[9],buf[10],buf[11]); 
		}
		*/

		offset += 12;
	} 
/*
	printf("xsize : %d, ysize: %d, nbits: %d\n", xsize, ysize, nbits);
	printf("TileWidth : %d, TileLength: %d\n", TileWidth, TileLength);
	printf("tile offset: %d, tile count %d, type %d \n", tileOff, tileCount, type); 
*/
		if (xsize < TileWidth) {
		   printf("image size is less than tile size: cannot convert\n");
		 //  exit(1);
		}


	inb = 0;
	if (nbits == 8) inb = 1;
	if (nbits == 16) inb = 2;
	if (nbits == 12) inb = 2;
	if (inb == 0)
	    { printf("BitsPerSample is not a valid number: %d\n", nbits); exit(1); }

	//if (nbits == 8) dm = (dtmax - dtmin)/ 255.0;
	//if (nbits == 16) dm = (dtmax - dtmin)/ 4095.0;
	if (inb == 1) dm = (dtmax - dtmin)/ 255.0;
	if (inb == 2) dm = (dtmax - dtmin)/ 4095.0;
	db = dtmin;
//	printf("dtmin %lf, dtmax %lf, dm %lf, db %lf\n", dtmin, dtmax, dm, db);



	if (tileOff != 0) {
		if ( (TileOffsets = (uint32 *) malloc(tileCount * sizeof(uint32)) ) == NULL)
			{ printf("unable to allocate TileOffsets\n"); exit(1); }


		fromget = get(fp, tileOff, buf, tileCount*4L);
		for (i = 0; i < tileCount; i++) {
	 TileOffsets[i] = fourbytes(buf[4*i],buf[4*i+1],buf[4*i+2],buf[4*i+3]);
		}

		if (fromget != tileCount*4L)
			printf("get tile offsets: asked for %ld, read %ld\n", tileCount, fromget);

	}


/* ************* create spider header ************* */

	lenbyt = xsize * 4;
	labrec = 1024 / lenbyt;
	if (1024%lenbyt != 0) labrec++;
	labbyt = labrec * lenbyt;
	nrec = labbyt / 4;

	if ( (sphdr = (float *) malloc(nrec * sizeof(float)) ) == NULL)
		{ printf("unable to allocate spihdr\n"); exit(1); }
	for (i = 0; i < nrec; i++) sphdr[i] = 0;
/* note the C indices are 1 less than the actual (Fortran) header locations */
	sphdr[0] = 1;
	sphdr[1] = ysize;
	sphdr[4] = 1;
	sphdr[5] = 0;
	sphdr[11] = xsize;
	sphdr[12] = labrec;  /* there are labrec records in the header */
	sphdr[21] = labbyt;  /* total number of bytes in the header */
	sphdr[22] = lenbyt;  /* the record length in bytes */
	sphdr[nrec-1] = overview;   /* NOT a STANDARD SPIDER HEADER VALUE! */
/*
	printf("lenbyt %d, labrec %d, labbyt %d, nrec %d\n", lenbyt, labrec, labbyt, nrec); */


/* ************* write output ************* */

	fwrite(sphdr, sizeof(float), nrec, ofp);

	imsize = TileWidth * TileLength;
	nxtiles = xsize / TileWidth;
	
	xrem = xsize%TileWidth;
	yrem = ysize%TileLength;

	nytiles = ysize / TileLength;
	n = imsize * (nxtiles + 1);

/*	printf("xrem %ld, yrem %ld, nxtiles %ld, nytiles %ld\n", xrem, yrem, nxtiles, nytiles); 
*/

/* 8 bit data: inb = 1 ; 16 bit data: inb = 2  */

	if ( (tline = (uint8 *) malloc(TileWidth * inb * sizeof(uint8)) ) == NULL)
		{ printf("unable to allocate tline\n"); exit(1); }

	if (nbits == 8) {
	  if ( (line = (uint8 *) malloc(xsize * sizeof(uint8)) ) == NULL)
		{ printf("unable to allocate line\n"); exit(1); }
	}
	else if (nbits == 16 || nbits == 12) {
	  if ( (line16 = (uint16 *) malloc(xsize * sizeof(uint16)) ) == NULL)
		{ printf("unable to allocate line16\n"); exit(1); }
	}

	if ( (fline = (float *) malloc(xsize * sizeof(float)) ) == NULL)
		{ printf("unable to allocate fline\n"); exit(1); }

	x = TileWidth;


  for (ti = 0; ti < nytiles; ti++) {

	for (j = 0; j < TileLength; j++) {

		for (k = 0; k < nxtiles; k++) {

			if (xrem !=0) tilenum = ti * (nxtiles + 1L) + k;
			else tilenum = (ti * nxtiles) + k;

			offset = TileOffsets[tilenum] + (j * x*inb);

			fromget = get(fp, offset, tline, x*inb);

			if (fromget != x*inb) 
				printf("get: tile %ld, asked for %ld, read %d\n", tilenum, x, fromget);
			
			linoff = k * TileWidth;

			if (nbits == 8) {
			  for (i = 0; i < TileWidth; i++) 
			    { line[linoff + i] = tline[i]; }
			}
			else {
			  for (i = 0; i < TileWidth; i++)
			   { line16[linoff + i] = ig2bytes(tline[i*inb],tline[i*inb+1]); }
			}

		
		}   // for k

		if (xrem != 0) {
			tilenum = ti * (nxtiles + 1L) + nxtiles;
			offset = TileOffsets[tilenum] + (j * x*inb);
			fromget = get(fp, offset, tline, xrem*inb);
			/*
			if (fromget != xrem*inb) 
				printf("get: tile %ld, asked for %ld, read %d\n", tilenum, xrem*inb, fromget);
			*/
			linoff = nxtiles * TileWidth;

			if (nbits == 8) {
			  for (i = 0; i < xrem; i++) 
			    { line[linoff + i] = tline[i]; }
			}
			else {
			  for (i = 0; i < xrem; i++) 
			    { line16[linoff + i] = ig2bytes(tline[i*inb],tline[i*inb+1]); }
			}


		}  /* end xrem */

/*		convert to transmissivities, and then to densities  */

		nvalue = TileWidth*nxtiles+xrem;
		for (i = 0; i < nvalue; i++) {
			if (nbits == 8) dt = line[i]; 
			else dt = line16[i];

			dt = 5 + log10( (dm * dt) + db);
  		    fline[i] = dt;
		}

		fwrite(fline, sizeof(float), nvalue, ofp);

	
	}   // for j  

/*	printf("tilenum %ld offset %ld\n", tilenum, offset);   */

  } // for ti	

  if (yrem != 0) {
	for (j = 0; j < yrem; j++) {
		for (k = 0; k < nxtiles; k++) {
			if (xrem !=0)
			   offset = TileOffsets[nytiles*(nxtiles+1)+k] + (j * x*inb);
			else 
			   offset = TileOffsets[(nytiles*nxtiles)+k] + (j * x*inb);

			get(fp, offset, tline, x*inb);
			linoff = k * TileWidth;
			if (nbits == 8) {
			  for (i = 0; i < TileWidth; i++) 
			    { line[linoff + i] = tline[i]; }
			}
			else {
			  for (i = 0; i < TileWidth; i++) 
			    { line16[linoff + i] = ig2bytes(tline[i*inb],tline[i*inb+1]); }
			}


		}   // for k

		if (xrem != 0) {
			offset = TileOffsets[nytiles*(nxtiles+1) + nxtiles] + (j * x*inb);
			get(fp, offset, tline, xrem*inb);
			linoff = 1 * nxtiles * TileWidth;
			if (nbits == 8) {
			for (i = 0; i < xrem; i++) 
			  { line[linoff + i] = tline[i]; }
			}
			else {
			  for (i = 0; i < xrem; i++) 
			   { line16[linoff + i] = ig2bytes(tline[i*inb],tline[i*inb+1]); }
			}

		}

/*		fwrite(line, sizeof(uint8), TileWidth*nxtiles+xrem, ofp);   */

		nvalue = TileWidth*nxtiles+xrem;
		for (i = 0; i < nvalue; i++) {
		   if (nbits == 8) dt = line[i];
		   else dt = line16[i];
		   dt = 5 + log10( (dm * dt) + db);
		   fline[i] = dt;

		}

		fwrite(fline, sizeof(float), nvalue, ofp);

	} // for j
  }

	fclose(fp); fclose(ofp);
}