#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* compile this file with cc -o adsc2sgi adsc2sgi.c -lm */

#define TRUE 1
#define FALSE 0

/* global notifier of frame byte=order */
int swap_bytes = FALSE;

/* frame handling routines */
FILE *GetFrame(char *filename);

typedef struct _SGIHeader
{
  short Magic;          /* Identification number (474) */
  char storage;         /* Compression flag */
  char bpc;             /* Bytes per pixel */
  unsigned short int dimension;       /* Number of image dimensions */
  unsigned short int xsize;           /* Width of image in pixels */
  unsigned short int ysize;           /* Height of image in pixels */
  unsigned short int zsize;           /* Number of bit planes */
  long pixmin;          /* Smallest pixel value */
  long pixmax;          /* Largest pixel value */
  char dummy1[4];       /* Not used */
  char name[80];	/* Name of image */
  long colormap;        /* Format of pixel data */
  char dummy2[404];     /* Not used */
} SGIHEAD;



main(int argv, char** argc)
{
    /* exit status */
    int return_code = 0;

    /* utility variables */
    int i;
    char *filename;
    FILE *frame = NULL;
    unsigned short int **image;

    /* frame information variables */
    int record_size;
    int num_records;
    int pixels_per_line;
    int image_width;
    int image_height;
    float eight;

    /* indicies */
    int x, z;
    int overloads = 0;
    
    /* limits */
    int high = 0;
    int low = 0xFFFF;
    double sum = 0;
    double sigma = 0;
    double scale = 0;
    double zoom = 1;
    int negate  = 0;
    int start_x = 0;
    int start_y = 0;
    int  stop_x = 0;
    int  stop_y = 0;
    
    /* storage */
    unsigned short int *record;
    unsigned short int *swapped;
    unsigned short int *temp;

    
    /* output RGB stuff */
    SGIHEAD        SGIhead;
    unsigned char *SGIdata;
    FILE *rgbout = NULL;
    char outfilename[256] = "new.rgb";

    /* begin progam here */
    printf("ADSC2SGI 1.0.1\t\tby James Holton 3-4-99\n\n");

    /* check argument list */
    for(i=1; i<argv; ++i)
    {
	if(strstr(argc[i], ".img"))
	{
	    /* filename specified */
	    if(frame == NULL)
	    {
		/* mimic input filename for output */
		filename = argc[i];
		while(strchr(filename, '/') != NULL)
		{
		    /* how come basename doesn't work? */
		    filename = (char *) strchr(filename, '/') + 1;
		}

		strncpy(outfilename, filename, strlen(filename)-4);
		strcat(outfilename, ".sgi");

		filename = argc[i];
		frame = GetFrame(filename);
	    }
	}
	if(strstr(argc[i], ".rgb") || strstr(argc[i], ".sgi"))
	{
	    strcpy(outfilename, argc[i]);
	}
	if(argc[i][0] == '-')
	{
	    /* option specified */
	    if(strstr(argc[i], "-negate")) negate = 1;
	    if(strstr(argc[i], "-invert")) negate = 1;
	    if(strstr(argc[i], "-scale") && (argv >= (i+1)))
	    {
		scale = atof(argc[i+1]);
	    }
	    if(strstr(argc[i], "-zoom") && (argv >= (i+1)))
	    {
		zoom = atof(argc[i+1]);
	    }
	    if(strstr(argc[i], "-box") && (argv >= (i+4)))
	    {
		start_x = atol(argc[i+1]);
		start_y = atol(argc[i+2]);
		 stop_x = atol(argc[i+3]);
		 stop_y = atol(argc[i+4]);
	    }
	}
    }
    
    if(frame != NULL)
    {
	/* Print intentions */
	printf("converting %s to %s\n", filename, outfilename);
	printf("scaling size by %g\n", zoom);
	if(scale == 0)
	{
	    printf("autoscaling intensity \n");
	}
	else
	{
	    printf("intensity scale set to %g\n", scale);
	}
	printf("\n");	

	/* Filename */
	printf("reading %s ", filename);

	record_size = 4608;
	num_records = 2304;
	pixels_per_line = 2304;

	if((stop_x == 0)||(stop_x > pixels_per_line)) stop_x = pixels_per_line;
	if((stop_y == 0)||(stop_y > num_records))     stop_y = num_records;
	if(stop_x < start_x) {i = start_x; start_x = stop_x; stop_x = i;}
	if(stop_y < start_y) {i = start_y; start_y = stop_y; stop_y = i;}
	image_width  = stop_x - start_x;
	image_height = stop_y - start_y;
	
	printf("box from (%d,%d) to (%d,%d)\n", start_x, start_y, stop_x, stop_y);
    
	/* allocate space for memeory-resident image */
	image = (unsigned short int **) calloc(sizeof(unsigned short int *), image_height);

	/* need a little extra memory for swapping bytes */
	if(swap_bytes == TRUE)
	{
	    swapped = calloc(1, image_width*2);
	    printf("swapping bytes...\n");
	}

	/* load appropriate records (line) into memory */
	for(z = 0; z < image_height; ++z)
	{
	    /* allocate space for this line */
	    record = calloc(2, image_width);

	    /* position file at beginning of the desired scan line */
	    fseek(frame, ((z+start_y)*record_size)+(start_x*2)+512, SEEK_SET);
	
	    /* read the record into memory */
	    fread(record, 2, image_width, frame);
	
	    /* Account for byte swapping */
	    if(swap_bytes == TRUE)
	    {
		/* preserve location of read-in pixels */
		temp = record;
	    
		/* swab() can only COPY swapped bytes */
		swab((char*)record, (char*)swapped, image_width*2);
		record = swapped;
	    
		/* exchange buffers for next time */
		swapped = temp;
	    }

	    /* enter this record into the image */
	    image[z] = record;
	    
	    /* do statistics on pixel values */
	    for(x = 0; x < image_width; ++x)
	    {
/*		i = start_x + x;
		/* initialize from words in file */
		if(record[x] == 0xFFFF)
		{
		    ++overloads;
		}
		else
		{
		    sigma += ( (double) image[z][x]) * ((double) image[z][x]);
		}

		/* debug *
		printf("reading (%d,%d) ", z, x);
		printf("from (%d+%d)*%d+%d*2+512 ", z, start_y, record_size, start_x );
		printf(" = %d  \n", image[z][x]);  /*  */
	    }
	}
	
	/* image has now been completely read into memory */
	fclose(frame);
	
	/* do some statistics */
	sigma /= (double) (image_width * image_height);
	sigma = sqrt(sigma);
    
	/* autoscale output image to 5-sigma cuttoff */
	if(scale == 0)
	{
	    scale = 256.0/(5*sigma);
	    printf("intensity scale set to %g\n", scale);
	}
	
	
	
	
	
	/* construct RGB image header */
	SGIhead.Magic = 474; /* always 474 for RGB */
	SGIhead.storage = 0; /* uncompressed */
	SGIhead.bpc = 1; /* one byte per pixel (channel) */
	SGIhead.dimension = 2; /* 2D image */
	SGIhead.xsize = image_width*zoom; /*  */
	SGIhead.ysize = image_height*zoom; /*  */
	SGIhead.zsize = 1; /* one z-pixel (* channels) */
	SGIdata = (unsigned char *) calloc(SGIhead.xsize*SGIhead.ysize*SGIhead.zsize, SGIhead.bpc);
	strncpy(SGIhead.dummy1, SGIdata, 4);
	strncpy(SGIhead.dummy2, SGIdata, 404);
	SGIhead.pixmin = 0; /* minimum pixel value */
	SGIhead.pixmax = 255; /* maximum pixel value */
	strcpy(SGIhead.name, "no name"); /* set to origional image name */
	SGIhead.colormap = 0; /* grayscale/RGB */
    
	/* convert and copy all image data into new RGB file */
	printf("transforming data ...\n");
	for(z = 0; z < SGIhead.ysize; ++z)
	{
	    for(x = 0; x < SGIhead.xsize; ++x)
	    {
		/* debug *
		printf("xform (%d,%d) -> (%d,%d) = ", (int)(z/zoom), (int)(x/zoom), z, x);
		printf("%d\n", image[(int)(z*zoom)][(int)(x/zoom)]); /*  */

		/* transform frame image */
		sum = scale*image[(int)(z/zoom)][(int)(x/zoom)];
		if(sum > 255) sum = 255;
		if(negate) sum = 255-sum;
		
		SGIdata[(z*SGIhead.xsize)+x] = (unsigned char) sum;
	    }
	}
    
	/* write out the RGB image */
	rgbout = fopen(outfilename, "wb");
	if(rgbout != NULL)
	{
	    printf("writing %s\n", outfilename);
	    fwrite(&SGIhead, 512, 1, rgbout);
	    fwrite(SGIdata, SGIhead.bpc, SGIhead.xsize*SGIhead.ysize*SGIhead.zsize, rgbout);
	    fclose(rgbout);
	}
	else
	{
	    printf("ERROR: could not open %s\n", outfilename);
	}

    }
    else
    {
	printf("usage: %s framefile.img [outfile.rgb] [-zoom factor] [-box x1 y1  x2 y2] [-negate]\n", argc[0]);
		
	printf("\n\t framefile.img - the ADSC CCD image file you want to convert\n");
	printf(  "\t   outfile.rgb - name of output SGI RGB file.\n");
	printf(  "\t        factor - zoom in by this factor.\n");
	printf(  "\t  x1 y1  x2 y2 - corners of box to convert (original image pixel coordinates).\n");
	printf(  "\t       -negate - black spots on white background\n");
	printf("\n");
	printf(  "NOTE:\t denzo/mosflm (X,Y) detector coordinates are (Y,X) here.");
	printf("\n");
		
	/* report unhappy ending */
	return_code = 1;
    }
    
    return return_code;
}


FILE *GetFrame(char *filename)
{
    FILE *frame;
    unsigned char header[512];
    unsigned char *string;
    unsigned char *byte_order;
    
    typedef union
    {
	unsigned char string[2];
	unsigned short integer;
    } TWOBYTES;
    TWOBYTES twobytes;
    twobytes.integer = 24954;
    
    frame = fopen(filename, "rb");

    if(frame != NULL)
    {
	fread(header, 1, 512, frame);

	/* What kind of file is this? */
	string = header;

	if(0!=strncmp(string, "{\nHEADER_BYTES=  512;\nDIM=2;\nBYTE_ORDER=", 40))
	{
	    /* probably not an ADSC frame */

	    /* inform the user */
	    printf("ERROR: %s does not look like an ADSC frame!\n", filename);
	    /* skip this file */
	    fclose(frame);
	    
	    frame = NULL;
	}
	else
	{
	    /* determine byte order on this machine */
	    if(0==strncmp(twobytes.string, "az", 2))
	    {
		byte_order = "big_endian";
	    }
	    else
	    {
		byte_order = "little_endian";
	    }
	    
	    /* see if we will need to swap bytes */
	    string = (char *) strstr(header, "BYTE_ORDER=")+11;
	    if(0==strncmp(byte_order, string, 10))
	    {
		swap_bytes = FALSE;
	    }
	    else
	    {
		swap_bytes = TRUE;
	    }
	}
    }
    else
    {
	/* fopen() failed */
	perror("adsc2sgi");
    }
    
    return frame;
}