/*
circ3: test 3-D dot-drawing speed on SGI with Opengl.  Uses z-buffering.

compile with:
    cc -O2 -mips2 -o circ3 circ3.c -lGLU -lGL -lX11 -lm

NOTE: behavior differences for antialiased points (style 'P'):
    on indy, P same as indigo2, but circles no bigger than 10 pixel diameter
    on crimson RE, P draws black squares around antialiased (aa) circles
    on indigo2 max impact, P draws thin black outlines around aa circles

Paul Heckbert	21 May 1996
ph@cs.cmu.edu
*/

#include <assert.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <GL/gl.h>	/* Opengl graphics library subroutine prototypes */
#include <GL/glx.h>	/* includes Xlib.h and other X stuff too */

#define UNIT_CIRCLE 1	/* SGI display list number for unit circle */
#define SPOT 2		/* SGI display list number for spot bitmap */

#define RAD_TO_DEG(x) ((x)*(180./M_PI))
#define FRAND() ((rand()&32767)/32767.)
#define RFRAC .4

#ifndef sgi
    #define sqrtf sqrt
    #define cosf cos
    #define sinf sin
    #define atan2f atan2
    #define asinf asin
#endif

typedef struct {float x, y, z;} vec3;

static int m, nsides, style;
static float r;
static Display *display;
static Window window;		/* window number */

#define VEC_SET(xx, yy, zz, a) ((a)->x = xx, (a)->y = yy, (a)->z = zz)
#define VEC_LEN(a) sqrtf((a)->x*(a)->x+(a)->y*(a)->y+(a)->z*(a)->z)

void vec_cross(vec3 *a, vec3 *b, vec3 *c) {
     c->x = a->y * b->z - b->y * a->z;
     c->y = a->z * b->x - b->z * a->x;
     c->z = a->x * b->y - b->x * a->y;
}

void vec_smul(float s, vec3 *a, vec3 *b) {
     b->x = s*a->x;
     b->y = s*a->y;
     b->z = s*a->z;
}


void simple_circle(vec3 *cen, vec3 *n, float r, int nsides) {
    /* doesn't assume n is a unit vector */

    int i;
    float ang, cosang, sinang;
    vec3 u, v, pt;

    if (n->x*n->x < .33) VEC_SET(0, n->z, -n->y, &u);	/* n X (1,0,0) */
    else VEC_SET(-n->z, 0, n->x, &u);			/* n X (0,1,0) */
    vec_smul(r/VEC_LEN(&u), &u, &u);
    vec_cross(n, &u, &v);			/* v is another tangent vector*/
    vec_smul(r/VEC_LEN(&v), &v, &v);
    /* u, v, and n are mutually orthogonal, u and v have length r */

    /* draw a disk */
    glBegin(GL_POLYGON);
    for (i=0; i<nsides; i++) {
	ang = 2.*M_PI*i/nsides;
	cosang = cosf(ang);
	sinang = sinf(ang);
	pt.x = cen->x + cosang*u.x + sinang*v.x;
	pt.y = cen->y + cosang*u.y + sinang*v.y;
	pt.z = cen->z + cosang*u.z + sinang*v.z;
	glVertex3fv(&pt.x);
    }
    glEnd();
}

#define NMAX 100
typedef struct {float cos, sin;} circle;
static circle circtab[NMAX];		/* lookup table for circle-drawing */

void table_init(int n) {
    int i;
    float ang;

    assert(nsides<=NMAX);
    nsides = n;
    for (i=0; i<nsides; i++) {
	ang = 2.*M_PI*i/nsides;
	circtab[i].cos = cosf(ang);
	circtab[i].sin = sinf(ang);
    }
}

void table_circle(vec3 *cen, vec3 *n, float r) {
    /* doesn't assume n is a unit vector */
    register int i;
    register circle *tab;
    vec3 u, v, pt;

    if (n->x*n->x < .33) VEC_SET(0, n->z, -n->y, &u);	/* n X (1,0,0) */
    else VEC_SET(-n->z, 0, n->x, &u);			/* n X (0,1,0) */
    vec_smul(r/VEC_LEN(&u), &u, &u);
    vec_cross(n, &u, &v);			/* v is another tangent vector*/
    vec_smul(r/VEC_LEN(&v), &v, &v);
    /* u, v, and n are mutually orthogonal, u and v have length r */

    glBegin(GL_POLYGON);
    for (tab=circtab, i=nsides; --i>=0; tab++) {
	pt.x = cen->x + tab->cos*u.x + tab->sin*v.x;
	pt.y = cen->y + tab->cos*u.y + tab->sin*v.y;
	pt.z = cen->z + tab->cos*u.z + tab->sin*v.z;
	glVertex3fv(&pt.x);
    }
    glEnd();
}

void list_init(int nsides) {
    int i;
    float ang;
    vec3 pt;

    glNewList(UNIT_CIRCLE, GL_COMPILE);
    /* draw a unit radius circle in xy plane centered on origin */
    glBegin(GL_POLYGON);
    for (i=0; i<nsides; i++) {
	ang = 2.*M_PI*i/nsides;
	pt.x = cosf(ang);
	pt.y = sinf(ang);
	glVertex2fv(&pt.x);
    }
    glEnd();
    glEndList();
}

void list_circle(vec3 *cen, vec3 *n, float r) {
    /* assumes n is a unit vector */
    float tx, ty;

    tx = -asinf(n->y);
    ty = atan2f(n->x, n->z);
    glPushMatrix();
    glTranslatef(cen->x, cen->y, cen->z);
    glRotatef(RAD_TO_DEG(ty), 0., 1., 0.);
    glRotatef(RAD_TO_DEG(tx), 1., 0., 0.);
    glScalef(r, r, r);
    glCallList(UNIT_CIRCLE);
    glPopMatrix();
}

void line_init(int nsides) {
    int i;
    float ang;
    vec3 pt;

    glNewList(UNIT_CIRCLE, GL_COMPILE);
    /* draw a unit radius circle in xy plane centered on origin */
    glBegin(GL_LINE_LOOP);
    for (i=0; i<nsides; i++) {
	ang = 2.*M_PI*i/nsides;
	pt.x = cosf(ang);
	pt.y = sinf(ang);
	glVertex2fv(&pt.x);
    }
    glEnd();
    glEndList();
}

void anti_point_init(float diam) {
    /* this doesn't work well with z-buffering */
    float diamrange[2];

    glGetFloatv(GL_POINT_SIZE_RANGE, diamrange);
    if (diam<diamrange[0] || diam>diamrange[1])
	printf("warning: point diameter range=[%g,%g], you asked for %g\n",
	    diamrange[0], diamrange[1], diam);

    glEnable(GL_POINT_SMOOTH);		/* turn on antialiasing of points */
    glHint(GL_POINT_SMOOTH_HINT, GL_NICEST);
    /* glHint(GL_POINT_SMOOTH_HINT, GL_FASTEST); */
    glEnable(GL_BLEND);			/* turn on blending */
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);	/* weighted avg blend */
    glPointSize(diam);			/* set point diameter */
}

#define FONT_INDEX 1

void font_init(float radius) {
    /* create a font bitmap for a circle */
    /* (this technique for fast dot-drawing suggested by Thad Beier) */

    int x, y, r = radius, d = 2*r+1, nw = (d+7)/8;
    GLubyte *bitmap;
    
    assert(bitmap = (GLubyte *)calloc(nw*d, sizeof bitmap[0]));
	/* calloc zeros the array */
#   define SETBIT(x, y) bitmap[(y)*nw+(x)/8] |= 1<<(~(x)&7)
    for (y= -r; y<=r; y++)
	for (x= -r; x<=r; x++)
	    if (x*x+y*y<=radius*radius) SETBIT(x+r, y+r);
    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
    glNewList(SPOT, GL_COMPILE);
    glBitmap(/*width&height*/ d, d, /*origin*/ -r, -r,
	/*increment for next char*/ d, 0., bitmap);
    glEndList();
    free(bitmap);
}

void font_circle(vec3 *cen) {
    glRasterPos3f(cen->x, cen->y, cen->z);
    glCallList(SPOT);
}

void redisplay() {
    static int wild = 0;
    int lat, lon;
    float nxy, theta, phi;
    vec3 n;

    glClearColor(0., 0., 0., 0.);
    glClear(GL_DEPTH_BUFFER_BIT | GL_COLOR_BUFFER_BIT);
    wild = !wild;
    if (!wild) glColor3f(1., 1., 1.);

    switch (style) {
	case 'p': case 'P': glBegin(GL_POINTS); break;
    }
    for (lat=0; lat<=m; lat++) {
	theta = M_PI*lat/m;
	n.z = cos(theta);
	nxy = sin(theta);
	for (lon=0; lon<2*m; lon++) {
	    phi = 2.*M_PI*lon/(2*m);
	    n.x = sin(phi)*nxy;
	    n.y = cos(phi)*nxy;
	    if (wild) glColor3f(FRAND(), FRAND(), FRAND());

	    switch (style) {
		case 's': simple_circle(&n, &n, r, nsides); break;
		case 't': table_circle(&n, &n, r); break;
		case 'd':
		case 'l': list_circle(&n, &n, r); break;
		case 'p':
		case 'P': glVertex3fv(&n.x); break;
		case 'f': font_circle(&n); break;
	    }
	}
    }
    switch (style) {
	case 'p': case 'P': glEnd(); break;
    }
    glXSwapBuffers(display, window);
}

/* Attributes for a normal, 24-bit, z-buffered and double-buffered window */
int gl_double_attrs[] =
    {GLX_RGBA, GLX_DEPTH_SIZE, 12, GLX_DOUBLEBUFFER, None};

static int is_map_notify(Display *display, XEvent *event, char *arg)
{
    return event->type==MapNotify && event->xmap.window==(Window)arg;
}

void gl_create_window(int x, int y, int width, int height) {
    /* Create a sub-window of the parent window, for Opengl use */

    XSetWindowAttributes winattr;
    XEvent event;
    XVisualInfo *visinfo;
    Colormap colormap;		/* colormap for window */
    GLXContext context;

    display = XOpenDisplay(0);

    /* decide between a single-buffered or double-buffered "Visual" */
    visinfo = glXChooseVisual(display, DefaultScreen(display), gl_double_attrs);
    if (!visinfo) {
	fprintf(stderr, "gl_create_window: Can't create window with the requested attributes\n");
	exit(1);
    }

    colormap = XCreateColormap(display,
	RootWindow(display, visinfo->screen),
	visinfo->visual, AllocNone);

    /*
     * note: if you try to create a window with a 24 bit visual when the
     * default visual is 8 bits (or more precisely, when the parent window
     * has a different visual or depth) and you don't explicitly set the
     * colormap and border, you will die with a BadMatch error.
     * You can thank the brain-damaged designers of X Windows for this.
     */
    winattr.colormap = colormap;
    winattr.border_pixel = 0;
    winattr.event_mask = StructureNotifyMask;

    /* create the subwindow */
    window = XCreateWindow(display,
	/*parent*/ RootWindow(display, DefaultScreen(display)),
	/*origin*/ x, y, /*width,height*/ width, height,
	/*borderwidth*/ 0, /*depth*/ visinfo->depth, /*class*/ InputOutput,
	visinfo->visual,
	/*window attribute mask*/ CWBorderPixel | CWColormap | CWEventMask,
	&winattr);

    if (x>=0) {
	XSizeHints size;
	size.x = x;
	size.y = y;
	size.flags = USPosition;
	/* window manager obeys USPosition requests, but not PPosition */
	XSetNormalHints(display, window, &size);
    }
    /* tell window manager to draw the window */
    XMapWindow(display, window);

    /* Wait for MapNotify event when the window is first drawn */
    XIfEvent(display, &event, is_map_notify, (char *)window);

    /* Create a graphics context if it hasn't been done before */
    /* Bind the X window to an OpenGL context */
    context = glXCreateContext(display, visinfo,
	/*context share list*/ 0, /*direct rendering?*/ GL_TRUE);
    assert(context);
    assert(glXMakeCurrent(display, window, context));

    glEnable(GL_DEPTH_TEST);	/* turn on z-buffering */
}

static char usage[] = "\n\
usage: circ3 [NSIDES [STYLE [NREPEAT [NCIRCLES]]]]\n\
default: circ3 10 t 50 14\n\
styles:\n\
	s	n-sided polygons, simple (slow method)\n\
	t	n-sided polygons, using table lookup\n\
	d	n-sided polygons, in display list\n\
	l	line drawing, in display list\n\
	p	aliased points\n\
	P	antialiased points\n\
	f	font (bitmap)\n\
note: in styles p,P,f, NSIDES is actually spot radius\n\
";

main(int ac, char **av) {
    static char opts[] = "stdlpPf";
    int rep, i;
    float ns;

    if (ac==2 && av[1][0]=='-') {
	fprintf(stderr, usage);
	fprintf(stderr, "style is one of [%s]\n", opts);
	exit(1);
    }
    ns = ac>1 ? atof(av[1]) : 10;	/* number of sides per circle */
    style = ac>2 ? av[2][0] : 't';	/* display style */
    rep = ac>3 ? atoi(av[3]) : 50;	/* number of times to redisplay */
    m = ac>4 ? atoi(av[4]) : 14;	/* number of circles vertically */

    nsides = (int)ns;
    r = .5*M_PI/m*RFRAC;
    if (!strchr(opts, style)) {
	printf("style must be [%s]\n", opts);
	exit(1);
    }

    gl_create_window(200, 0, 1024, 1024);
    switch (style) {
	case 't': table_init(nsides); break;
	case 'd': list_init(nsides); break;
	case 'l': line_init(nsides); break;
	    /* kludge: ns is used as radius for point & font modes */
	case 'p': glPointSize(ns*2.); break;
	case 'P': anti_point_init(ns*2.); break;
	case 'f': font_init(ns); break;
    }

    glMatrixMode(GL_PROJECTION);
    gluPerspective(40., 1., .1, 20.);
    glMatrixMode(GL_MODELVIEW);
    glTranslatef(0., 0., -4.);
    glRotatef(-80., 1., 0., 0.);
    for (i=0; i<rep; i++) {
	redisplay();
	glRotatef(1., 0., 0., 1.);
    }
}