Gelsight Sensor 重建 height-map —— 使用 DST-II 算子的快速泊松求解

Gelsight Sensor 重建 height-map —— 使用 DST-II 算子的快速泊松求解

gelsight

gelsight传感器具有软弹性体的接触表面，可以精确的重建物体形状还可以从传感器的变形推断接触力和滑动.

使用 DST-II 算子的快速泊松求解 height

文章：GelSight: High-Resolution Robot Tactile Sensors for Estimating Geometry and Force

文章最后求解的方程其实是: $p=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}$, u 是我们要重建的高度图. p是以知的，是u两个二阶导数的和.

在离散域下, 某个像素点处, $\frac{{\partial }^{2}u(i,j)}{\partial x^2}$ 可以写成 $\frac{u(i+1,j)-u(i,j)-[u(i,j)-u(i-1,j)]}{△x^2}$，同样 $\frac{{\partial }^{2}u(i,j)}{\partial y^2}$ 可以写成 $\frac{u(i,j+1)-u(i,j)-[u(i,j)-u(i,j-1)]}{△y^2}$ 。

故而可以有: $p_{i,j}=\frac{u(i+1,j)-u(i,j)-[u(i,j)-u(i-1,j)]}{△x^2}+\frac{u(i,j+1)-u(i,j)-[u(i,j)-u(i,j-1)]}{△y^2}$

引入DST-II: 单个方向的 DST 这样定义 $Y_k=2{\sum }_{i=0}^{N-1}{x}_{i}sin(\frac{\pi (i+1/2)(k+1)}{N})$

将DST带入上一条式子可以有: $\tilde{u}$ 是 $u$ 的DST

$u(i,j)=4{\sum }_{k_i=0}^{N_x-1}{\sum }_{k_j=0}^{N_y-1}\tilde{u}(k_i,k_j)sin(\frac{\pi (i+1/2)(k_i+1)}{N_x})sin(\frac{\pi (j+1/2)(k_j+1)}{N_y})$

$u(i-1,j)=4{\sum }_{k_i=0}^{N_x-1}{\sum }_{k_j=0}^{N_y-1}\tilde{u}(k_i,k_j)sin(\frac{\pi (i-1+1/2)(k_i+1)}{N_x})sin(\frac{\pi (j+1/2)(k_j+1)}{N_y})$

$u(i+1,j)=4{\sum }_{k_i=0}^{N_x-1}{\sum }_{k_j=0}^{N_y-1}\tilde{u}(k_i,k_j)sin(\frac{\pi (i+1+1/2)(k_i+1)}{N_x})sin(\frac{\pi (j+1/2)(k_j+1)}{N_y})$

将上面三个结果代入 $\frac{u(i+1,j)-u(i,j)-[u(i,j)-u(i-1,j)]}{△x^2}$, 得到:

$\frac{4}{△x^2}{\sum }_{k_i=0}^{N_x-1}{\sum }_{k_j=0}^{N_y-1}\tilde{u}(k_i,k_j){\lambda }_{x_i}sin(\frac{\pi (i+1/2)(k_i+1)}{N_x})sin(\frac{\pi (j+1/2)(k_j+1)}{N_y})$    (1)

其中 ${\lambda }_{x_i}=2\ast [\cos (\frac{\pi (k_i+1)}{N_x})-1]$

注：$$\begin{gathered} sin(\frac{\pi (i+1+1/2)(k_i+1)}{N_x})+sin(\frac{\pi (i-1+1/2)(k_i+1)}{N_x})-2\ast sin(\frac{\pi (i+1/2)(k_i+1)}{N_x}) \\ =\sin (\frac{\pi (k_i+1)}{N_x})\cos (\frac{\pi (i+1/2)(k_i+1)}{N_x})+\cos (\frac{\pi (k_i+1)}{N_x})\sin (\frac{\pi (i+1/2)(k_i+1)}{N_x}) \\ -\sin (\frac{\pi (k_i+1)}{N_x})\cos (\frac{\pi (i+1/2)(k_i+1)}{N_x})+\cos (\frac{\pi (k_i+1)}{N_x})\sin (\frac{\pi (i+1/2)(k_i+1)}{N_x}) \\ -2\ast sin(\frac{\pi (i+1/2)(k_i+1)}{N_x}) \\ =2\ast [\cos (\frac{\pi (k_i+1)}{N_x})-1]\ast \sin (\frac{\pi (i+1/2)(k_i+1)}{N_x}) \\ ={\lambda }_{x_i}\ast \sin (\frac{\pi (i+1/2)(k_i+1)}{N_x}) \end{gathered}$$

同样的可以得到 $\frac{u(i,j+1)-u(i,j)-[u(i,j)-u(i,j-1)]}{△y^2}$ 对应的 $\frac{4}{△y^2}{\sum }_{k_i=0}^{N_x-1}{\sum }_{k_j=0}^{N_y-1}\tilde{u}(k_i,k_j){\lambda }_{y_j}sin(\frac{\pi (i+1/2)(k_i+1)}{N_x})sin(\frac{\pi (j+1/2)(k_j+1)}{N_y})$    (2)

其中 ${\lambda }_{y_j}=2\ast [\cos (\frac{\pi (k_j+1)}{N_y})-1]$

对于 $p_{i,j}$ ，也可以使用DST:
$p(i,j)=4{\sum }_{k_i=0}^{N_x-1}{\sum }_{k_j=0}^{N_y-1}\tilde{p}(k_i,k_j)sin(\frac{\pi (i+1/2)(k_i+1)}{N_x})sin(\frac{\pi (j+1/2)(k_j+1)}{N_y})$    (3)

由 (3)=(1)+(2) 可以有: $\tilde{p}(k_i,k_j)=[\frac{{\lambda }_{x_i}}{△x^2}+\frac{{\lambda }_{y_j}}{△y^2}]\tilde{u}(k_i,k_j)$

所以 $\tilde{u}(k_i,k_j)=\tilde{p}(k_i,k_j)/[\frac{{\lambda }_{x_i}}{△x^2}+\frac{{\lambda }_{y_j}}{△y^2}]$, 得到 $\tilde{u}(k_i,k_j)$ 后在反求 IDST 就可以得到 $u(k_i,k_j)$

实际做法

文章求解高度时构建了像素强度I到放射率R的映射关系，称为lookup table，用于存放图像的两个一阶梯度 (还有一些其他的参数), 一阶梯度相减得到二阶梯度

	% Laplacian
	gyy(j+1,k) = gy(j+1,k) - gy(j,k);
	gxx(j,k+1) = gx(j,k+1) - gx(j,k);

得到p

    p = gxx + gyy; 

p的DST-II

    tt = dst(p2); 
	p2sin = dst(tt')';

$\tilde{p}(k_i,k_j)/[\frac{{\lambda }_{x_i}}{△x^2}+\frac{{\lambda }_{y_j}}{△y^2}]$

    denom = (2*cos(pi*x/(xdim-1))-2) + (2*cos(pi*y/(ydim-1)) - 2);

	p3 = p2sin./denom; 

反求 IDST 就可以得到 $u(k_i,k_j)$

    hh = idst(p3); img_tt = idst(hh')'; 

参考：
Fast-Poisson-Equation-Solver-using-DCT