图像预处理之warpaffine与双线性插值及其高性能实现

视频讲解：https://www.bilibili.com/video/BV1ZU4y1A7EG

代码Repo：https://github.com/shouxieai/tensorRT_Pro

本文为视频讲解的个人笔记。

warpaffine矩阵变换

对于坐标点的变换，我们通常考虑的是旋转、缩放、平移这三种变换。例如将点 $P(x,y)$ 旋转 $\theta$ 度，缩放 $scale$ 倍，平移 $ox,oy$ 。warpaffine 将坐标点的旋转、缩放、平移三种操作集成为一个矩阵乘法运算。

旋转变换

我们先来看旋转，如图所示，我们要将点 $P(x,y)$ 旋转到点 $P^{\prime }(x^{\prime },y^{\prime })$ ，推导的过程很简单，我们要求的就是 $x^{\prime },y^{\prime }$ 两点的坐标，将其转换为 $m\times cos(\theta +\alpha )$ 和 $m\times sin(\theta +\alpha )$ ，再用公式展开，即得结果（详见图中公式）：

$$\left\{\begin{array}{r}{x}^{\prime }\\ y^{\prime }\end{array}\right\}=\left\{\begin{array}{rc}cos(\theta )& -sin(\theta )\\ sin(\theta )& cos(\theta )\end{array}\right\}\left\{\begin{array}{r}x\\ y\end{array}\right\}$$
再考虑到我们在图像处理时的坐标系（如在 OpenCV 中的坐标系、常见目标检测的坐标系等）通常是原点在左上角，因此应该为：

$$\left\{\begin{array}{r}{x}^{\prime }\\ y^{\prime }\end{array}\right\}=\left\{\begin{array}{rc}cos(\theta )& sin(\theta )\\ -sin(\theta )& cos(\theta )\end{array}\right\}\left\{\begin{array}{r}x\\ y\end{array}\right\}$$
将旋转变换的矩阵记为 $R$ ，则 $P^{\prime }=RP$

缩放变换

缩放变换比较简单，两坐标直接乘以缩放系数 $scale$ 即可：

$$\begin{gathered} x^{\prime }=x\times scale \\ {y}^{\prime }=y\times scale \end{gathered}$$
写成矩阵形式即：

$$\left\{\begin{array}{r}{x}^{\prime }\\ y^{\prime }\end{array}\right\}=\left\{\begin{array}{rc}scale& 0\\ 0& scale\end{array}\right\}\left\{\begin{array}{r}x\\ y\end{array}\right\}$$
将缩放变换的变换矩阵记为 $S$，则：
$$P^{\prime }=SP$$
则旋转+缩放可以通过矩阵相乘写到同一个矩阵中：
$$\left\{\begin{array}{r}{x}^{\prime }\\ y^{\prime }\end{array}\right\}=\left\{\begin{array}{rc}cos(\theta )\times scale& sin(\theta )\times scale\\ -sin(\theta )\times scale& cos(\theta )\times scale\end{array}\right\}\left\{\begin{array}{r}x\\ y\end{array}\right\}$$
即：$P^{\prime }=SRP$

注意旋转和缩放顺序是随意的，不影响结果，这也可以通过代码来验证：

import numpy as nptheta = 0.8
scale = 2
rot = np.array([[np.cos(theta), np.sin(theta)],[-np.sin(theta), np.cos(theta)]
])sca = np.array([[scale, 0],[0, scale]
])print(np.allclose(rot @ sca, sca @ rot))
# 输出：True

平移变换

平移变换可以表示为：
$$\begin{gathered} x^{\prime }=x+ox \\ {y}^{\prime }=y+oy \end{gathered}$$
矩阵形式：
$$\left\{\begin{array}{r}{x}^{\prime }\\ y^{\prime }\end{array}\right\}=\left\{\begin{array}{rc}1& 0\\ 0& 1\end{array}\right\}\left\{\begin{array}{r}x\\ y\end{array}\right\}+\left\{\begin{array}{r}ox\\ oy\end{array}\right\}$$
可以发现，平移变换直接写成矩阵形式，已经不是单纯的矩阵相乘了，而是多了一个很麻烦的相加的操作。这就很难与我们之前的缩放+旋转的操作合并到一起，该怎么办呢？

我们可以增加一个维度，将二维的非齐次的形式转换为三维的齐次的形式，即这个知乎回答中所提到的：增加一个维度之后，就可以在高维度通过线性变换来完成低维度的放射变换。（该回答将放射变换讲的很形象，推荐阅读）。

那么我们增加一维 $(x,y,w)$，从而将点 $P$ 表示为 $P(\frac{x}{w},\frac{y}{w},1)$ ，这样平移变换就也可以表示为齐次矩阵乘的形式：

$$\left\{\begin{array}{r}{x}^{\prime }\\ y^{\prime }\\ w\end{array}\right\}=\left\{\begin{array}{rcc}1& 0& ox\\ 0& 1& oy\\ 0& 0& 1\end{array}\right\}\left\{\begin{array}{r}x\\ y\\ 1\end{array}\right\}$$
最后我们得到缩放+旋转+平移变换的矩阵表示（注意平移与缩放、旋转的顺序是不能随意调换的）：

$$\left\{\begin{array}{r}{x}^{\prime }\\ y^{\prime }\\ w\end{array}\right\}=\left\{\begin{array}{rcc}cos(\theta )\times scale& sin(\theta )\times scale& ox\\ -sin(\theta )\times scale& cos(\theta )\times scale& oy\\ 0& 0& 1\end{array}\right\}\left\{\begin{array}{r}x\\ y\\ 1\end{array}\right\}$$
将平移变换的变换矩阵记为 $R$ ，则：$P^{\prime }=TSRP$ ，可以将整个 warpaffine 三个变换操作的矩阵记为 $M$ ，即：$M=TSR,P^{\prime }=MP$ 。

warpaffine矩阵变换的反变换

• 旋转矩阵的逆矩阵，即是其转置：$R^{-1}=R^T$
• 整个 warp affine 的三个变换求反变换，对整个变换矩阵求逆即可：$P^{\prime }=MP,P=M^{-1}{P}^{\prime }$

目标检测中的常用预处理

在目标检测中，我们的预处理通常是先对图像进行等比缩放，然后居中，多余部分填充，就类似下图所展示的。

我们将这个过程分为三个步骤：

1. 等比缩放，矩阵 $S$ 实现

2. 将图片中心平移到左上坐标原点，矩阵 $O$ 实现

3. 将图片平移到目标位置的重心，矩阵 $T$ 实现

三步拆分法，看似麻烦了一点，实际上可以方便我们后续可能会需要到的更复杂的变换（比如在 $O$ 平移后加入旋转变换），并且便于记忆。

三步拆分法的矩阵表达：$P^{\prime }=TOSP$ 。

我们直接写出具体的矩阵：
$$\begin{gathered} scale=min(\frac{Dst.width}{Origin.width},\frac{Dst.height}{Origin.height}) \\ M=\left\{\begin{array}{llc}scale& 0& -\frac{scale\times Origin.width}{2}+\frac{Dst.width}{2}\\ 0& scale& -\frac{scale\times Origin.height}{2}+\frac{Dst.height}{2}\end{array}\right\} \end{gathered}$$

$$\left\{\begin{array}{l}{x}^{\prime }\\ y^{\prime }\end{array}\right\}=\left\{\begin{array}{llc}scale& 0& -\frac{scale\times Origin.width}{2}+\frac{Dst.width}{2}\\ 0& scale& -\frac{scale\times Origin.height}{2}+\frac{Dst.height}{2}\end{array}\right\}\left\{\begin{array}{l}x\\ y\\ 1\end{array}\right\}$$

逆变换：
$$\begin{gathered} k=scale \\ b1=-\frac{scale\times Origin.width}{2}+\frac{Dst.width}{2} \\ b2=-\frac{scale\times Origin.height}{2}+\frac{Dst.height}{2} \\ {x}^{\prime }=kx+b1 \\ {y}^{\prime }=ky+b2 \\ x=\frac{x^{\prime }-b1}{k}=x^{\prime }\times \frac{1}{k}+(-\frac{b1}{k}) \\ y=\frac{y^{\prime }-b2}{k}=y^{\prime }\times \frac{1}{k}+(-\frac{b2}{k}) \\ {M}^{-1}=\left\{\begin{array}{llc}\frac{1}{k}& 0& -\frac{b1}{k}\\ 0& \frac{1}{k}& -\frac{b2}{k}\end{array}\right\} \end{gathered}$$

warpaffine正逆变换代码实验

TODO

双线性插值

线性插值

距离目标点越远，影响就越小，因此权重是对面的距离占比。

如目标点距离冷水 0.6，距离热水 0.4，则冷水权重为 0.4 ，热水权重为 0.6 。

p0 = 20    # 冷水
p1 = 100   # 热水
pos = 0.6  # 应该多少度value = (1 - pos) * p0 + pos * p1
print(value)

双线性插值

线性插值的二维版本，原理一直，只是权重从计算长度占比改为计算面积占比。

调色板，红点对目标点（紫点）的影响权重即为对面的面积（红框面积）占总面积的比例。

高性能实现

为什么高性能？

• 我们在操作每个像素的过程中，可以将模型需要的像素级预处理（如减均值除标准差、除以255、BGR通道转换等）一并做了，避免多个操作分开来反复对每个像素进行循环访问这种低效行为。
• warpaffine 极其适合通过 cuda 核函数进行 GPU 加速。可以参考 repo 中的 preprocess_kernel.cu 。完整代码比较长这里就不放了。
• 以下是 warpaffine 双线性插值的 Python 实现，供参考：

def pyWarpAffine(image, M, dst_size, constant=(0, 0, 0)):# 注意输入的M矩阵格式，是Origin->Dst# 而这里需要的是Dst->Origin，所以要取逆矩阵M = cv2.invertAffineTransform(M)constant = np.array(constant)ih, iw   = image.shape[:2]dw, dh   = dst_sizedst      = np.full((dh, dw, 3), constant, dtype=np.uint8)irange   = lambda p: p[0] >= 0 and p[0] < iw and p[1] >= 0 and p[1] < ihfor y in range(dh):for x in range(dw):homogeneous = np.array([[x, y, 1]]).Tox, oy = M @ homogeneouslow_ox = int(np.floor(ox))low_oy = int(np.floor(oy))high_ox = low_ox + 1high_oy = low_oy + 1# p0     p1#      o# p2     p3pos = ox - low_ox, oy - low_oyp0_area = (1 - pos[0]) * (1 - pos[1])p1_area = pos[0] * (1 - pos[1])p2_area = (1 - pos[0]) * pos[1]p3_area = pos[0] * pos[1]p0 = low_ox, low_oyp1 = high_ox, low_oyp2 = low_ox, high_oyp3 = high_ox, high_oyp0_value = image[p0[1], p0[0]] if irange(p0) else constantp1_value = image[p1[1], p1[0]] if irange(p1) else constantp2_value = image[p2[1], p2[0]] if irange(p2) else constantp3_value = image[p3[1], p3[0]] if irange(p3) else constantdst[y, x] = p0_area * p0_value + p1_area * p1_value + p2_area * p2_value + p3_area * p3_value# 交换bgr  rgb# normalize ->  -mean /std# 1行代码实现normalize , /255.0# bgr bgr bgr -> bbb ggg rrr# focus# focus offset, 1行代码实现focusreturn dstcat1 = cv2.imread("cat1.png")
#acat1_cv, M, inv = align(cat1, (100, 100))
M = cv2.getRotationMatrix2D((0, 0), 30, 0.5)
acat1_cv = cv2.warpAffine(cat1, M, (100, 100))
acat1_py = pyWarpAffine(cat1, M, (100, 100))plt.figure(figsize=(10, 10))
plt.subplot(1, 2, 1)
plt.title("OpenCV")
plt.imshow(acat1_cv[..., ::-1])plt.subplot(1, 2, 2)
plt.title("PyWarpAffine")
plt.imshow(acat1_py[..., ::-1])