/***************************************************************************
 File Name  : image.c
 Purpose    : provides basic functionality to read, write, display, enhance
	      and correlate images ( SUN raster or BIN format )
 Release    : Version 1.0
 Date	    : Aug 31,1995

GSNAKE API is jointly developed by the Information Technology Institute (ITI), Singapore, and the School of Applied Science, Nanyang Technological
University (NTU), Singapore. 

These software programs are available to the user without any license or royalty fees. Permission is hereby granted to use, copy, modify, and distribute this software and its documentation for any purpose. ITI and NTU gives no warranty, express, implied, or statuary for the software and/or documentation provided, including, without limitation, waranty of merchantibility and warranty of fitness for a particular purpose. The software provided hereunder is on an "as is"  basis, and ITI and NTU has no obligation to provide maintenance, support, updates, enhaucements, or modifications.

GSNAKE API is available for any UNIX compatible system. User feedback, bugs, or software and manual suggestions should be sent via electronic mail to one of the following, who may or may not act on them as he/she desires :

		asschan@ntu.ac.sg
		kflai@iti.gov.sg
**************************************************************************/

#include "gsnake.h"
#include "xwindow.h"


/* Constructor and destructor methods */

IMAGE::IMAGE(void)
{
	data = NULL;
	ximg = NULL;
	row = col = 0;
}

IMAGE::~IMAGE(void)
{
	reset();
	
}

/* reset _iMemFrees previously allocated memory. Allows reuse of IMAGE */
void IMAGE::reset(void)
{

	_iMemFree(data);
	if(ximg) XDestroyImage(ximg);
	data=NULL;
	ximg=NULL;
	row = col = 0;
	
}


/* Method init allocates memory for matrix. It returns an integer 1/0 
   informing of memory allocation error if any */

int IMAGE::init(int _row, int _col)
{
	row = _row;
	col = _col;
	
	if( !(data = (float *) _iMemAlloc(_col*_row*sizeof(float))) ) {
                fprintf(stderr, "<IMAGE::init> : memory allocation error\n");
                return MEMORYERROR;
        }
        memset(data, 0, _col*_row*sizeof(float)) ;

	return NOERROR;
}


/* Method init_gauss33 performs memory allocation and initialises a gaussian
   3x3 conditioning matrix. */

void IMAGE::initAsGauss()
{	
	float term[3];
	int m,n;

	
	init(3,3);		/* initialise new matrix */	
	term[0] = term[2] = 0.2236;
	term[1] = 0.5477;
	
	for ( m = 0 ; m < row ; m++)
		for ( n = 0 ; n < col ; n++)
			put(m,n,term[m]*term[n]);
}
 

/* Method print prints out all values in matrix */

void IMAGE::print()
{
	register short i, j;

	printf("row : %d  col : %d\n", row, col);

	for (i = 0; i<row; i++) {
		for (j = 0; j<col; j++)
			printf("%f ", get(i,j));
		printf("\n");
	}
}

/* Method read reads data in char form from filename. File type can be of 
   image _bin file or sun raster file image. Returns error if otherwise. The    
   image data is then converted to a floating point form. Read returns
   an integer value indicating error if any .
   Modified to be able to read raster file images -- 8/2	
*/

IMAGE::read(char *filename)
{
	FILE *bin ;
	short format=_bin;
	unsigned char *buf ;
	unsigned char word[2];
	register int i, j ;
	unsigned char *cptr ;
	float *fptr ;
	int filecheck ;
	

	bin = fopen(filename, "rb") ;
	if( !bin ) {
		fprintf(stderr, "<IMAGE::read> : %s cannot be read\n", 
				filename) ;
		return FILEIOERROR ;
	}
       
	fread(&filecheck,sizeof(int),1,bin);
	rewind(bin);
	if (filecheck==0x59a66a95) (format= _ras);
	/* read file in BIN format */
	
	if (format == _bin)
	{
		/* read columns and rows */
		fread(word, sizeof(char), 2, bin) ;
		col = word[1]*256 + word[0] ;
		fread(word, sizeof(char), 2, bin) ;
		row = word[1]*256 + word[0] ;

		if( (init(row, col)) || !( buf = (unsigned char *)
				 _iMemAlloc(row*col*sizeof(unsigned char)))) {
                	fprintf(stderr, 
			    "<IMAGE::init> : memory allocation error\n");
			return MEMORYERROR ;
		}

		fread(buf, sizeof(char), row*col, bin);
	
		/* casting char to float */
		for(i=0, cptr = buf, fptr = data; i<row; i++)
			for(j=0; j<col; j++, cptr++, fptr++)
				*fptr = (float) *cptr ;
	}


	/* Reading raster file format */

	if (format == _ras)
	{
		int header[8];
		int length ;
		int depth ;
	
		for (i=0;i<8;i++) header[i]=0;

		/* read the raster file header = 8 integers */
		fread(header, sizeof(int), 8, bin) ;  
		fseek(bin, header[7], 1);
		col = header[1];
		row = header[2];
		length = header[4] ;
		depth = header[3];
		
		if ( (init(row, col)) || 
		    !(buf=(unsigned char*) _iMemAlloc(length) ) ) {

                	fprintf(stderr, 
			    "<IMAGE::init> : memory allocation error\n");
		    	return MEMORYERROR ;
		}
		fread(buf, 1, length, bin) ;

		if( depth == 8 ) {
		    for (i=0,cptr=buf, fptr=data; i<row; i++) {
			for (j=0;j<col;j++, cptr++, fptr++)
				*fptr= (float) *cptr;
			/* Raster file rounds up to nearest 16 bits.
			   For odd cols, must increase by 1 byte */
			if (col%2) cptr++;
		    }
		}
		else if (depth == 1) {
		    unsigned char c  ;
		    for(i=0, cptr=buf, fptr=data; i<row; i++) 
		        for(j=0, c=128; j<col; fptr++) {
		        	*fptr = ( *cptr & c ) ? GRAYLEVELRANGE-1 : 0.0 ;
		        	++j ;
		        	if( !(j%8) ) {
		        		c = 128 ;
		        		cptr++;
		        	}
		        	else ( c = c >> 1 ) ;
		        }
		}
		else 
		   fprintf(stderr, 
		   "<IMAGE::init> : do not know how to handle depth = %d\n", 
		   	depth) ;
	};

	_iMemFree(buf);
	fclose(bin);
	return NOERROR;
}

/* Write saves the image data onto file, specified by file type format. either
   bin or ras. Returns FILEIOERROR and MEMORYERROR if either occurs, other
   wise returns NOERROR */

IMAGE::write(char *filename, short format)
{
	FILE *bin;
	short row, col;
	register int i,j;
	unsigned char word[2];
	float min,max,prange;
	float *fptr;
	unsigned char *cptr, *buf;

	
	bin = fopen(filename, "wb") ;
	if( !bin ) {
		fprintf(stderr, "<IMAGE::write> : %s cannot be opened\n", 
				filename) ;
		return FILEIOERROR ;
	}
	
	col = getCol();
	row = getRow();

	if (!( buf = (unsigned char*) 
			_iMemAlloc(row*col*sizeof(unsigned char) )))
		return MEMORYERROR;

	/* Perform linear stretching of pixel values to full range of 255 */
	max=*data;
	min=*data;

	/* Get Max and Min pixel values in Image */
	for(i=0, fptr=data; i<row; i++)
		for(j=0; j<col; j++,fptr++)
		{
			if (*fptr > max) max = *fptr;
			else if (*fptr < min) min =*fptr;
		}

	/* Map to full range */	
	prange = max - min;	
	for(i=0, fptr=data; i<row; i++)
		for(j=0; j<col; j++,fptr++)
			*fptr = (*fptr-min)/prange * 255 ;	


	/* casting float to char */
	for(i=0, cptr=buf,fptr = data; i<row; i++)
		for(j=0; j<col; j++, fptr++,cptr++)
			*cptr = (unsigned) *fptr  ;
							
	if (format == _bin)
	{
		/* Write row and col info */
		printf("\nWriting to binary file ...\n");
		word[1] = col / 256;
		word[0] = col - word[1]*256;
		fwrite(word, sizeof(char), 2, bin);

		word[1] = row /256;
		word[0] = row - word[1]*256;
		fwrite(word, sizeof(char), 2, bin);
		
		/* write rest of data */
		fwrite(buf, sizeof(unsigned char), row*col, bin);
	}

	else 

	{	int header[8];
		char odd = col%2 ;	
	
		for (i=0;i<8;i++) header[i]=0;
		printf("\nWriting to sun raster file ... \n");
		header[0] = 0x59a66a95;	    /* Magic number */
		header[1] = col;	    /* ras_width */	
		header[2] = row;	    /* ras_height */	
		header[3] = 8;		    /* ras_depth */
		header[4] = row * col;      /* ras_length */		    
		header[5] = 1;		    /* ras_type */
		header[6] = 1; 		    /* ras_maptype */
		header[7] = 0;		    /* ras_maplength */
	

		fwrite(header, sizeof(int),8,bin);
		for (i=0,cptr=buf;i<row;i++,cptr+=col) {
			fwrite(cptr, sizeof(unsigned char),col,bin);
			if (odd) fputc(0,bin);
		}

	};		
	_iMemFree(buf);
	fclose(bin);
	return NOERROR;
}

	

/* Method show creates an Xwindow to show the image.
   Requires parameter blowup as magnification factor */

void IMAGE::show(unsigned char blowup,int num,int h_offset, int v_offset)
{

	/* raise or resize the window to fit */
	if (!h_offset && !v_offset) 
		xwin_raiseWindow(blowup*col*num, blowup*row);
		
	generateX(blowup);
	xwin_drawImg(ximg,h_offset*blowup,v_offset*blowup);
}

void IMAGE::generateX(unsigned char blowup)
{
 
	if( ximg ) XDestroyImage(ximg);
        ximg = xwin_creatXimage(data, row, col, blowup);

}

IMAGE *IMAGE::correlate(
        IMAGE *Template,
        int RowStep,
        int ColStep,
	char verbose)
{
        IMAGE *Output;
        register short im, in, om, on;
        short newrow, newcol;

	if(verbose) {
		fprintf(stdout, "\t\tCorrelating Image : ");
		fflush(stdout);
	}

        /* create a new output buffer */
        newrow = ROUNDOFF( (float) (getRow()) / (float) (RowStep) );
        newcol = ROUNDOFF( (float) (getCol()) / (float) (ColStep) );
        Output = new IMAGE;
	Output->init(newrow, newcol);

        if( !Output ) {
                fprintf(stderr, "correlate : Cannot create new image\n");
                return NULL;
        }

        /* main correlation loop */
        for (om = im = 0 ; om < newrow ; om++, im += RowStep) {
	    if( verbose && !(om%10)) {
		fprintf(stdout, "#");
		fflush(stdout);
	    }
            for (on =in = 0 ; on < newcol ; on++, in += ColStep) {
                /* correlation is centered at (im, in) at input image */
                Output->put(om, on, innerProd(Template, im, in));
            }
        }

	if(verbose)
		fprintf(stdout, "\n");

        return Output;
}

/*
 * Calculate the inner product (element by element product) between
 * template and portion of input. This is the primary operation
 * used in correlation and convolution
 */
double IMAGE::innerProd(
        IMAGE *Template,
        int InpRowIndex,        /* center of box in input matrix  */
        int InpColIndex)
{
        register short i, j;
        double RESULTS;
        short tm_off, tn_off, _row, _col;
        short m_ind, n_ind;
        short temp_m, temp_n;

        tm_off = (_row = Template->getRow())/ 2;
        tn_off = (_col = Template->getCol())/ 2;

        RESULTS = 0;
        for( i = 0  ; i < _row ; i++) /* correlates */
            for( j = 0 ; j < _col ;  j++) {

                m_ind = InpRowIndex + i - tm_off ;      /* index into input  */
                n_ind = InpColIndex + j - tn_off ;

                /* cap location to within input image range */
                temp_m = CAP(0, m_ind, row-1);
                temp_n = CAP(0, n_ind, col-1);

                RESULTS += Template->get(i,j) * get(temp_m, temp_n);

              }

        return RESULTS;
}


/* Method condition performs image conditioning based on histogram equalisation
   Parameters required are  */
   

void IMAGE::condition(double low_pct, double high_pct, double low_val, double high_val, double imm_pow) 	

{
	int histo[GRAYLEVELRANGE+1] ;      /* histogram for edge conditioning */
        double mapping[GRAYLEVELRANGE+1], cu[GRAYLEVELRANGE+1] ;
        double size = row*col ;
        register int i, j ;
	
        memset(histo, 0, (GRAYLEVELRANGE+1)*sizeof(int) ) ;
	for (i=0; i<row; i++)
	    for (j=0; j<col; j++)
		++histo[(int) CAP(0, *(data + i*col + j), GRAYLEVELRANGE)];

	/* get cumulative probability */
        memset(cu, 0, (GRAYLEVELRANGE+1)*sizeof(double)) ;
        for(i=0; i<GRAYLEVELRANGE+1; i++) {
                for(j=0; j<=i; j++)
                        cu[i] += histo[j] ;
                cu[i] /= size ;
        }

	/* obtain mapping table */
        mapping[0] = 0 ;
        for(i=1; i<GRAYLEVELRANGE+1; i++)  {

           if( cu[i] <= low_pct )
                mapping[i] = cu[i]/low_pct*low_val ;
           else if( cu[i] <= high_pct )
                mapping[i] = pow((cu[i]-low_pct)/(high_pct-low_pct), imm_pow)
                                        *(high_val-low_val)+low_val ;
           else
                mapping[i] = CAP(0, (cu[i]-high_pct)/(1-high_pct)*(1-high_val)
                                        + high_val, 1) ;
        }

        for(i=0; i<row; i++)
             for(j=0; j<col; j++)
        	*(data + i*col + j) =
                     mapping[(int)CAP(0, *(data + i*col + j), GRAYLEVELRANGE)];
}

/*
 * cut portion of give image. If memory allocation error will return a 
   NULL pointer , otherwise returns pointer to new image. If sx and sy 
   are illegal, will return NULL. Otherwise, it will size the box to the
   bounds of the image
 */

 
IMAGE *IMAGE::cut(int sx,int sy,int length, int height)
{
	
	IMAGE *segment;
        float *tptr,*sptr;
        int i,j;

	/* return NULL if illegal parameter */
	if (  (sx>col) || (sy>row)  ) return NULL;
	
	/* Trap cut to limit to bounds of image size */
	if (sx+length > col) length = col - sx;
	if (sy+height > row) height = row - sy;

	if ( !( segment = new IMAGE ) ) {
		fprintf(stderr,
			"<IMAGE::cut>: memory allocation error\n");
		return NULL;
	}
	 
        if (!(segment->init(height, length)) ) {

               	sptr = data+col*sy+sx;
               	tptr = segment->data;

               	for ( i=0;i<height ;i++) {

		   for (j=0;j<length;j++,sptr++,tptr++) 
                         	*tptr = *sptr;
                       	sptr += col - length;
                }
                return(segment);
        }
}


/* Method fill fills a specified square area with a value */
/* If all co-ordinate parameters are 0, will fill the whole image */

int IMAGE:: fill(float val,int sx, int sy, int length, int height)

{int i,j;

	if (!length && !height && !sx && !sy )
	{
		length = col;
		height = row;
	};
	
	if ( (row<(sx+length)) || (col<(sy+height)) )
		return PARAMERROR;
	
	for (i=sy;i<height;i++)
		for(j=sx;j<length;j++)
			put(i,j,val);

	return NOERROR;
};