从图片复原高度信息

本篇文章主要解释如何从图片信息中恢复物体的三维高度信息，并借以巩固自己的知识。作为第一篇CV相关的东西，稍微玩玩。

反射与摄像机模型

在单一平行光源且无环境光的条件下，模型认为某面点p处的亮度与光源照度$|\vec S_1|$、光源方向$\vec S_1$、面法线$\vec N(p)$、面反射度$\rho(p)$相关，满足公式

$B(p)=\rho(p) \vec N(p) \cdot \vec S_1$

认为摄像机响应为线性，即摄像机对亮度的实际成像与亮度线性相关。

$C(p)=kB(p)=k\rho(p) \vec N(p) \cdot \vec S_1$

而后，整理该式为

$\begin{cases} C(p)&=\vec V \cdot \vec g(p) \cr \vec V&=k\vec S_1 \cr \vec g(p)&=\rho(p)\vec N(p) \end{cases}$

其中，$V$是能够提前确定的环境因素。问题现在即变为已知$C(p)$和$\mathcal V$，求$\vec g(p)$。而后，设法线为单位法线，即有

$\begin{aligned} \vec N(p)&=\frac{\vec g(p)}{\left\|\vec g(p)\right\|} \cr \rho(p)&=\left\|\vec g(p)\right\| \end{aligned}$

现在，我们求得了面的法线。假设求得的单位法线为$(a(p),b(p),c(p))$。

求解$g(p)$

事实上，我们至少需要三个不同光源角度的同姿态照片，才可以求解。假设三次对p点测量结果分别为

$\mathcal C= \begin{pmatrix} C_1(p) \cr C_2(p) \cr C_3(p)\end{pmatrix}$

三次环境因素为

$\mathcal V= \begin{pmatrix} V_1^T \cr V_2^T \cr V_3^T \end{pmatrix}$

由此，有等式

$\mathcal C=\mathcal Vg(p)$

因$g(p)$需要确定的未知量有$3$个，所以至少需要三个光源数据，三光源必须在不同角度照射。

构建面

设我们将要复原的面为$(x,y,f(x,y))$。这个东西的单位法线可以求得，即

$\vec N(x,y)=\frac{1}{\sqrt{1+\frac{\partial f}{\partial x}^2+\frac{\partial f}{\partial y}^2}} \begin{pmatrix} \frac{\partial f}{\partial x}\cr \frac{\partial f}{\partial y}\cr 1 \end{pmatrix}$

于是，我们就找到了两个对应关系

$\begin{aligned} \frac{\partial f}{\partial x}&=\frac{a(p)}{c(p)} \cr \frac{\partial f}{\partial y}&=\frac{b(p)}{c(p)} \end{aligned}$

由此，我们就可以使用积分来复原一类$f(x,y)$了，具体来讲即

$f(x,y)=\oint_C \begin{pmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \end{pmatrix} \cdot dl+C_1$

其中$C$是一条从某个定点开始，到$(x,y)$结束的有向曲线。

细节

数据检查

由

$\frac{\partial^2 f}{\partial x\partial y}=\frac{\partial^2 f}{\partial y\partial x}$

可以得到一个检查数据的办法

$\frac{\partial (\frac{a(p)}{c(p)})}{\partial y}=\frac{\partial (\frac{b(p)}{c(p)})}{\partial x}$

因此，二者实际测量的差应该趋近于$0$，这保证了函数的可积分性。对于出现意外差值的数据，应该及时予以修正。

坐标系

$S_1$是在$x$点上指向光源的向量，其坐标在世界坐标系下。因此最终求出的法向量也是在世界坐标系下。这里隐藏了一个额外的问题，即如何确定平面姿态。然而这个问题看起来不好解决……所以应该尽可能以垂直角度拍摄待建物体，不要整太多花活。

实现

因为只是简单写写玩玩，也没有足够好的照片，所以就瞎糊点代码了。

img = Import["<path>"];
img = ImageResize[img, 100];

得到图片

随便假设光源在$(1,1,1)$附近。实际情况我也不知道该算它在哪，反正现在只是玩玩。代码求解出所有采样点的单位法向量（世界坐标系），并依据公式计算两个偏导数的数值，而后累加以重建$f(x,y)$。

normals =
 Map[Normalize@{x, y, z} /.
    Solve[{{.9, 1, 1} . {x, y, z} == #[[
         1]], {1, .9, 1} . {x, y, z} == #[[
         2]], {1, 1, .9} . {x, y, z} == #[[3]]}, {x, y, z}][[1]] &,
  ImageData[img], {2}];
deltas = Map[{#[[1]]/#[[3]], #[[2]]/#[[3]]} &, normals, {2}];

height = Accumulate /@ deltas[[;; , ;; , 2]];
vert = Accumulate@deltas[[;; , 1, 1]];
height = Table[height[[x]] + vert[[x]], {x, 1, Length[height]}];

输出结果，高处为黑色……emmmm……毕竟数据和图片都是瞎拍的，结果一看就不咋样。不过意思达到了，先这样吧。