PBRT_V2 总结记录 <110> SHTerms 和 SHTerms 和 SHEvaluate

1.

inline int SHTerms(int lmax) {
    return (lmax + 1) * (lmax + 1);
}

作用：

（根据 lmax 的，计算 basis functin 的数目，参考：《PBRT_V2 总结记录 <109> Spherical Harmonic (SH) Basis Functions》就可以知道， Lmax = 1 的时候，basis function 的个数是 Lmax = 0 和 Lmax = 1 的总和）

SHTerms() gives the total number of basis functions (or equivalently, the total number
of coefficients) that there are for a given lmax value. This function is used frequently
for tasks like computing the size of arrays needed to hold SH coefficients.

 

2.

inline int SHIndex(int l, int m) {
    return l*l + l + m;
}

作用：

（参考：《PBRT_V2 总结记录 <109> Spherical Harmonic (SH) Basis Functions》，根据 l,m 来 计算 basis function 的 index）

Given a 1D array to be used for storing SH basis function coefficients, the SHIndex()
computes the corresponding index into the array for a particular (l , m) basis function.

 

3.

// Spherical Harmonics Definitions
void SHEvaluate(const Vector &w, int lmax, float *out) {
    if (lmax > 28) {
        Error("SHEvaluate() runs out of numerical precision for lmax > 28. "
               "If you need more bands, try recompiling using doubles.");
        exit(1);
    }

    // Compute Legendre polynomial values for $\cos\theta$
    Assert(w.Length() > .995f && w.Length() < 1.005f);
    legendrep(w.z, lmax, out);

    // Compute $K_l^m$ coefficients
    float *Klm = ALLOCA(float, SHTerms(lmax));
    for (int l = 0; l <= lmax; ++l)
        for (int m = -l; m <= l; ++m)
            Klm[SHIndex(l, m)] = K(l, m);

    // Compute $\sin\phi$ and $\cos\phi$ values
    float *sins = ALLOCA(float, lmax+1), *coss = ALLOCA(float, lmax+1);
    float xyLen = sqrtf(max(0.f, 1.f - w.z*w.z));
    if (xyLen == 0.f) {
        for (int i = 0; i <= lmax; ++i) sins[i] = 0.f;
        for (int i = 0; i <= lmax; ++i) coss[i] = 1.f;
    }
    else
        sinCosIndexed(w.y / xyLen, w.x / xyLen, lmax+1, sins, coss);

    // Apply SH definitions to compute final $(l,m)$ values
    static const float sqrt2 = sqrtf(2.f);
    for (int l = 0; l <= lmax; ++l) {
        for (int m = -l; m < 0; ++m)
        {
            out[SHIndex(l, m)] = sqrt2 * Klm[SHIndex(l, m)] *
                out[SHIndex(l, -m)] * sins[-m];
            Assert(!isnan(out[SHIndex(l,m)]));
            Assert(!isinf(out[SHIndex(l,m)]));
        }
        out[SHIndex(l, 0)] *= Klm[SHIndex(l, 0)];
        for (int m = 1; m <= l; ++m)
        {
            out[SHIndex(l, m)] *= sqrt2 * Klm[SHIndex(l, m)] * coss[m];
            Assert(!isnan(out[SHIndex(l,m)]));
            Assert(!isinf(out[SHIndex(l,m)]));
        }
    }
}

作用：

（这里是 通过一个向量w 和 a maximum number of SH bands lmax 来 计算 Ylm，并保存 Ylm 到 out 参数中，Ylm 其实就是

《PBRT_V2 总结记录 <109> Spherical Harmonic (SH) Basis Functions》中的17.1 公式）

Now we can define the SH evaluation function, SHEvaluate(). SHEvaluate() takes a
direction vector w and a maximum number of SH bands lmax. It evaluates the real SH
functions up to l = lmax, storing the values of the functions in the provided out array,
using the indexing scheme defined by SHIndex(). The length of out must be at least
SHTerms(lmax).

 

细节

a. 

 legendrep(w.z, lmax, out);

作用：

（这里计算 ，把计算好的一系列的  都存到 out中，后面要使用）

The implementation starts by computing the associated Legendre polynomial values  using the legendrep() utility routine.

It stores the results in the out array;
after these values have been stored, the implementation below will multiply by the remaining
factors from Equation (17.7). Recall from Section 5.5.2 that in spherical coordinates,
cos θ = z. Thus, we can just pass the z component of w in to the legendrep()
function.

 

b.

static void legendrep(float x, int lmax, float *out) {
#define P(l,m) out[SHIndex(l,m)]
    // Compute $m=0$ Legendre values using recurrence
    P(0,0) = 1.f;
    P(1,0) = x;
    for (int l = 2; l <= lmax; ++l)
    {
        P(l, 0) = ((2*l-1)*x*P(l-1,0) - (l-1)*P(l-2,0)) / l;
        Assert(!isnan(P(l, 0)));
        Assert(!isinf(P(l, 0)));
    }

    // Compute $m=l$ edge using Legendre recurrence
    float neg = -1.f;
    float dfact = 1.f;
    float xroot = sqrtf(max(0.f, 1.f - x*x));
    float xpow = xroot;
    for (int l = 1; l <= lmax; ++l) {
        P(l, l) = neg * dfact * xpow;
        Assert(!isnan(P(l, l)));
        Assert(!isinf(P(l, l)));
        neg *= -1.f;      // neg = (-1)^l
        dfact *= 2*l + 1; // dfact = (2*l-1)!!
        xpow *= xroot;    // xpow = powf(1.f - x*x, float(l) * 0.5f);
    }

    // Compute $m=l-1$ edge using Legendre recurrence
    for (int l = 2; l <= lmax; ++l)
    {
        P(l, l-1) = x * (2*l-1) * P(l-1, l-1);
        Assert(!isnan(P(l, l-1)));
        Assert(!isinf(P(l, l-1)));
    }

    // Compute $m=1, \ldots, l-2$ values using Legendre recurrence
    for (int l = 3; l <= lmax; ++l)
        for (int m = 1; m <= l-2; ++m)
        {
            P(l, m) = ((2 * (l-1) + 1) * x * P(l-1,m) -
                       (l-1+m) * P(l-2,m)) / (l - m);
            Assert(!isnan(P(l, m)));
            Assert(!isinf(P(l, m)));
        }
    
    
#undef P
}

作用：

There are four main steps in the implementation of legendrep(); each one uses a different
recurrence to compute a different subset of the associated Legendre polynomial values for
the given value. Figure 17.4 illustrates which subset of (l , m) values each step computes.
The code here fills in the given out array where index SHIndex(l,m) holds the value of
the associated Legendre polynomial for the given x value. It only computes the values of

for m ≥ 0; the code in SHEvaluate() doesn’t need the values for m<0 in the end,
so they aren’t computed here

（不需要 m < 0, 参考 《PBRT_V2 总结记录 <109> Spherical Harmonic (SH) Basis Functions》有说明）

思路：

 

Figure 17.4: The Steps of Computing the Values of the Associated Legendre Polynomials. A
series of recurrences is used to compute their values in the function legendrep(). At each step, the
newly computed values have dark shading and the previously computed values have light shading.
(a) The values along m = 0 are computed for all l bands, (b) the m = l edge is computed from l = 1 to
l = lmax, (c) the m = l − 1 edge is computed for l = 2 to l = lmax, (d) finally, remaining points m = 1 to
m = l − 2 are computed for bands l = 3 and higher.

（a 图 计算 m = 0 ，L = 0~max 的 情况，b 图  m = L， L = 1~max 的情况， c 图 m = L - 1, L = 2~max 的情况， d 图 m = 1~L-2, L = 3 的情况，这样就可以把所有的P值计算出来，）

 

a图：

    // Compute $m=0$ Legendre values using recurrence
    P(0,0) = 1.f;
    P(1,0) = x;
    for (int l = 2; l <= lmax; ++l)
    {
        P(l, 0) = ((2*l-1)*x*P(l-1,0) - (l-1)*P(l-2,0)) / l;
        Assert(!isnan(P(l, 0)));
        Assert(!isinf(P(l, 0)));
    }

The m = 0，The first two of them are 

Subsequent ones are given by the recurrence

 

b 图：

// Compute $m=l$ edge using Legendre recurrence
    float neg = -1.f;
    float dfact = 1.f;
    float xroot = sqrtf(max(0.f, 1.f - x*x));
    float xpow = xroot;
    for (int l = 1; l <= lmax; ++l) {
        P(l, l) = neg * dfact * xpow;
        Assert(!isnan(P(l, l)));
        Assert(!isinf(P(l, l)));
        neg *= -1.f;      // neg = (-1)^l
        dfact *= 2*l + 1; // dfact = (2*l-1)!!
        xpow *= xroot;    // xpow = powf(1.f - x*x, float(l) * 0.5f);
    }

where m = l. For these values we use following recurrence:

In this expression, !! denotes the double factorial, which is defined as

 

c 图：

// Compute $m=l-1$ edge using Legendre recurrence
    for (int l = 2; l <= lmax; ++l)
    {
        P(l, l-1) = x * (2*l-1) * P(l-1, l-1);
        Assert(!isnan(P(l, l-1)));
        Assert(!isinf(P(l, l-1)));
    }

作用：

Next, the slice where m = l − 1 is computed. These values are found using the recurrence

which is equivalent to

 

d 图：

    for (int l = 3; l <= lmax; ++l)
        for (int m = 1; m <= l-2; ++m)
        {
            P(l, m) = ((2 * (l-1) + 1) * x * P(l-1,m) -
                       (l-1+m) * P(l-2,m)) / (l - m);
            Assert(!isnan(P(l, m)));
            Assert(!isinf(P(l, m)));
        }

作用：

（ L 偏移 -1 就可以得到 上面代码中的公式）

We can now finally fill in the remainder of the values for 0 < m < l − 1 using one last
recurrence,

Once again we offset l by −1 and divide the right-hand side by the resulting (l − m)
factor in the implementation below to compute the value of 

Note that the loop
starting and ending points are chosen carefully so that they don’t recompute values that
the previous recurrences have already covered.

 

细节

c

    // Compute $K_l^m$ coefficients
    float *Klm = ALLOCA(float, SHTerms(lmax));
    for (int l = 0; l <= lmax; ++l)
        for (int m = -l; m <= l; ++m)
            Klm[SHIndex(l, m)] = K(l, m);





static inline float K(int l, int m) {
    return sqrtf((2.f * l + 1.f) * INV_FOURPI * divfact(l, m));
}

static inline float divfact(int a, int b) {
    if (b == 0) return 1.f;
    float fa = a, fb = fabsf(b);
    float v = 1.f;
    for (float x = fa-fb+1.f; x <= fa+fb; x += 1.f)
        v *= x;
    return 1.f / v;
}

作用：

（计算 ）

Next, the SHEvaluate() function needs to compute the  coefficients using Equation
(17.6).

 

细节

d.

    // Compute $\sin\phi$ and $\cos\phi$ values
    float *sins = ALLOCA(float, lmax+1), *coss = ALLOCA(float, lmax+1);
    float xyLen = sqrtf(max(0.f, 1.f - w.z*w.z));
    if (xyLen == 0.f) {
        for (int i = 0; i <= lmax; ++i) sins[i] = 0.f;
        for (int i = 0; i <= lmax; ++i) coss[i] = 1.f;
    }
    else
        sinCosIndexed(w.y / xyLen, w.x / xyLen, lmax+1, sins, coss);






static void sinCosIndexed(float s, float c, int n,
                          float *sout, float *cout) {
    float si = 0, ci = 1;
    for (int i = 0; i < n; ++i) {
        // Compute $\sin{}i\phi$ and $\cos{}i\phi$ using recurrence
        *sout++ = si;
        *cout++ = ci;
        float oldsi = si;
        si = si * c + ci * s;
        ci = ci * c - oldsi * s;
    }
}

作用：

SHEvaluate() next computes the sine and cosine factors of Equation (17.7). From this
equation, the real spherical harmonics have factors of the form sin(−mφ) and cos(mφ);
thus, we need the values of sin −φ, sin(−2φ), and so forth up to sin(−mφ) for m = lmax
and similarly for the cosine factors.

（看 Equation (17.7)会发现，有会 sin(−mφ) and cos(mφ);也就是说，我们需要计算 例如：

sin −φ, sin(−2φ),sin(−3φ)，sin(−4φ).....,对于 cos 也是一样的道理）

 

（，

）

sins 和coss 会保存一系列的 sin(mφ)  和 cos(mφ)，

 

There are recurrence formulas that can be used to give the values of sin(iφ) and cos(iφ)
in terms of sin((i − 1)φ) and cos((i − 1)φ):

（这里 其实是可以 表示为

sin( iφ）= sin( (i-1)φ + φ ) = sin( (i-1)φ ) cos φ + cos( (i-1)φ ) sin (φ)

cos(iφ) = cos( (i-1)φ + φ ) = cos( (i-1)φ ) cos( φ) - sin ( (i-1)φ ) sin (φ)

 

细节

d

    // Apply SH definitions to compute final $(l,m)$ values
    static const float sqrt2 = sqrtf(2.f);
    for (int l = 0; l <= lmax; ++l) {
        for (int m = -l; m < 0; ++m)
        {
            out[SHIndex(l, m)] = sqrt2 * Klm[SHIndex(l, m)] *
                out[SHIndex(l, -m)] * sins[-m];
            Assert(!isnan(out[SHIndex(l,m)]));
            Assert(!isinf(out[SHIndex(l,m)]));
        }
        out[SHIndex(l, 0)] *= Klm[SHIndex(l, 0)];
        for (int m = 1; m <= l; ++m)
        {
            out[SHIndex(l, m)] *= sqrt2 * Klm[SHIndex(l, m)] * coss[m];
            Assert(!isnan(out[SHIndex(l,m)]));
            Assert(!isinf(out[SHIndex(l,m)]));
        }
    }

作用：

（整合 17.7 公式）

Given all of these components, the SHEvaluate() function can now finish computing the
values from Equation (17.7). Going into this fragment, the out array holds the values of
the associated Legendre polynomials for the corresponding (l , m) indices. For each term,
then, we need to multiply in the normalization factor and the sine or cosine factor.
The implementation below uses the fact that sin −x =−sin x. For them<0 terms, note
that the associated Legendre polynomial values for m>0 are used as per the definition
in Equation (17.7).