#ifndef lint
static char *RCSid = "$Id: binary.c,v 1.6 1994/09/13 16:11:18 alex Exp $";
#endif

/*
 * The addition of gnubin and binary, along with a small patch
 * to command.c, will permit gnuplot to plot binary files.
 * gnubin  - contains the code that relies on gnuplot include files
 *                     and other definitions
 * binary      - contains those things that are independent of those 
 *                     definitions and files
 *
 * With these routines, hidden line removal of your binary data is possible!
 *
 * Last update:  3/29/92 memory allocation bugs fixed. jvdwoude@hut.nl
 *               3/09/92 spelling errors, general cleanup, use alloc with no
 *                       nasty fatal errors
 *               3/03/92 for Gnuplot 3.24.
 * Created from code for written by RKC for gnuplot 2.0b.
 *
 * Copyright (c) 1991,1992 Robert K. Cunningham, MIT Lincoln Laboratory
 *
 */

#include <math.h>


#include "plot.h"   /* We have to get TRUE and FALSE */
#include "stdfn.h"
#include "binary.h"

/* 
 * This routine scans the first block of the file to see if the file is a 
 * binary file.  A file is considered binary if 10% of the characters in it 
 * are not in the ascii character set. (values < 128), or if a NUL is found.
 * I hope this doesn't break when used on the bizzare PC's.
 */
int
  is_binary_file(fp)
register FILE *fp;
{
  register int i,len;
  register int odd;                /* Contains a count of the odd characters */
  long where;
  register unsigned char *c;
  unsigned char buffer[512];

  if((where = ftell(fp)) == -1){ /* Find out where we start */
    fprintf(stderr,"Notice: Assuming unseekable data is not binary\n");
    return(FALSE);
  }
  else {
    rewind(fp);

    len = fread(buffer,sizeof(char),512,fp);
    if (len <= 0)	                  /* Empty file is declared ascii */
      return(FALSE);

    c = buffer;

    /* now scan buffer to look for odd characters */
    odd = 0;
    for (i=0; i<len; i++,c++) {
      if (!*c) {			  /* NUL _never_ allowed in text */
	odd += len;
	break;
      }
      else if ((*c & 128) ||/* Meta-characters--we hope it's not formatting */
	       (*c == 127)|| /* DEL */
	       (*c < 32 && 
		*c != '\n' && *c != '\r' && *c != '\b' &&
		*c != '\t' && *c != '\f' && *c != 27 /*ESC*/))
	odd++;
    }
  
    fseek(fp,where,0); /* Go back to where we started */

    if (odd * 10 > len)             /* allow 10% of the characters to be odd */
      return(TRUE);
    else
      return(FALSE);
  }
}
/*========================= I/O Routines ================================
  These may be useful for situations other than just gnuplot.  Note that I 
  have included the reading _and_ the writing routines, so others can create 
  the file as well as read the file.
*/

/*
  This function reads a matrix from a stream

  This routine never returns anything other than vectors and arrays
  that range from 0 to some number.  

*/
#define START_ROWS 100/* Each of these must be at least 1 */
#define ADD_ROWS 50
int
  fread_matrix(fin,ret_matrix,nr,nc,row_title,column_title)
FILE *fin;
float GPFAR * GPFAR * GPFAR *ret_matrix,GPFAR * GPFAR * row_title, GPFAR * GPFAR *column_title;
int *nr,*nc;
{
  float  GPFAR * GPFAR *m, GPFAR *rt, GPFAR *ct;
  register int num_rows = START_ROWS;
  register int num_cols;
  register int current_row = 0;
  register float  GPFAR * GPFAR *temp_array;
  float fdummy;
  
  fread(&fdummy,sizeof(fdummy),1,fin);
  num_cols = (int)fdummy;
  
  /* 
    Choose a reasonable number of rows,
    allocate space for it and continue until this space
    runs out, then extend the matrix as necessary.
    */
  ct = vector(0,num_cols-1);
  fread(ct,sizeof(*ct),num_cols,fin);

  rt = vector(0,num_rows-1);
  m = matrix(0,num_rows-1,0,num_cols-1);

  while(fread(&rt[current_row], sizeof(rt[current_row]), 1, fin)==1){ 
    /* We've got another row */
    if(fread(m[current_row],sizeof(*(m[current_row])),num_cols,fin)!=num_cols)
      return(FALSE);      /* Not a True matrix */

    current_row++;
    if(current_row>=num_rows){ /* We've got to make a bigger rowsize */
      temp_array = extend_matrix(m,0,num_rows-1,0,num_cols-1,
				 num_rows+ADD_ROWS-1,num_cols-1);
      rt = extend_vector(rt,0,num_rows-1,num_rows+ADD_ROWS-1);
      
      num_rows+= ADD_ROWS;
      m = temp_array;
    }
  }
  /*  finally we force the matrix to be the correct row size */
  /*  bug fixed. procedure called with incorrect 6th argument. jvdwoude@hut.nl */
  temp_array = retract_matrix(m,0,num_rows-1,0,num_cols-1,current_row-1,num_cols-1);
  /* Now save the things that change */
  *ret_matrix = temp_array;
  *row_title = retract_vector(rt, 0, num_rows-1, current_row-1);
  *column_title = ct;
  *nr = current_row;/* Really the total number of rows */
  *nc = num_cols;
  return(TRUE);
}

/* This writes a matrix to a stream 
   Note that our ranges are inclusive ranges--and we can specify subsets.
   This behaves similarly to the xrange and yrange operators in gnuplot
   that we all are familiar with.
*/
int
  fwrite_matrix(fout,m,nrl,nrh,ncl,nch,row_title,column_title)
register FILE *fout;
register float  GPFAR * GPFAR *m, GPFAR *row_title, GPFAR *column_title;
register int nrl,nrh,ncl,nch;
{
  register int j;
  float length;
  register int col_length;
  register int status;
  float  GPFAR *title = NULL;

  length = col_length = nch-ncl+1;

  if((status = fwrite((char*)&length,sizeof(float),1,fout))!=1){
    fprintf(stderr,"fwrite 1 returned %d\n",status);
    return(FALSE);
  }
  
  if(!column_title){
    column_title = title = vector(ncl,nch);
    for(j=ncl; j<=nch; j++)
      title[j] = j;
  }
  fwrite((char*)column_title,sizeof(float),col_length,fout);
  if(title){
    free_vector(title,ncl,nch);
    title = NULL;
  }

  if(!row_title){
    row_title = title = vector(nrl,nrh);
    for(j=nrl; j<=nrh; j++)
      title[j] = j;
  }
    
  for(j=nrl; j<=nrh; j++){
    fwrite((char*)&row_title[j],sizeof(float),1,fout);
    fwrite((char*)(m[j]+ncl),sizeof(float),col_length,fout);
  }
  if(title)
    free_vector(title,nrl,nrh);

  return(TRUE);
}

/*===================== Support routines ==============================*/

/******************************** VECTOR *******************************
 *       The following routines interact with vectors.
 *
 *   If there is an error we don't really return - int_error breaks us out.
 *
 *   This subroutine based on a subroutine listed in "Numerical Recipies in C",
 *   by Press, Flannery, Teukoilsky and Vetterling (1988).
 *
 */
float GPFAR *vector(nl,nh)
     register int nl,nh;
{
  register float GPFAR *vec;

  if (!(vec = (float GPFAR *)alloc((unsigned long) (nh-nl+1)*sizeof(float), NULL))){
    int_error("not enough memory to create vector",NO_CARET);
    return NULL;/* Not reached */
  }
  return (vec-nl);
}
/* 
 *  Free a vector allocated above
 *
 *   This subroutine based on a subroutine listed in "Numerical Recipies in C",
 *   by Press, Flannery, Teukoilsky and Vetterling (1988).
 *
 */
void 
  free_vector(vec,nl,nh)
float  GPFAR *vec;
int nl,nh;
{
  free(vec+nl);
}
/************ Routines to modify the length of a vector ****************/  
float  GPFAR *
  extend_vector(vec,old_nl,old_nh,new_nh)
float  GPFAR *vec;
register int old_nl,old_nh,new_nh;
{
  register float  GPFAR *new_v;
  new_v = (float GPFAR *)ralloc((void*)(vec+old_nl),
                       (unsigned long)(new_nh-old_nl+1)*sizeof(float),
		       "extend vector");
  return new_v-old_nl;
}

float  GPFAR *
  retract_vector(v,old_nl,old_nh,new_nh)
float  GPFAR *v;
register int old_nl,old_nh,new_nh;
{
  register float GPFAR *new_v;
  new_v = (float GPFAR *)ralloc((void*)(v+old_nl),
   (unsigned long)(new_nh-old_nl+1)*sizeof(float), "retract vector");
  return new_v-old_nl;
}
/***************************** MATRIX ************************
 *
 * 	  The following routines work with matricies
 *
 * 	 I always get confused with this, so here I write it down:
 * 			  for nrl<= nri <=nrh and
 * 			  for ncl<= ncj <=nch
 *  
 *   This matrix is accessed as:
 *   
 *     matrix[nri][ncj];
 *     where nri is the offset to the pointer to a vector where the
 *     ncjth element lies.
 * 
 *   If there is an error we don't really return - int_error breaks us out.
 *
 *   This subroutine based on a subroutine listed in "Numerical Recipies in C",
 *   by Press, Flannery, Teukoilsky and Vetterling (1988).
 *
 */
float 
   GPFAR * GPFAR *matrix(nrl,nrh,ncl,nch)
register int nrl,nrh,ncl,nch;
{
  register int i;
  register float GPFAR * GPFAR *m;

  m = (float GPFAR * GPFAR *)alloc((unsigned long)(nrh-nrl+1)*sizeof(float GPFAR *), "matrix");
  m -= nrl;

  for (i=nrl; i<=nrh; i++)
    {
      if(!(m[i] = (float GPFAR *) alloc((unsigned long)(nch-ncl+1)*sizeof(float), NULL))){
	free_matrix(m,nrl,i-1,ncl,nch);
	int_error("not enough memory to create matrix",NO_CARET);
	return NULL;
      }
      m[i] -= ncl;
    }
  return m;
}
/* 
 * Free a matrix allocated above
 *
 *
 *   This subroutine based on a subroutine listed in "Numerical Recipies in C",
 *   by Press, Flannery, Teukoilsky and Vetterling (1988).
 *
 */
void 
  free_matrix(m,nrl,nrh,ncl,nch)
float  GPFAR * GPFAR *m;
unsigned nrl,nrh,ncl,nch;
{
  register int i;

  for (i=nrl; i<=nrh; i++) 
    free((char GPFAR *) (m[i]+ncl));
  free((char GPFAR *) (m+nrl));
}
/*
  This routine takes a sub matrix and extends the number of rows and 
  columns for a new matrix
*/
float GPFAR * GPFAR *extend_matrix(a,nrl,nrh,ncl,nch,srh,sch)
     register float  GPFAR * GPFAR *a;
     register int nrl,nrh,ncl,nch;
     register int srh,sch;
{
  register int i;
  register float GPFAR * GPFAR *m;

  m = (float GPFAR * GPFAR *)ralloc((void*)(a+nrl),(unsigned long)(srh-nrl+1)*sizeof(float GPFAR *), "extend matrix");

  m -= nrl;

  if(sch != nch){
    for(i=nrl; i<=nrh; i++)
      {/* Copy and extend rows */
	if(!(m[i] = extend_vector(m[i],ncl,nch,sch))){
	  free_matrix(m,nrl,nrh,ncl,sch);
	  int_error("not enough memory to extend matrix",NO_CARET);
	  return NULL;
	}
      }
  }
  for(i=nrh+1; i<=srh; i++)
    {
      if(!(m[i] = (float GPFAR *) alloc((unsigned long) (nch-ncl+1)*sizeof(float), NULL))){
	free_matrix(m,nrl,i-1,nrl,sch);
	int_error("not enough memory to extend matrix",NO_CARET);
	return NULL;
      }
      m[i] -= ncl;
    }
  return m;
}
/*
  this routine carves a large matrix down to size
*/
float GPFAR * GPFAR *retract_matrix(a,nrl,nrh,ncl,nch,srh,sch)
     register float  GPFAR * GPFAR *a;
     register int nrl,nrh,ncl,nch;
     register int srh,sch;
{
  register int i;
  register float  GPFAR * GPFAR *m;

  for(i=srh+1; i<=nrh; i++) {
    free_vector(a[i],ncl,nch);
  }

  m = (float GPFAR * GPFAR *)ralloc((void*)(a+nrl), (unsigned long)(srh-nrl+1)*sizeof(float GPFAR *), "retract matrix");

  m -= nrl;

  if(sch != nch){
    for(i=nrl; i<=srh; i++)       
	if(!(m[i] = retract_vector(m[i],ncl,nch,sch))){ {/* Shrink rows */
	  free_matrix(m,nrl,srh,ncl,sch);
	  int_error("not enough memory to retract matrix",NO_CARET);
	  return NULL;
	}
      }
  }

  return m;
}

float 
   GPFAR * GPFAR *convert_matrix(a,nrl,nrh,ncl,nch)
float GPFAR *a;
register int nrl,nrh,ncl,nch;

/* allocate a float matrix m[nrl...nrh][ncl...nch] that points to the
matrix declared in the standard C manner as a[nrow][ncol], where 
nrow=nrh-nrl+1, ncol=nch-ncl+1.  The routine should be called with
the address &a[0][0] as the first argument.  This routine does
not free the memory used by the original array a but merely assigns
pointers to the rows. */

{
  register int i,j,ncol,nrow;
  register float GPFAR * GPFAR *m;

  nrow=nrh-nrl+1;
  ncol=nch-ncl+1;
  m = (float GPFAR * GPFAR *)alloc((unsigned long)(nrh-nrl+1)*sizeof(float GPFAR *), "convert_matrix");
  m -= nrl;

  m[nrl]=a-ncl;
  for(i=1,j=nrl+1;i<=nrow-1;i++,j++) m[j]=m[j-1]+ncol;
  return m;
}

void free_convert_matrix(b,nrl,nrh,ncl,nch)
float GPFAR* GPFAR *b;
register int nrl,nrh,ncl,nch;
{
	free((char*) (b+nrl));
}