prt1 | tinf

prt1：precompute part Diffuse Case

Outline

渲染方程（漫反射）的积分项可以拆成"光"和"传输"两部分，各自投影到球谐基（SH），渲染时只做点积：

Precompute lighting and light transport for each individual shading point.

L_{o} (x) = \frac{ρ}{π} \int_{Ω} L_{i} (ω) T (ω) d ω \overset{SH 正交性}{\to} \sum_{l, m} L_{l m} T_{l m}

光系数 $L_{l m}$ ：环境贴图投影到 SH（全局一份）。
传输系数 $T_{l m}$ ：每个顶点的 $T (ω) = max (N \cdot ω, 0) \cdot V (ω)$ 投影到 SH（每顶点一份）。

系数就是"函数在基函数轴上的投影"：

c_{l m} = \int_{Ω} f (ω) Y_{l m} (ω) d ω

系数与投影与点积

PrecomputeCubemapSH光系数

离散

积分无法解析计算，逐像素近似求和（cubemap 的每个像素贡献一份）：

c_{l m} \approx \sum_{px} f (ω_{p}) Y_{l m} (ω_{p}) Δ ω_{p}

量	对应
$f (ω_{p})$ （光函数值）	`Le`（该像素 RGB）
$Y_{l m} (ω_{p})$ （基函数值）	`sh::EvalSH(l, m, dirD)`
$Δ ω_{p}$ （立体角）	`CalcArea(x, y, width, height)`

像素方向

每个面的方向由三个轴向量拼出（cubemapFaceDirections）：

cpp

Eigen::Vector3f dir = (faceDirX * u + faceDirY * v + faceDirZ).normalized();
// u, v ∈ [-1, 1]，取自像素中心 (x + 0.5) / width

求和sum

cpp

for (int i = 0; i < 6; i++)                 // 6 个面
    for (int y = 0; y < height; y++)
        for (int x = 0; x < width; x++) {
            Eigen::Vector3f dir = cubemapDirs[i * width * height + y * width + x];
            int index = (y * width + x) * channel;
            Eigen::Array3f Le(images[i][index + 0], images[i][index + 1],
                              images[i][index + 2]);

            float angle = CalcArea(x, y, width, height);          // 注意签名 (u_, v_, width, height)...某种比较经典的求球面积分图...？
            Eigen::Vector3d dirD = dir.cast<double>();            // EvalSH 要 double

            for (int l = 0; l <= SHOrder; ++l)                    
                for (int m = -l; m <= l; ++m)                     
                    SHCoeffiecents[sh::GetIndex(l, m)] +=
                        Le * sh::EvalSH(l, m, dirD) * angle;      // f × Y × dω
        }

GetIndex(l, m) = l*(l+1)+m
EvalSH 需要 Eigen::Vector3d，dir 是 Vector3f，要 cast<double>()。
2 阶 SH → 9 个系数，GetIndex 下标范围 [0, 8]。

EvalSH

cpp

double EvalSH(int l, int m, const Eigen::Vector3d& dir) {
  if (l <= kHardCodedOrderLimit) {
    // Validate l and m here (don't do it generally since EvalSHSlow also
    // checks it if we delegate to that function).
    CHECK(l >= 0, "l must be at least 0.");
    CHECK(-l <= m && m <= l, "m must be between -l and l.");
    CHECK(NearByMargin(dir.squaredNorm(), 1.0), "dir is not unit.");

    switch (l) {
      case 0:
        return HardcodedSH00(dir);
      case 1:
        switch (m) {
          case -1:
            return HardcodedSH1n1(dir);
          case 0:
            return HardcodedSH10(dir);
          case 1:
            return HardcodedSH1p1(dir);
        }
      case 2:
        switch (m) {
          case -2:
            return HardcodedSH2n2(dir);
          case -1:
            return HardcodedSH2n1(dir);
          case 0:
            return HardcodedSH20(dir);
          case 1:
            return HardcodedSH2p1(dir);
          case 2:
            return HardcodedSH2p2(dir);
        }
      case 3:
        switch (m) {
          case -3:
            return HardcodedSH3n3(dir);
          case -2:
            return HardcodedSH3n2(dir);
          case -1:
            return HardcodedSH3n1(dir);
          case 0:
            return HardcodedSH30(dir);
          case 1:
            return HardcodedSH3p1(dir);
          case 2:
            return HardcodedSH3p2(dir);
          case 3:
            return HardcodedSH3p3(dir);
        }
      case 4:
        switch (m) {
          case -4:
            return HardcodedSH4n4(dir);
          case -3:
            return HardcodedSH4n3(dir);
          case -2:
            return HardcodedSH4n2(dir);
          case -1:
            return HardcodedSH4n1(dir);
          case 0:
            return HardcodedSH40(dir);
          case 1:
            return HardcodedSH4p1(dir);
          case 2:
            return HardcodedSH4p2(dir);
          case 3:
            return HardcodedSH4p3(dir);
          case 4:
            return HardcodedSH4p4(dir);
        }
    }

    // This is unreachable given the CHECK's above but the compiler can't tell.
    return 0.0;
  } else {
    // Not hard-coded so use the recurrence relation (which will convert this
    // to spherical coordinates).
    return EvalSHSlow(l, m, dir);
  }
}

Filter

和"邻域像素取平均再存回"的空间域模糊不同(split sum等...？)。
是投影到 SH 基 + 截断到 2 阶：高频细节全部丢掉，等价于低通滤波。
依据（Ramamoorthi & Hanrahan 2001）：漫反射 BRDF 本身就是[低通滤波器]，光照与 cos 核卷积后只含低频，9 个系数近似误差约 1%。

shFunc + ProjectFunction物体系数

outline

shFunc(phi, theta)：给定一个方向，返回一个数——传输函数值 $T (ω)$ 。
ProjectFunction(order, func, sample_count)：内部做分层均匀球面采样 + 乘基函数累加 + 乘权重 $4 π / N$ ，返回系数向量。

投影（上面写过的循环等）在 ProjectFunction 里完成。

ProjectFunction

cpp

std::unique_ptr<std::vector<double>> ProjectFunction(
    int order, const SphericalFunction& func, int sample_count) {
  CHECK(order >= 0, "Order must be at least zero.");
  CHECK(sample_count > 0, "Sample count must be at least one.");

  // This is the approach demonstrated in [1] and is useful for arbitrary
  // functions on the sphere that are represented analytically.
  const int sample_side = static_cast<int>(floor(sqrt(sample_count)));
  std::unique_ptr<std::vector<double>> coeffs(new std::vector<double>());
  coeffs->assign(GetCoefficientCount(order), 0.0);

  // generate sample_side^2 uniformly and stratified samples over the sphere
  std::random_device rd;
  std::mt19937 gen(rd());
  std::uniform_real_distribution<> rng(0.0, 1.0);
  for (int t = 0; t < sample_side; t++) {
    for (int p = 0; p < sample_side; p++) {
      double alpha = (t + rng(gen)) / sample_side;
      double beta = (p + rng(gen)) / sample_side;
      // See http://www.bogotobogo.com/Algorithms/uniform_distribution_sphere.php
      double phi = 2.0 * M_PI * beta;
      double theta = acos(2.0 * alpha - 1.0);

      // evaluate the analytic function for the current spherical coords
      double func_value = func(phi, theta);

      // evaluate the SH basis functions up to band O, scale them by the
      // function's value and accumulate them over all generated samples
      for (int l = 0; l <= order; l++) {
        for (int m = -l; m <= l; m++) {
          double sh = EvalSH(l, m, phi, theta);
          (*coeffs)[GetIndex(l, m)] += func_value * sh;
        }
      }
    }
  }

  // scale by the probability of a particular sample, which is
  // 4pi/sample_side^2. 4pi for the surface area of a unit sphere, and
  // 1/sample_side^2 for the number of samples drawn uniformly.
  double weight = 4.0 * M_PI / (sample_side * sample_side);
  for (unsigned int i = 0; i < coeffs->size(); i++) {
     (*coeffs)[i] *= weight;
  }

  return coeffs;
}

unshadowed

T (ω) = max (N \cdot ω, 0)

cpp

auto shFunc = [&](double phi, double theta) -> double {
    Eigen::Array3d d = sh::ToVector(phi, theta);
    const auto wi = Vector3f(d.x(), d.y(), d.z());
    if (m_Type == Type::Unshadowed)
        return std::max(0.0f, n.dot(wi));     // 就是 H，一个数
    ...
};

shadowed：可见性 = 一次射线求交

T (ω) = max (N \cdot ω, 0) \cdot V (ω)

cpp

else  // Shadowed
{
    float H = n.dot(wi);
    if (H < 0.0) return 0.0f;               // 下半球直接 0，省一次求交

    Intersection its;
    Ray3f ray(v, wi);                       // 默认 mint = Epsilon，自动避免自相交
    if (scene->rayIntersect(ray, its))
        return 0.0f;                        // 有遮挡 → V = 0
    return H;                               // 无遮挡 → V = 1
}

scene->rayIntersect(const Ray3f&, Intersection&) const：命中返回 true。
Ray3f(o, d) 默认 mint = Epsilon（ray.h），不需要手动偏移起点。
先判 H 后打射线，省一半求交；m_SampleCount（默认 100）控制精度，每顶点要打几百根射线，这就是预计算的耗时来源。

ProjectFunction 内部的权重

cpp

auto shCoeff = sh::ProjectFunction(SHOrder, shFunc, m_SampleCount);
for (int j = 0; j < shCoeff->size(); j++)
    m_TransportSHCoeffs.col(i).coeffRef(j) = (*shCoeff)[j];

采样数会被取成不大于它的最大完全平方数。
权重 $4 π / N$ 来自均匀球面采样的概率密度（球面积 $4 π$ ），在库内部已乘好。

关于讲义里的 rgb offset

讲义伪代码：

result[j + red_offset]   += value;
result[j + green_offset] += value;
result[j + blue_offset]  += value;

传输项是"无色"标量，光系数是 RGB 三通道 → 把同一标量复制三份对齐 RGB。
本框架不需要：m_TransportSHCoeffs 是 SHCoeffLength × vertexCount 的标量矩阵，Li() 里用同一个向量分别和 rL/gL/bL 点积，天然对齐。

（Li）

cpp

Color3f c0 = Color3f(rL.dot(sh0), gL.dot(sh0), bL.dot(sh0));   // 光系数 · 传输系数

result

! shadowed

unshadowed

Cost

L_{o} (p, ω_{o}) = \int_{Ω^{+}} L_{i} (p, ω_{i}) f_{r} (p, ω_{i}, ω_{o}) \cos θ_{i} V (p, ω_{i}) d ω_{i}

= \sum_{p} \sum_{q} c_{p} c_{q} \int_{Ω^{+}} B_{p} (ω_{i}) B_{q} (ω_{i}) d ω_{i}

L (ω_{i}) \approx \sum_{p} c_{p} B_{p} (ω_{i})

T (ω_{i}) \approx \sum_{q} c_{q} B_{q} (ω_{i})

the complexity is still O（n）.

Glossy Case

The Diffuse Case

Outline ​

PrecomputeCubemapSH光系数 ​

离散 ​

像素方向 ​

求和sum ​

Filter ​

shFunc + ProjectFunction物体系数 ​

outline ​

unshadowed ​

shadowed：可见性 = 一次射线求交 ​

ProjectFunction 内部的权重 ​

关于讲义里的 rgb offset ​

（Li） ​

result ​

Cost ​

Glossy Case ​

Outline

PrecomputeCubemapSH光系数

离散

像素方向

求和sum

Filter

shFunc + ProjectFunction物体系数

outline

unshadowed

shadowed：可见性 = 一次射线求交

ProjectFunction 内部的权重

关于讲义里的 rgb offset

（Li）

result

Cost

Glossy Case