00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "reflection.h"
00026 #include "color.h"
00027 #include "mc.h"
00028 #include "sampling.h"
00029 #include <stdarg.h>
00030
00031 COREDLL Spectrum FrDiel(float cosi, float cost,
00032 const Spectrum &etai,
00033 const Spectrum &etat) {
00034 Spectrum Rparl = ((etat * cosi) - (etai * cost)) /
00035 ((etat * cosi) + (etai * cost));
00036 Spectrum Rperp = ((etai * cosi) - (etat * cost)) /
00037 ((etai * cosi) + (etat * cost));
00038 return (Rparl*Rparl + Rperp*Rperp) / 2.f;
00039 }
00040 COREDLL Spectrum FrCond(float cosi,
00041 const Spectrum &eta,
00042 const Spectrum &k) {
00043 Spectrum tmp = (eta*eta + k*k) * cosi*cosi;
00044 Spectrum Rparl2 = (tmp - (2.f * eta * cosi) + 1) /
00045 (tmp + (2.f * eta * cosi) + 1);
00046 Spectrum tmp_f = eta*eta + k*k;
00047 Spectrum Rperp2 =
00048 (tmp_f - (2.f * eta * cosi) + cosi*cosi) /
00049 (tmp_f + (2.f * eta * cosi) + cosi*cosi);
00050 return (Rparl2 + Rperp2) / 2.f;
00051 }
00052 COREDLL Spectrum FresnelApproxEta(const Spectrum &Fr) {
00053 Spectrum reflectance = Fr.Clamp(0.f, .999f);
00054 return (Spectrum(1.) + reflectance.Sqrt()) /
00055 (Spectrum(1.) - reflectance.Sqrt());
00056 }
00057 COREDLL Spectrum FresnelApproxK(const Spectrum &Fr) {
00058 Spectrum reflectance = Fr.Clamp(0.f, .999f);
00059 return 2.f * (reflectance /
00060 (Spectrum(1.) - reflectance)).Sqrt();
00061 }
00062
00063 Spectrum BRDFToBTDF::f(const Vector &wo,
00064 const Vector &wi) const {
00065 return brdf->f(wo, otherHemisphere(wi));
00066 }
00067 Spectrum BRDFToBTDF::Sample_f(const Vector &wo, Vector *wi,
00068 float u1, float u2, float *pdf) const {
00069 Spectrum f = brdf->Sample_f(wo, wi, u1, u2, pdf);
00070 *wi = otherHemisphere(*wi);
00071 return f;
00072 }
00073 Fresnel::~Fresnel() { }
00074 Spectrum FresnelConductor::Evaluate(float cosi) const {
00075 return FrCond(fabsf(cosi), eta, k);
00076 }
00077 Spectrum FresnelDielectric::Evaluate(float cosi) const {
00078
00079 cosi = Clamp(cosi, -1.f, 1.f);
00080
00081 bool entering = cosi > 0.;
00082 float ei = eta_i, et = eta_t;
00083 if (!entering)
00084 swap(ei, et);
00085
00086 float sint = ei/et * sqrtf(max(0.f, 1.f - cosi*cosi));
00087 if (sint > 1.) {
00088
00089 return 1.;
00090 }
00091 else {
00092 float cost = sqrtf(max(0.f, 1.f - sint*sint));
00093 return FrDiel(fabsf(cosi), cost, ei, et);
00094 }
00095 }
00096 Spectrum SpecularReflection::Sample_f(const Vector &wo,
00097 Vector *wi, float u1, float u2, float *pdf) const {
00098
00099 *wi = Vector(-wo.x, -wo.y, wo.z);
00100 *pdf = 1.f;
00101 return fresnel->Evaluate(CosTheta(wo)) * R /
00102 fabsf(CosTheta(*wi));
00103 }
00104 Spectrum SpecularTransmission::Sample_f(const Vector &wo,
00105 Vector *wi, float u1, float u2, float *pdf) const {
00106
00107 bool entering = CosTheta(wo) > 0.;
00108 float ei = etai, et = etat;
00109 if (!entering)
00110 swap(ei, et);
00111
00112 float sini2 = SinTheta2(wo);
00113 float eta = ei / et;
00114 float sint2 = eta * eta * sini2;
00115
00116 if (sint2 > 1.) return 0.;
00117 float cost = sqrtf(max(0.f, 1.f - sint2));
00118 if (entering) cost = -cost;
00119 float sintOverSini = eta;
00120 *wi = Vector(sintOverSini * -wo.x,
00121 sintOverSini * -wo.y,
00122 cost);
00123 *pdf = 1.f;
00124 Spectrum F = fresnel.Evaluate(CosTheta(wo));
00125 return (et*et)/(ei*ei) * (Spectrum(1.)-F) * T /
00126 fabsf(CosTheta(*wi));
00127 }
00128 Spectrum Lambertian::f(const Vector &wo,
00129 const Vector &wi) const {
00130 return RoverPI;
00131 }
00132 Spectrum OrenNayar::f(const Vector &wo,
00133 const Vector &wi) const {
00134 float sinthetai = SinTheta(wi);
00135 float sinthetao = SinTheta(wo);
00136
00137 float maxcos = 0.f;
00138 if (sinthetai > 1e-4 && sinthetao > 1e-4) {
00139 float sinphii = SinPhi(wi), cosphii = CosPhi(wi);
00140 float sinphio = SinPhi(wo), cosphio = CosPhi(wo);
00141 float dcos = cosphii * cosphio + sinphii * sinphio;
00142 maxcos = max(0.f, dcos);
00143 }
00144
00145 float sinalpha, tanbeta;
00146 if (fabsf(CosTheta(wi)) > fabsf(CosTheta(wo))) {
00147 sinalpha = sinthetao;
00148 tanbeta = sinthetai / fabsf(CosTheta(wi));
00149 }
00150 else {
00151 sinalpha = sinthetai;
00152 tanbeta = sinthetao / fabsf(CosTheta(wo));
00153 }
00154 return R * INV_PI *
00155 (A + B * maxcos * sinalpha * tanbeta);
00156 }
00157 Microfacet::Microfacet(const Spectrum &reflectance,
00158 Fresnel *f,
00159 MicrofacetDistribution *d)
00160 : BxDF(BxDFType(BSDF_REFLECTION | BSDF_GLOSSY)),
00161 R(reflectance), distribution(d), fresnel(f) {
00162 }
00163 Spectrum Microfacet::f(const Vector &wo,
00164 const Vector &wi) const {
00165 float cosThetaO = fabsf(CosTheta(wo));
00166 float cosThetaI = fabsf(CosTheta(wi));
00167 Vector wh = Normalize(wi + wo);
00168 float cosThetaH = Dot(wi, wh);
00169 Spectrum F = fresnel->Evaluate(cosThetaH);
00170 return R * distribution->D(wh) * G(wo, wi, wh) * F /
00171 (4.f * cosThetaI * cosThetaO);
00172 }
00173 Lafortune::Lafortune(const Spectrum &r, u_int nl,
00174 const Spectrum *xx,
00175 const Spectrum *yy,
00176 const Spectrum *zz,
00177 const Spectrum *e, BxDFType t)
00178 : BxDF(t), R(r) {
00179 nLobes = nl;
00180 x = xx;
00181 y = yy;
00182 z = zz;
00183 exponent = e;
00184 }
00185 Spectrum Lafortune::f(const Vector &wo,
00186 const Vector &wi) const {
00187 Spectrum ret = R * INV_PI;
00188 for (u_int i = 0; i < nLobes; ++i) {
00189
00190 Spectrum v = x[i] * wo.x * wi.x + y[i] * wo.y * wi.y +
00191 z[i] * wo.z * wi.z;
00192 ret += v.Pow(exponent[i]);
00193 }
00194 return ret;
00195 }
00196 FresnelBlend::FresnelBlend(const Spectrum &d,
00197 const Spectrum &s,
00198 MicrofacetDistribution *dist)
00199 : BxDF(BxDFType(BSDF_REFLECTION | BSDF_GLOSSY)),
00200 Rd(d), Rs(s) {
00201 distribution = dist;
00202 }
00203 Spectrum FresnelBlend::f(const Vector &wo,
00204 const Vector &wi) const {
00205 Spectrum diffuse = (28.f/(23.f*M_PI)) * Rd *
00206 (Spectrum(1.) - Rs) *
00207 (1 - powf(1 - .5f * fabsf(CosTheta(wi)), 5)) *
00208 (1 - powf(1 - .5f * fabsf(CosTheta(wo)), 5));
00209 Vector H = Normalize(wi + wo);
00210 Spectrum specular = distribution->D(H) /
00211 (4.f * AbsDot(wi, H) *
00212 max(fabsf(CosTheta(wi)), fabsf(CosTheta(wo)))) *
00213 SchlickFresnel(Dot(wi, H));
00214 return diffuse + specular;
00215 }
00216 Spectrum BxDF::Sample_f(const Vector &wo, Vector *wi,
00217 float u1, float u2, float *pdf) const {
00218
00219 *wi = CosineSampleHemisphere(u1, u2);
00220 if (wo.z < 0.) wi->z *= -1.f;
00221 *pdf = Pdf(wo, *wi);
00222 return f(wo, *wi);
00223 }
00224 float BxDF::Pdf(const Vector &wo, const Vector &wi) const {
00225 return
00226 SameHemisphere(wo, wi) ? fabsf(wi.z) * INV_PI : 0.f;
00227 }
00228 float BRDFToBTDF::Pdf(const Vector &wo,
00229 const Vector &wi) const {
00230 return brdf->Pdf(wo, -wi);
00231 }
00232 Spectrum Microfacet::Sample_f(const Vector &wo, Vector *wi,
00233 float u1, float u2, float *pdf) const {
00234 distribution->Sample_f(wo, wi, u1, u2, pdf);
00235 if (!SameHemisphere(wo, *wi)) return Spectrum(0.f);
00236 return f(wo, *wi);
00237 }
00238 float Microfacet::Pdf(const Vector &wo,
00239 const Vector &wi) const {
00240 if (!SameHemisphere(wo, wi)) return 0.f;
00241 return distribution->Pdf(wo, wi);
00242 }
00243 void Blinn::Sample_f(const Vector &wo, Vector *wi,
00244 float u1, float u2, float *pdf) const {
00245
00246 float costheta = powf(u1, 1.f / (exponent+1));
00247 float sintheta = sqrtf(max(0.f, 1.f - costheta*costheta));
00248 float phi = u2 * 2.f * M_PI;
00249 Vector H = SphericalDirection(sintheta, costheta, phi);
00250 if (!SameHemisphere(wo, H)) H.z *= -1.f;
00251
00252 *wi = -wo + 2.f * Dot(wo, H) * H;
00253
00254 float blinn_pdf = ((exponent + 2.f) *
00255 powf(costheta, exponent)) /
00256 (2.f * M_PI * 4.f * Dot(wo, H));
00257 *pdf = blinn_pdf;
00258 }
00259 float Blinn::Pdf(const Vector &wo, const Vector &wi) const {
00260 Vector H = Normalize(wo + wi);
00261 float costheta = fabsf(H.z);
00262
00263 float blinn_pdf = ((exponent + 2.f) *
00264 powf(costheta, exponent)) /
00265 (2.f * M_PI * 4.f * Dot(wo, H));
00266 return blinn_pdf;
00267 }
00268 void Anisotropic::Sample_f(const Vector &wo, Vector *wi,
00269 float u1, float u2, float *pdf) const {
00270
00271 float phi, costheta;
00272 if (u1 < .25f) {
00273 sampleFirstQuadrant(4.f * u1, u2, &phi, &costheta);
00274 } else if (u1 < .5f) {
00275 u1 = 4.f * (.5f - u1);
00276 sampleFirstQuadrant(u1, u2, &phi, &costheta);
00277 phi = M_PI - phi;
00278 } else if (u1 < .75f) {
00279 u1 = 4.f * (u1 - .5f);
00280 sampleFirstQuadrant(u1, u2, &phi, &costheta);
00281 phi += M_PI;
00282 } else {
00283 u1 = 4.f * (1.f - u1);
00284 sampleFirstQuadrant(u1, u2, &phi, &costheta);
00285 phi = 2.f * M_PI - phi;
00286 }
00287 float sintheta = sqrtf(max(0.f, 1.f - costheta*costheta));
00288 Vector H = SphericalDirection(sintheta, costheta, phi);
00289 if (Dot(wo, H) < 0.f) H = -H;
00290
00291 *wi = -wo + 2.f * Dot(wo, H) * H;
00292
00293 float anisotropic_pdf = D(H) / (4.f * Dot(wo, H));
00294 if (Dot(wo, H) < 0.f) anisotropic_pdf = 0.f;
00295 *pdf = anisotropic_pdf;
00296 }
00297 void Anisotropic::sampleFirstQuadrant(float u1, float u2,
00298 float *phi, float *costheta) const {
00299 if (ex == ey)
00300 *phi = M_PI * u1 * 0.5f;
00301 else
00302 *phi = atanf(sqrtf((ex+1)/(ey+1)) *
00303 tanf(M_PI * u1 * 0.5f));
00304 float cosphi = cosf(*phi), sinphi = sinf(*phi);
00305 *costheta = powf(u2, 1.f/(ex * cosphi * cosphi +
00306 ey * sinphi * sinphi + 1));
00307 }
00308 float Anisotropic::Pdf(const Vector &wo,
00309 const Vector &wi) const {
00310 Vector H = Normalize(wo + wi);
00311
00312 float anisotropic_pdf = D(H) / (4.f * Dot(wo, H));
00313 if (Dot(wo, H) < 0.f) anisotropic_pdf = 0.f;
00314 return anisotropic_pdf;
00315 }
00316 Spectrum FresnelBlend::Sample_f(const Vector &wo,
00317 Vector *wi, float u1, float u2, float *pdf) const {
00318 if (u1 < .5) {
00319 u1 = 2.f * u1;
00320
00321 *wi = CosineSampleHemisphere(u1, u2);
00322 if (wo.z < 0.) wi->z *= -1.f;
00323 }
00324 else {
00325 u1 = 2.f * (u1 - .5f);
00326 distribution->Sample_f(wo, wi, u1, u2, pdf);
00327 if (!SameHemisphere(wo, *wi)) return Spectrum(0.f);
00328 }
00329 *pdf = Pdf(wo, *wi);
00330 return f(wo, *wi);
00331 }
00332 float FresnelBlend::Pdf(const Vector &wo,
00333 const Vector &wi) const {
00334 if (!SameHemisphere(wo, wi)) return 0.f;
00335 return .5f * (fabsf(wi.z) * INV_PI +
00336 distribution->Pdf(wo, wi));
00337 }
00338 Spectrum BxDF::rho(const Vector &w, int nSamples,
00339 float *samples) const {
00340 if (!samples) {
00341 samples =
00342 (float *)alloca(2 * nSamples * sizeof(float));
00343 LatinHypercube(samples, nSamples, 2);
00344 }
00345 Spectrum r = 0.;
00346 for (int i = 0; i < nSamples; ++i) {
00347
00348 Vector wi;
00349 float pdf = 0.f;
00350 Spectrum f =
00351 Sample_f(w, &wi, samples[2*i], samples[2*i+1], &pdf);
00352 if (pdf > 0.) r += f * fabsf(wi.z) / pdf;
00353 }
00354 return r / nSamples;
00355 }
00356 Spectrum BxDF::rho(int nSamples, float *samples) const {
00357 if (!samples) {
00358 samples =
00359 (float *)alloca(4 * nSamples * sizeof(float));
00360 LatinHypercube(samples, nSamples, 4);
00361 }
00362 Spectrum r = 0.;
00363 for (int i = 0; i < nSamples; ++i) {
00364
00365 Vector wo, wi;
00366 wo = UniformSampleHemisphere(samples[4*i], samples[4*i+1]);
00367 float pdf_o = INV_TWOPI, pdf_i = 0.f;
00368 Spectrum f =
00369 Sample_f(wo, &wi, samples[4*i+2], samples[4*i+3],
00370 &pdf_i);
00371 if (pdf_i > 0.)
00372 r += f * fabsf(wi.z * wo.z) / (pdf_o * pdf_i);
00373 }
00374 return r / (M_PI*nSamples);
00375 }
00376
00377 Spectrum BSDF::Sample_f(const Vector &wo, Vector *wi, BxDFType flags,
00378 BxDFType *sampledType) const {
00379 float pdf;
00380 Spectrum f = Sample_f(wo, wi, RandomFloat(), RandomFloat(),
00381 RandomFloat(), &pdf, flags, sampledType);
00382 if (!f.Black() && pdf > 0.) f /= pdf;
00383 return f;
00384 }
00385 Spectrum BSDF::Sample_f(const Vector &woW, Vector *wiW,
00386 float u1, float u2, float u3, float *pdf,
00387 BxDFType flags, BxDFType *sampledType) const {
00388
00389 int matchingComps = NumComponents(flags);
00390 if (matchingComps == 0) {
00391 *pdf = 0.f;
00392 return Spectrum(0.f);
00393 }
00394 int which = min(Floor2Int(u3 * matchingComps),
00395 matchingComps-1);
00396 BxDF *bxdf = NULL;
00397 int count = which;
00398 for (int i = 0; i < nBxDFs; ++i)
00399 if (bxdfs[i]->MatchesFlags(flags))
00400 if (count-- == 0) {
00401 bxdf = bxdfs[i];
00402 break;
00403 }
00404 Assert(bxdf);
00405
00406 Vector wi;
00407 Vector wo = WorldToLocal(woW);
00408 *pdf = 0.f;
00409 Spectrum f = bxdf->Sample_f(wo, &wi, u1, u2, pdf);
00410 if (*pdf == 0.f) return 0.f;
00411 if (sampledType) *sampledType = bxdf->type;
00412 *wiW = LocalToWorld(wi);
00413
00414 if (!(bxdf->type & BSDF_SPECULAR) && matchingComps > 1) {
00415 for (int i = 0; i < nBxDFs; ++i) {
00416 if (bxdfs[i] != bxdf &&
00417 bxdfs[i]->MatchesFlags(flags))
00418 *pdf += bxdfs[i]->Pdf(wo, wi);
00419 }
00420 }
00421 if (matchingComps > 1) *pdf /= matchingComps;
00422
00423 if (!(bxdf->type & BSDF_SPECULAR)) {
00424 f = 0.;
00425 if (Dot(*wiW, ng) * Dot(woW, ng) > 0)
00426
00427 flags = BxDFType(flags & ~BSDF_TRANSMISSION);
00428 else
00429
00430 flags = BxDFType(flags & ~BSDF_REFLECTION);
00431 for (int i = 0; i < nBxDFs; ++i)
00432 if (bxdfs[i]->MatchesFlags(flags))
00433 f += bxdfs[i]->f(wo, wi);
00434 }
00435 return f;
00436 }
00437 float BSDF::Pdf(const Vector &woW, const Vector &wiW,
00438 BxDFType flags) const {
00439 if (nBxDFs == 0.) return 0.;
00440 Vector wo = WorldToLocal(woW), wi = WorldToLocal(wiW);
00441 float pdf = 0.f;
00442 int matchingComps = 0;
00443 for (int i = 0; i < nBxDFs; ++i)
00444 if (bxdfs[i]->MatchesFlags(flags)) {
00445 ++matchingComps;
00446 pdf += bxdfs[i]->Pdf(wo, wi);
00447 }
00448 return matchingComps > 0 ? pdf / matchingComps : 0.f;
00449 }
00450 BSDF::BSDF(const DifferentialGeometry &dg,
00451 const Normal &ngeom, float e)
00452 : dgShading(dg), eta(e) {
00453 ng = ngeom;
00454 nn = dgShading.nn;
00455 sn = Normalize(dgShading.dpdu);
00456 tn = Cross(nn, sn);
00457 nBxDFs = 0;
00458 }
00459 Spectrum BSDF::f(const Vector &woW,
00460 const Vector &wiW, BxDFType flags) const {
00461 Vector wi = WorldToLocal(wiW), wo = WorldToLocal(woW);
00462 if (Dot(wiW, ng) * Dot(woW, ng) > 0)
00463
00464 flags = BxDFType(flags & ~BSDF_TRANSMISSION);
00465 else
00466
00467 flags = BxDFType(flags & ~BSDF_REFLECTION);
00468 Spectrum f = 0.;
00469 for (int i = 0; i < nBxDFs; ++i)
00470 if (bxdfs[i]->MatchesFlags(flags))
00471 f += bxdfs[i]->f(wo, wi);
00472 return f;
00473 }
00474 Spectrum BSDF::rho(BxDFType flags) const {
00475 Spectrum ret(0.);
00476 for (int i = 0; i < nBxDFs; ++i)
00477 if (bxdfs[i]->MatchesFlags(flags))
00478 ret += bxdfs[i]->rho();
00479 return ret;
00480 }
00481 Spectrum BSDF::rho(const Vector &wo, BxDFType flags) const {
00482 Spectrum ret(0.);
00483 for (int i = 0; i < nBxDFs; ++i)
00484 if (bxdfs[i]->MatchesFlags(flags))
00485 ret += bxdfs[i]->rho(wo);
00486 return ret;
00487 }
00488 MemoryArena BSDF::arena;