/*
 * National Center for SuperComputing Applications, University of Illinois
 *
 * This NCSA software product is public domain software.  Permission
 * is hereby granted to do whatever you like with it. Should you wish
 * to make a contribution towards the production of this software, please
 * send us your comments about your experience with the software,  why
 * you liked or disliked it, how you use it, and most importantly, how it
 * helps your work. We will receive your comments at softdev@ncsa.uiuc.edu.
 *
 * Please send bug reports to bugs@ncsa.uiuc.edu
 *
 * Author: Eng-Whatt Toh, National Center for Supercomputing Applications
 *         ewtoh@ncsa.uiuc.edu
 */
#if ! defined(lint) && ! defined(LINT)
static char rcs_id[] = "$Id: td_hdf.c,v 1.3 1993/10/28 15:23:09 gbourhis Exp $";
#endif

/*
 *	File:		td_hdf.c
 *	Contents:	HDF file calls
 *
 */

#include "td_com.h"
#include "gr_def.h"

extern ApplicationData gr_Data;


/*
 *	Return 1 if file is a HDF file
 */
int
td_FileIsHdf(pathname)
char *pathname;
{
	if (DFishdf(pathname) == 0)
		return(1);
	else
		return(0);

}


/*
 *	Reset current pointer to SDS in HDF file
 */
int
td_HdfSDrestart()
{
	return(DFSDrestart());	/* make sure we always start at the first SDS */
}


/*
 *	Reset current pointer to R8 in HDF file
 */
int
td_HdfR8restart()
{
	return(DFR8restart());
}


/*
 *	Reset current pointer to palette in HDF file
 */
int
td_HdfPalrestart()
{
	return(DFPrestart());
}


/*
 *	Return number of SDSs in file.  If it is not a HDF file, return -1
 */
int
td_HdfCheck(fileName)
char	*fileName;
{
	int		ret;
	int		rank;
	int32	dims[3];
	char	pathName[STRNG160];

	sprintf(pathName,"%s/%s",td_dirName,fileName);

	td_HdfSDrestart();
	dims[0]=dims[1]=dims[2] = 1; /* init them myself */
	ret = DFSDgetdims(pathName,&rank,dims,3);
	if (ret != -1)
	{
		ret = DFSDnumber(pathName);
		td_HdfSDrestart();
	}
	
	return(ret);
}

/*
 *	Read data and if it isn't float32, convert it.
 *	After that, transpose the data if asked.
 *	Written by Gilles Bourhis. Jan 93.
 */
#define CONVERSION()				\
  do {						\
      switch (numtype & ~DFNT_NATIVE)		\
	{					\
	case DFNT_INT8:				\
	  CONVFROM(int8);			\
	  break;				\
	case DFNT_UINT8:			\
	  CONVFROM(uint8);			\
	  break;				\
	case DFNT_INT16:			\
	  CONVFROM(int16);			\
	  break;				\
	case DFNT_UINT16:			\
	  CONVFROM(uint16);			\
	  break;				\
	case DFNT_INT32:			\
	  CONVFROM(int32);			\
	  break;				\
	case DFNT_UINT32:			\
	  CONVFROM(uint32);			\
	  break;				\
	case DFNT_FLOAT64:			\
	  CONVFROM(float64);			\
	  break;				\
	default:				\
	  gr_TextMsgOut("Unknown SDS data type");\
	  ret = -1;				\
	}					\
     } while (0)

#define CONVFROM(TYPE)					\
  do { 							\
        int i;						\
        for (i = 0; i < dims[0]*dims[1]*dims[2]; i++)	\
	  notr_data[i] = (float32) ((TYPE *)indata)[i];	\
     } while(0)

static int
td_getFloatData(pathName, rank, dims, data)
char *pathName;
int rank;
int32 dims[3];
float32 *data;
{
  int32 numtype;
  int ret;
  float32 *notr_data;		/* no transposed data */

  if (gr_Data.transpose)
    notr_data = (float32 *)malloc(sizeof(float32)*dims[0]*dims[1]*dims[2]);
  else
    notr_data = data;
  if (notr_data == NULL)
    {
      gr_TextMsgOut("malloc failed\n");
      return -1;
    }

  DFSDgetNT(&numtype);

  if ((numtype|DFNT_NATIVE) == DFNT_NFLOAT32)
    ret = DFSDgetdata(pathName, rank, dims, (VOIDP)notr_data);
  else
    {
      int32 *indata;
      indata = (int32 *)malloc(DFKNTsize(numtype|DFNT_NATIVE)*
			       dims[0]*dims[1]*dims[2]);
      if (indata == NULL)
	{
	  gr_TextMsgOut("malloc failed\n");
	  return -1;
	}
      else
	ret = DFSDgetdata(pathName, rank, dims, (VOIDP)indata);

      CONVERSION();
      free((void *)indata);
    }
  if (ret != -1 && gr_Data.transpose)
    {
      int
	newdims[3], i,
	eltsize = sizeof(float32),
	isfortran  = 0;

      ret = transar(&rank, dims, newdims, &eltsize, (char *)notr_data,
		    (char *)data, &isfortran);
      free((void *)notr_data);
      for (i=0; i<rank; i++) dims[i] = newdims[i];
    }

  return ret;
}

/*
 *	Load specified SDS into memory
 */
int
td_HdfLoad(hdf,num)
A_Data_t	*hdf;
int			num;
{
	int i,ret;

	td_HdfSDrestart();
	hdf->numSDs = DFSDnumber(hdf->pathName);
	if ((num<=0) || (num > hdf->numSDs))
		return(-1);


	hdf->dims[0] = hdf->dims[1] = hdf->dims[2] = 1; /* init them myself */

	i=1; ret=0;
	while ((i<=num) && (ret != -1))
	{
		ret = DFSDgetdims(hdf->pathName,&(hdf->rank),hdf->dims,3);
		i++;
	}
	/* if the rank < 3 in the last call to DFSDgetdims, we need
	   to fill the last dimsize to 1. gbourhis. Feb 93 */
	for (i=hdf->rank; i<3; i++)
	  hdf->dims[i] = 1;

	if (ret == -1)
		return(-1);
	else
	{
		int32 newdims[3];

		newdims[1] = newdims[2] = 1;
		for (i=0; i<hdf->rank; i++) /* if transpose is asked, we    */
		  if (gr_Data.transpose) /* need to call td_Malloc3Dfloat32 */
		    newdims[i] = hdf->dims[hdf->rank-1-i]; /* with dims     */
		  else			/* transposed. gbourhis Jan 93      */
		    newdims[i] = hdf->dims[i];

		hdf->data = (float32 ***)
			td_Malloc3Dfloat32(newdims[0],newdims[1],newdims[2]);
		if (hdf->data != NULL)
		{
		  /* gbourhis Jan 93: call my function instead
		     of DFSDgetdata */
		        ret = td_getFloatData(hdf->pathName, hdf->rank,
					      hdf->dims, hdf->data[0][0]);

			if (ret == -1)
			{
				sprintf(msg,
				"***ERROR: Cannot read 3D SDS from %s.\n",
					hdf->pathName);
				gr_TextMsgOut(msg);
				gr_TextMsgOut
				("Are you sure this is a 3D 32-bit floating point dataset?\n");
				td_Free3d(hdf->data);
			}
			else
				hdf->format = HDF;
		}
		else
		{
			td_Free3d(hdf->data);
			ret = -1;
			sprintf(msg,
			"***ERROR: Not enough memory to load %s SDS #%d.\n",
					hdf->pathName,num);
			gr_TextMsgOut(msg);
			gr_TextMsgOut
			("I suggest you quit NCSA X DataSlice and retry the loading.\n");
		}

		return(ret);
	}
}


/*
 *	Calculate minimum and maximum values in SDS
 */
void
td_HdfMinMax(hdf)
A_Data_t	*hdf;
{
	int32 x,y,z;

	hdf->min = hdf->max = hdf->data[0][0][0];

	for (x=0;x<hdf->dims[0];x++)
		for (y=0;y<hdf->dims[1];y++)
			for (z=0;z<hdf->dims[2];z++)
			{
				if (hdf->data[x][y][z] < hdf->min)
					hdf->min = hdf->data[x][y][z];
				else
				if (hdf->data[x][y][z] > hdf->max)
					hdf->max = hdf->data[x][y][z];
			}
	gr_TextMsgOut("Automatically calculated min and max values.\n");
}

/*
 *	Get a scale and convert it to float32
 *	Written by Gilles Bourhis. Feb 93.
 */
#undef CONVFROM
#define CONVFROM(TYPE)				\
  do {						\
       int i;					\
       for (i=0; i<size; i++)			\
	 scale[i] = (float32) ((TYPE *)indata)[i];\
     } while (0)

static int
td_Hdfgetdimscale(dim, size, scale)
int dim;
int32 size;
float32 scale[];
{
  int32 numtype;
  int ret;

  DFSDgetNT(&numtype);
  if ((numtype|DFNT_NATIVE) == DFNT_NFLOAT32)
    return DFSDgetdimscale(dim, size, (VOIDP)scale);
  else
    {
      int32 *indata;
      indata = (int32 *)malloc(DFKNTsize(numtype|DFNT_NATIVE)*size);
      if (indata == NULL)
	{
	  gr_TextMsgOut("malloc failed\n");
	  return -1;
	}
      else
	ret = DFSDgetdimscale(dim, size, (VOIDP)indata);
      if (ret != -1) CONVERSION();
      free((void *)indata);
    }
  return ret;
}

/*
 *	Return min and max in float, if found by DFSDgetrange
 *	Written by Gilles Bourhis. Jan 93.
 */
#undef CONVFROM
#define CONVFROM(TYPE)				\
  do {						\
       TYPE minr, maxr;				\
       ret = DFSDgetrange((VOIDP)&maxr, (VOIDP)&minr);	\
       if (ret != -1) *max = (float32)maxr, 	\
	   *min = (float32)minr;		\
     } while(0)
       
static int
td_HdfgetRange(max, min)
float32 *max, *min;
{
  int32 numtype;
  int ret;

  DFSDgetNT(&numtype);
  if ((numtype|DFNT_NATIVE) == DFNT_NFLOAT32)
    return DFSDgetrange((VOIDP)max, (VOIDP)min);

  CONVERSION();
  return ret;
}

/*
 *	Return attributes of SDS in a single string
 *	calling routine must deallocate string when done
 */
char *
td_HdfgetStats(hdf,toRet,nColors,nSplitColors,format)
A_Data_t	*hdf;
Boolean		toRet;
int			nColors,nSplitColors;
A_FileFormat_t format;
{
	int	 i,j,ret;
	char *strng;
	char tempstr[256];
	char tlabel[256], tunits[256], tfmt[256];


	if (toRet == TRUE)
	{
		hdf->dataName[0] = hdf->dataUnit[0] = hdf->dataFmt[0] = '\0';

		DFSDgetdatastrs(hdf->dataName,hdf->dataUnit,hdf->dataFmt,msg);

		strng = td_Malloc1D(5000,1,sizeof(char),"td_Malloc1D:string");
		sprintf(strng,"Name: %s, Units: %s, Fmt: %s\n",
				hdf->dataName,hdf->dataUnit,hdf->dataFmt);
	}
	if (format == HDF)
	{
		int k;
		for (i=0;i<3;i++)
		{
			/* gbourhis Feb 93: if transpose is asked, we
			   need to use the original index to get
			   the scale and strings info from the HDF file */
			if (gr_Data.transpose && i < hdf->rank)
			  k = hdf->rank-1-i; /* transpose back */
			else
			  k = i;
			hdf->scaleStr[i] = (float32 *)td_Malloc1D(1,hdf->dims[i],
				sizeof(float32), "scaleStr");
			/* gbourhis Feb 93: call my function instead
			   of DFSDgetdimscale */
			ret = td_Hdfgetdimscale(k+1,hdf->dims[i],hdf->scaleStr[i]);
			if (ret == -1)
			{
				/* Cannot get scale, so implement own scale */
				for (j=0;j<hdf->dims[i];j++)
					hdf->scaleStr[i][j] = (float32)j;
			}
			hdf->label[i][0] = hdf->units[i][0] = hdf->fmt[i][0] = '\0';
			ret = DFSDgetdimstrs(k+1,tlabel, tunits, tfmt);

			hdf->label[i][0] = '\0';
			hdf->units[i][0] = '\0';
			hdf->fmt[i][0] = '\0';

			strncat (hdf->label[i],tlabel,40);
			strncat (hdf->units[i],tunits,40);
			strncat (hdf->fmt[i],tfmt,40);

			if ((ret != -1) && (toRet == TRUE))
			  {
				sprintf(tempstr,"%s: dim=%ld, label= %s, units=%s\n",
						axesLabels[i],hdf->dims[i],hdf->label[i],hdf->units[i]);
				strcat(strng, tempstr);
			  }
		}
		ret = td_HdfgetRange(&(hdf->max),&(hdf->min));
		if (ret == -1)
			td_HdfMinMax(hdf);
	}
	else
	if (format == MIRIAD)
		td_HdfMinMax(hdf);

	hdf->range = hdf->max - hdf->min;
	hdf->rangeFrac = (float32)((float32)nColors/hdf->range);
	hdf->rangeFracSplit = (float32)((float32)nSplitColors/hdf->range);

	if (toRet == TRUE)
	{
		sprintf(tempstr,"Max= %8.5f, Min=%8.5f\n", hdf->max,hdf->min);
		strcat(strng, tempstr);
		return(strng);
	}
	else
		return(NULL);
}


/*
 *	Write processed SDS to a file.  Return -1 if failed.
 */
int
td_HdfPutSDS(hdf)
A_Data_t *hdf;
{
	int i,ret;

	ret = DFSDsetdims(3, hdf->dims);
	if (ret == -1) return(ret);
	for (i=0;i<3;i++)
	{
		ret = DFSDsetdimstrs(i+1, hdf->label[i], hdf->units[i], 
			hdf->fmt[i]);
		if (ret == -1) return(ret);
	}
	ret = DFSDsetrange((VOIDP)&hdf->max,(VOIDP)&hdf->min);
	if (ret == -1) return(ret);

	ret = DFSDadddata(hdf->pathName,3,hdf->dims,(VOIDP)hdf->data);

	return(ret);
}


/*
 *	Write R8 image to a file.  Return -1 if failed.
 */
int
td_HdfPutImage(pathname,pic,palette,xsize,ysize)
char *pathname;
char *pic;
unsigned char *palette;
int	xsize,ysize;
{
	int	ret;

	ret = DFR8setpalette(palette);
	ret = DFR8putimage(pathname,pic,xsize,ysize,DFTAG_RLE);

	return(ret);
}


/*
 *	Add R8 image to a file.  Return -1 if failed.
 */
int
td_HdfAddImage(pathname,pic,palette,xsize,ysize)
char *pathname;
char *pic;
unsigned char *palette;
int	xsize,ysize;
{
	int	ret;

	ret = DFR8setpalette(palette);
	ret = DFR8addimage(pathname,pic,xsize,ysize,DFTAG_RLE);

	return(ret);
}


/*
 *	Get dimensions of R8 image in a HDF file.  Return -1 if failed.
 *	If successful, return number of images in HDF file.
 */
int
td_HdfgetRasDims(pathname,xsize,ysize,ispal)
char	*pathname;
int		*xsize,*ysize,*ispal;
{
	int	ret;
	int32   xs, ys;

	ret = td_HdfR8restart();

	if (ret != -1)
		ret = DFR8getdims(pathname,&xs,&ys,ispal);
	*xsize = xs;
	*ysize = ys;
	if (ret != -1)
		ret = DFR8nimages(pathname);

	td_HdfR8restart();

	return(ret);
}


/*
 *	Get data of R8 image in a HDF file.
 */
char
*td_HdfgetRasData(pathname,palette,xsize,ysize,ispal,toRead)
char	*pathname;
char	*palette;
int		*xsize,*ysize,*ispal;
Boolean toRead;
{
	int	ret;
	int32   xs, ys;
	unsigned char	*data=NULL;

	ret = DFR8getdims(pathname,&xs,&ys,ispal);
	*xsize = xs;
	*ysize = ys;

	if ((ret != -1) && (toRead))
	{
		data = (uint8 *)td_Malloc1D((int)(*xsize),(int)(*ysize),
			(long)(sizeof(unsigned char)),"td_Malloc1D");

		DFR8getimage(pathname,data,*xsize,*ysize,(unsigned char *)palette);
	}

	return (char *)data;
}


/*
 *	Get data of R8 image in a HDF file given its frameID
 */
char
*td_HdfgetRaster(pathname,data,palette,xsize,ysize,ispal,frameID)
char	*pathname;
char	*data;
char	*palette;
int		*xsize,*ysize,*ispal;
int		frameID;
{
	int32   xs, ys;

	if (DFR8readref(pathname,frameID) == -1)
		return(NULL);
	if (DFR8getdims(pathname,&xs,&ys,ispal) == -1)
		return(NULL);
	*xsize = xs;
	*ysize = ys;
	if (DFR8getimage(pathname,(uint8 *)data,*xsize,*ysize,
			 (uint8 *)palette) == -1)
		return(NULL);
	else
		return(data);
}


/*
 *	Get palette data of R8 image in a HDF file.
 */
int
td_HdfPalLoad(pathname,palette)
char *pathname;
char *palette;
{
	td_HdfPalrestart();

	return(DFPgetpal(pathname,palette));
}
