Meta-Ball Rendering anyone?
Paul Haeberli
paul at manray.asd.sgi.com
Thu Feb 7 07:03:59 AEST 1991
/*
* Hi -
* I've been trying to get this Meta-Ball renderer
* to work for me with no success. Does anyone out there
* know of a public domain Meta-Ball renderer, or does
* anyone know how to make this program stop core dumping?
*
* Meta-ball ray tracing program from japanese PIXEL magazine
* Number 2, 1989. The only comments available were in Kanji.
* If you can figure out how to make this work, please let me know.
*
* To compile on SGI machine:
* cc blob.c -o blob -lgl_s -lm
*
* paul at sgi.com
* paul haeberli
* 415-962-3665
*
* Here is the data file provided in the article:
*
{
{
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
100 300 220 420
}
*
*
*/
#include "stdio.h"
#include "math.h"
#include "string.h"
#include "ctype.h"
#include "assert.h"
#define PI (3.1415265)
#define THD (1.0/3.0)
#define PMAX (200)
#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;
printf("line %d\n",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];
hx = 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;
if(t[pt->spl+1]<np->t0)
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]);
fscanf(f,"%d%d%d",&cr[pmax],&cg[pmax],&cb[pmax]);
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;
}
More information about the Comp.sys.sgi
mailing list