打败opencv ,哦，是快了3倍

大家好，本文转自我一个读者朋友Homio的文章，推荐给大家，希望对做这方便的同学有所帮助。

程序员，哦！不！软件工程师们都对opencv很熟悉，它在工作学习研究中起到了不可或缺的作用。但是它臃肿的身躯和迟缓的动作让人苦不堪言。

opencv初体验

今天的主角叫图像去畸变算法。先隆重请它出场🎉🎉🎉🎉🎉

图像去畸变效果

为了更好的成像质量引入的透镜导致了图像的畸变，无论是在摄影还是在AI图像识别方面都有还原事物真实面目的需求。那么我们首先体验一下在opencv中去畸变的快感吧。（偷偷说一下上面那张图是在ps中去畸变的👂🏻）

opencv中去畸变有两种形式，一种直接调用cv::undistort()函数，函数参数如下：

@param src Input (distorted) image.
@param dst Output (corrected) image that has the same size and type as src .
@param cameraMatrix Input camera matrix \f$A = \vecthreethree{f_x}{0}{c_x}{0}{f_y}{c_y}{0}{0}{1}\f$ .
@param distCoeffs Input vector of distortion coefficients
\f$(k_1, k_2, p_1, p_2[, k_3[, k_4, k_5, k_6[, s_1, s_2, s_3, s_4[, \tau_x, \tau_y]]]])\f$
of 4, 5, 8, 12 or 14 elements. If the vector is NULL/empty, the zero distortion coefficients are assumed.
@param newCameraMatrix Camera matrix of the distorted image. By default, it is the same as
cameraMatrix but you may additionally scale and shift the result by using a different matrix.*/
CV_EXPORTS_W void undistort( InputArray src, OutputArray dst,InputArray cameraMatrix,InputArray distCoeffs,InputArray newCameraMatrix = noArray() );

此法对每帧图像去畸变都要计算复杂的公式，耗时非常多，因此在实时视频处理方面，没人用它。

另一种方式就是调用cv::initUndistortRectifyMap()和cv::remap()。该方法是先根据相机内参获取到两个映射表，（相机内参是用张正友标定法获取的，这里就不解析了😭）再用remap()函数找到非畸变图像对映畸变图像的像素，直接把这个像素拷贝过来就🆗了。因此去畸变的过程就变成了查表的过程，复杂的运算只需要做一次就够了。

那么我们测试的代码如下：

int getInternalParm(const char *datFileName, Mat &mapx, Mat &mapy)
{Size imageSize;Matx33d intrinsic_matrix;Vec4d distortion_coeffs;mapx = Mat(imageSize, CV_32FC1);mapy = Mat(imageSize, CV_32FC1);FILE *camParam = fopen(datFileName, "rb");if (NULL == camParam){std::cout << "fopen fail:" << datFileName << std::endl;return -1;}fread(&intrinsic_matrix, sizeof(cv::Matx33d), 1, camParam);fread(&distortion_coeffs, sizeof(cv::Vec4d), 1, camParam);fread(&imageSize, sizeof(Size), 1, camParam);fclose(camParam);Mat R = Mat::eye(3, 3, CV_32F);fisheye::initUndistortRectifyMap(intrinsic_matrix, distortion_coeffs, R, intrinsic_matrix, imageSize, CV_32FC1, mapx, mapy);return 0;
}int main(int argc, char *argv[])
{Mat mapx;Mat mapy;if (-1 == getInternalParm("camParam.dat", mapx, mapy)){printf("getInternalParm error!\n");return 0;}Mat frame_bgr = imread("fisheys_modle.jpg");Mat image_undistort;remap(frame_bgr, image_undistort, mapx, mapy, INTER_NEAREST);    //去畸变return 0;
}

庖丁解牛

要想把opencv扔进垃圾桶，光完成上面的代码是没有任何用的，你完全不知道去畸变是怎么回事，只能看着opencv跑出来的图像大呼：好牛逼啊。然并卵。

径向畸变切向畸变

引起图像畸变的原因有两个，一种是由于透镜形状引起的畸变，称之为径向畸变，可以看到，镜像畸变基本上都是沿着半径方向分布。另一种是由于安装导致透镜和成像平面不平行导致的畸变，称之为切向畸变。

先来了解下相机的内参矩阵：

其中fx = f/dx f为相机焦距，单位mm，dx为像素x方向宽度，单位mm，1/dx表示x方向上1mm内有多少像素。fy同理。cx=u0，表示图像中心在像素坐标系下的位置。cy同理。

四个重要的坐标系

在去畸变算法中，有四个重要的坐标系：世界坐标系、相机坐标系、成像坐标系、像素坐标系。

世界坐标系->成像坐标系

完整的坐标变换过程

真实世界上的一点M' 经过旋转平移后达到成像平面上，其坐标我们用小写的x, y, z表示。然后除以z坐标映射到归一化坐标系，最后再经过内参矩阵A转换到像素坐标系，坐标用u，v表示。其中s表示可以对图像进行缩放，我们直接不缩放，令其为1即可。

由于我们不关心世界坐标系的点，直接将上述公式1-1归一化后代入公式1-2可得：

其中x/z我们称之为归一化坐标x' ，y/z为归一化坐标y' 。那么将上式展开可得归一化坐标系和像素坐标系的转换关系：

至此，我们完成了图像坐标系的转换，此后去畸变的过程就是在归一化坐标系下完成的。比如我们拿到一张畸变图像，可以依次取出像素点，转换到归一化坐标系下。归一化坐标系下有公式对应畸变图像和非畸变图像的坐标，那么就能获取到非畸变图像在归一化坐标系下的坐标，再根据式1-4求得非畸变图像在像素坐标系下的坐标。

那么归一化坐标系下畸变图像和非畸变图像是如何对应的呢？学者们给出了两种畸变模型：

针孔模型

鱼眼模型

其中，针孔模型适用于畸变比较小的镜头，鱼眼模型适用于畸变比较大的镜头。上式中k和p参数为径向畸变系数和切向畸变系数，在标定时已经获取。只有畸变坐标和非畸变坐标未知。因此可以根据一方求得另一方。求解过程列解如下：

针孔模型

鱼眼模型

九层垒土

那么我们将上述去畸变的过程自己用C语言写一下吧！

#include <sys/time.h>
#include <string>
#include <math.h>#define DST_WIDTH 640
#define DST_HEIGHT 480
#define SRC_VER_DEV 0 /**> the source image deviation in vertical*/
#define SRC_UPR_DEV 0  /**> the source image deviation in upright*/
#define SRC_WIDTH 640
#define SRC_HEIGHT 480typedef struct distort_coor_t
{uint32_t u_distorted;uint32_t v_distorted;
}distort_coor;void NV12_to_YUV444packed(uint8_t *src_date, uint8_t *dst_data, int width, int height)
{uint8_t *src_y = src_date;uint8_t *src_uv = src_date + width * height;uint32_t tmp = 0;int k = 0;for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){dst_data[k++] = src_y[i * width + j];tmp = 2 * ((i / 2) * (width / 2) + (j / 2));dst_data[k++] = src_uv[tmp];dst_data[k++] = src_uv[tmp + 1];}}return;
}void YUV444packed_to_NV12(uint8_t *src_data, uint8_t *dst_data, int width, int height)
{uint8_t *dst_y = dst_data;uint8_t *dst_u = dst_data + width * height;uint32_t tmp = 0;int k = 0;for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){dst_y[i * width + j] = src_data[k++];tmp = 2 * ((i / 2) * (width / 2) + (j / 2));dst_u[tmp] = src_data[k++];dst_u[tmp + 1] = src_data[k++];}}return;
}/************************************ function: get the mapping table between the undistorted picture and distorted picture* input: *cam_param:   *the intrinsic param and distortion coeffs file of camera. (using opencv data file)*          mapping_table[DST_HEIGHT][DST_WIDTH]:    the mapping table*          width:   the destination picture width*          height:    the destination picture height*          type: the model of distortion   1: pinhole, 2: fisheye* return: 0:success, -1 fail. * ********************************/
int get_mapping_table(const char *cam_param, distort_coor mapping_table[DST_HEIGHT][DST_WIDTH], const uint16_t width, const uint16_t height, uint8_t mode_type)
{/**> read the intrinsic_matrix and distortion_coeffs of camera */FILE *fp = fopen(cam_param, "rb");if (NULL == fp){printf("fopen fail: %s, errno=%d\n", cam_param, errno);return -1;}/**> 畸变参数(k1, k2为径向畸变系数，p1, p2为切向畸变系数) */double k1 = 0.0, k2 = 0.0, p1 = 0.0, p2 = 0.0;/**> 内参 fx, fy是以像素为单位的焦距，cx, cy是基准点(通常在图像的中心)*/double fx = 0.0, fy = 0.0, cx = 0.0, cy = 0.0;fread(&fx, sizeof(double), 1, fp);fseek(fp, sizeof(double), SEEK_CUR);fread(&cx, sizeof(double), 1, fp);fseek(fp, sizeof(double), SEEK_CUR);fread(&fy, sizeof(double), 1, fp);fread(&cy, sizeof(double), 1, fp);fseek(fp, 3 * sizeof(double), SEEK_CUR);fread(&k1, sizeof(double), 1, fp);fread(&k2, sizeof(double), 1, fp);fread(&p1, sizeof(double), 1, fp);fread(&p2, sizeof(double), 1, fp);fclose(fp);int rows = height, cols = width;/**> 计算去畸变后图像的内容(v, u为非畸变图像的高和宽) */for (int v = 0; v < rows; v++){for (int u = 0; u < cols; u++){uint32_t u_distorted = 0, v_distorted = 0;    /**> 畸变后的图像坐标 */float x1, y1, x2, y2; /**> 其中x1, y1为正常图像在归一化坐标系的坐标，x2, y2是畸变图像在归一化坐标系的坐标 */if (1 == mode_type)   /**> 针孔畸变模型 */{/**> 将image_undistort的坐标通过内参转换到归一化坐标系下，此时得到的归一化坐标是对的 */x1 = (u - cx) / fx;y1 = (v - cy) / fy;/**> 根据畸变模型，求得在归一化坐标系下畸变图像的坐标 */float r2;r2 = pow(x1, 2) + pow(y1, 2);x2 = x1 * (1 + k1 * r2 + k2 * pow(r2, 2)) + 2 * p1 * x1 * y1 + p2 * (r2 + 2 * x1 * x1);y2 = y1 * (1 + k1 * r2 + k2 * pow(r2, 2)) + p1 * (r2 + 2 * y1 * y1) + 2 * p2 * x1 * y1;/**> 求得畸变图像在像素坐标系下的坐标 */u_distorted = (uint32_t)fx * x2 + cx;v_distorted = (uint32_t)fy * y2 + cy;}else   /**> 鱼眼畸变模型 (圆形鱼眼图像) */{/**> 将image_undistort的坐标通过内参转换到归一化坐标系下，此时得到的是正常图像的归一化坐标 */x1 = (u - (cx + SRC_VER_DEV)) / fx;y1 = (v - (cy + SRC_UPR_DEV)) / fy;/**> 根据畸变模型，求得在归一化坐标系下畸变图像的坐标 */float r2;float r;float theta;float theta1;r2 = pow(x1, 2) + pow(y1, 2);r = sqrt(r2);theta = atan(r);theta1 = theta * (1 + k1 * pow(theta, 2) + k2 * pow(theta, 4) + p1 * pow(theta, 6) + p2 * pow(theta, 8));x2 = theta1 * x1 / r;y2 = theta1 * y1 / r;/**> 求得畸变图像在像素坐标系下的坐标 */u_distorted = (uint32_t)fx * x2 + cx;v_distorted = (uint32_t)fy * y2 + cy;}            mapping_table[v][u].u_distorted = u_distorted;mapping_table[v][u].v_distorted = v_distorted;}}return 0;
}void undistort_img(uint8_t *yuv_raw, uint8_t *yuv_undistort, const distort_coor mapping_table[DST_HEIGHT][DST_WIDTH], uint16_t width, uint16_t height)
{uint32_t tmp2 = 0;for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){tmp2 = mapping_table[i][j].v_distorted * SRC_WIDTH * 3 + (int)mapping_table[i][j].u_distorted * 3;memcpy(yuv_undistort, yuv_raw + tmp2, 3);yuv_undistort += 3;}}return;
}int main(int argc, char *argv[])
{static distort_coor mapping_table[DST_HEIGHT][DST_WIDTH] = {0};uint8_t *yuvImg_NV12 = (uint8_t *)calloc(rows * cols * 3, sizeof(uint8_t));uint8_t *yuvImg_444packed = (uint8_t *)calloc(rows * cols * 3, sizeof(uint8_t));uint8_t *yuvImg_undistort = (uint8_t *)calloc(DST_WIDTH * DST_HEIGHT * 3, sizeof(uint8_t));uint8_t *yuvImg_undistort_nv12 = (uint8_t *)calloc(DST_WIDTH * DST_HEIGHT * 3, sizeof(uint8_t));if (-1 == get_mapping_table("camParam.dat", mapping_table, DST_WIDTH, DST_HEIGHT, 2)){printf("get_mapping_table error!\n");return 0;}for (int i = 0; i < 1000; i++){NV12_to_YUV444packed(yuvImg_NV12, yuvImg_444packed, DST_WIDTH, DST_HEIGHT);undistort_img(yuvImg_444packed, yuvImg_undistort, mapping_table, DST_WIDTH, DST_HEIGHT);YUV444packed_to_NV12(yuvImg_undistort, yuvImg_undistort_nv12, DST_WIDTH, DST_HEIGHT); }free(yuvImg_I420);free(yuvImg_NV12);free(yuvImg_undistort);free(yuvImg_undistort_i420);free(yuvImg_undistort_nv12);return 0;
}

上述程序只提供大概思路，保证百分之百不能运行哈，比方说相机内参文件和图片来源都没有。我的思路是，先将YUV图像转成yuv444packed的格式（不清楚yuv格式的需要自行百度喽），因为这样像RGB似的好处理。去完畸变后再转回来。（opencv同样要从yuv转到rgb格式再转回来，因为后续编码和OSD还要在yuv图像的基础上做。）去畸变的思路参考了opencv，用了查表的思想。可以看到获取映射表的函数其中计算非常复杂，进行了大量的幂计算和乘法，竟然还有8次幂。可以想象如果每帧图像都这样计算那得是多么恐怖。查表就简单多了。可以看到整个undistort_img函数不过10行，找到对应的像素计算好偏移后直接拷贝就行了。

那么我们将opencv和自己写的函数各运行1000次（两者均用最近邻插值算法），对比一下，看看效率怎么样：

测试平台： imx8q

处理器： Cortex-A35

核心数： 4

最大主频： 1.2GHz

opencv和自写c语言去畸变比较

上面是opencv处理1000张640x480图片总共耗时22.16秒，下面是自写函数c代码处理1000张640x480图片耗时76.13秒。

磨刀不误砍柴工

别急，看一下进度条，故事还长，听我娓娓道来。。。

先来了解一下arm的neon协处理器。他是用来执行simd指令的。simd叫做单指令多数据操作。具体能干啥，且看下图：

上面描绘了sisd指令和simd指令的对比，简单来讲，就是simd可以把你4条指令执行完的事我一条指令就给搞定。我理解为就是把寄存器位数扩宽了，一个寄存器里可以存好几个数。

neon中有32个这样的寄存器，每个寄存器128位。可以存1个128位的数，或者2个64位的。。。和16个8位的。好了，simd了解这些就行了。关于指令问题，我们结合代码看。对了，还需要了解下内联汇编（）：

asm asm-qualifiers ( AssemblerTemplate : OutputOperands: InputOperands: Clobbers: GotoLabels)

内联汇编的基本格式如上，被4个冒号分为5部分，第一部分是汇编语句，然后是输出部分、后面是输入部分，接着是破坏描述部分（我不知道咋翻译），最后是什么标签？（我没用过）。不需要的部分可以省略。这里只做引入，详细的使用方法请参考gnu网站。

https://gcc.gnu.org/onlinedocs/gcc/Using-Assembly-Language-with-C.html

图穷匕首现

那么，我们开始表演吧！🎭

根据上面我们自己写的代码，需要优化三个函数。先把结果代码放上去，然后我们逐句讲解：

void NV12_to_YUV444packed(uint8_t *src_date, uint8_t *dst_data, int width, int height)
{#ifdef __ARM_NEON__uint8_t *src_y = src_date;uint8_t *src_uv = src_date + width * height;asm volatile ("mov x1, %[hei] \n" //x1 for height loop"mov x2, %[src_y] \n" //x2 for src_y pointer"mov x3, %[src_uv] \n" //x3 for src_uv pointer"mov x4, %[dst_data] \n" //x4 for dst pointer"height_loop_nv12_444pack:""mov x0, %[wid] \n" //x0 for width loop"and x5, x1, #1 \n"     //get if height is odd"mul x5, x5, %[wid] \n""sub x3, x3, x5 \n""width_loop_nv12_444pack:""ld1 {v0.16b}, [x2], #16 \n" //load 16byte yuv_y to v0 register "ld2 {v3.8b, v4.8b}, [x3], #16 \n" //load 16byte yuv_uv to v3/v4 register "zip1 v1.16b, v3.16b, v3.16b \n"    //duplicate each byte in v1"zip1 v2.16b, v4.16b, v4.16b \n""st3 {v0.16b, v1.16b, v2.16b}, [x4], #48 \n"    //store v0, v1, v2 to memory"sub x0, x0, #16 \n""cmp x0, #0 \n""bgt width_loop_nv12_444pack \n""sub x1, x1, #1 \n""cmp x1, #0 \n""bgt height_loop_nv12_444pack \n": [dst_data] "+r" (dst_data): [src_y] "r" (src_y), [src_uv] "r" (src_uv), [wid] "r" (width), [hei] "r" (height): "memory", "x0", "x1", "x2", "x3", "x4", "x5", "v0", "v1", "v2", "v3", "v4");#elseuint8_t *src_y = src_date;uint8_t *src_uv = src_date + width * height;uint32_t tmp = 0;int k = 0;for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){dst_data[k++] = src_y[i * width + j];tmp = 2 * ((i / 2) * (width / 2) + (j / 2));dst_data[k++] = src_uv[tmp];dst_data[k++] = src_uv[tmp + 1];}}#endifreturn;
}void YUV444packed_to_NV12(uint8_t *src_data, uint8_t *dst_data, int width, int height)
{#ifdef __ARM_NEON__uint8_t *dst_y = dst_data;uint8_t *dst_uv = dst_data + width * height;asm volatile ("mov x1, %[hei] \n" //x1 for height loop"mov x2, %[src] \n" //x2 for src pointer"mov x3, %[dst_y] \n" //x3 for dst_y pointer"mov x4, %[dst_uv] \n" //x4 for dst_uv pointer"height_loop_444pack_nv12:""mov x0, %[wid] \n" //x0 for width loop"and x5, x1, #1 \n"     //get if height is odd"mul x5, x5, %[wid] \n""sub x4, x4, x5 \n""width_loop_444pack_nv12:""ld3 {v0.16b, v1.16b, v2.16b}, [x2], #48 \n" //load 48byte yuv_y to v0,v1,v2 register "uzp1 v1.16b, v1.16b, v1.16b \n"    //store half of u in high section of v1 "uzp1 v2.16b, v2.16b, v2.16b \n"    //store half of v in high section of v2 "st1 {v0.16b}, [x3], #16 \n"    //store v0 to dst_y memory"st2 {v1.8b, v2.8b}, [x4], #16 \n"    //store v1, v2 to dst_uv memory"sub x0, x0, #16 \n""cmp x0, #0 \n""bgt width_loop_444pack_nv12 \n""sub x1, x1, #1 \n""cmp x1, #0 \n""bgt height_loop_444pack_nv12 \n": [dst_y] "+r" (dst_y), [dst_uv] "+r" (dst_uv): [src] "r" (src_data), [wid] "r" (width), [hei] "r" (height): "memory", "x0", "x1", "x2", "x3", "x4", "x5", "v0", "v1", "v2");#elseuint8_t *dst_y = dst_data;uint8_t *dst_u = dst_data + width * height;uint32_t tmp = 0;int k = 0;for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){dst_y[i * width + j] = src_data[k++];tmp = 2 * ((i / 2) * (width / 2) + (j / 2));dst_u[tmp] = src_data[k++];dst_u[tmp + 1] = src_data[k++];}}#endifreturn;
}void undistort_img(uint8_t *yuv_raw, uint8_t *yuv_undistort, const distort_coor mapping_table[DST_HEIGHT][DST_WIDTH], uint16_t width, uint16_t height)
{#ifdef __ARM_NEON__uint32_t tmp[4] = {0};tmp[0] = 3;tmp[1] = SRC_WIDTH;asm volatile ("mov x1, %[hei] \n" //x1 for height loop"mov x2, %[map] \n" //x2 for map pointer"mov x3, %[yuv444_dst] \n" //x3 for dst pointer"ld1 {v2.4s}, [%[tmp]] \n" "height_loop:""mov x0, %[wid] \n" //x0 for width loop"width_loop:""ld2 {v0.4s, v1.4s}, [x2], #32 \n" //v0 is u_distort, v1 is v_distort "mul v0.4s, v0.4s, v2.4s[0] \n""mul v1.4s, v1.4s, v2.4s[0] \n""mul v1.4s, v1.4s, v2.4s[1] \n""add v0.4s, v0.4s, v1.4s \n""mov w4, v0.4s[0] \n"   //get the offset of the first pixel"ldr x4, [%[yuv_raw], x4]\n""str x4, [x3], #3 \n"   //store the first pixel's yuv into memory"mov w4, v0.4s[1] \n""ldr x4, [%[yuv_raw], x4]\n""str x4, [x3], #3 \n""mov w4, v0.4s[2] \n""ldr x4, [%[yuv_raw], x4]\n""str x4, [x3], #3 \n""mov w4, v0.4s[3] \n""ldr x4, [%[yuv_raw], x4]\n""str x4, [x3], #3 \n""sub x0, x0, #4 \n""cmp x0, #0 \n""bgt width_loop \n""sub x1, x1, #1 \n""cmp x1, #0 \n""bgt height_loop \n": [yuv444_dst] "+r" (yuv_undistort): [yuv_raw] "r" (yuv_raw), [wid] "r" (width), [hei] "r" (height), [map] "r" (mapping_table), [tmp] "r" (tmp): "memory", "x0", "x1", "x2", "x3", "x4", "v0", "v1", "v2");#elseuint32_t tmp2 = 0;for (int i = 0; i < height; i++){for (int j = 0; j < width; j++){tmp2 = mapping_table[i][j].v_distorted * SRC_WIDTH * 3 + (int)mapping_table[i][j].u_distorted * 3;memcpy(yuv_undistort, yuv_raw + tmp2, 3);yuv_undistort += 3;}}#endifreturn;
}

我们来看undistort_img这个函数，竟然把10行代码写了这么多，哈啊哈

分析C代码，就是通过映射表找到非畸变图像对应畸变图像的坐标，然后计算出非畸变图像对于基地址的偏移，将yuv三个分量的数据拷贝到非畸变图像即可。用个图表示如下：

查表过程

从图中可知，畸变图像和非畸变图像的像素不是一一对应的，实际上是毫无规律的，这个得查映射表，根据映射表的坐标值算出偏移才能执行拷贝。（不然我直接memcpy宽x高x3不就行了嘛）。看程序流程，每次循环过程中都需要计算一个非畸变图像的偏移值。虽然这个计算量和创建映射表时进行的高次幂计算没法比。但是架不住每个像素都要计算一次啊，有乘有加的。这就轮到simd显摆了。

可以看到，虽然非畸变图像对应畸变图像的内存不是连续的。但是映射表的内存是连续的啊。非畸变（0，0）像素的映射表信息肯定存在映射表的前8个字节。（0，1）像素的往后排。这样我们就可以一次读取好几个像素的映射表。一次计算出好几个偏移，然后再分别拷贝呗。

说干就干！

输出部分[yuv444_dst] "+r" (yuv_undistort)表示作为输出结果，左边的变量在汇编代码中用，右边是c代码。+号表示可读写。r表示数据作为通用寄存器。输入部分[yuv_raw] "r" (yuv_raw), [wid] "r" (width), [hei] "r" (height), [map] "r" (mapping_table), [tmp] "r" (tmp)表示输入的原始数据。没有修饰符号修饰表示只读。破坏描述部分"memory", "x0", "x1", "x2", "x3", "x4", "v0", "v1", "v2"我理解为只要是你使用到的寄存器和内存，都要在这里显式的指明以便编译器知道。

汇编循环

图中我圈出来的功能很好理解，是完成c语言中的for循环。里面用到了mov、sub、cmp和bgt指令。我们一个一个说。

"mov x1, %[hei] \n"表示将变量height的值放入x1寄存器中。%[]跟上输入部分的变量，就相当于使用这个变量。在内联汇编中，每一行汇编语句都要跟个换行。

"sub x0, x0, #4 \n"表示将x0寄存器减4再放到x0中。这里可以看出来我们一次处理了4个像素。

"cmp x0, #0 \n"表示将x0中的值和立即数0比较，比较结果会影响spsr寄存器中的标志位N、S、C、V，关于SPSR想要了解更多，请自行翻阅arm手册。

SPSR寄存器

"bgt width_loop \n"就好理解了，就是如果上步比较的结果大于0，就跳到标签width_loop ，在内联汇编中，同一进程内，标签不允许重名，即使在不同函数、不同文件内也不行。

好了，循环的就介绍完了。

这两句是将映射表和目的地址的指针放到x2和x3寄存器中，没啥好说的吧。

介绍这两句需要先看学一下ld1和ld2指令的意义。

LD2 {V0.8H, V1.8H}, [X0]

ld2是simd指令。可以看到ld2将内存中数据取出，交叉放到v0和v1寄存器中。那么也可以理解ld3是取出放到3个v寄存器中。这样可以区分想要的数据。上述代码"ld2 {v0.4s, v1.4s}, [x2], #32 \n"，就是把映射表中的数据取出，放到v0和v1寄存器中。至于.4s表示的是这个寄存器存的是4个字长的数据（就是4个4字节的数据，刚好128bit）。刚好映射表中存的也是4字节的数据，u，v坐标以此存储。那么v0中存的就是4个像素的u坐标，v1中存的就是4个像素的v坐标。这样就实现u，v坐标的分离，方便单独操作。至于v2中存的是我想要用来当乘数的3和图像宽度。因为simd的乘法指令貌似不允许操作立即数。ld2最后面跟的立即数32表示的是读取完内存后，将指针往后偏移32个字节，以方便下次读取。因为有两个寄存器，每个寄存器存放16字节，所以要偏移32字节呀。刚才说x2中存放的是指针。那么[x2]就相当于解引用了，找到地址对应的数据。

后面的乘法指令和加法指令完成的就是c语言中tmp2的功能，将偏移量最后放到了v0寄存器中。其中v0.4s[0]表示v0寄存器的第一个通道，存储的也就是第一个像素的偏移值。同理还有二、三、四通道。

哦对了，ld1由于只有一个寄存器，那就不交叉存放啦，就直接读取16字节放到寄存器中就行啦。

这一部分是将刚才计算的偏移量存放到w4寄存器中，w4是x4的低半部分，因为是64位cpu，因此x4是64bit。所以w4是使用的低32bit，高位置零。

"ldr x4, [%[yuv_raw], x4]\n"是将原图中相对于基地址偏移x4地址上的数据加载到x4寄存器中。也就是具体像素的yuv数值啦。虽然yuv只占3个字节，但是我们取出来了8个字节。没关系，我们只要看前3个字节就好了。因为后面的数据会被下一次执行时覆盖掉。

"str x4, [x3], #3 \n"表示的是，将上面取出来的3字节yuv数值放到目的地址中去，然后将目的地址往后偏移3个字节，以便下次存放。（看吧，我就说会被覆盖掉吧）

再往后就是重复的操作了，因为我们刚刚计算了4个像素的。因此我们后面重复操作了3次，只是改了v0的通道而已。

至此，汇编代码就分析完啦。由于篇幅问题，关于yuv色彩空间模型的转换就不分析啦，感兴趣的可以留个作业，自己分析分析

那么我们看下速度吧！

opencv和自写汇编语言去畸变比较

同样的上面是opencv处理1000张的结果，耗时基本没变。下面是用汇编+SIMD加速的处理1000张的耗时。可以看到速度是opencv的三倍多！是自身c语言的11倍多！！！

后记

软件性能调优从来不是一件容易的事，汇编代码现在用的人很少了，又因为和cpu紧密相关，因此即使你认识人懂汇编的，他也不一定能帮到你。比方说我在色彩空间模型转换中用到的这三句：

"and x5, x1, #1 \n" 
"mul x5, x5, %[wid] \n"
"sub x4, x4, x5 \n"

其中x1存放的是循环中控制高度循环次数的数据，我想做的效果就是当x1为奇数时，x4就减去一个宽度值（内存指针回退），本来使用tst指令加subeq（条件指令）两句代码就能完成了。但是我的编译怎么也不支持subeq。后来找了半天才在arm开发手册中看到这句：

The A64 instruction set does not support conditional execution for every instruction. Predicated execution of instructions does not offer sufficient benefit to justify its significant use of opcode space.

原来在v7架构中，绝大部分指令都支持条件指令，这还是它们宣传的一大卖点。但是v8中就给取消了。真真给跪了。