Xref: utzoo comp.sys.sgi:8221 comp.graphics:15910
Path: utzoo!mnetor!tmsoft!torsqnt!news-server.csri.toronto.edu!bonnie.concordia.ca!uunet!olivea!decwrl!sgi!paul@manray.asd.sgi.com
From: paul@manray.asd.sgi.com (Paul Haeberli)
Newsgroups: comp.sys.sgi,comp.graphics
Subject: Working Meta-Ball renderer
Message-ID: <84792@sgi.sgi.com>
Date: 10 Feb 91 22:13:37 GMT
Sender: guest@sgi.sgi.com
Organization: Silicon Graphics, Inc., Mountain View, CA
Lines: 594

Here is a working Meta-Ball renderer from PIXEL magazine.  Thanks
to everyone how helped find the bugs and typos in the program to finally
make it work. . . . . .

/*
 *	blob -
 *		Working Meta-ball ray tracing program from japanese PIXEL 
 *	magazine Number 2, 1989.  The only comments available were in Kanji. 
 *
 *	To compile on SGI machine: 
 *		cc blob.c -o blob -lgl_s -lm
 *
 *	paul@sgi.com				
 *	paul haeberli
 *	415-962-3665
 *
 *	Here is the data file provided in the article:
 *
 *	The format of each blob is:
 *
 *	p
 *	78.0 63.0 63.0 		- x y z size
 *	0.0 70.0 0.0 		- x y z position 
 *	0.0 0.0 -90.0		- x y z rotation 
 *	100 200 250		- rgb color
 *	1.5			- weight
 *
 *	The format of the camera is:
 *
 *	c
 *	0.0 0.0 -150.0		- x y z camera position
 *	0.0 0.0 0.0		- x y z camera rotation
 *	0.0 0.0 -1.0		- x y z light direction vector
 *	0.5 			- a ambient
 *	210 120 180		- r g b back color
 *	100 300 220 420		- y y x x range of pixels to render
 *
{
{
p
78.0 63.0 63.0	
0.0 70.0 0.0	
0.0 0.0 -90.0
100 200 250
1.5

p
31.5 12.369 12.369
0.0 39.5 0.0
0.0 0.0 -90.0
200 130 100
1.5

p
67.98 77.28 77.28
-1.0 -14.5 0.0
0.0 0.0 -90.0
200 220 60
7.34

p
117.61 28.30 28.30
-74.0 33.0 0.0
0.0 0.0 37.74
20 200 100
1.5

p
123.54 26.83 26.83
77.0 33.0 0.0
0.0 0.0 -29.05
20 100 200
1.5

p
16.02 14.5 24.65 
-1.0 -19.0 -49.0
0.0 0.0 -95.19
20 20 200
6.9
}

c
0.0 0.0 -150.0
0.0 0.0 0.0
0.0 0.0 -1.0
0.5 
210 120 180
50 350 170 470

}
 *
 *
 */
#include "stdio.h"
#include "math.h"
#include "string.h"
#include "ctype.h"
#include "assert.h"

#define FIXBUGS			/* this fixes bugs in the printed article */

#define PI	(3.1415265)
#define THD	(1.0/3.0)
#define PMAX	(200)		/* The ultimate number 200 of meta-balls */
#define INFIN	(1000000.0)	

struct tlist {
    int mel, spl;
    float t0, t1, t2;
    struct tlist *next;
};

struct vect {
    float x, y, z;
};

struct vect eye, ray;

float nx, ny, nz, px, py, pz;
float sqrx, sqry, sqrz;
float dg = PI/180.0;
float rt_x, rt_y, rt_z;

FILE *fp;
float a[PMAX], b[PMAX], c[PMAX];
float cx[PMAX], cy[PMAX], cz[PMAX];
float tx[PMAX], ty[PMAX], tz[PMAX];
float cr[PMAX], cg[PMAX], cb[PMAX];
float wgt[PMAX];
float cc[PMAX];
float la[PMAX], lb[PMAX], lc[PMAX];
float lx, ly, lz, amb;
float daen(), sgn(), sgr(), ei(), wi();
float pix_r, pix_g, pix_b;
int jx, iy, pmax;
int b_r, b_g, b_b;
int y_b, y_e, x_b, x_e;
float r00, r01, r02, r10, r11, r12, r20, r21, r22;

main(argc,argv)
int argc;
char **argv;
{
    float sz = 100.0;
    float sx, sy, op;
    float lcx, lcy, lcz;
    float tmp_x, tmp_y, tmp_z;
    int i, j, k;

    if(argc<2) {
	fprintf(stderr,"usage: blob blob.dat\n");
	exit(1);
    }
    fp = fopen(argv[1],"r");
    if(readfile(fp,0) != 0) {
		inerror();
    }
    for(k=0; k<pmax; k++) {
	a[k] = sgn(a[k])/a[k]/a[k];
	b[k] = sgn(b[k])/b[k]/b[k];
	c[k] = sgn(c[k])/c[k]/c[k];
	lcx = eye.x-cx[k];
	lcy = eye.y-cy[k];
	lcz = eye.z-cz[k];

	myrot(tx[k],ty[k],tz[k]);
	tmp_x = (a[k]*r00*r00+b[k]*r10*r10+c[k]*r20*r20);
	tmp_y = (a[k]*r01*r01+b[k]*r11*r11+c[k]*r21*r21);
	tmp_z = (a[k]*r02*r02+b[k]*r12*r12+c[k]*r22*r22);

	tx[k] = 2.0*(a[k]*r01*r02+b[k]*r11*r12+c[k]*r21*r22);
	ty[k] = 2.0*(a[k]*r00*r02+b[k]*r10*r12+c[k]*r20*r22);
	tz[k] = 2.0*(a[k]*r00*r01+b[k]*r10*r11+c[k]*r20*r21);

	a[k] = tmp_x;
	b[k] = tmp_y;
	c[k] = tmp_z;

	cc[k] = lcx*lcx*a[k]+lcy*lcy*b[k]+lcz*lcz*c[k]-1.0;
	cc[k] = cc[k]+lcy*lcz*tx[k]+lcx*lcz*ty[k]+lcx*lcy*tz[k];
	la[k] = lcx*a[k];
	lb[k] = lcy*b[k];
	lc[k] = lcz*c[k];
    }

    myrot(rt_x,rt_y,rt_z);

    prefsize(640,400);
    winopen("blob");
    RGBmode();
    gconfig();
    RGBcolor(0,0,255);
    clear();

    for(i=y_b; i<y_e; i++) {
	sy = 200-i;
	for(j=x_b; j<x_e; j++) {
	    sx = 320-j;
	    jx = j;
	    iy = i;
	    op = sqrt(sx*sx+sy*sy+sz*sz);

	    tmp_x = sx/op;
	    tmp_y = sy/op;
	    tmp_z = sz/op;

	    ray.x = tmp_x*r00+tmp_y*r01+tmp_z*r02;
	    ray.y = tmp_x*r10+tmp_y*r11+tmp_z*r12;
	    ray.z = tmp_x*r20+tmp_y*r21+tmp_z*r22;
	    if(trace())
		dsp();
	    else
		back();
	}
    }
    while(1) 
	sleep(10);
}

trace()
{
    static int k, kk, flg;
    static float bright, td, tm, hi, da, db, dc, dd, tt, hx, hy, hz;
    static struct tlist list[50];
    static struct tlist *pt, *np, root;
    static float ds, cmp, op, tmpa;
    static float t[] = {0.0, 0.0, 0.0, 0.0, -INFIN};
    static float cos_a, phong, rfx, rfy, rfz;

    sqrx = ray.x*ray.x;
    sqry = ray.y*ray.y;
    sqrz = ray.z*ray.z;

    for(k=kk=0; k<pmax; k++) {
	if( (list[kk].t1 = daen(k,&list[kk].t2)) != INFIN ) {
	    list[kk].t0 = list[kk].t1;
	    list[kk].mel = k;
	    list[kk].spl = 0;
	    kk++;
	}
    }

    if(kk != 0) {
	da = db = dc = 0.0;

	flg = 0;
	for(;;) {
	    if(flg==0)
		flg = 1;
	    else {
		++(pt->spl);
		pt->t0 = t[pt->spl]; 
	    }

	    root.t0 = -100000.0;
	    mklist(&root,list,kk);
	    for(k=0; k<kk; k++) 
		mklist(&list[k],list,kk);
	    pt = root.next;
	    np = pt->next;

	    if(np == NULL && pt->spl>2)
		return 0;
	    td = ((pt->t2)-(pt->t1))*0.5;
	    tm = pt->t1+td;

	    hx = eye.x+tm*ray.x-cx[pt->mel];
	    hy = eye.y+tm*ray.y-cy[pt->mel];
	    hz = eye.z+tm*ray.z-cz[pt->mel];

	    hi = sqrt(hx*hx*a[pt->mel]+hy*hy*b[pt->mel]+hz*hz*c[pt->mel]
		+tx[pt->mel]*hy*hz+ty[pt->mel]*hx*hz+tz[pt->mel]*hx*hy);
	    ds = td*(1.0-ei(hi));

	    t[0] = pt->t1;
	    t[1] = pt->t1 + ds;
	    t[2] = pt->t2 - ds;
	    t[3] = pt->t2;

	    switch(pt->spl) {
		case 0:	
		    tmpa = wi(hi,pt->mel)/(ds*td);
		    break;
		case 1:	
		    tmpa = wi(hi,pt->mel)/(ei(hi)*-ds*td);
		    break;
		case 2:	
		    tmpa = wi(hi,pt->mel)/(ei(hi)*ds*td);
		    break;
		case 3:
		    tmpa = wi(hi,pt->mel)/(-ds*td);
		    break;
	    }

	    da += tmpa;
	    db += -2.0*tmpa*t[pt->spl];
	    dc += tmpa*t[pt->spl]*t[pt->spl];

	    if(da == 0.0)
		continue;

	    dd = db*db-4.0*da*(dc-1.0);
	    if(dd>0.0)
		tt = (-db+sqrt(dd))*0.5/da;
	    else
		continue;

	    if(tt<0)
		continue;

#ifdef FIXBUGS
	    if(np == 0 || t[pt->spl+1]<np->t0) 
#else
	    if(t[pt->spl+1]<np->t0)
#endif
		cmp = t[pt->spl+1];
	    else
		cmp = np->t0;

	    if(cmp == -INFIN)
		cmp = np->t0;

	    
	    if(pt->t0<tt && tt<cmp) {
		px = eye.x+ray.x*tt;
		py = eye.y+ray.y*tt;
		pz = eye.z+ray.z*tt;
		getnorm(root.next);

		ref(&rfx,&rfy,&rfz);
		cos_a = rfx*lx+rfy*ly+rfz*lz;
		if(cos_a<0.0)
		    cos_a = 0.0;
		phong = cos_a*cos_a*cos_a*cos_a*cos_a*cos_a;
		bright = nx*lx+ny*ly+nz*lz+phong+amb;
		if(bright<0.0)
		    bright = 0.0;

		pix_r = pix_r*bright;
		pix_g = pix_g*bright;
		pix_b = pix_b*bright;

		return 1;
	    }
	}
    } else
	return 0;
}

float daen(k,l)
int k;
float *l;
{
    float aa, bb, d;
    float lcx, lcy, lcz;

    lcx = eye.x-cx[k];
    lcy = eye.y-cy[k];
    lcz = eye.z-cz[k];
    aa = sqrx*a[k]+sqry*b[k]+sqrz*c[k];
    if(aa == 0.0)
	return INFIN;
    bb = 2.0*(ray.x*la[k]+ray.y*lb[k]+ray.z*lc[k]);

    aa = aa+ray.y*ray.z*tx[k]+ray.x*ray.z*ty[k]+ray.x*ray.y*tz[k];
    bb = bb+(lcz*ray.y+lcy*ray.z)*tx[k]+(lcx*ray.z+lcz*ray.x)*ty[k]
	   +(lcx*ray.y+lcy*ray.x)*tz[k];

    d = bb*bb-4.0*aa*cc[k];
    if(d<0.0)
	return INFIN;
    *l = (-bb+sqrt(d))*0.5/aa;
    return (-bb-sqrt(d))*0.5/aa;
}

dsp()
{
    int rr, gg ,bb;

    rr = pix_r;
    gg = pix_g;
    bb = pix_b;

    if(rr>255)
	rr = 255;
    if(gg>255)
	gg = 255;
    if(bb>255)
	bb = 255;
    cgpset(jx,iy,rr,gg,bb);
}

back()
{
    cgpset(jx,iy,b_r,b_g,b_b);
}

float sgn(v)
float v;
{
    if(v != 0.0)
	return ((v<0.0) ? -1.0 : 1.0);
    else
	return 0;
}

skip(f)
FILE *f;
{
    int ch;

    while(isspace(ch=getc(f)) && ch != EOF)
        ;
    return ch;
}

readfile(f,chr) 
FILE *f;
int chr;
{
    int ch, i, j, n;
    FILE *f2;

    pmax = 0;
    while(1) {
	switch(ch = skip(f)) {
	    case ';':
		while((ch=getc(f)) != '\n' && ch != EOF)
		    ;
		break;
	    case '{':
		if(readfile(f,chr) != 1)
		    inerror();
		break;
	    case '}':
		return 1;
	    case 'p':
		fscanf(f,"%f%f%f",&a[pmax],&b[pmax],&c[pmax]);
		fscanf(f,"%f%f%f",&cx[pmax],&cy[pmax],&cz[pmax]);
		fscanf(f,"%f%f%f",&tx[pmax],&ty[pmax],&tz[pmax]);
#ifdef FIXBUGS
		fscanf(f,"%f%f%f",&cr[pmax],&cg[pmax],&cb[pmax]);
#else
		fscanf(f,"%d%d%d",&cr[pmax],&cg[pmax],&cb[pmax]);
#endif
		fscanf(f,"%f",&wgt[pmax]);
		pmax++;
		break;
	    case 'c':
		fscanf(f,"%f%f%f",&eye.x,&eye.y,&eye.z);
		fscanf(f,"%f%f%f",&rt_x,&rt_y,&rt_z);
		fscanf(f,"%f%f%f",&lx,&ly,&lz);
		fscanf(f,"%f",&amb);
		fscanf(f,"%d%d%d",&b_r,&b_g,&b_b);
		fscanf(f,"%d%d%d%d",&y_b,&y_e,&x_b,&x_e);
		break;
	    case EOF:
		return 0;
	    default:
		inerror();
		break;
	}
    }
}

inerror()
{
    fprintf(stderr,"input read error\n");
    exit(1);
}

mklist(p,q,r)
struct tlist *p, *q;
int r;
{
    int k;
    float dt;
    float lst = INFIN;

    for(k=0; k<r; k++) {
	if((q+k) != p) {
	    dt = ((q+k)->t0)-(p->t0);
	    if((dt<lst) && (dt>=0.0) && !(dt==0.0 && q+k<p)) {
		lst = dt;
		p->next = q+k;
	    }
	}
    }	
    if(lst==INFIN)
	p->next = NULL;
}

float wi(di,i)
float di;
int i;
{
    if((0.0<=di) && (di<=THD))
	return (wgt[i]*(1.0-3.0*di*di));
    else if((di>=THD) && (di<=1.0))
	return (1.5*wgt[i]*(1.0-di)*(1.0-di));
    else
	return 0.0;
}

float ei(r)
float r;
{

    if(r<0.22)
	return (-0.5*r+0.33);
    return (0.346*r+0.243846);
}

getnorm(p)
struct tlist *p;
{
    float tpx, tpy, tpz;
    float r, pw, nxs, nys, nzs, op;
    int k;

    pix_r = pix_g = pix_b = 0.0;
    nxs = nys = nzs = 0.0;

    do {
	k = p->mel;
	tpx = px-cx[k];
	tpy = py-cy[k];
	tpz = pz-cz[k];
	r = tpx*tpx*a[k]+tpy*tpy*b[k]+tpz*tpz*c[k]
		+tx[k]*tpy*tpz+ty[k]*tpx*tpz+tz[k]*tpx*tpy;
	if(r<0)
	    r = 0;
	pw = wi(sqrt(r),k);
	nxs += (a[k]*tpx+0.5*(tz[k]*tpy+ty[k]*tpz))*pw;
	nys += (b[k]*tpy+0.5*(tz[k]*tpx+tx[k]*tpz))*pw;
	nzs += (c[k]*tpz+0.5*(ty[k]*tpx+tx[k]*tpy))*pw;

	pix_r += (float)cr[k]*pw;
	pix_g += (float)cg[k]*pw;
	pix_b += (float)cb[k]*pw;
    } while( (p=p->next) != NULL);

    nx = nxs;
    ny = nys;
    nz = nzs;
    op = sqrt(nx*nx+ny*ny+nz*nz);
    if(op == 0.0)
	op = 1.0;
    nx = nx/op;
    ny = ny/op;
    nz = nz/op;
}

cgpset(v,w,r,g,b)
int v,w,r,g,b;
{
    RGBcolor(r,g,b);
    pnt2i(v,w);
}

ref(rx,ry,rz)
float *rx, *ry, *rz;
{
    float w;

    w = -2.0*(ray.x*nx+ray.y*ny+ray.z*nz);
    *rx = ray.x+nx*w;
    *ry = ray.y+ny*w;
    *rz = ray.z+nz*w;
}

myrot(x,y,z)
float x,y,z;
{
    float snx, cnx, sny, cny, snz, cnz;

    cnx = cos(x*dg);
    snx = sin(x*dg);
    cny = cos(y*dg);
    sny = sin(y*dg);
    cnz = cos(z*dg);
    snz = sin(z*dg);

    r00 = cny*cnz;
    r01 = snx*sny*cnz-cnx*snz;
    r02 = cnx*sny*cnz+snx*snz;
    r10 = cny*snz;
    r11 = snx*sny*snz+cnx*cnz;
    r12 = cnx*sny*snz-snx*cnz;
    r20 = -sny;
    r21 = snx*cny;
    r22 = cnx*cny;
}