C++ BSDF::Sample_f方法代码示例

static Spectrum L(const Scene *scene, const Renderer *renderer,
        const Camera *camera, MemoryArena &arena, RNG &rng, int maxDepth,
        bool ignoreDirect, const MLTSample &sample) {
    // Generate camera ray from Metropolis sample
    RayDifferential ray;
    float cameraWeight = camera->GenerateRayDifferential(sample.cameraSample,
                                                         &ray);
    Spectrum pathThroughput = cameraWeight, L = 0.;
    bool specularBounce = false, allSpecular = true;
    for (int pathLength = 0; pathLength < maxDepth; ++pathLength) {
        // Find next intersection in Metropolis light path
        Intersection isect;
        if (!scene->Intersect(ray, &isect)) {
            bool includeLe = ignoreDirect ? (specularBounce && !allSpecular) :
                                            (pathLength == 0 || specularBounce);
            if (includeLe)
                for (uint32_t i = 0; i < scene->lights.size(); ++i)
                   L += pathThroughput * scene->lights[i]->Le(ray);
            break;
        }
        if (ignoreDirect ? (specularBounce && !allSpecular) :
                           (specularBounce || pathLength == 0))
            L += pathThroughput * isect.Le(-ray.d);
        BSDF *bsdf = isect.GetBSDF(ray, arena);
        const Point &p = bsdf->dgShading.p;
        const Normal &n = bsdf->dgShading.nn;
        Vector wo = -ray.d;
        const PathSample &ps = sample.pathSamples[pathLength];
        // Sample direct illumination for Metropolis path vertex
        if (!ignoreDirect || pathLength > 0) {
            LightSample lightSample(ps.lightDir0, ps.lightDir1, ps.lightNum0);
            BSDFSample bsdfSample(ps.bsdfLightDir0, ps.bsdfLightDir1,
                                  ps.bsdfLightComponent);
            uint32_t lightNum = Floor2Int(ps.lightNum1 * scene->lights.size());
            lightNum = min(lightNum, (uint32_t)(scene->lights.size()-1));
            const Light *light = scene->lights[lightNum];
            L += pathThroughput *
                 EstimateDirect(scene, renderer, arena, light, p, n, wo,
                     isect.rayEpsilon, sample.cameraSample.time, bsdf, rng,
                     lightSample, bsdfSample);
        }

        // Sample direction for outgoing Metropolis path direction
        BSDFSample outgoingBSDFSample(ps.bsdfDir0, ps.bsdfDir1,
                                      ps.bsdfComponent);
        Vector wi;
        float pdf;
        BxDFType flags;
        Spectrum f = bsdf->Sample_f(wo, &wi, outgoingBSDFSample,
                                    &pdf, BSDF_ALL, &flags);
        if (f.IsBlack() || pdf == 0.)
            break;
        specularBounce = (flags & BSDF_SPECULAR) != 0;
        allSpecular &= specularBounce;
        pathThroughput *= f * AbsDot(wi, n) / pdf;
        ray = RayDifferential(p, wi, ray, isect.rayEpsilon);
        //pathThroughput *= renderer->Transmittance(scene, ray, NULL, rng, arena);
    }
    return L;
}

// Metropolis Method Definitions
static uint32_t GeneratePath(const RayDifferential &r,
        const Spectrum &a, const Scene *scene, MemoryArena &arena,
        const vector<PathSample> &samples, PathVertex *path,
        RayDifferential *escapedRay, Spectrum *escapedAlpha) {
    PBRT_MLT_STARTED_GENERATE_PATH();
    RayDifferential ray = r;
    Spectrum alpha = a;
    if (escapedAlpha) *escapedAlpha = 0.f;
    uint32_t length = 0;
    for (; length < samples.size(); ++length) {
        // Try to generate next vertex of ray path
        PathVertex &v = path[length];
        if (!scene->Intersect(ray, &v.isect)) {
            // Handle ray that leaves the scene during path generation
            if (escapedAlpha) *escapedAlpha = alpha;
            if (escapedRay)   *escapedRay = ray;
            break;
        }

        // Record information for current path vertex
        v.alpha = alpha;
        BSDF *bsdf = v.isect.GetBSDF(ray, arena);
        v.bsdf = bsdf;
        v.wPrev = -ray.d;

        // Sample direction for outgoing Metropolis path direction
        float pdf;
        BxDFType flags;
        Spectrum f = bsdf->Sample_f(-ray.d, &v.wNext, samples[length].bsdfSample,
                                    &pdf, BSDF_ALL, &flags);
        v.specularBounce = (flags & BSDF_SPECULAR) != 0;
        v.nSpecularComponents = bsdf->NumComponents(BxDFType(BSDF_SPECULAR |
                                         BSDF_REFLECTION | BSDF_TRANSMISSION));
        if (f.IsBlack() || pdf == 0.f)
        {
            PBRT_MLT_FINISHED_GENERATE_PATH();
            return length+1;
        }

        // Terminate path with RR or prepare for finding next vertex
        const Point &p = bsdf->dgShading.p;
        const Normal &n = bsdf->dgShading.nn;
        Spectrum pathScale = f * AbsDot(v.wNext, n) / pdf;
        float rrSurviveProb = min(1.f, pathScale.y());
        if (samples[length].rrSample > rrSurviveProb)
        {
            PBRT_MLT_FINISHED_GENERATE_PATH();
            return length+1;
        }
        alpha *= pathScale / rrSurviveProb;
        //alpha *= renderer->Transmittance(scene, ray, NULL, rng, arena);
        ray = RayDifferential(p, v.wNext, ray, v.isect.rayEpsilon);
    }
    PBRT_MLT_FINISHED_GENERATE_PATH();
    return length;
}

Spectrum IrradianceCacheIntegrator::pathL(Ray &r, const Scene *scene,
        const Renderer *renderer, const Sample *sample, MemoryArena &arena) const {
    Spectrum L(0.f);
    Spectrum pathThroughput = 1.;
    RayDifferential ray(r);
    bool specularBounce = false;
    for (int pathLength = 0; ; ++pathLength) {
        // Find next vertex of path
        Intersection isect;
        if (!scene->Intersect(ray, &isect))
            break;
        if (pathLength == 0)
            r.maxt = ray.maxt;
        else if (pathLength == 1)
            pathThroughput *= renderer->Transmittance(scene, ray, sample, arena, NULL);
        else
            pathThroughput *= renderer->Transmittance(scene, ray, NULL, arena, sample->rng);
        // Possibly add emitted light at path vertex
        if (specularBounce)
            L += pathThroughput * isect.Le(-ray.d);
        // Evaluate BSDF at hit point
        BSDF *bsdf = isect.GetBSDF(ray, arena);
        // Sample illumination from lights to find path contribution
        const Point &p = bsdf->dgShading.p;
        const Normal &n = bsdf->dgShading.nn;
        Vector wo = -ray.d;
        L += pathThroughput *
            UniformSampleOneLight(scene, renderer, arena, p, n, wo, isect.rayEpsilon,
                                  bsdf, sample);
        if (pathLength+1 == maxIndirectDepth) break;
        // Sample BSDF to get new path direction
        // Get random numbers for sampling new direction, \mono{bs1}, \mono{bs2}, and \mono{bcs}
        Vector wi;
        float pdf;
        BxDFType flags;
        Spectrum f = bsdf->Sample_f(wo, &wi, BSDFSample(*sample->rng),
            &pdf, BSDF_ALL, &flags);
        if (f.IsBlack() || pdf == 0.)
            break;
        specularBounce = (flags & BSDF_SPECULAR) != 0;
        pathThroughput *= f * AbsDot(wi, n) / pdf;
        ray = RayDifferential(p, wi, ray, isect.rayEpsilon);
        // Possibly terminate the path
        if (pathLength > 2) {
            float rrProb = min(1.f, pathThroughput.y());
            if (sample->rng->RandomFloat() > rrProb)
                break;
            pathThroughput /= rrProb;
        }
    }
    return L;
}

Spectrum PathIntegrator::Li(const Scene *scene, const Renderer *renderer,
        const RayDifferential &r, const Intersection &isect,
        const Sample *sample, MemoryArena &arena) const {
    // Declare common path integration variables
    Spectrum pathThroughput = 1., L = 0.;
    RayDifferential ray(r);
    bool specularBounce = false;
    Intersection localIsect;
    const Intersection *isectp = &isect;
    for (int pathLength = 0; ; ++pathLength) {
        // Possibly add emitted light at path vertex
        if (pathLength == 0 || specularBounce)
            L += pathThroughput * isectp->Le(-ray.d);

        // Sample illumination from lights to find path contribution
        BSDF *bsdf = isectp->GetBSDF(ray, arena);
        const Point &p = bsdf->dgShading.p;
        const Normal &n = bsdf->dgShading.nn;
        Vector wo = -ray.d;
        if (pathLength < SAMPLE_DEPTH)
            L += pathThroughput *
                 UniformSampleOneLight(scene, renderer, arena, p, n, wo,
                     isectp->RayEpsilon, bsdf, sample, lightNumOffset[pathLength],
                     &lightSampleOffsets[pathLength], &bsdfSampleOffsets[pathLength]);
        else
            L += pathThroughput *
                 UniformSampleOneLight(scene, renderer, arena, p, n, wo,
                     isectp->RayEpsilon, bsdf, sample);

        // Sample BSDF to get new path direction

        // Get _outgoingBSDFSample_ for sampling new path direction
        BSDFSample outgoingBSDFSample;
        if (pathLength < SAMPLE_DEPTH)
            outgoingBSDFSample = BSDFSample(sample, pathSampleOffsets[pathLength], 0);
        else
            outgoingBSDFSample = BSDFSample(*sample->rng);
        Vector wi;
        float pdf;
        BxDFType flags;
        Spectrum f = bsdf->Sample_f(wo, &wi, outgoingBSDFSample,
                                    &pdf, BSDF_ALL, &flags);
        if (f.IsBlack() || pdf == 0.)
            break;
        specularBounce = (flags & BSDF_SPECULAR) != 0;
        pathThroughput *= f * AbsDot(wi, n) / pdf;
        ray = RayDifferential(p, wi, ray, isectp->RayEpsilon);

        // Possibly terminate the path
        if (pathLength > 3) {
            float continueProbability = min(.5f, pathThroughput.y());
            if (sample->rng->RandomFloat() > continueProbability)
                break;
            pathThroughput /= continueProbability;
        }
        if (pathLength == maxDepth)
            break;

        // Find next vertex of path
        if (!scene->Intersect(ray, &localIsect)) {
            if (specularBounce) {
                for (u_int i = 0; i < scene->lights.size(); ++i)
                   L += pathThroughput * scene->lights[i]->Le(ray);
            }
            break;
        }
        if (pathLength > 1)
            pathThroughput *= renderer->Transmittance(scene, ray, NULL, arena, sample->rng);
        isectp = &localIsect;
    }
    return L;
}

// WhittedIntegrator Method Definitions
Spectrum WhittedIntegrator::Li(const Scene *scene,
		const RayDifferential &ray, const Sample *sample,
		float *alpha) const {
	Intersection isect;
	Spectrum L(0.);
	bool hitSomething;
	// Search for ray-primitive intersection
	hitSomething = scene->Intersect(ray, &isect);
	if (!hitSomething) {
		// Handle ray with no intersection
		if (alpha) *alpha = 0.;
		for (u_int i = 0; i < scene->lights.size(); ++i)
			L += scene->lights[i]->Le(ray);
		if (alpha && !L.Black()) *alpha = 1.;
		return L;
	}
	else {
		// Initialize _alpha_ for ray hit
		if (alpha) *alpha = 1.;
		// Compute emitted and reflected light at ray intersection point
		// Evaluate BSDF at hit point
		BSDF *bsdf = isect.GetBSDF(ray);
		// Initialize common variables for Whitted integrator
		const Point &p = bsdf->dgShading.p;
		const Normal &n = bsdf->dgShading.nn;
		Vector wo = -ray.d;
		// Compute emitted light if ray hit an area light source
		L += isect.Le(wo);
		// Add contribution of each light source
		Vector wi;
		for (u_int i = 0; i < scene->lights.size(); ++i) {
			VisibilityTester visibility;
			Spectrum Li = scene->lights[i]->Sample_L(p, &wi, &visibility);
			if (Li.Black()) continue;
			Spectrum f = bsdf->f(wo, wi);
			if (!f.Black() && visibility.Unoccluded(scene))
				L += f * Li * AbsDot(wi, n) * visibility.Transmittance(scene);
		}
		if (rayDepth++ < maxDepth) {
			// Trace rays for specular reflection and refraction
			Spectrum f = bsdf->Sample_f(wo, &wi,
				BxDFType(BSDF_REFLECTION | BSDF_SPECULAR));
			if (!f.Black() && AbsDot(wi, n) > 0.f) {
				// Compute ray differential _rd_ for specular reflection
				RayDifferential rd(p, wi);
				rd.hasDifferentials = true;
				rd.rx.o = p + isect.dg.dpdx;
				rd.ry.o = p + isect.dg.dpdy;
				// Compute differential reflected directions
				Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx +
					bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
				Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy +
					bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
				Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
				float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
				float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
				rd.rx.d = wi -
				          dwodx + 2 * Vector(Dot(wo, n) * dndx +
						  dDNdx * n);
				rd.ry.d = wi -
				          dwody + 2 * Vector(Dot(wo, n) * dndy +
						  dDNdy * n);
				L += scene->Li(rd, sample) * f * AbsDot(wi, n);
			}
			f = bsdf->Sample_f(wo, &wi,
				BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));
			if (!f.Black() && AbsDot(wi, n) > 0.f) {
				// Compute ray differential _rd_ for specular transmission
				RayDifferential rd(p, wi);
				rd.hasDifferentials = true;
				rd.rx.o = p + isect.dg.dpdx;
				rd.ry.o = p + isect.dg.dpdy;
				
				float eta = bsdf->eta;
				Vector w = -wo;
				if (Dot(wo, n) < 0) eta = 1.f / eta;
				
				Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx + bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
				Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy + bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
				
				Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
				float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
				float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
				
				float mu = eta * Dot(w, n) - Dot(wi, n);
				float dmudx = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdx;
				float dmudy = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdy;
				
				rd.rx.d = wi + eta * dwodx - Vector(mu * dndx + dmudx * n);
				rd.ry.d = wi + eta * dwody - Vector(mu * dndy + dmudy * n);
				L += scene->Li(rd, sample) * f * AbsDot(wi, n);
			}
		}
		--rayDepth;
	}
	return L;
}

Spectrum IGIIntegrator::Li(const Scene *scene, const Renderer *renderer,
        const RayDifferential &ray, const Intersection &isect,
        const Sample *sample, RNG &rng, MemoryArena &arena) const {
    Spectrum L(0.);
    Vector wo = -ray.d;
    // Compute emitted light if ray hit an area light source
    L += isect.Le(wo);

    // Evaluate BSDF at hit point
    BSDF *bsdf = isect.GetBSDF(ray, arena);
    const Point &p = bsdf->dgShading.p;
    const Normal &n = bsdf->dgShading.nn;
    L += UniformSampleAllLights(scene, renderer, arena, p, n,
                    wo, isect.rayEpsilon, ray.time, bsdf, sample, rng,
                    lightSampleOffsets, bsdfSampleOffsets);
    // Compute indirect illumination with virtual lights
    uint32_t lSet = min(uint32_t(sample->oneD[vlSetOffset][0] * nLightSets),
                        nLightSets-1);
    for (uint32_t i = 0; i < virtualLights[lSet].size(); ++i) {
        const VirtualLight &vl = virtualLights[lSet][i];
        // Compute virtual light's tentative contribution _Llight_
        float d2 = DistanceSquared(p, vl.p);
        Vector wi = Normalize(vl.p - p);
        float G = AbsDot(wi, n) * AbsDot(wi, vl.n) / d2;
        G = min(G, gLimit);
        Spectrum f = bsdf->f(wo, wi);
        if (G == 0.f || f.IsBlack()) continue;
        Spectrum Llight = f * G * vl.pathContrib / nLightPaths;
        RayDifferential connectRay(p, wi, ray, isect.rayEpsilon,
                                   sqrtf(d2) * (1.f - vl.rayEpsilon));
        Llight *= renderer->Transmittance(scene, connectRay, NULL, rng, arena);

        // Possibly skip virtual light shadow ray with Russian roulette
        if (Llight.y() < rrThreshold) {
            float continueProbability = .1f;
            if (rng.RandomFloat() > continueProbability)
                continue;
            Llight /= continueProbability;
        }

        // Add contribution from _VirtualLight_ _vl_
        if (!scene->IntersectP(connectRay))
            L += Llight;
    }
    if (ray.depth < maxSpecularDepth) {
        // Do bias compensation for bounding geometry term
        int nSamples = (ray.depth == 0) ? nGatherSamples : 1;
        for (int i = 0; i < nSamples; ++i) {
            Vector wi;
            float pdf;
            BSDFSample bsdfSample = (ray.depth == 0) ?
                BSDFSample(sample, gatherSampleOffset, i) : BSDFSample(rng);
            Spectrum f = bsdf->Sample_f(wo, &wi, bsdfSample,
                                        &pdf, BxDFType(BSDF_ALL & ~BSDF_SPECULAR));
            if (!f.IsBlack() && pdf > 0.f) {
                // Trace ray for bias compensation gather sample
                float maxDist = sqrtf(AbsDot(wi, n) / gLimit);
                RayDifferential gatherRay(p, wi, ray, isect.rayEpsilon, maxDist);
                Intersection gatherIsect;
                Spectrum Li = renderer->Li(scene, gatherRay, sample, rng, arena,
                                           &gatherIsect);
                if (Li.IsBlack()) continue;

                // Add bias compensation ray contribution to radiance sum
                float Ggather = AbsDot(wi, n) * AbsDot(-wi, gatherIsect.dg.nn) /
                    DistanceSquared(p, gatherIsect.dg.p);
                if (Ggather - gLimit > 0.f && !isinf(Ggather)) {
                    float gs = (Ggather - gLimit) / Ggather;
                    L += f * Li * (AbsDot(wi, n) * gs / (nSamples * pdf));
                }
            }
        }
    }
    if (ray.depth + 1 < maxSpecularDepth) {
        Vector wi;
        // Trace rays for specular reflection and refraction
        L += SpecularReflect(ray, bsdf, rng, isect, renderer, scene, sample,
                             arena);
        L += SpecularTransmit(ray, bsdf, rng, isect, renderer, scene, sample,
                              arena);
    }
    return L;
}

Spectrum ExPhotonIntegrator::Li(const Scene *scene,
		const RayDifferential &ray, const Sample *sample,
		float *alpha) const {
	// Compute reflected radiance with photon map
	Spectrum L(0.);
	Intersection isect;
	if (scene->Intersect(ray, &isect)) {
		if (alpha) *alpha = 1.;
		Vector wo = -ray.d;
		// Compute emitted light if ray hit an area light source
		L += isect.Le(wo);
		// Evaluate BSDF at hit point
		BSDF *bsdf = isect.GetBSDF(ray);
		const Point &p = bsdf->dgShading.p;
		const Normal &n = bsdf->dgShading.nn;
		L += UniformSampleAllLights(scene, p, n,
			wo, bsdf, sample,
			lightSampleOffset, bsdfSampleOffset,
			bsdfComponentOffset);

		// Compute indirect lighting for photon map integrator
		L += LPhoton(causticMap, nCausticPaths, nLookup, bsdf,
			isect, wo, maxDistSquared);
		if (finalGather) {
#if 1
			// Do one-bounce final gather for photon map
			BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
				BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
			if (bsdf->NumComponents(nonSpecular) > 0) {
				// Find indirect photons around point for importance sampling
				u_int nIndirSamplePhotons = 50;
				PhotonProcess proc(nIndirSamplePhotons, p);
				proc.photons = (ClosePhoton *)alloca(nIndirSamplePhotons *
					sizeof(ClosePhoton));
				float searchDist2 = maxDistSquared;
				while (proc.foundPhotons < nIndirSamplePhotons) {
					float md2 = searchDist2;
					proc.foundPhotons = 0;
					indirectMap->Lookup(p, proc, md2);
					searchDist2 *= 2.f;
				}
				// Copy photon directions to local array
				Vector *photonDirs = (Vector *)alloca(nIndirSamplePhotons *
					sizeof(Vector));
				for (u_int i = 0; i < nIndirSamplePhotons; ++i)
					photonDirs[i] = proc.photons[i].photon->wi;
				// Use BSDF to do final gathering
				Spectrum Li = 0.;
				static StatsCounter gatherRays("Photon Map", // NOBOOK
					"Final gather rays traced"); // NOBOOK
				for (int i = 0; i < gatherSamples; ++i) {
					// Sample random direction from BSDF for final gather ray
					Vector wi;
					float u1 = sample->twoD[gatherSampleOffset[0]][2*i];
					float u2 = sample->twoD[gatherSampleOffset[0]][2*i+1];
					float u3 = sample->oneD[gatherComponentOffset[0]][i];
					float pdf;
					Spectrum fr = bsdf->Sample_f(wo, &wi, u1, u2, u3,
						&pdf, BxDFType(BSDF_ALL & (~BSDF_SPECULAR)));
					if (fr.Black() || pdf == 0.f) continue;
					// Trace BSDF final gather ray and accumulate radiance
					RayDifferential bounceRay(p, wi);
					++gatherRays; // NOBOOK
					Intersection gatherIsect;
					if (scene->Intersect(bounceRay, &gatherIsect)) {
						// Compute exitant radiance using precomputed irradiance
						Spectrum Lindir = 0.f;
						Normal n = gatherIsect.dg.nn;
						if (Dot(n, bounceRay.d) > 0) n = -n;
						RadiancePhotonProcess proc(gatherIsect.dg.p, n);
						float md2 = INFINITY;
						radianceMap->Lookup(gatherIsect.dg.p, proc, md2);
						if (proc.photon)
							Lindir = proc.photon->Lo;
						Lindir *= scene->Transmittance(bounceRay);
						// Compute MIS weight for BSDF-sampled gather ray
						// Compute PDF for photon-sampling of direction _wi_
						float photonPdf = 0.f;
						float conePdf = UniformConePdf(cosGatherAngle);
						for (u_int j = 0; j < nIndirSamplePhotons; ++j)
							if (Dot(photonDirs[j], wi) > .999f * cosGatherAngle)
								photonPdf += conePdf;
						photonPdf /= nIndirSamplePhotons;
						float wt = PowerHeuristic(gatherSamples, pdf,
							gatherSamples, photonPdf);
						Li += fr * Lindir * AbsDot(wi, n) * wt / pdf;
					}
				}
				L += Li / gatherSamples;
				// Use nearby photons to do final gathering
				Li = 0.;
				for (int i = 0; i < gatherSamples; ++i) {
					// Sample random direction using photons for final gather ray
					float u1 = sample->oneD[gatherComponentOffset[1]][i];
					float u2 = sample->twoD[gatherSampleOffset[1]][2*i];
					float u3 = sample->twoD[gatherSampleOffset[1]][2*i+1];
					int photonNum = min((int)nIndirSamplePhotons - 1,
						Floor2Int(u1 * nIndirSamplePhotons));
					// Sample gather ray direction from _photonNum_
					Vector vx, vy;
//.........这里部分代码省略.........

void IGIIntegrator::Preprocess(const Scene *scene) {
	if (scene->lights.size() == 0) return;
	// Compute samples for emitted rays from lights
	float *lightNum = new float[nLightPaths * nLightSets];
	float *lightSamp0 = new float[2 * nLightPaths *	nLightSets];
	float *lightSamp1 = new float[2 * nLightPaths * nLightSets];
	LDShuffleScrambled1D(nLightPaths, nLightSets, lightNum);
	LDShuffleScrambled2D(nLightPaths, nLightSets, lightSamp0);
	LDShuffleScrambled2D(nLightPaths, nLightSets, lightSamp1);
	// Precompute information for light sampling densities
	int nLights = int(scene->lights.size());
	float *lightPower = (float *)alloca(nLights * sizeof(float));
	float *lightCDF = (float *)alloca((nLights+1) * sizeof(float));
	for (int i = 0; i < nLights; ++i)
		lightPower[i] = scene->lights[i]->Power(scene).y();
	float totalPower;
	ComputeStep1dCDF(lightPower, nLights, &totalPower, lightCDF);
	for (u_int s = 0; s < nLightSets; ++s) {
		for (u_int i = 0; i < nLightPaths; ++i) {
			// Follow path _i_ from light to create virtual lights
			int sampOffset = s*nLightPaths + i;
			// Choose light source to trace path from
			float lightPdf;
			int lNum = Floor2Int(SampleStep1d(lightPower, lightCDF,
				totalPower, nLights, lightNum[sampOffset], &lightPdf) * nLights);
//			fprintf(stderr, "samp %f -> num %d\n", lightNum[sampOffset], lNum);
			Light *light = scene->lights[lNum];
			// Sample ray leaving light source
			RayDifferential ray;
			float pdf;
			Spectrum alpha =
				light->Sample_L(scene, lightSamp0[2*sampOffset],
						lightSamp0[2*sampOffset+1],
						lightSamp1[2*sampOffset],
						lightSamp1[2*sampOffset+1],
						&ray, &pdf);
			if (pdf == 0.f || alpha.Black()) continue;
			alpha /= pdf * lightPdf;
//			fprintf(stderr, "initial alpha %f, light # %d\n", alpha.y(), lNum);
			Intersection isect;
			int nIntersections = 0;
			while (scene->Intersect(ray, &isect) && !alpha.Black()) {
				++nIntersections;
				alpha *= scene->Transmittance(ray);
				Vector wo = -ray.d;
				BSDF *bsdf = isect.GetBSDF(ray);
				// Create virtual light at ray intersection point
				static StatsCounter vls("IGI Integrator", "Virtual Lights Created"); //NOBOOK
				++vls; //NOBOOK
				Spectrum Le = alpha * bsdf->rho(wo) / M_PI;
//				fprintf(stderr, "\tmade light with le y %f\n", Le.y());
				virtualLights[s].push_back(VirtualLight(isect.dg.p, isect.dg.nn, Le));
				// Sample new ray direction and update weight
				Vector wi;
				float pdf;
				BxDFType flags;
				Spectrum fr = bsdf->Sample_f(wo, &wi, RandomFloat(),
								 RandomFloat(), RandomFloat(),
								 &pdf, BSDF_ALL, &flags);
				if (fr.Black() || pdf == 0.f)
					break;
				Spectrum anew = alpha * fr * AbsDot(wi, bsdf->dgShading.nn) / pdf;
				float r = anew.y() / alpha.y();
//				fprintf(stderr, "\tr = %f\n", r);
				if (RandomFloat() > r)
					break;
				alpha = anew / r;
//				fprintf(stderr, "\tnew alpha %f\n", alpha.y());
				ray = RayDifferential(isect.dg.p, wi);
			}
			BSDF::FreeAll();
		}
	}
	delete[] lightNum; // NOBOOK
	delete[] lightSamp0; // NOBOOK
	delete[] lightSamp1; // NOBOOK
}

//.........这里部分代码省略.........
								}
								progress.Update();
							}
						}
						else {
							// Process indirect lighting photon intersection
							if (!indirectDone) {
								indirectPhotons.push_back(photon);
								if (indirectPhotons.size() == nIndirectPhotons) {
									indirectDone = true;
									nIndirectPaths = (int)nshot;
									indirectMap = new KdTree<Photon, PhotonProcess>(indirectPhotons);
								}
								progress.Update();
							}
						}
					}
					if (finalGather && RandomFloat() < .125f) {
						// Store data for radiance photon
						static StatsCounter rp("Photon Map", "Radiance photons created"); // NOBOOK
						++rp; // NOBOOK
						Normal n = photonIsect.dg.nn;
						if (Dot(n, photonRay.d) > 0.f) n = -n;
						radiancePhotons.push_back(RadiancePhoton(photonIsect.dg.p, n));
						Spectrum rho_r = photonBSDF->rho(BSDF_ALL_REFLECTION);
						rpReflectances.push_back(rho_r);
						Spectrum rho_t = photonBSDF->rho(BSDF_ALL_TRANSMISSION);
						rpTransmittances.push_back(rho_t);
					}
				}
				// Sample new photon ray direction
				Vector wi;
				float pdf;
				BxDFType flags;
				// Get random numbers for sampling outgoing photon direction
				float u1, u2, u3;
				if (nIntersections == 1) {
					u1 = RadicalInverse((int)nshot+1, 13);
					u2 = RadicalInverse((int)nshot+1, 17);
					u3 = RadicalInverse((int)nshot+1, 19);
				}
				else {
					u1 = RandomFloat();
					u2 = RandomFloat();
					u3 = RandomFloat();
				}

				// Compute new photon weight and possibly terminate with RR
				Spectrum fr = photonBSDF->Sample_f(wo, &wi, u1, u2, u3,
				                                   &pdf, BSDF_ALL, &flags);
				if (fr.Black() || pdf == 0.f)
					break;
				Spectrum anew = alpha * fr *
					AbsDot(wi, photonBSDF->dgShading.nn) / pdf;
				float continueProb = min(1.f, anew.y() / alpha.y());
				if (RandomFloat() > continueProb || nIntersections > 10)
					break;
				alpha = anew / continueProb;
				specularPath = (nIntersections == 1 || specularPath) &&
					((flags & BSDF_SPECULAR) != 0);
				photonRay = RayDifferential(photonIsect.dg.p, wi);
			}
		}
		BSDF::FreeAll();
	}

	progress.Done(); // NOBOOK

	// Precompute radiance at a subset of the photons
	KdTree<Photon, PhotonProcess> directMap(directPhotons);
	int nDirectPaths = nshot;
	if (finalGather) {
		ProgressReporter p2(radiancePhotons.size(), "Computing photon radiances"); // NOBOOK
		for (u_int i = 0; i < radiancePhotons.size(); ++i) {
			// Compute radiance for radiance photon _i_
			RadiancePhoton &rp = radiancePhotons[i];
			const Spectrum &rho_r = rpReflectances[i];
			const Spectrum &rho_t = rpTransmittances[i];
			Spectrum E;
			Point p = rp.p;
			Normal n = rp.n;
			if (!rho_r.Black()) {
				E = estimateE(&directMap,  nDirectPaths,   p, n) +
					estimateE(indirectMap, nIndirectPaths, p, n) +
					estimateE(causticMap,  nCausticPaths,  p, n);
				rp.Lo += E * INV_PI * rho_r;
			}
			if (!rho_t.Black()) {
				E = estimateE(&directMap,  nDirectPaths,   p, -n) +
					estimateE(indirectMap, nIndirectPaths, p, -n) +
					estimateE(causticMap,  nCausticPaths,  p, -n);
				rp.Lo += E * INV_PI * rho_t;
			}
			p2.Update(); // NOBOOK
		}
		radianceMap = new KdTree<RadiancePhoton,
			RadiancePhotonProcess>(radiancePhotons);
		p2.Done(); // NOBOOK
	}
}

Spectrum IGIIntegrator::Li(const Scene *scene,
		const RayDifferential &ray, const Sample *sample,
		   float *alpha) const {
	Spectrum L(0.);
	Intersection isect;
	if (scene->Intersect(ray, &isect)) {
		if (alpha) *alpha = 1.;
		Vector wo = -ray.d;
		// Compute emitted light if ray hit an area light source
		L += isect.Le(wo);
		// Evaluate BSDF at hit point
		BSDF *bsdf = isect.GetBSDF(ray);
		const Point &p = bsdf->dgShading.p;
		const Normal &n = bsdf->dgShading.nn;
		L += UniformSampleAllLights(scene, p, n,
					    wo, bsdf, sample,
					    lightSampleOffset, bsdfSampleOffset,
					    bsdfComponentOffset);
		// Compute indirect illumination with virtual lights
		u_int lSet = min(u_int(sample->oneD[vlSetOffset][0] * nLightSets),
		                 nLightSets-1);
		for (u_int i = 0; i < virtualLights[lSet].size(); ++i) {
			const VirtualLight &vl = virtualLights[lSet][i];
			// Add contribution from _VirtualLight_ _vl_
			// Ignore light if it's too close to current point
			float d2 = DistanceSquared(p, vl.p);
			//if (d2 < .8 * minDist2) continue;
			float distScale = SmoothStep(.8 * minDist2, 1.2 * minDist2, d2);
			// Compute virtual light's tentative contribution _Llight_
			Vector wi = Normalize(vl.p - p);
			Spectrum f = distScale * bsdf->f(wo, wi);
			if (f.Black()) continue;
			float G = AbsDot(wi, n) * AbsDot(wi, vl.n) / d2;
			Spectrum Llight = indirectScale * f * G * vl.Le /
				virtualLights[lSet].size();
			Llight *= scene->Transmittance(Ray(p, vl.p - p));
			// Possibly skip shadow ray with Russian roulette
			if (Llight.y() < rrThreshold) {
				float continueProbability = .1f;
				if (RandomFloat() > continueProbability)
					continue;
				Llight /= continueProbability;
			}
			static StatsCounter vlsr("IGI Integrator", "Shadow Rays to Virtual Lights"); //NOBOOK
			++vlsr; //NOBOOK
			if (!scene->IntersectP(Ray(p, vl.p - p, RAY_EPSILON,
					1.f - RAY_EPSILON)))
				L += Llight;
		}
		// Trace rays for specular reflection and refraction
		if (specularDepth++ < maxSpecularDepth) {
			Vector wi;
			// Trace rays for specular reflection and refraction
			Spectrum f = bsdf->Sample_f(wo, &wi,
						BxDFType(BSDF_REFLECTION | BSDF_SPECULAR));
			if (!f.Black()) {
				// Compute ray differential _rd_ for specular reflection
				RayDifferential rd(p, wi);
				rd.hasDifferentials = true;
				rd.rx.o = p + isect.dg.dpdx;
				rd.ry.o = p + isect.dg.dpdy;
				// Compute differential reflected directions
				Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx +
					bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
				Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy +
					bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
				Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
				float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
				float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
				rd.rx.d = wi -
					dwodx + 2 * Vector(Dot(wo, n) * dndx +
						 dDNdx * n);
				rd.ry.d = wi -
					dwody + 2 * Vector(Dot(wo, n) * dndy +
						 dDNdy * n);
				L += scene->Li(rd, sample) * f * AbsDot(wi, n);
			}
			f = bsdf->Sample_f(wo, &wi,
					   BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));
			if (!f.Black()) {
				// Compute ray differential _rd_ for specular transmission
				RayDifferential rd(p, wi);
				rd.hasDifferentials = true;
				rd.rx.o = p + isect.dg.dpdx;
				rd.ry.o = p + isect.dg.dpdy;

				float eta = bsdf->eta;
				Vector w = -wo;
				if (Dot(wo, n) < 0) eta = 1.f / eta;

				Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx + bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
				Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy + bsdf->dgShading.dndv * bsdf->dgShading.dvdy;

				Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
				float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
				float dDNdy = Dot(dwody, n) + Dot(wo, dndy);

				float mu = eta * Dot(w, n) - Dot(wi, n);
				float dmudx = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdx;
				float dmudy = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdy;
//.........这里部分代码省略.........

Spectrum PhotonIntegrator::Li(const Scene *scene,
		const RayDifferential &ray, const Sample *sample,
		float *alpha) const {
	// Compute reflected radiance with photon map
	Spectrum L(0.);
	Intersection isect;
	if (scene->Intersect(ray, &isect)) {
		if (alpha) *alpha = 1.;
		Vector wo = -ray.d;
		// Compute emitted light if ray hit an area light source
		L += isect.Le(wo);
		// Evaluate BSDF at hit point
		BSDF *bsdf = isect.GetBSDF(ray);
		const Point &p = bsdf->dgShading.p;
		const Normal &n = bsdf->dgShading.nn;
		// Compute direct lighting for photon map integrator
		if (directWithPhotons)
			L += LPhoton(directMap, nDirectPaths, nLookup,
				bsdf, isect, wo, maxDistSquared);
		else
			L += UniformSampleAllLights(scene, p, n,
				wo, bsdf, sample,
				lightSampleOffset, bsdfSampleOffset,
				bsdfComponentOffset);
		
		// Compute indirect lighting for photon map integrator
		L += LPhoton(causticMap, nCausticPaths, nLookup, bsdf,
			isect, wo, maxDistSquared);
		if (finalGather) {
			// Do one-bounce final gather for photon map
			Spectrum Li(0.);
			for (int i = 0; i < gatherSamples; ++i) {
				// Sample random direction for final gather ray
				Vector wi;
				float u1 = sample->twoD[gatherSampleOffset][2*i];
				float u2 = sample->twoD[gatherSampleOffset][2*i+1];
				float u3 = sample->oneD[gatherComponentOffset][i];
				float pdf;
				Spectrum fr = bsdf->Sample_f(wo, &wi, u1, u2, u3,
					&pdf, BxDFType(BSDF_ALL & (~BSDF_SPECULAR)));
				if (fr.Black() || pdf == 0.f) continue;
				RayDifferential bounceRay(p, wi);
				static StatsCounter gatherRays("Photon Map", // NOBOOK
					"Final gather rays traced"); // NOBOOK
				++gatherRays; // NOBOOK
				Intersection gatherIsect;
				if (scene->Intersect(bounceRay, &gatherIsect)) {
					// Compute exitant radiance at final gather intersection
					BSDF *gatherBSDF = gatherIsect.GetBSDF(bounceRay);
					Vector bounceWo = -bounceRay.d;
					Spectrum Lindir =
						LPhoton(directMap, nDirectPaths, nLookup,
							gatherBSDF, gatherIsect, bounceWo, maxDistSquared) +
						LPhoton(indirectMap, nIndirectPaths, nLookup,
							gatherBSDF, gatherIsect, bounceWo, maxDistSquared) +
						LPhoton(causticMap, nCausticPaths, nLookup,
							gatherBSDF, gatherIsect, bounceWo, maxDistSquared);
					Lindir *= scene->Transmittance(bounceRay);
					Li += fr * Lindir * AbsDot(wi, n) / pdf;
				}
			}
			L += Li / float(gatherSamples);
		}
		else
			L += LPhoton(indirectMap, nIndirectPaths, nLookup,
				bsdf, isect, wo, maxDistSquared);
		if (specularDepth++ < maxSpecularDepth) {
			Vector wi;
			// Trace rays for specular reflection and refraction
			Spectrum f = bsdf->Sample_f(wo, &wi,
				BxDFType(BSDF_REFLECTION | BSDF_SPECULAR));
			if (!f.Black()) {
				// Compute ray differential _rd_ for specular reflection
				RayDifferential rd(p, wi);
				rd.hasDifferentials = true;
				rd.rx.o = p + isect.dg.dpdx;
				rd.ry.o = p + isect.dg.dpdy;
				// Compute differential reflected directions
				Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx +
					bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
				Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy +
					bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
				Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
				float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
				float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
				rd.rx.d = wi -
				          dwodx + 2 * Vector(Dot(wo, n) * dndx +
						  dDNdx * n);
				rd.ry.d = wi -
				          dwody + 2 * Vector(Dot(wo, n) * dndy +
						  dDNdy * n);
				L += scene->Li(rd, sample) * f * AbsDot(wi, n);
			}
			f = bsdf->Sample_f(wo, &wi,
				BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));
			if (!f.Black()) {
				// Compute ray differential _rd_ for specular transmission
				RayDifferential rd(p, wi);
				rd.hasDifferentials = true;
				rd.rx.o = p + isect.dg.dpdx;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		if (pdf == 0.f || alpha.Black()) continue;
		alpha /= pdf * lightPdf;
		if (!alpha.Black()) {
			// Follow photon path through scene and record intersections
			bool specularPath = false;
			Intersection photonIsect;
			int nIntersections = 0;
			while (scene->Intersect(photonRay, &photonIsect)) {
				++nIntersections;
				// Handle photon/surface intersection
				alpha *= scene->Transmittance(photonRay);
				Vector wo = -photonRay.d;
				BSDF *photonBSDF = photonIsect.GetBSDF(photonRay);
				BxDFType specularType = BxDFType(BSDF_REFLECTION |
					BSDF_TRANSMISSION | BSDF_SPECULAR);
				bool hasNonSpecular = (photonBSDF->NumComponents() >
					photonBSDF->NumComponents(specularType));
				if (hasNonSpecular) {
					// Deposit photon at surface
					Photon photon(photonIsect.dg.p, alpha, wo);
					if (nIntersections == 1) {
						// Process direct lighting photon intersection
						if (!directDone) {
							directPhotons.push_back(photon);
							if (directPhotons.size() == nDirectPhotons) {
								directDone = true;
								nDirectPaths = (int)nshot;
								directMap =
									new KdTree<Photon,
											   PhotonProcess>(directPhotons);
							}
							progress.Update(); // NOBOOK
						}
					}
					else if (specularPath) {
						// Process caustic photon intersection
						if (!causticDone) {
							causticPhotons.push_back(photon);
							if (causticPhotons.size() == nCausticPhotons) {
								causticDone = true;
								nCausticPaths = (int)nshot;
								causticMap =
									new KdTree<Photon,
										       PhotonProcess>(causticPhotons);
							}
							progress.Update();
						}
					}
					else {
						// Process indirect lighting photon intersection
						if (!indirectDone) {
							indirectPhotons.push_back(photon);
							if (indirectPhotons.size() == nIndirectPhotons) {
								indirectDone = true;
								nIndirectPaths = (int)nshot;
								indirectMap =
									new KdTree<Photon,
											   PhotonProcess>(indirectPhotons);
							}
							progress.Update();
						}
					}
				}
				// Sample new photon ray direction
				Vector wi;
				float pdf;
				BxDFType flags;
				// Get random numbers for sampling outgoing photon direction
				float u1, u2, u3;
				if (nIntersections == 1) {
					u1 = (float)RadicalInverse((int)nshot+1, 13);
					u2 = (float)RadicalInverse((int)nshot+1, 17);
					u3 = (float)RadicalInverse((int)nshot+1, 19);
				}
				else {
					u1 = RandomFloat();
					u2 = RandomFloat();
					u3 = RandomFloat();
				}
				Spectrum fr = photonBSDF->Sample_f(wo, &wi, u1, u2, u3,
					&pdf, BSDF_ALL, &flags);
				if (fr.Black() || pdf == 0.f)
					break;
				specularPath = (nIntersections == 1 || specularPath) &&
					((flags & BSDF_SPECULAR) != 0);
				alpha *= fr * AbsDot(wi, photonBSDF->dgShading.nn) / pdf;
				photonRay = RayDifferential(photonIsect.dg.p, wi);
				// Possibly terminate photon path
				if (nIntersections > 3) {
					float continueProbability = .5f;
					if (RandomFloat() > continueProbability)
						break;
					alpha /= continueProbability;
				}
			}
		}
		BSDF::FreeAll();
	}
	progress.Done(); // NOBOOK
}

Spectrum IrradianceCache::Li(const Scene *scene, const RayDifferential &ray,
		const Sample *sample, float *alpha) const {
	Intersection isect;
	Spectrum L(0.);
	if (scene->Intersect(ray, &isect)) {
		if (alpha) *alpha = 1.;
		// Evaluate BSDF at hit point
		BSDF *bsdf = isect.GetBSDF(ray);
		Vector wo = -ray.d;
		const Point &p = bsdf->dgShading.p;
		const Normal &n = bsdf->dgShading.nn;
		// Compute direct lighting for irradiance cache
		L += isect.Le(wo);
		L += UniformSampleAllLights(scene, p, n, wo, bsdf, sample,
			lightSampleOffset, bsdfSampleOffset,
			bsdfComponentOffset);
		// Compute indirect lighting for irradiance cache
		if (specularDepth++ < maxSpecularDepth) {
			Vector wi;
			// Trace rays for specular reflection and refraction
			Spectrum f = bsdf->Sample_f(wo, &wi,
				BxDFType(BSDF_REFLECTION | BSDF_SPECULAR));
			if (!f.Black()) {
				// Compute ray differential _rd_ for specular reflection
				RayDifferential rd(p, wi);
				rd.hasDifferentials = true;
				rd.rx.o = p + isect.dg.dpdx;
				rd.ry.o = p + isect.dg.dpdy;
				// Compute differential reflected directions
				Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx +
					bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
				Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy +
					bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
				Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
				float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
				float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
				rd.rx.d = wi -
				          dwodx + 2 * Vector(Dot(wo, n) * dndx +
						  dDNdx * n);
				rd.ry.d = wi -
				          dwody + 2 * Vector(Dot(wo, n) * dndy +
						  dDNdy * n);
				L += scene->Li(rd, sample) * f * AbsDot(wi, n);
			}
			f = bsdf->Sample_f(wo, &wi,
				BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));
			if (!f.Black()) {
				// Compute ray differential _rd_ for specular transmission
				RayDifferential rd(p, wi);
				rd.hasDifferentials = true;
				rd.rx.o = p + isect.dg.dpdx;
				rd.ry.o = p + isect.dg.dpdy;
				
				float eta = bsdf->eta;
				Vector w = -wo;
				if (Dot(wo, n) < 0) eta = 1.f / eta;
				
				Normal dndx = bsdf->dgShading.dndu * bsdf->dgShading.dudx + bsdf->dgShading.dndv * bsdf->dgShading.dvdx;
				Normal dndy = bsdf->dgShading.dndu * bsdf->dgShading.dudy + bsdf->dgShading.dndv * bsdf->dgShading.dvdy;
				
				Vector dwodx = -ray.rx.d - wo, dwody = -ray.ry.d - wo;
				float dDNdx = Dot(dwodx, n) + Dot(wo, dndx);
				float dDNdy = Dot(dwody, n) + Dot(wo, dndy);
				
				float mu = eta * Dot(w, n) - Dot(wi, n);
				float dmudx = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdx;
				float dmudy = (eta - (eta*eta*Dot(w,n))/Dot(wi, n)) * dDNdy;
				
				rd.rx.d = wi + eta * dwodx - Vector(mu * dndx + dmudx * n);
				rd.ry.d = wi + eta * dwody - Vector(mu * dndy + dmudy * n);
				L += scene->Li(rd, sample) * f * AbsDot(wi, n);
			}
		}
		--specularDepth;
		// Estimate indirect lighting with irradiance cache
		Normal ng = isect.dg.nn;
		if (Dot(wo, ng) < 0.f) ng = -ng;
		BxDFType flags = BxDFType(BSDF_REFLECTION |
		                          BSDF_DIFFUSE |
								  BSDF_GLOSSY);
		L += IndirectLo(p, ng, wo, bsdf, flags, sample, scene);
		flags = BxDFType(BSDF_TRANSMISSION |
		                 BSDF_DIFFUSE |
						 BSDF_GLOSSY);
		L += IndirectLo(p, -ng, wo, bsdf, flags, sample, scene);
	}
	else {
		// Handle ray with no intersection
		if (alpha) *alpha = 0.;
		for (u_int i = 0; i < scene->lights.size(); ++i)
			L += scene->lights[i]->Le(ray);
		if (alpha && !L.Black()) *alpha = 1.;
		return L;
	}
	return L;
}

Spectrum GlossyPRTIntegrator::Li(const Scene *scene, const Renderer *,
        const RayDifferential &ray, const Intersection &isect,
        const Sample *sample, RNG &rng, MemoryArena &arena) const {
    Spectrum L = 0.f;
    Vector wo = -ray.d;
    // Compute emitted light if ray hit an area light source
    L += isect.Le(wo);

    // Evaluate BSDF at hit point
    BSDF *bsdf = isect.GetBSDF(ray, arena);
    const Point &p = bsdf->dgShading.p;
    // Compute reflected radiance with glossy PRT at point

    // Compute SH radiance transfer matrix at point and SH coefficients
    Spectrum *c_t = arena.Alloc<Spectrum>(SHTerms(lmax));
    Spectrum *T = arena.Alloc<Spectrum>(SHTerms(lmax)*SHTerms(lmax));
    SHComputeTransferMatrix(p, isect.rayEpsilon, scene, rng, nSamples,
                            lmax, T);
    SHMatrixVectorMultiply(T, c_in, c_t, lmax);

    // Rotate incident SH lighting to local coordinate frame
    Vector r1 = bsdf->LocalToWorld(Vector(1,0,0));
    Vector r2 = bsdf->LocalToWorld(Vector(0,1,0));
    Normal nl = Normal(bsdf->LocalToWorld(Vector(0,0,1)));
    Matrix4x4 rot(r1.x, r2.x, nl.x, 0,
                  r1.y, r2.y, nl.y, 0,
                  r1.z, r2.z, nl.z, 0,
                     0,    0,    0, 1);
    Spectrum *c_l = arena.Alloc<Spectrum>(SHTerms(lmax));
    SHRotate(c_t, c_l, rot, lmax, arena);
    #if 0

    // Sample BSDF and integrate against direct SH coefficients
    float *Ylm = ALLOCA(float, SHTerms(lmax));
    int ns = 1024;
    for (int i = 0; i < ns; ++i) {
        Vector wi;
        float pdf;
        Spectrum f = bsdf->Sample_f(wo, &wi, BSDFSample(rng), &pdf);
        if (pdf > 0.f && !f.IsBlack() && !scene->IntersectP(Ray(p, wi))) {
            f *= fabsf(Dot(wi, n)) / (pdf * ns);
            SHEvaluate(bsdf->WorldToLocal(wi), lmax, Ylm);
    
            Spectrum Li = 0.f;
            for (int j = 0; j < SHTerms(lmax); ++j)
                Li += Ylm[j] * c_l[j] * f;
            L += Li.Clamp();
        }
    }
    #else

    // Compute final coefficients _c\_o_ using BSDF matrix
    Spectrum *c_o = arena.Alloc<Spectrum>(SHTerms(lmax));
    SHMatrixVectorMultiply(B, c_l, c_o, lmax);

    // Evaluate outgoing radiance function for $\wo$ and add to _L_
    Vector woLocal = bsdf->WorldToLocal(wo);
    float *Ylm = ALLOCA(float, SHTerms(lmax));
    SHEvaluate(woLocal, lmax, Ylm);
    Spectrum Li = 0.f;
    for (int i = 0; i < SHTerms(lmax); ++i)
        Li += Ylm[i] * c_o[i];
    L += Li.Clamp();
    #endif
    return L;
}

Spectrum IrradianceCache::IndirectLo(const Point &p,
		const Normal &n, const Vector &wo, BSDF *bsdf,
		BxDFType flags, const Sample *sample,
		const Scene *scene) const {
	if (bsdf->NumComponents(flags) == 0)
		return Spectrum(0.);
	Spectrum E;
	if (!InterpolateIrradiance(scene, p, n, &E)) {
		// Compute irradiance at current point
		u_int scramble[2] = { RandomUInt(), RandomUInt() };
		float sumInvDists = 0.;
		for (int i = 0; i < nSamples; ++i) {
			// Trace ray to sample radiance for irradiance estimate
			// Update irradiance statistics for rays traced
			static StatsCounter nIrradiancePaths("Irradiance Cache",
				"Paths followed for irradiance estimates");
			++nIrradiancePaths;
			float u[2];
			Sample02(i, scramble, u);
			Vector w = CosineSampleHemisphere(u[0], u[1]);
			RayDifferential r(p, bsdf->LocalToWorld(w));
			if (Dot(r.d, n) < 0) r.d = -r.d;
			Spectrum L(0.);
			// Do path tracing to compute radiance along ray for estimate
			{
			// Declare common path integration variables
			Spectrum pathThroughput = 1.;
			RayDifferential ray(r);
			bool specularBounce = false;
			for (int pathLength = 0; ; ++pathLength) {
				// Find next vertex of path
				Intersection isect;
				if (!scene->Intersect(ray, &isect))
					break;
				if (pathLength == 0)
					r.maxt = ray.maxt;
				pathThroughput *= scene->Transmittance(ray);
				// Possibly add emitted light at path vertex
				if (specularBounce)
					L += pathThroughput * isect.Le(-ray.d);
				// Evaluate BSDF at hit point
				BSDF *bsdf = isect.GetBSDF(ray);
				// Sample illumination from lights to find path contribution
				const Point &p = bsdf->dgShading.p;
				const Normal &n = bsdf->dgShading.nn;
				Vector wo = -ray.d;
				L += pathThroughput *
					UniformSampleOneLight(scene, p, n, wo, bsdf, sample);
				if (pathLength+1 == maxIndirectDepth) break;
				// Sample BSDF to get new path direction
				// Get random numbers for sampling new direction, \mono{bs1}, \mono{bs2}, and \mono{bcs}
				float bs1 = RandomFloat(), bs2 = RandomFloat(), bcs = RandomFloat();
				Vector wi;
				float pdf;
				BxDFType flags;
				Spectrum f = bsdf->Sample_f(wo, &wi, bs1, bs2, bcs,
					&pdf, BSDF_ALL, &flags);
				if (f.Black() || pdf == 0.)
					break;
				specularBounce = (flags & BSDF_SPECULAR) != 0;
				pathThroughput *= f * AbsDot(wi, n) / pdf;
				ray = RayDifferential(p, wi);
				// Possibly terminate the path
				if (pathLength > 3) {
					float continueProbability = .5f;
					if (RandomFloat() > continueProbability)
						break;
					pathThroughput /= continueProbability;
				}
			}
			}
			E += L;
			float dist = r.maxt * r.d.Length();
			sumInvDists += 1.f / dist;
		}
		E *= M_PI / float(nSamples);
		// Add computed irradiance value to cache
		// Update statistics for new irradiance sample
		static StatsCounter nSamplesComputed("Irradiance Cache",
			"Irradiance estimates computed");
		++nSamplesComputed;
		// Compute bounding box of irradiance sample's contribution region
		static float minMaxDist =
			.001f * powf(scene->WorldBound().Volume(), 1.f/3.f);
		static float maxMaxDist =
			.125f * powf(scene->WorldBound().Volume(), 1.f/3.f);
		float maxDist = nSamples / sumInvDists;
		if (minMaxDist > 0.f)
			maxDist = Clamp(maxDist, minMaxDist, maxMaxDist);
		maxDist *= maxError;
		BBox sampleExtent(p);
		sampleExtent.Expand(maxDist);
		octree->Add(IrradianceSample(E, p, n, maxDist),
			sampleExtent);
	}
	return .5f * bsdf->rho(wo, flags) * E;
}