// fijiart
// Allows interactive manipulation of images-rendered-as-terrain
//
// Adaptation of CLIP fiji, which was adaptation of a very old
// program I wrote c. 1988.
//
// For full usage information, just type
// fijiart
// on the command line.
// Usual usage is
// fijiart inpic
//
// John Robinson 2004.

// Here's the usagetext displayed in the interaction windows
const char usagetext[]="3D CONTROL:\n\
1. Right drag in \"Fiji picture\"\n\
to shift the terrain picture and\n\
centre the desired rotation axis.\n\
2. Left drag in \"Fiji picture\" to\n\
rotate the terrain.\n\
3. Left drag in \"Framing\" to shift\n\
the camera over the tilted terrain.\n\
4. Right vertical drag in \"Framing\"\n\
to push and pull the camera.\n\
5. Right horizontal drag in \n\
\"Framing\" to roll the camera.\n\n\
EASY USAGE:\n\
In all 3 windows, left drag only.\n\n\
SAVING:\n\
Close \"Fiji picture\" to finish and\n\
save result in file \"lastfijipic\".";

const char gainzoomtext[]="Left drag in here to adjust\n\
terrain gain (surface height).\n\
Right drag to adjust focal length.";

#include "picture.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

const float LUMINANCE_MULTIPLIER = 230;
const float AMBIENT = 20.0;
/* Biggest image we will handle */
const int BIGIMAGE = 40000;

// float versions of command line arguments
float d_focallength = 200.0;
float d_camdistance = 200.0;
float d_ZGAIN = 40.0;
float d_XPLUS = 0.0;
float d_YPLUS = 0.0;
float d_XCENT = 0.0;
float d_YCENT = 0.0;

// int versions of command line parameters
int ZGAIN = 40;
int focallength = 200;
int camdistance = 200;
int camdistanceset = 0;
int output_xsiz, output_ysiz;
int XPLUS = 0;
int XPLUSset = 0;
int YPLUS = 0;
int YPLUSset = 0;
int XCENT = 0;
int YCENT = 0;
float camtheta = 0.0;

int **outbuffer;
int outimage;
int	**inbuffer;	/* pointer to buffer of input image */
int	inimage;	/* input image number */
int maxval, minval;
picture_of_int *inpic;
picture_of_int *outpic;

/* Parameters to define angle of view */
const float DEGRAD = (float)3.141592/180;
float thetax = 0.0;
float thetaz = 0.0;
float light_thetax = -40.0;
float light_thetaz = -30.0;

int perspective_flag = 1;
int texture_flag = 1;
int interactive_flag = 1;
char in_filename[64];
char texture_filename[64];


// My implementation of getopt that can run on Windows too
int optind = 1;
char *optarg = 0;
int my_getopt(int argc, char **argv, char *str) {
	if (optind == argc)
		return EOF;
	if (argv[optind][0] == '-') {
		char *p = str;
		while (*p) {
			if (*p == ':') { p++; continue; }
			if (*p == argv[optind][1]) {
				if (*(p+1) == ':') {
					optarg = argv[optind+1];
					optind += 2;
					return (int) *p;
					}
				else {
					optarg = 0; 
					optind++;
					return (int) *p;
					}
				}
			p++;
			}
		}
	return EOF;
	}

void handle_args(int argc, char *argv[])
    {
    int c;
    while ((c = my_getopt(argc, argv, "bpTg:u:v:w:l:m:d:f:x:y:X:Y:")) != EOF)
        {
        switch (c)
            {
            case 'b':
                interactive_flag = 0;
                break;
            case 'p':
                perspective_flag = 0;
                break;
            case 'u':
                camtheta = (float)atof(optarg);
                break;
            case 'v':
                thetax = (float)atof(optarg);
                break;
            case 'w':
                thetaz = (float)atof(optarg);
                break;
	    case 'x':
		d_XCENT = (float)atof(optarg);
		XCENT = -1;
		break;
	    case 'y':
		d_YCENT = (float)atof(optarg);
		YCENT = -1;
		break;
	    case 'X':
		d_XPLUS = (float)atof(optarg);
		XPLUS = -1;
		XPLUSset = 1;
		break;
	    case 'Y':
		d_YPLUS = (float)atof(optarg);
		YPLUS = -1;
		YPLUSset = 1;
		break;
	    case 'g':
		d_ZGAIN = (float)atof(optarg);
		ZGAIN = -1;
		break;
	    case 'f':
		d_focallength = (float)atof(optarg);
		focallength = -1;
		if (!camdistanceset)
			camdistance = focallength;
		// So set cameradistance to focallength unless explicitly
		// set otherwise
		break;
	    case 'd':
		d_camdistance = (float)atof(optarg);
		camdistance = -1;
		camdistanceset = 1;
		break;
            case 'l':
                light_thetax = (float)atof(optarg);
                break;
            case 'm':
                light_thetaz = (float)atof(optarg);
		break;
	    case 'T':
		texture_flag = 0;
		break;
            }
        }
    if (((c = (argc - optind)) != 1)&&(c != 2))
	{
	fprintf(stderr,"Usage: fijiart [options] depth_file [texture_file]\n");
	fprintf(stderr,"Options:\n");
	fprintf(stderr,"\t[-b] DON'T do interactively (b = batch)\n");
	fprintf(stderr,"\t[-p] DON'T render in perspective\n");
	fprintf(stderr,"\t[-T] DON'T use texture\n");
	fprintf(stderr,"[Following 7 all as proportion of picture width]\n");
	fprintf(stderr,"\t[-X X framing offset] \n");
	fprintf(stderr,"\t[-Y Y framing offset] \n");
	fprintf(stderr,"\t[-x X rotation centre] \n");
	fprintf(stderr,"\t[-y Y rotation centre] \n");
	fprintf(stderr,"\t[-f focal length (and default initial distance)]\n");
	fprintf(stderr,"\t[-d camera distance]\n");
	fprintf(stderr,"\t[-g amplitude_gain]\n");
	fprintf(stderr,"\t[-u camtheta]\n");
	fprintf(stderr,"\t[-v thetax]\n");
	fprintf(stderr,"\t[-w thetaz] \n");
	fprintf(stderr,"\t[-l light_thetax]\n"); 
	fprintf(stderr,"\t[-m light_thetaz]  \n");
	exit(1);
	}
    strcpy(in_filename,argv[optind++]);
    if (c == 2)
	strcpy(texture_filename,argv[optind]);
    else
	strcpy(texture_filename,in_filename);
    }

void create_clip_im() {
	outpic = new picture_of_int(output_ysiz, output_xsiz, 0);
	outbuffer = outpic->bufptr();
	*outpic = 0;
	}


int	nrows, ncols;	/* number of image rows and columns */
int get_image_and_range(char *infile) {
	colour_picture colinpic(infile);
	inpic = new picture_of_int(colinpic.nrows(), colinpic.ncols());
	inpic->loadascolour(colinpic);
/*
	if (colinpic.ismonochrome())
		inpic = new picture_of_int(colinpic.mono);
	else
		{
		// Stack up colours into one picture_of_int
		int minval = colinpic.min();
		if (minval < 0) {
			colinpic.red -= minval;
			colinpic.green -= minval;
			colinpic.blue -= minval;
			}
		int maxval = colinpic.max();
		if (maxval > 255) {
			colinpic.red *= 255;
			colinpic.green *= 255;
			colinpic.blue *= 255;
			colinpic.red /= maxval;
			colinpic.green /= maxval;
			colinpic.blue /= maxval;
			}
		colinpic.red *= 256*256;
		colinpic.green *= 256;
		inpic = new picture_of_int(colinpic.blue);
		*inpic += colinpic.red;
		*inpic += colinpic.green;
		}
	nrows = inpic->nrows();
	ncols = inpic->ncols();
	inbuffer = inpic->bufptr();
	maxval = inpic->max();
	minval = inpic->min();
	return colinpic.ismonochrome();
*/
	nrows = inpic->nrows();
	ncols = inpic->ncols();
	inbuffer = inpic->bufptr();
	maxval = inpic->max();
	minval = inpic->min();
	return 0;
}

void scaleup_params() {
	if (XCENT < 0)
		XCENT = (int) (d_XCENT*ncols); 
	if (YCENT < 0)
		YCENT = (int) (d_YCENT*ncols);
	if (XPLUS < 0)
		XPLUS = (int) (d_XPLUS*ncols); 
	if (YPLUS < 0)
		YPLUS = (int) (d_YPLUS*ncols);
	if (focallength < 0)
		focallength = (int) (d_focallength*ncols); 
	if (camdistance < 0)
		camdistance = (int) (d_camdistance*ncols);
	if (ZGAIN < 0)
		ZGAIN = (int) (d_ZGAIN*ncols);
	}
void scaledown_params() {
	d_XCENT = ((float) XCENT)/ncols; d_YCENT = ((float) YCENT)/ncols;
	d_XPLUS = ((float) XPLUS)/ncols; d_YPLUS = ((float) YPLUS)/ncols;
	d_focallength = ((float) focallength)/ncols; 
	d_camdistance = ((float) camdistance)/ncols;
	d_ZGAIN = ((float) ZGAIN)/ncols;
	}
	
void remove_image() {
	delete inpic;
	}


struct gridquad
	{
	int i, j;	/* Indices in array, for use after sorting */
	float depth;
	} *gridquad, *pgq;
struct grid
	{
	float coords[3];
	int overlayfunction;
	int overlayline;
	float transcoords[3];
	float luminance;
	} *grid, *pgrid;

/*
** Comparison routine for quicksorting
*/
int lowfirst(const void *p1, const void *p2){
	const struct gridquad *pg1 = (const struct gridquad *)p1;
	const struct gridquad *pg2 = (const struct gridquad *)p2;
	if (pg1->depth > pg2->depth)
		return(1);
	else if (pg2->depth > pg1->depth)
		return(-1);
	return(0);
	}



float newaxes[3][3];
/* Transform coordinates by rotation */
void transform(float locthetaz,float locthetax,int perspective)
	{
	// cout << "transform(" << thetaz << "," << thetax << "," << perspective << ")\n";
	float sinthetaz, costhetaz;
	float sinthetax, costhetax;
	float sincamtheta, coscamtheta;
	float xsines[BIGIMAGE];
	float xcosines[BIGIMAGE];
	float ysines[BIGIMAGE];
	float ycosines[BIGIMAGE];
	float newx,newy,newz;
	float outx,outy,outz;
	int i,j;
	sinthetaz = sin(locthetaz);
	costhetaz = cos(locthetaz);
	sinthetax = sin(locthetax);
	costhetax = cos(locthetax);
	sincamtheta = sin(camtheta);
	coscamtheta = cos(camtheta);
	/*
	** First do rotation about z. Use fact that x and y start life as
	** a regular grid
	*/
	pgrid = grid;
	for (i = 0; i < ncols; i++)
		{
		xsines[i] = (pgrid->coords[0]+XCENT)*sinthetaz;
		xcosines[i] = (pgrid->coords[0]+XCENT)*costhetaz;
		pgrid++;
		}
	pgrid = grid;
	for (i = 0; i < nrows; i++)
		{
		ysines[i] = (pgrid->coords[1]+YCENT)*sinthetaz;
		ycosines[i] = (pgrid->coords[1]+YCENT)*costhetaz;
		pgrid += ncols;
		}
	pgrid = grid;
	for (i = 0; i < nrows; i++)
	    {
	    for (j = 0; j < ncols; j++, pgrid++)
		{
		newx = xcosines[j] - ysines[i];
		newy = xsines[j] + ycosines[i];
		newz = pgrid->coords[2];
		outx = newx;
		outy = newz*sinthetax + newy*costhetax;
		outz = newz*costhetax - newy*sinthetax;
		if (perspective) {
			outx *= focallength/(camdistance-outz);
			outy *= focallength/(camdistance-outz);
			}
		if (camtheta) {
			pgrid->transcoords[0] = outx*coscamtheta - 
							outy*sincamtheta;
			pgrid->transcoords[1] = outx*sincamtheta +
							outy*coscamtheta;
			}
		else {
			pgrid->transcoords[0] = outx;
			pgrid->transcoords[1] = outy;
			}
		pgrid->transcoords[2] = outz;
		}
	    }
	// Calculate projection of axes onto new space
	// origin maps to origin
	for (i = 0; i < 3; i++) {
		switch (i) {
			case 0:
				// x axis: mapping of (1,0,0)
				newx = costhetaz; newy = sinthetaz; newz = 0;
				break;
			case 1:
				// y axis: mapping of (0,1,0)
				newx = sinthetaz; newy = costhetaz; newz =0;
				break;
			case 2:
				// z axis: mapping of (0,0,1)
				newx = 0; newy = 0; newz = 1;
				break;
			}
		outx = newx;
		outy = newz*sinthetax + newy*costhetax;
		outz = newz*costhetax - newy*sinthetax;
		if (perspective) {
			outx *= focallength/(camdistance-outz);
			outy *= focallength/(camdistance-outz);
			}
		if (camtheta) {
			newaxes[i][0] = outx*coscamtheta - 
							outy*sincamtheta;
			newaxes[i][1] = outx*sincamtheta +
							outy*coscamtheta;
			}
		else {
			newaxes[i][0] = outx;
			newaxes[i][1] = outy;
			}
		newaxes[i][2] = outz;
		}
/*
	printf("Projection of axes onto displayed picture:\n");
	printf("x: (%f,%f,%f)\n",newaxes[0][0],newaxes[0][1],newaxes[0][2]);
	printf("y: (%f,%f,%f)\n",newaxes[1][0],newaxes[1][1],newaxes[1][2]);
	printf("z: (%f,%f,%f)\n",newaxes[2][0],newaxes[2][1],newaxes[2][2]);
*/
	}

struct {
	int line;
	int xs, xe;
	int ls, le;
	} polyrast[BIGIMAGE];

		
int polyplot(int x,int y,int l)
	{
	int i;
	for (i = 0; i < BIGIMAGE; i++)
		{
		if (polyrast[i].line == y)
			{
			if (x < polyrast[i].xs)
				{
				polyrast[i].xs = x;
				polyrast[i].ls = l;
				}
			else if (x > polyrast[i].xe)
				{
				polyrast[i].xe = x;
				polyrast[i].le = l;
				}
			return(0);
			}
		else if (polyrast[i].line == -1)
			{
			polyrast[i].line = y;
			polyrast[i].xs = x;
			polyrast[i].xe = x;
			polyrast[i].ls = l;
			polyrast[i].le = l;
			return(0);
			}
		}
	return(-1);
	}
void xyplot(int **buf, int xs, int xe, int ys, int ye, int ls, int le) {
	register int x = (xs + xe)>>1;
	register int l;
	register int y = (ys + ye)>>1;
	if ((y == ys)||(y == ye))
		return;
	if (ls == le)
		l = ls;
	else 
		// l = (ls>>1) + (le>>1);	// Works because 7 bit precision
					// means bottom bits are blanked
		l = (ls+le)>>1;
	buf[y][x] = l;
	l &= 0xfefefefeL;
	polyplot(x,y,l);
	xyplot(buf, xs, x, ys, y, ls, l);
	xyplot(buf, x, xe, y, ye, l, le);
	}
		
void xplot(int *buf, int xs, int xe, int ls, int le) {
	if (ls == le) {
		register int *p = buf + xs + 1;
		register int cnt = xe - xs - 1;
		while(cnt--)
			*p++ = ls;
		return;
		}
	int x = (xs + xe)>>1;
	// int l = (ls>>1) + (le>>1);	// Works because 7 bit precision
					// means bottom bits are blanked
	int l = (ls+le)>>1;
	buf[x] = l;
	l &= 0xfefefefeL;		// Continue 7 bit precision
	if (x > xs+1)
		xplot(buf, xs, x, ls, l);
	if (xe > x+1)
		xplot(buf, x, xe, l, le);
	}

void polygoninit() {
	for (int i = 0; i < BIGIMAGE; i++)
		polyrast[i].line = -1;	
	}

void polygon (picture_of_int& pic, int *x, int *y, int *l)
	{
	int i;
	int ll[4];
	// If all of polygon is outside limits, don't draw
	int numbad = 0;
	for (i = 0; i < 4; i++) {
		int thisbad = 0;
		if (x[i] < 0) { x[i] = 0; thisbad = 1; }
		else if (x[i] >= pic.ncols()) { x[i] = pic.ncols()-1; thisbad = 1;}
		if (y[i] < 0) { y[i] = 0; thisbad = 1; }
		else if (y[i] >= pic.nrows()) { y[i] = pic.nrows()-1; thisbad = 1;}
		numbad += thisbad;
		}
	if (numbad == 4)
		return;
	for (i = 0; i < BIGIMAGE; i++) {
		if (polyrast[i].line < 0)
			break;
		polyrast[i].line = -1;	
		}
	for (i = 0; i < 4; i++)
		ll[i] = l[i]&0xfefefefeL;	// All done to 7 bit precision
	for (i = 0; i < 4; i++) {
		pic[y[i]][x[i]] = ll[i];
		polyplot(x[i],y[i],ll[i]);
		}
	xyplot(pic.bufptr(),x[0],x[1],y[0],y[1],ll[0],ll[1]);
	xyplot(pic.bufptr(),x[1],x[2],y[1],y[2],ll[1],ll[2]);
	xyplot(pic.bufptr(),x[2],x[3],y[2],y[3],ll[2],ll[3]);
	xyplot(pic.bufptr(),x[3],x[0],y[3],y[0],ll[3],ll[0]);
	for (i = 0; i < BIGIMAGE; i++) {
		if (polyrast[i].line < 0)
			break;
		if (polyrast[i].xe > polyrast[i].xs + 1)
			xplot(pic.rowptr(polyrast[i].line),
				polyrast[i].xs, polyrast[i].xe,
				polyrast[i].ls, polyrast[i].le);
		}
	}
			
void get_normal(float *c1,float *c2,float *c3,float *c4,float *result)
	{
	float nx,ny,nz;
	float ux,uy,uz;
	float vx,vy,vz;
	float modulus;
	ux = c1[0] - c3[0];
	uy = c1[1] - c3[1];
	uz = c1[2] - c3[2];

	vx = c2[0] - c4[0];
	vy = c2[1] - c4[1];
	vz = c2[2] - c4[2];

	/* Cross product to find normal */
	nx = uy*vz - uz*vy;
	ny = - (ux*vz - uz*vx);
	nz = ux*vy - uy*vx;
	modulus = sqrt(nx*nx + ny*ny + nz*nz);
	if (nz < 0)
		nz = -nz;
	*result = AMBIENT + (LUMINANCE_MULTIPLIER*nz)/modulus;
	}

int main(int argc,char *argv[])
	{
	int i, j;
	int texture_is_monochrome = 1;

    	handle_args(argc,argv);
    	thetax *= DEGRAD;
    	thetaz *= DEGRAD;
	camtheta *= DEGRAD;
    	light_thetax *= DEGRAD;
    	light_thetaz *= DEGRAD;
    	get_image_and_range(in_filename);
	scaleup_params();
    	grid = new struct grid[nrows*ncols];
    	gridquad = new struct gridquad[(nrows-1)*(ncols-1)];
	polygoninit();

/* Generate values for z at grid points */

	pgrid = grid;
	for (i = 0; i < nrows; i++)
		{
		for (j = 0; j < ncols; j++)
			{
			pgrid->coords[0] =  (j - (float)((ncols-1)/2.0));
			pgrid->coords[1] =  (i - (float)((nrows-1)/2.0));
     pgrid->coords[2] = ((float)(inbuffer[i][j]-minval)*ZGAIN)/(maxval-minval);
			pgrid++;
			}
		}
	remove_image( );
	output_xsiz = ncols;
	output_ysiz = nrows;
	if (texture_flag)
		texture_is_monochrome = get_image_and_range(texture_filename);
	else	
	/* Deal with lighting */
	/*  transform all grid points */
		transform(light_thetaz,light_thetax,0);
	/* Now work out normals */
	pgrid = grid;
	for (i = 0; i < nrows; i++)
		{
		for (j = 0; j < ncols; j++)
			{
			if ((i==0)||(j==0)||(i==nrows-1)||(j==ncols-1))
				{
				pgrid->luminance = 0.0;
				pgrid++;
				continue;
				}
			if (texture_flag) 
				pgrid->luminance = (float) inbuffer[i][j];
			else
				get_normal((pgrid-1)->transcoords,
				(pgrid-ncols)->transcoords,
				(pgrid+1)->transcoords,
				(pgrid+ncols)->transcoords,
				&pgrid->luminance);
			pgrid++;
			}
		}
	create_clip_im(); // So this sets size of output pic
	if (!XPLUSset)
		XPLUS = ncols/2;
	if (!YPLUSset)
		YPLUS = nrows/2;
	if (interactive_flag)
		outpic->show("Fiji picture");
	picture_of_int frame(256,256);
	picture_of_int *framingpic = &frame;
	picture_of_int textpic(usagetext,255,0);
	frame.paste_centred_on(textpic,128,128);
	picture_of_int gainzoom(40,256);
	picture_of_int *gainzoompic = &gainzoom;
	picture_of_int textpic2(gainzoomtext,255,0);
	gainzoom.paste_centred_on(textpic2,128,20);
	if (interactive_flag) {
		framingpic->show("Framing");
		gainzoompic->show("L:gain, R:zoom");
		}
	
	int pointerx = -1;
	int pointery = -1;
	int thistime = 1;
	int pointingleft = 0;
	int pointingright = 0;
	int tophalf = 0;	// For left point and x movement
	while(1) {
		if (!thistime) {
			if (framingpic->pointleftbutton()) {
			    if (pointingleft && ((framingpic->pointx() != pointerx)||(framingpic->pointy() != pointery))) {
				int pointxoffset = pointerx - framingpic->pointx();
				int pointyoffset = pointery - framingpic->pointy();
				pointerx = framingpic->pointx();
				pointery = framingpic->pointy();
				XPLUS -= pointxoffset;
				YPLUS -= pointyoffset;
				thistime = 2;
				}
			    else if (!pointingleft) {
				pointingleft = 1;
				pointerx = framingpic->pointx();
				pointery = framingpic->pointy();
				}
			    }
			else if (framingpic->pointrightbutton()) {
			    if (pointingright && ((framingpic->pointx() != pointerx)||(framingpic->pointy() != pointery))) {
				int pointxoffset = pointerx - framingpic->pointx();
				int pointyoffset = pointery - framingpic->pointy();
				pointerx = framingpic->pointx();
				pointery = framingpic->pointy();
				camtheta += pointxoffset*180*DEGRAD/(output_xsiz);
				camdistance -= pointyoffset*200/output_ysiz;
				thistime = 1;
				}
			    else if (!pointingright) {
				pointingright = 1;
				pointerx = framingpic->pointx();
				pointery = framingpic->pointy();
				}
			    }
			else if (outpic->pointrightbutton()) {
			    if (pointingright && ((outpic->pointx() != pointerx)||(outpic->pointy() != pointery))) {
				int pointxoffset = pointerx - outpic->pointx();
				int pointyoffset = pointery - outpic->pointy();
				pointerx = outpic->pointx();
				pointery = outpic->pointy();
				XCENT -= (int)(pointxoffset*newaxes[0][0]
						+ pointyoffset*newaxes[1][0]);
				YCENT -= (int)(-pointxoffset*newaxes[0][1]
						+ pointyoffset*newaxes[1][1]);
				camdistance -=(int)( 
				  (pointyoffset*newaxes[2][1]*200)/output_ysiz);
				thistime = 1;
				}
			    else if (!pointingright) {
				pointingright = 1;
				pointerx = outpic->pointx();
				pointery = outpic->pointy();
				}
			    }
			else if (outpic->pointleftbutton()) {
			    if (pointingleft && ((outpic->pointx() != pointerx)||(outpic->pointy() != pointery))) {
				int pointxoffset = pointerx - outpic->pointx();
				int pointyoffset = pointery - outpic->pointy();
				pointerx = outpic->pointx();
				pointery = outpic->pointy();
				thetax -= pointyoffset*180*DEGRAD/(output_ysiz);
				if (tophalf)
				  thetaz-=pointxoffset*360*DEGRAD/(output_xsiz);
				else 
				  thetaz+=pointxoffset*360*DEGRAD/(output_xsiz);
				thistime = 1;
				}
			    else if (!pointingleft) {
				pointingleft = 1;
				pointerx = outpic->pointx();
				pointery = outpic->pointy();
				if (pointery < output_ysiz/2)
					tophalf = 1;
				else
					tophalf = 0;
				}
			    }
			else if (gainzoompic->pointrightbutton()) {
			    if (pointingright && ((gainzoompic->pointx() != pointerx)||(gainzoompic->pointy() != pointery))) {
				int pointxoffset = pointerx - gainzoompic->pointx();
				int pointyoffset = pointery - gainzoompic->pointy();
				pointerx = gainzoompic->pointx();
				pointery = gainzoompic->pointy();
				focallength -= pointxoffset*200/output_xsiz;
				thistime = 1;
				}
			    else if (!pointingright) {
				pointingright = 1;
				pointerx = gainzoompic->pointx();
				pointery = gainzoompic->pointy();
				}
			    }
			else if (gainzoompic->pointleftbutton()) {
			    if (pointingleft && ((gainzoompic->pointx() != pointerx)||(gainzoompic->pointy() != pointery))) {
				int pointxoffset = pointerx - gainzoompic->pointx();
				int pointyoffset = pointery - gainzoompic->pointy();
				pointerx = gainzoompic->pointx();
				pointery = gainzoompic->pointy();
				if (pointxoffset) {
					ZGAIN -= (int) (pointxoffset*200/(output_xsiz));
					pgrid = grid;
					for (i = 0; i < nrows; i++)
						{
						for (j = 0; j < ncols; j++)
							{
							pgrid->coords[0] =  (j - (float)((ncols-1)/2.0));
							pgrid->coords[1] =  (i - (float)((nrows-1)/2.0));
     				pgrid->coords[2] = ((float)(inbuffer[i][j]-minval)*ZGAIN)/(maxval-minval);
							pgrid++;
							}
						}
					}
				thistime = 1;
				}
			    else if (!pointingleft) {
				pointingleft = 1;
				pointerx = gainzoompic->pointx();
				pointery = gainzoompic->pointy();
				}
			    }
			else {
				pointingleft = 0;
				pointingright = 0;
				}
			}
		if ((!thistime)&&(!outpic->closerequested())&&interactive_flag)
			continue;	// Don't refresh
		thistime--;
		if (!thistime) { 	// A change that requires retransform
		    /* Retransform all grid points */
		    transform(thetaz,thetax,perspective_flag);
		    /* Put into gridquad arrays, then order, lowest first */
		    pgq = gridquad;
		    for (i = 0; i < nrows-1; i++)
	    	    {
	    	    for (j = 0; j < ncols-1; j++)
			{
			pgrid = grid + i*ncols + j;
			pgq->i = i;
			pgq->j = j;
			pgq->depth = pgrid->transcoords[2] +
					(pgrid+1)->transcoords[2] +
					(pgrid+ncols+1)->transcoords[2] +
					(pgrid+ncols)->transcoords[2];
			pgq->depth /= 4.0;
			pgq++;
			}
	    	    }
		    qsort(gridquad,(ncols-1)*(nrows-1),sizeof(struct
			gridquad), lowfirst);
		    }
		else
		    thistime--;
		/* Now have ordered list of quads. Simply draw them in order */
		pgq = gridquad;
		for (i = 0; i < nrows-1; i++)
	    	{
	    	for (j = 0; j < ncols-1; j++)
			{
			int l[4];	/* Luminances at the vertices */
			int vx[4];
			int vy[4];
			struct grid *p;
			p = grid + (ncols*pgq->i) + pgq->j;
			l[0] = (int) p->luminance;
			vx[0] = (int)(p->transcoords[0] + XPLUS);
			vy[0] = (int)(p->transcoords[1] + YPLUS);
			p++;
			l[1] = (int) p->luminance;
			vx[1] = (int)(p->transcoords[0] + XPLUS);
			vy[1] = (int)(p->transcoords[1] + YPLUS);
			p += ncols;
			l[2] = (int) p->luminance;
			vx[2] = (int)(p->transcoords[0] + XPLUS);
			vy[2] = (int)(p->transcoords[1] + YPLUS);
			p--;
			l[3] = (int) p->luminance;
			vx[3] = (int)(p->transcoords[0] + XPLUS);
			vy[3] = (int)(p->transcoords[1] + YPLUS);
			polygon(*outpic,vx,vy,l);
			pgq++;
			}
	    	}
		if (outpic->closerequested()||!interactive_flag) {
			colour_picture writeit(nrows,ncols);
			for (i = 0; i < nrows; i++)
			for (j = 0; j < ncols; j++) {
				writeit.red[i][j] = ((*outpic)[i][j]>>16)&0xff;
				writeit.green[i][j] = ((*outpic)[i][j]>>8)&0xff;
				writeit.blue[i][j] = (*outpic)[i][j]&0xff;
				}
			writeit.write("lastfijipic");
			scaledown_params();
			printf("Output written to \"lastfijipic\"\n");
			printf("Pastable commandline:\nfijiart -b ");
			printf("-f %f ",d_focallength);
			printf("-d %f ",d_camdistance);
			printf("-v %f -w %f -u %f ",
				thetax/DEGRAD,thetaz/DEGRAD,camtheta/DEGRAD); 
			printf("-X %f -Y %f ",
				d_XPLUS, d_YPLUS);
			printf("-x  %f -y %f ",
				d_XCENT, d_YCENT);
			printf("-g %f ", d_ZGAIN);
			printf("%s\n", in_filename);
			break;
			}
		if (texture_is_monochrome)
			outpic->reshow();
		else
			outpic->reshowascolour();
		*outpic = 0;
		}
		return 0;
	}