自学教程：C++ BxDFType函数代码示例

Spectrum FresnelSpecular::Sample_f(const Vector3f &wo, Vector3f *wi,                                   const Point2f &u, Float *pdf,                                   BxDFType *sampledType) const {    Float F = FrDielectric(CosTheta(wo), etaA, etaB);    if (u[0] < F) {        // Compute specular reflection for _FresnelSpecular_        // Compute perfect specular reflection direction        *wi = Vector3f(-wo.x, -wo.y, wo.z);        if (sampledType)            *sampledType = BxDFType(BSDF_SPECULAR | BSDF_REFLECTION);        *pdf = F;        return F * R / AbsCosTheta(*wi);    } else {        // Compute specular transmission for _FresnelSpecular_        // Figure out which $/eta$ is incident and which is transmitted        bool entering = CosTheta(wo) > 0;        Float etaI = entering ? etaA : etaB;        Float etaT = entering ? etaB : etaA;        // Compute ray direction for specular transmission        if (!Refract(wo, Faceforward(Normal3f(0, 0, 1), wo), etaI / etaT, wi))            return 0;        Spectrum ft = T * (1 - F);        // Account for non-symmetry with transmission to different medium        if (mode == TransportMode::Radiance)            ft *= (etaI * etaI) / (etaT * etaT);        if (sampledType)            *sampledType = BxDFType(BSDF_SPECULAR | BSDF_TRANSMISSION);        *pdf = 1 - F;        return ft / AbsCosTheta(*wi);    }}

Spectrum ExPhotonIntegrator::LPhoton(		KdTree<Photon, PhotonProcess> *map,		int nPaths, int nLookup, BSDF *bsdf,		const Intersection &isect, const Vector &wo,		float maxDistSquared) {	Spectrum L(0.);	if (!map) return L;	BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |		BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);	if (bsdf->NumComponents(nonSpecular) == 0)		return L;	static StatsCounter lookups("Photon Map", "Total lookups"); // NOBOOK	// Initialize _PhotonProcess_ object, _proc_, for photon map lookups	PhotonProcess proc(nLookup, isect.dg.p);	proc.photons =		(ClosePhoton *)alloca(nLookup * sizeof(ClosePhoton));	// Do photon map lookup	++lookups;  // NOBOOK	map->Lookup(isect.dg.p, proc, maxDistSquared);	// Accumulate light from nearby photons	static StatsRatio foundRate("Photon Map", "Photons found per lookup"); // NOBOOK	foundRate.Add(proc.foundPhotons, 1); // NOBOOK	// Estimate reflected light from photons	ClosePhoton *photons = proc.photons;	int nFound = proc.foundPhotons;	Normal Nf = Dot(wo, bsdf->dgShading.nn) < 0 ? -bsdf->dgShading.nn :		bsdf->dgShading.nn;	if (bsdf->NumComponents(BxDFType(BSDF_REFLECTION |			BSDF_TRANSMISSION | BSDF_GLOSSY)) > 0) {		// Compute exitant radiance from photons for glossy surface		for (int i = 0; i < nFound; ++i) {			const Photon *p = photons[i].photon;			BxDFType flag = Dot(Nf, p->wi) > 0.f ?				BSDF_ALL_REFLECTION : BSDF_ALL_TRANSMISSION;			float k = kernel(p, isect.dg.p, maxDistSquared);			L += (k / nPaths) * bsdf->f(wo, p->wi, flag) * p->alpha;		}	}	else {		// Compute exitant radiance from photons for diffuse surface		Spectrum Lr(0.), Lt(0.);		for (int i = 0; i < nFound; ++i) {			if (Dot(Nf, photons[i].photon->wi) > 0.f) {				float k = kernel(photons[i].photon, isect.dg.p,					maxDistSquared);				Lr += (k / nPaths) * photons[i].photon->alpha;			}			else {				float k = kernel(photons[i].photon, isect.dg.p,					maxDistSquared);				Lt += (k / nPaths) * photons[i].photon->alpha;			}		}		L += Lr * bsdf->rho(wo, BSDF_ALL_REFLECTION) * INV_PI +			Lt * bsdf->rho(wo, BSDF_ALL_TRANSMISSION) * INV_PI;	}	return L;}

Spectrum BSDF::Sample_f(const Vector3f &woWorld, Vector3f *wiWorld,                        const Point2f &u, Float *pdf, BxDFType type,                        BxDFType *sampledType) const {    ProfilePhase pp(Prof::BSDFEvaluation);    // Choose which _BxDF_ to sample    int matchingComps = NumComponents(type);    if (matchingComps == 0) {        *pdf = 0;        if (sampledType) *sampledType = BxDFType(0);        return Spectrum(0);    }    int comp =        std::min((int)std::floor(u[0] * matchingComps), matchingComps - 1);    // Get _BxDF_ pointer for chosen component    BxDF *bxdf = nullptr;    int count = comp;    for (int i = 0; i < nBxDFs; ++i)        if (bxdfs[i]->MatchesFlags(type) && count-- == 0) {            bxdf = bxdfs[i];            break;        }    Assert(bxdf);    // Remap _BxDF_ sample _u_ to $[0,1)^2$    Point2f uRemapped(u[0] * matchingComps - comp, u[1]);    // Sample chosen _BxDF_    Vector3f wi, wo = WorldToLocal(woWorld);    *pdf = 0;    if (sampledType) *sampledType = bxdf->type;    Spectrum f = bxdf->Sample_f(wo, &wi, uRemapped, pdf, sampledType);    if (*pdf == 0) {        if (sampledType) *sampledType = BxDFType(0);        return 0;    }    *wiWorld = LocalToWorld(wi);    // Compute overall PDF with all matching _BxDF_s    if (!(bxdf->type & BSDF_SPECULAR) && matchingComps > 1)        for (int i = 0; i < nBxDFs; ++i)            if (bxdfs[i] != bxdf && bxdfs[i]->MatchesFlags(type))                *pdf += bxdfs[i]->Pdf(wo, wi);    if (matchingComps > 1) *pdf /= matchingComps;    // Compute value of BSDF for sampled direction    if (!(bxdf->type & BSDF_SPECULAR) && matchingComps > 1) {        bool reflect = Dot(*wiWorld, ng) * Dot(woWorld, ng) > 0;        f = 0.;        for (int i = 0; i < nBxDFs; ++i)            if (bxdfs[i]->MatchesFlags(type) &&                ((reflect && (bxdfs[i]->type & BSDF_REFLECTION)) ||                 (!reflect && (bxdfs[i]->type & BSDF_TRANSMISSION))))                f += bxdfs[i]->f(wo, wi);    }    return f;}

Spectrum LPhoton(KdTree<Photon> *map, int nPaths, int nLookup,      MemoryArena &arena, BSDF *bsdf, RNG &rng, const Intersection &isect,      const Vector &wo, float maxDistSquared) {    Spectrum L(0.);    BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |        BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);    if (map && bsdf->NumComponents(nonSpecular) > 0) {        PBRT_PHOTON_MAP_STARTED_LOOKUP(const_cast<DifferentialGeometry *>(&isect.dg));        // Do photon map lookup at intersection point        PhotonProcess proc(nLookup, arena.Alloc<ClosePhoton>(nLookup));        map->Lookup(isect.dg.p, proc, maxDistSquared);        // Estimate reflected radiance due to incident photons        ClosePhoton *photons = proc.photons;        int nFound = proc.nFound;        Normal Nf = Faceforward(bsdf->dgShading.nn, wo);        if (bsdf->NumComponents(BxDFType(BSDF_REFLECTION |                BSDF_TRANSMISSION | BSDF_GLOSSY)) > 0) {            // Compute exitant radiance from photons for glossy surface            for (int i = 0; i < nFound; ++i) {                const Photon *p = photons[i].photon;                float k = kernel(p, isect.dg.p, maxDistSquared);                L += (k / (nPaths * maxDistSquared)) * bsdf->f(wo, p->wi) *                     p->alpha;            }        }        else {            // Compute exitant radiance from photons for diffuse surface            Spectrum Lr(0.), Lt(0.);            for (int i = 0; i < nFound; ++i) {                if (Dot(Nf, photons[i].photon->wi) > 0.f) {                    float k = kernel(photons[i].photon, isect.dg.p,                        maxDistSquared);                    Lr += (k / (nPaths * maxDistSquared)) * photons[i].photon->alpha;                }                else {                    float k = kernel(photons[i].photon, isect.dg.p,                        maxDistSquared);                    Lt += (k / (nPaths * maxDistSquared)) * photons[i].photon->alpha;                }            }            const int sqrtRhoSamples = 4;            float rhoRSamples[2*sqrtRhoSamples*sqrtRhoSamples];            StratifiedSample2D(rhoRSamples, sqrtRhoSamples, sqrtRhoSamples, rng);            float rhoTSamples[2*sqrtRhoSamples*sqrtRhoSamples];            StratifiedSample2D(rhoTSamples, sqrtRhoSamples, sqrtRhoSamples, rng);            L += Lr * bsdf->rho(wo, sqrtRhoSamples*sqrtRhoSamples, rhoRSamples,                                BSDF_ALL_REFLECTION) * INV_PI +                 Lt * bsdf->rho(wo, sqrtRhoSamples*sqrtRhoSamples, rhoTSamples,                                BSDF_ALL_TRANSMISSION) * INV_PI;        }        PBRT_PHOTON_MAP_FINISHED_LOOKUP(const_cast<DifferentialGeometry *>(&isect.dg),            proc.nFound, proc.nLookup, &L);    }    return L;}

Spectrum BSDF::Sample_f(const Vector &woW, Vector *wiW,		float u1, float u2, float u3, float *pdf,		BxDFType flags, BxDFType *sampledType) const {	// Choose which _BxDF_ to sample	int matchingComps = NumComponents(flags);	if (matchingComps == 0) {		*pdf = 0.f;		return Spectrum(0.f);	}	int which = min(Floor2Int(u3 * matchingComps),		matchingComps-1);	BxDF *bxdf = NULL;	int count = which;	for (int i = 0; i < nBxDFs; ++i)		if (bxdfs[i]->MatchesFlags(flags))			if (count-- == 0) {				bxdf = bxdfs[i];				break;			}	Assert(bxdf); // NOBOOK	// Sample chosen _BxDF_	Vector wi;	Vector wo = WorldToLocal(woW);	*pdf = 0.f;	Spectrum f = bxdf->Sample_f(wo, &wi, u1, u2, pdf);	if (*pdf == 0.f) return 0.f;	if (sampledType) *sampledType = bxdf->type;	*wiW = LocalToWorld(wi);	// Compute overall PDF with all matching _BxDF_s	if (!(bxdf->type & BSDF_SPECULAR) && matchingComps > 1) {		for (int i = 0; i < nBxDFs; ++i) {			if (bxdfs[i] != bxdf &&			    bxdfs[i]->MatchesFlags(flags))				*pdf += bxdfs[i]->Pdf(wo, wi);		}	}	if (matchingComps > 1) *pdf /= matchingComps;	// Compute value of BSDF for sampled direction	if (!(bxdf->type & BSDF_SPECULAR)) {		f = 0.;		if (Dot(*wiW, ng) * Dot(woW, ng) > 0)			// ignore BTDFs			flags = BxDFType(flags & ~BSDF_TRANSMISSION);		else			// ignore BRDFs			flags = BxDFType(flags & ~BSDF_REFLECTION);		for (int i = 0; i < nBxDFs; ++i)			if (bxdfs[i]->MatchesFlags(flags))				f += bxdfs[i]->f(wo, wi);	}	return f;}

Spectrum BSDF::f(const Vector &woW,		const Vector &wiW, BxDFType flags) const {	Vector wi = WorldToLocal(wiW), wo = WorldToLocal(woW);	if (Dot(wiW, ng) * Dot(woW, ng) > 0)		// ignore BTDFs		flags = BxDFType(flags & ~BSDF_TRANSMISSION);	else		// ignore BRDFs		flags = BxDFType(flags & ~BSDF_REFLECTION);	Spectrum f = 0.;	for (int i = 0; i < nBxDFs; ++i)		if (bxdfs[i]->MatchesFlags(flags))			f += bxdfs[i]->f(wo, wi);	return f;}

// Skin Method DefinitionsBSDF *Skin::GetBSDF(const DifferentialGeometry &dgGeom, const DifferentialGeometry &dgShading) const {	// Declare skin coefficients	static float diffuse[3] = {  0.428425f,  0.301341f,  0.331054f};	static float xy0[3] =     { -1.131747f, -1.016939f, -0.966018f};	static float z0[3] =      { -1.209182f, -1.462488f, -1.222419f};	static float e0[3] =      {  6.421658f,  3.699932f,  3.524889f};	static float xy1[3] =     { -0.546570f, -0.643533f, -0.638934f};	static float z1[3] =      {  0.380123f,  0.410559f,  0.437367f};	static float e1[3] =      {  3.685044f,  4.266495f,  4.539742f};	static float xy2[3] =     { -0.998888f, -1.020153f, -1.027479f};	static float z2[3] =      {  0.857998f,  0.703913f,  0.573625f};	static float e2[3] =      { 64.208486f, 63.919687f, 43.809866f};	static Spectrum xy[3] = { Spectrum(xy0), Spectrum(xy1), Spectrum(xy2) };	static Spectrum z[3] = { Spectrum(z0), Spectrum(z1), Spectrum(z2) };	static Spectrum e[3] = { Spectrum(e0), Spectrum(e1), Spectrum(e2) };	// Allocate _BSDF_, possibly doing bump-mapping with _bumpMap_	DifferentialGeometry dgs;	if (bumpMap)		Bump(bumpMap, dgGeom, dgShading, &dgs);	else		dgs = dgShading;	BSDF *bsdf = BSDF_ALLOC(BSDF)(dgs, dgGeom.nn);	bsdf->Add(BSDF_ALLOC(Lafortune)(Spectrum(diffuse), 3, xy, xy, z, e,		BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)));	return bsdf;}

Heitz::Heitz(const Spectrum &reflectance, Fresnel *f,             const HeitzDistribution &d)    : BxDF(BxDFType(BSDF_REFLECTION | BSDF_GLOSSY)),      mUseUniformSampling(false),   // TODO: test      R(reflectance), mDistribution(d), fresnel(f){}

// commented, dpl 10 august 2005FresnelBlend::FresnelBlend(const Spectrum &d,                           const Spectrum &s,						   MicrofacetDistribution *dist)	: BxDF(BxDFType(BSDF_REFLECTION | BSDF_GLOSSY)),	  Rd(d), Rs(s) {	distribution = dist;}

void Gen_Sample_f(BSDF* bsdf,    const Vector & wo, Vector* wi, float* pdf, Spectrum* f){    // only glossy or diffuse reflections (no specular reflections)    BxDFType inflags = BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE | BSDF_GLOSSY);    BxDFType outflags;    BSDFSample sample(rng);    *f = bsdf->Sample_f(wo, wi, sample, pdf, inflags, &outflags);    // double check bsdf->Pdf() gives us the same answer    Vector wiL = bsdf->WorldToLocal(*wi);    float wiCosTheta = wiL.z;    bool validSample = (wiCosTheta > 1e-7);        if (validSample) {        float verifyPdf = bsdf->Pdf(wo, *wi, inflags);        if (fabs(verifyPdf - *pdf) > 1e-4) {            fprintf(stderr, "BSDF::Pdf() doesn't match BSDF::Sample_f() !/n"                "  Sample_f pdf %.3f, Pdf pdf %.3f/n"                "  wo %.3f %.3f %.3f, wi %.3f %.3f %.3f/n",                *pdf, verifyPdf, wo[0], wo[1], wo[2], (*wi)[0], (*wi)[1], (*wi)[2]);            fprintf(stderr, "blah! validSample %d, wiCosTheta %.3f, wiL %.3f %.3f %.3f/n",                validSample, wiCosTheta, wiL[0], wiL[1], wiL[2]);        }    }}

Spectrum UniformSampleOneLight(const Scene *scene,        const Renderer *renderer, MemoryArena &arena, const Point &p,        const Normal &n, const Vector &wo, float rayEpsilon, float time,        BSDF *bsdf, const Sample *sample, RNG &rng, int lightNumOffset,        const LightSampleOffsets *lightSampleOffset,        const BSDFSampleOffsets *bsdfSampleOffset) {    // Randomly choose a single light to sample, _light_    int nLights = int(scene->lights.size());    if (nLights == 0) return Spectrum(0.);    int lightNum;    if (lightNumOffset != -1)        lightNum = Floor2Int(sample->oneD[lightNumOffset][0] * nLights);    else        lightNum = Floor2Int(rng.RandomFloat() * nLights);    lightNum = min(lightNum, nLights-1);    Light *light = scene->lights[lightNum];    // Initialize light and bsdf samples for single light sample    LightSample lightSample;    BSDFSample bsdfSample;    if (lightSampleOffset != NULL && bsdfSampleOffset != NULL) {        lightSample = LightSample(sample, *lightSampleOffset, 0);        bsdfSample = BSDFSample(sample, *bsdfSampleOffset, 0);    }    else {        lightSample = LightSample(rng);        bsdfSample = BSDFSample(rng);    }    return (float)nLights *        EstimateDirect(scene, renderer, arena, light, p, n, wo,                       rayEpsilon, time, bsdf, rng, lightSample,                       bsdfSample, BxDFType(BSDF_ALL & ~BSDF_SPECULAR));}

// BrushedMetal Method DefinitionsBSDF *BrushedMetal::GetBSDF(const DifferentialGeometry &dgGeom, const DifferentialGeometry &dgShading) const {	// Declare brushedmetal coefficients	static float diffuse[3] = { 0, 0, 0 };	static float xy0[3] =     {  -1.11854f, -1.11845f, -1.11999f  };	static float z0[3] =      {   1.01272f,  1.01469f,  1.01942f  };	static float e0[3] =      {  15.8708f,  15.6489f,  15.4571f   };	static float xy1[3] =     {  -1.05334f, -1.06409f, -1.08378f  };	static float z1[3] =      {   0.69541f,  0.662178f, 0.626672f };	static float e1[3] =      { 111.267f,   88.9222f,  65.2179f   };	static float xy2[3] =     {  -1.01684f, -1.01635f, -1.01529f  };	static float z2[3] =      {   1.00132f,  1.00112f,  1.00108f  };	static float e2[3] =      { 180.181f,  184.152f,  195.773f    };	static Spectrum xy[3] = { Spectrum(xy0), Spectrum(xy1), Spectrum(xy2) };	static Spectrum z[3] = { Spectrum(z0), Spectrum(z1), Spectrum(z2) };	static Spectrum e[3] = { Spectrum(e0), Spectrum(e1), Spectrum(e2) };	// Allocate _BSDF_, possibly doing bump-mapping with _bumpMap_	DifferentialGeometry dgs;	if (bumpMap)		Bump(bumpMap, dgGeom, dgShading, &dgs);	else		dgs = dgShading;	BSDF *bsdf = BSDF_ALLOC(BSDF)(dgs, dgGeom.nn);	bsdf->Add(BSDF_ALLOC(Lafortune)(Spectrum(*diffuse), 3, xy, xy, z, e,		BxDFType(BSDF_REFLECTION | BSDF_GLOSSY)));	return bsdf;}

// Integrator Utility FunctionsSpectrum UniformSampleAllLights(const Scene *scene,        const Renderer *renderer, MemoryArena &arena, const Point &p,        const Normal &n, const Vector &wo, float rayEpsilon,        float time, BSDF *bsdf, const Sample *sample, RNG &rng,        const LightSampleOffsets *lightSampleOffsets,        const BSDFSampleOffsets *bsdfSampleOffsets) {    Spectrum L(0.);    for (uint32_t i = 0; i < scene->lights.size(); ++i) {        Light *light = scene->lights[i];        int nSamples = lightSampleOffsets ?                       lightSampleOffsets[i].nSamples : 1;        // Estimate direct lighting from _light_ samples        Spectrum Ld(0.);        for (int j = 0; j < nSamples; ++j) {            // Find light and BSDF sample values for direct lighting estimate            LightSample lightSample;            BSDFSample bsdfSample;            if (lightSampleOffsets != NULL && bsdfSampleOffsets != NULL) {                lightSample = LightSample(sample, lightSampleOffsets[i], j);                bsdfSample = BSDFSample(sample, bsdfSampleOffsets[i], j);            }            else {                lightSample = LightSample(rng);                bsdfSample = BSDFSample(rng);            }            Ld += EstimateDirect(scene, renderer, arena, light, p, n, wo,                rayEpsilon, time, bsdf, rng, lightSample, bsdfSample,                BxDFType(BSDF_ALL & ~BSDF_SPECULAR));        }        L += Ld / nSamples;    }    return L;}

SpecularReflection::SpecularReflection( const Spectrum& R , const Fresnel* fresnel )	: BxDF( BxDFType( REFLECTION | SPECULAR ) )	, R( R )	, fresnel( fresnel ){}

// Felt Method DefinitionsBSDF *Felt::GetBSDF(const DifferentialGeometry &dgGeom, const DifferentialGeometry &dgShading) const {	// Declare felt coefficients	static float diffuse[3] = {  0.025865f,  0.025865f,  0.025865f};	static float xy0[3] =     { -0.304075f, -0.304075f, -0.304075f};	static float z0[3] =      { -0.065992f, -0.065992f, -0.065992f};	static float e0[3] =      {  3.047892f,  3.047892f,  3.047892f};	static float xy1[3] =     { -0.749561f, -0.749561f, -0.749561f};	static float z1[3] =      { -1.167929f, -1.167929f, -1.167929f};	static float e1[3] =      {  6.931827f,  6.931827f,  6.931827f};	static float xy2[3] =     {  1.004921f,  1.004921f,  1.004921f};	static float z2[3] =      { -0.205529f, -0.205529f, -0.205529f};	static float e2[3] =      { 94.117332f, 94.117332f, 94.117332f};	static Spectrum xy[3] = { Spectrum(xy0), Spectrum(xy1), Spectrum(xy2) };	static Spectrum z[3] = { Spectrum(z0), Spectrum(z1), Spectrum(z2) };	static Spectrum e[3] = { Spectrum(e0), Spectrum(e1), Spectrum(e2) };	// Allocate _BSDF_, possibly doing bump-mapping with _bumpMap_	DifferentialGeometry dgs;	if (bumpMap)		Bump(bumpMap, dgGeom, dgShading, &dgs);	else		dgs = dgShading;	BSDF *bsdf = BSDF_ALLOC(BSDF)(dgs, dgGeom.nn);	bsdf->Add(BSDF_ALLOC(Lafortune)(Spectrum(diffuse), 3, xy, xy, z, e,		BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)));	return bsdf;}

// Primer Method DefinitionsBSDF *Primer::GetBSDF(const DifferentialGeometry &dgGeom, const DifferentialGeometry &dgShading) const {	// Declare primer coefficients	static float diffuse[3] = {  0.118230f,  0.121218f,  0.133209f};	static float xy0[3] =     { -0.399286f, -1.033473f, -1.058104f};	static float z0[3] =      {  0.167504f,  0.009545f, -0.068002f};	static float e0[3] =      {  2.466633f,  7.637253f,  8.117645f};	static float xy1[3] =     { -1.041861f, -1.100108f, -1.087779f};	static float z1[3] =      {  0.014375f, -0.198147f, -0.053605f};	static float e1[3] =      {  7.993722f, 29.446268f, 41.988990f};	static float xy2[3] =     { -1.098605f, -0.379883f, -0.449038f};	static float z2[3] =      { -0.145110f,  0.159127f,  0.173224f};	static float e2[3] =      { 31.899719f,  2.372852f,  2.636161f};	static Spectrum xy[3] = { Spectrum(xy0), Spectrum(xy1), Spectrum(xy2) };	static Spectrum z[3] = { Spectrum(z0), Spectrum(z1), Spectrum(z2) };	static Spectrum e[3] = { Spectrum(e0), Spectrum(e1), Spectrum(e2) };	// Allocate _BSDF_, possibly doing bump-mapping with _bumpMap_	DifferentialGeometry dgs;	if (bumpMap)		Bump(bumpMap, dgGeom, dgShading, &dgs);	else		dgs = dgShading;	BSDF *bsdf = BSDF_ALLOC(BSDF)(dgs, dgGeom.nn);	bsdf->Add(BSDF_ALLOC(Lafortune)(Spectrum(diffuse), 3, xy, xy, z, e,		BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)));	return bsdf;}

// BluePaint Method DefinitionsBSDF *BluePaint::GetBSDF(const DifferentialGeometry &dgGeom, const DifferentialGeometry &dgShading) const {	// Declare bluepaint coefficients	static float diffuse[3] = {  0.3094f,    0.39667f,   0.70837f  };	static float xy0[3] =     {  0.870567f,  0.857255f,  0.670982f };	static float z0[3] =      {  0.803624f,  0.774290f,  0.586674f };	static float e0[3] =      { 21.820103f, 18.597755f,  7.472717f };	static float xy1[3] =     { -0.451218f, -0.406681f, -0.477976f };	static float z1[3] =      {  0.023123f,  0.017625f,  0.227295f };	static float e1[3] =      {  2.774499f,  2.581499f,  3.677653f };	static float xy2[3] =     { -1.031545f, -1.029426f, -1.026588f };	static float z2[3] =      {  0.706734f,  0.696530f,  0.687715f };	static float e2[3] =      { 66.899060f, 63.767912f, 57.489181f };	static Spectrum xy[3] = { Spectrum(xy0), Spectrum(xy1), Spectrum(xy2) };	static Spectrum z[3] = { Spectrum(z0), Spectrum(z1), Spectrum(z2) };	static Spectrum e[3] = { Spectrum(e0), Spectrum(e1), Spectrum(e2) };	// Allocate _BSDF_, possibly doing bump-mapping with _bumpMap_	DifferentialGeometry dgs;	if (bumpMap)		Bump(bumpMap, dgGeom, dgShading, &dgs);	else		dgs = dgShading;	BSDF *bsdf = BSDF_ALLOC(BSDF)(dgs, dgGeom.nn);	bsdf->Add(BSDF_ALLOC(Lafortune)(Spectrum(diffuse), 3, xy, xy, z, e,		BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)));	return bsdf;}

RGB BSDF::f(const Vector &woWorld, const Vector &wiWorld,		BxDFType flags) const {	Vector wo = WorldToLocal(woWorld);	Vector wi = WorldToLocal(wiWorld);	if (Dot(woWorld, mNG) * Dot(wiWorld, mNG) > 0) { //相乘大于零说明在同一半球		flags = BxDFType(flags & ~BSDF_TRANSMISSION); //去除折射	} else {		flags = BxDFType(flags & ~BSDF_REFLECTION); //去除反射	}	RGB f = 0;	for (int i = 0; i < mNumBxdf; ++i) {		if (mBxdfs[i]->MatchesFlag(flags)) {			f += mBxdfs[i]->f(wo, wi);		}	}	return f;}

SpecularTransmission::SpecularTransmission(SpectrumCoef_d i_transmittance, double i_refractive_index_inner, double i_refractive_index_outer):    BxDF(BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR)), m_transmittance(i_transmittance),    m_refractive_index_inner(i_refractive_index_inner), m_refractive_index_outer(i_refractive_index_outer),    m_fresnel(i_refractive_index_inner, i_refractive_index_outer){    ASSERT(InRange(i_transmittance,0.0,1.0));    ASSERT(i_refractive_index_inner>0.0 && i_refractive_index_outer>0.0);}

// Metropolis Method Definitionsstatic 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;}

OrenNayar::OrenNayar(SpectrumCoef_d i_reflectance, double i_sigma):BxDF(BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE)), m_reflectance(i_reflectance)  {  ASSERT(InRange(i_reflectance,0.0,1.0));  ASSERT(i_sigma>=0.0 && i_sigma<=1.0);  double sigma_sqr = i_sigma*i_sigma;  m_A = 1.0 - (0.5 * sigma_sqr / (sigma_sqr + 0.33));  m_B = 0.45 * sigma_sqr / (sigma_sqr + 0.09);  }

Spectrum DipoleSubsurfaceIntegrator::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;    Spectrum rho_dr = Spectrum(1.0f);    // Evaluate BSSRDF and possibly compute subsurface scattering    BSSRDF *bssrdf = isect.GetBSSRDF(ray, arena);    if (bssrdf && octree) {	        Spectrum sigma_a  = bssrdf->sigma_a();        Spectrum sigmap_s = bssrdf->sigma_prime_s();        Spectrum sigmap_t = sigmap_s + sigma_a;        if (!sigmap_t.IsBlack()) {            // Use hierarchical integration to evaluate reflection from dipole model            PBRT_SUBSURFACE_STARTED_OCTREE_LOOKUP(const_cast<Point *>(&p));            DiffusionReflectance Rd(sigma_a, sigmap_s, bssrdf->eta());            Spectrum Mo = octree->Mo(octreeBounds, p, Rd, maxError);            FresnelDielectric fresnel(1.f, bssrdf->eta());            Spectrum Ft = Spectrum(1.f) - fresnel.Evaluate(AbsDot(wo, n));            float Fdt = 1.f - Fdr(bssrdf->eta());	    // modulate SSS contribution by rho_dr            //L += (INV_PI * Ft) * (Fdt * Mo);	    rho_dr = wet->integrate_BRDF(bsdf, ray.d, 10, BxDFType(BSDF_REFLECTION | BSDF_GLOSSY));	    L += (INV_PI * Ft) * (Fdt * Mo) * (Spectrum(1.0f) - rho_dr);	    //L += (INV_PI * Ft) * (Fdt * Mo) * (Spectrum(0.0f));	                PBRT_SUBSURFACE_FINISHED_OCTREE_LOOKUP();        }    }    L += UniformSampleAllLights(scene, renderer, arena, p, n,        wo, isect.rayEpsilon, ray.time, bsdf, sample, rng, lightSampleOffsets,        bsdfSampleOffsets);    if (ray.depth < maxSpecularDepth) {        // Trace rays for specular reflection and refraction.      //TODO: this has no effect?        L += SpecularReflect(ray, bsdf, rng, isect, renderer, scene,			     sample, arena);        L += SpecularTransmit(ray, bsdf, rng, isect,			      renderer, scene, sample, arena);    }    return L;}

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,        BSDF *bsdf, RNG &rng, const LightSample &lightSample,        const BSDFSample &bsdfSample) {    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);        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);                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 flags = BxDFType(BSDF_ALL & ~BSDF_SPECULAR);        Spectrum f = bsdf->Sample_f(wo, &wi, bsdfSample, &bsdfPdf, flags);        if (!f.IsBlack() && bsdfPdf > 0.) {            lightPdf = light->Pdf(p, 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, rayEpsilon, INFINITY, time);                if (scene->Intersect(ray, &lightIsect)) {                    if (lightIsect.primitive->GetAreaLight() == light)                        Li = lightIsect.Le(-wi);                }                else                    Li = light->Le(ray);                if (!Li.IsBlack()) {                    Li *= renderer->Transmittance(scene, ray, NULL, rng, arena);                    Ld += f * Li * AbsDot(wi, n) * weight / bsdfPdf;                }            }        }    }    return Ld;}

Spectrum IrradianceCacheIntegrator::Li(const Scene *scene,        const Renderer *renderer, const RayDifferential &ray, const Intersection &isect,        const Sample *sample, RNG &rng, MemoryArena &arena) const {    Spectrum L(0.);    // Evaluate BSDF at hit point    BSDF *bsdf = isect.GetBSDF(ray, arena);    Vector wo = -ray.d;    const Point &p = bsdf->dgShading.p;    const Normal &n = bsdf->dgShading.nn;    L += isect.Le(wo);    // Compute direct lighting for irradiance cache    L += UniformSampleAllLights(scene, renderer, arena, p, n, wo,             isect.rayEpsilon, ray.time, bsdf, sample, rng,             lightSampleOffsets, bsdfSampleOffsets);    // Compute indirect lighting for irradiance cache    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);    }    // Estimate indirect lighting with irradiance cache    Normal ng = isect.dg.nn;    ng = Faceforward(ng, wo);    // Compute pixel spacing in world space at intersection point    float pixelSpacing = sqrtf(Cross(isect.dg.dpdx, isect.dg.dpdy).Length());    BxDFType flags = BxDFType(BSDF_REFLECTION | BSDF_DIFFUSE | BSDF_GLOSSY);    L += indirectLo(p, ng, pixelSpacing, wo, isect.rayEpsilon,                    bsdf, flags, rng, scene, renderer, arena);    flags = BxDFType(BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);    L += indirectLo(p, -ng, pixelSpacing, wo, isect.rayEpsilon,                    bsdf, flags, rng, scene, renderer, arena);    return L;}

Spectrum_d PhotonLTEIntegrator::_LookupCausticRadiance(const BSDF *ip_bsdf, const DifferentialGeometry &i_dg, const Vector3D_d &i_direction, ThreadSpecifics i_ts) const  {  ASSERT(ip_bsdf);  ASSERT(i_direction.IsNormalized());  MemoryPool *p_pool = i_ts.mp_pool;  if (mp_photon_maps->GetCausticMap() == NULL)    return Spectrum_d();  BxDFType non_specular = BxDFType(BSDF_REFLECTION | BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);  if (ip_bsdf->GetComponentsNum(non_specular) == 0)    return Spectrum_d();  // Allocate array for the nearest photons.  NearestPhoton *p_nearest_photons = (NearestPhoton*)p_pool->Alloc(m_params.m_caustic_lookup_photons_num * sizeof(NearestPhoton));  PhotonFilter filter(i_dg.m_point, i_dg.m_geometric_normal, MAX_NORMAL_DEVIATION_COS);  size_t photons_found = mp_photon_maps->GetCausticMap()->GetNearestPoints(i_dg.m_point, m_params.m_caustic_lookup_photons_num, p_nearest_photons, filter, m_params.m_max_caustic_lookup_dist);  if (photons_found == 0)    return Spectrum_d();  double max_dist_sqr = 0.0;  Spectrum_d radiance;  for(size_t i=0;i<photons_found;++i)    {    Point3D_d photon_position = Convert<double>( p_nearest_photons[i].mp_point->m_point );    Vector3D_d photon_direction = p_nearest_photons[i].mp_point->m_incident_direction.ToVector3D<double>();    Spectrum_d photon_weight = Convert<double>( p_nearest_photons[i].mp_point->m_weight );    double tmp_dist_sqr = Vector3D_d(photon_position - i_dg.m_point).LengthSqr();    if (tmp_dist_sqr > max_dist_sqr) max_dist_sqr = tmp_dist_sqr;    double kernel = _PhotonKernel( Vector3D_d(photon_position-i_dg.m_point).LengthSqr(), max_dist_sqr);    radiance += ip_bsdf->Evaluate(photon_direction, i_direction) * photon_weight * kernel;    }  if (photons_found<m_params.m_caustic_lookup_photons_num || max_dist_sqr==0.0)    max_dist_sqr = m_params.m_max_caustic_lookup_dist * m_params.m_max_caustic_lookup_dist;  else    /*    Since the max_dist_sqr is exactly equal to the squared distance to the farthest photon we need to multiply the area by the correcting factor.    The easy way to understand it is the following. Think of what will happen if we decrease the radius a little bit.    The farthest photon will drop out while the area won't change significantly. Thus the resulting radiance value would change by the value brought by the farthest photon.    We need to increase the total area by a half of a single photon's area.    */    max_dist_sqr *= (photons_found) / (photons_found-0.5);  return radiance / (mp_photon_maps->GetNumberOfCausticPaths() * max_dist_sqr);  }

Spectrum MetropolisRenderer::Lpath(const Scene *scene,        const PathVertex *cameraPath, int cameraPathLength,        MemoryArena &arena, const vector<LightingSample> &samples,        RNG &rng, float time, const Distribution1D *lightDistribution,        const RayDifferential &eRay, const Spectrum &eAlpha) const {    PBRT_MLT_STARTED_LPATH();    Spectrum L = 0.;    bool previousSpecular = true, allSpecular = true;    for (int i = 0; i < cameraPathLength; ++i) {        // Initialize basic variables for camera path vertex        const PathVertex &vc = cameraPath[i];        const Point &pc = vc.bsdf->dgShading.p;        const Normal &nc = vc.bsdf->dgShading.nn;        // Add emitted light from vertex if appropriate        if (previousSpecular && (directLighting == NULL || !allSpecular))            L += vc.alpha * vc.isect.Le(vc.wPrev);        // Compute direct illumination for Metropolis path vertex        Spectrum Ld(0.f);        if (directLighting == NULL || !allSpecular) {            // Choose light and call _EstimateDirect()_ for Metropolis vertex            const LightingSample &ls = samples[i];            float lightPdf;            uint32_t lightNum = lightDistribution->SampleDiscrete(ls.lightNum,                                                                  &lightPdf);            const Light *light = scene->lights[lightNum];            PBRT_MLT_STARTED_ESTIMATE_DIRECT();                        Ld = vc.alpha *                 EstimateDirect(scene, this, arena, light, pc, nc, vc.wPrev,                                vc.isect.rayEpsilon, time, vc.bsdf, rng, NULL,                                ls.lightSample, ls.bsdfSample,                                BxDFType(BSDF_ALL & ~BSDF_SPECULAR)) / lightPdf;            PBRT_MLT_FINISHED_ESTIMATE_DIRECT();        }        previousSpecular = vc.specularBounce;        allSpecular &= previousSpecular;        L += Ld;    }    // Add contribution of escaped ray, if any    if (!eAlpha.IsBlack() && previousSpecular &&        (directLighting == NULL || !allSpecular))        for (uint32_t i = 0; i < scene->lights.size(); ++i)           L += eAlpha * scene->lights[i]->Le(eRay);    PBRT_MLT_FINISHED_LPATH();    return L;}

Spectrum SamplerIntegrator::SpecularTransmit(    const RayDifferential &ray, const SurfaceInteraction &isect,    const Scene &scene, Sampler &sampler, MemoryArena &arena, int depth) const {    Vector3f wo = isect.wo, wi;    Float pdf;    const Point3f &p = isect.p;    const Normal3f &ns = isect.shading.n;    const BSDF &bsdf = *isect.bsdf;    Spectrum f = bsdf.Sample_f(wo, &wi, sampler.Get2D(), &pdf,                               BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));    Spectrum L = Spectrum(0.f);    if (pdf > 0.f && !f.IsBlack() && AbsDot(wi, ns) != 0.f) {        // Compute ray differential _rd_ for specular transmission        RayDifferential rd = isect.SpawnRay(wi);        if (ray.hasDifferentials) {            rd.hasDifferentials = true;            rd.rxOrigin = p + isect.dpdx;            rd.ryOrigin = p + isect.dpdy;            Float eta = bsdf.eta;            Vector3f w = -wo;            if (Dot(wo, ns) < 0) eta = 1.f / eta;            Normal3f dndx = isect.shading.dndu * isect.dudx +                            isect.shading.dndv * isect.dvdx;            Normal3f dndy = isect.shading.dndu * isect.dudy +                            isect.shading.dndv * isect.dvdy;            Vector3f dwodx = -ray.rxDirection - wo,                     dwody = -ray.ryDirection - wo;            Float dDNdx = Dot(dwodx, ns) + Dot(wo, dndx);            Float dDNdy = Dot(dwody, ns) + Dot(wo, dndy);            Float mu = eta * Dot(w, ns) - Dot(wi, ns);            Float dmudx =                (eta - (eta * eta * Dot(w, ns)) / Dot(wi, ns)) * dDNdx;            Float dmudy =                (eta - (eta * eta * Dot(w, ns)) / Dot(wi, ns)) * dDNdy;            rd.rxDirection =                wi + eta * dwodx - Vector3f(mu * dndx + dmudx * ns);            rd.ryDirection =                wi + eta * dwody - Vector3f(mu * dndy + dmudy * ns);        }        L = f * Li(rd, scene, sampler, arena, depth + 1) * AbsDot(wi, ns) / pdf;    }    return L;}

Spectrum SpecularTransmit(const RayDifferential &ray, BSDF *bsdf,        RNG &rng, const Intersection &isect, const Renderer *renderer,        const Scene *scene, const Sample *sample, MemoryArena &arena) {    Vector wo = -ray.d, wi;    float pdf;    const Point &p = bsdf->dgShading.p;    const Normal &n = bsdf->dgShading.nn;    Spectrum f = bsdf->Sample_f(wo, &wi, BSDFSample(rng), &pdf,                               BxDFType(BSDF_TRANSMISSION | BSDF_SPECULAR));    Spectrum L = 0.f;    if (pdf > 0.f && !f.IsBlack() && AbsDot(wi, n) != 0.f) {        // Compute ray differential _rd_ for specular transmission        RayDifferential rd(p, wi, ray, isect.rayEpsilon);        if (ray.hasDifferentials) {            rd.hasDifferentials = true;            rd.rxOrigin = p + isect.dg.dpdx;            rd.ryOrigin = 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.rxDirection - wo, dwody = -ray.ryDirection - 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.rxDirection = wi + eta * dwodx - Vector(mu * dndx + dmudx * n);            rd.ryDirection = wi + eta * dwody - Vector(mu * dndy + dmudy * n);        }        PBRT_STARTED_SPECULAR_REFRACTION_RAY(const_cast<RayDifferential *>(&rd));        Spectrum Li = renderer->Li(scene, rd, sample, rng, arena);        L = f * Li * AbsDot(wi, n) / pdf;        PBRT_FINISHED_SPECULAR_REFRACTION_RAY(const_cast<RayDifferential *>(&rd));    }    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;}