Newsgroups: comp.sys.sgi.graphics
From: ph+@cs.cmu.edu (Paul Heckbert)
Subject: drawing dots fast (C code)
Organization: School of Computer Science, Carnegie Mellon
Date: Fri, 8 Jul 1994 20:11:18 GMT

I've spent quite a bit of time figuring out how to draw dots (circles)
in 3-D quickly on Silicon Graphics machines.
For the conclusions of my experiments, see the comments just below.
To try things for yourself, simply compile and run.

Paul Heckbert						    ph@cs.cmu.edu
Computer Science Dept., Carnegie Mellon University
5000 Forbes Ave, Pittsburgh PA 15213-3891, USA
World Wide Web:  http://www.cs.cmu.edu:8001/afs/cs/user/ph/www/heckbert.html



# to unpack, cut here and run the following shell archive through sh
# contents: circ2.c demo.csh
#
echo extracting circ2.c
sed 's/^_//' <<'EOF11868' >circ2.c
_/*
_circ2: test 3-D dot-drawing speed on SGI.  Uses z-buffering.
_
_compile with:
_    cc -O2 -o circ2 circ2.c -lgl_s -lm
_
_tried several drawing modes:
_   c  circf - REALLY SLOW, don't use it if speed matters
_	(always draws an 80-sided polygon!)
_   s  simple polygonal method, unoptimized - faster than 'c'
_   t  polygonal method with table lookup optimization - quite fast **
_   o  callobj method - also fast, but less flexible than 't',
_	and no faster except when nsides huge
_   l  line drawing - no faster than 't' or 'o'
_   p  points - squares in screen space, fast, but a bit ugly
_	("in screen space" means the squares always face you)
_	This method and method 'P' don't work on low end machines,
_	where pntsizef() does little or nothing.
_   P  antialiased points - don't work well in combination with z-buffering
_   f  fonts - fastest way to do dots in screen space **
_
_winning methods are marked with "**"
_
_CONCLUSIONS:
_
_circf is ridiculously slow (11 to 29 times slower than methods 't' or 'o').
_Whoever wrote circf should be spanked!
_
_With methods 't' and 'o', there's not much speed difference between 3-sided
_and 10-sided polygons, but beyond nsides=10, it starts to slow down.
_The speed of polygon-drawing is an irregular function of the number of sides
_(it is constant for a while, then jumps).  The locations of the jumps
_in the time graph changes with the size of the polygon in screen space.
_
_I've done most of my testing on a Reality Engine, but most of the above
_holds on Indigo2 Extremes and Personal Irises also.
_
_A good rule of thumb for choosing the number of sides to your polygons for
_methods 's' and 't' is that the number of sides required when displaying
_a regular polygon of radius r pixels to make it look like a circle is
_proportional to sqrt(r).  Use nsides=6 for the smallest circles and
_about nsides=90 for the largest circles (radius 512 pixels).
_
_To draw spheres, check out "man libsphere".
_
_Paul Heckbert	6 July 1994
_ph@cs.cmu.edu
_*/
_
_#include <assert.h>
_#include <string.h>
_#include <math.h>
_#include <gl.h>		/* iris graphics library subroutine prototypes */
_
_#define UNIT_CIRCLE 1	/* SGI object number */
_#define RAD_TO_DEG(x) ((x)*(180./M_PI))
_
_#define RFRAC .4
_
_typedef struct {float x, y, z;} vec3;
_
_static int m, nsides, style;
_static float r;
_
_#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 circf_circle(vec3 *cen, vec3 *n, float r) {
_    /* cen is center, n is normal, r is radius */
_    /* assumes n is a unit vector */
_
_    float tx, ty;
_
_    tx = -asinf(n->y);
_    ty = atan2f(n->x, n->z);
_    pushmatrix();
_    translate(cen->x, cen->y, cen->z);
_    rot(RAD_TO_DEG(ty), 'y');
_    rot(RAD_TO_DEG(tx), 'x');
_    circf(0., 0., r);
_    popmatrix();
_}
_
_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 */
_    bgnpolygon();
_    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;
_	v3f(&pt.x);
_    }
_    endpolygon();
_}
_
_#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 */
_
_    bgnpolygon();
_    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;
_	v3f(&pt.x);
_    }
_    endpolygon();
_}
_
_void object_init(int nsides) {
_    int i;
_    float ang;
_    vec3 pt;
_
_    makeobj(UNIT_CIRCLE);
_    /* draw a unit radius circle in xy plane centered on origin */
_    bgnpolygon();
_    for (i=0; i<nsides; i++) {
_	ang = 2.*M_PI*i/nsides;
_	pt.x = cosf(ang);
_	pt.y = sinf(ang);
_	v2f(&pt.x);
_    }
_    endpolygon();
_    closeobj();
_}
_
_void object_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);
_    pushmatrix();
_    translate(cen->x, cen->y, cen->z);
_    rot(RAD_TO_DEG(ty), 'y');
_    rot(RAD_TO_DEG(tx), 'x');
_    scale(r, r, r);
_    callobj(UNIT_CIRCLE);
_    popmatrix();
_}
_
_void line_init(int nsides) {
_    int i;
_    float ang;
_    vec3 pt;
_
_    makeobj(UNIT_CIRCLE);
_    /* draw a unit radius circle in xy plane centered on origin */
_    bgnclosedline();
_    for (i=0; i<nsides; i++) {
_	ang = 2.*M_PI*i/nsides;
_	pt.x = cosf(ang);
_	pt.y = sinf(ang);
_	v2f(&pt.x);
_    }
_    endclosedline();
_    closeobj();
_}
_
_void anti_point_init(float diam) {
_    /* this doesn't work well with z-buffering */
_
_    subpixel(TRUE);			/* don't round points' coordinates */
_    /* pntsmooth(SMP_ON);		/ draw antialiased squares */
_    pntsmooth(SMP_ON|SMP_SMOOTHER);	/* draw antialiased circles */
_    blendfunction(BF_SA, BF_MSA);	/* do weighted average blending */
_    /* blendfunction(BF_SA, BF_ONE);	/ accumulation to saturation */
_    pntsizef(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+15)/16;
_    unsigned short *bitmap;
_    Fontchar chars[2];
_	/* \000 is unused (because it ends string) */
_	/* \001 is circle */
_    
_    assert(bitmap = (unsigned short *)calloc(nw*d, sizeof bitmap[0]));
_	/* calloc zeros the array */
_#   define SETBIT(x, y) bitmap[(y)*nw+(x)/16] |= 1<<(~(x)&15)
_    for (y= -r; y<=r; y++)
_	for (x= -r; x<=r; x++)
_	    if (x*x+y*y<=radius*radius) SETBIT(x+r, y+r);
_    chars[1].offset = 0;
_    chars[1].w = d;
_    chars[1].h = d;
_    chars[1].xoff = -r;
_    chars[1].yoff = -r;
_    chars[1].width = d+1;
_    defrasterfont(FONT_INDEX, /*height*/ d,
_	sizeof chars/sizeof chars[0], /*character info*/ chars,
_	/*length of bitmap array*/ nw*d, bitmap);
_    free(bitmap);
_    font(FONT_INDEX);
_}
_
_void font_circle(vec3 *cen) {
_    cmov(cen->x, cen->y, cen->z);
_    charstr("\001");
_}
_
_void redisplay() {
_    static int wild = 0;
_    int lat, lon;
_    float nxy, theta, phi;
_    vec3 n;
_
_    RGBcolor(0, 0, 0);
_    clear();
_    zclear();
_    wild = !wild;
_    if (!wild) RGBcolor(255, 255, 255);
_
_    switch (style) {
_	case 'p': case 'P': bgnpoint(); 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) RGBcolor(rand()&255, rand()&255, rand()&255);
_	    switch (style) {
_		case 'c': circf_circle(&n, &n, r); break;
_		case 's': simple_circle(&n, &n, r, nsides); break;
_		case 't': table_circle(&n, &n, r); break;
_		case 'o':
_		case 'l': object_circle(&n, &n, r); break;
_		case 'p':
_		case 'P': v3f(&n.x); break;
_		case 'f': font_circle(&n); break;
_	    }
_	}
_    }
_    switch (style) {
_	case 'p': case 'P': endpoint(); break;
_    }
_    swapbuffers();
_}
_
_main(int ac, char **av) {
_    static char opts[] = "cstolpPf";
_    int rep, i;
_    float ns;
_
_    if (ac==2 && av[1][0]=='-') {
_	printf("usage: circ2 [NSIDES [STYLE [REP [NCIRCLES]]]]\n");
_	printf("default: circ2 10 t 50 14\n");
_	printf("style is one of [%s]\n", opts);
_	printf("note: in style p,P,f, NSIDES is actually spot radius\n");
_	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);
_    }
_
_    prefposition(200, 1223, 0, 1023);
_    foreground();		/* keeps process from being forked */
_    winopen("circ");		/* create a screen window */
_    RGBmode();			/* 24 bits per pixel */
_    zbuffer(1);			/* turn on z-buffering */
_    doublebuffer();
_    gconfig();			/* call this if you're using RGBmode */
_    reshapeviewport();		/* map viewport to screen window */
_
_    switch (style) {
_	case 'c': glcompat(GLC_OLDPOLYGON, 0); break;	/* see note below */
_	case 't': table_init(nsides); break;
_	case 'o': object_init(nsides); break;
_	case 'l': line_init(nsides); break;
_	    /* kludge: ns is used as radius for point & font modes */
_	case 'p': pntsizef(ns*2); break;
_	case 'P': anti_point_init(ns*2); break;
_	case 'f': font_init(ns); break;
_    }
_    perspective(400, 1., .1, 20.);
_    translate(0., 0., -4.);
_    rot(-80., 'x');
_    for (i=0; i<rep; i++) {
_	redisplay();
_	rot(1., 'z');
_    }
_}
_
_/*
_note: "man circf" says:
_
_     By default circf draws a polygon using the old-style scan conversion
_     algorithm, rather than the point sampling algorithm used by bgnpolygon
_     and endpolygon.  Point sampling is preferred over old-style scan
_     conversion because it has higher performance, it results in higher-
_     quality images when blending, logical operations, or stenciling are used,
_     and it supports new rendering features such as texture mapping, fog, and
_     antialiasing.  Call glcompat(GLC_OLDPOLYGON,0) to enable point sampling
_     of polygons generated by circf.
_*/
EOF11868
echo extracting demo.csh
sed 's/^_//' <<'EOF11869' >demo.csh
_#!/bin/csh -f
_# shell script to give a quick demo
_foreach n (3 20)
_    foreach x (c s t o l p P f)
_	echo circ2 $n $x 20 6
_	circ2 $n $x 20 6
_    end
_end
EOF11869
exit