造轮子之仿射变换

有人说，我们不应该再造轮子；也有人说，学习怎么造轮子可以带来更深的理解。我说，用轮子有用轮子的乐趣，造轮子有造轮子的乐趣。总之，开心就好。

仿射变换，可能是图像处理中最为基本的变换之一了吧。OpenCV中就有现成的函数，大可以一调了之。然而，如果不能使用OpenCV呢？比如，在DSP上实现。有人会说，TI的DSP直接有一个医学图像库可以调用仿射变换啊。好吧，我被打败了，但TI的医学图像库的仿射变换有一篇文档。那篇文档写得实在是太感人肺腑了，难道我们不该顺便学习一下吗？所以有了这篇文章（其实大多数就是翻译一下啦，不用那么认真。直接看原文吧：））

简介

图像中要使用仿射变换的目的其实就是为了把一个本来已经偏了的东西放正。比如，我之前做的身份证扫描的项目。身份证通过CIS的扫描，变成了一个背景是黑色，前景当然就是身份证的东西了。（透露这么多已经很危险了，就不放图了）可以想象一下，身份证不一定就很正吧，会有些许的倾斜。那问题来了，客户想要的是标标准准的正置着的裁好边的身份证图像。我该怎么做呢？很简单，就是分成两个步骤：

1. 找到身份证的四个顶点的坐标。
2. 利用仿射变换对这四个顶点坐标进行仿射变换。（其实只要三个顶点坐标就OK了。）

举这么个例子，想象力比较丰富的可能就知道仿射变换是要干什么了。仿射变换其实就是图像变换的一列组合，包括：旋转(rotating)、移位(shifting)、缩放(scaling)、切变(shearing)。

图像的变换究竟是怎么变换的呢？事实上，图像的变换其实就是坐标点位置的变换。正变换是这样：知道了原图，知道了变换方式，我们就可以推测出目标图像是哪样。但其实我们使用的往往是逆变换。也就是：我们首先对目标图像有一个大小的期望(分辨率为$W\times H$)，那么每一个像素该填什么值呢？就可以通过逆变换，找到其在原图像上的对应位置。此时，看看原图上的周围像素点是多少，就可以用双线性插值等算法计算出目标图像上该点的像素值了。

基本操作

我们可以写出下面的式子来表示仿射坐标变换的通式：

$$p_T=A\cdot p_S$$

其中，

$$p_S=\left[\begin{array}{c}{x}_S\\ y_S\\ 1\end{array}\right]$$

$$p_T=\left[\begin{array}{c}{x}_T\\ y_T\\ 1\end{array}\right]$$

分别是原点位置和目标点位置。

$$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ 0& 0& 1\end{array}\right]$$

是变换矩阵。这个矩阵可以拆分为如下几个矩阵：

平移矩阵：

$$A_T=\left[\begin{array}{ccc}1& 0& d_x\\ 0& 1& d_y\\ 0& 0& 1\end{array}\right]$$

这个平移矩阵的作用就是将图像平移。只要代入原式，就可以看到，$x_T=x_S+d_x$，$y_T=y_S+d_y$。

对图像放缩的变换矩阵如下：

$$A_S=\left[\begin{array}{ccc}{\beta }_x& 0& 0\\ 0& {\beta }_y& 0\\ 0& 0& 1\end{array}\right]$$

对图像旋转的变换矩阵如下：

$$A_R=\left[\begin{array}{ccc}\cos \theta & -\sin \theta & 0\\ \sin \theta & \cos \theta & 0\\ 0& 0& 1\end{array}\right]$$

其中，$\theta$为旋转的角度

切变也是一种线性变换：

$$A_Q=\left[\begin{array}{ccc}1& s_x& 0\\ s_y& 1& 0\\ 0& 0& 1\end{array}\right]$$

其中，$s_x$和$s_y$分别是$x$方向和$y$方向的切变系数。

原图

沿x方向切变

旋转

综合

算法细节

一般地，我们会去查看目标上的点在原图像上的位置，并用原图像上对应点的像素值插值后得到目标图像上该点的像素值。通过目标点的位置来计算原点的位置使用的是逆变换：

$$p_S=A^{-1}\cdot p_T$$
$$A^{-1}=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ 0& 0& 1\end{array}\right]$$

在下图中，蓝色格子表示目标图像，红色格子表示原图像。

目标图像的(0,0)点对应的是原图像的(xshift,yshift)点。从而：

$$\left[\begin{array}{c}\text{xshift}\\ \text{yshift}\\ 1\end{array}\right]=\left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$$

因此，$b_{13}=xshift,b_{23}=yshift$。

令[xstep_c,ystep_c]为目标图像移动一列后，原图像的坐标发生的相应的变化。则：

KaTeX parse error: Expected 'EOF', got '_' at position 37: …{c} \text{xstep_̲c} \\ \text{yst…

因此，$b_{11}=xste{p}_c,b_{21}=yste{p}_c$。

类似地，令[xstep_r,ystep_r]为目标图像移动一行后，原图像的坐标发生的相应的变化。则：

KaTeX parse error: Expected 'EOF', got '_' at position 37: …{c} \text{xstep_̲r} \\ \text{yst…

这意味着，$b_{12}=xste{p}_r,b_{22}=yste{p}_r$。

因此，原图像的位置可以这样计算得到：

KaTeX parse error: Expected 'EOF', got '_' at position 106: …cc} \text{xstep_̲c} & \text{xste…

实现

以上描述的关系如下图所示：

伪代码如下：

双线性插值

源代码

下面是在C语言中实现的源代码片段

void affine_warp_8q(const uint8_t * RESTRICT in_img, uint16_t src_rows, uint16_t src_cols, 
				    const AffineWarpParamq_t *pWarp, 
				    uint8_t* RESTRICT dst_img, short dst_rows, short dst_cols)
{
    uint16_t count_row, count_col; // loop counters
    int xstart_r;
    int ystart_r;
    int xf, yf;
    int ind;
    uint16_t cxf, cyf;
    uint16_t oneminuscxf, oneminuscyf;
    uint32_t s01f, s23f, soutf;
    uint8_t s0, s1, s2, s3;

    uint16_t xi, yi;
    int xstep_cq = pWarp->xstep_c;
    int ystep_cq = pWarp->ystep_c;
    int xstep_rq = pWarp->xstep_r;
    int ystep_rq = pWarp->ystep_r;

    /* Set the starting point in source corresponding to
    * top left corner of target block.
    */
    xstart_r = pWarp->xshift;
    ystart_r = pWarp->yshift;

    for (count_row = 0; count_row < dst_rows; count_row++)
    {
        // Set the starting point for the given row
        xf = xstart_r;
        yf = ystart_r;
        for (count_col = 0; count_col < dst_cols; count_col++)
        {
            xi = (uint16_t)(xf >> AFFINEWARP_ISHIFT);
            yi = (uint16_t)(yf >> AFFINEWARP_ISHIFT);

            /* extract most significant 16 bits  */
            cxf = (uint16_t)xf;
            cyf = (uint16_t)yf;
            oneminuscxf = 0xFFFF - cxf;
            oneminuscyf = 0xFFFF - cyf;

            // Get four input pixels for interpolation
            ind = yi * (src_cols) + xi;
            s0 = in_img[ind++];
            s1 = in_img[ind];
            ind += src_cols - 1;
            s2 = in_img[ind++];
            s3 = in_img[ind];

            // Interpolate
            s01f = oneminuscxf * s0 + cxf * s1;
            s23f = oneminuscxf * s2 + cxf * s3;

            s01f >>= 8;
            s23f >>= 8;
            soutf = oneminuscyf * s01f + cyf * s23f;
            soutf >>= 24;

            // set output pixel
            dst_img[count_row * dst_cols + count_col] = (uint8_t)soutf;

            // Add source increments for advancing by one column
            xf += xstep_cq;
            yf += ystep_cq;
        }
        // Add source increment for advancing by one row.
        xstart_r += xstep_rq;
        ystart_r += ystep_rq;
    }
}

/* affine_warp_8f
 * warp affine the pSource image to pTarget image
 * pSource - pointer to source image
 * pTarget - pointer to dest image
 * pWarp   - affine parameters
 */
void affine_warp_8f(const ImgBlock8_t *pSource, const AffineWarpParamf_t *pWarp, ImgBlock8_t *pTarget)
{
    ushort count_row, count_col; // loop counters
    float xstart_r;
    float ystart_r;
    float xf, yf;
    float cxf, cyf;
    float s01f, s23f, soutf;
    uchar s0, s1, s2, s3;
    uchar * restrict out_val;

    ushort xi, yi;

    int *restrict pf;
    int i_soutf, i_xi, i_yi, E;
    const int zero = 0;

    // Set the starting point in source corresponding to
    // top left corner of target block.
    xstart_r = pWarp->xshift;
    ystart_r = pWarp->yshift;

    #pragma MUST_ITERATE(400, 400, 4)
    for (count_row = 0; count_row < pTarget->num_rows; ++count_row)
    {
        // Set the starting point for the given row
        xf = xstart_r;
        yf = ystart_r;

        #pragma MUST_ITERATE(800, 800, 4)
        for (count_col = 0; count_col < pTarget->num_cols; ++count_col)
        {
            if (xf < 1.0)
                xi = 0;
            else
            {
                pf = (int*)&xf;
                E = ((*pf >> 23) & 0xff) - 127;
                i_xi = *pf & 0x7fffff;
                i_xi >>= 23 - E;
                i_xi |= (1 << E);

                xi = (ushort)(i_xi);
            }

            if (yf < 1.0)
                yi = 0;
            else
            {
                pf = (int*)&yf;
                E = ((*pf >> 23) & 0xff) - 127;
                i_yi = *pf & 0x7fffff;
                i_yi >>= 23 - E;
                i_yi |= (1 << E);

                yi = (ushort)(i_yi);
            }

            cxf = xf - xi;
            cyf = yf - yi;

            s0 = GET_PIXEL(pSource, yi, xi);
            s1 = GET_PIXEL(pSource, yi, xi + 1);
            s2 = GET_PIXEL(pSource, yi + 1, xi);
            s3 = GET_PIXEL(pSource, yi + 1, xi + 1);

            s01f = (1.0 - cxf) * s0 + cxf * s1;
            s23f = (1.0 - cxf) * s2 + cxf * s3;
            soutf = (1.0 - cyf) * s01f + cyf * s23f;

            soutf = soutf + 0.5 + 1;

            pf = (int*)&soutf;
            E = ((*pf >> 23) & 0xff) - 127;
            *pf &= (~zero - ((1 << (23 - E)) - 1));
            soutf -= 1.;
            if (soutf > 255.0) soutf = 255.0;

            if (soutf < 1.0)
                out_val = 0;
            else
            {
                pf = (int*)&soutf;
                E = ((*pf >> 23) & 0xff) - 127;
                i_soutf = *pf & 0x7fffff;
                i_soutf >>= 23 - E;
                i_soutf |= (1 << E);

                out_val = (uchar*)i_soutf;
            }

            pTarget->pBase[count_row * pTarget->row_stride + count_col] = *out_val;

            xf += pWarp->xstep_c;
            yf += pWarp->ystep_c;
        }

        xstart_r += pWarp->xstep_r;
        ystart_r += pWarp->ystep_r;
    }
}

其中，warp_affine_8q是定点数的版本；warp_affine_8f是浮点数的版本。

性能说明

最后有必要说明一下性能。其实也就两句话：

1. 近邻域插值方法比双线性插值方法快；
2. 定点比浮点快。（只要是DSP，大概都这样吧）