轨迹相似度度量方法总结

发表于 2019-06-10 | 更新于 2019-06-11

轨迹相似度度量有广泛的用途，如语音识别分类、模板匹配、信息检索等。下面介绍各种轨迹相似度度量方法。

轨迹定义

轨迹可由时间域到空间域的映射函数表示，如：

$$t\stackrel{F}{\longrightarrow}R^d, d>1$$

度量方法

轨迹相似度度量方法主要有：

基于点方法： EDR，LCSS，DTW等
基于形状的方法： Frechet， Hausdorff
基于分段的方法：One Way Distance, LIP distance
基于特定任务的方法：TRACLUS， Road Network，grid等

基于点的方法

DTW

DTW（Dynamic Time Warping, 动态时间规整）可以计算两个时间序列的相似度，尤其适用于不同长度、不同节奏的时间序列（比如不同的人读同一个词的音频序列）。DTW将自动warping扭曲时间序列（即在时间轴上进行局部的缩放），使得两个序列的形态尽可能的一致，得到最大可能的相似度。

Dynamic Time Warping（DTW）诞生有一定的历史了（日本学者Itakura提出），它出现的目的也比较单纯，是一种衡量两个长度不同的时间序列的相似度的方法。应用也比较广，主要是在模板匹配中，比如说用在孤立词语音识别（识别两段语音是否表示同一个单词），手势识别，数据挖掘和信息检索等中。

设 $P=<p_1,p_2,…,p_m>$ 和 $Q=<q_1,q_2,…,q_n>$ 是两个时间序列，则 $P$ 和 $Q$的距离 $DTW(P,Q)$ 定义如下：

$$DTW(P,Q)= \left{
\begin{array}{lcl}
0 & &if\ m=n=0\
\infty & &if\ m=0\ or\ n=0\
dist(p_1,q_1)+min
\left{
\begin{aligned}
DTW(Rest(P),Rest(Q))\
DTW(Rest(P),Q)\
DTW(P,Rest(Q))\
\end{aligned}
\right} & & otherwise
\end{array}
\right.
$$

参考链接

经典控制与现代控制理论的区别与联系

发表于 2019-06-09 | 更新于 2021-01-22

控制理论是工程学与数学的跨领域分支，主要处理在有输入信号的动力系统的行为。系统的外部输入称为“参考值”，系统中的一个或多个变量需随着参考值变化，控制器处理系统的输入，使系统输出得到预期的效果。

控制理论一般的目的是借由控制器的动作让系统稳定，也就是系统维持在设定值，而且不会在设定值附近晃动。

控制论简介

控制理论是

一个研究如何调整动态系统特性的理论。
科学中跨学科的领域，起源于工程及数学，逐渐的应用在许多社会科学中，例如心理学、社会学（社会学中的控制理论）、犯罪学及金融系统。

控制系统可以视为具有四种机能的系统：量测、比较、计算及修正。这四个机能可以用五种元素来实现：感测器、换能器、发送器、控制器及最终控制元件。量测机能是由感测器、换能器及发送器执行，在实务应用上，这三个元素会整合在一个单体内，像是电阻温度计。比较和计算的机能是由控制器执行，可能是电子式的比例控制（P控制）、PI控制、PID控制、双稳态的迟滞控制，也可能是可编程逻辑控制器（PLC）。早期的控制器也可能是机械式的，像是离心式调速器或是化油器。修正机能是由最终控制元件执行，最终控制元件改变系统的输出，因此影响操纵或控制的变量。

范例

车辆的巡航定速系统是让车辆维持在由驾驶者设定的固定参考速度。此时控制器为巡航定速系统，车辆为受控体（plant），而系统是由控制器和车辆所组成，而控制变量是引擎节流阀的位置．会决定引擎可以产生的功率。

一种最单纯的作法是当驾驶者启动巡航定速系统时，固定引擎节流阀的位置。但是若驾驶者在平坦的路面启动巡航定速系统，车辆在上坡时速度会较慢，车辆在下坡时速度又会较快。这种的控制器称为开环控制器，因为没有去量测系统输出（车辆速度）并且影响控制变量（节流阀位置），因此此系统无法去针对车辆遇到的变化（像路面坡度的变动）去进行调整。

在闭环控制系统中，利用感测器量测系统输出（车辆速度），并将资料送入控制器中，控制器依资料调整控制变量（节流阀位置），来达到维持理想系统输出（使车辆速度和驾驶者设定的参考速度一致）。此时若车辆在上坡时，感测器会量到车辆的速度变慢，因此会调整节流阀位置，加大引擎输出功率，使马达加速。因为有量测车辆速度的回授，因此控制器可以配合车辆速度的变化进行动态调整。因此产生了控制系统中的“环”范式：控制变量影响系统输出，而再根据量测到的系统输出去调整控制变量。

经典控制理论

为了克服开环控制器的限制，在控制理论中导入了反馈。闭环的控制器利用回授来控制动态系统的状态或输出。其名称来自系统中的讯息路径：程序输入（例如马达的电压）影响程序输出（例如马达的电流或转矩），利用感测器量测输出，再将量测资料送到控制器中处理，结果送回控制器作为输入信号之一，因此成为一闭环。

相对于开环控制器，闭环控制器有以下的优点：

噪声抑制能力（像巡航定速中的路面坡度）。
即使在数学模型有一些不确定性的情形下（如模型结构和实际系统不是完全符合，或是模型参数和实际数值不是完全一致），仍有一定程度的性能。
可以稳定不稳定的系统
减少对于参数变动的灵敏度
提升命令追随（命令变化时，系统配合命令变化）的性能
有些系统中，同时出现开环及闭环的控制，此时的开环会称为前馈，目的是为了提升命令追随的性能。

PID控制器是常见的闭回路控制器架构。

闭环传递函数

系统的输出y(t)借由感测器F量测后，和参考值r(t)相减，控制器C根据参考值和输出值的误差e调整受控体P的输入u，如图1所示，这类的控制器称为闭环控制器。

由于只有一个输入和输出，此系统会称为SISO（单一输入单一输出）控制系统。MIMO（多重输入多重输出）控制系统是指输入或输出不只一个，在实际应用上也很常见，其输入变量和输出变量会用向量表示，而不是单一数值的标量。在分布参量系统中，向量可能是无限维的，即一般的函数。

闭环传递函数示意图

图1 闭环传递函数示意图

若假设控制器C、受控体P及感测器F都是线性及非时变的（各模组输入和输出的关系不随时间改变），可以将上述系统用拉普拉斯转换来分析，因此可以得到以下的关系：

$$Y(s)=P(s)U(s),!$$
$$U(s)=C(s)E(s),!$$
$$E(s)=R(s)-F(s)Y(s),!$$

其中 s为拉普拉斯转换中的复变量，若要求解Y(s)用R(s)表示，可得：

$$Y(s)=\left({\frac {P(s)C(s)}{1+F(s)P(s)C(s)}}\right)R(s)=H(s)R(s),!$$

表示式 $H(s)={\frac {P(s)C(s)}{1+F(s)P(s)C(s)}},!$即为系统的闭环传递函数，分子是从r到y的前馈（开环）增益，分母是1加上经过反馈环的增益．即闭环增益，若$|P(s)C(s)|\gg 1,!$，，也就是说在各s下，其范数都很大，且 $|F(s)|\approx 1,!$，则Y(s)近似于R(s)，此时输出会紧密的追随参考输入。

PID

请参考链接[4]。

现代控制理论

经典控制理论以频域分析为主，而现代控制理论利用时域的状态空间表示法，将系统中的输入、输出及状态变量之间的关系用一阶的微分方程表示。为了抽象化输入、输出及状态变量的数量，这些变量一般会用向量来表示，而微分方程或代数方程（当系统是线性时）则会以矩阵形式表示。状态空间表示法也称为时域分析，提供一个方便且简洁的方式针对多重输入及输出的系统建模及分析，在有输入和输出时，也可以利用拉氏转换，将系统所有的资料包括在其中。现代控制理论不同于频域分析，可以分析非线性或不是零初始条件的系统。状态空间就是指坐标轴为状态变量的空间，系统的状态可以表示为状态空间中的一个向量。

控制理论主题

稳定性

在控制理论中的稳定性是指控制系统的状态在特定条件下，可以维持在一定的范围内，不会发散，而在什么范围内才算是稳定则依系统种类而不同。

没有输入信号的动力系统，其稳定性是用李雅普诺夫稳定性来描述，也就是任何初始条件在 $x_{0}$ 附近的轨迹均能维持在 $x_{0}$ 附近
有输入信号的线性系统，其稳定性是用有界输入有界输出稳定性（BIBO 稳定性）来描述，针对任何有界的输入信号，其输出也是有界。
有输入信号的非线性系统，其稳定性是用输入-状态稳定性（input-to-state stability），结合了李雅普诺夫稳定性及类似有界输入有界输出稳定性的表示方式。

可控制性及可观测性

可控制性和可观测性分别是输入和状态，输出和状态之间的性质。是在分析控制系统，决定控制策略或判断是否可以使系统稳定时所需要的重要性质。

可控制性是指是可以用适当的控制信号作为输入，使特定状态变量的数值变成０，和利用输入调整状态变量的能力有关，若一个状态变量是不可控制的，表示没有输入可以调整这一个状态，若一系统中所有不可控制的状态变量，其动态特性都是稳定的，则此系统称为可稳定的（stabilizable）。

可观测性是指可以用输出的量测及计算得到状态变量的值，若一个状态变量是不可观测的，表示无法确认此一状态是否稳定，也就无法用此状态来稳定整个系统。若一系统中所有不可控制的观测变量都是稳定的，则此系统称为可检测的（detectable）。

控制规格

在控制原理的基础下，已发展出许多不同的控制策略，从非常通用的（PID控制器），到针对特殊系统的控制，尤其是机器人或是航空器的巡航定速控制。

一个控制问题会有许多的规格，其中稳定性是必要条件的，不论系统开环稳定性如何，控制器需确保在闭环下是稳定的。性能不佳或是调整不当的控制器可能使系统变的不稳定，甚至可能比开环还要不稳定，这是应尽量要避免的。

模型识别及鲁棒性

控制系统一定会有一定程度的鲁棒性。控制器一般是依照一个假设的受控系统模组再进行设计，鲁棒性是指一控制器配合的受控系统和原来假设的系统有一点不同，控制器的特性不会有太大的变化。这个规格在实际的控制器中相当重要，因为很少实际系统会完全符合描述它的微分方程，在选择系统数学模型时，一般会进行简化，否则数学模型会非常复杂，甚至无法求得一个完整的模型。

系统分类

线性系统控制

针对MIMO的系统，极点的指定可以用开环系统的状态空间，再将极点放在指定位置，计算对应的回授矩阵。若在复杂的系统中，上述的程序需要用电脑辅助计算才能达到，而且不保证其鲁棒性。而且一般而言无法量到所有的系统状态，在极点指定的设计时需加入观测器（observer）的设计。

非线性系统控制

像机器人学及航天产业中的程序一般都有高度非线性的动态，在控制理论中有时可以用线性化的方式转换为线性系统，再依线性系统的方式控制。但有时需要用一些可以配合非线性系统使用的非线性控制理论，例如回授线性化、反推控制、滑动模式控制等。轨迹线性化控制一般利用李亚普诺夫稳定性的基础。微分几何用做为一数学工具，将许多广为人知的线性控制概念扩展到非线性控制中，但其中又有其微妙之处，因此变成一个更有挑战性的问题。

分散式系统

分散控制系统是指一个系统由多个控制器来控制。分散控制有几个好处，例如可以控制一个位在广大地理区域的系统，各控制器之间可以用通讯网络彼此交换资料，并协调彼此的行动。

参考链接

ROS中Package安装方法

发表于 2019-06-08

ROS（机器人操作系统，Robot Operating System），是专为机器人软件开发所设计出来的一套电脑操作系统架构。它是一个开源的元级操作系统（后操作系统），提供类似于操作系统的服务，包括硬件抽象描述、底层驱动程序管理、共用功能的执行、程序间消息传递、程序发行包管理，它也提供一些工具和库用于获取、建立、编写和执行多机融合的程序。

ROS的Package资源非常丰富，官方库中就有两千多Package，这对扩充ROS的功能十分重要。下面即介绍ROS中Package的安装方法，主要分成两种方法：

Deb安装方式

deb方式安装方法十分简单，根据ROS版本，直接运行apt-get命令，例如：

1	$ sudo apt-get install ros-kinetic-camera-calibration

源码安装方式

源码安装方式稍微复杂，安装方法如下：

创建catkin工作空间
在catkin工作空间的src文件夹下，下载ROS的Package源代码
使用catkin build命令编译安装

参考链接

安装ROS软件包,by ferstar.

ROS编译命令catkin简析

发表于 2019-06-06

目前编译ROS的Package有两种方法：

catkin_make
catkin build

catkin_make

catkin_make 是一个命令行工具，它简化了catkin的标准工作流程。你可以认为catkin_make是在CMake标准工作流程中依次调用了cmake 和 make。

使用方法如下：

1 2	# 在catkin工作空间下 $ catkin_make [make_targets] [-DCMAKE_VARIABLES=...]

catkin

catkin是一个用于处理catkin元构建系统和catkin工作区的命令行工具。其用法如下：

`catkin VERB -h` for help on each verb listed below:
  
    build	Builds a catkin workspace.
    clean	Deletes various products of the build verb.
    config	Configures a catkin workspace's context.
    create	Creates catkin workspace resources like packages.
    env	        Run an arbitrary command in a modified environment. 
    init	Initializes a given folder as a catkin workspace.
    list	Lists catkin packages in the workspace or other arbitray folders.
    locate	Get the paths to various locations in a workspace.
    profile	Manage config profiles for a catkin workspace.

同样可使用catkin build命令编译ROS的package。

catkin_make与catkin build的区别

与catkin_make不同，catkin命令行工具不仅仅是围绕cmake和make命令的瘦包装器。 catkin build命令隔离地在工作空间的源空间中构建每个包，以防止构建时串扰。因此，在其最简单的用法中，catkin构建的行为类似于catkin_make_isolated的并行化版本。

参考链接

Migrating from catkin_make,by catkin_tools homepage.
编译ROS程序包,by ros wiki.

卡尔曼滤波入门

发表于 2019-06-04 | 更新于 2023-01-31

卡尔曼滤波（Kalman filter）是一种高效率的递归滤波器（自回归滤波器），它能够从一系列的不完全及包含噪声的测量中，估计动态系统的状态。卡尔曼滤波会根据各测量量在不同时间下的值，考虑各时间下的联合分布，再产生对未知变数的估计，因此会比只以单一测量量为基础的估计方式要准。卡尔曼滤波得名自主要贡献者之一的鲁道夫·卡尔曼。

卡尔曼滤波在技术领域有许多的应用。常见的有飞机及太空船的导引、导航及控制。卡尔曼滤波也广为使用在时间序列的分析中，例如信号处理及计量经济学中。卡尔曼滤波也是机器人运动规划及控制的重要主题之一，有时也包括在轨迹最佳化。卡尔曼滤波也用在中轴神经系统运动控制的建模中。因为从给与运动命令到收到感觉神经的回授之间有时间差，使用卡尔曼滤波有助于建立符合实际的系统，估计运动系统的目前状态，并且更新命令。

卡尔曼滤波原理解析

卡尔曼滤波器的状态矩阵方程如图1所示。

卡尔曼滤波器的状态矩阵方程

图1 卡尔曼滤波器的状态矩阵方程

其中，下角标上的k是状态。此处我们将其视为离散时间间隔，比如说k=1 代表 1ms, k=2 代表2ms。

我们的目的是找到信号 $x$ 的估值 $\hat{x} $，并且希望能对所有的k值都能找到对应的估值。

另外此处的 $Z_k$ 是实际测量值，记住我们对该值并不完全信任，否则我们也不用费这么多事了。 $K_k$ 称为卡尔曼增益（也是最重要的量）， $\hat{x}_k$ 是前一状态下的信号估值。

现在我们有了测量值，前一状态的信号估值。该方程中唯一未知的量就是卡尔曼增益 $K_k$ 了。对于每个状态，我们都需要计算对应的。这事不简单，但好在我们有所需的计算工具。

另一方面，假设 $K_k$ 等于0.5，我们会发现该式变成了一个简单的求平均值公式。换句话说，随着状态的变化，我们的 $K_k$ 值将越来越“聪明”。

卡尔曼滤波器的构造

建立模型

此步最为关键，你必须确保卡尔曼滤波器适用于你要解决的问题。

卡尔曼滤波器的两个方程如下:

$$x_k=Ax_{k-1}+Bu_k+w_{k-1} \tag{1}$$
$$z_k=Hx_k+v_k \tag{2}$$

式(1)表达的是每个 $x_k$ 都可以通过一个线性随机方程估计出来。任意 $x_k$ 都是其前一时刻的值与过程噪音的线性组合（这个很难概念化）。请记住，大部分情况下该式没有控制信号 $u_k$ 项。

式(2)告诉我们任何测量值 $z_k$ （无法确定精确与否的测量值）都是信号值与测量噪声的线性组合。这两个分量符合高斯分布。

过程噪声与测量噪声互相统计独立。

$A, B, H$ 是一般形式的矩阵。但在大多数信号处理问题中，这些量仅为数值。而且虽然这些值在状态变换时会改变，大多数情况下我们都可以假设他们为定值。

如果我们十分确定我们的系统符合此模型，那么唯一剩下要做的事就是估计噪音函数 $w_{k-1}$ 和 $v_k$ 的平均值以及标准差。我们知道，在实际生活中没有信号满足高斯分布，但我们可以近似其为高斯分布。

该近似问题不大，因为我们将看到卡尔曼滤波器算法会逐渐向正确的（噪音函数的）估计值收敛，即使高斯噪声参数估计不佳。

唯一需要记住的是：你估计出来的噪音参数越好（越接近实际），你估计的（输出真实值）就越好。

开始卡尔曼滤波

如果你的模型适用于卡尔曼滤波器，那么接下来的步骤就是决定一些必要的参数以及初始值。

卡尔曼滤波器包含的方程可分为两个方程集：时间更新方程组（用于预测）以及测量更新方程组（用于修正）。这两个方程组在滤波器运行的每一步（每个状态）下都会执行，如图2所示。

卡尔曼滤波的两个步骤

图2 卡尔曼滤波的两个步骤

建模部分已经在步骤一完成了，所以矩阵A，B和H已知。这些矩阵很可能是一个常数，而且大部分情况下会等于1。

剩下的最让人难受的部分就是决定R和Q的值了。R的值还是很容易找的，因为一般情况下我们对环境中的噪音还是能够确认的。（起码能用仪器测一下）。但是找Q的值就没那么直观了。

为了使滤波器能够运行，我们需要知道 $x_0$ 和 $P_0$ 的估计值。

迭代

在获得了滤波器运行所需的所有信息后，我们就可以估值迭代了。记住：前一状态的估值将成为当前状态的输入。

卡尔曼滤波的迭代运行

图3 卡尔曼滤波的迭代运行

此处 $\hat{x}_k^-$ 是预估值，从某种角度来说是第二部分运行前对 x 的一个粗略估计值。

同时 $P_k^-$ 叫做预估误差协方差。在第二步“测量更新”中我们将会用到这两个预估值。

$\hat{x}_k$为在时间 k 时的 x 的估计值。（也是我们最想获得的值）。同时，我们得到了用于k+1时刻计算的 $P_k$ 值。

下一次迭代不会用到我们求得的卡尔曼增益 $K_k$ 的值，该值隐藏而神秘，并且是这些方程集的最重要的部分。

我们在第二步“测量更新”中求得的值也叫做后部值(posterior values)。这个名称也很说得通。

卡尔曼滤波器的应用示例

参考链接

求解射线与三角形交点的算法

发表于 2019-06-04

求解射线与三角形的交点，在光线追踪、碰撞检测、目标拾取等场景中经常使用，是计算机图形学中最基本的操作。下面介绍常用的求解射线与三角形交点的算法。

问题定义

求解射线与三角形交点示意图如图1所示。

求解射线与三角形交点示意图

图1 求解射线与三角形交点示意图

射线的参数方程如下，其中O是射线的起点，D是射线的方向，t是常数。

$$O+Dt$$

该方程的含义是一个点从起点O开始，沿着方向D移动任意长度，得到终点R，根据t值的不同，得到的R值也不同，所有这些不同的R值便构成了整条射线，比如下面的射线，起点是P0，方向是u，p0 + tu也就构成了整条射线。

射线方程示意图

图2 射线方程示意图

三角形的参数方程如下，其中$V_0$，$V_1$和$V_2$是三角形的三个点，$u, v$是$V_1$和$V_2$的权重，$1-u-v$是$V_0$的权重，并且满足$u>=0, v >= 0,u+v<=1$。

$$(1-u-v)V_0+uV_1+vV_2$$

三角形方程示意图

图3 三角形方程示意图

直观方法

求解射线与三角形的交点最直观的方法如下：

判断射线是否与平面相交
判断点是否在三角形内

但该方法需要额外计算三角形所在平面，效率不高。

Moller-Trumbore方法(Journal of Graphic Tools, 1997)

Moller-Trumbore方法中，求射线与三角形的交点即求解如下方程:

$$O+Dt=(1-u-v)V_0+uV_1+vV_2$$

其中t,u,v是未知数，其他都是已知的。

移项并整理，将t,u,v提取出来作为未知数，得到下面的线性方程组：

$$\begin{bmatrix}
-D& V_1-V_0 &V_2-V_0
\end{bmatrix}\begin{bmatrix}
t\u\v
\end{bmatrix}=O-V_0$$

现在开始解这个方程组，这里要用到两个知识点，一是克莱姆法则，二是向量的混合积。

令$E_1 = V_1 - V_0，E_2 = V_2 - V_0，T = O - V_0$上式可以改写成:

$$\begin{bmatrix}
-D& E_1 & E_2
\end{bmatrix}\begin{bmatrix}
t\u\v
\end{bmatrix}=T$$

根据克莱姆法则，可得到t,u,v的解为：

$$
\begin{bmatrix}
t\u\v
\end{bmatrix}
=\frac{1}{\begin{vmatrix}
-D & E_1 & E_2
\end{vmatrix} }
\begin{vmatrix}
T&E_1&E_2\
-D& T& E_2\
-D & E_1& T
\end{vmatrix}
$$

根据混合积公式:

$$\begin{vmatrix}
a&b&c
\end{vmatrix}
=a\times{b}\cdot{c}$$

上式改写为：

$$
\begin{bmatrix}
t\u\v
\end{bmatrix}
=\frac{1}{\begin{bmatrix}
-D \times E_2 \cdot E_1
\end{bmatrix} }
\begin{vmatrix}
T \times E_1 \cdot E_2\
D \times E_2 \cdot T\
T \times E_1 \cdot D
\end{vmatrix}
$$

令$P=D \times E_2$，$Q=T \times E_1$，得到最终的公式:

$$\begin{bmatrix}
t\u\v
\end{bmatrix}
=\frac{1}{\begin{bmatrix}
P \cdot E_1
\end{bmatrix} }
\begin{vmatrix}
Q \cdot E_2\
P \cdot T\
Q \cdot D
\end{vmatrix}
$$

之所以提炼出P和Q是为了避免重复计算。

参考链接

射线和三角形的相交检测（ray triangle intersection test）,by zdd.
光线-三角形求交测试算法[译], by PKUWWT.
克莱姆法则,by wikipedia.
混合积,by wikipedia.
重心坐标,by wikipedia.
重心坐标（Barycentric coordinates）,by 杨超.

PLY文件格式分析

发表于 2019-06-04

PLY文件是一种存储3D模型的文件格式，全名为多边形档案（Polygon File Format）或史丹佛三角形档案（Stanford Triangle Format）。

该格式主要用以储存立体扫描结果的三维数值，透过多边形片面的集合描述三维物体，与其他格式相较之下这是较为简单的方法。它可以储存的资讯包含颜色、透明度、表面法向量、材质座标与资料可信度，并能对多边形的正反两面设定不同的属性。

在档案内容的储存上PLY有两种版本，分别是纯文字（ASCII）版本与二元码（binary）版本，其差异在储存时是否以ASCII编码表示元素资讯。

通过分析PLY文件格式，我们可以进而初窥存储3D模型的文件奥秘。

文件格式

各种文件格式通常分成文件头和文件内容，PLY格式也不例外。

每个PLY档都包含档头（header），用以设定网格模型的“元素”与“属性”，以及在档头下方接着一连串的元素“数值资料”。一般而言，网格模型的“元素”就是顶点（vertices）、面（faces），另外还可能包含有边（edges）、深度图样本（samples of range maps）与三角带（triangle strips）等元素。无论是纯文字与二元码的PLY档，档头资讯都是以ASCII编码编写，接续其后的数值资料才有编码之分。

PLY档案以此行作为开头，以识别PLY格式：

ply

接着第二行是版本资讯，目前有三种写法：

1
2
3

format ascii 1.0
format binary_little_endian 1.0
format binary_big_endian 1.0

其中ascii, binary_little_endian, binary_big_endian是档案储存的编码方式，而1.0是遵循的标准版本（现阶段仅有PLY 1.0版）。在档头中可使用’comment’作为一行的开头以编写注解，例如：

1	comment This is a comment!

描述元素及属性，必须使用’element’及’property’的关键字，一般的格式为element下方接着属性列表，例如:

element <element name> <number in file>
property <data_type> <property name 1>
property <data_type> <property name 2>
property <data_type> <property name 3>

‘property’不仅定义了资料的型态，其出现顺序亦定义了资料的顺序。内定的资料形态有两种写法：一种是char uchar short ushort int uint float double,另外一种是具有位元长度的int8 uint8 int16 uint16 int32 uint32 float32 float64。例如，描述一个包含12个顶点的物体，每个顶点使用3个单精度浮点数 (x,y,z）代表点的座标，使用3个unsigned char代表顶点颜色，颜色顺序为 (B, G, R),则档头的写法为：

element vertex 12
property float x
property float y
property float z
property uchar blue
property uchar green
property uchar red

其中vertex是内定的元素类型，接续的6行property描述构成vertex元素的数值字段顺序代表的意义，及其资料形态。

另一个常使用的元素是面。由于一个面是由3个以上的顶点所组成，因此使用一个“顶点列表”即可描述一个面, PLY格式使用一个特殊关键字’property list’定义之。例如，一个具有10个面的物体，其PLY档头可能包含：

1 2	element face 10 property list uchar int vertex_indices

‘property list’表示该元素face的特性是由一行的顶点列表来描述。列表开头以uchar型态的数值表示列表的项目数，后面接着资料型态为int的顶点索引值（vertex_indices），顶点索引值从0开始。

最后，标头必须以此行结尾：

1	end_header

档头后接着的是元素资料（端点座标、拓朴连结等）。在ASCII格式中各个端点与面的资讯都是以独立的一行描述，而二元编码格式则连续储存这些资料，载入时须以’element’定义的元素数目以及’property’中设定的资料形态计算各笔字段的长度。

示例

一个典型的PLY文件结构分成三部分：

1
2
3

文件头 (从ply开始到end_header）
顶点元素列表
面元素列表

其中的顶点元素列表一般以x y z方式排列，形态如文件头所定义；而面元素列表是以下列格式表示。

1	<组成面的端点數N> <端点#1的索引> <端点#2的索引> ... <端点#N的索引>

以存储一个立方体模型的PLY文件为例，其内容为：

ply
format ascii 1.0
comment made by anonymous
comment this file is a cube
element vertex 8
property float32 x
property float32 y
property float32 z
element face 12
property list uint8 int32 vertex_index
end_header
0 0 0
0 25.8 0
18.9 0 0
18.9 25.8 0
0 0 7.5
0 25.8 7.5
18.9 0 7.5
18.9 25.8 7.5
3 5 1 0
3 5 4 0
3 4 0 2
3 4 6 2
3 7 5 4
3 7 6 4
3 3 2 1
3 1 2 0
3 5 7 1
3 7 1 3
3 7 6 3
3 6 3 2

PLY文件可用blender软件打开。

参考链接

PLY,by wikipedia.
Opengl学习笔记：（一）.Ply文件文件格式和文件读取,by 为何走到这里.

图像特征检测与描述

发表于 2019-06-02

特征检测（英语：Feature detection）是计算机视觉和图像处理中的一个概念。它指的是使用计算机提取图像信息，决定每个图像的点是否属于一个图像特征。特征检测的结果是把图像上的点分为不同的子集，这些子集往往属于孤立的点、连续的曲线或者连续的区域。

图像特征分类

特征的精确定义往往由问题或者应用类型决定。特征是一个数字图像中“有趣”的部分，它是许多计算机图像分析算法的起点。因此一个算法是否成功往往由它使用和定义的特征决定。因此特征检测最重要的一个特性是“可重复性”：同一场景的不同图像所提取的特征应该是相同的。

特征检测是图象处理中的一个初级运算，也就是说它是对一个图像进行的第一个运算处理。它检查每个像素来确定该像素是否代表一个特征。假如它是一个更大的算法的一部分，那么这个算法一般只检查图像的特征区域。作为特征检测的一个前提运算，输入图像一般通过高斯模糊核在尺度空间中被平滑。此后通过局部导数运算来计算图像的一个或多个特征。

常用的图像特征分成以下四类：

边缘

边缘指组成两个图像区域之间边界（或边缘）的像素。一般一个边缘的形状可以是任意的，还可能包括交叉点。在实践中边缘一般被定义为图像中拥有大的梯度的点组成的子集。一些常用的算法还会把梯度高的点联系起来来构成一个更完善的边缘的描写。

角指图像中点似的特征，在局部它有两维结构。早期的算法首先进行边缘检测，然后分析边缘的走向来寻找边缘突然转向（角）。后来发展的算法不再需要边缘检测这个步骤，而是可以直接在图像梯度中寻找高度曲率。

区域

与角不同的是区域描写一个图像中的一个区域性的结构，但是区域也可能仅由一个像素组成，因此许多区域检测也可以用来监测角。一个区域监测器检测图像中一个对于角监测器来说太平滑的区域。

区域检测可以被想象为把一张图像缩小，然后在缩小的图像上进行角检测。

长条形的物体被称为脊。在实践中脊可以被看作是代表对称轴的一维曲线，此外局部针对于每个脊像素有一个脊宽度。从灰梯度图像中提取脊要比提取边缘、角和区域困难。在空中摄影中往往使用脊检测来分辨道路，在医学图像中它被用来分辨血管。

边缘检测

边缘检测（英语：Edge detection）是图像处理和计算机视觉中的基本问题，边缘检测的目的是标识数字图像中亮度变化明显的点。图像属性中的显著变化通常反映了属性的重要事件和变化。这些包括（i）深度上的不连续、（ii）表面方向不连续、（iii）物质属性变化和（iv）场景照明变化。边缘检测是图像处理和计算机视觉中，尤其是特征检测中的一个研究领域。

图像边缘检测大幅度地减少了数据量，并且剔除了可以认为不相关的信息，保留了图像重要的结构属性。有许多方法用于边缘检测，它们的绝大部分可以划分为两类：基于查找一类和基于零穿越的一类。基于查找的方法通过寻找图像一阶导数中的最大和最小值来检测边界，通常是将边界定位在梯度最大的方向。基于零穿越的方法通过寻找图像二阶导数零穿越来寻找边界，通常是Laplacian过零点或者非线性差分表示的过零点。

边缘检测的方法

有许多用于边缘检测的方法，他们大致可分为两类：

基于搜索的边缘检测

基于搜索的边缘检测方法首先计算边缘强度，通常用一阶导数表示，例如梯度模;然后，用计算估计边缘的局部方向，通常采用梯度的方向，并利用此方向找到局部梯度模的最大值.

基于零交叉的边缘检测

基于零交叉的方法找到由图像得到的二阶导数的零交叉点来定位边缘.通常用拉普拉斯算子或非线性微分方程的零交叉点，我们将在后面的小节中描述.

滤波做为边缘检测的预处理通常是必要的，通常采用高斯滤波.

已发表的边缘检测方法应用计算边界强度的度量，这与平滑滤波有本质的不同.正如许多边缘检测方法依赖于图像梯度的计算，他们用不同种类的滤波器来估计x-方向和y-方向的梯度.

角检测

角检测（英语：Corner detection）或兴趣点检测（interest point detection），是计算机视觉系统中用来提取特征以及推测图像内容的一种方法.角检测的应用很广，经常用在运动检测，跟踪，图像镶嵌（image mosaicing），全景图缝合（panorama stiching），三维建模以及物体识别中.

问题定义

两条边的交点形成一个角（点）。而图像的要点（也称为受关注点）是指图像中具有代表性以及稳健性（即指该点能够在有噪声干扰的情况下也能稳定的被定位，在大陆亦被称为：鲁棒性）的点。也就是说，要点可以是角（点），也可以不是，例如局部亮点或暗点，线段终点，或者曲线上的曲率最大值点。在实际应用中，很多所谓的（角）点检测算法其实是检测要点，而不仅仅是角（点）。所以，如果我们只想检测角的话，还需要对检测出的要点进一步分析。

角检测的方法

在现实世界中，角点对应于物体的拐角，道路的十字路口、丁字路口等。从图像分析的角度来定义角点可以有以下两种定义：

角点可以是两个边缘的角点；
角点是邻域内具有两个主方向的特征点；
前者往往需要对图像边缘进行编码，这在很大程度上依赖于图像的分割与边缘提取，具有相当大的难度和计算量，且一旦待检测目标局部发生变化，很可能导致操作的失败。早期主要有Rosenfeld和Freeman等人的方法，后期有CSS等方法。

基于图像灰度的方法通过计算点的曲率及梯度来检测角点，避免了第一类方法存在的缺陷，此类方法主要有Moravec算子、Forstner算子、Harris算子、SUSAN算子等。

角检测方法发展历程

图1 角检测方法发展历程

Harris算法

Harris角点检测基本原理如图2所示。人眼对角点的识别通常是在一个局部的小区域或小窗口完成的。如果在各个方向上移动这个特征的小窗口，窗口内区域的灰度发生了较大的变化，那么就认为在窗口内遇到了角点。如果这个特定的窗口在图像各个方向上移动时，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点；如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生变化，那么，窗口内的图像可能就是一条直线的线段。

Harris角点检测基本原理

图2 Harris角点检测基本原理

SIFT算法

SIFT(Scale-invariant features transform, 尺度不变特征变换)是一种检测局部特征的算法，该算法通过求一幅图中的特征点（interest points,or corner points）及其有关scale 和 orientation 的描述子得到特征并进行图像特征点匹配，获得了良好效果。SIFT特征不只具有尺度不变性，即使改变旋转角度，图像亮度或拍摄视角，仍然能够得到好的检测效果。该算法由 David Lowe在1999年所发表，2004年完善总结。

SIFT算子是把图像中检测到的特征点用一个128维的特征向量进行描述，因此一幅图像经过SIFT算法后表示为一个128维的特征向量集，该特征向量集具有对图像缩放，平移，旋转不变的特征，对于光照、仿射和投影变换也有一定的不变性，是一种非常优秀的局部特征描述算法。

SIFT算法的流程分别为：

尺度空间极点检测
关键点精确定位
关键点的方向确定
特征向量的生成

ORB算法

ORB的全称是ORiented Brief，是文章ORB: an efficient alternative to SIFT or SURF中提出的一种新的角点检测与特征描述算法。实际上，ORB算法是将FAST角点检测与BRIEF特征描述结合并进行了改进。

ORB特征由关键点和描述子两部分组成。它的关键点称为“Oriented Fast”，是一种改进的FAST角点。描述子称为“BRIEF(Binary Robust Independent Elementary Feature)”。因此提取ORB特征分为如下两个步骤：

提取FAST角点(相较于原版FAST角点，ORB中计算了特征点的主方向，为后续的BRIEF描述子增加了旋转不变性)
计算BRIEF描述子

参考链接

特征检测,by wikipedia.
边缘检测,by wikipedia.
角检测,by wikipedia.
Harris角点,by Ronny.
特征点匹配——SIFT算法详解,by lhanchao.
ORB特征提取、匹配及实现,by zhaoxuhui.
传统计算机视觉中图像特征匹配方法的原理介绍（SIFT 和 ORB）,by Zhang Bin.

相机位姿估计入门

发表于 2019-05-29 | 更新于 2019-06-01

相机位姿估计是指给定若干图像，估计其中相机运动的问题。求解方法通常分特征点法和直接法两种。下面主要介绍特征点法。

特征点法的思路是先从图像当中提取许多特征，然后在图像间进行特征匹配，这样就得到许多匹配好的点，再根据这些点进行相机位姿的求解。根据传感器形式的不同，可以分成三种情况：

2D-2D，单目相机获取的影像，只能获得像素坐标
3D-3D配对点，RGBD或双目相机，可以获取深度信息
3D-2D，已知一张图中的3D信息，另一张图只有2D信息

问题分析

在两帧各自的相机坐标系中，设P点的相机坐标系的坐标分别为$P_1$、$P_2$，如图1所示。

相机位姿估计示意图

图1 相机位姿估计示意图

其中：

$$P_{1}=\begin{pmatrix}
X_{1}\
Y_{1}\
Z_{1}
\end{pmatrix},P_{2}=\begin{pmatrix}
X_{2}\
Y_{2}\
Z_{2}
\end{pmatrix}$$

相机从第一帧移动到第二帧，其旋转矩阵设为R，平移向量设为t。这里的R表示的旋转是相对于第一帧的姿态改变量。那么有:

$$P_{2}=RP_{1}+t \tag{1}$$

相机位姿估计的最终目的就是要根据这个运动方程，求解出相机的运动，也就是R、t。注意这里的R、t并不是相机外参！

由于$P_1$、$P_2$的坐标未知，但其在像素坐标系的坐标$p_1$、$p_2$已知，因此，根据小孔成像模型寻找$p_1$、$p_2$和$P_1$、$P_2$之间的关系。

小孔成像模型

小孔成像模型的示意图如图2所示。设P点在世界坐标系下坐标为$P_w$在相机坐标系下坐标为$P_c$，物理成像平面对应坐标为 $P’$, 像素平面对应坐标为p, 有：

$$p=\begin{pmatrix}
u\
v
\end{pmatrix},P^{‘}=\begin{pmatrix}
X^{‘}\
Y^{‘}\
Z^{‘}
\end{pmatrix},P_c=\begin{pmatrix}
X_c\
Y_c\
Z_c
\end{pmatrix},P_w=\begin{pmatrix}
X_w\
Y_w\
Z_w
\end{pmatrix}$$

小孔成像模型示意图

图2 小孔成像模型示意图

根据小孔成像模型，有如下公式：

$$
\begin{pmatrix}
u\
v\
1
\end{pmatrix}=\frac{1}{Z_{c}}\begin{pmatrix}
f_{x} & 0 & c_{x}\
0 & f_{y} & c_{y}\
0 & 0 & 1
\end{pmatrix}\begin{pmatrix}
X_{c}\
Y_{c}\
Z_{c}
\end{pmatrix} \tag{2}$$

该公式给出了点P相机坐标系下坐标与像素平面坐标之间的关系。

2D-2D估计

设相机内参矩阵为K，即

$$
K=\begin{pmatrix}
f_{x} & 0 & c_{x}\
0 & f_{y} & c_{y}\
0 & 0 & 1
\end{pmatrix}
$$

将公式(2)代入公式(1)，则

$$K^{-1}Z_{2}p_{2}=RK^{-1}Z_{1}p_{1}+t \tag{3}$$

内参矩阵K短期内不变，但$Z_1$很明显一般情况下不等于$Z_2$，除非是相机绕着P点为圆心旋转。现在已知的是K、p1、p2，待求R，t。然而这里还有$Z_1$、$Z_2$是未知的，因此必须想办法将$Z$消去。

将相机坐标系下坐标$P_c$归一化，即

$$
P_{c}^{‘}=\begin{pmatrix}
{X_{c}}/{Z_{c}}\
{Y_{c}}/{Z_{c}}\
1
\end{pmatrix}
$$

由公式(2)可得：

$$
\begin{pmatrix}
u\
v\
1
\end{pmatrix}=\begin{pmatrix}
f_{x} & 0 & c_{x}\
0 & f_{y} & c_{y}\
0 & 0 & 1
\end{pmatrix}\begin{pmatrix}
{X_{c}}/{Z_{c}}\
{Y_{c}}/{Z_{c}}\
1
\end{pmatrix} =K\begin{pmatrix}
{X_{c}}/{Z_{c}}\
{Y_{c}}/{Z_{c}}\
1
\end{pmatrix} \tag{4}$$

将公式(4)代入公式(1)，则有：

$$K^{-1}p_{2}=RK^{-1}p_{1}+t$$

这样便能求解了。但是也导致了一个问题，因为不同帧对应的相机坐标系中P的Z值并不相等。在实际操作中分别除以其本身从而将Z分量归一化为1，这其实就丢失了真实的位置信息，不同帧的缩放是不等的！从而导致2D-2D估计两两帧之间，每次估计的尺度都是不同的。这也就是单目SLAM的尺度不确定性。

通过2D-2D，最终可获得P点像素坐标 $p$ 对应相机坐标系下归一化坐标$P_{c}^{‘}$。

三角测量

三角测量在三角学与几何学上是一借由测量目标点与固定基准线的已知端点的角度，测量目标距离的方法。而不是直接测量特定位置的距离（三边量测法）。当已知一个边长及两个观测角度时，观测目标点可以被标定为一个三角形的第三个点。

基于两固定角度之距离量测

假设一量测目标点及两个已知座标的参考点可形成一个三角形，则借由计算三角形其中参考边的长度，量测两参考点与目标点形成的角度，即可找出目标点的距离及座标。

基于两固定角度之距离量测示意图

图3 基于两固定角度之距离量测示意图

三角测量可用来计算岸边与船只之间的距离及座标。A顶点的观察者测量岸边与船只之间的角度α，B点的观察者则依同理测量出角度β，由长度l或已知的A及B点座标，则可由正弦定理取得在C点船只的座标及距离d。

计算过程如下：

$$ \ell = \frac{d}{\tan \alpha} + \frac{d}{\tan \beta}$$

根据三角恒等式${\displaystyle \tan \alpha ={\frac {\sin \alpha }{\cos \alpha }}}$和${\displaystyle \sin \left(\alpha +\beta \right)=\sin \alpha \cos \beta +\cos \alpha \sin \beta }$，此式可等于：

$${\displaystyle \ell =d\left({\frac {\cos \alpha }{\sin \alpha }}+{\frac {\cos \beta }{\sin \beta }}\right)}$$

$${\displaystyle \ell =d\ {\frac {\sin(\alpha +\beta )}{\sin \alpha \sin \beta }}}$$

因此，

$${\displaystyle d=\ell \ {\frac {\sin \alpha \sin \beta }{\sin(\alpha +\beta )}}}$$

由此便可简单定义出一未知点与观察点间的距离，以及与观察点往东西、南北向相差的位移量，终得完整座标。

2D-2D估计获得特征点在相机坐标系下的归一化3D坐标，结合三角测量，可获得特征点在相机坐标系下的深度，两者结合即获得特征点在相机坐标系下的3D坐标。

3D-2D估计

3D-2D估计本质是PnP一个问题，即给定世界坐标系中n个3D点及其在图像中的相应2D投影的情况下，估计校准相机的姿势的问题。相机姿势由6个自由度（DOF）组成，其由旋转（滚动，俯仰和偏航）以及相机相对于世界的3D平移构成。该问题源于相机校准，并且在计算机视觉和其他领域中具有许多应用，包括3D姿态估计、机器人和增强现实。对于n = 3，存在一个常用的问题解决方案，称为P3P，并且许多解决方案适用于n≥3的一般情况。

P3P

P3P仅需要使用三对匹配点，就可以完成相机的位姿估计。

假设空间中有A,B,C三点，投影到成像平面中有a,b,c三点，在PnP问题中，A,B,C在世界坐标系下的坐标是已知的，但是在相机坐标系下的坐标是未知的。a,b,c的坐标是已知的。PnP的目的就是要求解A,B,C在相机坐标系下的坐标值。如下图所示。需要注意的是三角形abc和三角形ABC不一定是平行的。

P3P问题示意图

图4 P3P问题示意图

根据余弦定理有：

$$OA^2 + OB^2 - 2OA \cdot OB \cdot \cos(a,b) = AB^2 \ OB^2 + OC^2 - 2OB \cdot OC \cdot \cos(b,c) = BC^2 \ OA^2 + OC^2 - 2OA \cdot OC \cdot \cos(a,c) = AC^2$$

记$x=\dfrac{OA}{OC}$,$y=\dfrac{OB}{OC}$，因为A,B,C在相机坐标系中的坐标未知，因此x，y是未知的。

另记$u=\dfrac{BC^2}{AB^2}$,$w=\dfrac{AC}{AB}$, 根据A,B,C的世界坐标，u,w是可以求出的。

通过一系列的转化可以得到两个等式：

$$(1-u)y^2-ux^2-\cos(b,c)y+2uxy \cos(a,b) +1 = 0 \ (1-w)x^2-wy^2-\cos(a,c)x+2wxy \cos(a,b) +1 = 0$$

该方程组是关于x,y的一个二元二次方程，可以通过吴消元法求解。最多可能得到四个解，因此在三个点之外还需要一组匹配点进行验证。

至此，通过x和y就可以求得A，B，C在相机坐标下的坐标值。因此3D-2D问题转变成了3D-3D的位姿估计问题。而带有匹配信息的3D-3D位姿求解非常容易。

Bundle Adjustment

假设某空间点坐标为$P_i = [X_i, Y_i, Z_i]$, 其投影的像素坐标为$p_i=[u_i,v_i]$。这些在PnP问题里都是已知的。在相机坐标系下有$c=[x_i, y_i, z_i]$，这个坐标通过P3P或者其他解法有了粗略的估计。根据针孔相机模型可得：

$$z_i p_i = KTP_i = K \exp([\xi]_{\times})P_i$$

根据这个等式可以构造出一个最小二乘问题：

$$\xi^* = \arg \min \limits {\xi} \dfrac{1}{2} \sum\limits _{i=1} ^n \begin{Vmatrix} p_i - \dfrac{1}{z_i} K \exp([\xi]{\times})P_i \end{Vmatrix} _2 ^2$$

该问题的误差项，是将像素坐标与3D点按照当前估计的位姿进行投影得到的位置相比较得到的误差，所以称之为重投影误差。如图5所示。

重投影误差示意图

图5 重投影误差示意图

这个最小二乘问题主要优化两个变量，第一是对相机位姿的优化，也就是对李代数的优化，第二是对空间点P的优化，也就是P点的优化。

3D-3D估计

3D-3D的位姿估计问题是指，对于空间中的某一点，我们知道这个点在两个相机坐标系中的三维坐标，如何利用这两个三维坐标来求解这两个相机坐标系的运动就是3D-3D的位姿评估问题。这个问题通常用迭代最近点(Iterative Closest Point,ICP)求解。

假设空间中的一系列点在第一个相机坐标系下的三维坐标为$C={c_1,…,c_n}$，在第二个相机坐标系下匹配的三维坐标为$C’={c_1’,…,c’_n}$。则有：

$$\forall i, \ c_i=Rc_i’+t$$

对于ICP的求解主要分为两种方式：利用线性代数的求解和利用非线性优化方式求解。

线性代数求解

构造误差项：

$$e_i = c_i - (Rc_i’+t)$$

将这个误差项构造成一个最小二乘问题：

$$\min \limits _{R,t} J= \dfrac{1}{2} \sum \limits _{i=1} ^n \begin{Vmatrix} c_i - (Rc_i’+t) \end{Vmatrix} _2 ^2$$

通过求解这个最小二乘问题，我们可以得到R和t。

总结

相机位姿估计最终目标是获得表征相机运动的旋转矩阵R和平移向量t，可分成两种方法：一种是2D-2D估计加三角测量，另一种是3D-2D估计加3D-3D估计。

参考链接

SLAM相机位姿估计(1),by zhaoxuhui.
2D-2D相机位姿估计,by 金戈大王.
单目相机中的对极几何,by 一索哥传奇.
三角测量,by 金戈大王.
三角测量,by wikipedia.
单目相机中的三角化测量,by 一索哥传奇.
3D-2D相机位姿估计,by 金戈大王.
3D-3D相机位姿估计,by 金戈大王.
相机位姿求解问题？, by zhihu.
3D-2D的运动估计,by 一索哥传奇.
3D-3D的运动估计,by 一索哥传奇.
相机位姿求解——P3P问题,by 达达MFZ.
图像二维坐标转世界三维坐标,by 橙子.

光线追踪基本概念入门

发表于 2019-05-26 | 更新于 2022-05-11

光线追踪（Ray tracing）是三维计算机图形学中的特殊渲染算法，跟踪从眼睛发出的光线而不是光源发出的光线，通过这样一项技术生成编排好的场景的数学模型显现出来。光线追踪的优点可以提供更为真实的光影效果，缺点是计算量巨大。

基本概念

光线追踪与光栅化渲染作为相对的两个概念，理解光栅化渲染更能解释光线追踪的概念。

光栅化渲染是将向量图形格式表示的图像转换成位图以用于显示器或者打印机输出的过程，如图1所示。

光线追踪示意图

图1 光栅化渲染示意图

光线追踪的示意图如图2所示。[Whitted 1980]提出了使用光线跟踪来在计算机上生成图像的方法，这一方法后来也被称为经典光线跟踪方法或Whitted-style 光线跟踪方法。其主要思想是从视点向成像平面上的像素发射光线，找到与该光线相交的最近物体的交点，如果该点处的表面是散射面，则计算光源直接照射该点产生的颜色；如果该点处表面是镜面或折射面，则继续向反射或折射方向跟踪另一条光线，如此递归下去，直到光线逃逸出场景或达到设定的最大递归深度。这种经典的方法可以产生镜面反射、折射、阴影等效果，不过不能实现其他的全局光照的效果。

光线追踪示意图

图2 光线追踪示意图

辐射度学基本量

图形学模拟可见光与各种材质的交互，这个过程涉及到能量的传输。辐射度学（Radiometry）是度量电磁辐射能量传输的学科，也是基于物理着色模型的基础。

能量

能量（Energy），用符号Q表示，单位焦耳（J），每个光子都具有一定量的能量，和频率相关，频率越高，能量也越高。

功率

功率（Power），单位瓦特（Watts），或者焦耳／秒（J/s）。辐射度学中，辐射功率也被称为辐射通量（Radiant Flux）或者通量（Flux），指单位时间内通过表面或者空间区域的能量的总量，用符号\Phi 表示，定义 $\Phi = \frac{ dQ}{dt}$。

辐照度和辐出度

辐照度（Irradiance），指单位时间内到达单位面积的辐射能量，或到达单位面积的辐射通量，也就是通量对于面积的密度。用符号E表示，单位 $W / m^{2}$ 。定义为 $E = \frac{d\Phi }{dA}$。

辐出度（Radiant Existance），也称为辐射出射度、辐射度（Radiosity），用符号M表示。辐出度与辐照度类似，唯一的区别在辐出度衡量的是离开表面的通量密度，辐照度衡量的是到达表面的通量密度。辐照度和辐出度都可以称为辐射通量密度（Radiant Flux Density）。

辐射强度

立体角则是度量三维角度的量，用符号 $\omega$ 表示，单位为立体弧度（也叫球面度，Steradian，简写为sr），等于立体角在单位球上对应的区域的面积（实际上也就是在任意半径的球上的面积除以半径的平方 $\omega = \frac {s}{r^{2} }$ ），单位球的表面积是 $4\pi$ ，所以整个球面的立体角也是$4\pi$ 。

辐射强度（Radiant Intensity），指通过单位立体角的辐射通量。用符号I表示，单位 $W / sr$，定义为 $I = \frac{d \Phi }{d \omega }$ 。之所以引入辐射强度，是因为有时候要度量通过一个点的通量的密度，但因为点的面积是0，无法使用辐照度，所以引入辐射强度。辐射强度不会随距离变化而变化，不像点光源的辐照度会随距离增大而衰减，这是因为立体角不会随距离变化而变化。

辐射率

辐射率（Radiance），指每单位面积每单位立体角的辐射通量密度。用符号 $L$ 表示，单位 $W/m^{2} sr$ ，定义为 $L = \frac{d \Phi }{d\omega d A^{\bot } }$ 。其中 $dA^{\bot}$是微分面积dA在垂直于光线方向的投影。

渲染方程（The Rendering Equation）

[Kajiya 1986]第一次将渲染方程引入图形学，使用它来解释光能传输的产生的各种现象。这一方程描述了场景中光能传输达到稳定状态以后，物体表面某个点在某个方向上的辐射亮度（Radiance）与入射辐射亮度等的关系。

$$L_o(x,w_o)=L_e(x,w_0)+\int_\Omega{f_r(x,w_i,w_0)L_i(x,w_i)cos\theta_idw_i}$$

其中，$L_o(x,w_o)$ 表示物体表面点 $x$ 处在方向 $𝜔_𝑜$ 上出射的辐射亮度，$𝐿𝑒(𝑥, 𝜔_0)$ 表示在该点该方向上自辐射的亮度。$𝐿𝑖(𝑥,𝜔_𝑖)$ 表示该点处 $𝜔_𝑖$ 方向入射的辐射亮度，$𝑓𝑟(𝑥, 𝜔_𝑖, 𝜔_𝑜)$ 是双向反射分布函数（BRDF），描述的是入射方向的辐射亮度对出射方向的贡献，$cos 𝜃𝑖$ 是$𝜔_𝑖$与表面法向的点积。在这一方程的基础上，辐射度方法和蒙特卡罗光线跟踪的方法就可以看成是对方程中积分的不同的数值求解方法。

BRDF描述的是表面本身的性质，比如它的光滑程度，导电程度等等。由于四面八方的光线都会作用在这个表面，所以我们需要对所有方向进行积分，也就是一个球面上的积分，考虑到积分项中的 $(w_i\cdot n)$ ，那么只有位于正半空间的方向才会对最终积分有贡献，所以最后这个球面的积分就变成了一个半球的积分，如图3所示。

光线追踪示意图

图3 BRDF示意图

双向反射分布函数（bidirectional reflectance distribution function, BRDF）

在计算机图形学领域，着色（Shading）是指根据表面或者多边形相对光源和相机的角度和距离来计算它的颜色的过程。不同的用途可以使用不同的着色算法，CAD等追求响应速度的交互式图形领域可以使用简单快速的着色算法，卡通油画等艺术效果可以使用非真实感（Nonphotorealistic）着色算法，而追求真实感的CG电影或游戏则可以使用基于物理建模的着色算法。而BRDF是基于物理建模的着色算法的理论基础。

我们看到一个表面，实际上是周围环境的光照射到表面上，然后表面将一部分光反射到我们眼睛里。双向反射分布函数BRDF（Bidirectional Reflectance Distribution Function）就是描述表面入射光和反射光关系的。

对于一个方向的入射光，表面会将光反射到表面上半球的各个方向，不同方向反射的比例是不同的，我们用BRDF来表示指定方向的反射光和入射光的比例关系，BRDF定义为：

$$f(l,v)=\frac{dL_o(v)}{dE(l)}$$

其中，$f$就是BRDF，$l$是入射光方向，$v$是观察方向，也就是我们关心的反射光方向。$d L_o(v)$ 是表面反射到$v$方向的反射光的微分辐射率。表面反射到$v$方向的反射光的辐射率为$L_o(v)$，来自于表面上半球所有方向的入射光线的贡献，而微分辐射率 $d L_o(v)$ 特指来自方向 $l$ 的入射光贡献的反射辐射率。$dE(l)$是表面上来自入射光方向 $l$ 的微分辐照度。表面接收到的辐照度为 $E$ ，来自上半球所有方向的入射光线的贡献，而微分辐照度 $dE(l)$ 特指来自于方向 $l$ 的入射光。

光照模型（illumination model）

当光照射到物体表面时,物体对光会发生反射、透射、吸收、衍射、折射、和干涉,其中被物体吸收的部分转化为热,反射、透射的光进入人的视觉系统,使我们能看见物体。为模拟这一现象,我们建立一些数学模型来替代复杂的物理模型，这些模型就称为明暗效应模型或者光照明模型。

局部光照模型

在真实感图形学中,仅处理光源直接照射物体表面的光照明模型被称为局部光照明模型。局部光照明模型的分类如图4所示。

局部光照模型分类

图4 局部光照模型分类

局部光照模型是一种比较简单的光照模型，它是与光栅化渲染算法相适应的，光栅化算法一次只考虑一个像素的光照强度，因此局部光照模型不能计算某像素受其他像素影响的光照强度部分。也就是说，局部光照模型只对物体进行直接光照的计算，而不考虑其他的间接影响。

全局光照模型

全局光照模型是基于光学物理原理的，光照强度的计算依赖于光能在现实世界中的传播情况，考虑光线与整个场景中各物体表面及物体表面间的相互影响，包括多次反射、透射、散射等。因此，与局部光照模型相比，全局光照模型需要相当大的计算量，但同时也能取得非常逼真的真实效果。全局光照模型分类如图5所示。

全局光照模型分类

图5 全局光照模型分类

求交检测

在光线追踪过程中，从眼睛发出的光线与3D模型的三角面求交是一个复杂问题。通常精致的3D模型可能由几十万至上百万三角面构成，如果采用穷举法求交点，其时间复杂度将是O(n), 过于复杂。为了减少不必要的求交检测，应采用空间划分技术，最常用的是平衡kdtree算法，提高求交检测的效率。

k-d tree

在计算机科学里，k-d树（ k-维树的缩写）是在k维欧几里德空间组织点的数据结构。k-d树可以使用在多种应用场合，如多维键值搜索（例：范围搜寻及最邻近搜索）。k-d树是空间二分树（Binary space partitioning ）的一种特殊情况。

k-d树是每个节点都为k维点的二叉树。所有非叶子节点可以视作用一个超平面把空间分割成两个半空间。节点左边的子树代表在超平面左边的点，节点右边的子树代表在超平面右边的点。选择超平面的方法如下：每个节点都与k维中垂直于超平面的那一维有关。因此，如果选择按照x轴划分，所有x值小于指定值的节点都会出现在左子树，所有x值大于指定值的节点都会出现在右子树。这样，超平面可以用该x值来确定，其法线为x轴的单位向量。

参考链接

光线追踪,by wikipedia.
栅格化,by wikipedia.
渲染,by wikipedia.
光线追踪基本概念与代码实现,by 鹅城惊喜师爷.
一篇光线追踪的入门,by 洛城.
基于蒙特卡罗的光线跟踪绘制方法,by 严俊.
基于物理着色：BRDF,by Maple.
Monte-Carlo Ray Tracing System (一)原理以及设计,by 已退逼乎.
冯氏光照模型–镜面光的计算,by MooAiFighting.
什么是光照模型,by 黄琦.
蒙特卡洛光线追踪,by sunacmer.
k-d树,by wikipedia.
辐射强度、辐亮度、辐照度——一文搞定,by 三眼二郎.
从光栅化到光线追踪,by CrazyEngineCo.
windows下没有srand48和drand48的解决方法,by 查志强.
辐射照度、辐射强度、光照度、发光强度（差异以及如何相互转换）（易懂讲解）,by 三眼二郎.