原论文：https://people.csail.mit.edu/kaiming/sig13/index.html

摘要

拍摄全景照片的时候，因为人手在抖，所以可能拍出来的图片边界就不是一个矩形框。在本文中，作者提出了一种可以感知图片内容的弯曲（warp）算法，来将这种形状非矩形的图片转化为矩形图片。分为两步，第一步是没有mesh的，初步把一张图片弯曲成矩形。第二步就是在矩形中建立mesh，然后在其中优化图片，使得图片中的各种形状保持原样（例如直线保持直线的样子）。

简介

同上介绍的，关于拍摄全景照片，因为手抖，无法拉出一个矩形图片。

一个简单的解决办法是：剪裁。从这个图片中剪裁出一个矩形区域。缺点是会损失很多信息。还有一种方法就是在原来的图片边缘的外面进行图片合成，缺陷是，对于复杂的语义信息无法处理。

作者这里提出了一种只会引入“可接受”的、“不容易注意到”的扭曲，并且保护图片原本内容的方法。

前人工作

Image alignment and stitching

全景照片实际上是将多张照片缝合（stitching）在一起的。通常用某种图片上的特征，来将原照片全部投影到一个相同的坐标系下。

Projections

众所周知，3D场景投影到2D会不可避免地导致一些折叠的扭曲。例如透视投影，虽然可以保持直线的形状，但是其他物体的形状则会被拉伸。圆柱投影和球面投影可以保持物体的形状，但会弯曲直线。

也有研究者提出了一种在自适应流形（adaptive manifold）上进行投影和缝合的方法。不过它是适应相机的移动，而不是适应图片的内容。

之前也有大量的工作，考虑将这些折叠后的扭曲，通过某种方式，让它们出现在不容易注意到的位置。虽然视觉上可以接受，但是无法移除固有的扭曲。

Image retargeting and warping

Image retargeting是一种基于图片内容的缩放图片的方式。但是它们都假设输入的图片具有矩形形状，而且其中一部分必须还要定义一个grid mesh。

Image completion

即前面提到的，在图片的边缘增加内容，使其变成一个矩形图片。一般的方法是把图片里面的内容，复制到需要合成的地方。但是对于需要生成语义内容来说，表现不好。

算法

算法分为两部分，一部分是改进的Seam Carving，是一个Local Warp的步骤，其将图片初步弯曲成矩形。然后另一部分（global warp）在前一步的基础上，建立一个grid mesh，在mesh上进行优化，来保证原图的物体形状和直线形状不变。

Local Warp

原版的Seam Carving算法是在图片中插入水平或者垂直的接缝（seam），来穿过整张图片，然后再从接缝中水平或者垂直地扩展一个像素。作者这里想到，接缝不穿过整张图片，这样就可以改变图片边缘的形状，从而变成矩形。

原版的算法可见https://faculty.runi.ac.il/arik/scweb/imret/index.html

如上图，其中蓝色的部分是有效的图片像素，白色的部分是需要我们填充的像素。

首先定义boundary segment是待填充像素中的，在图片顶部/底部/左侧/右侧的一行/一列像素。例如上图。我们找到目前最长的一个boundary segment，然后从它扩展得到一个sub-image。

这里左侧右侧的boundary segment，扩展出来的sub-image的高度保持，而宽度则等于整张图。同理，顶部和底部的boundary segment，扩展出来的sub-image的宽度保持，高度等于整张图。

然后我们在这个sub-image上利用seam carving的方法，去进行图像扩展。方法很简单，首先根据图片颜色梯度信息，计算每一个像素的能量。然后使用动态规划，求解最短路径。

以Figure 3(ii)为例，计算得到每个像素的能量。因为我们需要向右边填入像素，所以我们需要一个竖直方向的seam。于是我们用动态规划的方式，求底部像素到顶部像素的最短路（8邻居路径），然后所有底部像素中距离顶部像素最近的像素，及其到顶部的路径，成为我们的seam。

注意这里使用8邻居是源于原算法的研究，另外，我们需要给待填入的像素的能量赋值为$10^8$或者更大的数，来防止seam穿过这些像素。

填充像素的时候，我们把seam及其右边的像素都往右移动一位。原本seam的位置则求取左右两个像素的平均。

原版算法中提到，为了防止反复在相同的位置插入seam、导致如下的情况

之前的作者建议同时采集多根seam，然后依次填入像素。

但是这在我们这篇文章中是行不通的。我们的Local Warp要求我们每次填充完一根seam之后，都要重新计算最长的bound segment，我们不能反复对同一个sub-image填入像素。

我这里的实现采用了一个比较trick的方法，算出seam后，对原本的seam的位置的能量设置为$10^5$这样的数，比一般的大，比非法像素的小。之后让像素和能量都右移。现在原本的seam位置的能量和其右边一个像素的能量都为$10^5$了，这样就有更少的可能会在这里再次插入seam。

总而言之，重复以上过程，直到填充完所有像素。例如我实现的下例

其原图为

可以看到，虽然填充是填充了，但是歪歪扭扭。于是还需要下面的步骤。

Global Warp

Mesh Placement

首先，我们在Local Warp之后的图片的长宽上，将其分成$(40-1)*(10-9)$个网格，其中$40*10$是网格的顶点数量。网格是正交的。

在Local Warp中，我们对像素进行了移动，我们需要记录下每个像素移动的偏移量。然后，我们使用网格顶点所在像素的偏移量，将各个网格顶点移动回去，就完成了mesh placement这一步。之后local warp的图像就可以抛弃了。移动后的结果例如：

Energy Function

这里设计能量的目的有三个，对于目前的网格与最终的网格结果之间：

1. 保持网格的形状尽量不变，亦即Shape Preservation
2. 保持网格内的线段的笔直，并且与其平行的其他线段，在处理后依然平行，亦即Line Preservation
3. 网格上下左右贴合目标矩形图像框，亦即Boundary Constraints

Shape Preservation

最终输出的mesh的顶点（vertex）集记作$V$，其中的每一个元素$v_i=(x_i, y_i)$，可以用坐标表示出来。我们进行完Mesh Placement后得到的初始mesh记作$\hat V$

保型的能量定义为

$$E_S(V) = \dfrac{1}{N}\sum_q||(A_q(A_q^TA_q)^{-1}A_q^T-I)V_q||^2$$

这里$N$是quad的个数。其中

$$A_q = \begin{bmatrix} \hat x_0 & -\hat y_0 & 1 & 0\\ \hat y_0 & \hat x_0 & 0 & 1\\ \vdots & \vdots & \vdots & \vdots\\ \hat x_3 & -\hat y_3 & 1 & 0\\ \hat y_3 & \hat x_3 & 0 & 1 \end{bmatrix}$$

$$V_q = \begin{bmatrix} x_0 \\ y_0 \\ \vdots \\ x_3 \\ y_3 \end{bmatrix}$$

这里也就是说第$q$个quad的四个顶点的坐标。其中$A_q$中的是初始mesh的坐标，$V_q$是目标mesh的坐标。

Line Preservation

我们目前计算过的东西是不包含“线段”这样的信息的。我们首先要做的事是，使用任何一种线段检测算法，来算出图片中的线段。

根据作者的说法，我们需要将这些线段切分开来，放到mesh的各个quad里。然后将所有线段的方向角（范围$[-\dfrac{\pi}{2},\dfrac{\pi}{2})$），离散化为$M=50$个角度，每个离散化角度称为一个bin。为了保持直线性和平行性，设计的能量需要让每个bin里的线段拥有同样的角度。

给定一个线段，其方向向量$e$可以由线段的两个端点相减得到。如果我们将线段的两个端点用其所在的quad的$V_q$进行双线性插值表示，那么$e$就可以表示为$V_q$的线性函数。

同样的$e$是我们的目标，而$\hat e$则是最初的对应线段的方向向量。对于某一条线段，给定旋转操作的目标角度$\theta_m$，它的能量为

$$||sR\hat e-e||^2$$

其中

$$R = \begin{bmatrix} \cos\theta_m & -\sin\theta_m\\ \sin\theta_m & \cos\theta_m \end{bmatrix}$$

是旋转矩阵，而

$$s=(\hat e^T\hat e)^{-1}\hat e^TR^Te$$

是放缩因子。能量可以重新写为

$$||Ce||^2$$

其中

$$C = R\hat e(\hat e^T\hat e)^{-1}\hat e^TR^T-I$$

因为$\hat e$和$e$都可以写作$\hat V_q$和$V_q$的线性函数，所以能量就可以写成二次函数。

总的能量为

$$E_L(V, \{\theta_m\}) = \dfrac{1}{N_L}\sum_j||C_j(\theta_{m(j)})e_{q(j)}||^2$$

这里$N_L$为线段总数。对于第$j$条线段，$\theta_{m(j)}$代表其所属的bin中的目标旋转角度，$q(j)$代表其所属的quad。

Boundary Constraints

直白点说，就是把mesh的上下左右拉到图片框边缘上。

$$E_B(V) = \sum_{v_i\in L} x_i^2 + \sum_{v_i\in R} (x_i-w)^2 + \sum_{v_i\in T} y_i^2 + \sum_{v_i\in B} (y_i-h)^2$$

$L/R/T/B$分别是左右上下边界的mesh顶点，$w, h$则是图片框的宽高。

最终能量

将前面的几个加和，添加权重

$$E(V, \{\theta_m\}) = E_S(V)+\lambda_LE_L(V, \{\theta_m\}) + \lambda_B E_B(V)$$

显然，我们需要图片至少能够贴合图片框，我们设置$\lambda_B=10^8$来确保这一点。然后，关于$\lambda_L$的选择，作者说选取$>10$比较好，原文采取了$\lambda_L=100$。

优化算法

作者指出，可以分别优化$V$和$\{\theta_m\}$。每各优化一次算迭代一次，总共迭代10次。

固定$\{\theta_m\}$优化$V$

此时能量函数是$V$的二次函数，可以求解线性方程来计算。总共就几百个顶点，所以总体上来说，这部分的开销很小。

这里的求解是通过最小二乘法来算的，具体见我自己的实现部分。

固定$V$优化$\{\theta_m\}$

因为另外两项和$\theta$无关，所以我们只用优化$E_L$这一部分。

作者指出，可以通过牛顿迭代等方法去优化。但是，几何意义上，$E_L$的目标是找到一个$\theta_m$，使得所有在第$m$个bin中的线段$e$和其初始线段$\hat e$的夹角，都近似等于$\theta_m$。所以，我们可以简单的计算所有这个夹角，然后加和，求平均，作为新的$\theta_m$。

防止拉伸和后处理

如果我们的目标图相框的长宽比设置的不是很好，那么我们的结果图就会被拉伸。为了减少这种拉伸问题，我们在global warping之后更新目标框。

对于每个quad，我们计算它的$x$方向的伸缩因子为

$$s_x=(x_{max}-x_{min})/(\hat x_{max}-\hat x_{min})$$

所有quad的平均伸缩因子就为$\bar s_x$。同理可以计算出$\bar s_y$。之后我们的目标框就设置为$(w/\bar s_x, h/\bar s_y)$。

之后我们再在这个新的目标框上跑一次global warping

获得最终的mesh之后，我们就可以通过一些方法来获得最终的图像了。以图形学的思想为例，最终的mesh构成了我们要绘制的三角面。三角面空间坐标为最终mesh的坐标，uv坐标则是初始mesh的坐标。texture则是原图。

作者指出，这样的实现有一个小问题，因为mesh不会只包住有效像素，例如

这里的quad有一部分就包含了外部的无效像素。进行纹理采样时，作者简单地将无效像素替换为最近的有效像素。

实现细节

由于mesh是很“光滑”的，作者指出我们可以在缩小的图片上进行mesh的计算。首先将图片缩小到一个固定的大小，然后再在上面进行local/global warping。然后我们将mesh的尺度放大到原来的图片大小，再用这个mesh去进行纹理采样。

本文局限

如果图片缺失的部分太多，那么处理结果不好

本文的算法可能会导致边缘凹进去

作者指出，可以手动添加一个透明区域，当作已知像素。

将其当作已知像素来warp，之后在本文的算法处理完之后，再在这个透明区域中使用image completion等方法。

另外，本文无法处理未识别到的线段。这可以通过用户交互来缓解。

最后，本文无法应对线段数量过多的情况，这是warping方式的共同挑战。

我的实现

https://github.com/kegalas/Rectangling_Panoramic_Images

一些细节

首先是我之前也提到过的，为了防止反复在同一个地方插入seam，我们添加了额外的一个像素类别SEAM_PIXEL，让他的能量为$10^5$。比一般的像素大，比无效像素小。

在Global Warping主要要说明的部分则是关于如何优化能量。从$E_S$开始

$$E_S(V) = \dfrac{1}{N}\sum_q||(A_q(A_q^TA_q)^{-1}A_q^T-I)V_q||^2$$

我们需要使其最小。众所周知，使向量的范数平方最小，只用使其每一维都最小即可。或者说，使向量尽可能成为零向量。

这里的$(A_q(A_q^TA_q)^{-1}A_q^T-I)$是一个$8\times 8$矩阵（记作$A$好了），而$V_q$是一个8维向量。我们可以将所有的$V_q$收尾接在一起。如果我们的mesh是$40\times 10$的，那么新的拼接向量$V_Q$的维度就是$(40-1)\times(10-1)\times 8$

$$V_Q = \begin{bmatrix} V_1 \\ V_2 \\ \vdots \\ V_{N-1} \\ V_{N} \end{bmatrix}$$

同理，我们也可以把$A$拼接。如下

$$A_Q = \begin{bmatrix} A_1 & & & \\ & A_2 & & \\ & & \ddots & \\ & & & A_N \end{bmatrix}$$

我们的目标是使得

$$E_S(V) = \dfrac{1}{N}A_QV_Q = \overrightarrow{0}$$

这里$A_Q$是用$\hat V$算的，而$V_Q$是我们最终mesh的坐标。也就是说我们需要解方程算出来$V_Q$。这其实就是最小二乘法。我们使用Eigen的QR分解对$A_Q$分解，然后再求解$V_Q$。这里建议使用稀疏矩阵。

另外，每一个mesh的顶点只能有一个实例，不能单纯地把所有quad的四个顶点放进去同一列。每个顶点只能出现在$V_Q$里一次。我们直接存储$V_Q$为$40\times 10\times 2$维向量。即每个顶点都只存一份（包含$x,y$坐标）。然后我们去修改$A_Q$的顺序就可以在向量乘法的时候一一对应。具体见我写的代码。调整后写为

$$E_S(V) = \dfrac{1}{N}AV = \overrightarrow{0}$$

接下来比较简单的是$E_B$

$$E_B(V) = \sum_{v_i\in L} x_i^2 + \sum_{v_i\in R} (x_i-w)^2 + \sum_{v_i\in T} y_i^2 + \sum_{v_i\in B} (y_i-h)^2$$

对于mesh上的一个顶点$v_i$，其有四种情况，在左侧、右侧、上部，下部。当然有四个顶点同时属于两种情况。也就是说，有一些顶点需要限制$x$坐标，有一些需要限制$y$坐标，有一些两个都需要限制，而有一些则两个都不需要限制。

设$B_i$如下

$$B_i = \begin{bmatrix} a & 0 \\ 0 & b \end{bmatrix}$$

一个点如果需要限制$x$坐标，则$a=1$，如果需要限制$y$坐标，则$b=1$。否则都等于$0$。

例如左上角的顶点$(x, y)$，我们的目标是

$$\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} 0\\ 0 \end{bmatrix}$$

对于右下角的顶点，目标是

$$\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} w\\ h \end{bmatrix}$$

对于右侧的顶点，目标是

$$\begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} w\\ 0 \end{bmatrix}$$

对于中间的无限制顶点，目标是

$$\begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} x\\ y \end{bmatrix} = \begin{bmatrix} 0\\ 0 \end{bmatrix}$$

以此类推。

同样的，我们把$v_i$合并成一个$40\times 10\times 2$向量（同前）。将$B_i$也合并。

$$B = \begin{bmatrix} B_1 & & & \\ & B_2 & & \\ & & \ddots & \\ & & & B_N \end{bmatrix}$$

这里就不用重新调整$B$内容的顺序了。接下来合并目标$v_{ti}$

$$V_T = \begin{bmatrix} v_{t1} \\ v_{t2} \\ \vdots \\ v_{tn} \end{bmatrix}$$

最终的目标是

$$E_B(V) = BV = V_T$$

同样，Eigen解方程可得$V_Q$。

接下来是很复杂的$E_L$

设原始的一条线段为$\hat L=[\hat x_{l0}, \hat y_{l0}, \hat x_{l1}, \hat y_{l1}]^T$。根据我们前面的说法，一条线段会被mesh切割，然后放进各个quad中。我们用quad的四个顶点来进行双线性插值，来表示目标的$L$。方式如下，首先用双线性插值从$\hat V_q$计算出$V_q$

$$\begin{bmatrix} \hat x_0 & \hat y_0 & \hat x_0\hat y_0 & 1 & & & &\\ & & & & \hat x_0 & \hat y_0 & \hat x_0\hat y_0 & 1\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ \hat x_3 & \hat y_3 & \hat x_3\hat y_3 & 1 & & & &\\ & & & & \hat x_3 & \hat y_3 & \hat x_3\hat y_3 & 1 \end{bmatrix}F_L = \begin{bmatrix} x_0\\ y_0\\ \vdots\\ x_3\\ y_3 \end{bmatrix}$$

我们把上式左边的$8\times 8$矩阵记作$W$，上式就为$WF_L=V_q$

然后用$\hat L$表示$L$，有

$$\begin{bmatrix} \hat x_{l0} & \hat y_{l0} & \hat x_{l0}\hat y_{l0} & 1 & & & &\\ & & & & \hat x_{l0} & \hat y_{l0} & \hat x_{l0}\hat y_{l0} & 1\\ \hat x_{l1} & \hat y_{l1} & \hat x_{l1}\hat y_{l1} & 1 & & & &\\ & & & & \hat x_{l1} & \hat y_{l1} & \hat x_{l1}\hat y_{l1} & 1 \end{bmatrix}F_L = \begin{bmatrix} x_{l0}\\ y_{l0}\\ x_{l1}\\ y_{l1} \end{bmatrix}$$

上式左侧$4\times 8$矩阵记作$H$，上式就为$HF_L=L$

于是

$$L = HF_L = HW^{-1}V_q$$

$H$和$W$都为已知，所以$L$表示为$V_q$的线性函数。计算$L$的向量$e$使用两个端点相减，矩阵为

$$D = \begin{bmatrix} -1 & 0 & 1 & 0\\ 0 & -1 & 0 & 1 \end{bmatrix}$$

有

$$e = DL = \begin{bmatrix} -1 & 0 & 1 & 0\\ 0 & -1 & 0 & 1 \end{bmatrix} \begin{bmatrix} x_{l0}\\ y_{l0}\\ x_{l1}\\ y_{l1} \end{bmatrix}= \begin{bmatrix} x_{l1}-x_{l0}\\ y_{l1}-y_{l0} \end{bmatrix}$$

$E_L$中的$||Ce||^2$中的$e$最终就表示为

$$\underset{(2, 1)}{e} = \underset{(2, 4)}{D}\underset{(4, 8)}{H}\underset{(8,8)}{W^{-1}}\underset{(8, 1)}{V_q}$$

而

$$C = R\hat e(\hat e^T\hat e)^{-1}\hat e^TR^T-I$$

其中$\hat e$已知，$R$只和$\theta_m$有关。可以直接算，与$V_q$无关。

同样的，我们把所有线段的$e$和其$C$拼接（具体见代码），有

$$E_L(V,\{\theta_m\}) = \dfrac{1}{N_L}Ce = \overrightarrow{0}$$

同样的，在优化$V_q$时，上式是$V_q$的线性函数，只需要用Eigen解方程即可。

最终，我们的目标是使得

$$E = E_S+\lambda_LE_L+\lambda_BE_B$$

最小，即

$$\begin{bmatrix} \dfrac{1}{N}A\\ \dfrac{\lambda_L}{N_L}CDHW^{-1}\\ \lambda_BB \end{bmatrix}V = \begin{bmatrix} \overrightarrow{0}\\ \overrightarrow{0}\\ V_T \end{bmatrix}$$

用Eigen求解出V即可。注意其中的$CDHW^{-1}$指的是每一个线段的矩阵拼接后的结果，具体拼接见代码。