C++ VisibilityTester类代码示例

Spectrum EstimateIrradianceDirect(const Scene* scene,
								  const Light* light, const Point& p,
								  const Normal& n,
								  int lightSamp,
								  u_int sampleNum)
{
	Spectrum Ed(0.); // direct irradiance
	// Find light and BSDF sample values for direct lighting estimate
	float ls1, ls2;
	
	ls1 = RandomFloat();
	ls2 = RandomFloat();
	
	// Sample light source with multiple importance sampling
	Vector wi;
	float lightPdf;
	VisibilityTester visibility;
	Spectrum Li = light->Sample_L(p, n, ls1, ls2, &wi, &lightPdf, &visibility);
	//printf("got light sample ");
	//Li.printSelf();
	if (lightPdf > 0. && !Li.Black()) {
		if (visibility.Unoccluded(scene)) {
			// Add light's contribution to reflected radiance
			Li *= visibility.Transmittance(scene);
			Ed += Li * AbsDot(wi, n) / lightPdf;
		}
	}
	return Ed;
}

Spectrum EstimateDirect(const Scene *scene, const Renderer *renderer,
        MemoryArena &arena, const Light *light, const Point &p,
        const Normal &n, const Vector &wo, float rayEpsilon, float time,
        const BSDF *bsdf, RNG &rng, const LightSample &lightSample,
        const BSDFSample &bsdfSample, BxDFType flags) {
    Spectrum Ld(0.);
    // Sample light source with multiple importance sampling
    Vector wi;
    float lightPdf, bsdfPdf;
    VisibilityTester visibility;
    light->Sample_L(p, rayEpsilon, lightSample, time,
                                  &wi, &lightPdf, &visibility);
    Spectrum Li(1.);
    if (lightPdf > 0. && !Li.IsBlack()) {
        Spectrum f = bsdf->f(wo, wi, flags);
        if (!f.IsBlack() && visibility.Unoccluded(scene)) {
            // Add light's contribution to reflected radiance
            Li *= visibility.Transmittance(scene, renderer, NULL, rng, arena);
            if (light->IsDeltaLight())
                Ld += f * (AbsDot(wi, n) / lightPdf);
            else {
                bsdfPdf = bsdf->Pdf(wo, wi, flags);
                float weight = PowerHeuristic(1, lightPdf, 1, bsdfPdf);
                Ld += f * (AbsDot(wi, n) * weight / lightPdf);
            }
        }
    }

       return Ld;
}

void Light::SHProject(const PbrtPoint &p, float pEpsilon, int lmax,
        const Scene *scene, bool computeLightVisibility, float time,
        RNG &rng, Spectrum *coeffs) const {
    for (int i = 0; i < SHTerms(lmax); ++i)
        coeffs[i] = 0.f;
    uint32_t ns = RoundUpPow2(nSamples);
    uint32_t scramble1D = rng.RandomUInt();
    uint32_t scramble2D[2] = { rng.RandomUInt(), rng.RandomUInt() };
    float *Ylm = ALLOCA(float, SHTerms(lmax));
    for (uint32_t i = 0; i < ns; ++i) {
        // Compute incident radiance sample from _light_, update SH _coeffs_
        float u[2], pdf;
        Sample02(i, scramble2D, u);
        LightSample lightSample(u[0], u[1], VanDerCorput(i, scramble1D));
        Vector wi;
        VisibilityTester vis;
        Spectrum Li = Sample_L(p, pEpsilon, lightSample, time, &wi, &pdf, &vis);
        if (!Li.IsBlack() && pdf > 0.f &&
            (!computeLightVisibility || vis.Unoccluded(scene))) {
            // Add light sample contribution to MC estimate of SH coefficients
            SHEvaluate(wi, lmax, Ylm);
            for (int j = 0; j < SHTerms(lmax); ++j)
                coeffs[j] += Li * Ylm[j] / (pdf * ns);
        }
    }
}

Spectrum EstimateDirect(const Scene *scene, const Renderer *renderer,
        MemoryArena &arena, const Light *light, const Point &p,
        const Normal &n, const Vector &wo, float rayEpsilon, float time,
        const BSDF *bsdf, RNG &rng, const LightSample &lightSample,
        const BSDFSample &bsdfSample, BxDFType flags) {
    Spectrum Ld(0.);
    // Sample light source with multiple importance sampling
    Vector wi;
    float lightPdf, bsdfPdf;
    VisibilityTester visibility;
    Spectrum Li = light->Sample_L(p, rayEpsilon, lightSample, time,
                                  &wi, &lightPdf, &visibility);
    if (lightPdf > 0. && !Li.IsBlack()) {
        Spectrum f = bsdf->f(wo, wi, flags);
        if (!f.IsBlack() && visibility.Unoccluded(scene)) {
            // Add light's contribution to reflected radiance
            Li *= visibility.Transmittance(scene, renderer, NULL, rng, arena);
            if (light->IsDeltaLight())
                Ld += f * Li * (AbsDot(wi, n) / lightPdf);
            else {
                bsdfPdf = bsdf->Pdf(wo, wi, flags);
                float weight = PowerHeuristic(1, lightPdf, 1, bsdfPdf);
                Ld += f * Li * (AbsDot(wi, n) * weight / lightPdf);
            }
        }
    }

    // Sample BSDF with multiple importance sampling
    if (!light->IsDeltaLight()) {
        BxDFType sampledType;
        Spectrum f = bsdf->Sample_f(wo, &wi, bsdfSample, &bsdfPdf, flags,
                                    &sampledType);
        if (!f.IsBlack() && bsdfPdf > 0.) {
            float weight = 1.f;
            if (!(sampledType & BSDF_SPECULAR)) {
                lightPdf = light->Pdf(p, wi);
                if (lightPdf == 0.)
                    return Ld;
                weight = PowerHeuristic(1, bsdfPdf, 1, lightPdf);
            }
            // Add light contribution from BSDF sampling
            Intersection lightIsect;
            Spectrum Li(0.f);
            RayDifferential ray(p, wi, rayEpsilon, INFINITY, time);
            if (scene->Intersect(ray, &lightIsect)) {
                if (lightIsect.primitive->GetAreaLight() == light)
                    Li = lightIsect.Le(-wi);
            }
            else
                //Li = light->Le(ray,);
            	Li = Spectrum(0.);
            if (!Li.IsBlack()) {
                Li *= renderer->Transmittance(scene, ray, NULL, rng, arena);
                Ld += f * Li * AbsDot(wi, n) * weight / bsdfPdf;
            }
        }
    }
    return Ld;
}

LightInfo ProjectionLight::Sample_L(const Point &p, float pEpsilon,
    const LightSample &ls, float time) const {
    auto wi = Normalize(lightPos - p);
    VisibilityTester visibility;
    visibility.SetSegment(p, pEpsilon, lightPos, 0., time);
    auto L = Intensity * Projection(-wi) / DistanceSquared(lightPos, p);
    return LightInfo{L,wi,1.f,visibility};
}

  csColor ProxyLight::SampleLight (const csVector3& point, const csVector3& n,
    float u1, float u2, csVector3& lightVec, float& pdf, VisibilityTester& vistest,
    const csPlane3* visLimitPlane)
  {
    // Setup clipped visibility ray
    const csVector3 lightPos = GetLightSamplePosition (u1, u2);
    csSegment3 visSegment (lightPos, point);
    if (!lightFrustum.Contains (point-lightFrustum.GetOrigin ()) ||
      !lightFrustum.Intersect (visSegment))
    {
      pdf = 0.0f;
      return csColor (0,0,0);
    }

    if (visLimitPlane)
      csIntersect3::SegmentPlane (*visLimitPlane, visSegment);

    vistest.AddSegment (ownerSector->kdTree, visSegment.Start (), visSegment.End ());

    csVector3 parentLightVec;
    csPlane3 transformedPlane;
    transformedPlane = proxyTransform.Other2This (portalPlane);

    const csColor parentLight = parent->SampleLight (point, n, u1, u2, 
      parentLightVec, pdf, vistest, &transformedPlane);
    lightVec = proxyTransform.Other2ThisRelative (parentLightVec);

    return parentLight;
  }

// WhittedIntegrator Method Definitions
Spectrum WhittedIntegrator::Li(const RayDifferential &ray, const Scene &scene,
                               Sampler &sampler, MemoryArena &arena,
                               int depth) const {
    Spectrum L(0.);
    // Find closest ray intersection or return background radiance
    SurfaceInteraction isect;
    if (!scene.Intersect(ray, &isect)) {
        for (const auto &light : scene.lights) L += light->Le(ray);
        return L;
    }

    // Compute emitted and reflected light at ray intersection point

    // Initialize common variables for Whitted integrator
    const Normal3f &n = isect.shading.n;
    Vector3f wo = isect.wo;

    // Compute scattering functions for surface interaction
    isect.ComputeScatteringFunctions(ray, arena);
    if (!isect.bsdf)
        return Li(isect.SpawnRay(ray.d), scene, sampler, arena, depth);

    // Compute emitted light if ray hit an area light source
    L += isect.Le(wo);

    // Add contribution of each light source
    for (const auto &light : scene.lights) {
        Vector3f wi;
        Float pdf;
        VisibilityTester visibility;
        Spectrum Li =
            light->Sample_Li(isect, sampler.Get2D(), &wi, &pdf, &visibility);
        if (Li.IsBlack() || pdf == 0) continue;
        Spectrum f = isect.bsdf->f(wo, wi);
        if (!f.IsBlack() && visibility.Unoccluded(scene))
            L += f * Li * AbsDot(wi, n) / pdf;
    }
    if (depth + 1 < maxDepth) {
        // Trace rays for specular reflection and refraction
        L += SpecularReflect(ray, isect, scene, sampler, arena, depth);
        L += SpecularTransmit(ray, isect, scene, sampler, arena, depth);
    }
    return L;
}

// WhittedIntegrator Method Definitions
Spectrum WhittedIntegrator::Li(const Scene *scene,
        const Renderer *renderer, const RayDifferential &ray,
        const Intersection &isect, const Sample *sample,
        MemoryArena &arena) const {
    Spectrum L(0.);
    // Compute emitted and reflected light at ray intersection point

    // Evaluate BSDF at hit point
    BSDF *bsdf = isect.GetBSDF(ray, arena);

    // 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;
        float pdf;
        Spectrum Li = scene->lights[i]->Sample_L(p, isect.rayEpsilon,
            LightSample(*sample->rng), sample->time, &wi, &pdf, &visibility);
        if (Li.IsBlack() || pdf == 0.f) continue;
        Li /= pdf;
        Spectrum f = bsdf->f(wo, wi);
        if (!f.IsBlack() && visibility.Unoccluded(scene))
            L += f * Li * AbsDot(wi, n) *
                visibility.Transmittance(scene, renderer,
                                         sample, NULL, arena);
    }
    if (ray.depth + 1 < maxDepth) {
        // Trace rays for specular reflection and refraction
        L += SpecularReflect(ray, bsdf, *sample->rng, isect, renderer,
                             scene, sample, arena);
        L += SpecularTransmit(ray, bsdf, *sample->rng, isect, renderer,
                              scene, sample, arena);
    }
    return L;
}

Spectrum VolumePatIntegrator::EstimateDirectLight(const Scene *scene,
        const Renderer *renderer, MemoryArena &arena, const Light *light,
        const Point &p, const Normal &n, const Vector &wo, float rayEpsilon,
        float time, RNG &rng, const LightSample &lightSample) const {
    VolumeRegion *vr = scene->volumeRegion;
    if (!vr) Spectrum(0.);
    Spectrum Ld(0.);
    // Sample light source.
    Vector wi;
    float lightPdf;
    VisibilityTester visibility;
    Spectrum Li = light->Sample_L(p, rayEpsilon, lightSample, time,
                                  &wi, &lightPdf, &visibility);
    if (lightPdf > 0. && !Li.IsBlack() && visibility.Unoccluded(scene)) {
        // Add light's contribution to reflected radiance
        Li *= visibility.Transmittance(scene, renderer, NULL, rng, arena);
        // Li *= PowerHeuristic(1, lightPdf, 1, (1/M_PI_4));
        Ld += vr->p(p, -wi, wo, time) * vr->Sigma_s(p, wo, time) * Li / lightPdf;
    }
    return Ld;
}

Spectrum SingleScatteringFluorescenceRWLIntegrator::Li(const Scene *scene,
        const Renderer *renderer, const RayDifferential &ray,
        const Sample *sample, RNG &rng, Spectrum *T, MemoryArena &arena) const {
    VolumeRegion *vr = scene->volumeRegion;
    float t0, t1;
    if (!vr || !vr->IntersectP(ray, &t0, &t1) || (t1-t0) == 0.f) {
        *T = 1.f;
        return 0.f;
    }
    // Do single scattering volume integration in _vr_
    Spectrum Lv(0.);

    // Prepare for volume integration stepping
    int nSamples = Ceil2Int((t1-t0) / stepSize);
    float step = (t1 - t0) / nSamples;
    Spectrum Tr(1.f);
    Point p = ray(t0), pPrev;
    Vector w = -ray.d;
    t0 += sample->oneD[scatterSampleOffset][0] * step;

    // Compute sample patterns for single scattering samples
    float *lightNum = arena.Alloc<float>(nSamples);
    LDShuffleScrambled1D(1, nSamples, lightNum, rng);
    float *lightComp = arena.Alloc<float>(nSamples);
    LDShuffleScrambled1D(1, nSamples, lightComp, rng);
    float *lightPos = arena.Alloc<float>(2*nSamples);
    LDShuffleScrambled2D(1, nSamples, lightPos, rng);
    uint32_t sampOffset = 0;
    for (int i = 0; i < nSamples; ++i, t0 += step) {
        // Advance to sample at _t0_ and update _T_
        pPrev = p;
        p = ray(t0);

        Ray tauRay(pPrev, p - pPrev, 0.f, 1.f, ray.time, ray.depth);
        Spectrum stepTau = vr->tau(tauRay, 0.5f * stepSize, rng.RandomFloat());
        Tr *= Exp(-stepTau);

        // Possibly terminate ray marching if transmittance is small
        if (Tr.y() < 1e-3) {
            const float continueProb = .5f;
            if (rng.RandomFloat() > continueProb) {
                Tr = 0.f;
                break;
            }
            Tr /= continueProb;
        }

        // Compute fluorescence emission
        Spectrum sigma = vr->Mu(p, w, ray.time);
        if (!sigma.IsBlack() && scene->lights.size() > 0) {
            int nLights = scene->lights.size();
            int ln = min(Floor2Int(lightNum[sampOffset] * nLights),
                         nLights-1);
            Light *light = scene->lights[ln];

            // Add contribution of _light_ due to the in-scattering at _p_
            float pdf;
            VisibilityTester vis;
            Vector wo;
            LightSample ls(lightComp[sampOffset], lightPos[2*sampOffset],
                           lightPos[2*sampOffset+1]);
            Spectrum L = light->Sample_L(p, 0.f, ls, ray.time, &wo, &pdf, &vis);
            if (!L.IsBlack() && pdf > 0.f && vis.Unoccluded(scene)) {
                Spectrum Ld = L * vis.Transmittance(scene, renderer, NULL, rng,
                                                    arena);
                int lambdaExcIndex = light->GetLaserWavelengthIndex();
                float Lpower = Ld.GetLaserEmissionPower(lambdaExcIndex);
                float yield = vr->Yeild(Point());
                Spectrum fEx = vr->fEx(Point());
                Spectrum fEm = vr->fEm(Point());
                float scale = fEx.GetSampleValueAtWavelengthIndex(lambdaExcIndex);
                Lv += Lpower * Tr * sigma * vr->p(p, w, -wo, ray.time) *
                        scale * fEm * yield * float(nLights) / pdf;
            }
        }
        ++sampOffset;
    }
    *T = Tr;
    return Lv * step;
}

Spectrum EstimateDirect(const Scene *scene,
        const Light *light, const Point &p,
		const Normal &n, const Vector &wo,
		BSDF *bsdf, const Sample *sample, int lightSamp,
		int bsdfSamp, int bsdfComponent, u_int sampleNum) {
  Spectrum Ld = Spectrum(0.f);

	// Find light and BSDF sample values for direct lighting estimate
	float ls1, ls2, bs1, bs2, bcs;
	if (lightSamp != -1 && bsdfSamp != -1 &&
		sampleNum < sample->n2D[lightSamp] &&
		sampleNum < sample->n2D[bsdfSamp]) {
		ls1 = sample->twoD[lightSamp][2*sampleNum];
		ls2 = sample->twoD[lightSamp][2*sampleNum+1];
		bs1 = sample->twoD[bsdfSamp][2*sampleNum];
		bs2 = sample->twoD[bsdfSamp][2*sampleNum+1];
		bcs = sample->oneD[bsdfComponent][sampleNum];
	}
	else {
		ls1 = RandomFloat();
		ls2 = RandomFloat();
		bs1 = RandomFloat();
		bs2 = RandomFloat();
		bcs = RandomFloat();
	}
	// Sample light source with multiple importance sampling
	Vector wi;
	float lightPdf, bsdfPdf;
	VisibilityTester visibility;
	Spectrum Li = light->Sample_L(p, n, ls1, ls2, &wi, &lightPdf, &visibility);
	if (lightPdf > 0. && !Li.Black()) {
		
		if(bsdf->NumComponents(BSDF_FLUORESCENT) > 0)
		{
			Bispectrum * fluoro = (Bispectrum*)bsdf->f_ptr(wo, wi, BxDFType(BSDF_FLUORESCENT));
			if (visibility.Unoccluded(scene)) {
				// Add light's contribution to reflected radiance
				Li *= visibility.Transmittance(scene);
				if (light->IsDeltaLight())
					Ld += fluoro->output(Li, true, false)  * AbsDot(wi, n);
				else {
					bsdfPdf = bsdf->Pdf(wo, wi);
					//float weight = PowerHeuristic(1, lightPdf, 1, bsdfPdf);
					Ld += fluoro->output(Li, true, false) * AbsDot(wi, n)*bsdfPdf;
				}
			}
		}else{
		
			Spectrum f = bsdf->f(wo, wi, BxDFType(BSDF_ALL & ~BSDF_FLUORESCENT));
			if (!f.Black() && visibility.Unoccluded(scene)) {
				// Add light's contribution to reflected radiance
				Li *= visibility.Transmittance(scene);
				if (light->IsDeltaLight())
					Ld += f * Li * AbsDot(wi, n);
				else {
					bsdfPdf = bsdf->Pdf(wo, wi);
					//float weight = PowerHeuristic(1, lightPdf, 1, bsdfPdf);
					Ld += f * Li * AbsDot(wi, n)*bsdfPdf;
				}
			}
		}
		
		
	}/*
	// Sample BSDF with multiple importance sampling
	if (!light->IsDeltaLight()) {
		BxDFType flags = BxDFType(BSDF_ALL & ~BSDF_SPECULAR & ~BSDF_FLUORESCENT);
		if(bsdf->NumComponents(BSDF_FLUORESCENT) > 0)
		{
			Bispectrum * fluoro = (Bispectrum*)bsdf->Sample_f_ptr(wo, &wi, bs1, bs2, bcs, &bsdfPdf, BxDFType(BSDF_FLUORESCENT));
			if (bsdfPdf > 0.) {
				lightPdf = light->Pdf(p, n, wi);
				if (lightPdf > 0.) {
					// Add light contribution from BSDF sampling
					float weight = PowerHeuristic(1, bsdfPdf, 1, lightPdf);
					Intersection lightIsect;
					Spectrum Li(0.f);
					RayDifferential ray(p, wi);
					if (scene->Intersect(ray, &lightIsect)) {
						if (lightIsect.primitive->GetAreaLight() == light)
							Li = lightIsect.Le(-wi);
					}
					else
						Li = light->Le(ray);
					if (!Li.Black()) {
						Li *= scene->Transmittance(ray);
						Ld += fluoro->output(Li, true, false) * AbsDot(wi, n)* weight / bsdfPdf;
					}
				}
			}
		}else{
		
			Spectrum f = bsdf->Sample_f(wo, &wi, bs1, bs2, bcs, &bsdfPdf, flags);
			if (!f.Black() && bsdfPdf > 0.) {
				lightPdf = light->Pdf(p, n, wi);
				if (lightPdf > 0.) {
					// Add light contribution from BSDF sampling
					float weight = PowerHeuristic(1, bsdfPdf, 1, lightPdf);
					Intersection lightIsect;
					Spectrum Li(0.f);
//.........这里部分代码省略.........

Spectrum SingleScatteringIntegrator::Li(const Scene *scene, const Renderer *renderer,
        const RayDifferential &ray, const Sample *sample,
        Spectrum *T, MemoryArena &arena) const {
    VolumeRegion *vr = scene->volumeRegion;
    float t0, t1;
    if (!vr || !vr->IntersectP(ray, &t0, &t1)) {
        *T = 1.f;
        return 0.f;
    }
    // Do single scattering volume integration in _vr_
    Spectrum Lv(0.);

    // Prepare for volume integration stepping
    int nSamples = Ceil2Int((t1-t0) / stepSize);
    float step = (t1 - t0) / nSamples;
    Spectrum Tr(1.f);
    Point p = ray(t0), pPrev;
    Vector w = -ray.d;
    t0 += sample->oneD[scatterSampleOffset][0] * step;

    // Compute sample patterns for single scattering samples
    float *lightNum = arena.Alloc<float>(nSamples);
    LDShuffleScrambled1D(1, nSamples, lightNum, *sample->rng);
    float *lightComp = arena.Alloc<float>(nSamples);
    LDShuffleScrambled1D(1, nSamples, lightComp, *sample->rng);
    float *lightPos = arena.Alloc<float>(2*nSamples);
    LDShuffleScrambled2D(1, nSamples, lightPos, *sample->rng);
    u_int sampOffset = 0;
    for (int i = 0; i < nSamples; ++i, t0 += step) {
        // Advance to sample at _t0_ and update _T_
        pPrev = p;
        p = ray(t0);
        Ray tauRay(pPrev, p - pPrev, 0.f, 1.f, ray.time, ray.depth);
        Spectrum stepTau = vr->tau(tauRay,
                                   .5f * stepSize, sample->rng->RandomFloat());
        Tr *= Exp(-stepTau);

        // Possibly terminate ray marching if transmittance is small
        if (Tr.y() < 1e-3) {
            const float continueProb = .5f;
            if (sample->rng->RandomFloat() > continueProb) break;
            Tr /= continueProb;
        }

        // Compute single-scattering source term at _p_
        Lv += Tr * vr->Lve(p, w, ray.time);
        Spectrum ss = vr->sigma_s(p, w, ray.time);
        if (!ss.IsBlack() && scene->lights.size() > 0) {
            int nLights = scene->lights.size();
            int ln = min(Floor2Int(lightNum[sampOffset] * nLights),
                         nLights-1);
            Light *light = scene->lights[ln];
            // Add contribution of _light_ due to scattering at _p_
            float pdf;
            VisibilityTester vis;
            Vector wo;
            LightSample ls(lightComp[sampOffset], lightPos[2*sampOffset],
                           lightPos[2*sampOffset+1]);
            Spectrum L = light->Sample_L(p, 0.f, ls, ray.time, &wo, &pdf, &vis);
            if (!L.IsBlack() && pdf > 0.f && vis.Unoccluded(scene)) {
                Spectrum Ld = L * vis.Transmittance(scene, renderer, NULL, sample->rng, arena);
                Lv += Tr * ss * vr->p(p, w, -wo, ray.time) * Ld * float(nLights) / pdf;
            }
        }
        ++sampOffset;
    }
    *T = Tr;
    return Lv * step;
}

Spectrum PhotonVolumeIntegrator::Li(const Scene *scene, const Renderer *renderer,
        const RayDifferential &ray, const Sample *sample, RNG &rng,
        Spectrum *T, MemoryArena &arena) const {
 	
 	VolumeRegion *vr = scene->volumeRegion;
    RainbowVolume* rv = dynamic_cast<RainbowVolume*>(vr);
 	KdTree<Photon>* volumeMap = photonShooter->volumeMap; 

 	float t0, t1;
 	if (!vr || !vr->IntersectP(ray, &t0, &t1) || (t1-t0) == 0.f){
 		*T = 1.f;
 	 	return 0.f;
 	 }
 	// Do single scattering & photon multiple scattering volume integration in _vr_
 	Spectrum Lv(0.);


 	// Prepare for volume integration stepping
 	int nSamples = Ceil2Int((t1-t0) / stepSize);
 	float step = (t1 - t0) / nSamples;
 	Spectrum Tr(1.f);
 	Point p = ray(t0), pPrev;
 	Vector w = -ray.d;
 	t0 += sample->oneD[scatterSampleOffset][0] * step;

 	float *lightNum = arena.Alloc<float>(nSamples);
    LDShuffleScrambled1D(1, nSamples, lightNum, rng);
    float *lightComp = arena.Alloc<float>(nSamples);
    LDShuffleScrambled1D(1, nSamples, lightComp, rng);
    float *lightPos = arena.Alloc<float>(2*nSamples);
    LDShuffleScrambled2D(1, nSamples, lightPos, rng);
 	int sampOffset = 0;

 	ClosePhoton *lookupBuf = new ClosePhoton[nSamples];

 	for (int i = 0; i < nSamples; ++i, t0 += step) {
 		// Advance to sample at _t0_ and update _T_
 		pPrev = p;
 		p = ray(t0);

 		Ray tauRay(pPrev, p - pPrev, 0.f, 1.f, ray.time, ray.depth);

 		Spectrum stepTau = vr->tau(tauRay,.5f * stepSize, rng.RandomFloat());
 		Tr = Exp(-stepTau);

 		// Possibly terminate raymarching if transmittance is small.
 		if (Tr.y() < 1e-3) {
 			const float continueProb = .5f;
 			if (rng.RandomFloat() > continueProb){
 				Tr = 0.f;
 				break;
 			}
 			Tr /= continueProb;
 		}
		
		
 		// Compute single-scattering source term at _p_ & photon mapped MS
 		Spectrum L_i(0.);
 		Spectrum L_d(0.);
 		Spectrum L_ii(0.);
 		
 		// Lv += Tr*vr->Lve(p, w, ray.time);
 		Spectrum ss = vr->sigma_s(p, w, ray.time);
 		Spectrum sa = vr->sigma_a(p, w, ray.time);

 		if (!ss.IsBlack() && scene->lights.size() > 0) {
 			int nLights = scene->lights.size();
 			int ln =
 				min(Floor2Int(lightNum[sampOffset] * nLights),
 				    nLights-1);
 			Light *light = scene->lights[ln];
 			// Add contribution of _light_ due to scattering at _p_
 			float pdf;
 			VisibilityTester vis;
 			Vector wo;

 			LightSample ls(lightComp[sampOffset], lightPos[2*sampOffset],
                           lightPos[2*sampOffset+1]);
            Spectrum L = light->Sample_L(p, 0.f, ls, ray.time, &wo, &pdf, &vis);
            

 			if (!L.IsBlack() && pdf > 0.f && vis.Unoccluded(scene)) {

                Spectrum Ld = L * vis.Transmittance(scene,renderer, NULL, rng, arena);
                if(rv){
                    L_d = rv->rainbowReflection(Ld, ray.d, wo);
                }
                else {
                    L_d = vr->p(p, w, -wo, ray.time) * Ld * float(nLights)/pdf;
                }
 			}
 		}
		// Compute 'indirect' in-scattered radiance from photon map
        if(!rv){
            L_ii += LPhoton(volumeMap, nUsed, lookupBuf, w, p, vr, maxDistSquared, ray.time);
        }
        
		// Compute total in-scattered radiance
		if (sa.y()!=0.0 || ss.y()!=0.0)
			L_i = L_d + (ss/(sa+ss))*L_ii;
//.........这里部分代码省略.........

void DipoleSubsurfaceIntegrator::Preprocess(const Scene *scene,
        const Camera *camera, const Renderer *renderer) {
    if (scene->lights.size() == 0) return;
    vector<SurfacePoint> pts;
    // Get _SurfacePoint_s for translucent objects in scene
    if (filename != "") {
        // Initialize _SurfacePoint_s from file
        vector<float> fpts;
        if (ReadFloatFile(filename.c_str(), &fpts)) {
            if ((fpts.size() % 8) != 0)
                Error("Excess values (%d) in points file \"%s\"", int(fpts.size() % 8),
                      filename.c_str());
            for (u_int i = 0; i < fpts.size(); i += 8)
                pts.push_back(SurfacePoint(Point(fpts[i], fpts[i+1], fpts[i+2]),
                                           Normal(fpts[i+3], fpts[i+4], fpts[i+5]),
                                           fpts[i+6], fpts[i+7]));
        }
    }
    if (pts.size() == 0) {
        Point pCamera = camera->CameraToWorld(camera->shutterOpen,
                                              Point(0, 0, 0));
        FindPoissonPointDistribution(pCamera, camera->shutterOpen,
                                     minSampleDist, scene, &pts);
    }

    // Compute irradiance values at sample points
    RNG rng;
    MemoryArena arena;
    PBRT_SUBSURFACE_STARTED_COMPUTING_IRRADIANCE_VALUES();
    ProgressReporter progress(pts.size(), "Computing Irradiances");
    for (uint32_t i = 0; i < pts.size(); ++i) {
        SurfacePoint &sp = pts[i];
        Spectrum E(0.f);
        for (uint32_t j = 0; j < scene->lights.size(); ++j) {
            // Add irradiance from light at point
            const Light *light = scene->lights[j];
            Spectrum Elight = 0.f;
            int nSamples = RoundUpPow2(light->nSamples);
            uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
            uint32_t compScramble = rng.RandomUInt();
            for (int s = 0; s < nSamples; ++s) {
                float lpos[2];
                Sample02(s, scramble, lpos);
                float lcomp = VanDerCorput(s, compScramble);
                LightSample ls(lpos[0], lpos[1], lcomp);
                Vector wi;
                float lightPdf;
                VisibilityTester visibility;
                Spectrum Li = light->Sample_L(sp.p, sp.rayEpsilon,
                    ls, camera->shutterOpen, &wi, &lightPdf, &visibility);
                if (Dot(wi, sp.n) <= 0.) continue;
                if (Li.IsBlack() || lightPdf == 0.f) continue;
                Li *= visibility.Transmittance(scene, renderer, NULL, rng, arena);
                if (visibility.Unoccluded(scene))
                    Elight += Li * AbsDot(wi, sp.n) / lightPdf;
            }
            E += Elight / nSamples;
        }
        if (E.y() > 0.f)
        {
            irradiancePoints.push_back(IrradiancePoint(sp, E));
            PBRT_SUBSURFACE_COMPUTED_IRRADIANCE_AT_POINT(&sp, &E);
        }
        arena.FreeAll();
        progress.Update();
    }
    progress.Done();
    PBRT_SUBSURFACE_FINISHED_COMPUTING_IRRADIANCE_VALUES();

    // Create octree of clustered irradiance samples
    octree = octreeArena.Alloc<SubsurfaceOctreeNode>();
    for (uint32_t i = 0; i < irradiancePoints.size(); ++i)
        octreeBounds = Union(octreeBounds, irradiancePoints[i].p);
    for (uint32_t i = 0; i < irradiancePoints.size(); ++i)
        octree->Insert(octreeBounds, &irradiancePoints[i], octreeArena);
    octree->InitHierarchy();
}

Spectrum EstimateDirect(const Interaction &it, const Point2f &uScattering,
                        const Light &light, const Point2f &uLight,
                        const Scene &scene, Sampler &sampler,
                        MemoryArena &arena, bool handleMedia, bool specular) {
    BxDFType bsdfFlags =
        specular ? BSDF_ALL : BxDFType(BSDF_ALL & ~BSDF_SPECULAR);
    Spectrum Ld(0.f);
    // Sample light source with multiple importance sampling
    Vector3f wi;
    Float lightPdf = 0, scatteringPdf = 0;
    VisibilityTester visibility;
    Spectrum Li = light.Sample_Li(it, uLight, &wi, &lightPdf, &visibility);
    if (lightPdf > 0 && !Li.IsBlack()) {
        // Compute BSDF or phase function's value for light sample
        Spectrum f;
        if (it.IsSurfaceInteraction()) {
            // Evaluate BSDF for light sampling strategy
            const SurfaceInteraction &isect = (const SurfaceInteraction &)it;
            f = isect.bsdf->f(isect.wo, wi, bsdfFlags) *
                AbsDot(wi, isect.shading.n);
            scatteringPdf = isect.bsdf->Pdf(isect.wo, wi, bsdfFlags);
        } else {
            // Evaluate phase function for light sampling strategy
            const MediumInteraction &mi = (const MediumInteraction &)it;
            Float p = mi.phase->p(mi.wo, wi);
            f = Spectrum(p);
            scatteringPdf = p;
        }
        if (!f.IsBlack()) {
            // Compute effect of visibility for light source sample
            if (handleMedia)
                Li *= visibility.Tr(scene, sampler);
            else if (!visibility.Unoccluded(scene))
                Li = Spectrum(0.f);

            // Add light's contribution to reflected radiance
            if (!Li.IsBlack()) {
                if (IsDeltaLight(light.flags))
                    Ld += f * Li / lightPdf;
                else {
                    Float weight =
                        PowerHeuristic(1, lightPdf, 1, scatteringPdf);
                    Ld += f * Li * weight / lightPdf;
                }
            }
        }
    }

    // Sample BSDF with multiple importance sampling
    if (!IsDeltaLight(light.flags)) {
        Spectrum f;
        bool sampledSpecular = false;
        if (it.IsSurfaceInteraction()) {
            // Sample scattered direction for surface interactions
            BxDFType sampledType;
            const SurfaceInteraction &isect = (const SurfaceInteraction &)it;
            f = isect.bsdf->Sample_f(isect.wo, &wi, uScattering, &scatteringPdf,
                                     bsdfFlags, &sampledType);
            f *= AbsDot(wi, isect.shading.n);
            sampledSpecular = sampledType & BSDF_SPECULAR;
        } else {
            // Sample scattered direction for medium interactions
            const MediumInteraction &mi = (const MediumInteraction &)it;
            Float p = mi.phase->Sample_p(mi.wo, &wi, uScattering);
            f = Spectrum(p);
            scatteringPdf = p;
        }
        if (!f.IsBlack() && scatteringPdf > 0) {
            // Account for light contributions along sampled direction _wi_
            Float weight = 1;
            if (!sampledSpecular) {
                lightPdf = light.Pdf_Li(it, wi);
                if (lightPdf == 0) return Ld;
                weight = PowerHeuristic(1, scatteringPdf, 1, lightPdf);
            }

            // Find intersection and compute transmittance
            SurfaceInteraction lightIsect;
            Ray ray = it.SpawnRay(wi);
            Spectrum Tr(1.f);
            bool foundSurfaceInteraction =
                handleMedia ? scene.IntersectTr(ray, sampler, &lightIsect, &Tr)
                            : scene.Intersect(ray, &lightIsect);

            // Add light contribution from material sampling
            Spectrum Li(0.f);
            if (foundSurfaceInteraction) {
                if (lightIsect.primitive->GetAreaLight() == &light)
                    Li = lightIsect.Le(-wi);
            } else
                Li = light.Le(ray);
            if (!Li.IsBlack()) Ld += f * Li * Tr * weight / scatteringPdf;
        }
    }
    return Ld;
}