/*
 * anti-aliased isoplot into a file
 * does this really obey window?
 */
#include <u.h>
#include <libc.h>
#include <libg.h>
#include <fb.h>
#include <stdio.h>
float **data;
int ndata;
typedef struct Pixel{
	unsigned char m, a;
}Pixel;
Pixel **pic;
int wid=640, hgt=512;
typedef struct Acc{
	float accm, acca;
	float m, a;
}Acc;
Acc *col;
int lo, hi;
float *line, *linealpha, *stripe, *stripealpha;
#define	NLIGHT	16
int nlight=1;
struct light{
	float lx, ly, lz;	/* direction to light */
	float lum;		/* brightness of light */
	float sx, sy, sz;	/* direction to hilite */
}light[NLIGHT]={
	2, 3, 9, 1, 0, 0
};
float spec=0;			/* amount of specular reflection */
float bump=80;			/* exponent of specular bump */
float diff=.98;			/* amount of diffuse reflection */
float ambient=.02;		/* amount of ambient light */
double range;
void erread(int f, char *p, int n){
	if(read(f, p, n)!=n){
		fprint(2, "aplot: eof!\n");
		exits("eof");
	}
}
#define	cvt(name, type)\
void name(int f, float *fp, int n)\
{\
	type line[8192];	/* bug, could be too small! */\
	type *lp=line;\
	int i;\
	erread(f, (char *)line, n*sizeof(type));\
	for(i=0;i!=n;i++)\
		*fp++=*lp++;\
}
cvt(cint, int);
cvt(cshort, short);
cvt(clong, long);
cvt(cdouble, double);
cvt(cchar, char);
cvt(cuchar, unsigned char);
void cfloat(int f, float *fp, int n){
	erread(f, (char *)fp, n*sizeof(float));
}
void getlights(char *);
int size(int, int);
void fixlights(void);
void adraw(int);
void *emalloc(long n){
	void *v=malloc(n);
	if(v==0){
		fprint(2, "aplot: no space\n");
		exits("no space");
	}
	return v;
}
main(int argc, char *argv[]){
	float hi=0, scale;
	register i, j, f;
	char *picargv[10];
	char window[1024];
	char *line;
	void (*crack)(int, float *, int);
	int x0=0, y0=0;
	PICFILE *outf;
	/* Crack flags */
	if(getflags(argc, argv, "al:1[lightfile]t:1[type]r:1[range]w:4[x0 y0 x1 y1]")!=2)
		usage("infile");
	if(flag['l']) getlights(flag['l'][0]);
	if(flag['r']) range=atof(flag['r'][0]);
	if(flag['w']){
		x0=atoi(flag['w'][0]);
		y0=atoi(flag['w'][1]);
		wid=atoi(flag['w'][2]-x0);
		hgt=atoi(flag['w'][3]-y0);
	}
	if(wid<=0 || hgt<=0){
		fprint(2, "%s: output window too small (%dx%d)\n", argv[0], wid, hgt);
		exits("window too small");
	}
	pic=emalloc(hgt*sizeof(Pixel *));
	pic[0]=emalloc(wid*hgt*sizeof(Pixel));
	memset(pic[0], 0, wid*hgt*sizeof(Pixel));
	for(i=1;i!=hgt;i++) pic[i]=pic[i-1]+wid;
	line=emalloc(2*wid*sizeof(char));
	col=emalloc(hgt*sizeof(Acc));
	line=emalloc(hgt*sizeof(float));
	linealpha=emalloc(hgt*sizeof(float));
	stripe=emalloc(hgt*sizeof(float));
	stripealpha=emalloc(hgt*sizeof(float));
	lo=hgt;
	hi=0;
	switch(flag['t']?flag['t'][0][0]:'f'){
	case 's': i=sizeof(short);	crack=cshort;	break;
	case 'i': i=sizeof(int);	crack=cint;	break;
	case 'l': i=sizeof(long);	crack=clong;	break;
	case 'f': i=sizeof(float);	crack=cfloat;	break;
	case 'd': i=sizeof(double);	crack=cdouble;	break;
	case 'c': i=sizeof(char);	crack=cchar;	break;
	case 'u': i=sizeof(unsigned char);crack=cuchar;	break;
	default: fprint(2, "bad type %c\n", flag['t'][0][0]); exits("bad type");
	}
	/* Check and read input */
	if((f=open(argv[1], 0))<0){
		perror(argv[1]);
		exits("open input");
	}
	ndata=size(f, i);
	if(ndata<=0){
		fprint(2, "%s: data file too small (ndata=%d)\n", argv[0], ndata);
		exits("too small");
	}
	data=emalloc(ndata*sizeof(float *));
	data[0]=emalloc(ndata*ndata*sizeof(float));
	for(i=1;i!=ndata;i++) data[i]=data[i-1]+ndata;
	fixlights();
	for(i=0;i!=ndata;i++){
		(*crack)(f, data[i], ndata);
		for(j=0;j!=ndata;j++)
			if(fabs(data[i][j])>hi)hi=fabs(data[i][j]);
	}
	fprint(2, "size %dx%d, range %g\n", ndata, ndata, hi);
	/* Create output */
	outf=picopen_w("OUT", "runcode", x0, y0, wid, hgt, flag['a']?"m":"ma", argv, 0);
	if(outf==0){
		perror("OUT");
		exits("creat output");
	}
	/* Rescale input data */
	scale=flag['r']?range:hi;
	if(scale!=0.){
		for(i=0;i!=ndata;i++)for(j=0;j!=ndata;j++)
			data[i][j]*=ndata/scale;
	}
	/* Draw image */
	for(j=0;j!=wid;j++)
		adraw(j);
	/* Write output */
	for(i=0;i!=hgt;i++){
		if(flag['a']){
			for(j=0;j!=wid;j++)
				line[j]=pic[i][j].m;
		}
		else{
			for(j=0;j!=wid;j++){
				line[2*j]=pic[i][j].m;
				line[2*j+1]=pic[i][j].a;
			}
		}
		picwrite(outf, line);
	}
	exits("");
}
int size(int f, int l){
	Dir stbuf;
	register i;
	dirfstat(f, &stbuf);
	for(i=0;i*i*l<stbuf.length;i++);
	if(i*i*l!=stbuf.length)
		fprint(2, "Non-square file padded to %dx%d\n", i, i);
	return i;
}
float sample(int x, int y){
	if(x<0) x=0; else if(ndata<=x) x=ndata-1;
	if(y<0) y=0; else if(ndata<=y) y=ndata-1;
	return data[y][x];
}
double fmin(double a, double b){ return a<b?a:b; }
int min(int a, int b){ return a<b?a:b; }
double fmax(double a, double b){ return a>b?a:b; }
int max(int a, int b){ return a>b?a:b; }
int clamp(double a){
	return a<0?0:a<255?a:255;
}
float intens(int x, int y){
	struct light *l, *el=&light[nlight];
	float nx, ny, nz, v, len;
	nx=sample(x+1, y)-sample(x-1, y);
	ny=sample(x, y+1)-sample(x, y-1);
	nz=2;
	len=sqrt(nx*nx+ny*ny+nz*nz);
	nx/=len;
	ny/=len;
	nz/=len;
	v=ambient;
	for(l=light;l!=el;l++)
		v+=l->lum*(diff*fmax(nx*l->lx+ny*l->ly+nz*l->lz, 0.)
			+spec*pow(fmax(nx*l->sx+ny*l->sy+nz*l->sz, 0.), bump));
	return v;
}
typedef struct{
	int x, y, z;
	float tx, tz, pv;
}Point3d;
int inside(Point3d p){
	return(0<=p.x && p.x<ndata && 0<=p.y && p.y<ndata);
}
Point3d add3(Point3d p, Point3d q){
	p.x+=q.x;
	p.y+=q.y;
	p.z+=q.z;
	return(p);
}
Point3d Pt3(int x, int y, int z){
	Point3d p;
	p.x=x;
	p.y=y;
	p.z=z;
	return(p);
}
/*
 * samples transformed into screen coordinates
 */
void tsample(Point3d *p){
	p->tz=sample(p->x, p->y)+(p->x+p->y)/2.+hgt/4;
	p->tx=p->y-p->x+wid/2;
	p->pv=intens(p->x, p->y);
}
void flushline(void);
void accumulate(Point3d, Point3d, Point3d);
void flushstripe(int);
/*
 * draw an antialised vertical stripe of the picture in pic[.][j].
 * This stripe maps into the line y=j-wid/2+x in the data array.
 * Points [x, y, f(x, y)] map to [(x+y)/2+f(x, y)+hgt/4, y-x+wid/2]
 * In front to back order, we accumulate sub-stripes into an array that
 * gets composited into the ultimate vertical scanline and zeroed every time the
 * stripe changes direction from up to down or vice versa.
 */
#define	trap(a, b) ((b.tx-a.tx)*(a.tz+b.tz))
void adraw(int j){
	int dir=0, plo;
	Point3d p, q, r;
	p.x=0;
	p.y=j-ndata;
	tsample(&p);
	q=add3(p, Pt3(0, 1, 0));
	if(!inside(p)){
		p=add3(p, Pt3(-p.y, -p.y, 0));
		tsample(&p);
		if(!inside(p))
			return;
	}
	tsample(&q);
	if(!inside(q)){
		q=add3(q, Pt3(-q.y, -q.y, 0));
		tsample(&q);
		if(!inside(q))
			return;
	}
	for(;;){
		plo=p.x+p.y<q.x+q.y;
		if(plo)
			r=add3(p, Pt3(1, 1, 0));
		else
			r=add3(q, Pt3(1, 1, 0));
		if(!inside(r))
			break;
		tsample(&r);
		if((trap(p, q)+trap(q, r)<trap(p, r))!=dir){
			flushline();
			dir=!dir;
		}
		accumulate(p, q, r);
		if(plo)
			p=r;
		else
			q=r;
	}
	flushstripe(j);
}
void flushline(void){
	int i;
	for(i=lo;i<=hi;i++){
		if(stripe[i]<1. && fabs(col[i].a-col[i].acca)>1e-3){
			col[i].m+=col[i].accm*(1-col[i].a);
			col[i].a+=col[i].acca*(1-col[i].a);
		}
		col[i].accm=0;
		col[i].acca=0;
	}
	lo=hgt;
	hi=0;
}
void flushstripe(int j){
	int i;
	flushline();
	for(i=0;i!=hgt;i++){
		pic[hgt-1-i][j].m=clamp(col[i].m*256);
		pic[hgt-1-i][j].a=clamp(col[i].a*256);
		col[i].m=0;
		col[i].a=0;
	}
}
void accumulate(Point3d p, Point3d q, Point3d r){
	Point3d s, t;
	register float a, newa, da, y, z, dz, z0, z1, pv;
	pv=r.pv;
	/* sort so that p.tz<=q.tz<=r.tz */
	if(p.tz>q.tz){ t=p; p=q; q=t; }
	if(q.tz>r.tz){ t=q; q=r; r=t; }
	if(p.tz>q.tz){ t=p; p=q; q=t; }
	if(r.tz==p.tz)
		return;
	s.tz=q.tz;
	s.tx=p.tx+(r.tx-p.tx)*(q.tz-p.tz)/(r.tz-p.tz);
	if(q.tx>s.tx){ t=q; q=s; s=t; }
	/* p.tz<=q.tz==s.tz<=r.tz and q.tx<=s.tx */
	if(p.tz!=q.tz){
		y=(s.tx-q.tx)/(q.tz-p.tz)/2.;
		a=0.;
		z0=ceil(p.tz);
		z1=ceil(q.tz);
		for(z=z0;z<=z1;z+=1){
			dz=fmin(z, q.tz)-p.tz;
			newa=dz*dz*y;
			da=newa-a;
			if(0<z && z<=hgt){
				col[(int)z-1].accm+=da*pv;
				col[(int)z-1].acca+=da;
			}
			a=newa;
		}
	}
	if(r.tz!=q.tz){
		y=(s.tx-q.tx)/(r.tz-q.tz)/2.;
		a=0.;
		z0=floor(r.tz);
		z1=floor(q.tz);
		for(z=z0;z>=z1;z-=1){
			dz=r.tz-fmax(z, q.tz);
			newa=dz*dz*y;
			da=newa-a;
			if(0<=z && z<hgt){
				col[(int)z].accm+=da*pv;
				col[(int)z].acca+=da;
			}
			a=newa;
		}
	}
	lo=max(0, min((int)floor(p.tz), lo));
	hi=min(hgt-1, max((int)ceil(r.tz), hi));
}
void fixlights(void){
	register struct light *lp, *el=&light[nlight];
	float l;
	for(lp=light;lp!=el;lp++){
		l=sqrt(lp->lx*lp->lx+lp->ly*lp->ly+lp->lz*lp->lz);
		if(l!=0.){
			lp->lx/=l;
			lp->ly/=l;
			lp->lz/=l;
		}
		lp->sx=lp->lx+.5;
		lp->sy=lp->ly+.5;
		lp->sz=lp->lz+sqrt(.5);
		l=sqrt(lp->sx*lp->sx+lp->sy*lp->sy+lp->sz*lp->sz);
		if(l!=0.){
			lp->sx/=l;
			lp->sy/=l;
			lp->sz/=l;
		}
	}
}
void getlights(char *name){
	register struct light *lp=light;
	register FILE *fd;
	double x, y, z, lum;
	char line[512];
	if((fd=fopen(name, "r"))==NULL){
		perror(name);
		exits("open lights");
	}
	while(fgets(line, sizeof line, fd)!=NULL){
		if(sscanf(line, "light%lf%lf%lf%lf", &x, &y, &z, &lum)==4){
			if(lp==&light[NLIGHT]){
				fprint(2, "aplot: can only do %d lights\n",
					NLIGHT);
				exits("too many lights");
			}
			lp->lx=x;
			lp->ly=y;
			lp->lz=z;
			lp->lum=lum;
			lp++;
		}else if(sscanf(line, "ambient%lf", &x)==1)
			ambient=x;
		else if(sscanf(line, "spec%lf", &x)==1)
			spec=x;
		else if(sscanf(line, "diff%lf", &x)==1)
			diff=x;
		else if(sscanf(line, "bump%lf", &x)==1)
			bump=x;
		else{
			if(line[0])
				line[strlen(line-1)]='\0';
			fprint(2, "aplot: what means `%s'?\n", line);
		}
	}
	nlight=lp-light;
	fclose(fd);		
}