这个系列是Physically Based Rendering: From Theory To Implementation的读书笔记，本节主要介绍反射模型。
Basic Terminology
我们能够看到物体的本质原因是物体反射了光，因此反射模型对于渲染出正确的图像至关重要。在Surface Reflection一节中我们介绍过BRDF的概念，它表示在$\omega_o$方向上反射光与$\omega_i$方向上接收光的比值。除了BRDF外，对于透明材质光线还会出现透射的情况，描述这种现象的函数称为BTDF。在PBRT中我们使用散射(scattering)将反射和透射统一起来，这样描述光线和材料表面相互作用的函数称为BSDF。
Reflection Categories
首先来考虑反射的情况。在PBRT中我们将反射分成4种不同的类型：漫反射(diffuse)、光泽镜面反射(glossy specular)、完美镜面反射(perfect specular)以及反向反射(retro-reflective)。其中漫反射表示物体接收光线后再各个方向上都具有相同的反射；光泽镜面反射表示光线的反射会聚集在一定的范围上；完美镜面反射则更加严格，它表示物体反射的光线只会出现在一个确定的方向上；而反向反射则表示反射的方向与入射方向在一个区域内。现实中的反射现象可以看作以上四种反射类型的组合，它们的示意图如下：
同时还需要注意的是反射材料的方向性。现实生活中大部分材料是各向同性(isotropic)的，对于各向同性材料反射只和入射光线与法线的夹角有关，当入射光线绕法线旋转时反射的能量不变；而对于各项异性(anisotropic)材料，当入射光线绕法线旋转时反射的能量也会同时发生变化，在渲染时要格外注意。生活中常见的各项异性材料包括拉丝后的金属(brushed metal)、布料等等。
Geometric Setting
PBRT在考虑反射时会在光线与材料交点的局部坐标系中进行计算，这个局部坐标系称为着色坐标系(shading coordinate system)。这样在进行计算时首先会把入射光线放在着色坐标系中并计算反射光线，然后再把反射光线转换到世界坐标系中进行后续的计算。
着色坐标系可以通过三个正交的向量$(\mathbf{s}, \mathbf{t}, \mathbf{n})$来定义，其中$\mathbf{n}$为交点处的法线方向，另外两个方向则是局部坐标系的x轴和y轴。
在计算方向时我们往往会使用球坐标来代替直角坐标系，此时球面上的方向可以使用两个角度$(\theta, \phi)$来表示。其中$\theta$表示该方向与z轴的夹角，$\phi$则是方向投影后和x轴的夹角。
还需要注意的是PBRT在计算反射时有一些假定：
- 入射光线$\omega_i$和反射光线$\omega_o$在计算时都需要进行单位化，而且它们都射向光线交点外；
- PBRT中默认法向$\mathbf{n}$指向物体外侧，这样便于判断光线是射入还是射出透射的物体：光线方向与$\mathbf{n}$在同一半球时表示射入，在相反的半球上时表示射出；
- 求交计算出的坐标系与着色坐标系可能是不同的，在计算凹凸贴图时会修改求交计算的坐标系；
- BRDF和BTDF不会假定入射光线$\omega_i$和反射光线$\omega_o$在同一半球上。
Specular Reflection and Transmission
对于完美镜面反射的情况我们可以通过几何光学来描述入射光线和出射光线的关系：
\[\theta_o = \theta_i\] \[\phi_o = \phi_i + \pi\]对于透射的情况，$\phi_t$与$\phi_i$的关系与反射一致，但$\theta_t$则需要满足折射定律(Snell’s law)：
\[\eta_i \sin \theta_i = \eta_t \sin \theta_t\]其中$\eta$为折射率(refractive index)。在通常情况下不同波长的光线具有不同的折射率，这就导致入射光线发生折射时会散射到不同的方向上，这种现象称为色散(dispersion)。不过在图形学中一般会忽略这种现象，在大多数情况下我们会直接假定不同波长的光线都具有相同的折射率。
Fresnel Reflectance
除了方向外我们还需要计算反射和透射光线的能量占入射光线能量的比例。对于反射的情况这个比例由Fresnel方程来描述，它是Maxwell方程组在光滑表面上的解。对于给定的折射率和入射角度，Fresnel方程的解给出了入射光线在两个极化方向上入射光的反射率。在大多数情况下光线的极性不会产生人眼可见的视觉效果，因此在PBRT中我们忽略光的极性直接假设光线不具有极性，这样Fresnel反射可以认为是垂直和平行方向反射率的平方平均值。
在计算Fresnsl方程时需要考虑不同的材料：
- 电介质(dielectrics)：电介质是不导电的材料，它们的折射率为实数(一般介于1~3之间)
- 导体(conductors)：导体可以导电，它们的折射率为复数$\bar{\eta} = \eta + \text{i} k$
- 半导体(Semiconductors)：在PBRT中不考虑半导体这种材质
对于实数折射率的情况，反射率计算公式为：
\[r_\parallel = \frac{\eta_t \cos \theta_i - \eta_i \cos \theta_t}{\eta_t \cos \theta_i + \eta_i \cos \theta_t}\] \[r_\perp = \frac{\eta_i \cos \theta_i - \eta_t \cos \theta_t}{\eta_i \cos \theta_i + \eta_t \cos \theta_t}\]$r_\parallel$和$r_\perp$为光线在两个极化方向上的反射率。对于非极性的光线，Fresnel反射率计算公式为：
\[F_r = \frac{1}{2} (r_\parallel^2 + r_\perp^2)\]根据能量守恒，电介质透射的能量比例为$1 - F_r$。对于透射的情况还需要计算透射的角度，这可以通过两种介质中的折射率来计算。
还需要注意的是当光线从高折射率射向低折射率的介质时，根据折射定律可能会出现$\sin \theta_t \gt 1$的情况。发生这种现象最大的入射角度称为临界角度(critical angle)，当入射角度$\theta_i$大于临界角度时所有的入射能量会反射出去，因此这样的现象也称为全内反射(total internal reflection)。显然，当发生全内反射时就不再需要求解Fresnel方程来计算反射率了。
对于折射率是复数的情况，一部分入射能量会被材料吸收并转化成热能。此时Fresnel方程会依赖于折射率的虚部$k$，它也称为材料的吸收系数(absorption coefficient)。一般情况下我们不需要考虑导体的折射情况，因此对于电介质与导体的交界处的反射有
\[r_\perp = \frac{a^2 + b^2 - 2a \cos \theta + \cos^2 \theta}{a^2 + b^2 + 2a \cos \theta + \cos^2 \theta}\] \[r_\parallel = r_\perp \frac{\cos^2 \theta \ (a^2 + b^2) - 2 a \cos \theta \sin^2 \theta + \sin^4 \theta}{\cos^2 \theta \ (a^2 + b^2) + 2 a \cos \theta \sin^2 \theta + \sin^4 \theta}\] \[a^2 + b^2 = \sqrt{(\eta^2 - k^2 - \sin^2 \theta)^2 + 4 \eta^2 k^2}\]其中$\eta + \text{i} k = \bar{\eta}_t / \bar{\eta}_i$为两种介质的相对折射率。
Specular Reflection
在此基础上我们就可以推导出镜面反射的计算公式：
\[L_o(\omega_o) = \int f_r(\omega_o, \omega_i) L_i(\omega_i) \vert \cos \theta \vert \ d \omega = F_r(\omega_r) L_i(\omega_r)\]其中$\omega_r = R(\omega_o, \mathbf{n})$为$\omega_o$经过法线为$\mathbf{n}$的平面反射后的方向。我们可以利用$\delta(x)$函数来构造出镜面反射的BRDF：
\[f_r(\omega_o, \omega_i) = F_r(\omega_r) \frac{\delta (\omega_i - \omega_r)}{\vert \cos \theta_r \vert}\]Specular Transmission
接下来我们通过折射定律来推导镜面折射情况下的BTDF。假设光线在两种材料的边界上发生折射，材料的折射率分别为$\eta_i$和$\eta_o$，折射能量的比例为$\tau = 1 - F_r(\omega_i)$。因此可以得到入射和折射能量之间的关系：
\[d \Phi_o = \tau d \Phi_i\]利用能量通量与radiance的关系可以得到：
\[L_o \cos \theta_o \ d A \ d \omega_o = \tau (L_i \cos \theta_i \ d A \ d \omega_i)\]将微分立体角展开得到：
\[L_o \cos \theta_o \ d A \ \sin \theta_o \ d \theta_o \ d \phi_o = \tau (L_i \cos \theta_i \ d A \ \sin \theta_i \ d \theta_i \ d \phi_i)\]
对折射定律两边同时进行微分可以得到：
\[\eta_o \cos \theta_o \ d \theta_o = \eta_i \cos \theta_i \ d \theta_i\] \[\frac{\cos \theta_o \ d \theta_o}{\cos \theta_i \ d \theta_i} = \frac{\eta_i}{\eta_o}\]整理上面的式子得到：
\[L_o \eta_i^2 \ d \phi_o = \tau L_i \eta_o^2 \ d \phi_i\]由于$\phi_o = \phi_i + \pi$，因此得到$L_o$和$L_i$的关系：
\[L_o = \tau L_i \frac{\eta_o^2}{\eta_i^2}\]这样就得到了镜面折射情况下的BTDF：
\[f_r(\omega_o, \omega_i) = \frac{\eta_o^2}{\eta_i^2} (1 - F_r(\omega_i)) \frac{\delta (\omega_i - T(\omega_o, \mathbf{n}) )}{\vert \cos \theta_i \vert}\]其中$T(\omega_o, \mathbf{n})$为$\omega_o$在法线为$\mathbf{n}$的平面上折射后的方向。
Lambertian Reflection
Lambertian反射模型是最简单的BRDF。对于漫反射材质，Lambertian模型假设各个方向上反射的能量是相等的。根据能量守恒可以得到Lambertian模型的BRDF：
\[f_r(\omega_o, \omega_i) = \frac{R}{\pi}\]其中$R$表示入射能量被反射的比例。
Microfacet Models
微表面模型(microfacet models)在基于物理的材质建模中有着大量的应用。它的基本思想是把粗糙材料的表面看作是各个方向光滑平面的集合，这样尽管在宏观尺度上材料表面的法向是唯一的，但在微观尺度上则对应一个法向的分布。对于越粗糙的表面法向分布变化越剧烈，越光滑的表面法向分布的变化越平缓。
微表面模型的BRDF通过对法向分布的建模来描述光线在材料表面发生散射的现象。假设面积微元$d A$和单个微表面的面积相比足够大，这样材质的BRDF就可以看做是这些微表面相互作用的平均。
显然微表面模型的BRDF有两个重要组成部分：微表面的分布以及光线和单个微表面的BRDF。对于第二项，在大多数微表面模型中都假设单个微表面满足光滑镜面反射。当然也有一些模型会考虑镜面的透射，而后面介绍的Oren–Nayar模型则是使用Lambertian模型来考虑微表面反射。
另一方面微表面模型需要考虑不同方向微表面之间的相对关系，根据入射光线和观察的方向微表面之间有三种几何关系：masking、shadowing以及interreflection。masking是指观测方向会被其它微表面遮挡，shadowing是指观测方向上的微表面位于其它微表面的阴影中，而interreflection则是说光线经过反射后出射方向与入射方向相同。
Oren–Nayar Diffuse Reflection
Oren和Nayar发现真实世界中的材质不会遵循Lambertian模型。特别地，对于粗糙表面的材质当观测方向接近入射方向时物体会比Lambertian模型显得更亮一些。在此基础上他们开发了一种基于V形微表面分布的材质模型，该模型包含一个参数$\sigma$表示微表面方向分布的标准差。Oren–Nayar模型的BRDF为：
\[f_r(\omega_i, \omega_o) = \frac{R}{\pi} \bigg(A + B \max \big(0, \cos (\phi_i - \phi_o) \big) \sin \alpha \tan \beta \bigg)\] \[A = 1 - \frac{\sigma^2}{2 (\sigma^2 + 0.33)}\] \[B = \frac{0.45 \sigma^2}{\sigma^2 + 0.09}\] \[\alpha = \max (\theta_i, \theta_o)\] \[\beta = \min (\theta_i, \theta_o)\]Microfacet Distribution Functions
在微表面模型中法向分布函数(normal distribution function)$D(\omega_h)$描述了法向$\omega_h$在微表面处对应的面积微元。在PBRT中使用着色坐标来计算法向分布函数，它需要满足归一化条件，即各个方向投影后的面积微元积分为$d A$：
\[\int_{H^2(\mathbf{n})} D(\omega_h) \cos \theta_h \ d \omega_h = 1\]
Beckmann–Spizzichino模型将NDF定义为：
\[D(\omega_h) = \frac{e^{-\tan^2 \theta_h / \alpha^2}}{\pi \alpha^2 \cos^4 \theta_h}\]其中$\sigma$为微表面的RMS slope，$\alpha = \sqrt{2} \sigma$。
对于各向异性材料，它的NDF定义为：
\[D(\omega_h) = \frac{e^{-\tan^2 \theta_h (\cos^2 \phi_h / \alpha_x^2 + \sin^2 \phi_h / \alpha_y^2)}}{\pi \alpha_x \alpha_y \cos^4 \theta_h}\]其中$\alpha_x$和$\alpha_y$分别为x轴和y轴上的分布参数。同时不难发现各向同性材料的NDF相当于$\alpha_x = \alpha_y$情况下的特例。
另一种常用的NDF模型是Trowbridge–Reitz模型，它的定义为：
\[D(\omega_h) = \frac{1}{\pi \alpha_x \alpha_y \cos^4 \theta_h \big( 1 + \tan^2 \theta_h (\cos^2 \phi_h / \alpha_x^2 + \sin^2 \phi_h / \alpha_y^2) \big)^2}\]和Beckmann–Spizzichino模型相比Trowbridge–Reitz模型是一个长尾分布，也更接近于现实生活中材料表面的性质。这两种模型的对比可见下图：
Masking and Shadowing
仅有NDF并不能够完整描述微表面的几何性质，对于masking和shadowing的问题我们需要一个几何函数$G_1(\omega, \omega_h)$来进行描述，称为Smith’s masking-shadowing function。$G_1(\omega, \omega_h)$表示朝向为$\omega_h$的微表面从$\omega$方向进行观察时的可见性，在大多数情况下我们可以假定各个方向的微表面可见性与它们的朝向无关，因此可以将它简记为$G_1(\omega)$。
除此之外根据可见微表面的面积微元与$d A$之间的几何关系，$G_1$需要满足约束条件：
\[\cos \theta = \int_{H^2(\mathbf{n})} G_1(\omega, \omega_h) \max (0, \omega \cdot \omega_h) D(\omega_h) \ d \omega_h\]
由于微表面可以理解为一个高度场，我们可以按照观测的方向将微表面分为面向$\omega$的部分和背对$\omega$的部分，它们的投影面积分别为$A^+$和$A^-$。结合上式有：
\[\cos \theta = A^+ (\omega) - A^- (\omega)\] \[G_1 (\omega) = \frac{A^+ (\omega) - A^- (\omega)}{A^+ (\omega)}\]通常情况下masking-shadowing function还可以通过一个辅助函数$\Lambda (\omega)$来表示，它定义为被遮挡的面积占可见面积的比例：
\[\Lambda (\omega) = \frac{A^- (\omega)}{A^+ (\omega) - A^- (\omega)} = \frac{A^- (\omega)}{\cos \theta}\]这样我们就可以用$\Lambda (\omega)$来表示$G_1 (\omega)$：
\[G_1 (\omega) = \frac{1}{1 + \Lambda (\omega)}\]需要说明的是仅有NDF并不能推导出唯一的$\Lambda (\omega)$，实际上对于同一个分布函数可能存在若干个$\Lambda$满足约束条件。如果我们进一步假定微表面上每一点的高度与相邻点的高度无关则可以推导出很多分布函数唯一对应的$\Lambda$。尽管这个假定在现实中可能是过于强的，但推导出的$\Lambda$和真实测量的$\Lambda$相比也是足够准确的。
假设微表面上每一点的高度分布与相邻点之间没有关联，各向同性Beckmann–Spizzichino分布函数对应的$\Lambda$为：
\[\Lambda (\omega) = \frac{1}{2} \bigg( \text{erf} (a) - 1 + \frac{e^{-a^2}}{a \sqrt{\pi}} \bigg)\]其中$a = \frac{1}{\alpha \tan \theta}$，$\text{erf} (x) = \frac{2}{\sqrt{\pi}} \int_0^\pi e^{-x’^2} d x’$。对于各向异性的材质，它的$\Lambda$函数则不能简单将$\alpha$推广成$\alpha_x$和$\alpha_y$，通常情况下需要插值出一个等效的$\alpha$来进行计算。
Trowbridge–Reitz分布对应的$\Lambda$则十分简单：
\[\Lambda (\omega) = \frac{-1 + \sqrt{1 + \alpha^2 \tan^2 \theta}}{2}\]不同$\alpha$对应的$G_1$函数的形式可见下图。显然越粗糙的表面($\alpha$越大)$G_1$函数衰减得越快。
与微表面几何性质相关的最后一个函数为$G(\omega_o, \omega_i)$，它表示在$\omega_i$和$\omega_o$两个方向都可见的微表面面积比例。如果我们假设微表面在各个方向上的可见性是相互独立的，有：
\[G(\omega_o, \omega_i) = G_1(\omega_o) G_1(\omega_i)\]在实践中发现各个方向上的可见性是不完全独立的，显然观察的角度越接近$G_1(\omega_i)$和$G_1(\omega_o)$之间的关联就越大。因此更精确的模型为：
\[G(\omega_o, \omega_i) = \frac{1}{1 + \Lambda(\omega_o) + \Lambda(\omega_i)}\]The Torrance–Sparrow Model
Torrance和Sparrow通过假设微表面满足光滑镜面来对金属材质进行建模，这样就只需要考虑法向为半程向量(half-angle vector)$\omega_h$的微表面：
\[\omega_h = \widehat{\omega_i + \omega_o}\]
接下来我们开始推导Torrance–Sparrow模型的BRDF。首先考虑法向为$\omega_h$的微表面接收的能量微分为：
\[d \Phi_h = L_i(\omega_i) \ d \omega \ d A^\perp (\omega_h) = L_i(\omega_i) \ d \omega \cos \theta_h \ d A (\omega_h)\]
根据NDF有
\[d A(\omega_h) = D(\omega_h) \ d \omega_h \ d A\]将它代回到能量微分中得到
\[d \Phi_h = L_i(\omega_i) \ d \omega \cos \theta_h \ D(\omega_h) \ d \omega_h \ d A\]根据Fresnel定律，反射出的能量通量为：
\[d \Phi_o = F_r (\omega_o) \ d \Phi_h\]再利用radiance的定义，可以得到反射出的radiance为：
\[L_o (\omega_o) = \frac{d \Phi_o}{d \omega_o \ d A \ \cos \theta_o}\]把上面的式子结合起来可以得到：
\[L_o (\omega_o) = \frac{F_r (\omega_o) L_i(\omega_i) \ d \omega_i \ \cos \theta_h \ D(\omega_h) \ d \omega_h}{\cos \theta_o \ d \omega_o}\]对于完美镜面反射，微分立体角满足：
\[d \omega_h = \frac{d \omega_o}{4 \cos \theta_h}\]这样得到化简后的反射radiance：
\[L_o (\omega_o) = \frac{F_r (\omega_o) L_i(\omega_i) D(\omega_h) \ d \omega_i}{4 \cos \theta_o}\]再加入几何可见项后就可以得到Torrance–Sparrow模型的BRDF：
\[f_r (\omega_o, \omega_i) = \frac{D(\omega_h) \ G(\omega_o, \omega_i) \ F_r(\omega_o)}{4 \cos \theta_o \ \cos \theta_i}\]Torrance–Sparrow模型的优势在于它没有假设材料的NDF或是Fresnel反射率，因此我们可以使用任意形式的NDF也可以直接使用电介质或是导体的Fresnel反射率计算公式。
我们同样可以推导出Torrance–Sparrow模型的BTDF。假设光线从折射率为$\eta_i$的材质射向折射率为$\eta_o$的材质，此时微分立体角的关系为：
\[d \omega_h = \frac{\eta_o^2 \vert \omega_o \cdot \omega_h \vert \ d \omega_o}{\big( \eta_i (\omega_i \cdot \omega_h) + \eta_o (\omega_o \cdot \omega_h) \big)^2}\]这样BTDF的表达式为：
\[f_r (\omega_o, \omega_i) = \frac{D(\omega_h) \ G(\omega_o, \omega_i) \ (1 - F_r(\omega_o))}{\big( (\omega_o \cdot \omega_h) + \eta ((\omega_i \cdot \omega_h)) \big)^2} \ \frac{\vert \omega_i \cdot \omega_h \vert \ \vert \omega_o \cdot \omega_h \vert}{\cos \theta_o \ \cos \theta_i}\]其中$\eta = \frac{\eta_i}{\eta_o}$，折射的半程向量满足$\omega_h = \omega_o + \eta \omega_i$。
Fresnel Incidence Effects
很多BRDF模型没有考虑到材料的表面性质对反射能量的影响。假设我们对木质材料进行抛光或是对墙面进行粉刷，材质本身没有发生任何变化但对表面的处理会改变光线在材料表面的行为。在这种情况下就不能按照通用的BRDF来考虑着色了。
针对这种现象Ashikhmin和Shirley开发了一种多层材质的BRDF，它假设光线在材料表面发生光泽镜面反射但在下层则发生漫反射。显然下层的漫反射材料会因为光线已经在表面发生镜面反射而接收到更少的能量。
这样整体的BRDF就是上层光泽表面与下层漫反射材料BRDF的加权平均。经过推导，光泽表面的BRDF为：
\[f_r (p, \omega_o, \omega_i) = \frac{D(\omega_h) F(\omega_o)}{4 (\omega_h \cdot \omega_i) \ \max(\mathbf{n} \cdot \omega_o, \mathbf{n} \cdot \omega_i)}\]而对于下层漫反射材料则使用了Schlick近似来逼近Fresnel反射率：
\[F_r(\cos \theta) = R + (1 - R) (1 - \cos \theta)^5\]其中$R$为光线垂直入射时材料的反射率。这样漫反射材质的BRDF为：
\[f_r (p, \omega_o, \omega_i) = \frac{28 R_d}{23 \pi} (1 - R_s) \bigg( 1 - \big(1 - \frac{(\mathbf{n} \cdot \omega_i)}{2} \big)^5 \bigg) \bigg( 1 - \big(1 - \frac{(\mathbf{n} \cdot \omega_o)}{2} \big)^5 \bigg)\]最后将两个BRDF相加即为整体材料的BRDF。
Fourier Basis BSDFs
尽管Torrance–Sparrow模型和Oren–Nayar已经能够表达很多常见的材质，但仍然由很多材质的BSDF是难以通过解析的形式来进行表示。对于这部分材质我们可以通过专业设备进行采样并建立一个3D或是4D的表格，这样在进行计算时可以通过查表的方式获取BSDF。在PBRT中使用FourierBSDF类来表示这种需要进行采样和查表的BSDF模型。
首先我们需要将BSDF展开成入射方向和出射方向的函数：
\[f(\omega_i, \omega_o) = f(\mu_i, \phi_i, \mu_o, \phi_o)\]其中$\mu_i$和$\mu_o$分别是入射光线和出射光线天顶角(zenith angles)的余弦，而$\phi_i$和$\phi_o$则是对应光线的方位角(azimuth angles)。
对于各向同性材料，BSDF仅与方位角的相对关系有关：
\[f(\omega_i, \omega_o) = f(\mu_i, \mu_o, \phi_i - \phi_o) = f(\mu_i, \mu_o, \phi)\]同时，BSDF也是相对方位角的偶函数：
\[f(\mu_i, \mu_o, \phi) = f(\mu_i, \mu_o, -\phi)\]因此我们可以用一个仅包含余弦项的Fourier级数来表示BSDF：
\[f(\mu_i, \mu_o, \phi) = \sum_{k=0}^{m-1} a_k(\mu_i, \mu_o) \cos \big(k (\phi_i - \phi_o) \big)\]对于每一项系数$a_k (\mu_i, \mu_o)$则可以用一个$n \times n$的矩阵来表示。这样整个BSDF就可以通过$m$个矩阵来表示，每个矩阵对应BSDF在对应频率上的响应。显然要存储整个BSDF就需要$m \times n \times n$个数字作为存储空间。
同时在计算余弦项的时候还可以使用递推来减少计算和存储的复杂度：
\[\cos (k \phi) = (2 \cos \phi) \cos((k-1) \phi) - \cos ((k-2) \phi)\]Spline Interpolation
在查表时需要通过插值来获得Fourier级数的每一项系数，在PBRT中使用Catmull–Rom样条来实现插值。以一维函数为例，假设已知函数$f$和它的导数$f’$在离散点$x_0, x_1, …, x_k$处的取值，则在区间[$x_i$, $x_{i+1}$]上可以使用一个三次函数来近似原函数：
\[p_i(x) = a x^3 + b x^2 + c x + d\]为简化讨论，我们假设[$x_i$, $x_{i+1}$] = [0, 1]。带入边界处的函数值和导数作为约束条件，可以解出样条函数的4个系数：
\[a = f'(x_0) + f'(x_1) + 2 f(x_0) - 2 f(x_1)\] \[b = 3 f(x_1) - 3 f(x_0) - 2 f'(x_0) - f'(x_1)\] \[c = f'(x_0)\] \[d = f(x_0)\]把插值函数写成端点函数(导数)值的组合，有
\[\begin{aligned} p(x) &= (2 x^3 - 3 x^2 + 1) f(x_0) \\ &+ (-2 x^3 + 3 x^2) f(x_1) \\ &+ (x^3 - 2 x^2 + x) f'(x_0) \\ &+ (x^3 - x^2) f'(x_1) \end{aligned}\]在大多数情况下我们很难直接得到函数的导数，此时可以使用差分来代替导数：
\[f'(x_0) = \frac{f(x_1) - f(x_{-1})}{1 - x_{-1}}\] \[f'(x_1) = \frac{f(x_2) - f(x_0)}{x_2}\]这样就可以把插值函数写成端点函数值的组合：
\[\begin{aligned} p(x) &= \frac{x^3 - 2 x^2 + x}{x_{-1} - 1} f(x_{-1}) \\ &+ \bigg( 2 x^3 - 3 x^2 + 1 - \frac{x^3 - x^2}{x_2} \bigg) f(x_0) \\ &+ \bigg(-2 x^3 + 3 x^2 + \frac{x^3 - 2 x^2 + x}{1 - x_{-1}} \bigg) f(x_1) \\ &+ \frac{x^3 - x^2}{x_2} f(x_2) \end{aligned}\]注意到4个系数之间存在一定的关系，因此可以将上式改写为：
\[p(x) = w_0 f(x_{-1}) + w_1 f(x_0) + w_2 f(x_1) + w_3 f(x_2)\]其中
\[w_0 = \frac{x^3 - 2 x^2 + x}{x_{-1} - 1}\] \[w_1 = 2 x^3 - 3 x^2 + 1 - \frac{x^3 - x^2}{x_2} = (2 x^3 - 3 x^2 + 1) - w_3\] \[w_2 = -2 x^3 + 3 x^2 + \frac{x^3 - 2 x^2 + x}{1 - x_{-1}} = (-2 x^3 + 3 x^2) - w_0\] \[w_3 = \frac{x^3 - x^2}{x_2}\]对于更一般的情况则需要将区间[$x_i$, $x_{i+1}$]进行规范化，我们引入$t$作为缩放后的变量来代替$x$：
\[t = \frac{x - x_0}{x_1 - x_0}\] \[w_0 = (t^3 - 2 t^2 + t) \frac{x_1 - x_0}{x_{-1} - x_1}\] \[w_1 = 2 t^3 - 3 t^2 + 1 - w_3\] \[w_2 =-2 t^3 + 3 t^2 - w_0\] \[w_3 = (t^3 - t^2) \frac{x_1 - x_0}{x_2 - x_0}\]