00001
00002
00003 #include <stdio.h>
00004 #include <stdlib.h>
00005
00006 #include "pbrt.h"
00007 #include "reflection.h"
00008 #include "montecarlo.h"
00009 #include "memory.h"
00010 #include "api.h"
00011 #include "paramset.h"
00012 #include "shapes/disk.h"
00013
00014 static MemoryArena arena;
00015 static RNG rng;
00016
00017
00018 double spectrumRedValue(const Spectrum & s)
00019 {
00020 return ((float*) (& s))[0];
00021 }
00022
00023 typedef void (*CreateBSDFFunc)(BSDF* bsdf);
00024
00025 void createBlinn0(BSDF* bsdf);
00026 void createBlinn05(BSDF* bsdf);
00027 void createBlinn2(BSDF* bsdf);
00028 void createBlinn30and0(BSDF* bsdf);
00029 void createAniso0_0(BSDF* bsdf);
00030 void createAniso10_10(BSDF* bsdf);
00031 void createAniso30_30(BSDF* bsdf);
00032 void createLambertian(BSDF* bsdf);
00033 void createOrenNayar0(BSDF* bsdf);
00034 void createOrenNayar20(BSDF* bsdf);
00035 void createFresnelBlend0(BSDF* bsdf);
00036 void createFresnelBlend30(BSDF* bsdf);
00037 void createPlastic(BSDF* bsdf);
00038 void createSubstrate(BSDF* bsdf);
00039
00040
00041 typedef void (*GenSampleFunc)(BSDF* bsdf,
00042 const Vector & wo, Vector* wi, float* pdf, Spectrum* f);
00043
00044 void Gen_Sample_f(BSDF* bsdf,
00045 const Vector & wo, Vector* wi, float* pdf, Spectrum* f);
00046 void Gen_CosHemisphere(BSDF* bsdf,
00047 const Vector & wo, Vector* wi, float* pdf, Spectrum* f);
00048 void Gen_UniformHemisphere(BSDF* bsdf,
00049 const Vector & wo, Vector* wi, float* pdf, Spectrum* f);
00050
00051
00052 int main(int argc, char *argv[])
00053 {
00054 Options opt;
00055 pbrtInit(opt);
00056
00057
00058
00059 const int estimates = 100000000;
00060
00061
00062 const double environmentRadiance = 1.0;
00063
00064
00065 fprintf(stderr, "outgoing radiance from a surface viewed\n"
00066 "straight on with uniform lighting\n\n"
00067 " uniform incoming radiance = %.3f\n"
00068 " monte carlo samples = %d\n\n\n",
00069 environmentRadiance, estimates);
00070
00071
00072 CreateBSDFFunc BSDFFuncArray[] = {
00073 createBlinn0,
00074 createBlinn05,
00075 createBlinn2,
00076 createBlinn30and0,
00077 createAniso0_0,
00078 createAniso10_10,
00079 createAniso30_30,
00080
00081
00082
00083
00084
00085
00086
00087 };
00088
00089 const char* BSDFFuncDescripArray[] = {
00090 "Blinn (exponent 0)",
00091 "Blinn (exponent 0.5)",
00092 "Blinn (exponent 2)",
00093 "Blinn (exponent 30 and 0)",
00094 "Anisotropic (exponent 0, 0)",
00095 "Anisotropic (exponent 10, 10)",
00096 "Anisotropic (exponent 30, 30)",
00097
00098
00099
00100
00101
00102
00103
00104 };
00105
00106 GenSampleFunc SampleFuncArray[] = {
00107 Gen_Sample_f,
00108 Gen_CosHemisphere,
00109
00110 };
00111
00112 const char* SampleFuncDescripArray[] = {
00113 "BSDF Importance Sampling",
00114 "Cos Hemisphere",
00115
00116 };
00117
00118 int numModels = sizeof(BSDFFuncArray) / sizeof(BSDFFuncArray[0]);
00119 int numModelsDescrip = sizeof(BSDFFuncDescripArray) /
00120 sizeof(BSDFFuncDescripArray[0]);
00121 int numGenerators = sizeof(SampleFuncArray) / sizeof(SampleFuncArray[0]);
00122 int numGeneratorsDescrip = sizeof(SampleFuncDescripArray) /
00123 sizeof(SampleFuncDescripArray[0]);
00124
00125 if (numModels != numModelsDescrip) {
00126 fprintf(stderr, "BSDFFuncArray and BSDFFuncDescripArray out of sync!\n");
00127 exit(1);
00128 }
00129
00130 if (numGenerators != numGeneratorsDescrip) {
00131 fprintf(stderr, "SampleFuncArray and SampleFuncDescripArray out of sync!\n");
00132 exit(1);
00133 }
00134
00135
00136 for (int model = 0; model < numModels; model++) {
00137
00138 BSDF* bsdf;
00139
00140
00141
00142
00143 {
00144 Transform t = RotateX(-90);
00145 bool reverseOrientation = false;
00146 ParamSet p;
00147
00148 Reference<Shape> disk = new Disk(new Transform(t), new Transform(Inverse(t)),
00149 reverseOrientation, 0., 1., 0, 360.);
00150 if (!disk) {
00151 fprintf(stderr, "Could not load disk plugin\n"
00152 " make sure the PBRT_SEARCHPATH environment variable is set\n");
00153 exit(1);
00154 }
00155
00156 Point origin(0.1, 1, 0);
00157 Vector direction(0, -1, 0);
00158 float tHit, rayEps;
00159 Ray r(origin, direction, 1e-3, INFINITY);
00160 DifferentialGeometry* dg = BSDF_ALLOC(arena, DifferentialGeometry)();
00161 disk->Intersect(r, &tHit, &rayEps, dg);
00162
00163 bsdf = BSDF_ALLOC(arena, BSDF)(*dg, dg->nn);
00164 (BSDFFuncArray[model])(bsdf);
00165 }
00166
00167
00168
00169 Vector woL = Normalize(Vector(0, 0, 1));
00170 Vector wo = bsdf->LocalToWorld(woL);
00171 const Normal &n = bsdf->dgShading.nn;
00172
00173
00174 for (int gen = 0; gen < numGenerators; gen++) {
00175 double redSum = 0.0;
00176
00177 const int numHistoBins = 10;
00178 double histogram[numHistoBins][numHistoBins];
00179 for (int i = 0; i < numHistoBins; i++) {
00180 for (int j = 0; j < numHistoBins; j++) {
00181 histogram[i][j] = 0;
00182 }
00183 }
00184 int badSamples = 0;
00185 int outsideSamples = 0;
00186
00187 int warningTarget = 1;
00188 for (int sample = 0; sample < estimates; sample++) {
00189 Vector wi;
00190 float pdf;
00191 Spectrum f;
00192
00193
00194 (SampleFuncArray[gen])(bsdf, wo, & wi, & pdf, & f);
00195
00196 double redF = spectrumRedValue(f);
00197
00198
00199 Vector wiL = bsdf->WorldToLocal(wi);
00200 float x = Clamp(wiL.x, -1.f, 1.f);
00201 float y = Clamp(wiL.y, -1.f, 1.f);
00202 float wiPhi = (atan2(y, x) + M_PI) / (2.0 * M_PI);
00203 float wiCosTheta = wiL.z;
00204 bool validSample = (wiCosTheta > 1e-7);
00205 if (wiPhi < -0.0001 || wiPhi > 1.0001 || wiCosTheta > 1.0001) {
00206
00207 fprintf(stderr, "bad wi! %.3f %.3f %.3f, (%.3f %.3f)\n",
00208 wiL[0], wiL[1], wiL[2], wiPhi, wiCosTheta);
00209 } else if (validSample) {
00210 int histoPhi = (int) (wiPhi * numHistoBins);
00211 int histoCosTheta = (int) (wiCosTheta * numHistoBins);
00212 histogram[histoCosTheta][histoPhi] += 1.0 / pdf;
00213 }
00214
00215 if (!validSample) {
00216 outsideSamples++;
00217 } else if (pdf == 0.f || isnan(pdf) || redF < 0 || isnan(redF)) {
00218 if (badSamples == warningTarget) {
00219 fprintf(stderr, "warning %d, bad sample %d! "
00220 "pdf: %.3f, redF: %.3f\n",
00221 warningTarget, sample, pdf, redF);
00222 warningTarget *= 10;
00223 }
00224 badSamples++;
00225 } else {
00226
00227
00228 redSum += redF * environmentRadiance * AbsDot(wi, n) / pdf;
00229 }
00230 }
00231 int goodSamples = estimates - badSamples;
00232
00233
00234 fprintf(stderr, "*** BRDF: '%s', Samples: '%s'\n\n"
00235 "wi histogram showing the relative weight in each bin\n"
00236 " all entries should be close to 2pi = %.5f:\n"
00237 " (%d bad samples, %d outside samples)\n\n"
00238 " cos(theta) bins\n",
00239 BSDFFuncDescripArray[model], SampleFuncDescripArray[gen],
00240 M_PI * 2.0, badSamples, outsideSamples);
00241 double totalSum = 0.0;
00242 for (int i = 0; i < numHistoBins; i++) {
00243 fprintf(stderr, " phi bin %02d:", i);
00244 for (int j = 0; j < numHistoBins; j++) {
00245 fprintf(stderr, " %5.2f", histogram[i][j] *
00246 numHistoBins * numHistoBins / goodSamples);
00247 totalSum += histogram[i][j];
00248 }
00249 fprintf(stderr, "\n");
00250 }
00251 fprintf(stderr, "\n final average : %.5f (error %.5f)\n\n"
00252 " radiance = %.5f\n\n",
00253 totalSum / goodSamples, totalSum / goodSamples - M_PI * 2.0,
00254 redSum / goodSamples);
00255 }
00256 }
00257
00258 pbrtCleanup();
00259 return 0;
00260 }
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270 void Gen_Sample_f(BSDF* bsdf,
00271 const Vector & wo, Vector* wi, float* pdf, Spectrum* f)
00272 {
00273
00274 BxDFType inflags = BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE | BSDF_GLOSSY);
00275 BxDFType outflags;
00276 BSDFSample sample(rng);
00277 *f = bsdf->Sample_f(wo, wi, sample, pdf, inflags, &outflags);
00278
00279
00280 Vector wiL = bsdf->WorldToLocal(*wi);
00281 float wiCosTheta = wiL.z;
00282 bool validSample = (wiCosTheta > 1e-7);
00283
00284 if (validSample) {
00285 float verifyPdf = bsdf->Pdf(wo, *wi, inflags);
00286 if (fabs(verifyPdf - *pdf) > 1e-4) {
00287 fprintf(stderr, "BSDF::Pdf() doesn't match BSDF::Sample_f() !\n"
00288 " Sample_f pdf %.3f, Pdf pdf %.3f\n"
00289 " wo %.3f %.3f %.3f, wi %.3f %.3f %.3f\n",
00290 *pdf, verifyPdf, wo[0], wo[1], wo[2], (*wi)[0], (*wi)[1], (*wi)[2]);
00291 fprintf(stderr, "blah! validSample %d, wiCosTheta %.3f, wiL %.3f %.3f %.3f\n",
00292 validSample, wiCosTheta, wiL[0], wiL[1], wiL[2]);
00293 }
00294 }
00295 }
00296
00297 void Gen_CosHemisphere(BSDF* bsdf,
00298 const Vector & wo, Vector* wi, float* pdf, Spectrum* f)
00299 {
00300 float u1 = rng.RandomFloat();
00301 float u2 = rng.RandomFloat();
00302 Vector wiL = CosineSampleHemisphere(u1, u2);
00303 *wi = bsdf->LocalToWorld(wiL);
00304 float cosTheta = wiL.z;
00305 *pdf = CosineHemispherePdf(cosTheta, u2);
00306
00307 *f = bsdf->f(wo, *wi);
00308 }
00309
00310 void Gen_UniformHemisphere(BSDF* bsdf,
00311 const Vector & wo, Vector* wi, float* pdf, Spectrum* f)
00312 {
00313 float u1 = rng.RandomFloat();
00314 float u2 = rng.RandomFloat();
00315 Vector wiL = UniformSampleHemisphere(u1, u2);
00316 *wi = bsdf->LocalToWorld(wiL);
00317 *pdf = UniformHemispherePdf();
00318
00319 *f = bsdf->f(wo, *wi);
00320 }
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330 void createBlinn0(BSDF* bsdf)
00331 {
00332 const float blinnExponent = 0.0;
00333 Spectrum Ks(1);
00334 MicrofacetDistribution* distribution = BSDF_ALLOC(arena, Blinn)(blinnExponent);
00335 Fresnel* fresnel = BSDF_ALLOC(arena, FresnelNoOp)();
00336 BxDF* bxdf = BSDF_ALLOC(arena, Microfacet)(Ks, fresnel, distribution);
00337 bsdf->Add(bxdf);
00338 }
00339
00340
00341 void createBlinn05(BSDF* bsdf)
00342 {
00343 const float blinnExponent = 0.5;
00344 Spectrum Ks(1);
00345 MicrofacetDistribution* distribution = BSDF_ALLOC(arena, Blinn)(blinnExponent);
00346 Fresnel* fresnel = BSDF_ALLOC(arena, FresnelNoOp)();
00347 BxDF* bxdf = BSDF_ALLOC(arena, Microfacet)(Ks, fresnel, distribution);
00348 bsdf->Add(bxdf);
00349 }
00350
00351
00352 void createBlinn2(BSDF* bsdf)
00353 {
00354 const float blinnExponent = 2.0;
00355 Spectrum Ks(1);
00356 MicrofacetDistribution* distribution = BSDF_ALLOC(arena, Blinn)(blinnExponent);
00357 Fresnel* fresnel = BSDF_ALLOC(arena, FresnelNoOp)();
00358 BxDF* bxdf = BSDF_ALLOC(arena, Microfacet)(Ks, fresnel, distribution);
00359 bsdf->Add(bxdf);
00360 }
00361
00362
00363 void createBlinn30and0(BSDF* bsdf)
00364 {
00365 const float blinnExponent1 = 30.0;
00366 Spectrum Ks(0.5);
00367 MicrofacetDistribution* distribution1 = BSDF_ALLOC(arena, Blinn)(blinnExponent1);
00368 Fresnel* fresnel = BSDF_ALLOC(arena, FresnelNoOp)();
00369 BxDF* bxdf1 = BSDF_ALLOC(arena, Microfacet)(Ks, fresnel, distribution1);
00370 bsdf->Add(bxdf1);
00371
00372 const float blinnExponent2 = 0.0;
00373 MicrofacetDistribution* distribution2 = BSDF_ALLOC(arena, Blinn)(blinnExponent2);
00374 BxDF* bxdf2 = BSDF_ALLOC(arena, Microfacet)(Ks, fresnel, distribution2);
00375 bsdf->Add(bxdf2);
00376 }
00377
00378
00379 void createAniso0_0(BSDF* bsdf)
00380 {
00381 const float aniso1 = 0.0;
00382 const float aniso2 = 0.0;
00383 Spectrum Ks(1);
00384 MicrofacetDistribution* distribution =
00385 BSDF_ALLOC(arena, Anisotropic(aniso1, aniso2));
00386 Fresnel* fresnel = BSDF_ALLOC(arena, FresnelNoOp)();
00387 BxDF* bxdf = BSDF_ALLOC(arena, Microfacet)(Ks, fresnel, distribution);
00388 bsdf->Add(bxdf);
00389 }
00390
00391
00392 void createAniso10_10(BSDF* bsdf)
00393 {
00394 const float aniso1 = 10.0;
00395 const float aniso2 = 10.0;
00396 Spectrum Ks(1);
00397 MicrofacetDistribution* distribution =
00398 BSDF_ALLOC(arena, Anisotropic(aniso1, aniso2));
00399 Fresnel* fresnel = BSDF_ALLOC(arena, FresnelNoOp)();
00400 BxDF* bxdf = BSDF_ALLOC(arena, Microfacet)(Ks, fresnel, distribution);
00401 bsdf->Add(bxdf);
00402 }
00403
00404 void createAniso30_30(BSDF* bsdf)
00405 {
00406 const float aniso1 = 30.0;
00407 const float aniso2 = 30.0;
00408 Spectrum Ks(1);
00409 MicrofacetDistribution* distribution =
00410 BSDF_ALLOC(arena, Anisotropic(aniso1, aniso2));
00411 Fresnel* fresnel = BSDF_ALLOC(arena, FresnelNoOp)();
00412 BxDF* bxdf = BSDF_ALLOC(arena, Microfacet)(Ks, fresnel, distribution);
00413 bsdf->Add(bxdf);
00414 }
00415
00416
00417 void createLambertian(BSDF* bsdf)
00418 {
00419 Spectrum Kd(1);
00420 BxDF* bxdf = BSDF_ALLOC(arena, Lambertian)(Kd);
00421 bsdf->Add(bxdf);
00422 }
00423
00424
00425 void createOrenNayar0(BSDF* bsdf)
00426 {
00427 Spectrum Kd(1);
00428 float sigma = 20.0;
00429 BxDF* bxdf = BSDF_ALLOC(arena, OrenNayar)(Kd, sigma);
00430 bsdf->Add(bxdf);
00431 }
00432
00433
00434 void createOrenNayar20(BSDF* bsdf)
00435 {
00436 Spectrum Kd(1);
00437 float sigma = 0.0;
00438 BxDF* bxdf = BSDF_ALLOC(arena, OrenNayar)(Kd, sigma);
00439 bsdf->Add(bxdf);
00440 }
00441
00442
00443 void createFresnelBlend0(BSDF* bsdf)
00444 {
00445 Spectrum d(0.5);
00446 Spectrum s(0.5);
00447 float exponent = 0.0;
00448
00449 MicrofacetDistribution* distribution = BSDF_ALLOC(arena, Blinn)(exponent);
00450 BxDF* bxdf = BSDF_ALLOC(arena, FresnelBlend)(d, s, distribution);
00451 bsdf->Add(bxdf);
00452 }
00453
00454
00455 void createFresnelBlend30(BSDF* bsdf)
00456 {
00457 Spectrum d(0.5);
00458 Spectrum s(0.5);
00459 float exponent = 30.0;
00460
00461 MicrofacetDistribution* distribution = BSDF_ALLOC(arena, Blinn)(exponent);
00462 BxDF* bxdf = BSDF_ALLOC(arena, FresnelBlend)(d, s, distribution);
00463 bsdf->Add(bxdf);
00464 }
00465
00466
00467 void createPlastic(BSDF* bsdf)
00468 {
00469
00470
00471 Spectrum kd(0.5);
00472 BxDF *diff = BSDF_ALLOC(arena, Lambertian)(kd);
00473 Fresnel *fresnel = BSDF_ALLOC(arena, FresnelDielectric)(1.5f, 1.f);
00474 Spectrum ks = (0.5);
00475 float rough = 0.1;
00476 BxDF *spec = BSDF_ALLOC(arena, Microfacet)(ks, fresnel,
00477 BSDF_ALLOC(arena, Blinn)(1.f / rough));
00478 bsdf->Add(diff);
00479 bsdf->Add(spec);
00480 }
00481
00482
00483
00484 void createSubstrate(BSDF* bsdf)
00485 {
00486
00487
00488 Spectrum d(0.5);
00489 Spectrum s(0.5);
00490 float u = 0.1;
00491 float v = 0.1;
00492
00493 bsdf->Add(BSDF_ALLOC(arena, FresnelBlend)(d, s,
00494 BSDF_ALLOC(arena, Anisotropic)(1.f/u, 1.f/v)));
00495 }
00496
00497
00498