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

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

#define TRUE 1
#define FALSE 0

typedef union
{
	float real;
	unsigned char string[4];
	unsigned int integer;
} FOURBYTES;

/* function prototypes */
FOURBYTES fix_float(FOURBYTES wrong_order);


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

/* frame handling routines */
FILE *GetFrame(char *filename, char *header);
unsigned char header[2048];



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

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

    /* frame information variables */
    int record_size;
    int num_records;
    int pixels_per_line;
    int image_width;
    int image_height;
    float eight = 8;
    
    /* output image variables */
    int xsize;
    int ysize;

    /* indicies */
    int x, z;
    int in_x, in_z;
    int overloads = 0;
    
    /* limits */
    int high = 0;
    int low = 0xFFFF;
    double sum = 0;
    double sigma = 0;
    int stat_pixels = 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;
    int box_offset_x = 0;
    int box_offset_y = 0;

    /* storage */
    unsigned short int *record;
    unsigned short int *swapped;
    unsigned short int *temp;

    
    /* output PGM stuff */
    unsigned char *PGMdata;
    FILE *pgmout = NULL;
    char outfilename[256] = "";

    /* begin progam here */
    printf("OSC2PGM 1.0\t\tby James Holton 2-1-01\n\n");
    
    /* check argument list */
    for(i=1; i<argv; ++i)
    {
	if(strstr(argc[i], ".osc") || strstr(argc[i], ".stl") )
	{
	    /* filename specified */
	    if(frame == NULL)
	    {
		/* mimic input filename for output */
		filename = argc[i];
		while((char *) strchr(filename, '/') != NULL)
		{
		    /* how come basename doesn't work? */
		    filename = (char *) strchr(filename, '/') + 1;
		}

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

		filename = argc[i];
		frame = GetFrame(filename, header);
	    }
	}
	if(strstr(argc[i], ".pgm"))
	{
	    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_y = atol(argc[i+1]);
		start_x = atol(argc[i+2]);
		 stop_y = atol(argc[i+3]);
		 stop_x = 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);

	/* read file sizes from header */
	memcpy(&record_size, &header[784], 4);
	memcpy(&num_records, &header[788], 4);
	memcpy(&pixels_per_line, &header[768], 4);
	memcpy(&eight, &header[800], 4);
	
	/* most probable record sizes *
	record_size = 4096;
	num_records = 1900;
	pixels_per_line = 1900;
	eight = 8;
	/*  */

	/* handle weird box definitions */
	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;}
	/* default to whole image */
	if(stop_x == 0) stop_x = pixels_per_line;
	if(stop_y == 0) stop_y = num_records;

	/* preserve output image size (for later) */
	xsize = (int) (zoom* (double) (stop_x - start_x));
	ysize = (int) (zoom* (double) (stop_y - start_y));
	box_offset_x = start_x;
	box_offset_y = start_y;
	
	/* clip edges (for reading only) */
	printf("box from (%d,%d) to (%d,%d)\n", start_x, start_y, stop_x, stop_y);
	if(stop_x > pixels_per_line) stop_x = pixels_per_line;
	if(stop_y > num_records)     stop_y = num_records;
	if(start_x < 0)              start_x = 0;
	if(start_y < 0)              start_y = 0;
	image_width  = stop_x - start_x;
	image_height = stop_y - start_y;
	box_offset_x -= start_x;
	box_offset_y -= start_y;
	
	/* allocate space for memeory-resident image */
	image = (unsigned long int **) calloc(sizeof(unsigned long int *), image_height);
	record = calloc(2, image_width);

	/* 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 */
	    image[z] = calloc(4, image_width);

	    /* position file at beginning of the desired scan line */
	    fseek(frame, ((z+start_y+1)*record_size)+(start_x*2), 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 */
	    /* promote pixel words to actual values and */
	    /* do statistics on pixel values */
	    for(x = 0; x < image_width; ++x)
	    {
		/* initialize from words in file */
		if(record[x] == 0xFFFF)
		{
		    ++overloads;
		    /* simulate burn-thru effect? */
		    record[x] = 0;
		}
		if(record[x] >= 0)
		{
		    image[z][x] = (unsigned long) record[x];
		}
		if(record[x] < -1)
		{
		    image[z][x] = (unsigned long int) eight*(0x7fff & record[x]);
	    
		    /* this should never happen *
		    if(record[x] < -28672)
		    printf("oddball pixel!: (%4d, %4d) = %4d\n", x, z-1, record[x]); */
		}

		/* debug *
		printf("reading (%d,%d) = ", z, x);
		printf("%d  \n", image[z][x]);  /*  */

		/* now collect statistics */
		if((record[x] != 0)&&(record[x] != 0xFFFF))
		{
		    sigma += ( (double) image[z][x]) * ((double) image[z][x]);
		    ++stat_pixels;
		}
	    }
	}
	
	/* image has now been completely read into memory */
	fclose(frame);
	
	/* finish statistics */
	if(stat_pixels == 0)
	{
	    /* no data read? */
	    stat_pixels=1; sigma=1;
	}

	sigma /= (double) stat_pixels;
	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);
	}
	
	
	
	/* convert and copy all image data into new PGM file */
	PGMdata = (unsigned char *) calloc(xsize, ysize);

	
	printf("transforming data ...\n");
	for(z = 0; z < ysize; ++z)
	{
	    for(x = 0; x < xsize; ++x)
	    {
		/* transform output coordinates to input image frame (pixel units) */
		in_x = (int) (x/zoom) + box_offset_x;
		in_z = (int) (z/zoom) + box_offset_y;
		
		/* check for out-of-bounds data */
		if((in_x<0)||(in_z<0)||(in_x>=image_width)||(in_z>=image_height))
		{
		    /* slam non-image output to zero */
		    sum = 0.0;
		}
		else
		{
		    /* transform frame image */
		    sum = scale*image[in_z][in_x];
		}
		/* clip output intensity */
		if(sum > 255) sum = 255;
		/* optional negative image */
		if(negate) sum = 255-sum;
				
		PGMdata[(z*xsize)+x] = (unsigned char) sum;
	    }
	}
    
	/* write out the PGM image */
	pgmout = fopen(outfilename, "wb");
	if(pgmout != NULL)
	{
	    printf("writing %s\n", outfilename);
	    fprintf(pgmout, "P2\n%d %d\n255\n", ysize, xsize);
	    fprintf(pgmout, "# (%d,%d) to (%d,%d) on %s\n", start_x, start_y, stop_x, stop_y, filename);
	    fprintf(pgmout, "# pixels scaled by %f and zoomed in %fX\n", scale, zoom);
	    
	    /* print out bytes as text integers */
	    i=0;
	    for(x = xsize-1; x >= 0; --x)
	    {
		for(z = 0; z < ysize; ++z)
		{
		    fprintf(pgmout, "%3d ", PGMdata[(z*xsize)+x]);
		    
		    /* wrap long lines as per PGM definition */
		    ++i;
		    if((i%20==0)||((i%xsize) == (xsize-1))) fprintf(pgmout, "\n");
		}
	    }
	    
	    fclose(pgmout);
	}
	else
	{
	    printf("ERROR: could not open %s\n", outfilename);
	}

    }
    else
    {
	printf("usage: %s framefile.osc [outfile.pgm] [-zoom factor] [-box x1 y1  x2 y2] [-negate]\n", argc[0]);
		
	printf("\n\t framefile.osc - the RAXIS II/IV image file you want to convert\n");
	printf(  "\t   outfile.pgm - name of output PGM (portable graymap) 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 = 2;
    }
    
    return return_code;
}


FILE *GetFrame(char *filename, char *header)
{
    FILE *frame;
    unsigned char *string;
    unsigned char *byte_order;
    unsigned long int test;
    
    typedef union
    {
	unsigned char string[2];
	unsigned short integer;
    } TWOBYTES;

    TWOBYTES twobytes;
    FOURBYTES fourbytes;

    twobytes.integer = 24954;
    
    frame = fopen(filename, "rb");

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

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

	if(0!=strncmp(string, "R-AXIS", 6))
	{
	    /* probably not an OSC frame */

	    /* inform the user */
	    printf("ERROR: %s does not look like an Raxis II 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 */
	    memcpy(&fourbytes, &header[784], 4);

	    if(fourbytes.integer < 10000)
	    {
		swap_bytes = FALSE;
	    }
	    else
	    {
		swap_bytes = TRUE;
		
		/* byte-swap the header */
		printf("swapping header\n");
		string = (char*) calloc(1, 2048);
		swab((char*)header, (char*)string, 2048);
		memcpy(&header, &string, 2048);
		free(string);
	    }
	}
    }
    else
    {
	/* fopen() failed */
	perror("osc2pgm");
    }
    
    return frame;
}