/* * blob.c * * Copyright (C) 1990, 1991, Mark Podlipec, Craig E. Kolb * All rights reserved. * * This software may be freely copied, modified, and redistributed * provided that this copyright notice is preserved on all copies. * * You may not distribute this software, in whole or in part, as part of * any commercial product without the express consent of the authors. * * There is no warranty or other guarantee of fitness of this software * for any purpose. It is provided solely "as is". * * $Id: blob.c,v 2.1 1992/06/02 14:08:51 mcgregor Exp mcgregor $ * * $Log: blob.c,v $ * Revision 2.1 1992/06/02 14:08:51 mcgregor * Installed the individual metaball surface support. A step up from the * individual color specification support. !NB! Attenuation adjustment may * look kludgy, but a clear model of how to "blend" attenuation values is * not available right now. * * Revision 1.1.1.1 1992/05/08 17:50:59 mcgregor * This baby is working!!!!! I placed the metaball_surface_contrib * call in BlobIntersect per Mark's suggestion. Next-- SURFACES! * * Revision 1.2 1992/05/07 20:23:45 mcgregor * This version works, again, at generating a multiple colored blob primitive * with full texturing. The bug pointed out by Mark Podlipec still exists: * a shadowed (or is it non-shadowed?) surface uses the ambient value from * the last trip through this code. It has to do with BlobNormal _NOT_ * being called for some things. I just want to save a working edition even * if it is only 90% working. * * Revision 1.1 1992/04/27 20:17:23 mcgregor * Initial revision * * Revision 4.0 91/07/17 14:36:07 kolb * Initial version. * */ #include "libcommon/common.h" #include "libsurf/surface.h" #include "geom.h" #include "blob.h" static Methods *iBlobMethods = NULL; static char blobName[] = "blob"; unsigned long BlobTests, BlobHits; Surface init_surf; /* Preserve the surface details as found in input file */ /* * Blob/Metaball Description * * In this implementation a Blob is made up of a threshold T, and a * group of 1 or more Metaballs. Each metaball 'i' is defined by * its position (xi,yi,zi), its radius ri, and its strength si. * Around each Metaball is a density function di(Ri) defined by: * * di(Ri) = c4i * Ri^4 + c2i * Ri^2 + c0i 0 <= Ri <= ri * di(Ri) = 0 ri < Ri * * where * c4i = si / (ri^4) * c2i = -(2 * si) / (ri^2) * c0i = si * Ri^2 is the distance from a point (x,y,z) to the center of * the Metaball - (xi,yi,zi). * * The density function looks sort of like a W (Float-u) with the * center hump crossing the d-axis at si and the two humps on the * side being tangent to the R-axis at -ri and +ri. Only the * range [0,ri] is being used. * I chose this function so that derivative = 0 at the points 0 and +ri. * This keeps the surface smooth when the densities are added. I also * wanted no Ri^3 and Ri terms, since the equations are easier to * solve without them. This led to the symmetry about the d-axis and * the derivative being equal to zero at -ri as well. * * The surface of the Blob is defined by: * * * N * --- * F(x,y,z) = \ di(Ri) = T * / * --- * i=0 * * where * * di(Ri) is given above * Ri^2 = (x-xi)^2 + (y-yi)^2 + (z-zi)^2 * N is number of Metaballs in Blob * T is the threshold * (xi,yi,zi) is the center of Metaball i * */ /***************************************************************************** * Create & return reference to a metaball. * gnm 4-30-92 added reference to blob's surface */ Blob * BlobCreate(surf, T, mlist, npoints) Surface *surf; /* surface reference for the blob */ Float T; /* Treshold */ MetaList *mlist; /* Pointer to linked list of metaballs */ int npoints; /* Number of metaballs in blob instance */ { Blob *blob; int i; MetaList *cur; /* * There has to be at least one metaball in the blob. * Note: if there's only one metaball, the blob * will be a sphere. */ if (npoints < 1) { RLerror(RL_WARN, "blob field not correct.\n"); return (Blob *)NULL; }; blob = (Blob *)Malloc(sizeof(Blob)); if((blob->surf = surf) != NULL) { /* save the original surface details ...er... oops this'll only save the most *recent* blob instance's surface. I have a plan to save each blob's (and metaballs' ??!) surface, but I'll implement that later. */ save_orig_master_surf(surf); }; blob->T = T; blob->list=(MetaVector *)Malloc((unsigned)(npoints*sizeof(MetaVector))); blob->num = npoints; /* blob will, initially be single color, not multicolored. */ blob->MColor_Flag = FALSE; /* * Initialize Metaball list */ for(i=0;imvec.c0 > -EPSILON) && (cur->mvec.c0 < EPSILON) ) || (cur->mvec.rs < EPSILON) ) { RLerror(RL_WARN,"Degenerate metaball in blob (sr).\n"); return (Blob *)NULL; }; /* store radius squared */ blob->list[i].rs = cur->mvec.rs * cur->mvec.rs; /* Calculate and store coefficients for each metaball */ blob->list[i].c0 = cur->mvec.c0; blob->list[i].c2 = -(2.0 * cur->mvec.c0) / blob->list[i].rs; blob->list[i].c4 = cur->mvec.c0 / (blob->list[i].rs * blob->list[i].rs); blob->list[i].x = cur->mvec.x; blob->list[i].y = cur->mvec.y; blob->list[i].z = cur->mvec.z; blob->list[i].surf = cur->mvec.surf; /* store metaball black body color info blob->list[i].mcolor.red = cur->mvec.mcolor.red; blob->list[i].mcolor.green = cur->mvec.mcolor.green; blob->list[i].mcolor.blue = cur->mvec.mcolor.blue; */ /* If any mcolor value is != 1.0 we need to assume that this blob will be multi colored. if((blob->list[i].mcolor.red != 1.0) || (blob->list[i].mcolor.green != 1.0) || (blob->list[i].mcolor.blue != 1.0)) */ /* if the surf spec is non-NULL then we have a multisurfaced blob */ if(blob->list[i].surf != NULL) blob->MColor_Flag = TRUE; mlist = mlist->next; free((voidstar)cur); }; /* * Allocate room for Intersection Structures and * Allocate room for an array of pointers to point to * the Intersection Structures. */ blob->ilist = (MetaInt *)Malloc( 2 * blob->num * sizeof(MetaInt)); blob->iarr = (MetaInt **)Malloc( 2 * blob->num * sizeof(MetaInt *)); return blob; } /* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: * * Function: save_orig_master_surf(surf) * Surface *surf; * * Saves the original, master surface description(if any). * * Errors: * * Notes: This needs some _SERIOUS_ reworking! * * -------------------------------------------------------------------- * who when what * -------------------------------------------------------------------- * GN McGregor * -------------------------------------------------------------------- */ int save_orig_master_surf(surf) Surface *surf; { init_surf.name = surf->name; /* Name */ init_surf.amb.r = surf->amb.r; /* Ambient 'curve' */ init_surf.amb.g = surf->amb.g; /* Ambient 'curve' */ init_surf.amb.b = surf->amb.b; /* Ambient 'curve' */ init_surf.diff.r = surf->diff.r; /* Diffuse reflection 'curve' */ init_surf.diff.g = surf->diff.g; /* Diffuse reflection 'curve' */ init_surf.diff.b = surf->diff.b; /* Diffuse reflection 'curve' */ init_surf.spec.r = surf->spec.r; /* Specular reflection 'curve' */ init_surf.spec.g = surf->spec.g; /* Specular reflection 'curve' */ init_surf.spec.b = surf->spec.b; /* Specular reflection 'curve' */ init_surf.translu.r = surf->translu.r; /* Diffuse transmission 'curve' */ init_surf.translu.g = surf->translu.g; /* Diffuse transmission 'curve' */ init_surf.translu.b = surf->translu.b; /* Diffuse transmission 'curve' */ init_surf.body.r = surf->body.r; /* Specular transmission 'curve' */ init_surf.body.g = surf->body.g; /* Specular transmission 'curve' */ init_surf.body.b = surf->body.b; /* Specular transmission 'curve' */ init_surf.srexp = surf->srexp; /* Specular reflection exponent */ init_surf.stexp = surf->stexp; /* Specular transmission exponent */ init_surf.statten = surf->statten; /* Specular transmission attenuation */ init_surf.index = surf->index; /* Index of refraction */ init_surf.reflect = surf->reflect; /* Specular reflectivity */ init_surf.transp = surf->transp; /* Specular transmittance */ init_surf.translucency = surf->translucency; /* Diffuse transmittance */ init_surf.noshadow = surf->noshadow; /* No shadowing? */ init_surf.next = surf->next; /* Next surface in list (if any) */ return; } /* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: * * Function: Methods *BlobMethods() * * * Errors: * * Notes: * * -------------------------------------------------------------------- * who when what * -------------------------------------------------------------------- * GN McGregor * -------------------------------------------------------------------- */ Methods * BlobMethods() { if (iBlobMethods == (Methods *)NULL) { iBlobMethods = MethodsCreate(); iBlobMethods->create = (GeomCreateFunc *)BlobCreate; iBlobMethods->methods = BlobMethods; iBlobMethods->name = BlobName; iBlobMethods->intersect = BlobIntersect; iBlobMethods->normal = BlobNormal; iBlobMethods->bounds = BlobBounds; iBlobMethods->stats = BlobStats; iBlobMethods->checkbounds = TRUE; iBlobMethods->closed = TRUE; }; return(iBlobMethods); } /***************************************************************************** * Function used by qsort() when sorting the Ray/Blob intersection list */ int MetaCompare(A,B) char *A,*B; { MetaInt **AA,**BB; AA = (MetaInt **) A; BB = (MetaInt **) B; if (AA[0]->bound == BB[0]->bound) return(0); if (AA[0]->bound < BB[0]->bound) return(-1); return(1); /* AA[0]->bound is > BB[0]->bound */ } /***************************************************************************** * Ray/metaball intersection test. */ int BlobIntersect(blob, ray, mindist, maxdist) Blob *blob; Ray *ray; Float mindist, *maxdist; { double c[5], root; Float dist; MetaInt *ilist,**iarr; register int i,j,inum; extern void qsort(); int metaball_surface_contrib(); /* Modify surface details accord. to di(Ri) */ BlobTests++; ilist = blob->ilist; iarr = blob->iarr; /* * The first step in calculating the Ray/Blob intersection is to * divide the Ray into intervals such that only a fixed set of * Metaballs contribute to the density function along that interval. * This is done by finding the set of intersections between the Ray * and each Metaball's Sphere/Region of influence, which has a * radius ri and is centered at (xi,yi,zi). * Intersection information is kept track of in the MetaInt * structure and consists of: * * type indicates whether this intersection is the start(R_START) * of a Region or the end(R_END) of one. * pnt the Metaball of this intersection * bound the distance from Ray origin to this intersection * * This list is then sorted by bound and used later to find the Ray/Blob * intersection. */ inum = 0; for(i=0; i < blob->num; i++) /* Check all metaballs for intersection */ { register Float xadj, yadj, zadj; register Float b, t, rs; register Float dmin,dmax; rs = blob->list[i].rs; xadj = blob->list[i].x - ray->pos.x; yadj = blob->list[i].y - ray->pos.y; zadj = blob->list[i].z - ray->pos.z; /* * Ray/Sphere of Influence intersection */ b = xadj * ray->dir.x + yadj * ray->dir.y + zadj * ray->dir.z; t = b * b - xadj * xadj - yadj * yadj - zadj * zadj + rs; /* * Don't accept imaginary or single roots. A single root is a ray * tangent to the Metaball's Sphere/Region. The Metaball's * contribution to the overall density function at this point is * zero anyway. */ if (t > 0.0) /* we have two roots and an intersection*/ { t = sqrt(t); dmin = b - t; /* * only interested in stuff in front of ray origin */ if (dmin < mindist + EPSILON) dmin = mindist + EPSILON; dmax = b + t; if (dmax > dmin) /* we have a valid Region */ { /* * Initialize min/start and max/end Intersections Structures * for this Metaball */ ilist[inum].type = R_START; ilist[inum].pnt = i; ilist[inum].bound = dmin; for(j=0;j<5;j++) ilist[inum].c[j] = 0.0; iarr[inum] = &(ilist[inum]); inum++; ilist[inum].type = R_END; ilist[inum].pnt = i; ilist[inum].bound = dmax; for(j=0;j<5;j++) ilist[inum].c[j] = 0.0; iarr[inum] = &(ilist[inum]); inum++; }; /* End of valid Region */ }; /* End of two roots */ }; /* End of loop through metaballs */ /* * If there are no Ray/Metaball intersections there will * not be a Ray/Blob intersection. Exit now. */ if (inum == 0) { return FALSE; } /* * Sort Intersection list. No sense using qsort if there's only * two intersections. * * Note: we actually aren't sorting the Intersection structures, but * an array of pointers to the Intersection structures. * This is faster than sorting the Intersection structures themselves. */ if (inum==2) { MetaInt *t; if (ilist[0].bound > ilist[1].bound) { t = iarr[0]; iarr[0] = iarr[1]; iarr[1] = t; }; } else qsort((voidstar)(iarr), (unsigned)inum, sizeof(MetaInt *), MetaCompare); /* * Finding the Ray/Blob Intersection * * The non-zero part of the density function for each Metaball is * * di(Ri) = c4i * Ri^4 + c2i * Ri^2 + c0i * * To find find the Ray/Blob intersection for one metaball * substitute for distance * * Ri^2 = (x-xi)^2 + (y-yi)^2 + (z-zi)^2 * * and then substitute the Ray equations: * * x = x0 + x1 * t * y = y0 + y1 * t * z = z0 + z1 * t * * to get a big mess :^). Actually, it's a Quartic in t and it's fully * listed farther down. Here's a short version: * * c[4] * t^4 + c[3] * t^3 + c[2] * t^2 + c[1] * t + c[0] = T * * Don't forget that the Blob is defined by the density being equal to * the threshold T. * To calculate the intersection of a Ray and two or more Metaballs, * the coefficients are calculated for each Metaball and then added * together. We can do this since we're working with polynomials. * The points of intersection are the roots of the resultant equation. * * The algorithm loops through the intersection list, calculating * the coefficients if an intersection is the start of a Region and * adding them to all intersections in that region. * When it detects a valid interval, it calculates the * roots from the starting intersection's coefficients and if any of * the roots are in the interval, the smallest one is returned. * */ { /* This curly brace is for variable scoping control. */ register Float *tmpc; MetaInt *strt,*tmp; register int istrt,itmp; register int num,exitflag,inside; int clear; /* * Start the algorithm outside the first interval */ inside = 0; istrt = 0; strt = iarr[istrt]; if (strt->type != R_START) RLerror(RL_WARN,"MetaInt sanity check FAILED!\n"); /* * Loop through intersection. If a root is found the code * will return at that point. */ while( istrt < inum ) /* inum ==> number of intersections aka iarr entries */ { /* * Check for multiple intersections at the same point. * This is also where the coefficients are calculated * and spread throughout that Metaball's sphere of * influence. */ do /* BEGIN: search for valid interval(s) */ { /* find out which metaball */ i = strt->pnt; /* only at starting boundaries do this */ if (strt->type == R_START) { register MetaVector *ml; register Float a1,a0; register Float xd,yd,zd; register Float c4,c2,c0; /* multi colored blob variables */ register Float Ri2; /* Ri^2 */ register Float Ri4; /* Ri^4 */ register Float DF; /* di(Ri) / T */ /* we're inside */ inside++; /* As promised, the full equations * * c[4] = c4*a2*a2; * c[3] = 4.0*c4*a1*a2; * c[2] = 4.0*c4*a1*a1 + 2.0*c4*a2*a0 + c2*a2; * c[1] = 4.0*c4*a1*a0 + 2.0*c2*a1; * c[0] = c4*a0*a0 + c2*a0 + c0; * * where * a2 = (x1*x1 + y1*y1 + z1*z1) = 1.0 because the ray * is normalized * a1 = (xd*x1 + yd*y1 + zd*z1) * a0 = (xd*xd + yd*yd + zd*zd) * xd = (x0 - xi) * yd = (y0 - yi) * zd = (z0 - zi) * (xi,yi,zi) is center of Metaball * (x0,y0,z0) is Ray origin * (x1,y1,z1) is normalized Ray direction * c4,c2,c0 are the coefficients for the * Metaball's density function * */ /* * Some compilers would recalculate * this each time its used. */ ml = &(blob->list[i]); xd = ray->pos.x - ml->x; yd = ray->pos.y - ml->y; zd = ray->pos.z - ml->z; a1 = (xd * ray->dir.x + yd * ray->dir.y + zd * ray->dir.z); a0 = (xd * xd + yd * yd + zd * zd); c0 = ml->c0; c2 = ml->c2; c4 = ml->c4; c[4] = c4; c[3] = 4.0*c4*a1; c[2] = 2.0*c4*(2.0*a1*a1 + a0) + c2; c[1] = 2.0*a1*(2.0*c4*a0 + c2); c[0] = c4*a0*a0 + c2*a0 + c0; /* * Since we've gone through the trouble of calculating the * coefficients, put them where they'll be used. * Starting at the current intersection(It's also the start of * a region), add the coefficients to each intersection until * we reach the end of this region. */ itmp = istrt; tmp = strt; while( (tmp->pnt != i) || (tmp->type != R_END) ) { tmpc = tmp->c; for(j=0;j<5;j++) *tmpc++ += c[j]; itmp++; tmp = iarr[itmp]; } } /* End of start of a Region */ /* * Since the intersection wasn't the start * of a region, it must the end of one. */ else inside--; /* * If we are inside a region(or regions) and the next * intersection is not at the same place, then we have * the start of an interval. Set the exitflag. */ if ((inside > 0) && (strt->bound != iarr[istrt+1]->bound) ) exitflag = 1; else /* move to next intersection */ { istrt++; strt = iarr[istrt]; exitflag = 0; } /* * Check to see if we're at the end. If so then exit with * no intersection found. */ if (istrt >= inum) { return FALSE; } } while(!exitflag); /* END: of Search for valid interval (do-while) */ /* * Find Roots along this interval */ /* get coefficients from Intersection structure */ tmpc = strt->c; for(j=0;j<5;j++) c[j] = *tmpc++; /* Don't forget to put in threshold */ c[0] -= blob->T; /* Use Graphic Gem's Quartic Root Finder */ num = FindRoot(c,4,strt->bound,iarr[istrt+1]->bound,&root,-1); if ( num && (root > mindist) && (root < *maxdist) ) { *maxdist = root; BlobHits++; if(blob->MColor_Flag) metaball_surface_contrib(blob,ray,root); return TRUE; } /* * Move to next intersection */ istrt++; strt = iarr[istrt]; }; /* End of loop through Intersection List */ }; /* End of Solving for Ray/Blob Intersection */ /* * return negative */ return FALSE; } /*********************************************** * Find the Normal of a Blob at a given point * * The normal of a surface at a point (x0,y0,z0) is the gradient of that * surface at that point. The gradient is the vector * * Fx(x0,y0,z0) , Fy(x0,y0,z0) , Fz(x0,y0,z0) * * where Fx(),Fy(),Fz() are the partial dervitives of F() with respect * to x,y and z respectively. Since the surface of a Blob is made up * of the sum of one or more polynomials, the normal of a Blob at a point * is the sum of the gradients of the individual polynomials at that point. * The individual polynomials in this case are di(Ri) i = 0,...,num . * * To find the gradient of the contributing polynomials * take di(Ri) and substitute U = Ri^2; * * di(U) = c4i * U^2 + c2i * U + c0i * * dix(U) = 2 * c4i * Ux * U + c2i * Ux * * U = (x-xi)^2 + (y-yi)^2 + (z-zi)^2 * * Ux = 2 * (x-xi) * * Rearranging we get * * dix(x0,y0,z0) = 4 * c4i * xdi * disti + 2 * c2i * xdi * diy(x0,y0,z0) = 4 * c4i * ydi * disti + 2 * c2i * ydi * diz(x0,y0,z0) = 4 * c4i * zdi * disti + 2 * c2i * zdi * * where * xdi = x0-xi * ydi = y0-yi * zdi = y0-yi * disti = xdi*xdi + ydi*ydi + zdi*zdi * (xi,yi,zi) is the center of Metaball i * c4i,c2i,c0i are the coefficients of Metaball i's density function */ int BlobNormal(blob, pos, nrm, gnrm) Blob *blob; Vector *pos, *nrm, *gnrm; { register int i; int clear; /* * Initialize normals to zero */ nrm->x = nrm->y = nrm->z = 0.0; /* * Loop through Metaballs. If the point is within a Metaball's * Sphere of influence, calculate the gradient and add it to the * normals */ for(i=0;i < blob->num; i++) { register MetaVector *sl; register Float dist,xd,yd,zd; sl = &(blob->list[i]); xd = pos->x - sl->x; yd = pos->y - sl->y; zd = pos->z - sl->z; dist = xd*xd + yd*yd + zd*zd; if (dist <= sl->rs ) { register Float temp; register Float Ri4, Ri2, DF; /* temp is negative so normal points out of blob */ temp = -2.0 * (2.0 * sl->c4 * dist + sl->c2); nrm->x += xd * temp; nrm->y += yd * temp; nrm->z += zd * temp; }; }; (void)VecNormalize(nrm); *gnrm = *nrm; return FALSE; } /***************************************************************************** * Calculate the extent of the Blob */ void BlobBounds(blob, bounds) Blob *blob; Float bounds[2][3]; { Float r; Float minx,miny,minz,maxx,maxy,maxz; int i; /* * Loop through Metaballs to find the minimun and maximum in each * direction. */ for( i=0; i< blob->num; i++) { register Float d; r = sqrt(blob->list[i].rs); if (i==0) { minx = blob->list[i].x - r; miny = blob->list[i].y - r; minz = blob->list[i].z - r; maxx = blob->list[i].x + r; maxy = blob->list[i].y + r; maxz = blob->list[i].z + r; } else { d = blob->list[i].x - r; if (dlist[i].y - r; if (dlist[i].z - r; if (dlist[i].x + r; if (d>maxx) maxx = d; d = blob->list[i].y + r; if (d>maxy) maxy = d; d = blob->list[i].z + r; if (d>maxz) maxz = d; }; }; bounds[LOW][X] = minx; bounds[HIGH][X] = maxx; bounds[LOW][Y] = miny; bounds[HIGH][Y] = maxy; bounds[LOW][Z] = minz; bounds[HIGH][Z] = maxz; return; } char * BlobName() { return blobName; } void BlobStats(tests, hits) unsigned long *tests, *hits; { *tests = BlobTests; *hits = BlobHits; return; } void BlobMethodRegister(meth) UserMethodType meth; { if (iBlobMethods) iBlobMethods->user = meth; return; } /* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: * * Function: restore_orig_surf_descrip(blob) * restores the original surface details as they were passed * into Rayshade for the blob instance. * Errors: None at this time. * * Notes: This _SHOULD_ bullet proof the code, but where????? * * -------------------------------------------------------------------- * who when what * -------------------------------------------------------------------- * GN McGregor 6-May-92 development revisions * -------------------------------------------------------------------- */ int restore_orig_surf_descrip(blob) Blob *blob; { Surface *surf; /* As the routine name suggests we are restoring the original surface description as it was handed in to the program for the blob. While this may be a bit of overkill, it does make this part of the code a bit more bullet proof. */ if((surf = blob->surf)) { surf->name = init_surf.name; /* Name */ surf->amb.r = init_surf.amb.r; /* Ambient 'curve' */ surf->amb.g = init_surf.amb.g; /* Ambient 'curve' */ surf->amb.b = init_surf.amb.b; /* Ambient 'curve' */ surf->diff.r = init_surf.diff.r; /* Diffuse reflection 'curve' */ surf->diff.g = init_surf.diff.g; /* Diffuse reflection 'curve' */ surf->diff.b = init_surf.diff.b; /* Diffuse reflection 'curve' */ surf->spec.r = init_surf.spec.r; /* Specular reflection 'curve' */ surf->spec.g = init_surf.spec.g; /* Specular reflection 'curve' */ surf->spec.b = init_surf.spec.b; /* Specular reflection 'curve' */ surf->translu.r = init_surf.translu.r; /* Diffuse transmission 'curve' */ surf->translu.g = init_surf.translu.g; /* Diffuse transmission 'curve' */ surf->translu.b = init_surf.translu.b; /* Diffuse transmission 'curve' */ surf->body.r = init_surf.body.r; /* Specular transmission 'curve' */ surf->body.g = init_surf.body.g; /* Specular transmission 'curve' */ surf->body.b = init_surf.body.b; /* Specular transmission 'curve' */ surf->srexp = init_surf.srexp; /* Specular reflection exponent */ surf->stexp = init_surf.stexp; /* Specular transmission exponent */ surf->statten = init_surf.statten; /* Specular transmission attenuation */ surf->index = init_surf.index; /* Index of refraction */ surf->reflect = init_surf.reflect; /* Specular reflectivity */ surf->transp = init_surf.transp; /* Specular transmittance */ surf->translucency = init_surf.translucency; /* Diffuse transmittance */ surf->noshadow = init_surf.noshadow; /* No shadowing? */ surf->next = init_surf.next; /* Next surface in list (if any) */ }; return(GOOD); } /* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: * * Function: metaball_surface_contrib(blob, ray, dist) * * * Errors: * * Notes: * * -------------------------------------------------------------------- * who when what * -------------------------------------------------------------------- * GN McGregor * -------------------------------------------------------------------- */ int metaball_surface_contrib(blob, ray, dist) Blob *blob; Ray *ray; Float dist; { Vector intersect_pos; int i; /* We have the distance from the ray's origin to the point of intersection. We have to find the x, y, z values for the point equal to dist along the ray.The intersection x, y, z values can be drived by extending the unit vector like so... */ VecAddScaled(ray->pos, dist, ray->dir, &intersect_pos); /* Prepare to modify the surface */ tweak_blob_surf(1, blob->surf, 1.0, &(blob->list[0])); for(i=0;i < blob->num; i++) /* Scan metaballs for surface contributions */ { register MetaVector *sl; register Float range,xd,yd,zd; register Float Ri2; /* R(i)^2 in Mark's Equation up above. */ sl = &(blob->list[i]); xd = intersect_pos.x - sl->x; yd = intersect_pos.y - sl->y; zd = intersect_pos.z - sl->z; Ri2 = xd*xd + yd*yd + zd*zd; if (Ri2 <= sl->rs ) { register Float DF; DF = (((sl->c4*Ri2) + sl->c2)*Ri2 + sl->c0) / blob->T; /*DF = ((sl->c4*Ri4) + (sl->c2*Ri2) + sl->c0) / blob->T;*/ tweak_blob_surf(0, blob->surf, DF, sl->surf); }; /* END metaball density calc. */ }; return(GOOD); } /* ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: * * Function: tweak_blob_surf(clear, surf, df, mcolor) * * This routine will fiddle with those components of the surface that * have RGB values as part of their specs. The more sophisticated * concept and operation would involve allowing each metaball to possess * its own surface. After this exercise that possibility may prove easier * and more useful. * * * Errors: * * Notes: * * -------------------------------------------------------------------- * who when what * -------------------------------------------------------------------- * GN McGregor 6-May-92 development revisions * -------------------------------------------------------------------- */ int tweak_blob_surf(clear, surf, df, msurf) int clear; /* flag-signal to init surface spec prior to tweaking */ Surface *surf; /* pointer to surface spec. */ Float df; /* density fraction-> contribution weight of metaball */ Surface *msurf; /* metaball's surface spec. */ { if(surf) { if(clear) /* We want to prepare the surface for its dynamic update. */ { surf->amb.r = 0.0; surf->amb.g = 0.0; surf->amb.b = 0.0; surf->diff.r = 0.0; surf->diff.g = 0.0; surf->diff.b = 0.0; surf->spec.r = 0.0; surf->spec.g = 0.0; surf->spec.b = 0.0; surf->translu.r = 0.0; surf->translu.g = 0.0; surf->translu.b = 0.0; surf->body.r = 0.0; surf->body.g = 0.0; surf->body.b = 0.0; surf->srexp = 0.0; surf->stexp = 0.0; surf->statten = 0.0; surf->index = 0.0; surf->reflect = 0.0; surf->transp = 0.0; surf->translucency = 0.0; } else /* We are updating the surface components. */ { /* The new surface component color is the weighted sum of each of the contributing metaballs colors. The weighting factor is fractiobnal contribution of the metaball to the blob's density function. ( i.e. d(i) / T ). */ /* Trying out: use surface specs as a scale value. */ surf->amb.r += (msurf->amb.r * df) * init_surf.amb.r; surf->amb.g += (msurf->amb.g * df) * init_surf.amb.g; surf->amb.b += (msurf->amb.b * df) * init_surf.amb.b; surf->diff.r += (msurf->diff.r * df) * init_surf.diff.r; surf->diff.g += (msurf->diff.g * df) * init_surf.diff.g; surf->diff.b += (msurf->diff.b * df) * init_surf.diff.b; surf->spec.r += (msurf->spec.r * df) * init_surf.spec.r; surf->spec.g += (msurf->spec.g * df) * init_surf.spec.g; surf->spec.b += (msurf->spec.b * df) * init_surf.spec.b; surf->translu.r += (msurf->translu.r * df) * init_surf.translu.r; surf->translu.g += (msurf->translu.g * df) * init_surf.translu.g; surf->translu.b += (msurf->translu.b * df) * init_surf.translu.b; surf->body.r += (msurf->body.r * df) * init_surf.body.r; surf->body.g += (msurf->body.g * df) * init_surf.body.g; surf->body.b += (msurf->body.b * df) * init_surf.body.b; /* Model physical behavior as sum of weighted contributions from metaballs. */ surf->srexp += (msurf->srexp * df);/* + init_surf.srexp;*/ surf->stexp += (msurf->stexp * df);/* + init_surf.stexp;*/ /* Attenuation _is_ a problem... Let's try a method where the attenuation is equal to the maximum value we see. */ if(msurf->statten > init_surf.statten) surf->statten = (msurf->statten > surf->statten) ? msurf->statten : surf->statten; else surf->statten = (init_surf.statten > surf->statten) ? init_surf.statten : surf->statten; surf->index += (msurf->index * df);/* * init_surf.index;*/ surf->reflect += (msurf->reflect * df);/* * init_surf.reflect;*/ surf->transp += (msurf->transp * df);/* * init_surf.transp;*/ surf->translucency += (msurf->translucency * df);/* * init_surf.translucency;*/ }; }; return(GOOD); }