/* matlab2hdf.c - transform a MATLAB binary file to HDF format */

#include <stdio.h>
#include <string.h>
#include "df.h"
#ifdef MAC
#  include <CursorCtl.h>
#else
   extern char *malloc();
#  ifdef SUN
   extern free();
#  else
   extern void free(void*);
#  endif
#endif
extern int loadmat();

usage()
{
#	define pr(X) fprintf(stderr,X)
	pr("usage: matlab2hdf [-rn1tTw] [-R min max] file1 [file2 ..]\n");
	pr("    transform MATLAB binary files filei into HDF files.\n");
	pr("Naming: if filei is of the form anything.m or anything.mat,\n");
	pr("    the new file is anything.hdf. Otherwise, '.hdf' is simply\n");
	pr("    appended to the old name.\n");
	pr("Option flags: -r make also Raster images, -n do Not write float data,\n");
	pr("    -1 include vectors and scalars, -t write Tracing output to stderr,\n");
	pr("    -T transpose, -w possibly Wrap around in making raster images,\n");
	pr("    -R min max: give data Range for all raster images (min,max: float).\n");
	exit(-1);
}

#define real double
real min=0.,max=1.;
int PresetRange = 0; /*false*/
int Rasters = 0;
int SDSs = 1;
int Include1x1 = 0;
int Tracing = 0;
int Transposing = 0;
int DontWrap = 1;

mknewfilename(old,new)
char *old,*new;
{
	int n;
	n = strlen(old);
	if (old[n-2] == '.' && (old[n-1] == 'm' || old[n-1] == 'M')) {
		strcpy(new,old);
		new[n-1] = 'h'; new[n] = 'd'; new[n+1] = 'f'; new[n+2] = '\0';
	} else if (old[n-4] == '.' && (old[n-3] == 'm' || old[n-3] == 'M')
				&& (old[n-2] == 'a' || old[n-2] == 'A')
				&& (old[n-1] == 't' || old[n-1] == 'T') ) {
		strcpy(new,old);
		new[n-3] = 'h'; new[n-2] = 'd'; new[n-1] = 'f';
	} else {
		strcpy(new,old);
		strcat(new,".hdf");
	}
}

static chk(code)
{
	if (code) {
		fprintf(stderr,"*** Unexpected HDF error %d, DFerror=%d\n",code,DFerror);
		exit(code);
	}
}

static memerr() {fprintf(stderr,"*** Not enough memory\n"); exit(-1);}

float32* mkfloat(xp,n,m)
double *xp; int n,m;
{
	size_t nm = (size_t)n*(size_t)m;
	float32* result;
	size_t i,j;
	result = (float32*) malloc(nm*sizeof(float32));
	if (result==0) memerr();
	if (Transposing)
		for (i=0; i<m; i++) for (j=0; j<n; j++)
			result[j*m+i] = (float32)xp[i*n+j];
	else
		for (i=0; i<nm; i++) result[i] = (float32)xp[i];
	return result;
}

unsigned char *mkraster(xp,n,m)
double *xp; int n,m;
{
	size_t nm = (size_t)n*(size_t)m;
	size_t i,j;
	real maxval,minval,scaling,x;
	unsigned char *result;
	result = (unsigned char *)malloc(nm*sizeof(unsigned char));
	if (result==0) memerr();
	if (PresetRange) {maxval=max; minval=min;} else {
		maxval = -1.234E33; minval = -maxval;
		for (i=0; i<nm; i++) {
			if (xp[i]<minval) minval=xp[i];
			if (xp[i]>maxval) maxval=xp[i];
		}
	}
	scaling = 253.0/(maxval - minval);
	if (PresetRange && DontWrap)
		if (Transposing)
			for (i=0; i<m; i++) for (j=0; j<n; j++) {
				x = xp[i*n+j];
				if (x<minval) x=minval; if (x>maxval) x=maxval;
				result[j*m+i] = 1+(unsigned char)((x-minval)*scaling);
			}
		else
			for (i=0; i<nm; i++) {
				x = xp[i];
				if (x<minval) x=minval; if (x>maxval) x=maxval;
				result[i] = 1+(unsigned char)((x-minval)*scaling);
			}	
	else
		if (Transposing)
			for (i=0; i<m; i++) for (j=0; j<n; j++)
				result[j*m+i] = 1+(unsigned char)((xp[i*n+j]-minval)*scaling);
		else
			for (i=0; i<nm; i++)
				result[i] = 1+(unsigned char)((xp[i]-minval)*scaling);
	return result;
}

ProcessFile(fn)
char *fn;
{
	FILE *fp;
	char name[20]; char newfilename[80];
	int type, mrows, ncols, imagf, ret;
	uint16 ref;
	double *xr, *xi; float *xfr,*xfi; unsigned char *xrr, *xri;
	int32 dims[2];
	size_t i,j;
	fp = fopen(fn,"r");
	if (fp == 0) {fprintf(stderr,"*** File %s does not exists\n",fn); return 0;}
	mknewfilename(fn,newfilename);
	if (Tracing) {
		fprintf(stderr,"%s -> %s",fn,newfilename);
	}
	remove(newfilename);
	while (loadmat(fp,&type,name,&mrows,&ncols,&imagf,&xr,&xi)==0) {
#		ifdef MAC
		SpinCursor(-32);
#		endif
		if (Transposing)
			{dims[0] = mrows; dims[1] = ncols;}
		else
			{dims[1] = mrows; dims[0] = ncols;}
		if (SDSs && (Include1x1 || mrows>1 && ncols>1)) {
			if (imagf) {
				char name2[25];
				xfr = mkfloat(xr,mrows,ncols);
				xfi = mkfloat(xi,mrows,ncols);
				strcpy(name2,name);
				strcat(name2,"_re");
				chk(DFSDclear());
				chk(DFSDsetdims(2,dims));
				chk(DFSDsetdatastrs(name2,"","",""));
				chk(DFSDadddata(newfilename,2,dims,(float32*)xfr));
				if (Tracing) fprintf(stderr," :%s (SDS)",name2);
				ref = DFSDlastref();
				chk(DFANputlabel(newfilename,DFTAG_SDG,ref,name2));
				strcpy(name2,name);
				strcat(name2,"_im");
				chk(DFSDclear());
				chk(DFSDsetdims(2,dims));
				chk(DFSDsetdatastrs(name2,"","",""));
				chk(DFSDadddata(newfilename,2,dims,(float32*)xfi));
				if (Tracing) fprintf(stderr," :%s (SDS)",name2);
				ref = DFSDlastref();
				chk(DFANputlabel(newfilename,DFTAG_SDG,ref,name2));
				free(xfr); free(xfi);
			} else {
				xfr = mkfloat(xr,mrows,ncols);
				chk(DFSDclear());
				chk(DFSDsetdims(2,dims));
				chk(DFSDsetdatastrs(name,"","",""));
				chk(DFSDadddata(newfilename,2,dims,(float32*)xfr));
				if (Tracing) fprintf(stderr," :%s (SDS)",name);
				ref = DFSDlastref();
				chk(DFANputlabel(newfilename,DFTAG_SDG,ref,name));
				free(xfr);
			}
		}
		if (Rasters && (Include1x1 || mrows>1 && ncols>1)) {
			if (imagf) {
				char name2[25];
				xrr = mkraster(xr,mrows,ncols);
				xri = mkraster(xi,mrows,ncols);
				strcpy(name2,name);
				strcat(name2,"_re");
				if (Transposing)
					chk(DFR8addimage(newfilename,xrr,ncols,mrows,0));
				else
					chk(DFR8addimage(newfilename,xrr,mrows,ncols,0));
				if (Tracing) fprintf(stderr," :%s (RIG)",name2);
				ref = DFR8lastref();
				chk(DFANputlabel(newfilename,DFTAG_RIG,ref,name2));
				strcpy(name2,name);
				strcat(name2,"_im");
				if (Transposing)
					chk(DFR8addimage(newfilename,xri,ncols,mrows,0));
				else
					chk(DFR8addimage(newfilename,xri,mrows,ncols,0));
				if (Tracing) fprintf(stderr," :%s (RIG)",name2);
				ref = DFR8lastref();
				chk(DFANputlabel(newfilename,DFTAG_RIG,ref,name2));
				free(xrr); free(xri);
			} else {
				xrr = mkraster(xr,mrows,ncols);
				if (Transposing)
					chk(DFR8addimage(newfilename,xrr,ncols,mrows,0));
				else
					chk(DFR8addimage(newfilename,xrr,mrows,ncols,0));
				if (Tracing) fprintf(stderr," :%s (RIG)",name);
				ref = DFR8lastref();
				chk(DFANputlabel(newfilename,DFTAG_RIG,ref,name));
				free(xrr);
			}
		}
		free(xr);
		if (imagf) free(xi);
	}
	if (Tracing) fprintf(stderr,"\n");
	fclose(fp);
}

scanflags(s)
char *s;
{
	int i;
	for (i=0; s[i]; i++) switch (s[i]) {
		case 'r': Rasters = 1; break;
		case 'n': SDSs = 0; break;
		case '1': Include1x1 = 1; break;
		case 't': Tracing = 1; break;
		case 'T': Transposing = 1; break;
		case 'w': DontWrap = 0; break;
		default : fprintf(stderr,"*** Unknown flag %c\n",s[i]);
	}
}

main(argc,argv)
int argc; char *argv[];
{
	int i,start=1;
	float Min,Max;
#	ifdef MAC
	InitCursorCtl(0);
#	endif
	if (argc <= 1) usage();
	if (argv[1][0] == '-' && argv[1][1] != 'R') {scanflags(argv[1]+1); start++;}
	if (start>=argc) usage();
	if (argv[start][0] == '-' && argv[start][1] == 'R') {
		start+= 3; if (start>=argc) usage();
		if (sscanf(argv[start-2],"%f",&Min)!=1) usage();
		if (sscanf(argv[start-1],"%f",&Max)!=1) usage();
		if (min < max) {min = Min; max = Max; PresetRange = 1;}
		else fprintf(stderr,"*** min>=max, using default scaling...\n");
	}
	for (i=start; i<argc; i++) ProcessFile(argv[i]);
	return 0;
}