PBRT_V2 总结记录 <100> IrradianceCacheIntegrator

class IrradianceCacheIntegrator : public SurfaceIntegrator {
public:
    // IrradianceCacheIntegrator Public Methods
    IrradianceCacheIntegrator(float minwt, float minsp, float maxsp,
                              float maxang, int maxspec, int maxind, int ns) {
        minWeight = minwt;
        minSamplePixelSpacing = minsp;
        maxSamplePixelSpacing = maxsp;
        cosMaxSampleAngleDifference = cosf(Degrees(maxang));
        nSamples = ns;
        maxSpecularDepth = maxspec;
        maxIndirectDepth = maxind;
        mutex = RWMutex::Create();
        lightSampleOffsets = NULL;
        bsdfSampleOffsets = NULL;
    }
    ~IrradianceCacheIntegrator();
    Spectrum Li(const Scene *scene, const Renderer *renderer,
        const RayDifferential &ray, const Intersection &isect,
        const Sample *sample, RNG &rng, MemoryArena &arena) const;
    void RequestSamples(Sampler *sampler, Sample *sample, const Scene *scene);
    void Preprocess(const Scene *, const Camera *, const Renderer *);
private:
    // IrradianceCacheIntegrator Private Data
    float minSamplePixelSpacing, maxSamplePixelSpacing;
    float minWeight, cosMaxSampleAngleDifference;
    int nSamples, maxSpecularDepth, maxIndirectDepth;
    mutable RWMutex *mutex;

    // Declare sample parameters for light source sampling
    LightSampleOffsets *lightSampleOffsets;
    BSDFSampleOffsets *bsdfSampleOffsets;
    mutable Octree<IrradianceSample *> *octree;

    // IrradianceCacheIntegrator Private Methods
    Spectrum indirectLo(const Point &p, const Normal &ng, float pixelSpacing,
        const Vector &wo, float rayEpsilon,BSDF *bsdf, BxDFType flags, RNG &rng,
        const Scene *scene, const Renderer *renderer, MemoryArena &arena) const;
    bool interpolateE(const Scene *scene,
            const Point &p, const Normal &n, Spectrum *E, Vector *wi) const;
    Spectrum pathL(Ray &r, const Scene *scene, const Renderer *renderer,
        RNG &rng, MemoryArena &arena) const;
};

类的作用：

（Path Tracing 实现中，一个Pixel需要大量的Sample才可以消除noise，这样会涉及大量的计算，所以，出现其他的办法给去解决这些问题，这些方法一般就是重新使用之前计算过的结果，而不是每一次都是重新计算，这样就可以节省计算量。所以这些方法 属于 biased approaches）

With unbiased light transport algorithms like path tracing, some scenes can take a large
number of rays (and corresponding compute time) to generate images without objectionable
noise. One approach to this problem has been the development of biased approaches
to solving the LTE. These approaches generally reuse previously computed results
over multiple exitant radiance computations, even when the values used don’t give
the precise quantity that needs to be computed (for example, by reusing an illumination
value from a nearby point under the assumption that illumination is slowly changing).

The irradiance caching algorithm is based on the observation that, while direct lighting
often changes rapidly from point to point (e.g., consider a hard shadow edge), indirect
lighting is often much more slowly changing. Therefore, if we compute an accurate
representation of indirect light at a sparse(稀少) set of sample points in the scene and then
interpolate nearby samples to compute an approximate representation of indirect light
at any particular point being shaded, we can expect that the error introduced by not
recomputing indirect lighting everywhere shouldn’t be too bad, and we can achieve a
substantial computational savings by not recomputing this information at every point.

 

个人理解的算法思路：

（Irradiance caching algorithm 缓存一些交点的 E(iradiance) 和 wAvg(average incident direction，入射方向的平均值)，E 是利用 Path Tracing 算法来进行计算的，那么并不是所有的交点都需要计算这个 E 值和 wAvg 值，对于某一些 符合条件的交点，可以直接从 缓存中 取出 与它 相邻的 交点的 E 和 wAvg，拿缓存中的交点的 E和 wAvg来计算 其他交点的 E 和 wAvg，得到 E和 wAvg之后，就可以计算出这个E 对 方向 w0 所产生的 radiance，这个radiance就是这个点所产生的间接光）

In the irradiance caching algorithm implemented here, indirect lighting is computed at
a subset of the points in the scene that are shaded (as opposed to a predetermined set of
points) and stored in a spatial data structure. When exitant radiance at a point is being
computed, the cache is first searched for one or more acceptable nearby samples, using a
set of error metrics to determine if the already existing samples are acceptable. If not
enough acceptable samples are found, a new one is computed and added to the data
structure.

In order to have a compact representation of indirect light, this algorithm only stores the
irradiance and the average direction of incident irradiance at each point, rather than a
directionally varying radiance distribution, thus reducing the representation of incident
illumination to just a single Spectrum and a single Vector.

 

1.  Spectrum Li(const Scene *scene, const Renderer *renderer,
        const RayDifferential &ray, const Intersection &isect,
        const Sample *sample, RNG &rng, MemoryArena &arena) const;

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;
}

 细节：

a.

    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);

作用：

计算交点的直接光照。

 

b.

    // 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);

作用：

（这里就是处理 reflective 和 transmissive 的 间接光，如果计算 间接光，核心操作都放在了 indirectLo 函数中）

If the surface has both reflective and transmissive nonspecular components, then reflection
and transmission must be handled separately with two independent irradiance values,
since the irradiance values from the hemisphere that is needed for reflective surfaces
and the hemisphere on the opposite side of the surface that is needed for transmissive
surfaces can be completely different. The indirectLo() method handles either case, interpolating
a value using the cache or computing a new value and using it to compute
exitant radiance due to incident irradiance.

 

2.

 Spectrum indirectLo(const Point &p, const Normal &ng, float pixelSpacing,
        const Vector &wo, float rayEpsilon,BSDF *bsdf, BxDFType flags, RNG &rng,
        const Scene *scene, const Renderer *renderer, MemoryArena &arena) const;

Spectrum IrradianceCacheIntegrator::indirectLo(const Point &p,
        const Normal &ng, float pixelSpacing, const Vector &wo,
        float rayEpsilon, BSDF *bsdf, BxDFType flags, RNG &rng,
        const Scene *scene, const Renderer *renderer,
        MemoryArena &arena) const {
    if (bsdf->NumComponents(flags) == 0)
        return Spectrum(0.);
    Spectrum E;
    Vector wi;
    // Get irradiance _E_ and average incident direction _wi_ at point _p_
    if (!interpolateE(scene, p, ng, &E, &wi)) {
        // Compute irradiance at current point
        PBRT_IRRADIANCE_CACHE_STARTED_COMPUTING_IRRADIANCE(const_cast<Point *>(&p), const_cast<Normal *>(&ng));
        uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
        float minHitDistance = INFINITY;
        Vector wAvg(0,0,0);
        Spectrum LiSum = 0.f;
        for (int i = 0; i < nSamples; ++i) {
            // Sample direction for irradiance estimate ray
            float u[2];
            Sample02(i, scramble, u);
            Vector w = CosineSampleHemisphere(u[0], u[1]);
            RayDifferential r(p, bsdf->LocalToWorld(w), rayEpsilon);
            r.d = Faceforward(r.d, ng);

            // Trace ray to sample radiance for irradiance estimate
            PBRT_IRRADIANCE_CACHE_STARTED_RAY(&r);
            Spectrum L = pathL(r, scene, renderer, rng, arena);
            LiSum += L;
            wAvg += r.d * L.y();
            minHitDistance = min(minHitDistance, r.maxt);
            PBRT_IRRADIANCE_CACHE_FINISHED_RAY(&r, r.maxt, &L);
        }
        E = (M_PI / float(nSamples)) * LiSum;
        PBRT_IRRADIANCE_CACHE_FINISHED_COMPUTING_IRRADIANCE(const_cast<Point *>(&p), const_cast<Normal *>(&ng));

        // Add computed irradiance value to cache

        // Compute irradiance sample's contribution extent and bounding box
        float maxDist = maxSamplePixelSpacing * pixelSpacing;
        float minDist = minSamplePixelSpacing * pixelSpacing;
        float contribExtent = Clamp(minHitDistance / 2.f, minDist, maxDist);
        BBox sampleExtent(p);
        sampleExtent.Expand(contribExtent);
        PBRT_IRRADIANCE_CACHE_ADDED_NEW_SAMPLE(const_cast<Point *>(&p), const_cast<Normal *>(&ng), contribExtent, &E, &wAvg, pixelSpacing);

        // Allocate _IrradianceSample_, get write lock, add to octree
        IrradianceSample *sample = new IrradianceSample(E, p, ng, wAvg,
                                                        contribExtent);
        RWMutexLock lock(*mutex, WRITE);
        octree->Add(sample, sampleExtent);
        wi = wAvg;
    }

    // Compute reflected radiance due to irradiance and BSDF
    if (wi.LengthSquared() == 0.f) return Spectrum(0.);
    return bsdf->f(wo, Normalize(wi), flags) * E;
}

作用：

（indirectLo 主要就是计算BSDF 对应 flag的 间接光）

The indirectLo() method returns the reflected radiance for the BSDF components specified
by the flags parameter, computed by finding the irradiance at the given point in
the hemisphere by the given normal. A surprising number of additional parameters to
the method are needed; many of them are needed for cases where a new irradiance value
needs to be computed and added to the cache.

The indirectLo() method , interpolating
a value using the cache or computing a new value and using it to compute
exitant radiance due to incident irradiance.

 

细节：

a.

if (!interpolateE(scene, p, ng, &E, &wi))

作用：

（在 interpolateE 方法中，主要的作用就是在 缓存中 根据新的交点和新的法线 来寻找缓存中合适的 E 和 wAvg，找到合适的E 和 wAvg 之后，就进行 对应的插值操作，计算得到新的交点的 E 和 wAvg, 但是，目前先来看看，怎么往缓存添加 数据，）

The interpolateE() method attempts to interpolate the irradiance values already in
the octree to approximate the irradiance at the given point in the hemisphere with the
given normal. If suitable irradiance samples aren’t found in the cache, a new sample is
computed and added. In either case, E is initialized with an approximation of irradiance
for this point, wi is initialized to the average direction of incident radiance, and exitant
radiance can be computed with the approximation developed above.

 

b. 这里就是在缓存中找不到合适的交点，就计算目前交点的 E 和 wAvg。

        // Compute irradiance at current point
        
        uint32_t scramble[2] = { rng.RandomUInt(), rng.RandomUInt() };
        float minHitDistance = INFINITY;
        Vector wAvg(0,0,0);
        Spectrum LiSum = 0.f;
        for (int i = 0; i < nSamples; ++i) {
            // Sample direction for irradiance estimate ray
            float u[2];
            Sample02(i, scramble, u);
            Vector w = CosineSampleHemisphere(u[0], u[1]);
            RayDifferential r(p, bsdf->LocalToWorld(w), rayEpsilon);
            r.d = Faceforward(r.d, ng);

            // Trace ray to sample radiance for irradiance estimate
            PBRT_IRRADIANCE_CACHE_STARTED_RAY(&r);
            Spectrum L = pathL(r, scene, renderer, rng, arena);
            LiSum += L;
            wAvg += r.d * L.y();
            minHitDistance = min(minHitDistance, r.maxt);
            PBRT_IRRADIANCE_CACHE_FINISHED_RAY(&r, r.maxt, &L);
        }
        E = (M_PI / float(nSamples)) * LiSum;
        

        

作用：

（上面的代码计算交点的E 和wAvg， 通过 自己在半球面上发射 nSample 条Ray，每一条Ray都是 执行 pathL 方法，也就是 Path Tracing 来计算 Li，再根据

来计算这个交点最后的 E 值。而 wAvg的计算 也是在上面进行的 ： wAvg += r.d * L.y();）

 

If the interpolateE() method isn’t able to find enough nearby irradiance samples of
good enough quality, it computes a new irradiance value and adds it to the cache.

To compute a new irradiance estimate, we need to estimate the value of the integral

To compute the irradiance estimate from the radiance values along the sample rays, the
implementation uses the standard Monte Carlo estimator:

Because the rays are sampled from a cosine-weighted distribution, p(ω) = cos θ/π, so

 

c.

        // Add computed irradiance value to cache

        // Compute irradiance sample's contribution extent and bounding box
        float maxDist = maxSamplePixelSpacing * pixelSpacing;
        float minDist = minSamplePixelSpacing * pixelSpacing;
        float contribExtent = Clamp(minHitDistance / 2.f, minDist, maxDist);
        BBox sampleExtent(p);
        sampleExtent.Expand(contribExtent);
        PBRT_IRRADIANCE_CACHE_ADDED_NEW_SAMPLE(const_cast<Point *>(&p), const_cast<Normal *>(&ng), contribExtent, &E, &wAvg, pixelSpacing);

        // Allocate _IrradianceSample_, get write lock, add to octree
        IrradianceSample *sample = new IrradianceSample(E, p, ng, wAvg,
                                                        contribExtent);
        RWMutexLock lock(*mutex, WRITE);
        octree->Add(sample, sampleExtent);
        wi = wAvg;



struct IrradianceSample {
    // IrradianceSample Constructor
    IrradianceSample() { }
    IrradianceSample(const Spectrum &e, const Point &P, const Normal &N,
        const Vector &pd, float md) : E(e), n(N), p(P), wAvg(pd) {
        maxDist = md;
    }
    Spectrum E;
    Normal n;
    Point p;
    Vector wAvg;
    float maxDist;
};

作用：

（这里可以简单理解，这里就是 新的交点的信息(IrradianceSample) 加入到 缓存到，IrradianceSample 会保存这个交点的 坐标，法线，wAvg，这个交点和附近最近物体的距离的一半。(用于判断是否与新的交点相似)）

There are three criteria that determine the region in which an irradiance sample may be
valid. The first two criteria are based on pixel spacing constraints. As described above,
the pixel spacing criteria ensure that, no matter what, we don’t deposit samples at a
rate greater than minSamplePixelSpacing as measured on the image plane, and that,
no matter what, we place a sample at the minimum rate of maxSamplePixelSpacing.
Those parameters are converted to world-space distances here by multiplying by the
provided value of the length of the distance between pixels in world space at the point
pixelSpacing.

The third criterion is based on the minimum distance to the closest intersection found
when computing the irradiance sample. The idea is that if there is a nearby occluding
object, then the indirect irradiance is likely to be changing quickly around the sample.
For example, near the corner where a wall meets the ceiling, irradiance will in general be
changing more rapidly than it will be in the middle of the ceiling. This term captures this
effect, reducing the region of applicability for samples that are near occluders.

Given the sampleExtent bounding box that bounds the sample’s valid region, we can
allocate an irradiance sample object and add it to the octree. Note that we acquire a write
lock to the octree before updating it, thus ensuring that other threads aren’t concurrently
traversing the tree as it is being modified here.

 

d.

bool IrradianceCacheIntegrator::interpolateE(const Scene *scene,
        const Point &p, const Normal &n, Spectrum *E,
        Vector *wi) const {
    if (!octree) return false;
    
    IrradProcess proc(p, n, minWeight, cosMaxSampleAngleDifference);
    RWMutexLock lock(*mutex, READ);
    octree->Lookup(p, proc);
    
    if (!proc.Successful()) return false;
    *E = proc.GetIrradiance();
    *wi = proc.GetAverageDirection();
    return true;
}

作用：

（在interpolateE 的方法中，在缓存中寻找合适的交点，插值交点的数据，得出新的数据给新的交点用，那么怎么去寻找合适的交点呢，主要的工作都是在 Octee::Lookup中进行，Lookup 会执行 IrradProcess 的 operator()函数，比较缓存中的每一个交点，看看哪一些交点比较合适，从而计算出新的交点数据）

Much of the work in the interpolateE() method is done by the Octree::Lookup()
method, which traverses the nodes of the octree that the given point is inside and calls
the operator() method of the IrradProcess object for each IrradianceSample in each

of these nodes. IrradProcess decides if each sample is acceptable and accumulates the
value of the interpolated result. Note that a reader lock is acquired before the lookup is
performed; doing so ensures that another thread can’t insert new values into the octree
as it is being traversed in this thread.

 

e.

bool IrradProcess::operator()(const IrradianceSample *sample) {
    // Compute estimate error term and possibly use sample
    float perr = Distance(p, sample->p) / sample->maxDist;
    float nerr = sqrtf((1.f - Dot(n, sample->n)) /
                       (1.f - cosMaxSampleAngleDifference));
    float err = max(perr, nerr);
    
    if (err < 1.) {
        ++nFound;
        float wt = 1.f - err;
        E += wt * sample->E;
        wAvg += wt * sample->wAvg;
        sumWt += wt;
    }
    return true;
}

Spectrum GetIrradiance() const { return E / sumWt; }

Vector GetAverageDirection() const { return wAvg; }

作用：

（perr 其实就是 计算 新的交点和 缓存交点 的距离，与 缓存交点最近碰撞物体距离的一半 进行 比值，其实就是 利用距离 来确定两个交点的关系，只有相近缓存交点才会被挑选出来，

nerr 其实就是 比较 新的交代呢和 缓存交代的 法线的差距，只有法线相接近的缓存交点才会被挑选出来。

如果perr 和 nerr 都不大的情况下，那就是，这个缓存交点与 新的交点相似的情况下，就可以 利用加权平均来计算 这个新的交点的 E 和wAvg

）

Each time the IrradProcess callback method is called by Octree::Lookup(), it is passed a
pointer to an irradiance sample from the octree. The method computes an error term
that tries to estimate the error from using the sample at the shading point, which is
compared to user-supplied error limits; if the sample passes the test, it’s added to a
running weighted sum of samples that contribute to the point being shaded.

Two error terms are computed to determine whether a particular candidate irradiance
sample can reasonably be used at the lookup point. The first, perr, is based on the ratio
of the distance from the lookup point to the point where the sample was computed to
the sample’s maximum allowed distance value.

The second error term, nerr, is based on the angle between the surface normal of the
lookup point and the normal direction used for computing the irradiance sample. The
estimated error increases as the angle approaches the maximum allowed angle. For efficiency,
an error term based on the cosine of the angle is used, so that only a dot product
is needed to compare their orientations.

Samples with acceptable error are then added to the weighted averages. Each sample
is weighted by one minus its error, thus ensuring that low-error samples have higher
contribution to the final interpolated result.

The minimum weight and cosMaxSampleAngleDifference values control the error metric
used in determining if a particular irradiance sample can be used at a given point without
introducing excessive error. Their use will be explained in the code that uses them.

 

f.

// Compute reflected radiance due to irradiance and BSDF
    if (wi.LengthSquared() == 0.f) return Spectrum(0.);
    return bsdf->f(wo, Normalize(wi), flags) * E;

作用：

（根据公式，

利用 wAvg 和 E 就可以计算用 Lo，也就是可以计算出在 wo 方向上的 radiance,这里就顺利地计算出 间接光的 radiance）

Following the approach introduced by Tabellion and Lamorlette (2004), we can approximate
the incident irradiance by incident radiance along a single direction—in some sense
converting it to an equivalent directional light source. Doing so maintains a very rough
approximation of the directional distribution of incident radiance at the point being
shaded. The end result for perfectly diffuse surfaces is the same as if we just multiplied
the irradiance by the reflectance, but for moderately glossy surfaces this approximation
produces better results than assuming a uniform incident distribution of illumination.

Given an irradiance value E and an incident direction ωavg, we can compute the amount
of radiance L that must arrive along that direction to give the same irradiance:

so Lavg = E/|cos θavg|. Substituting into the scattering equation, we have

Thus, computing the reflected radiance is straightforward. The only complication is that
the length of the weighted average direction of incident radiance function may be zero;
this happens when there is no irradiance at the point. In this case, a black spectrum
is returned immediately to avoid calling Vector::Normalize() on a zero-length vector,
which would in turn return not-a-number values.