通过投影统计计算灰度图像朝向

图像处理中会遇到需要计算物体朝向的情况，前文计算了二值图物体朝向，本文进一步放宽条件，计算灰度图像朝向。

问题描述

• 对于一个灰度图像，需要求解图像的朝向

解决思路

• 按照根据投影计算二值图朝向的思路，在二值图中默认密度为1，因此在计算中忽略了图像的像素值，对于灰度图相当于薄板密度不再为1

• 于是可以通过计算加权距离，求解灰度图像体朝向

• 我们可以将该带密度的薄板理解为平面上的二维概率分布，即将整个图像每个像素值除以总质量，即得到了总和为1的 $x,y$ 的联合分布

  当图中物体部分密度相同（非零），背景部分密度为 0 时，灰度图像朝向即坍缩为二值图的朝向计算问题

• 那么我们的目标是找到某个方向，使得该联合分布投影到这个方向得到的边缘分布方差最小，我们认为这个方向为与物体朝向垂直的方向。

计算方法

• 考虑单通道非零灰度图像$I$，假设其中的点有 $N$ 个
• 第$i$个点坐标为 $(x_i,y_i),0<i<N+1$，合并为列向量

$$\mathbf{\text{x}}=[x_1,x_2,\dots ,x_N{]}^T,\mathbf{\text{y}}=[y_1,y_2,\dots ,y_N{]}^T$$

• 合并为维度 $N\times 2$ 的矩阵 $\mathbf{\text{M}}=[\mathbf{\text{x}},\mathbf{\text{y}}]$

• 对应 $N$ 个点的密度 $[{\rho }_1,{\rho }_2,\dots ,{\rho }_N{]}^T$

• 总质量为 $S=\sum _{n=1}^N{\rho }_n$

• 将原始图像转义为二维分布，$N$ 个点的概率密度为 $\mathbf{\text{p}}=[p_1,p_2,\dots ,p_N{]}^T$，其中 $p_i={\rho }_i/S$

• 给定一条直线，直线过点 $(x_t,y_t)$，直线的法向量 $\mathbf{\text{v}}=[a,b{]}^T$，该法向量为单位向量，有：

$$a^2+b^2=1$$

• 则直线可表示为：

$$[x-x_t,y-y_t]\mathbf{\text{v}}=0$$

• 二值图上的点 $(x,y)$ 到直线的有向距离为： 点到$(x_t,y_t)$ 的向量到直线法向量 $\mathbf{\text{v}}$ 的投影:

$$d=[x-x_t,y-y_t]\mathbf{\text{v}}$$

• 所有点的有向距离组成向量：

$$\begin{array}{cc}\mathbf{\text{d}}& =[\mathbf{\text{M}}-[x_t,y_t]]\mathbf{\text{v}}\\ & =\mathbf{\text{M}}\mathbf{\text{v}}-[x_t,y_t]\mathbf{\text{v}}\end{array}$$

• 设 $[x_t,y_t]\mathbf{\text{v}}$ 为常数 $c$，有：

$$d=[x,y]\mathbf{\text{v}}-c，\mathbf{\text{d}}=\mathbf{\text{M}}\mathbf{\text{v}}-c$$

• 我们的目标是找到所有点的有向距离组成的样本集合的方差统计量最小，有向距离的均值为：

$$\begin{gathered} \bar{d}=\sum _{i=1}^N{p}_i{\mathbf{\text{m}}}_i\mathbf{\text{v}}-c \\ =[\bar{x},\bar{y}]\mathbf{\text{v}}-c \end{gathered}$$

其中 $[\bar{x},\bar{y}]$ 为 $x,y$ 的期望向量

• 灰度图像，有向距离的方差可表示为：

$$\begin{array}{ccc}Var& =& \frac{1}{N}\sum _{i=1}^N(d_i-\bar{d}{)}^2{p}_i\\ & =& \frac{1}{N}\sum _{i=1}^N([x_i,y_i]\mathbf{\text{v}}-c-[\bar{x},\bar{y}]\mathbf{\text{v}}+c{)}^2{p}_i\\ & =& \frac{1}{N}\sum _{i=1}^N([x_i-\bar{x},y_i-\bar{y}]\mathbf{\text{v}}{)}^2{p}_i\end{array}$$

• 可见二值图到直线的有向距离的方差与直线的偏置无关，仅与法向量有关，设 ${\mathbf{\text{M}}}^{\prime }=\mathbf{\text{M}}-[\bar{x},\bar{y}]$，有：

$$\begin{array}{ccc}Var& =& \frac{1}{N}\sum _{i=1}^N([x_i-\bar{x},y_i-\bar{y}]\mathbf{\text{v}}{)}^2{p}_i\\ & =& \frac{1}{N}({\sqrt{\mathbf{\text{p}}}}^T{\mathbf{\text{M}}}^{\prime }\mathbf{\text{v}}{)}^T({\sqrt{\mathbf{\text{p}}}}^T{\mathbf{\text{M}}}^{\prime }\mathbf{\text{v}})\\ & =& \frac{1}{N}{\mathbf{\text{v}}}^T{\mathbf{\text{M}}}^{\prime T}\sqrt{\mathbf{\text{p}}}{\sqrt{\mathbf{\text{p}}}}^T{\mathbf{\text{M}}}^{\prime }\mathbf{\text{v}}\end{array}$$

• 可以看做 $Var$ 是关于 $\mathbf{\text{v}}$ 的二次型，系数矩阵为 $\frac{1}{N}{\mathbf{\text{M}}}^{\prime T}\sqrt{\mathbf{\text{p}}}{\sqrt{\mathbf{\text{p}}}}^T{\mathbf{\text{M}}}^{\prime }$
• 注意：我们需要该方差最小，但不能直接对 $\mathbf{\text{v}}$ 求导，令其取值为零：

$$\begin{array}{c}\frac{\partial Var}{\partial \mathbf{v}}=\frac{2}{N}{\mathbf{\text{M}}}^{\prime T}\sqrt{\mathbf{\text{p}}}{\sqrt{\mathbf{\text{p}}}}^T{\mathbf{\text{M}}}^{\prime }\mathbf{v}\mathbf{=}\mathbf{0}\\ {\mathbf{\text{M}}}^{\prime T}\sqrt{\mathbf{\text{p}}}{\sqrt{\mathbf{\text{p}}}}^T{\mathbf{\text{M}}}^{\prime }\mathbf{v}\mathbf{=}\mathbf{0}\end{array}$$

• 因为这样做会得到 $\mathbf{v}$ 为零向量的平凡解，但是事实上 $\mathbf{v}$ 是有模长为 1 的约束的，因此我们需要解的是带约束的优化方程：

  $$\begin{array}{c}minimizeVar=\frac{1}{N}{\mathbf{\text{v}}}^T{\mathbf{\text{M}}}^{\prime T}\sqrt{\mathbf{\text{p}}}{\sqrt{\mathbf{\text{p}}}}^T{\mathbf{\text{M}}}^{\prime }\mathbf{\text{v}}\\ subjectto{\mathbf{v}}^{\mathbf{T}}\mathbf{v}\mathbf{-}\mathbf{1}\mathbf{=}\mathbf{0}\end{array}$$

解方程

拉格朗日

• 至此之后的方程解法和之前完全相同，只是矩阵数值根据密度发生了一些变化
• 之后解该方程需要在 $|\mathbf{v}\mathbf{|}\mathbf{=}\mathbf{1}$ 的约束下进行，因此可以引入拉格朗日乘数法
• 可以参考 拉格朗日解该方程 的流程

特征向量

• 我们用 $\mathbf{A}$ 表示上述系数矩阵 $\frac{1}{N}{\mathbf{\text{M}}}^{\prime T}\sqrt{\mathbf{\text{p}}}{\sqrt{\mathbf{\text{p}}}}^T{\mathbf{\text{M}}}^{\prime }$ (此处为 $2\times 2$ 的矩阵)
• 因此我们的优化目标变为: ${\mathbf{v}}^{\mathbf{T}}\mathbf{A}\mathbf{v}$
• 拉格朗日方程为：

$$L={\mathbf{v}}^{\mathbf{T}}\mathbf{A}\mathbf{v}\mathbf{+}\lambda \mathbf{(}{\mathbf{v}}^{\mathbf{T}}\mathbf{v}\mathbf{-}\mathbf{1}\mathbf{)}$$

• 对 $\mathbf{v}$ 求导，令其为 0：

$$\begin{array}{c}2A\mathbf{v}\mathbf{+}\mathbf{2}\lambda \mathbf{v}\mathbf{=}\mathbf{0}\\ A\mathbf{v}\mathbf{+}\lambda \mathbf{v}\mathbf{=}\mathbf{0}\\ (A+\lambda I)\mathbf{v}\mathbf{=}\mathbf{0}\end{array}$$

• 哈，这不就是在求 $\mathbf{A}$ 的特征向量吗？
• 求就完了，会得到两个特征向量，特征值小的对应的特征向量为方差最小的投影方向。

参考资料

• https://www.zywvvd.com/notes/study/image-processing/robot-vision/chapter-3/binary-image-moment/binary-image-moment
• https://www.zywvvd.com/notes/study/image-processing/robot-vision/chapter-3/obj-orientation/obj-orientation/

文章链接：
https://www.zywvvd.com/notes/study/image-processing/robot-vision/chapter-3/obj-orientation/gray-obj-orientation/