Jinyu Li

Photo

2023-10-06T01:21:09+00:00

旅途中偶得～

class Life : public SLAM<Life>

2022-10-06T12:24:36+00:00

The long dormancy is over, now it’s the time for new adventures.

The excitement of learning separates youth from old age. As long as you're learning you're not old.

—Rosalyn S. Yalow

写出了一个有趣的 bug ……

2019-05-08T18:53:29+00:00

背景建筑图片来自 Wikimedia 。

微软从来不会数数

2018-10-22T10:00:00+00:00

如题，亮点自寻……

旋转变量求导

2018-07-10T11:35:00+00:00

在涉及空间三维运动的问题中经常会遇到含有旋转变量的能量，最优化这样的能量时，便需要对旋转变量进行求导。在当前旋转的切空间上，运用李代数的方法进行局部求导已经是一种标准的做法，但其中推导比较复杂，不熟悉的话容易出错。为此，这篇文章整理了一些常见的包含旋转变量的能量形式以及它们的导数。

基本方法

我们回避过于形式化的完整推导，采用一种仿照对偶数¹的简易推导方法。如果对于一个函数 $f$ 在点 $x$ 处有 $f(x)+\delta f = f(x+\delta x)$ 其中，$\delta x$ 和 $\delta f$ 看作是微元或者对偶量，那么我们就可以从中提取 $\delta f/\delta x$ 作为 $f(x)$ 的导数。从具体步骤上来说，首先在需要求导的点 $x$ 上添加微元 $\delta x$ 后应用函数 $f$ ，然后将其展开为 $f(x)$ 加上微元余项的形式，最后就可以提取导数了。

旋转小量

三维旋转只有三个自由度，然而它高度非线性，因而我们采用李代数的方法，在切空间中进行求导。这样一来，旋转小量便对应了三维的向量。旋转切空间的李代数到旋转的李群之间由指数映射和对数映射联系起来，对一个旋转 $R$ 添加小量微扰 $\delta R$ 可以表示为 $R\oplus \delta R = R\exp \delta R$ 。相对的，两个相近的旋转 $R_1$ 和 $R_2$ 之间的差可以表示为 $R_2 \ominus R_1 = \log(R_1^\top R_2)$ 。

此外，我们有一些有用的关系：

$\exp(\delta R) = I+[\delta R]_\times$
$\exp(\delta R)U = U\exp(U^\top\delta R)$
$\log(U\exp(\delta R)) = \log U + J_r^{-1}(\log U)\delta R$

注意由于 $\delta R$ 是微元，上面三个关系的等号可以认为成立，否则第一和第三个关系只在 $\delta R$ 比较小时近似成立。

常见基本形式的导数

1. 三维点的旋转

空间三维点的旋转对应于 $f_1(R) = Rx$ 或 $f_2(R) = R^\top x$ ，其中 $x$ 是三维常向量。以 $f_1$ 为例：

\[\begin{align} f_1(R\oplus\delta R) - f_1(R) &= (R\exp\delta R)x-Rx \nonumber \\ &= R(I+[\delta R]_\times)x - Rx \nonumber \\ &= R[\delta R]_\times x \nonumber \\ &= -R[x]_\times \delta R \end{align}\]

因此，对应的导数为 $-R[x]_\times$ 。

下表为上述两个形式对应的导数：

$i$	$f_i(R)$	$[f_i(R\oplus\delta R)-f_i(R)]/\delta R$
$1$	$Rx$	$-R[x]_\times$
$2$	$R^\top x$	$[R^\top x]_\times$

2. 旋转的复合

旋转的复合对应于 $f_3(R) = UR$ 、$f_4(R) = UR^\top$ 、$f_5(R) = RU$ 和 $f_6(R) = R^\top U$ 四种形式中的一种，其中 $U$ 为一个常量旋转。由于这四个形式最终结果是旋转，根据参数化方式的不同，结果的维度也不一样，我们转而在结果的切空间内求导，此时导数对应了 $3\times3$ 矩阵。结果的导数形式中，所有的旋转也相应采用其 $3\times3$ 矩阵形式表达，以 $f_5$ 为例：

\[\begin{align} f_5(R\oplus \delta R) \ominus f_5(R) &= (R\exp (\delta R) U) \ominus f_5(R) \nonumber \\ &= (RU \exp (U^\top \delta R)) \ominus (RU) \nonumber \\ &= \log((RU)^\top RU\exp(U^\top\delta R)) \nonumber \\ &= U^\top\delta R \end{align}\]

因此，对应的切空间内的导数为 $U^\top$ 。下表列出了上述四个形式对应的导数：

$i$	$f_i(R)$	$[f_i(R\oplus \delta R) \ominus f_i(R)]/\delta R$
$3$	$UR$	$I_{3\times 3}$
$4$	$UR^\top$	$-R$
$5$	$RU$	$U^\top$
$6$	$R^\top U$	$-U^\top R$

3. 旋转的差值

旋转的差值对应于 $f_7(R) = \log(UR)$ 、 $f_8(R) = \log(UR^\top)$ 、$f_9(R) = \log(RU)$ 和 $f_{10}=\log(R^\top U)$ 四种形式中的一种。以 $f_9$ 为例：

\[\begin{align} f_9(R\oplus\delta R)-f_9(R) &= \log(R\exp(\delta R) U) - f_9(R) \nonumber \\ &= \log(RU\exp(U^\top\delta R)) - f_9(R) \nonumber \\ &= \log(RU)+J_r^{-1}(\log(RU))U^\top \delta R-f_9(R) \nonumber \\ &= J_r^{-1}(\log(RU))U^\top\delta R = J_r^{-1}(f_9(R))U^\top\delta R \end{align}\]

因此，对应的导数为 $J_r^{-1}(f_9(R))U^\top$ 。

下表列出了这四个形式对应的导数：

$i$	$f_i(R)$	$[f_i(R\oplus \delta R)-f_i(R)]/\delta R$
$7$	$\log(UR)$	$J_r^{-1}(f_7(R))$
$8$	$\log(UR^\top)$	$-J_r^{-1}(f_8(R))R$
$9$	$\log(RU)$	$J_r^{-1}(f_9(R))U^\top$
$10$	$\log(R^\top U)$	$-J_r^{-1}(f_{10}(R))U^\top R$

在面对复合函数时，只要正确使用链式法则结合上面的导数，便可以得到想要的结果了。

Dual Number on Wikipedia ↩

OpenCV 中进行 YUV420 到 RGB 的转换

2018-03-28T10:46:00+00:00

最近的一个项目中，需要从远程接收 H.264 编码的视频流，然后在本地进行解码后使用图像。

本地从 H.264 码流解码出的图像采用了 YCbCr(I420) 格式，对于一帧 $W\times H$ 的图像，它将像素色彩的三个分量分别储存在三个平面上，也就是三个缓冲区，首尾相连。第一个平面对应 Y 色度平面，储存了亮度信息，它的大小也是 $W\times H$。第二个和第三个平面对应 Cb 和 Cr 色度平面，储存了蓝色和红色的色差信息，由于利用了人眼视觉上对色彩的空间分布相对不敏感的特点，这两个色度平面对原始图像的两个方向上各进行了$\frac12$的降采样，也就是说每$2\times 2$的像素共用一组 $(C_b,C_r)$ 色差信息。由于这个原因，这两个平面的大小都是 $\frac W2 \times \frac H2$。

理论虽然是这样，实际要用 OpenCV 来做就比较麻烦了（并不）。

先说说麻烦在那里：由于三个色度平面大小不一，在如何使用 cv::Mat 表示上就产生了许多的猜测。网上一番搜索，基本都是介绍原理，实现需要一个像素一个像素手工转换。对于我这种懒人，我是非常希望可以用一句 cv::cvtColor 搞定的。这样做有很多的好处，比如代码更精简，比如工作量小，比如可以更快~~，比如有问题了可以甩锅给 OpenCV~~…… 可稍微看了一下文档就会意识到， OpenCV 的文档根本就是垃圾，从文档是不可能知道该怎么用的！

那么真的要自己去实现一个么？不需要！ cv::cvtColor 究竟支不支持这样的转换呢？支持！

所以，怎么做呢？

有的时候，代码就是最好的文档（不过 OpenCV 又一次证明我错了）。经过浏览 OpenCV 中 cv::cvtColor 的代码实现，我注意到了这样的段落（项目中使用的是 OpenCV 2.4.13）：

// :::
case CV_YUV2BGR_YV12: case CV_YUV2RGB_YV12: case CV_YUV2BGRA_YV12: case CV_YUV2RGBA_YV12:
case CV_YUV2BGR_IYUV: case CV_YUV2RGB_IYUV: case CV_YUV2BGRA_IYUV: case CV_YUV2RGBA_IYUV:
    {
        //http://www.fourcc.org/yuv.php#YV12 == yuv420p -> It comprises an NxM Y plane followed by (N/2)x(M/2) V and U planes.
        //http://www.fourcc.org/yuv.php#IYUV == I420 -> It comprises an NxN Y plane followed by (N/2)x(N/2) U and V planes
// :::

从这里可以看到 CV_YUV2{BGR|RGB}[A]_{YV12|IYUV} 系列枚举值对应了我们需要的色度转换。通过进一步阅读后面的处理（朋友，我替你读过了，所以不要读了，太浪费生命），可以知道正确的方法是建立一个 cv::Mat ，直接按照标准 YCbCr(I420) 的平面大小分配空间并填充数据，然后调用 cv::cvtColor 。整个过程的代码如下：

cv::Mat YUV420_to_BGR888(
	int width, int height,
	const uchar * Y, const uchar * Cb, const uchar * Cr,
    int strideY, int strideCbCr
) {
	int uvwidth = width / 2;
	int uvheight = height / 2;
	int size = width * height;
	int uvsize = uvwidth * uvheight;
	
	// whole data size = Y + U + V = W*H + (W/2)*(H/2)*2 = W*(H+H/2)
	cv::Mat YCbCrData(cv::Size(width, height + uvheight), CV_8UC1);
	
	// get the pointer to the beginning of each plane.
	uchar* pY = YCbCrData.data;
	uchar* pCb = pY + size;
	uchar* pCr = pCb + uvsize;
	
	// copy Y channel for each line
	for (int i = 0; i < height; ++i) {
	    memcpy(pY + i*width, Y + i*strideY, width);
	}
	
	// copy Cb and Cr channel for each line
	for (int i = 0; i < uvheight; ++i) {
	    memcpy(pCb + i*uvwidth, Cb + i*strideCbCr, uvwidth);
	    memcpy(pCr + i*uvwidth, Cr + i*strideCbCr, uvwidth);
	}
	
	cv::Mat BGRData(cv::Size(width, height), CV_8UC3);
	
	cv::cvtColor(YCbCrData, BGRData, CV_YUV2BGR_IYUV);
	
	return BGRData;
}

在上面的代码里，我假定了 width 和 height 都是偶数，奇数的情形还请自行处理。此外，可以看到我这里对 YCbCr 和 YUV 不加区分的使用了，这两者根据场合还是存在不同的。fourcc 的网站上提供了比较详细的信息，可以参考 YUV pixel formats 和 YUV to RGB Conversion 。

在 vector 里使用 bool

2018-01-25T05:39:00+00:00

Q: 我想在 std::vector 里使用 bool 可以么？
A: 当然是不可以！

Q: 为什么不可以？
A: 其实也可以……

Q: 一会儿可以一会儿不可以，别逗我行不行？
A: 你试试 &v[0] ……

Q: 为啥别的类型都行，但 std::vector<bool> 就不让 &v[0] ？
A: 因为 std::vector 对 bool 进行了特化，把若干 bool 压缩存储在了整数类型的每个二进制位上，因此无法取地址。

Q: 好吧那怎么办……
A: 下面这段代码送给你，good luck and have fun！

struct boolean {
    boolean() = default;
    boolean(bool value) : value(value) {}
    operator bool&() { return value; }
    operator const bool&() const { return value; }
    bool* operator&() { return &value; }
    const bool* operator&() const { return &value; }
private:
    bool value;
};

最小二乘问题（十五）

2017-12-05T15:22:36+00:00

对于线性最小二乘问题 $\min\|Ax-b\|^2$ ，我们知道它的最优解同时是标准方程 $A^TAx = A^Tb$ 的解。

如果设 $M = A^TA$ ，在 $A$ 列满秩时，$M$ 是对称正定矩阵。我们知道，存在唯一的下三角矩阵 $L$ ，在对角元大于零的前提下满足 $LL^T = M$ 。这被称为矩阵 $M$ 的 Cholesky 分解。

与此同时，如果我们对 $A$ 进行 QR 分解，同样存在唯一的 $QR=A$ 满足 $R$ 的对角元大于零。将后者代入 $M=A^TA$ ，便可发现 $R^TR = R^TQ^TQR = A^TA = M = LL^T$ 。这意味着：矩阵 $A$ 的 QR 分解对应着矩阵 $M$ 的 Cholesky 分解。

回到最开始的线性最小二乘问题，若 $QR = A$ ，对原始问题的余项进行正交投影变换得到如下的新问题：

\[\min\|Q^T(Ax-b)\|^2 = \min\|Rx-Q^Tb\|^2\]

它的标准方程是 $R^TRx=R^TQ^Tb$，也就恰好是 $A^TAx=A^Tb$ 。这也就是说这个新问题具有与原问题相同的最优解。

前面的一系列分析基于了 $A$ 列满秩的假设，在 $A$ 的行数小于列数时，这是不成立的。但抛开矩阵分解结果唯一性，上面的 QR 分解与 Cholesky 分解的类推关系依旧成立。这意味着我们依旧可以用这种正交投影的方式来得到与原始问题具有相同解的新问题。

The Busy Indicator

2017-11-29T14:01:58+00:00

当手机从一天不充没问题变成每天两充还不够的时候，不是手机电池不行了，是接近年底了真的很忙啊……

“阈值”优化

2017-11-20T13:53:39+00:00

有时我们想求一系列数值 $f_1, f_2, \dots, f_n$ 中绝对值大于 $\lambda$ 的部分，如何非常有逼格地做这件事儿呢？

有趣的是，下面的优化解答了这一问题，其中 $f=(f_1, f_2, \dots, f_n)^T$ 是这些数值构成的列向量：

\[\min \frac12\|x-f\|_2^2+\lambda\|x\|_1\]

求解令上面问题最优的 $x$ ，则 $x$ 中非零项所在的维度便对应了想要的 $f_i$ 所在的维度。更具体地说，$x_i = \mathrm{sgn}(f_i)\cdot\max(|f_i|-\lambda, 0)$ 。这个操作等于将 $f$ 的所有元素向零方向收缩 $\lambda$，因此得名 shrinkage 。

从这个看似简单的优化，可以扩展出很多复杂的数学工具，包括著名的 LASSO 和 SVT 算法。