计算机视觉

目标：让机器能够“理解”图像的内容：

图像的语义：识别
图像的三维：重建

图形经过建模，利用计算机图形学渲染（绘制）得到图像；
图像经过图像处理，利用计算机视觉得到图形。

Introduction

图像表示

程序中的图像表示

struct MyImage
{
     int width, height; // 大小
     int type; // 类型，含通道数、位深度信息
		       /* CV_8UC3 : unsigned char [3]
		          CV_32SC1: int [1]
                  CV_32UC1: uint [1]
                  CV_32FC4: float [4]
               */
     void* data; // 图像数据
     int step; // 步长(每行所占用的字节数):用于数据对齐；表示子区域(ROI)
};

访问像素

img.type = CV_8UC3; // 8位无符号，3通道数据
uchar* get_pixel(const MyImage &img, int x, int y)
{
    return (uchar*)img.data + y * img.step + x * 3;
}

img.type=CV_32SC3; // 32位带符号，3通道数据
int*  get_pixel(const MyImage &img, int x, int y)
{
    return (int*)((char*)img.data + y * img.step + x * 3 * 4);
}

遍历像素

void scan_pixels(uchar *data, int width, int height, int step, int nc)
{
    uchar *row = data;
    for (int yi = 0; yi < height; ++yi, row += step)
    {
        uchar *px = row;
        for (int xi = 0; xi < width; ++xi, px += nc)
        {
            // px now address the pixel (xi, yi)
        }
    }
}

void scan_roi_pixels(MyImage &img, int x, int y, int roi_width, int roi_height)
{// 通道数nc = img.nc();
    scan_pixels(get_pixel(img, x, y), roi_width, roi_height, img.step, img.nc());
}

Image Processing

像素灰度变换
亮度变换
基本灰度变换函数
alt text

Log Transformation： $s=c\cdot\log(1+r)$ ， $c$ 为常数
Power-Law Transformations： $s=cr^\gamma$ $s = c r^{γ}$ ， $c,\gamma$ $c, γ$ 为正常数
- $s=cr^\gamma$ 变换函数 $c=1, \gamma$ 取不同值下的图像

对比度变换
增加灰度的变化范围
alt text

Quantize & Threshold 量化和阈值
分段函数
sigmoid函数

代数运算：多幅图像
xor or and
alpha混合： $C=\alpha F +(1-\alpha)B$
背景相减：两幅图像像素相减，差值大于某阈值的位置作为前景，即可恢复出前景。

基本几何处理
水平&垂直翻转 (flip)

// 垂直翻转
void vflip(const void *in, int width, int height, int istep, int pix_size, void *out, int ostep)
{
    out = (char*)out + (height - 1) * ostep;
    for(int yi = 0; yi < height; ++yi, in = (char*)in + istep, out = (char*)out - ostep)
    {
        memcpy(out, in, width * pix_size);
    }
}
// 水平翻转
void hflip(const void *in, int width, int height, int istep,  int pix_size, void *out, int ostep)
{
    char * _in = (char*)in;
    char *_out = (char*)out+(width-1)*pix_size;
    
    for(int yi = 0; yi < height; ++yi, _in += istep, _out += ostep)
    {
        char *in_x = _in, 
        char *out_x = _out;
        for(int xi = 0; xi < width; ++xi, in_x += px_size, out_x -= px_size)
            memcpy(out_x, in_x, px_size) ;
    }
}

缩放 (resize / zoom in /zoom out /scale)
实现方式：

Projection (from source to target)：对小图中的每个像素，计算其在大图中对应的像素，再拷贝小图的像素值到大图。
- 问题：大图中的某些像素可能无法得到填充

Lookup (from target to source)：对大图中的每个像素，计算其在小图中对应的像素，再拷贝小图的像素值到大图。

问题：计算得到的在小图中的位置不是整数

解决：

重采样Resampling

最近邻 (Nearest Neighbor)
- x=int(x + 0.5), y=int(y + 0.5);

双线性插值 (Bilinear Interpolation)

两水平一垂直

float bilinear(float a, float b, float c, float d, float dx, float dy)
{
    float  h1 = a + dx * (b - a); // = (1 - dx) * a + dx * b
    float  h2 = c + dx * (d - c);
    return h1 + dy * (h2 - h1);
}

两垂直一水平

float bilinear(float a, float b, float c, float d, float dx, float dy)
{
    float  w1 = a + dy * (c - a); // = (1 - dy) * a + dy * c
    float  w2 = b + dy * (d - b);
    return w1 + dx * (w2 - w1);
}

两种方法等价

双三次插值 (Bicubic Interpolation)

超分辨率Super-Resolution

变换矩阵

仿射变换
$\begin{bmatrix} x'\\ y'\\ 1 \end{bmatrix} = \begin{pmatrix} a & b & c\\ d & e & f\\ 0 & 0 & 1 \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$
- 平移 $\begin{bmatrix} x'\\ y'\\ 1 \end{bmatrix} = \begin{pmatrix} 1 & 0 & d_x\\ 0 & 1 & d_y\\ 0 & 0 & 1 \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$
- 缩放 $\begin{bmatrix} x'\\ y'\\ 1 \end{bmatrix} = \begin{pmatrix} s_x & 0 & 0\\ 0 & s_y & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$
- 旋转：左手系，顺时针旋转 $\theta$ $\begin{bmatrix} x'\\ y'\\ 1 \end{bmatrix} = \begin{pmatrix} \cos\theta & -\sin\theta & 0\\ \sin\theta & \cos\theta & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$
- x切变 $\begin{bmatrix} x'\\ y'\\ 1 \end{bmatrix} = \begin{pmatrix} 1 & s & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$
- y切变 $\begin{bmatrix} x'\\ y'\\ 1 \end{bmatrix} = \begin{pmatrix} 1 & 0 & 0\\ s & 1 & 0\\ 0 & 0 & 1 \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$
透视变换 $\begin{bmatrix} x'\\ y'\\ w \end{bmatrix} = \begin{pmatrix} a & b & c\\ d & e & f\\ g & h & 1 \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$

仿射变换 (Affine Transform)：线性变换+平移变换，在同一平面内变换

保持共线性、共线向量长度比、重心坐标
特殊的仿射变换
- 刚性变换(Rigid Transformation)：
  - 只包含平移和旋转
  - 保持物体的形状(保角)和尺寸
  - 相当于正交变换
- 相似变换(Similar Transformation)：
  - 只包含平移、旋转和等比缩放
  - 保持物体的形状

透视变换（Perspective Transform）：可表示不同视角观察到的同一平面，或同一视角观察到的不同平面之间的变换

仿射变换应用：图像匹配
不共线的3个平面点对决定一个二维仿射变换

\begin{bmatrix} u_i'\\ v_i' \end{bmatrix} = \begin{pmatrix} a & b & c\\ d & e & f \end{pmatrix} \begin{bmatrix} u_i\\ v_i\\ 1 \end{bmatrix} \implies \begin{cases} u_0a+v_0b+c=u_0'\\ u_0d+v_0e+cf=v_0'\\ u_1a+v_1b+c=u_1'\\ u_1d+v_1e+cf=v_1'\\ u_2a+v_2b+c=u_2'\\ u_2d+v_2e+cf=v_2'\\ \end{cases}

超过3个点对
$A=\argmin_A \sum_i\parallel Ap_i-p_i'\parallel^2$
基于仿射变换的图像匹配

在第t帧检测特征点（特征检测）
计算特征点在第t+1帧的对应（特征跟踪）
根据特征点对估计第t帧到第t+1帧的仿射变换A
利用A对第t帧的图像进行变换，将变换的结果作为与第t+1帧配准的图像： $Ap_t\leftrightarrow p_{t+1}$

正向查找/逆向查找应用：图像变形 (Image Warping)
坐标归一化，获取原图中对应位置的像素

Spatial Filtering

空间域滤波器

线性滤波器： $f$ $f$ 是线性函数
- $p'=\sum_iw_ip_i$
  $p'(x,y)=\sum_{s=-a}^a\sum_{s=-b}^b w(s,t)p(x+s,y+t)$ ， $w(s,t)$ 为滤波核
- 均值滤波器： $p'=\frac{1}{9}(p_0+p_1+p_2+\cdots+p_8)$
非线性滤波器： $f$ $f$ 是非线性函数
- 最大值滤波器： $p'=\max(p_0,p_1,p_2\dots p_8)$
- 最小值滤波器： $p'=\min(p_0,p_1,p_2\dots p_8)$
- 中值滤波器
- 双边滤波

线性滤波器

平滑滤波器(Smoothing Filter)(低通滤波器)

原理：加权平均（均值滤波，高斯滤波）
边界处理
- 原因1：对边界附近的像素，滤波核的部分可能落在图像区域外。
- 原因2：如果不进行边界处理，不同像素点作为滤波核中心的次数不一样，有些点不会作为滤波核的中心出现
- 实现：
  - 对边界外某个范围内的区域进行填充
    - 常数填充
    - 镜像填充（以边界为轴，取图像内对称点的像素值填充）
  - 在边界附近调整滤波核的大小

快速均值滤波

积分图：图像I的积分图S是与其大小相同的图像，S的每一像素S(u,v)存贮的是I(u,v)左上角所有像素的颜色值之和。（二维前缀和）

void integral_image(const uchar *src, int width, int height, int sstride, int *pint, int istride)
{
    int *prow = new int[width];

    memset(prow, 0, sizeof(int)*width);

    for(int yi = 0; yi < height; ++yi, src += sstride, pint += istride)
    {
        prow[0] += src[0];  pint[0] = prow[0]; // for the first pixel
        
        for(int xi = 1; xi < width; ++xi)
        {
            prow[xi] += src[xi];
            pint[xi] = pint[xi - 1] + prow[xi];
        }
    }
    delete[]prow;
}

高斯滤波
- 以高斯函数为滤波核 $G(x,y)=\frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}$
- 高斯滤波具有行列可分离性，核大小为M的二维图像高斯滤波，等价于同样核大小的一维高斯滤波在行列方向的叠加。所以在对二维图像进行高斯滤波时，为了减少运算时间，可以将行列单独进行滤波后叠加，将二维滤波转化为两个一维滤波。
- $\sigma$ 越大，核就应该越大， $size=[6\sigma-1]$

锐化滤波器(Sharpening Filter)(高通滤波器)
- 基本高通滤波器
  - 滤波器在中心有正系数，而在边缘上有负系数，且总和为0
  - 在平坦变化的区域很暗，在剧烈变化的区域很亮
  - 可以用来边缘检测，滤波之后，只有边缘区域是亮的
- 导数滤波器
  - 一阶导： $f'(x)=f(x+1)-f(x)$ $f^{'} (x) = f (x + 1) - f (x)$
    - 零：在数值不变的区域
    - 非零：在数值变化的区域
  - 二阶导： $f''(x)=f(x+1)+f(x-1)-2f(x)$ $f^{''} (x) = f (x + 1) + f (x - 1) - 2 f (x)$
    - 零：数值均匀变化的区域
    - 非零：数值剧变的区域
  - 二阶导在边缘处过零点，可以简单确定边缘。但二阶导对噪声敏感，因此需要先进行高斯滤波去噪
  - 一阶导的运用：梯度
    - 图像处理中最常见的差分方法就是梯度 $\nabla f_{x, y}=(G_x, G_y)=(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y})$ ，梯度就是一个向量，指向灰度变化最大的方向，其长度为： $mag(\nabla f)=[G_x^2+G_y^2]^{\frac{1}{2}}=[(\frac{\partial f}{\partial x})^2+(\frac{\partial f}{\partial y})^2]^\frac{1}{2}$
    - Roberts算子
    - Prewitt算子
    - Sobel算子
      
      Sobel算子优点：Sobel算子跟高斯核相似，在中间的像素有较大权重，在两侧的像素的权重较小，对于噪声的抗干扰能力比较强，可以有效地去除一些小的噪点。Sobel算子的计算量比较小，处理速度比较快。Sobel算子的检测结果比较清晰，边缘位置比较精准。
  - 二阶导的运用：拉普拉斯
    - Laplacian定义为x,y方向上的二阶导的和： $\nabla^2f=\frac{\partial^2 f}{\partial x^2}+\frac{\partial^2 f}{\partial y^2}$ ，离散表示： $\frac{\partial^2 f}{\partial x^2}=f(x+1,y)+f(x-1,y)-2f(x,y),\frac{\partial^2 f}{\partial y^2}=f(x,y+1)+f(x,y-1)-2f(x,y)$
    - 拉普拉斯滤波核可以中间是正数周围是负数，也可以中间是负数周围是正数
    - 可以通过将拉普拉斯检测出来的图像加上原图像来锐化原图像
    - 高斯拉普拉斯算子(Laplacian of Gaussian,LoG)：高斯函数的二阶导（先高斯滤波再求导等价于先对高斯函数求导再滤波）

图像的线性滤波是对图像的线性变换，可以表示成对图像的矩阵变换形式

非线性滤波器

次序统计滤波器(Order Statistics Filters)

中值滤波器：取邻域中灰度的中值作为该点的灰度值
- 尺寸越大图像越模糊
- 可用于消除椒盐噪声
最大值滤波器
最小值滤波器
快速最大/最小值滤波
1. 每次选取一个区域M，将M的中间位置当作c
2. 对位置i来说，若i在c左边，则位置i的值为从i到c的最大/最小值；若i在c右边，则位置i的值为从c到i的最大/最小值
3. 由2得到一个序列T，T从c往左往右都是单调不减/增的
4. 求一个范围内的最大/最小值，只需要比较该范围的左右端点即可

双边滤波器
高斯滤波的问题：消除噪音的同时也使边缘变得模糊。因为高斯滤波在滤波过程中只关注了位置信息，即在滤波窗口内，距离中心点越近的点权重越大。这在像素值出现剧变的边缘区域，会损失有用的边缘信息。
双边滤波：计算权重时同时考虑空间位置和像素颜色之差，在高斯滤波的基础上加入了像素值权重项。在图像的平坦区域，像素值变化很小，对应的像素范围域权重接近于1，此时空间域权重起主要作用，相当于进行高斯模糊；在图像的边缘区域，像素值变化很大，像素范围域权重变大，从而保持了边缘的信息。
alt text

Matching

特征检测

目标：稳定检测；易于匹配
挑战：对图像几何变换的稳定性；对颜色、光照变化的稳定性
角点检测(Corner Detection)
角点：在任意方向像素颜色都会有明显变化。等价于：在角点附近，图像梯度具有至少两个主方向。
Harris角点检测
将图像区域内的一个窗口进行小位移移动，观察窗口内的像素变化，用以下公式来衡量像素的变化： $E(u,v)=\sum_{x,y}w(x,y)[I(x+u,y+v)-I(x,y)]^2$ ， $w(x,y)$ 是窗口的权重函数（可以是高斯函数，也可以全部设置为1）， $I(x+u,y+v)$ 是移动后的像素值， $I(x,y)$ 是移动前的像素值。将 $I(x+u,y+v)$ 项泰勒展开，再将平方项展开写成矩阵形式，经过一系列推导，可以将 $E(x,y)$ 写成以下形式：

E(u,v)\approx\begin{bmatrix} u & v \end{bmatrix} M \begin{bmatrix} u\\ v \end{bmatrix}\\ M=\sum_{x,y}w(x,y) \begin{bmatrix} I_x^2 & I_xI_y\\ I_xI_y & I_y^2 \end{bmatrix}

这也是椭圆的矩阵表示形式，假设 $M$ 的特征值为 $\lambda_1,\lambda_2$ ，那么分为以下三种情况：

$\lambda_1,\lambda_2$ 都很小，且二者比值近似，说明 $E$ 在各个方向上几乎都没有变化，这表示图像中的平面
$\lambda_1,\lambda_2$ 一大一小，且二者相差较大，这表示图像中的边缘
$\lambda_1,\lambda_2$ 都很大，且二者比值近似，说明 $E$ 在各个方向上都有变化，这表示图像中的角点

为了简化计算，定义了响应函数 $R=\det (M)-\alpha \text{trace}(M)^2=\lambda_1\lambda_2-\alpha(\lambda_1+\lambda_2)^2$ ， $R$ 为响应值，角点的响应值较大，平坦区域的响应值较小，边缘的响应值为负数； $\alpha$ 为经验常数，通常为0.04~0.06。
得到角点之后进行非极大值抑制操作。

Harris对光照变化和几何变换的稳定性？
Harris具有旋转不变性和光照不变性，不具有尺度不变性

旋转不变性：图像旋转后，角点检测过程中得到了椭圆只是方向发生变化，但形状没变，响应函数R没变，故具有旋转不变性
不具有尺度不变性：尺度变化会将角点变成边缘
光照不变性：整体光照改变了像素值的大小，不改变梯度

斑点检测(Blob Detection)
利用图像的Hessian矩阵（图像某一点处的二阶导数），在两个正交方向上寻找强梯度。使用公式 $\det(H)=I_{xx}I_{yy}-I_{xy}^2=\lambda_1\lambda_2$ ，结果：响应主要在角点和纹理比较强的区域。

如何实现尺度不变性？
分辨率越高，尺度越小；分辨率越低，尺度越大
高斯金字塔
高斯核是实现尺度变换的唯一核
通过降采样构建高斯金字塔。降采样前需要先对图像进行平滑以避免走样。
alt text
每一层都通过对前一层的平滑和下采样来得到
拉普拉斯金字塔

拉普拉斯金字塔可以看作是是对高斯的差分(Difference of Gaussian,DoG)，拉普拉斯金字塔的第i层是高斯金字塔的第i层减去第i-1层的扩展（插值+滤波）得到的。

得到金字塔后，在不同尺度同时进行检测和匹配，就可以实现尺度不变性
实现：

穷举搜索：利用多个大小的窗口对图像进行扫描
- 缺点：
  - 效率低
  - 不利于识别
自动尺寸选择
- 思想：设计一个尺度无关的函数。取函数的局部最大值，最大值对应的区域对于图像尺寸来说是不变的，最后归一化到固定尺寸
- 实现：使用LoG（LoG的 $\sigma$ $σ$ 与检测框的尺寸正相关，因此可以使用不同 $\sigma$ $σ$ 的LoG进行不同尺度的特征检测），然后记录LoG在尺度空间上的局部最大值对应的位置。实际使用时常用DoG来代替LoG，计算更方便。具体流程是：
  - 在尺度空间上检测DoG的局部最大值
  - 非极大值抑制
  - 消除边缘响应（由于某些原因，在检测极值点时，边缘容易被检测出来）

特征匹配

将特征所在的局部图像块转化为向量表示的特征描述，然后进行匹配。
特征描述要具有不变性（对一张图片的两次检测结果应该相同）和区分性（不同特征之间应该是可区分的）。
特征匹配的方法：

SSD (Sum of Square Difference)： $ssd(P,Q)=\sum_{(x,y)}|P(x,y)-Q(x,y)|^2$
SAD (Sum of Absolute Difference)： $sad(P,Q)=\sum_{(x,y)}|P(x,y)-Q(x,y)|$
NCC (Normalized Cross Correlation)： $ncc(P,Q)=\langle \frac{P-\overline{P}}{\parallel P-\overline{P}\parallel}, \frac{Q-\overline{Q}}{\parallel Q-\overline{Q}\parallel}\rangle$ （两个向量点积）
Census Transform：对图像块的二值编码，比中间像素小的视为0，比中间像素大的视为1，然后将两个二值序列异或，异或值越小则认为越相似。 $ct(P,Q)=CT(P) xor CT(Q)$
BRIEF

考虑对齐误差
以上方法都是对对应位置的像素进行计算，这样做的问题是一些微小的变动可能对匹配得分产生重大影响。因此需要统计整个图像块内的信息作为描述。

SIFT(Scale Invariant Feature Transform)
- 流程：
  - 构造DoG尺度空间
  - 在各尺度上定位关键点
  - 为关键点分配方向角
  - 形成特征描述符
    - 描述符的计算：将一个区域分成 $4\times 4$ 的子区域块，对于每一块子区域计算梯度方向直方图（8个参考角度），总共 $4\times 4\times 8=128$ 维。
PCA-SIFT
- 对SIFT的描述符计算进行了改进
SURF(Speeded Up Robust Features)
- 在特征点周围取一个 $4\times 4$ $4 \times 4$ 矩形区域块，但是取的矩形方向是沿着特征点的主方向。每个子区域统计25个像素的水平方向和垂直方向的Haar小波特征，这里的水平方向和垂直方向都是相对主方向而言的。该Haar小波特征为水平方向值之和、垂直方向值之和、水平方向绝对值之和、垂直方向绝对值之和。总共 $4\times 4\times 4=64$ $4 \times 4 \times 4 = 64$ 维
  - Haar特征（哈尔特征）：本例中使用Haar特征的两个模板。模板是矩形框，由黑色部分和白色部分构成。使用方法为：利用白色部分像素和减去黑色部分像素和。Haar特征反映了图像的灰度变化情况
  - Haar特征在灰度分布均匀的区域特征值趋近于0

如何实现旋转不变性？
以图像块中的主方向作为定位（统计图像块中像素的梯度，哪个梯度方向的像素最多就作为主方向），通过主方向对图像快进行旋转，这将使图像块旋转到一个规定的方向。

SIFT主方向选择：基于梯度直方图的主方向和特征描述
- 统计梯度直方图
- 选择最显著的梯度方向
- 归一化：将图像旋转到指定方向
SURF主方向选择：基于小波梯度的主方向和特征描述
- 统计特征点圆形邻域内的Harr小波特征。即在特征点的圆形邻域内统计 $60\degree$ 扇形内的所有水平、垂直Harr小波特征总和，然后扇形以0.2弧度大小的间隔进行旋转并再次统计该区域内的小波特征，最后将值最大的扇形的方向作为该特征点的主方向。
ORB(Orented-FAST and Rotated BRIEF)主方向选择：基于亮度质心的主方向+基于BRIEF的特征描述
- 亮度质心的计算同时考虑了像素位置和像素值大小

特征匹配：向量搜索

BruteForce
Hash-based
- Local Sensitive Hash (LSH)
Tree-based
- Approximate Nearest Neighbor (ANN)
- FLANN library

特征匹配的基本原理

尺度不变
旋转不变
光照不变
对齐误差

全局特征表示：对图像或者目标的整体进行编码表示，主要面向高层语义相关的任务，如识别、检测、检索等。

颜色直方图：
- 统计图像的颜色直方图，并将直方图作为图像的特征表示
- 缺点：只能表达图像的整体颜色分布，不能表示语义信息。
梯度方向直方图(Histogram of Oriented Gradients, HoG)
- 将图像均匀划分成 $8\times8$ 像素的区域，并计算每个区域内的梯度方向分布。将所有区域的梯度方向直方图连接为一个向量作为图像的特征表示
- 缺点：利用了目标的形状信息，但是对形状变化敏感。
词袋模型(Bag of Words, BoW)
- 将目标表示为局部特征（字典单词）的分布
  - 从训练集提取局部特征，并聚类获得字典
  - 提取图像特征，并与字典中的特征进行匹配，计算分布直方图
- 优点：具有较强的语义表达能力；对目标的形变不敏感
- 缺点：丢失了局部特征的空间信息
DPM(deformable parts model)
- 局部可形变，兼有整体的空间关系约束
卷积神经网络
- 越靠后的卷积层，语义表达能力越强，可以将最后一层的输出作为图像的全局特征表示

局部特征表示：对局部图像块进行编码，是图像匹配、三维重建等的基础。
全局特征表示：对整幅图像进行编码，是图像识别、检索等的基础。

运动估计

忽略外点(Outlier)：

齐次坐标和变换矩阵
- 仿射变换：线性变换+平移，改变物体的位置和形状，但保持平直性
  常见的变换矩阵：平移、缩放、旋转、x切变，y切变（前文有提到）
  不共线的三个平面点对确定一个二维仿射变换；大于三个点对可以通过优化的方式估计
- 特殊仿射变换
  - 相似变换：只包含平移、旋转、缩放，保持物体的形状
    $\begin{bmatrix} x'\\ y' \end{bmatrix} = \begin{pmatrix} s\cos\theta & -s\sin\theta & t_x\\ s\sin\theta & s\cos\theta & t_y \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$
  - 刚性变换：只包含平移和旋转，保持物体的形状(保角)和尺寸。相当于正交变换。
    
    $\begin{bmatrix} x'\\ y' \end{bmatrix} = \begin{pmatrix} \cos\theta & -\sin\theta & t_x\\ \sin\theta & \cos\theta & t_y \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$
    - 如何求解：非线性最小二乘
- 透视变换
  - 平面到平面的保持直线性的映射
  - 任意一个 $3\times3$ 可逆矩阵都是透视变换
  - 任意透视变换都可表示为 $3\times3$ 可逆矩阵
  - 中心投影对应的平面映射是透视变换
    - 场景（三维）中任意平面到图像平面的映射
    - 场景中同一平面在不同视点下图像之间的对应点的映射
    - 旋转相机在不同角度得到的图像之间的映射
- 透视变换估计
  
  $\begin{bmatrix} x'\\ y'\\ w \end{bmatrix} = \begin{pmatrix} a & b & c\\ d & e & f\\ g & h & 1 \end{pmatrix} \begin{bmatrix} x\\ y\\ 1 \end{bmatrix}$
  
  8个参数，最少需要4个点对
  - 非线性最小二乘
  - 直接线性变换

考虑外点：

RANSAC(RANdom SAmple Consensus)
RANSAC，随机抽样一致算法。是一种在包含离群点在内的数据集里，通过迭代的方式估计模型的参数的方法。有一定概率得到一个合理的结果。
思路：规避外点的影响，只使用内点。如果一个离群点被选择来计算当前的拟合，那么这个结果对剩下的点就不会有很好的拟合结果。
流程：
- 随机选择一组种子点（随机选取的点默认是内点）来进行基本的变换估计
- 利用这一组种子点计算变换公式（利用随机选择的局内点拟合一个模型）
- 找到符合这一变换公式的点，将其标注为内点
- 如果内点的数量足够大，那么重新通过所有内点重新计算上面得到的模型的最小二乘估计，来评估拟合出来的模型
- 如果当前模型的效果比当前最好模型的效果好，则更新最好模型，否则，抛弃，重新开始迭代。
需要多少次采样？
- w是内点的比例
- n个点能定义一个模型（直线需要2个）
- 进行了k次采样，每次取n个点
一次取样中，选出n个点全为内点的概率为 $w^n$ ，k次取样都没有全部取完n个内点的概率是 $(1-w^n)^k$ 。因此选择足够大k使其低于期望失败率。

总结：
- RANSAC将数据划分为内点和外点，并从内点的最小集合进行估计。
- 通过对所有内点的估计优化初始估计（例如标准的最小二乘法）。
- 但这可能会改变内点，因此会重新拟合并分类内点和外点。
问题：
- 离群点的比例可能很大（90%以上）
- 离群点的比例未知

Reconstruction

三维感知

单张图像：不准确
基于双目/多目：物体的深度值越大，在两幅图之间的视觉差越小。我们可以基于这个来判断物体在三维空间中的深度。

三角化(Triangulation)
通过两张图片对应点与相机所在直线，来确定三维空间中物体所在的位置。
alt text
进行三角化时需要知道的信息：

外参：相机的位置、朝向
内参：焦距等参数
像素对应关系

针孔相机模型
alt text
$(x,y)$ 是成像面上成像点的坐标， $(X,Y,Z)$ 是空间中的三维点的坐标， $f$ 是成像面距离相机光心的距离。根据点在三维空间中的坐标计算出该点在屏幕上的坐标。矩阵形式的表达：

\begin{pmatrix} x\\ y\\ 1 \end{pmatrix} = \begin{bmatrix} f & & \\ & f & \\ & & 1 \end{bmatrix} \begin{bmatrix} 1 & & & 0\\ & 1 & & 0\\ & & 1 & 0 \end{bmatrix} \begin{pmatrix} X\\ Y\\ Z\\ 1 \end{pmatrix}

相机坐标系

主轴(Principal axis)：从中心出发，与图像平面垂直
主点(Principal point)：主轴与图像平面的交点，理想情况在图像中心 $p=(0, 0)$

主点偏移
主点不在图像中心
alt text

\begin{pmatrix} x\\ y\\ 1 \end{pmatrix} = \begin{bmatrix} f & & p_x \\ & f & p_y\\ & & 1 \end{bmatrix} \begin{bmatrix} 1 & & & 0\\ & 1 & & 0\\ & & 1 & 0 \end{bmatrix} \begin{pmatrix} X\\ Y\\ Z\\ 1 \end{pmatrix}

像素长宽比
CCD单元长宽比不为1：pixel size: $\frac{1}{m_x}\times\frac{1}{m_y}$

\begin{aligned} \begin{pmatrix} x\\ y\\ 1 \end{pmatrix} &= \begin{bmatrix} m_x & & \\ & m_y & \\ & & 1 \end{bmatrix} \begin{bmatrix} f & & p_x \\ & f & p_y\\ & & 1 \end{bmatrix} \begin{bmatrix} 1 & & & 0\\ & 1 & & 0\\ & & 1 & 0 \end{bmatrix} \begin{pmatrix} X\\ Y\\ Z\\ 1 \end{pmatrix}\\ &= \begin{bmatrix} f_x & & p_x \\ & f_y & p_y\\ & & 1 \end{bmatrix} \begin{bmatrix} 1 & & & 0\\ & 1 & & 0\\ & & 1 & 0 \end{bmatrix} \begin{pmatrix} X\\ Y\\ Z\\ 1 \end{pmatrix} \end{aligned}

像素不是矩形
CCD行与列不垂直

\begin{pmatrix} x\\ y\\ 1 \end{pmatrix} = \begin{bmatrix} f_x & s & p_x \\ & f_y & p_y\\ & & 1 \end{bmatrix} \begin{bmatrix} 1 & & & 0\\ & 1 & & 0\\ & & 1 & 0 \end{bmatrix} \begin{pmatrix} X\\ Y\\ Z\\ 1 \end{pmatrix}

从相机坐标系的三维点 $X=(X,Y,Z)$ 到像素坐标的变换 $x=K[I|0]X_{cam}$ ， $K$ 为内参矩阵
若空间坐标系不与相机坐标系重合：
已知相机在世界坐标系中的位姿：

相机光心在世界坐标系的坐标为C
三坐标轴相对于世界坐标系的旋转为R（旋转矩阵）

$X_{cam}=R(X-C)=[R|-RC] \begin{pmatrix} X\\ 1 \end{pmatrix}$

$x=K[I|0]X_{cam}=K[R|-RC] \begin{pmatrix} X\\ 1 \end{pmatrix} =K[R|t]\tilde{X}$
相机投影矩阵： $P=K[R|t],t=-RC$ ， $K$ 是相机内参矩阵， $[R|t]$ 是相机外参矩阵
由此，我们得到了三维重建的基本参数
alt text

运动推断结构（Structure from Motion，SFM）

运动：相机的运动
结构：场景的三维点云
SFM：从相机运动获取场景的三维点云

如果运动已知，那么只需要图像匹配+三角化，然而相机运动和三维点云都未知。

极线几何
alt text
左视图点 $p$ 在右视图的对应点 $p'$ 一定位于 $l'$ 上
右视图点 $p'$ 在左视图的对应点 $p$ 一定位于 $l$ 上

极线约束(Epipolar Constraint)
alt text
$\hat{x_i}=K^{-1}x_j,\parallel \hat{x_j}\parallel=1$
$x_j$ 是成像面上的二维坐标 $\hat{x_j}$ 是相机坐标系的空间坐标，因为 $x=K[I|0]X_{cam}=KX_{cam}$ ，所以 $X_{cam}=K^{-1}x$ ，然后归一化到单位长度，作为一个方向向量，其方向由相机光心指向该点
$d_1\hat{x_1}=p_1=Rp_0+t=R(d_0\hat{x_0})+t$
$d_0,d_1$ 是深度，是未知的， $R,t$ 是右边坐标系相对左边坐标系的变换参数。理论上来说 $d_0\hat{x_0},d_1\hat{x_1}$ 会相交于一点 $p$ ，但由于相机分辨率以及噪声的存在，实际情况中可能不存在这样的交点，那么可以用公垂线的中点来代替 $p$ 点。
等式两边同时与 $t$ 做叉乘，即两边同时乘一个 $t$ 的叉积矩阵：
$d_1[t]_{\times}\hat{x_1}=d_0[t]_{\times}R\hat{x_0}$
等式两边同时与 $\hat{x_1}$ 做点积：
$d_1\hat{x_1}^T[t]_{\times}\hat{x_1}=d_0\hat{x_1}^T[t]_{\times}R\hat{x_0}$
等式左侧， $[t]_{\times}\hat{x_1}$ 的结果是一个向量，其方向垂直于 $\hat{x_1}$ ，再与 $\hat{x_1}$ 做点积，结果是0，所以：
$d_0\hat{x_1}^T[t]_{\times}R\hat{x_0}=0\implies\hat{x_1}^TE\hat{x_0}=0,E=[t]_{\times}R$
由此，我们得到了本质矩阵 $E$ ，由 $R,t$ 决定。如果已知本质矩阵 $E$ 的话，就能得到 $R,t$ 然后就能得到深度 $d_0,d_1$ 了。
从 $E$ 到 $R,t$
任意本质矩阵都可以通过SVD分解为如下形式： $E=[t]_{\times}R=Udiag(1,1,0)V^T$
，相应的 $[R,t]$ 存在4种可能： $[R,t]=[UWV^T,u_3]or[UWV^T,-u_3]or[UW^TV^T,u_3]or[UW^TV^T,-u_3]$ ， $u_3$ 是 $U$ 的第3个（最小奇异值）奇异向量；
$W=\begin{pmatrix} 0 & -1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 1 \end{pmatrix}$ 。
对于以上四种情况，只有(a)是符合实际情况的，所以对于四种可能的解，选取三维点在相机前方最多的解作为结果。
alt text

如何计算本质矩阵E

方法1：
基础矩阵 $\begin{cases} \hat{x_i}=K^{-1}x_j,\parallel \hat{x_j}\parallel=1\\ \hat{x_1}^TE\hat{x_0}=0 \end{cases} \implies x_1^T(K_1^{-1})^TEK_0^{-1}x_0=0\implies x_1^TFx_0=0$ 基础矩阵 $F=(K_1^{-1})^TEK_0^{-1}$
因为 $x_1^TFx_0=0$ ，所以一个对应的点对，就能得到一条方程，只要有足够的点对，就能得到基础矩阵 $F$ ，从而解得本质矩阵 $E=K_1^TFK_0$
方法2：
直接通过 $\hat{x_1}^TE\hat{x_0}=0$ 来估计本质矩阵，同样通过多个对应点对来估计。

二视图SFM：只能求出相机的相对位置和姿态；位移t只能得到方向
多视图SFM：二视图+融合的方法可能存在冲突，因为同一三维点可能在多个视图同时被观察到
为了使公式计算出的同一特征点对应的投影点尽量与ground truth一致，采用集束调整(Bundle Adjustment)的方法：

$\min\sum_i\sum_j(\tilde{x}_i^j-K[R_i|t_i]X^j)^2$
非线性最小二乘
使用二视图结果初始化

参考资料