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
00026 #include "stdafx.h"
00027 #include "integrators/dipolesubsurface.h"
00028 #include "scene.h"
00029 #include "montecarlo.h"
00030 #include "sampler.h"
00031 #include "progressreporter.h"
00032 #include "intersection.h"
00033 #include "paramset.h"
00034 #include "reflection.h"
00035 #include "octree.h"
00036 #include "camera.h"
00037 #include "floatfile.h"
00038 struct DiffusionReflectance;
00039
00040
00041 struct SubsurfaceOctreeNode {
00042
00043 SubsurfaceOctreeNode() {
00044 isLeaf = true;
00045 sumArea = 0.f;
00046 for (int i = 0; i < 8; ++i)
00047 ips[i] = NULL;
00048 }
00049 void Insert(const BBox &nodeBound, IrradiancePoint *ip,
00050 MemoryArena &arena) {
00051 Point pMid = .5f * nodeBound.pMin + .5f * nodeBound.pMax;
00052 if (isLeaf) {
00053
00054 for (int i = 0; i < 8; ++i) {
00055 if (!ips[i]) {
00056 ips[i] = ip;
00057 return;
00058 }
00059 }
00060
00061
00062 isLeaf = false;
00063 IrradiancePoint *localIps[8];
00064 for (int i = 0; i < 8; ++i) {
00065 localIps[i] = ips[i];
00066 children[i] = NULL;
00067 }
00068 for (int i = 0; i < 8; ++i) {
00069 IrradiancePoint *ip = localIps[i];
00070
00071 int child = (ip->p.x > pMid.x ? 4 : 0) +
00072 (ip->p.y > pMid.y ? 2 : 0) + (ip->p.z > pMid.z ? 1 : 0);
00073 if (!children[child])
00074 children[child] = arena.Alloc<SubsurfaceOctreeNode>();
00075 BBox childBound = octreeChildBound(child, nodeBound, pMid);
00076 children[child]->Insert(childBound, ip, arena);
00077 }
00078
00079 }
00080
00081 int child = (ip->p.x > pMid.x ? 4 : 0) +
00082 (ip->p.y > pMid.y ? 2 : 0) + (ip->p.z > pMid.z ? 1 : 0);
00083 if (!children[child])
00084 children[child] = arena.Alloc<SubsurfaceOctreeNode>();
00085 BBox childBound = octreeChildBound(child, nodeBound, pMid);
00086 children[child]->Insert(childBound, ip, arena);
00087 }
00088 void InitHierarchy() {
00089 if (isLeaf) {
00090
00091 float sumWt = 0.f;
00092 uint32_t i;
00093 for (i = 0; i < 8; ++i) {
00094 if (!ips[i]) break;
00095 float wt = ips[i]->E.y();
00096 E += ips[i]->E;
00097 p += wt * ips[i]->p;
00098 sumWt += wt;
00099 sumArea += ips[i]->area;
00100 }
00101 if (sumWt > 0.f) p /= sumWt;
00102 E /= i;
00103 }
00104 else {
00105
00106 float sumWt = 0.f;
00107 uint32_t nChildren = 0;
00108 for (uint32_t i = 0; i < 8; ++i) {
00109 if (!children[i]) continue;
00110 ++nChildren;
00111 children[i]->InitHierarchy();
00112 float wt = children[i]->E.y();
00113 E += children[i]->E;
00114 p += wt * children[i]->p;
00115 sumWt += wt;
00116 sumArea += children[i]->sumArea;
00117 }
00118 if (sumWt > 0.f) p /= sumWt;
00119 E /= nChildren;
00120 }
00121 }
00122 Spectrum Mo(const BBox &nodeBound, const Point &p, const DiffusionReflectance &Rd,
00123 float maxError);
00124
00125
00126 Point p;
00127 bool isLeaf;
00128 Spectrum E;
00129 float sumArea;
00130 union {
00131 SubsurfaceOctreeNode *children[8];
00132 IrradiancePoint *ips[8];
00133 };
00134 };
00135
00136
00137 struct DiffusionReflectance {
00138
00139 DiffusionReflectance(const Spectrum &sigma_a, const Spectrum &sigmap_s,
00140 float eta) {
00141 A = (1.f + Fdr(eta)) / (1.f - Fdr(eta));
00142 sigmap_t = sigma_a + sigmap_s;
00143 sigma_tr = Sqrt(3.f * sigma_a * sigmap_t);
00144 alphap = sigmap_s / sigmap_t;
00145 zpos = Spectrum(1.f) / sigmap_t;
00146 zneg = zpos * (1.f + (4.f/3.f) * A);
00147 }
00148 Spectrum operator()(float d2) const {
00149 Spectrum dpos = Sqrt(Spectrum(d2) + zpos * zpos);
00150 Spectrum dneg = Sqrt(Spectrum(d2) + zneg * zneg);
00151 Spectrum Rd = (1.f / (4.f * M_PI)) *
00152 ((zpos * (dpos * sigma_tr + Spectrum(1.f)) *
00153 Exp(-sigma_tr * dpos)) / (dpos * dpos * dpos) -
00154 (zneg * (dneg * sigma_tr + Spectrum(1.f)) *
00155 Exp(-sigma_tr * dneg)) / (dneg * dneg * dneg));
00156 return Rd.Clamp();
00157 }
00158
00159
00160 Spectrum zpos, zneg, sigmap_t, sigma_tr, alphap;
00161 float A;
00162 };
00163
00164
00165
00166
00167 DipoleSubsurfaceIntegrator::~DipoleSubsurfaceIntegrator() {
00168 delete[] lightSampleOffsets;
00169 delete[] bsdfSampleOffsets;
00170 }
00171
00172
00173 void DipoleSubsurfaceIntegrator::RequestSamples(Sampler *sampler, Sample *sample,
00174 const Scene *scene) {
00175
00176 uint32_t nLights = scene->lights.size();
00177 lightSampleOffsets = new LightSampleOffsets[nLights];
00178 bsdfSampleOffsets = new BSDFSampleOffsets[nLights];
00179 for (uint32_t i = 0; i < nLights; ++i) {
00180 const Light *light = scene->lights[i];
00181 int nSamples = light->nSamples;
00182 if (sampler) nSamples = sampler->RoundSize(nSamples);
00183 lightSampleOffsets[i] = LightSampleOffsets(nSamples, sample);
00184 bsdfSampleOffsets[i] = BSDFSampleOffsets(nSamples, sample);
00185 }
00186 }
00187
00188
00189 void DipoleSubsurfaceIntegrator::Preprocess(const Scene *scene,
00190 const Camera *camera, const Renderer *renderer) {
00191 if (scene->lights.size() == 0) return;
00192 vector<SurfacePoint> pts;
00193
00194 if (filename != "") {
00195
00196 vector<float> fpts;
00197 if (ReadFloatFile(filename.c_str(), &fpts)) {
00198 if ((fpts.size() % 8) != 0)
00199 Error("Excess values (%d) in points file \"%s\"", int(fpts.size() % 8),
00200 filename.c_str());
00201 for (u_int i = 0; i < fpts.size(); i += 8)
00202 pts.push_back(SurfacePoint(Point(fpts[i], fpts[i+1], fpts[i+2]),
00203 Normal(fpts[i+3], fpts[i+4], fpts[i+5]),
00204 fpts[i+6], fpts[i+7]));
00205 }
00206 }
00207 if (pts.size() == 0) {
00208 Point pCamera = camera->CameraToWorld(camera->shutterOpen,
00209 Point(0, 0, 0));
00210 FindPoissonPointDistribution(pCamera, camera->shutterOpen,
00211 minSampleDist, scene, &pts);
00212 }
00213
00214
00215 RNG rng;
00216 MemoryArena arena;
00217 PBRT_SUBSURFACE_STARTED_COMPUTING_IRRADIANCE_VALUES();
00218 ProgressReporter progress(pts.size(), "Computing Irradiances");
00219 for (uint32_t i = 0; i < pts.size(); ++i) {
00220 SurfacePoint &sp = pts[i];
00221 Spectrum E(0.f);
00222 for (uint32_t j = 0; j < scene->lights.size(); ++j) {
00223
00224 const Light *light = scene->lights[j];
00225 Spectrum Elight = 0.f;
00226 int nSamples = RoundUpPow2(light->nSamples);
00227 uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
00228 uint32_t compScramble = rng.RandomUInt();
00229 for (int s = 0; s < nSamples; ++s) {
00230 float lpos[2];
00231 Sample02(s, scramble, lpos);
00232 float lcomp = VanDerCorput(s, compScramble);
00233 LightSample ls(lpos[0], lpos[1], lcomp);
00234 Vector wi;
00235 float lightPdf;
00236 VisibilityTester visibility;
00237 Spectrum Li = light->Sample_L(sp.p, sp.rayEpsilon,
00238 ls, camera->shutterOpen, &wi, &lightPdf, &visibility);
00239 if (Dot(wi, sp.n) <= 0.) continue;
00240 if (Li.IsBlack() || lightPdf == 0.f) continue;
00241 Li *= visibility.Transmittance(scene, renderer, NULL, rng, arena);
00242 if (visibility.Unoccluded(scene))
00243 Elight += Li * AbsDot(wi, sp.n) / lightPdf;
00244 }
00245 E += Elight / nSamples;
00246 }
00247 irradiancePoints.push_back(IrradiancePoint(sp, E));
00248 PBRT_SUBSURFACE_COMPUTED_IRRADIANCE_AT_POINT(&sp, &E);
00249 arena.FreeAll();
00250 progress.Update();
00251 }
00252 progress.Done();
00253 PBRT_SUBSURFACE_FINISHED_COMPUTING_IRRADIANCE_VALUES();
00254
00255
00256 octree = octreeArena.Alloc<SubsurfaceOctreeNode>();
00257 for (uint32_t i = 0; i < irradiancePoints.size(); ++i)
00258 octreeBounds = Union(octreeBounds, irradiancePoints[i].p);
00259 for (uint32_t i = 0; i < irradiancePoints.size(); ++i)
00260 octree->Insert(octreeBounds, &irradiancePoints[i], octreeArena);
00261 octree->InitHierarchy();
00262 }
00263
00264
00265 Spectrum DipoleSubsurfaceIntegrator::Li(const Scene *scene, const Renderer *renderer,
00266 const RayDifferential &ray, const Intersection &isect,
00267 const Sample *sample, RNG &rng, MemoryArena &arena) const {
00268 Spectrum L(0.);
00269 Vector wo = -ray.d;
00270
00271 L += isect.Le(wo);
00272
00273
00274 BSDF *bsdf = isect.GetBSDF(ray, arena);
00275 const Point &p = bsdf->dgShading.p;
00276 const Normal &n = bsdf->dgShading.nn;
00277
00278 BSSRDF *bssrdf = isect.GetBSSRDF(ray, arena);
00279 if (bssrdf && octree) {
00280 Spectrum sigma_a = bssrdf->sigma_a();
00281 Spectrum sigmap_s = bssrdf->sigma_prime_s();
00282 Spectrum sigmap_t = sigmap_s + sigma_a;
00283 if (!sigmap_t.IsBlack()) {
00284
00285 PBRT_SUBSURFACE_STARTED_OCTREE_LOOKUP(const_cast<Point *>(&p));
00286 DiffusionReflectance Rd(sigma_a, sigmap_s, bssrdf->eta());
00287 Spectrum Mo = octree->Mo(octreeBounds, p, Rd, maxError);
00288 FresnelDielectric fresnel(1.f, bssrdf->eta());
00289 Spectrum Ft = Spectrum(1.f) - fresnel.Evaluate(AbsDot(wo, n));
00290 float Fdt = 1.f - Fdr(bssrdf->eta());
00291 L += (INV_PI * Ft) * (Fdt * Mo);
00292 PBRT_SUBSURFACE_FINISHED_OCTREE_LOOKUP();
00293 }
00294 }
00295 L += UniformSampleAllLights(scene, renderer, arena, p, n,
00296 wo, isect.rayEpsilon, ray.time, bsdf, sample, rng, lightSampleOffsets,
00297 bsdfSampleOffsets);
00298 if (ray.depth < maxSpecularDepth) {
00299
00300 L += SpecularReflect(ray, bsdf, rng, isect, renderer, scene, sample,
00301 arena);
00302 L += SpecularTransmit(ray, bsdf, rng, isect, renderer, scene, sample,
00303 arena);
00304 }
00305 return L;
00306 }
00307
00308
00309 Spectrum SubsurfaceOctreeNode::Mo(const BBox &nodeBound, const Point &pt,
00310 const DiffusionReflectance &Rd, float maxError) {
00311
00312 float dw = sumArea / DistanceSquared(pt, p);
00313 if (dw < maxError && !nodeBound.Inside(pt))
00314 {
00315 PBRT_SUBSURFACE_ADDED_INTERIOR_CONTRIBUTION(const_cast<SubsurfaceOctreeNode *>(this));
00316 return Rd(DistanceSquared(pt, p)) * E * sumArea;
00317 }
00318
00319
00320 Spectrum Mo = 0.f;
00321 if (isLeaf) {
00322
00323 for (int i = 0; i < 8; ++i) {
00324 if (!ips[i]) break;
00325 PBRT_SUBSURFACE_ADDED_POINT_CONTRIBUTION(const_cast<IrradiancePoint *>(ips[i]));
00326 Mo += Rd(DistanceSquared(pt, ips[i]->p)) * ips[i]->E * ips[i]->area;
00327 }
00328 }
00329 else {
00330
00331 Point pMid = .5f * nodeBound.pMin + .5f * nodeBound.pMax;
00332 for (int child = 0; child < 8; ++child) {
00333 if (!children[child]) continue;
00334 BBox childBound = octreeChildBound(child, nodeBound, pMid);
00335 Mo += children[child]->Mo(childBound, pt, Rd, maxError);
00336 }
00337 }
00338 return Mo;
00339 }
00340
00341
00342 DipoleSubsurfaceIntegrator *CreateDipoleSubsurfaceIntegrator(const ParamSet ¶ms) {
00343 int maxDepth = params.FindOneInt("maxdepth", 5);
00344 float maxError = params.FindOneFloat("maxerror", .05f);
00345 float minDist = params.FindOneFloat("minsampledistance", .25f);
00346 string pointsfile = params.FindOneString("pointsfile", "");
00347 if (PbrtOptions.quickRender) { maxError *= 4.f; minDist *= 4.f; }
00348 return new DipoleSubsurfaceIntegrator(maxDepth, maxError, minDist, pointsfile);
00349 }
00350
00351