Python正面硬刚C语言,结果会怎样?

时间:2021-10-06 19:06:44

Python正面硬刚C语言,结果会怎样?

雷达数据可视化是雷达数据处理的最后阶段,通常是将一个二维数组的数据转换为扇形图像。这个二维数组的行数对应着雷达的扫描半径,扫描半径越大,行数越多;数据的列数和雷达的扫描角度相关,扫描角度越大,列数越多。

Python正面硬刚C语言,结果会怎样?

雷达扫描数据样例(扫描半径1km,扫描范围130°)

比如,上面这张图就是一个扫描半径1km、扫描范围130°的雷达二维数据的直观显示,下面这张图则是由这个数据转换得到的扇形图像。

Python正面硬刚C语言,结果会怎样?

原始数据转扇形图像(顺时针扫描,初始相位155°)

二维数据转扇形图像的原理很简单,就是将二维数据的每一列写到输出图像的对应像素上。如果输出图像的扇形弧长比原始数据的列数多,则需要插值。下图是二维数据转扇形图像的原理示意图,原始数据共有8列,而输出图像的圆弧长24个像素(弧长取决于雷达扫描角度和扫描半径),故每一列数据被重复操作了3次,每次的旋转角度各不相同。

Python正面硬刚C语言,结果会怎样?

二维数据转扇形图像原理示意图

近日有网友求援,要我帮忙优化一个用于雷达数据可视化的Python脚本。略作分析之后,基于二维数据转扇形图像的基本原理,我为求援的网友重写了一个新的脚本文件,全部代码大约50余行。

  1. # -*- coding:utf-8 -*- 
  2.  
  3. import os, time 
  4. import numpy as np 
  5. from PIL import Image 
  6.  
  7. def outimg(fn_squ, fn_fan, angle, r0=0, phase=180, cw=True): 
  8.     """将矩形图像转为环形 
  9.  
  10.     fn_squ      - 输入文件名 
  11.     fn_fan      - 输出文件名 
  12.     angle       - 环形夹角度数 
  13.     r0          - 环形内圆半径,默认r0为0,输出扇形 
  14.     phase       - 初始相位(原点在输出图像的中心,以指向右侧的水平线为0°,逆时针方向为正) 
  15.     cw          - 顺时针扫描 
  16.     ""
  17.  
  18.     im = np.array(Image.open(fn_squ)) # 读图像文件为NumPy数组 
  19.     h, w, d = im.shape # 矩形图像的高度、宽度和通道数 
  20.  
  21.     r1 = h + r0 # 扇形半径 
  22.     k = int(np.ceil(np.radians(angle)*r1/w)) # 插值系数(自动确定,无需修改) 
  23.  
  24.     xs = np.ones((2*r1-1, 2*r1-1), dtype=np.int32) * -1 # 列索引数组 
  25.     ys = np.ones((2*r1-1, 2*r1-1), dtype=np.int32) * -1 # 行索引数组 
  26.     rs = np.linspace(r1, r0, h) # 半径序列 
  27.     hs = range(h) # 行序列 
  28.  
  29.     if cw: # 顺时针扫描 
  30.         theta = np.radians(np.linspace(phase, phase-angle, k*w)) 
  31.     else: # 逆时针扫描 
  32.         theta = np.radians(np.linspace(phase, phase+angle, k*w)) 
  33.  
  34.     for i in range(k*w): 
  35.         x = np.int32(np.cos(theta[i])*rs) + r1 - 1 
  36.         y = -np.int32(np.sin(theta[i])*rs) + r1 + 1 
  37.  
  38.         xs[(y, x)] = i//k 
  39.         ys[(y, x)] = hs 
  40.  
  41.     im_fan = im[(ys, xs)] # 从原始数据得到扇形图像数据 
  42.     im_fan[np.where(xs == -1)] = (0,0,0,0) # 空白部分置为透明 
  43.     Image.fromarray(im_fan).save(fn_fan) # 保存为文件 
  44.  
  45.  
  46. if __name__ == '__main__'
  47.     fn_squ = 'res/raw_d130_1km.png' 
  48.     fn_fan = 'res/fan_d130_1km.png' 
  49.  
  50.     t0 = time.time() 
  51.     outimg(fn_squ, fn_fan, angle=130, r0=100, phase=155, cw=True
  52.     t1 = time.time() 
  53.  
  54.     print('图像已处理完并保存,耗时%d毫秒'%int((t1-t0)*1000)) 

使用上面展示的扫描半径1km、扫描范围130°的雷达二维数据(可直接下载图像文件作为测试数据),这段代码生成扇形图像大约耗时1.6秒钟。发给求援的网友之后,很快传来了反馈消息:新的脚本不但可以正常运行,速度更是提升了20倍左右。略带夸张的千恩万谢之后,这位网友又说,他们原本对优化没有抱多大期望,只想尝试一下;如果优化结果不理想的话,打算用C替换这个脚本的;现在好了,处理速度足可满足需求,无需再用C重写了。

帮忙的事情算是圆满结束了,但这位网友的话却让我萌生了一个想法:用C来实现同样的功能,究竟会比Python快多少呢?平时总听到很多人说,Python如何如何慢,何不借此问题,让Python和C来一个正面较量呢?

坐而论道,不如起而行之。几个小时之后,我写完了下面这段同样是实现二维数据转扇形图像的C代码。其中加载图像文件和保存图像文件,借用了GitHub上的一个C/C++图像库。这个名为stb的图像库,并非无名之辈,单是Contributors就有188人之多,持续开发近10年之久,圈内也算小有名气。若要运行下面的代码,请先去stb的GitHub(https://github.com/nothings/stb/)下载stb_image.h和stb_image_write.h两个头文件。

  1. #include <stdio.h> 
  2. #include <stdlib.h> 
  3. #include <windows.h> 
  4.  
  5. #define _USE_MATH_DEFINES 
  6. #include <math.h> 
  7.  
  8. #define STB_IMAGE_IMPLEMENTATION 
  9. #include "stb_image.h" 
  10.  
  11. #define STB_IMAGE_WRITE_IMPLEMENTATION 
  12. #include "stb_image_write.h" 
  13.  
  14. int main() { 
  15.     LARGE_INTEGER li; 
  16.     LONGLONG startTime, stopTime, freq; 
  17.     QueryPerformanceFrequency(&li); 
  18.     freq = li.QuadPart; 
  19.     QueryPerformanceCounter(&li); 
  20.     startTime = li.QuadPart; // 记录开始时间 
  21.  
  22.     char* rawFile = "D://MyCcode//RadoData2Image//res//raw_d130_1km.png"
  23.     char* outFile = "D://MyCcode//RadoData2Image//res//fan_d130_1km.png"
  24.     int w_raw = 0, h_raw = 0, chn = 0; 
  25.  
  26.     unsigned char* radoData = stbi_load(rawFile, &w_raw, &h_raw, &chn, 0); 
  27.  
  28.     int r0 = 100, cw = 1, r1 = h_raw + r0; 
  29.     double angle = 130.0, phase = 155.0; 
  30.     int w_out = 2*r1 - 1, h_out = 2*r1 - 1; 
  31.     int size_out = w_out * h_out * chn; 
  32.     int k = (int)(ceil((M_PI*angle/180.0)*w_out/w_raw)); 
  33.     int arc = k * w_raw; 
  34.     double step = angle/(arc-1); 
  35.  
  36.     char* fanData; 
  37.     fanData = (char*)malloc(size_out); // 生成保存转换结果的数组 
  38.     for (int i=0; i<size_out; ++i) { 
  39.         fanData[i] = 0; // 初始化像素,全部透明 
  40.     } 
  41.  
  42.     double theta, sinv, cosv; 
  43.     int x, y, col_raw, pos_raw, pos_out; 
  44.     for (int i=0; i<arc; i++) { 
  45.         if (cw == 1) 
  46.             theta = M_PI * (phase - i*step) / 180.0; 
  47.         else 
  48.             theta = M_PI * (phase + i*step) / 180.0; 
  49.  
  50.         sinv = sin(theta); 
  51.         cosv = cos(theta); 
  52.         col_raw = i/k; 
  53.  
  54.         for (int r=r0; r<r1; r++) { 
  55.             x = (int)(cosv*r) + r1 - 1; 
  56.             y = -(int)(sinv*r) + r1 + 1; 
  57.             pos_out = (y * h_out + x) * chn; 
  58.             pos_raw = ((h_raw - 1 - r + r0) * w_raw + col_raw) * chn; 
  59.             for (int j=0; j<4; j++) 
  60.                 fanData[pos_out+j] = radoData[pos_raw+j]; 
  61.         } 
  62.     } 
  63.  
  64.     // 将转换结果保存为文件 
  65.     stbi_write_png(outFile, w_out, h_out, chn, fanData, 0); 
  66.  
  67.     QueryPerformanceCounter(&li); 
  68.     stopTime = li.QuadPart; // 记录结束时间 
  69.     int costTime =(int)((stopTime-startTime)*1000/freq); 
  70.     printf("Time cost: %d ms\n", costTime); 
  71.  
  72.     return 0; 

激动人心的时刻终于到了,我迫不及待地点击了“构建并运行”按钮。看起来一切顺利,屏幕迅速滚动并最终定格。

Python正面硬刚C语言,结果会怎样?

C代码运行截图

什么?2668毫秒?竟然比Python慢了1000毫秒?不可能!!!直觉告诉我,一定是哪里出现了问题。接下来我又花了几个小时,反复检查验证,但结果和过程都没有发现问题。下表是10次运行结果的耗时记录,结果显示,在相同的测试条件下,Python平均耗时1660毫秒,C平均耗时2582毫秒,Python耗时大约是C的64%。

No. Python C
1 1635ms 2596ms
2 1652ms 2599ms
3 1696ms 2609ms
4 1673ms 2557ms
5 1633ms 2550ms
6 1632ms 2584ms
7 1626ms 2567ms
8 1729ms 2603ms
9 1642ms 2562ms
1 1691ms 2593ms
平均 1660ms 2582ms
 

尽管不可思议,但我现在开始尝试相信这个结果了。读者您呢?要是有疑问或建议,欢迎留言。如有更加高效的Python代码或着C代码,请发私信给我,让我们一起将这场Python和C语言的正面交锋延续下去、延伸开来。

原文链接:https://mp.weixin.qq.com/s/t1dLgdlg6Msm412SVS_DBA