paul@manray.asd.sgi.com (Paul Haeberli) (02/07/91)
/*
* 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@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;
}jk87377@korppi.tut.fi (Kouhia Juhana Krister) (02/07/91)
In article <84263@sgi.sgi.com> paul@manray.asd.sgi.com (Paul Haeberli) writes: > * anyone know how to make this program stop core dumping? > float sx, sy, op; [deleted] > for(i=y_b; i<y_e; i++) { > sy = 200-i; I'm not sure, but I would like to write sy = 200.0-(float)i; instead of what you have. > printf("line %d\n",i); > for(j=x_b; j<x_e; j++) { > sx = 320-j; The same. >trace() >{ I haven't been in C-language courses yet, so, this is a wild guess: Need you declare the type of the trace()? Have you an ANSI-C there? If the compiler assume that trace() is int-type, that's ok. But remember, C-language is for human - Fortran is what it is... I would like to write 'int trace()'. > if(tt<0) tt<0.0 I'm sorry about deleting all other code and not telling where all 'errors' are. Thank you for posting the code! Juhana Kouhia jk87377@tut.fi
drw900@anusf.anu.edu.au ("Drew R Whitehouse") (02/08/91)
In article <84263@sgi.sgi.com>, paul@manray.asd.sgi.com (Paul Haeberli) writes: |> /* |> * 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? Here's something to try out, it produces a picture - I'm not sure if this was what they looked like in the article !! /*---------------------------------------------------------------------------*/ /* Drew Whitehouse, E-mail: drw900@anusf.anu.edu.au */ /* Visualization Group, Fax : +61 (0)6 247 3425 */ /* Australian National University Phone : +61 (0)6 249 5985 */ /* Supercomputer Facility. */ /* GPO Box 4, Canberra ACT Australia 2601. */ /*---------------------------------------------------------------------------*/ /* * 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@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 = 700.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(1000,1000); /* foreground(); added for debugging, Drew */ 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; ...........buggy !!!!!! */ if(np == NULL || pt->spl>2) /* this maybe ??????? Drew */ 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\n",&cx[pmax],&cy[pmax],&cz[pmax]); fscanf(f,"%f %f %f\n",&tx[pmax],&ty[pmax],&tz[pmax]); /* had a simple scanf bug in the next line, Drew */ fscanf(f,"%f %f %f\n",&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; }
phs424l@vaxc.cc.monash.edu.au (02/13/91)
drw900@anusf.anu.edu.au ("Drew R Whitehouse") writes: > |> /* > |> * 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? > > > > Here's something to try out, it produces a picture - I'm not > sure if this was what they looked like in the article !! Is this for a Silicon Graphics machine? If not, what? Daf. -- ------------------------------------------------------------------------------- Davyd Norris : "Strange people physicists... In my opinion Physics Dept. : those that aren't already dead are in some Monash University, : way very, very sick." - Douglas Adams, Vic, Australia. : "The Long Dark Tea-time Of The Soul" -------------------------------------------------------------------------------
naughton@wind.Eng.Sun.COM (Patrick Naughton) (02/16/91)
phs424l@vaxc.cc.monash.edu.au writes: |> |> drw900@anusf.anu.edu.au ("Drew R Whitehouse") writes: |> > |> /* |> > |> * 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? |> > |> > |> > |> > Here's something to try out, it produces a picture - I'm not |> > sure if this was what they looked like in the article !! |> |> Is this for a Silicon Graphics machine? If not, what? |> Daf. |> -- |> ------------------------------------------------------------------------------- |> Davyd Norris : "Strange people physicists... In my opinion |> Physics Dept. : those that aren't already dead are in some |> Monash University, : way very, very sick." - Douglas Adams, |> Vic, Australia. : "The Long Dark Tea-time Of The Soul" |> ------------------------------------------------------------------------------- Yes it is for the SGI as is... but link it with this file and it will run on X11. (You'll need a 24 bit screen to make it work well...:-) #include <X11/Xlib.h> Display *dpy = 0; int screen; Colormap cmap; GC gc; Window win; int winw; int winh; void prefsize(w, h) { winw = w; winh = h; } void winopen(name) char *name; { dpy = XOpenDisplay(0); screen = DefaultScreen(dpy); cmap = DefaultColormap(dpy, screen); gc = DefaultGC(dpy, screen); win = XCreateSimpleWindow(dpy, RootWindow(dpy, screen), 0, 0, winw, winh, 1, 0); XSelectInput(dpy, win, ExposureMask); XStoreName(dpy, win, name); XMapWindow(dpy, win); while (1) { XEvent ev; XNextEvent(dpy, &ev); if (ev.type == Expose) break; } } void RGBmode() { } void gconfig() { } void RGBcolor(r, g, b) { XColor color; color.pixel = 0; color.red = r << 8; color.green = g << 8; color.blue = b << 8; color.flags = DoRed | DoGreen | DoBlue; XAllocColor(dpy, cmap, &color); XSetForeground(dpy, gc, color.pixel); } void clear() { XFillRectangle(dpy, win, gc, 0, 0, winw, winh); } pnt2i(v, w) int v; int w; { XDrawPoint(dpy, win, gc, v, w); } -- ______________________________________________________________________ Patrick J. Naughton email: naughton@sun.com Sun Labs vmail: (415) 336 - 1080