| 1 | greg | 2.1 | #ifndef lint | 
| 2 | greg | 2.7 | static const char RCSid[] = "$Id: m_wgmdf.c,v 2.6 2024/12/18 17:57:06 greg Exp $"; | 
| 3 | greg | 2.1 | #endif | 
| 4 |  |  | /* | 
| 5 |  |  | *  Shading function for programmable Ward-Geisler-Moroder-Duer material. | 
| 6 |  |  | */ | 
| 7 |  |  |  | 
| 8 |  |  | #include "copyright.h" | 
| 9 |  |  |  | 
| 10 |  |  | #include  "ray.h" | 
| 11 |  |  | #include  "ambient.h" | 
| 12 |  |  | #include  "otypes.h" | 
| 13 |  |  | #include  "rtotypes.h" | 
| 14 |  |  | #include  "source.h" | 
| 15 |  |  | #include  "func.h" | 
| 16 |  |  | #include  "random.h" | 
| 17 |  |  | #include  "pmapmat.h" | 
| 18 |  |  |  | 
| 19 |  |  | #ifndef  MAXITER | 
| 20 |  |  | #define  MAXITER        10              /* maximum # specular ray attempts */ | 
| 21 |  |  | #endif | 
| 22 |  |  | /* estimate of Fresnel function */ | 
| 23 |  |  | #define  FRESNE(ci)     (exp(-5.85*(ci)) - 0.00202943064) | 
| 24 |  |  | #define  FRESTHRESH     0.017999        /* minimum specularity for approx. */ | 
| 25 |  |  |  | 
| 26 |  |  | /* | 
| 27 |  |  | *      This routine implements the anisotropic Gaussian | 
| 28 |  |  | *  model described by Ward in a 1992 Siggraph article and updated by | 
| 29 |  |  | *  Geisler-Moroder and Duer in a 2010 article in High Performance Graphics. | 
| 30 |  |  | *      We do not reorient incoming ray, using side in part to determine | 
| 31 |  |  | *  reflectance values.  Most parameters are programmable with their own | 
| 32 |  |  | *  modifiers and/or value expressions. | 
| 33 |  |  | * | 
| 34 |  |  | *  Arguments for MAT_WGMDF are: | 
| 35 |  |  | *      13+     rs_mod  rs  rs_urough rs_vrough | 
| 36 |  |  | *              ts_mod  ts  ts_urough ts_vrough | 
| 37 |  |  | *              td_mod | 
| 38 |  |  | *              ux uy uz  funcfile  transform | 
| 39 |  |  | *      0 | 
| 40 |  |  | *      9+      rfdif gfdif bfdif | 
| 41 |  |  | *              rbdif gbdif bbdif | 
| 42 |  |  | *              rtdif gtdif btdif | 
| 43 |  |  | *              A10 .. | 
| 44 |  |  | * | 
| 45 |  |  | *  Where the rs_urough or rs_vrough expression yields zero, mirror-Fresnel | 
| 46 |  |  | *  effects are computed, similar to MAT_PLASTIC and MAT_METAL.  The | 
| 47 |  |  | *  rs* expressions should not vary with incident angle, or the material | 
| 48 |  |  | *  will not be physically valid.  Similarly, the ts* expressions should | 
| 49 |  |  | *  give the same value for coincident direction vectors from either side. | 
| 50 |  |  | *      There are independent modifiers for specular reflection, | 
| 51 |  |  | *  transmission, and diffuse transmission.  Diffuse reflection | 
| 52 |  |  | *  applies the material's main modifier, which doesn't apply to | 
| 53 |  |  | *  anything else by default.  However, any of the modifiers may be | 
| 54 |  |  | *  ALIASMOD, which will use the main material modifier, or VOIDID, | 
| 55 |  |  | *  which will just be white. | 
| 56 |  |  | *      Diffuse reflection and transmission colors and patterns add to | 
| 57 |  |  | *  the specular components, and are only adjusted with mirror-Fresnel | 
| 58 |  |  | *  reflection if specular reflection is greater than FRESHTHRESH.  The | 
| 59 |  |  | *  specular transmission is likewise adjusted in such cases.  Specified | 
| 60 |  |  | *  values for all components should sum to less than 1, but like other | 
| 61 |  |  | *  Radiance materials, this is not enforced, nor is a warning issued. | 
| 62 |  |  | */ | 
| 63 |  |  | /* specularity flags */ | 
| 64 |  |  | #define  SP_REFL        01              /* has reflected specular component */ | 
| 65 |  |  | #define  SP_TRAN        02              /* has transmitted specular */ | 
| 66 |  |  | #define  SP_RPURE       04              /* mirror reflection */ | 
| 67 |  |  | #define  SP_TPURE       010             /* has view component */ | 
| 68 |  |  | #define  SP_FLAT        020             /* flat reflecting surface */ | 
| 69 |  |  | #define  SP_RBLT        040             /* reflection below sample threshold */ | 
| 70 |  |  | #define  SP_TBLT        0100            /* transmission below threshold */ | 
| 71 |  |  |  | 
| 72 |  |  | typedef struct { | 
| 73 |  |  | char            *nam;           /* modifier name */ | 
| 74 |  |  | int             hastexture;     /* has a texture? */ | 
| 75 |  |  | FVECT           pnorm;          /* perturbed normal direction */ | 
| 76 |  |  | double          pdot;           /* perturbed dot product */ | 
| 77 |  |  | SCOLOR          pcol;           /* pattern color */ | 
| 78 |  |  | } MODVAL;               /* modifier-derived values */ | 
| 79 |  |  |  | 
| 80 |  |  | typedef struct { | 
| 81 |  |  | MODVAL          mo;             /* modifier parameters */ | 
| 82 |  |  | SCOLOR          scol;           /* modified diffuse color */ | 
| 83 |  |  | } DCOMP;                /* diffuse component parameters */ | 
| 84 |  |  |  | 
| 85 |  |  | typedef struct { | 
| 86 |  |  | MODVAL          mo;             /* modifier parameters */ | 
| 87 |  |  | SCOLOR          scol;           /* modified specular color */ | 
| 88 |  |  | FVECT           u, v;           /* u and v in-plane vectors */ | 
| 89 |  |  | double          u_alpha;        /* u roughness */ | 
| 90 |  |  | double          v_alpha;        /* v roughness */ | 
| 91 |  |  | } SCOMP;                /* specular component parameters */ | 
| 92 |  |  |  | 
| 93 |  |  | typedef struct { | 
| 94 |  |  | RAY             *rp;            /* ray pointer */ | 
| 95 |  |  | OBJREC          *mtp;           /* material pointer */ | 
| 96 |  |  | MFUNC           *mf;            /* pointer to expression list */ | 
| 97 |  |  | int             specfl;         /* specularity flags, defined above */ | 
| 98 |  |  | FVECT           ulocal;         /* u-vector in local coordinates */ | 
| 99 |  |  | DCOMP           rd, td;         /* diffuse component params */ | 
| 100 |  |  | SCOMP           rs, ts;         /* specular component params */ | 
| 101 |  |  | FVECT           prdir;          /* vector in transmitted direction */ | 
| 102 |  |  | } WGMDDAT;              /* WGMD material data */ | 
| 103 |  |  |  | 
| 104 |  |  | #define clr_comps(wp)   ((wp)->specfl = 0, \ | 
| 105 |  |  | (wp)->rd.mo.nam = (wp)->td.mo.nam = \ | 
| 106 |  |  | (wp)->rs.mo.nam = (wp)->ts.mo.nam = "") | 
| 107 |  |  |  | 
| 108 |  |  | /* assign modifier values */ | 
| 109 |  |  | static int | 
| 110 |  |  | set_modval(MODVAL *mp, OBJECT omod, const RAY *r) | 
| 111 |  |  | { | 
| 112 |  |  | RAY     tr; | 
| 113 |  |  |  | 
| 114 |  |  | if (!mp->nam[0]) | 
| 115 |  |  | mp->nam = (omod == OVOID) ? VOIDID : objptr(omod)->oname; | 
| 116 |  |  | else if (!strcmp(mp->nam, VOIDID)) | 
| 117 |  |  | omod = OVOID; | 
| 118 |  |  | else if (omod == OVOID) | 
| 119 |  |  | return(0); | 
| 120 |  |  | tr = *r;                        /* independent modifier */ | 
| 121 |  |  | raytexture(&tr, omod); | 
| 122 |  |  | if (DOT(tr.pert,tr.pert) > FTINY*FTINY) { | 
| 123 |  |  | mp->pdot = raynormal(mp->pnorm, &tr); | 
| 124 |  |  | mp->hastexture = 1; | 
| 125 |  |  | } else { | 
| 126 |  |  | VCOPY(mp->pnorm, tr.ron); | 
| 127 |  |  | mp->pdot = tr.rod; | 
| 128 |  |  | mp->hastexture = 0; | 
| 129 |  |  | } | 
| 130 |  |  | copyscolor(mp->pcol, tr.pcol); | 
| 131 |  |  | return(1); | 
| 132 |  |  | } | 
| 133 |  |  |  | 
| 134 |  |  | /* fill modifier values, using previous setting if found */ | 
| 135 |  |  | static int | 
| 136 |  |  | fill_modval(MODVAL *mp, const WGMDDAT *wp) | 
| 137 |  |  | { | 
| 138 |  |  | if (mp == &wp->rd.mo) {         /* special case (should be first) */ | 
| 139 |  |  | set_modval(mp, wp->mtp->omod, wp->rp); | 
| 140 |  |  | return(1); | 
| 141 |  |  | }                               /* use main modifier? */ | 
| 142 |  |  | if (!strcmp(mp->nam, ALIASMOD) || !strcmp(mp->nam, wp->rd.mo.nam)) { | 
| 143 |  |  | *mp = wp->rd.mo; | 
| 144 |  |  | return(1); | 
| 145 |  |  | }                               /* check others */ | 
| 146 |  |  | if (mp != &wp->td.mo && !strcmp(mp->nam, wp->td.mo.nam)) { | 
| 147 |  |  | *mp = wp->td.mo; | 
| 148 |  |  | return(1); | 
| 149 |  |  | } | 
| 150 |  |  | if (mp != &wp->rs.mo && !strcmp(mp->nam, wp->rs.mo.nam)) { | 
| 151 |  |  | *mp = wp->rs.mo; | 
| 152 |  |  | return(1); | 
| 153 |  |  | } | 
| 154 |  |  | if (mp != &wp->ts.mo && !strcmp(mp->nam, wp->ts.mo.nam)) { | 
| 155 |  |  | *mp = wp->ts.mo; | 
| 156 |  |  | return(1); | 
| 157 |  |  | }                               /* new modifier */ | 
| 158 |  |  | return(set_modval(mp, lastmod(objndx(wp->mtp), mp->nam), wp->rp)); | 
| 159 |  |  | } | 
| 160 |  |  |  | 
| 161 | greg | 2.6 | /* set calculation context for given component of MAT_WGMDF */ | 
| 162 | greg | 2.5 | static int | 
| 163 |  |  | setWGMDfunc(MODVAL *mp, const WGMDDAT *wp) | 
| 164 |  |  | { | 
| 165 | greg | 2.6 | static char     lastMod[MAXSTR]; | 
| 166 | greg | 2.5 | double          sf; | 
| 167 |  |  | FVECT           vec; | 
| 168 |  |  |  | 
| 169 |  |  | if (setfunc(wp->mtp, wp->rp) == 0 && | 
| 170 |  |  | !strcmp(mp->nam, lastMod)) | 
| 171 |  |  | return(0);      /* already set */ | 
| 172 | greg | 2.6 | strcpy(lastMod, mp->nam); | 
| 173 | greg | 2.5 | /* else (re)assign special variables */ | 
| 174 | greg | 2.6 | sf = 1 - 2*(wp->rp->rod < 0); | 
| 175 |  |  | varset("RdotP`", '=', mp->pdot*sf); | 
| 176 |  |  | multv3(vec, mp->pnorm, funcxf.xfm); | 
| 177 | greg | 2.5 | sf /= funcxf.sca; | 
| 178 |  |  | varset("NxP`", '=', vec[0]*sf); | 
| 179 |  |  | varset("NyP`", '=', vec[1]*sf); | 
| 180 |  |  | varset("NzP`", '=', vec[2]*sf); | 
| 181 |  |  | return(1); | 
| 182 |  |  | } | 
| 183 |  |  |  | 
| 184 | greg | 2.1 | /* assign indicated diffuse component (do !trans first) */ | 
| 185 |  |  | static void | 
| 186 |  |  | set_dcomp(WGMDDAT *wp, int trans) | 
| 187 |  |  | { | 
| 188 |  |  | DCOMP           *dp = trans ? &wp->td : &wp->rd; | 
| 189 |  |  | const int       offs = trans ? 6 : 3*(wp->rp->rod < 0); | 
| 190 |  |  |  | 
| 191 |  |  | if (trans) {                    /* transmitted diffuse? */ | 
| 192 |  |  | if (intens(wp->mtp->oargs.farg+offs) <= FTINY) { | 
| 193 |  |  | scolorblack(dp->scol); | 
| 194 |  |  | return; | 
| 195 |  |  | } | 
| 196 |  |  | dp->mo.nam = wp->mtp->oargs.sarg[8]; | 
| 197 |  |  | if (!fill_modval(&dp->mo, wp)) { | 
| 198 |  |  | sprintf(errmsg, | 
| 199 |  |  | "unknown diffuse transmission modifier '%s'", | 
| 200 |  |  | dp->mo.nam); | 
| 201 |  |  | objerror(wp->mtp, USER, errmsg); | 
| 202 |  |  | } | 
| 203 |  |  | } else                          /* no priors for main mod */ | 
| 204 |  |  | fill_modval(&dp->mo, wp); | 
| 205 |  |  |  | 
| 206 |  |  | setscolor(dp->scol, wp->mtp->oargs.farg[offs], | 
| 207 |  |  | wp->mtp->oargs.farg[offs+1], | 
| 208 |  |  | wp->mtp->oargs.farg[offs+2]); | 
| 209 |  |  | smultscolor(dp->scol, dp->mo.pcol); | 
| 210 |  |  | } | 
| 211 |  |  |  | 
| 212 |  |  | /* assign indicated specular component */ | 
| 213 |  |  | static void | 
| 214 |  |  | set_scomp(WGMDDAT *wp, int trans) | 
| 215 |  |  | { | 
| 216 | greg | 2.6 | SCOMP   *sp = trans ? &wp->ts : &wp->rs; | 
| 217 |  |  | EPNODE  **exa = wp->mf->ep + 3*(trans != 0); | 
| 218 |  |  | double  coef; | 
| 219 | greg | 2.5 | /* constant zero check */ | 
| 220 | greg | 2.7 | if (exa[0]->type == NUM && exa[0]->v.num <= FTINY) | 
| 221 |  |  | goto blackout; | 
| 222 |  |  | /* need modifier */ | 
| 223 | greg | 2.5 | sp->mo.nam = wp->mtp->oargs.sarg[4*(trans != 0)]; | 
| 224 |  |  | if (!fill_modval(&sp->mo, wp)) { | 
| 225 |  |  | sprintf(errmsg, "unknown specular %s modifier '%s'", | 
| 226 |  |  | trans ? "transmission" : "reflection", sp->mo.nam); | 
| 227 |  |  | objerror(wp->mtp, USER, errmsg); | 
| 228 |  |  | } | 
| 229 | greg | 2.7 | if (sintens(sp->mo.pcol) <= FTINY) | 
| 230 |  |  | goto blackout;          /* got black pattern */ | 
| 231 |  |  | setWGMDfunc(&sp->mo, wp);       /* else compute coefficient */ | 
| 232 | greg | 2.1 | errno = 0; | 
| 233 | greg | 2.6 | coef = evalue(exa[0]); | 
| 234 | greg | 2.1 | if ((errno == EDOM) | (errno == ERANGE)) { | 
| 235 |  |  | objerror(wp->mtp, WARNING, "specular compute error"); | 
| 236 | greg | 2.7 | goto blackout; | 
| 237 | greg | 2.1 | } | 
| 238 | greg | 2.7 | if (coef <= FTINY)              /* negligible value? */ | 
| 239 |  |  | goto blackout; | 
| 240 | greg | 2.1 | copyscolor(sp->scol, sp->mo.pcol); | 
| 241 |  |  | scalescolor(sp->scol, coef); | 
| 242 | greg | 2.5 | errno = 0;                      /* else get roughness */ | 
| 243 | greg | 2.6 | sp->u_alpha = evalue(exa[1]); | 
| 244 |  |  | sp->v_alpha = (sp->u_alpha > FTINY) ? evalue(exa[2]) : 0.0; | 
| 245 | greg | 2.1 | if ((errno == EDOM) | (errno == ERANGE)) { | 
| 246 |  |  | objerror(wp->mtp, WARNING, "roughness compute error"); | 
| 247 | greg | 2.7 | goto blackout; | 
| 248 | greg | 2.1 | }                               /* we have something... */ | 
| 249 |  |  | wp->specfl |= trans ? SP_TRAN : SP_REFL; | 
| 250 |  |  | if (sp->v_alpha <= FTINY) {     /* is it pure specular? */ | 
| 251 |  |  | wp->specfl |= trans ? SP_TPURE : SP_RPURE; | 
| 252 |  |  | sp->u_alpha = sp->v_alpha = 0.0; | 
| 253 |  |  | return; | 
| 254 | greg | 2.7 | }                               /* else get aniso coordinates */ | 
| 255 | greg | 2.1 | fcross(sp->v, sp->mo.pnorm, wp->ulocal); | 
| 256 |  |  | if (normalize(sp->v) == 0.0) {  /* orientation vector==normal? */ | 
| 257 |  |  | if (fabs(sp->u_alpha - sp->v_alpha) > 0.001) | 
| 258 |  |  | objerror(wp->mtp, WARNING, "bad orientation vector"); | 
| 259 |  |  | getperpendicular(sp->u, sp->mo.pnorm, 1);       /* punting */ | 
| 260 |  |  | fcross(sp->v, sp->mo.pnorm, sp->u); | 
| 261 |  |  | sp->u_alpha = sp->v_alpha = sqrt( 0.5 * | 
| 262 |  |  | (sp->u_alpha*sp->u_alpha + sp->v_alpha*sp->v_alpha) ); | 
| 263 |  |  | } else | 
| 264 |  |  | fcross(sp->u, sp->v, sp->mo.pnorm); | 
| 265 | greg | 2.7 | return; | 
| 266 |  |  | blackout: | 
| 267 |  |  | scolorblack(sp->scol);          /* zero out component */ | 
| 268 | greg | 2.1 | } | 
| 269 |  |  |  | 
| 270 |  |  | /* sample anisotropic Gaussian specular */ | 
| 271 |  |  | static void | 
| 272 |  |  | agaussamp(WGMDDAT *wp) | 
| 273 |  |  | { | 
| 274 |  |  | RAY     sr; | 
| 275 |  |  | FVECT   h; | 
| 276 |  |  | double  rv[2]; | 
| 277 |  |  | double  d, sinp, cosp; | 
| 278 |  |  | int     maxiter, ntrials, nstarget, nstaken; | 
| 279 |  |  | int     i; | 
| 280 |  |  | /* compute reflection */ | 
| 281 |  |  | if ((wp->specfl & (SP_REFL|SP_RPURE|SP_RBLT)) == SP_REFL && | 
| 282 |  |  | rayorigin(&sr, RSPECULAR, wp->rp, wp->rs.scol) == 0) { | 
| 283 |  |  | SCOLOR  scol; | 
| 284 |  |  | nstarget = 1; | 
| 285 |  |  | if (specjitter > 1.5) { /* multiple samples? */ | 
| 286 |  |  | nstarget = specjitter*wp->rp->rweight + .5; | 
| 287 |  |  | if (sr.rweight <= minweight*nstarget) | 
| 288 |  |  | nstarget = sr.rweight/minweight; | 
| 289 |  |  | if (nstarget > 1) { | 
| 290 |  |  | d = 1./nstarget; | 
| 291 |  |  | scalescolor(sr.rcoef, d); | 
| 292 |  |  | sr.rweight *= d; | 
| 293 |  |  | } else | 
| 294 |  |  | nstarget = 1; | 
| 295 |  |  | } | 
| 296 |  |  | scolorblack(scol); | 
| 297 |  |  | dimlist[ndims++] = (int)(size_t)wp->mtp; | 
| 298 |  |  | maxiter = MAXITER*nstarget; | 
| 299 |  |  | for (nstaken = ntrials = 0; (nstaken < nstarget) & | 
| 300 |  |  | (ntrials < maxiter); ntrials++) { | 
| 301 |  |  | if (ntrials) | 
| 302 |  |  | d = frandom(); | 
| 303 |  |  | else | 
| 304 |  |  | d = urand(ilhash(dimlist,ndims)+samplendx); | 
| 305 |  |  | multisamp(rv, 2, d); | 
| 306 |  |  | d = 2.0*PI * rv[0]; | 
| 307 |  |  | cosp = tcos(d) * wp->rs.u_alpha; | 
| 308 |  |  | sinp = tsin(d) * wp->rs.v_alpha; | 
| 309 |  |  | d = 1./sqrt(cosp*cosp + sinp*sinp); | 
| 310 |  |  | cosp *= d; | 
| 311 |  |  | sinp *= d; | 
| 312 |  |  | if ((0. <= specjitter) & (specjitter < 1.)) | 
| 313 |  |  | rv[1] = 1.0 - specjitter*rv[1]; | 
| 314 |  |  | d = (rv[1] <= FTINY) ? 1.0 : sqrt( -log(rv[1]) / | 
| 315 |  |  | (cosp*cosp/(wp->rs.u_alpha*wp->rs.u_alpha) + | 
| 316 |  |  | sinp*sinp/(wp->rs.v_alpha*wp->rs.v_alpha)) ); | 
| 317 |  |  | for (i = 0; i < 3; i++) | 
| 318 |  |  | h[i] = wp->rs.mo.pnorm[i] + | 
| 319 |  |  | d*(cosp*wp->rs.u[i] + sinp*wp->rs.v[i]); | 
| 320 |  |  | d = -2.0 * DOT(h, wp->rp->rdir) / (1.0 + d*d); | 
| 321 |  |  | VSUM(sr.rdir, wp->rp->rdir, h, d); | 
| 322 |  |  | /* sample rejection test */ | 
| 323 |  |  | d = DOT(sr.rdir, wp->rp->ron); | 
| 324 |  |  | if ((d > 0) ^ (wp->rp->rod > 0)) | 
| 325 |  |  | continue; | 
| 326 |  |  | checknorm(sr.rdir); | 
| 327 |  |  | if (nstarget > 1) {     /* W-G-M-D adjustment */ | 
| 328 |  |  | if (nstaken) rayclear(&sr); | 
| 329 |  |  | rayvalue(&sr); | 
| 330 |  |  | d = 2./(1. + wp->rp->rod/d); | 
| 331 |  |  | scalescolor(sr.rcol, d); | 
| 332 |  |  | saddscolor(scol, sr.rcol); | 
| 333 |  |  | } else { | 
| 334 |  |  | rayvalue(&sr); | 
| 335 |  |  | smultscolor(sr.rcol, sr.rcoef); | 
| 336 |  |  | saddscolor(wp->rp->rcol, sr.rcol); | 
| 337 |  |  | } | 
| 338 |  |  | ++nstaken; | 
| 339 |  |  | } | 
| 340 |  |  | if (nstarget > 1) {             /* final W-G-M-D weighting */ | 
| 341 |  |  | smultscolor(scol, sr.rcoef); | 
| 342 |  |  | d = (double)nstarget/ntrials; | 
| 343 |  |  | scalescolor(scol, d); | 
| 344 |  |  | saddscolor(wp->rp->rcol, scol); | 
| 345 |  |  | } | 
| 346 |  |  | ndims--; | 
| 347 |  |  | } | 
| 348 |  |  | /* compute transmission */ | 
| 349 |  |  | if ((wp->specfl & (SP_TRAN|SP_TPURE|SP_TBLT)) == SP_TRAN && | 
| 350 |  |  | rayorigin(&sr, TSPECULAR, wp->rp, wp->ts.scol) == 0) { | 
| 351 |  |  | nstarget = 1; | 
| 352 |  |  | if (specjitter > 1.5) { /* multiple samples? */ | 
| 353 |  |  | nstarget = specjitter*wp->rp->rweight + .5; | 
| 354 |  |  | if (sr.rweight <= minweight*nstarget) | 
| 355 |  |  | nstarget = sr.rweight/minweight; | 
| 356 |  |  | if (nstarget > 1) { | 
| 357 |  |  | d = 1./nstarget; | 
| 358 |  |  | scalescolor(sr.rcoef, d); | 
| 359 |  |  | sr.rweight *= d; | 
| 360 |  |  | } else | 
| 361 |  |  | nstarget = 1; | 
| 362 |  |  | } | 
| 363 |  |  | dimlist[ndims++] = (int)(size_t)wp->mtp; | 
| 364 |  |  | maxiter = MAXITER*nstarget; | 
| 365 |  |  | for (nstaken = ntrials = 0; (nstaken < nstarget) & | 
| 366 |  |  | (ntrials < maxiter); ntrials++) { | 
| 367 |  |  | if (ntrials) | 
| 368 |  |  | d = frandom(); | 
| 369 |  |  | else | 
| 370 |  |  | d = urand(ilhash(dimlist,ndims)+1823+samplendx); | 
| 371 |  |  | multisamp(rv, 2, d); | 
| 372 |  |  | d = 2.0*PI * rv[0]; | 
| 373 |  |  | cosp = tcos(d) * wp->ts.u_alpha; | 
| 374 |  |  | sinp = tsin(d) * wp->ts.v_alpha; | 
| 375 |  |  | d = 1./sqrt(cosp*cosp + sinp*sinp); | 
| 376 |  |  | cosp *= d; | 
| 377 |  |  | sinp *= d; | 
| 378 |  |  | if ((0. <= specjitter) & (specjitter < 1.)) | 
| 379 |  |  | rv[1] = 1.0 - specjitter*rv[1]; | 
| 380 |  |  | if (rv[1] <= FTINY) | 
| 381 |  |  | d = 1.0; | 
| 382 |  |  | else | 
| 383 |  |  | d = sqrt(-log(rv[1]) / | 
| 384 |  |  | (cosp*cosp/(wp->ts.u_alpha*wp->ts.u_alpha) + | 
| 385 |  |  | sinp*sinp/(wp->ts.v_alpha*wp->ts.v_alpha))); | 
| 386 |  |  | for (i = 0; i < 3; i++) | 
| 387 |  |  | sr.rdir[i] = wp->prdir[i] + | 
| 388 |  |  | d*(cosp*wp->ts.u[i] + sinp*wp->ts.v[i]); | 
| 389 |  |  | /* rejection test */ | 
| 390 |  |  | if ((DOT(sr.rdir,wp->rp->ron) > 0) == (wp->rp->rod > 0)) | 
| 391 |  |  | continue; | 
| 392 |  |  | normalize(sr.rdir);     /* OK, normalize */ | 
| 393 |  |  | if (nstaken)            /* multi-sampling? */ | 
| 394 |  |  | rayclear(&sr); | 
| 395 |  |  | rayvalue(&sr); | 
| 396 |  |  | smultscolor(sr.rcol, sr.rcoef); | 
| 397 |  |  | saddscolor(wp->rp->rcol, sr.rcol); | 
| 398 |  |  | ++nstaken; | 
| 399 |  |  | } | 
| 400 |  |  | ndims--; | 
| 401 |  |  | } | 
| 402 |  |  | } | 
| 403 |  |  |  | 
| 404 |  |  | /* compute source contribution for MAT_WGMDF */ | 
| 405 |  |  | static void | 
| 406 |  |  | dirwgmdf(SCOLOR scval, void *uwp, FVECT ldir, double omega) | 
| 407 |  |  | { | 
| 408 |  |  | WGMDDAT         *wp = (WGMDDAT *)uwp; | 
| 409 |  |  | const int       hitfront = (wp->rp->rod > 0); | 
| 410 |  |  | double          fresadj = 1.; | 
| 411 |  |  | double          ldot; | 
| 412 |  |  | double          dtmp, dtmp1, dtmp2; | 
| 413 |  |  | FVECT           h; | 
| 414 |  |  | double          au2, av2; | 
| 415 |  |  | SCOLOR          sctmp; | 
| 416 |  |  |  | 
| 417 |  |  | scolorblack(scval);             /* will add component coefficients */ | 
| 418 |  |  |  | 
| 419 |  |  | /* XXX ignores which side is lit */ | 
| 420 |  |  | if (wp->specfl & SP_RPURE && pbright(wp->rs.scol) >= FRESTHRESH) | 
| 421 |  |  | fresadj = 1. - FRESNE(fabs(DOT(wp->rs.mo.pnorm,ldir))); | 
| 422 |  |  |  | 
| 423 |  |  | if (sintens(wp->rd.scol) > FTINY && | 
| 424 |  |  | ((ldot = DOT(wp->rd.mo.pnorm,ldir)) > 0) == hitfront) { | 
| 425 |  |  | /* | 
| 426 |  |  | *  Compute diffuse reflection coefficient for source. | 
| 427 |  |  | */ | 
| 428 |  |  | copyscolor(sctmp, wp->rd.scol); | 
| 429 |  |  | dtmp = fabs(ldot) * omega * (1.0/PI) * fresadj; | 
| 430 |  |  | scalescolor(sctmp, dtmp); | 
| 431 |  |  | saddscolor(scval, sctmp); | 
| 432 |  |  | } | 
| 433 |  |  | if (sintens(wp->td.scol) > FTINY && | 
| 434 |  |  | ((ldot = DOT(wp->td.mo.pnorm,ldir)) > 0) ^ hitfront) { | 
| 435 |  |  | /* | 
| 436 |  |  | *  Compute diffuse transmission coefficient for source. | 
| 437 |  |  | */ | 
| 438 |  |  | copyscolor(sctmp, wp->td.scol); | 
| 439 |  |  | dtmp = fabs(ldot) * omega * (1.0/PI) * fresadj; | 
| 440 |  |  | scalescolor(sctmp, dtmp); | 
| 441 |  |  | saddscolor(scval, sctmp); | 
| 442 |  |  | } | 
| 443 |  |  | #if 0   /* XXX not yet implemented */ | 
| 444 |  |  | if (ambRayInPmap(wp->rp)) | 
| 445 |  |  | return;         /* specular accounted for in photon map */ | 
| 446 |  |  | #endif | 
| 447 |  |  | if ((wp->specfl & (SP_REFL|SP_RPURE)) == SP_REFL && | 
| 448 |  |  | ((ldot = DOT(wp->rs.mo.pnorm,ldir)) > 0) == hitfront) { | 
| 449 |  |  | /* | 
| 450 |  |  | *  Compute specular reflection coefficient for source using | 
| 451 |  |  | *  anisotropic Gaussian distribution model. | 
| 452 |  |  | */ | 
| 453 |  |  | /* add source width if flat */ | 
| 454 |  |  | if (wp->specfl & SP_FLAT) | 
| 455 |  |  | au2 = av2 = omega * (0.25/PI); | 
| 456 |  |  | else | 
| 457 |  |  | au2 = av2 = 0.0; | 
| 458 |  |  | au2 += wp->rs.u_alpha*wp->rs.u_alpha; | 
| 459 |  |  | av2 += wp->rs.v_alpha*wp->rs.v_alpha; | 
| 460 |  |  | /* half vector */ | 
| 461 |  |  | VSUB(h, ldir, wp->rp->rdir); | 
| 462 |  |  | /* ellipse */ | 
| 463 |  |  | dtmp1 = DOT(wp->rs.u, h); | 
| 464 |  |  | dtmp1 *= dtmp1 / au2; | 
| 465 |  |  | dtmp2 = DOT(wp->rs.v, h); | 
| 466 |  |  | dtmp2 *= dtmp2 / av2; | 
| 467 |  |  | /* W-G-M-D model */ | 
| 468 |  |  | dtmp = DOT(wp->rs.mo.pnorm, h); | 
| 469 |  |  | dtmp *= dtmp; | 
| 470 |  |  | dtmp1 = (dtmp1 + dtmp2) / dtmp; | 
| 471 |  |  | dtmp = exp(-dtmp1) * DOT(h,h) / | 
| 472 |  |  | (PI * dtmp*dtmp * sqrt(au2*av2)); | 
| 473 |  |  |  | 
| 474 |  |  | if (dtmp > FTINY) {             /* worth using? */ | 
| 475 |  |  | copyscolor(sctmp, wp->rs.scol); | 
| 476 |  |  | dtmp *= fabs(ldot) * omega; | 
| 477 |  |  | scalescolor(sctmp, dtmp); | 
| 478 |  |  | saddscolor(scval, sctmp); | 
| 479 |  |  | } | 
| 480 |  |  | } | 
| 481 |  |  | if ((wp->specfl & (SP_TRAN|SP_TPURE)) == SP_TRAN && | 
| 482 |  |  | ((ldot = DOT(wp->ts.mo.pnorm,ldir)) > 0) ^ hitfront) { | 
| 483 |  |  | /* | 
| 484 |  |  | *  Compute specular transmission coefficient for source. | 
| 485 |  |  | */ | 
| 486 |  |  | /* roughness + source */ | 
| 487 |  |  | au2 = av2 = omega * (1.0/PI); | 
| 488 |  |  | au2 += wp->ts.u_alpha*wp->ts.u_alpha; | 
| 489 |  |  | av2 += wp->ts.v_alpha*wp->ts.v_alpha; | 
| 490 |  |  | /* "half vector" */ | 
| 491 |  |  | VSUB(h, ldir, wp->prdir); | 
| 492 |  |  | dtmp = DOT(h,h); | 
| 493 |  |  | if (dtmp > FTINY*FTINY) { | 
| 494 |  |  | dtmp1 = DOT(h,wp->ts.mo.pnorm); | 
| 495 |  |  | dtmp = 1.0 - dtmp1*dtmp1/dtmp; | 
| 496 |  |  | } | 
| 497 |  |  | if (dtmp > FTINY*FTINY) { | 
| 498 |  |  | dtmp1 = DOT(h,wp->ts.u); | 
| 499 |  |  | dtmp1 *= dtmp1 / au2; | 
| 500 |  |  | dtmp2 = DOT(h,wp->ts.v); | 
| 501 |  |  | dtmp2 *= dtmp2 / av2; | 
| 502 |  |  | dtmp = (dtmp1 + dtmp2) / dtmp; | 
| 503 |  |  | dtmp = exp(-dtmp); | 
| 504 |  |  | } else | 
| 505 |  |  | dtmp = 1.0; | 
| 506 |  |  | /* Gaussian */ | 
| 507 |  |  | dtmp *= (1.0/PI) * sqrt(-ldot/(wp->ts.mo.pdot*au2*av2)); | 
| 508 |  |  |  | 
| 509 |  |  | if (dtmp > FTINY) {             /* worth using? */ | 
| 510 |  |  | copyscolor(sctmp, wp->ts.scol); | 
| 511 |  |  | dtmp *= omega; | 
| 512 |  |  | scalescolor(sctmp, dtmp); | 
| 513 |  |  | saddscolor(scval, sctmp); | 
| 514 |  |  | } | 
| 515 |  |  | } | 
| 516 |  |  | } | 
| 517 |  |  |  | 
| 518 |  |  | /* color a ray that hit a programmable WGMD material */ | 
| 519 |  |  | int | 
| 520 |  |  | m_wgmdf(OBJREC *m, RAY *r) | 
| 521 |  |  | { | 
| 522 |  |  | RAY             lr; | 
| 523 |  |  | WGMDDAT         wd; | 
| 524 |  |  | SCOLOR          sctmp; | 
| 525 |  |  | FVECT           anorm; | 
| 526 |  |  | int             i; | 
| 527 |  |  |  | 
| 528 |  |  | if (!backvis & (r->rod < 0.0)) { | 
| 529 |  |  | raytrans(r); | 
| 530 |  |  | return(1);              /* backside invisible */ | 
| 531 |  |  | } | 
| 532 |  |  | if ((m->oargs.nsargs < 13) | (m->oargs.nfargs < 9)) | 
| 533 |  |  | objerror(m, USER, "bad number of arguments"); | 
| 534 | greg | 2.2 |  | 
| 535 |  |  | if (r->crtype & SHADOW && !strcmp(m->oargs.sarg[5], "0")) | 
| 536 |  |  | return(1);              /* first shadow test */ | 
| 537 | greg | 2.1 | clr_comps(&wd); | 
| 538 |  |  | wd.rp = r; | 
| 539 |  |  | wd.mtp = m; | 
| 540 |  |  | wd.mf = getfunc(m, 12, 0xEEE, 1); | 
| 541 | greg | 2.5 | set_dcomp(&wd, 0);              /* gets main modifier */ | 
| 542 |  |  | setWGMDfunc(&wd.rd.mo, &wd);    /* get local u vector */ | 
| 543 | greg | 2.1 | errno = 0; | 
| 544 |  |  | for (i = 0; i < 3; i++) | 
| 545 |  |  | wd.ulocal[i] = evalue(wd.mf->ep[6+i]); | 
| 546 |  |  | if ((errno == EDOM) | (errno == ERANGE)) | 
| 547 |  |  | wd.ulocal[0] = wd.ulocal[1] = wd.ulocal[2] = 0.0; | 
| 548 |  |  | else if (wd.mf->fxp != &unitxf) | 
| 549 |  |  | multv3(wd.ulocal, wd.ulocal, wd.mf->fxp->xfm); | 
| 550 |  |  |  | 
| 551 |  |  | set_scomp(&wd, 1);              /* sets SP_TPURE */ | 
| 552 |  |  | if (r->crtype & SHADOW && !(wd.specfl & SP_TPURE)) | 
| 553 | greg | 2.2 | return(1);              /* second shadow test */ | 
| 554 |  |  | set_dcomp(&wd, 1); | 
| 555 | greg | 2.1 | set_scomp(&wd, 0); | 
| 556 |  |  | wd.specfl |= SP_FLAT*(r->ro != NULL && isflat(r->ro->otype)); | 
| 557 |  |  | /* apply Fresnel adjustments? */ | 
| 558 |  |  | if (wd.specfl & SP_RPURE && pbright(wd.rs.scol) >= FRESTHRESH) { | 
| 559 |  |  | const double    fest = FRESNE(fabs(wd.rs.mo.pdot)); | 
| 560 |  |  | for (i = NCSAMP; i--; ) | 
| 561 |  |  | wd.rs.scol[i] += fest*(1. - wd.rs.scol[i]); | 
| 562 | greg | 2.4 | scalescolor(wd.rd.scol, 1.-fest); | 
| 563 | greg | 2.1 | scalescolor(wd.ts.scol, 1.-fest); | 
| 564 |  |  | scalescolor(wd.td.scol, 1.-fest); | 
| 565 |  |  | } | 
| 566 |  |  | /* check specular thresholds */ | 
| 567 |  |  | wd.specfl |= SP_RBLT*((wd.specfl & (SP_REFL|SP_RPURE)) == SP_REFL && | 
| 568 |  |  | specthresh >= pbright(wd.rs.scol)-FTINY); | 
| 569 |  |  | wd.specfl |= SP_TBLT*((wd.specfl & (SP_TRAN|SP_TPURE)) == SP_TRAN && | 
| 570 |  |  | specthresh >= pbright(wd.ts.scol)-FTINY); | 
| 571 |  |  | /* get through direction */ | 
| 572 |  |  | if (wd.specfl & SP_TRAN && wd.ts.mo.hastexture && | 
| 573 |  |  | !(r->crtype & (SHADOW|AMBIENT))) { | 
| 574 |  |  | for (i = 0; i < 3; i++) /* perturb */ | 
| 575 |  |  | wd.prdir[i] = r->rdir[i] - wd.ts.mo.pnorm[i] + r->ron[i]; | 
| 576 |  |  | if ((DOT(wd.prdir,r->ron) > 0) ^ (r->rod > 0)) | 
| 577 |  |  | normalize(wd.prdir);    /* OK */ | 
| 578 |  |  | else                    /* too much */ | 
| 579 |  |  | VCOPY(wd.prdir, r->rdir); | 
| 580 |  |  | } else | 
| 581 |  |  | VCOPY(wd.prdir, r->rdir); | 
| 582 |  |  | /* transmitted view ray? */ | 
| 583 |  |  | if ((wd.specfl & (SP_TRAN|SP_TPURE|SP_TBLT)) == (SP_TRAN|SP_TPURE) && | 
| 584 |  |  | rayorigin(&lr, TRANS, r, wd.ts.scol) == 0) { | 
| 585 |  |  | VCOPY(lr.rdir, wd.prdir); | 
| 586 |  |  | rayvalue(&lr); | 
| 587 |  |  | smultscolor(lr.rcol, lr.rcoef); | 
| 588 |  |  | saddscolor(r->rcol, lr.rcol); | 
| 589 |  |  | if (scolor_mean(wd.ts.scol) >= 0.999) { | 
| 590 |  |  | /* completely transparent */ | 
| 591 |  |  | smultscolor(lr.mcol, lr.rcoef); | 
| 592 |  |  | copyscolor(r->mcol, lr.mcol); | 
| 593 |  |  | r->rmt = r->rot + lr.rmt; | 
| 594 |  |  | r->rxt = r->rot + lr.rxt; | 
| 595 |  |  | } else if (pbright(wd.ts.scol) > | 
| 596 |  |  | pbright(wd.td.scol) + pbright(wd.rd.scol)) | 
| 597 |  |  | r->rxt = r->rot + raydistance(&lr); | 
| 598 |  |  | } | 
| 599 |  |  | if (r->crtype & SHADOW) | 
| 600 |  |  | return(1);              /* the rest is shadow */ | 
| 601 |  |  | /* mirror ray? */ | 
| 602 |  |  | if ((wd.specfl & (SP_REFL|SP_RPURE|SP_RBLT)) == (SP_REFL|SP_RPURE) && | 
| 603 |  |  | rayorigin(&lr, REFLECTED, r, wd.rs.scol) == 0) { | 
| 604 |  |  | VSUM(lr.rdir, r->rdir, wd.rs.mo.pnorm, 2.*wd.rs.mo.pdot); | 
| 605 |  |  | /* fall back if would penetrate */ | 
| 606 |  |  | if (wd.rs.mo.hastexture && | 
| 607 |  |  | (DOT(lr.rdir,r->ron) > 0) ^ (r->rod > 0)) | 
| 608 |  |  | VSUM(lr.rdir, r->rdir, r->ron, 2.*r->rod); | 
| 609 |  |  | checknorm(lr.rdir); | 
| 610 |  |  | rayvalue(&lr); | 
| 611 |  |  | smultscolor(lr.rcol, lr.rcoef); | 
| 612 |  |  | copyscolor(r->mcol, lr.rcol); | 
| 613 |  |  | saddscolor(r->rcol, lr.rcol); | 
| 614 |  |  | r->rmt = r->rot; | 
| 615 |  |  | if (wd.specfl & SP_FLAT && | 
| 616 |  |  | !wd.rs.mo.hastexture | (r->crtype & AMBIENT)) | 
| 617 |  |  | r->rmt += raydistance(&lr); | 
| 618 |  |  | } | 
| 619 |  |  | if (wd.specfl & (SP_REFL|SP_TRAN))      /* specularly scattered rays */ | 
| 620 |  |  | agaussamp(&wd);         /* checks *BLT flags */ | 
| 621 |  |  |  | 
| 622 |  |  | if (sintens(wd.rd.scol) > FTINY) {      /* ambient from this side */ | 
| 623 |  |  | if (r->rod > 0) { | 
| 624 |  |  | VCOPY(anorm, wd.rd.mo.pnorm); | 
| 625 |  |  | } else { | 
| 626 |  |  | anorm[0] = -wd.rd.mo.pnorm[0]; | 
| 627 |  |  | anorm[1] = -wd.rd.mo.pnorm[1]; | 
| 628 |  |  | anorm[2] = -wd.rd.mo.pnorm[2]; | 
| 629 |  |  | } | 
| 630 |  |  | copyscolor(sctmp, wd.rd.scol); | 
| 631 |  |  | if (wd.specfl & SP_RBLT)        /* add in specular as well? */ | 
| 632 |  |  | saddscolor(sctmp, wd.rs.scol); | 
| 633 |  |  | multambient(sctmp, r, anorm); | 
| 634 |  |  | saddscolor(r->rcol, sctmp);     /* add to returned color */ | 
| 635 |  |  | } | 
| 636 |  |  | if (sintens(wd.td.scol) > FTINY) {      /* ambient from other side */ | 
| 637 |  |  | if (r->rod > 0) { | 
| 638 |  |  | anorm[0] = -wd.td.mo.pnorm[0]; | 
| 639 |  |  | anorm[1] = -wd.td.mo.pnorm[1]; | 
| 640 |  |  | anorm[2] = -wd.td.mo.pnorm[2]; | 
| 641 |  |  | } else { | 
| 642 |  |  | VCOPY(anorm, wd.td.mo.pnorm); | 
| 643 |  |  | } | 
| 644 |  |  | copyscolor(sctmp, wd.td.scol); | 
| 645 |  |  | if (wd.specfl & SP_TBLT)        /* add in specular as well? */ | 
| 646 |  |  | saddscolor(sctmp, wd.ts.scol) | 
| 647 |  |  | multambient(sctmp, r, anorm); | 
| 648 |  |  | saddscolor(r->rcol, sctmp); | 
| 649 |  |  | } | 
| 650 |  |  | direct(r, dirwgmdf, &wd);       /* add direct component last */ | 
| 651 |  |  | return(1); | 
| 652 |  |  | } |