mat opencv 修改roi_为什么opencv实现的图像几何变换这么快

想必所有学过数字图像处理的同学们当年都动手写过图像的几何变换，也就是resize，rotate，crop，warp affine和warp perpspective这些。也许，更爱学习的同学还实现过piecewise affine和image distortion。

即便今天，深度学习当道图像领域，这些方法依然被广泛用在image augmentation，cnn 后处理（crop bbox）中，甚至作为cnn中的某些层，比如upsampling2d，cropping2d甚至spatial transfomer。

当年在实现，也许好学的你还用openmp和sse写过这些算法的并行版本，比同学们的算法快了 三四倍以上，可是你却发现不管你怎么写，都快不过opencv的版本，总要慢两倍以上。今天我们就来探究一下opencv的内部实现。

通常的实现如下，通过逆变换遍历目标图像的每个像素点投影回原图像进行插值。使用opencv的数据结构Mat实现如下，为简单起见，这里只考虑单通道图像 ，

void warp_affine_f32(const cv::Mat &src_img, cv::Mat &dst_img, std::vector<float> trans_mat) {
	const uint8_t *src = src_img.ptr<uint8_t>();
	uint8_t *dst = dst_img.ptr<uint8_t>();

	int out_h = dst_img.rows;
	int out_w = dst_img.cols;
	int in_h = src_img.rows;
	int in_w = src_img.cols;
	std::memset(dst, 0, out_h * out_w);
	float m00 = trans_mat[0], m01 = trans_mat[1], m02 = trans_mat[2], m10 = trans_mat[3], m11 = trans_mat[4],
		m12 = trans_mat[5];
	for (int dy = 0; dy < out_h; dy++) {
		uint8_t *dst_row = dst + dy * out_w;
		for (int dx = 0; dx < out_w; dx++) {
		float xf = m00 * dx + m01 * dy + m02;
		float yf = m10 * dx + m11 * dy + m12;
		int ix = floor(xf);
		int iy = floor(yf);
		int ix2 = ix + 1;
		int iy2 = iy + 1;
		float b = xf - ix;
		float a = yf - iy;

		float f00 = (1.0f - a) * (1.0f - b);
		float f01 = a * (1.0f - b);
		float f10 = (1.0f - a) * b;
		float f11 = a * b;

		int idx1 = iy * in_w + ix;
		int idx2 = iy2 * in_w + ix;
		int idx3 = iy * in_w + ix2;
		int idx4 = iy2 * in_w + ix2;
		uint8_t p1 = 0, p2 = 0, p3 = 0, p4 = 0;
		if (ix < in_w && ix >= 0 && iy < in_h && iy >= 0)
			p1 = src[idx1];
		if (ix < in_w && ix >= 0 && iy2 < in_h && iy2 >= 0)
			p2 = src[idx2];
		if (ix2 < in_w && ix2 >= 0 && iy < in_h && iy >= 0)
			p3 = src[idx3];
		if (ix2 < in_w && ix2 >= 0 && iy2 < in_h && iy2 >= 0)
			p4 = src[idx4];
		float tmp = f00 * p1 + f10 * p2 + f01 * p3 + f11 * p4;
		dst_row[dx] = (uint8_t)(tmp);
		}
	}
}

投影回来的点要找离自己最近的四个点计算插值系数，这几个距离值就是a，b，1-a和1-b ，插值系数就分别是a*b，（1-a）*b，a*（1-b）和（1-a）*（1-b）

大家已经发现了中间的计算结果都是用float32 来表示的，然而我们最终的输出数据都是uint8，所以计算过程明显是精度过高了，一定有一些不影响最终结果的小数被round舍弃掉了。所以opencv做了这样一步优化，上面四个距离值都分布在0和1之间，那么就把这个区间分成256份，那四个插值系数的可能组合就有256*256种组合就行了，这种思想其实就是采样的过程。计算插值系数就变成了查表的过程，具体就是将插值系数存成一个二维查找表即可，距离值乘以256 然后floor就是索引值，另外很关键的一点，这张表只用算一次，完全可以声明为静态变量。双三次插值的做法是类似的，opencv对应源码如下，

const int INTER_REMAP_COEF_BITS=15;
const int INTER_REMAP_COEF_SCALE=1 << INTER_REMAP_COEF_BITS;

static uchar NNDeltaTab_i[INTER_TAB_SIZE2][2];

static float BilinearTab_f[INTER_TAB_SIZE2][2][2];
static short BilinearTab_i[INTER_TAB_SIZE2][2][2];


static inline void interpolateLinear( float x, float* coeffs ) {
    coeffs[0] = 1.f - x;
    coeffs[1] = x;
}


static void initInterTab1D(int method, float* tab, int tabsz) {
    float scale = 1.f/tabsz;
    if( method == INTER_LINEAR ) {
        for( int i = 0; i < tabsz; i++, tab += 2 )
            interpolateLinear( i*scale, tab );
    }
    else
        CV_Error( CV_StsBadArg, "Unknown interpolation method" );
}

所有查找表都是在整个opencv初始化的时候就建立好了，tabsz这个参数比较关键，可以理解为采样间隔，每种插值方法都有两张表，一个浮点的，另一个用于下面即将要讲到的顶点化方法。

这是第一步优化，第二步因为结果是8位无符号整数，一定不需要32位这么高的精度，16位即可，因为输出是整数(如果处理器支持float16，把所有float32全部变成float16也是可以work的)，也可以不需要浮点计算，int16即可，所以修改后的版本如下，大家可以自己试一试，

void warp_affine_i16(const cv::Mat &src_img, cv::Mat &dst_img, std::vector<float> trans_mat) {
  const uint8_t *src = src_img.ptr<uint8_t>();
  uint8_t *dst = dst_img.ptr<uint8_t>();

  int out_h = dst_img.rows;
  int out_w = dst_img.cols;
  int in_h = src_img.rows;
  int in_w = src_img.cols;

  const int INTER_BITS = 5;
  const int INTER_BITS2 = INTER_BITS * 2;
  const int INTER_TAB_SIZE = (1 << INTER_BITS);

  std::memset(dst, 0, out_h * out_w);

  float m00 = trans_mat[0], m01 = trans_mat[1], m02 = trans_mat[2], m10 = trans_mat[3], m11 = trans_mat[4],
        m12 = trans_mat[5];
  int16_t m00_16 = m00 * INTER_TAB_SIZE, m01_16 = m01 * INTER_TAB_SIZE, m02_16 = m02 * INTER_TAB_SIZE,
          m10_16 = m10 * INTER_TAB_SIZE, m11_16 = m11 * INTER_TAB_SIZE, m12_16 = m12 * INTER_TAB_SIZE;
  for (int dy = 0; dy < out_h; dy++) {
    uint8_t *dst_row = dst + dy * out_w;
    for (int dx = 0; dx < out_w; dx++) {
      int16_t xf_16 = m00_16 * dx + m01_16 * dy + m02_16;
      int16_t yf_16 = m10_16 * dx + m11_16 * dy + m12_16;
      int ix = (xf_16 >> INTER_BITS);
      int iy = (yf_16 >> INTER_BITS);
      int ix2 = ix + 1;
      int iy2 = iy + 1;
      int16_t b_16 = xf_16 - (ix << INTER_BITS);
      int16_t a_16 = yf_16 - (iy << INTER_BITS);
      int16_t f00_16 = ((INTER_TAB_SIZE - a_16) * (INTER_TAB_SIZE - b_16));
      int16_t f01_16 = (a_16 * (INTER_TAB_SIZE - b_16));
      int16_t f10_16 = ((INTER_TAB_SIZE - a_16) * b_16);
      int16_t f11_16 = (a_16 * b_16);
      int idx1 = iy * in_w + ix;
      int idx2 = iy2 * in_w + ix;
      int idx3 = iy * in_w + ix2;
      int idx4 = iy2 * in_w + ix2;
      uint8_t p1 = 0, p2 = 0, p3 = 0, p4 = 0;
      if (ix < in_w && ix >= 0 && iy < in_h && iy >= 0)
        p1 = src[idx1];
      if (ix < in_w && ix >= 0 && iy2 < in_h && iy2 >= 0)
        p2 = src[idx2];
      if (ix2 < in_w && ix2 >= 0 && iy < in_h && iy >= 0)
        p3 = src[idx3];
      if (ix2 < in_w && ix2 >= 0 && iy2 < in_h && iy2 >= 0)
        p4 = src[idx4];
      int32_t tmp = f00_16 * p1 + f10_16 * p2 + f01_16 * p3 + f11_16 * p4;
      dst_row[dx] = (uint8_t)(tmp >> INTER_BITS2);
    }
  }
}

这是个比较粗糙的版本，在这个基础上，使用多线程和simd，仔细调寄存器的使用，在用上前面的trick，最终性能应该能和opencv一战。

最后，这个算法体现了数字图像和数字信号处理算法的一个特点，输出输出都是可以量化的，理应中间计算也可以量化，而这算是最早的一个例子，让人不禁感慨，image processing is designed for hardware。

希望spatial transfomer有一天也能成为cnn的标准组件，设计真的很精巧，让人看到了传统图像算法的精巧。