图像预处理之warpaffine与双线性插值及其高性能实现
warpaffine矩阵变换
对于坐标点的变换，我们通常考虑的是旋转、缩放、平移这三种变换。例如将点 P ( x , y ) P(x,y) P(x,y) 旋转 θ \theta θ 度，缩放 s c a l e scale scale 倍，平移 o x , o y ox,oy ox,oy 。warpaffine 将坐标点的旋转、缩放、平移三种操作集成为一个矩阵乘法运算。
旋转变换
我们先来看旋转，如图所示，我们要将点 P ( x , y ) P(x,y) P(x,y) 旋转到点 P ′ ( x ′ , y ′ ) P'(x',y') P′(x′,y′) ，推导的过程很简单，我们要求的就是 x ′ , y ′ x',y' x′,y′ 两点的坐标，将其转换为 m × c o s ( θ + α ) m\times cos(\theta+\alpha) m×cos(θ+α) 和 m × s i n ( θ + α ) m\times sin(\theta+\alpha) m×sin(θ+α) ，再用公式展开，即得结果（详见图中公式）：
{
x
′
y
′
}
=
{
c
o
s
(
θ
)
−
s
i
n
(
θ
)
s
i
n
(
θ
)
c
o
s
(
θ
)
}
{
x
y
}
\left\{ \begin{array}{rc} x' \\ y' \end{array} \right\}= \left\{ \begin{array}{rc} cos(\theta) & -sin(\theta) \\ sin(\theta) & cos(\theta) \end{array} \right\} \left\{ \begin{array}{rc} x \\ y \end{array} \right\}
{x′y′}={cos(θ)sin(θ)−sin(θ)cos(θ)}{xy}
再考虑到我们在图像处理时的坐标系（如在 OpenCV 中的坐标系、常见目标检测的坐标系等）通常是原点在左上角，因此应该为：
{
x
′
y
′
}
=
{
c
o
s
(
θ
)
s
i
n
(
θ
)
−
s
i
n
(
θ
)
c
o
s
(
θ
)
}
{
x
y
}
\left\{ \begin{array}{rc} x' \\ y' \end{array} \right\}= \left\{ \begin{array}{rc} cos(\theta) & sin(\theta) \\ -sin(\theta) & cos(\theta) \end{array} \right\} \left\{ \begin{array}{rc} x \\ y \end{array} \right\}
{x′y′}={cos(θ)−sin(θ)sin(θ)cos(θ)}{xy}
将旋转变换的矩阵记为
R
R
R ，则
P
′
=
R
P
P'=RP
P′=RP
缩放变换
缩放变换比较简单，两坐标直接乘以缩放系数 s c a l e scale scale 即可：
x
′
=
x
×
s
c
a
l
e
y
′
=
y
×
s
c
a
l
e
x'=x\times scale \\ y'=y\times scale
x′=x×scaley′=y×scale
写成矩阵形式即：
{
x
′
y
′
}
=
{
s
c
a
l
e
0
0
s
c
a
l
e
}
{
x
y
}
\left\{ \begin{array}{rc} x' \\ y' \end{array} \right\}= \left\{ \begin{array}{rc} scale & 0 \\ 0 & scale \end{array} \right\} \left\{ \begin{array}{rc} x \\ y \end{array} \right\}
{x′y′}={scale00scale}{xy}
将缩放变换的变换矩阵记为
S
S
S，则：
P
′
=
S
P
P'=SP
P′=SP
则旋转+缩放可以通过矩阵相乘写到同一个矩阵中：
{
x
′
y
′
}
=
{
c
o
s
(
θ
)
×
s
c
a
l
e
s
i
n
(
θ
)
×
s
c
a
l
e
−
s
i
n
(
θ
)
×
s
c
a
l
e
c
o
s
(
θ
)
×
s
c
a
l
e
}
{
x
y
}
\left\{ \begin{array}{rc} x' \\ y' \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}{rc} x \\ y \end{array} \right\}
{x′y′}={cos(θ)×scale−sin(θ)×scalesin(θ)×scalecos(θ)×scale}{xy}
即：
P
′
=
S
R
P
P'=SRP
P′=SRP
注意旋转和缩放顺序是随意的，不影响结果，这也可以通过代码来验证：
import numpy as np
theta = 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
平移变换
平移变换可以表示为：
x
′
=
x
+
o
x
y
′
=
y
+
o
y
x'=x+ox\\ y'=y+oy
x′=x+oxy′=y+oy
矩阵形式：
{
x
′
y
′
}
=
{
1
0
0
1
}
{
x
y
}
+
{
o
x
o
y
}
\left\{ \begin{array}{rc} x' \\ y' \end{array} \right\}= \left\{ \begin{array}{rc} 1 & 0 \\ 0 & 1 \end{array} \right\} \left\{ \begin{array}{rc} x \\ y \end{array} \right\} + \left\{ \begin{array}{rc} ox \\ oy \end{array} \right\}
{x′y′}={1001}{xy}+{oxoy}
可以发现，平移变换直接写成矩阵形式，已经不是单纯的矩阵相乘了，而是多了一个很麻烦的相加的操作。这就很难与我们之前的缩放+旋转的操作合并到一起，该怎么办呢？
我们可以增加一个维度，将二维的非齐次的形式转换为三维的齐次的形式，即这个知乎回答中所提到的：增加一个维度之后，就可以在高维度通过线性变换来完成低维度的放射变换。（该回答将放射变换讲的很形象，推荐阅读）。
那么我们增加一维 ( x , y , w ) (x,y,w) (x,y,w)，从而将点 P P P 表示为 P ( x w , y w , 1 ) P(\frac{x}{w},\frac{y}{w},1) P(wx,wy,1) ，这样平移变换就也可以表示为齐次矩阵乘的形式：
{
x
′
y
′
w
}
=
{
1
0
o
x
0
1
o
y
0
0
1
}
{
x
y
1
}
\left\{ \begin{array}{rc} x' \\ y' \\ w \\ \end{array} \right\}= \left\{ \begin{array}{rc} 1 & 0 & ox \\ 0 & 1 & oy \\ 0 & 0 & 1 \end{array} \right\} \left\{ \begin{array}{rc} x \\ y \\ 1 \end{array} \right\}
⎩⎨⎧x′y′w⎭⎬⎫=⎩⎨⎧100010oxoy1⎭⎬⎫⎩⎨⎧xy1⎭⎬⎫
最后我们得到缩放+旋转+平移变换的矩阵表示（注意平移与缩放、旋转的顺序是不能随意调换的）：
{
x
′
y
′
w
}
=
{
c
o
s
(
θ
)
×
s
c
a
l
e
s
i
n
(
θ
)
×
s
c
a
l
e
o
x
−
s
i
n
(
θ
)
×
s
c
a
l
e
c
o
s
(
θ
)
×
s
c
a
l
e
o
y
0
0
1
}
{
x
y
1
}
\left\{ \begin{array}{rc} x' \\ y' \\ w \end{array} \right\}= \left\{ \begin{array}{rc} 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}{rc} x \\ y \\ 1 \end{array} \right\}
⎩⎨⎧x′y′w⎭⎬⎫=⎩⎨⎧cos(θ)×scale−sin(θ)×scale0sin(θ)×scalecos(θ)×scale0oxoy1⎭⎬⎫⎩⎨⎧xy1⎭⎬⎫
将平移变换的变换矩阵记为
R
R
R ，则：
P
′
=
T
S
R
P
P'=TSRP
P′=TSRP ，可以将整个 warpaffine 三个变换操作的矩阵记为
M
M
M ，即：
M
=
T
S
R
,
P
′
=
M
P
M=TSR,\ \ P'=MP
M=TSR, P′=MP 。
warpaffine矩阵变换的反变换
- 旋转矩阵的逆矩阵，即是其转置： R − 1 = R T R^{-1}=R^T R−1=RT
- 整个 warp affine 的三个变换求反变换，对整个变换矩阵求逆即可： P ′ = M P , P = M − 1 P ′ P'=MP,\ \ P=M^{-1}P' P′=MP, P=M−1P′
目标检测中的常用预处理
在目标检测中，我们的预处理通常是先对图像进行等比缩放，然后居中，多余部分填充，就类似下图所展示的。
我们将这个过程分为三个步骤：
- 等比缩放，矩阵 S S S 实现
- 将图片中心平移到左上坐标原点，矩阵 O O O 实现
- 将图片平移到目标位置的重心，矩阵 T T T 实现
三步拆分法，看似麻烦了一点，实际上可以方便我们后续可能会需要到的更复杂的变换（比如在 O O O 平移后加入旋转变换），并且便于记忆。
三步拆分法的矩阵表达： P ′ = T O S P P'=TOSP P′=TOSP 。
我们直接写出具体的矩阵：
s
c
a
l
e
=
m
i
n
(
D
s
t
.
w
i
d
t
h
O
r
i
g
i
n
.
w
i
d
t
h
,
D
s
t
.
h
e
i
g
h
t
O
r
i
g
i
n
.
h
e
i
g
h
t
)
M
=
{
s
c
a
l
e
0
−
s
c
a
l
e
×
O
r
i
g
i
n
.
w
i
d
t
h
2
+
D
s
t
.
w
i
d
t
h
2
0
s
c
a
l
e
−
s
c
a
l
e
×
O
r
i
g
i
n
.
h
e
i
g
h
t
2
+
D
s
t
.
h
e
i
g
h
t
2
}
scale = min(\frac{Dst.width}{Origin.width}, \frac{Dst.height}{Origin.height}) \\ \\ M = \left\{ \begin{array}{ll} 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\}
scale=min(Origin.widthDst.width,Origin.heightDst.height)M={scale00scale−2scale×Origin.width+2Dst.width−2scale×Origin.height+2Dst.height}
{ x ′ y ′ } = { s c a l e 0 − s c a l e × O r i g i n . w i d t h 2 + D s t . w i d t h 2 0 s c a l e − s c a l e × O r i g i n . h e i g h t 2 + D s t . h e i g h t 2 } { x y 1 } \left\{ \begin{array}{ll} x' \\ y' \\ \end{array} \right\}= \left\{ \begin{array}{ll} 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}{ll} x \\ y \\ 1 \end{array} \right\} {x′y′}={scale00scale−2scale×Origin.width+2Dst.width−2scale×Origin.height+2Dst.height}⎩⎨⎧xy1⎭⎬⎫
逆变换：
k
=
s
c
a
l
e
b
1
=
−
s
c
a
l
e
×
O
r
i
g
i
n
.
w
i
d
t
h
2
+
D
s
t
.
w
i
d
t
h
2
b
2
=
−
s
c
a
l
e
×
O
r
i
g
i
n
.
h
e
i
g
h
t
2
+
D
s
t
.
h
e
i
g
h
t
2
x
′
=
k
x
+
b
1
y
′
=
k
y
+
b
2
x
=
x
′
−
b
1
k
=
x
′
×
1
k
+
(
−
b
1
k
)
y
=
y
′
−
b
2
k
=
y
′
×
1
k
+
(
−
b
2
k
)
M
−
1
=
{
1
k
0
−
b
1
k
0
1
k
−
b
2
k
}
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' = kx + b1 \\ y' = ky + b2 \\ x = \frac{x' - b1}{k} = x'\times \frac{1}{k} + (-\frac{b1}{k}) \\ y = \frac{y' - b2}{k} = y'\times \frac{1}{k} + (-\frac{b2}{k}) \\ M^{-1} = \left\{ \begin{array}{ll} \frac{1}{k} & 0 & -\frac{b1}{k} \\ 0 & \frac{1}{k} & -\frac{b2}{k} \\ \end{array} \right\}
k=scaleb1=−2scale×Origin.width+2Dst.widthb2=−2scale×Origin.height+2Dst.heightx′=kx+b1y′=ky+b2x=kx′−b1=x′×k1+(−kb1)y=ky′−b2=y′×k1+(−kb2)M−1={k100k1−kb1−kb2}
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_size
dst = 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] < ih
for y in range(dh):
for x in range(dw):
homogeneous = np.array([[x, y, 1]]).T
ox, oy = M @ homogeneous
low_ox = int(np.floor(ox))
low_oy = int(np.floor(oy))
high_ox = low_ox + 1
high_oy = low_oy + 1
# p0 p1
# o
# p2 p3
pos = ox - low_ox, oy - low_oy
p0_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_oy
p1 = high_ox, low_oy
p2 = low_ox, high_oy
p3 = high_ox, high_oy
p0_value = image[p0[1], p0[0]] if irange(p0) else constant
p1_value = image[p1[1], p1[0]] if irange(p1) else constant
p2_value = image[p2[1], p2[0]] if irange(p2) else constant
p3_value = image[p3[1], p3[0]] if irange(p3) else constant
dst[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行代码实现focus
return dst
cat1 = 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])
