Python正面硬刚C语言，结果会怎样？

水深无声 2022-09-08 11:57 182阅读 0赞

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_Q1NETiBA5aSp5YWD5rWq5a2Q_size_118_color_FFFFFF_t_70_g_se_x_16_pic_center]

雷达扫描数据样例（扫描半径1km，扫描范围130°）

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

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_Q1NETiBA5aSp5YWD5rWq5a2Q_size_85_color_FFFFFF_t_70_g_se_x_16_pic_center]

原始数据转扇形图像（顺时针扫描，初始相位155°）

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

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

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

# -*- coding:utf-8 -*-
    
    import os, time
    import numpy as np
    from PIL import Image
    
    def outimg(fn_squ, fn_fan, angle, r0=0, phase=180, cw=True):
        """将矩形图像转为环形 fn_squ - 输入文件名 fn_fan - 输出文件名 angle - 环形夹角度数 r0 - 环形内圆半径，默认r0为0，输出扇形 phase - 初始相位（原点在输出图像的中心，以指向右侧的水平线为0°，逆时针方向为正） cw - 顺时针扫描 """
        
        im = np.array(Image.open(fn_squ)) # 读图像文件为NumPy数组
        h, w, d = im.shape # 矩形图像的高度、宽度和通道数
        
        r1 = h + r0 # 扇形半径
        k = int(np.ceil(np.radians(angle)*r1/w)) # 插值系数（自动确定，无需修改）
        
        xs = np.ones((2*r1-1, 2*r1-1), dtype=np.int32) * -1 # 列索引数组
        ys = np.ones((2*r1-1, 2*r1-1), dtype=np.int32) * -1 # 行索引数组
        rs = np.linspace(r1, r0, h) # 半径序列
        hs = range(h) # 行序列
        
        if cw: # 顺时针扫描
            theta = np.radians(np.linspace(phase, phase-angle, k*w))
        else: # 逆时针扫描
            theta = np.radians(np.linspace(phase, phase+angle, k*w))
        
        for i in range(k*w):
            x = np.int32(np.cos(theta[i])*rs) + r1 - 1
            y = -np.int32(np.sin(theta[i])*rs) + r1 + 1
            
            xs[(y, x)] = i//k
            ys[(y, x)] = hs
        
        im_fan = im[(ys, xs)] # 从原始数据得到扇形图像数据
        im_fan[np.where(xs == -1)] = (0,0,0,0) # 空白部分置为透明
        Image.fromarray(im_fan).save(fn_fan) # 保存为文件
        
        
    if __name__ == '__main__':
        fn_squ = 'res/raw_d130_1km.png'
        fn_fan = 'res/fan_d130_1km.png'
        
        t0 = time.time()
        outimg(fn_squ, fn_fan, angle=130, r0=100, phase=155, cw=True)
        t1 = time.time()
        
        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][stb_GitHub]下载stb\_image.h和stb\_image\_write.h两个头文件。

#include <stdio.h>
    #include <stdlib.h>
    #include <windows.h>
    
    #define _USE_MATH_DEFINES
    #include <math.h>
    
    #define STB_IMAGE_IMPLEMENTATION
    #include "stb_image.h"
    
    #define STB_IMAGE_WRITE_IMPLEMENTATION
    #include "stb_image_write.h"
    
    int main() { 
        LARGE_INTEGER li;
        LONGLONG startTime, stopTime, freq;
        QueryPerformanceFrequency(&li);
        freq = li.QuadPart;
        QueryPerformanceCounter(&li);
        startTime = li.QuadPart; // 记录开始时间
    
        char* rawFile = "D://MyCcode//RadoData2Image//res//raw_d130_1km.png";
        char* outFile = "D://MyCcode//RadoData2Image//res//fan_d130_1km.png";
        int w_raw = 0, h_raw = 0, chn = 0;
    
        unsigned char* radoData = stbi_load(rawFile, &w_raw, &h_raw, &chn, 0);
    
        int r0 = 100, cw = 1, r1 = h_raw + r0;
        double angle = 130.0, phase = 155.0;
        int w_out = 2*r1 - 1, h_out = 2*r1 - 1;
        int size_out = w_out * h_out * chn;
        int k = (int)(ceil((M_PI*angle/180.0)*w_out/w_raw));
        int arc = k * w_raw;
        double step = angle/(arc-1);
    
        char* fanData;
        fanData = (char*)malloc(size_out); // 生成保存转换结果的数组
        for (int i=0; i<size_out; ++i) { 
            fanData[i] = 0; // 初始化像素，全部透明
        }
    
        double theta, sinv, cosv;
        int x, y, col_raw, pos_raw, pos_out;
        for (int i=0; i<arc; i++) { 
            if (cw == 1)
                theta = M_PI * (phase - i*step) / 180.0;
            else
                theta = M_PI * (phase + i*step) / 180.0;
    
            sinv = sin(theta);
            cosv = cos(theta);
            col_raw = i/k;
    
            for (int r=r0; r<r1; r++) { 
                x = (int)(cosv*r) + r1 - 1;
                y = -(int)(sinv*r) + r1 + 1;
                pos_out = (y * h_out + x) * chn;
                pos_raw = ((h_raw - 1 - r + r0) * w_raw + col_raw) * chn;
                for (int j=0; j<4; j++)
                    fanData[pos_out+j] = radoData[pos_raw+j];
            }
        }
    
        // 将转换结果保存为文件
    	stbi_write_png(outFile, w_out, h_out, chn, fanData, 0);
    
        QueryPerformanceCounter(&li);
        stopTime = li.QuadPart; // 记录结束时间
        int costTime =(int)((stopTime-startTime)*1000/freq);
        printf("Time cost: %d ms\n", costTime);
    
        return 0;
    }

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

![在这里插入图片描述][183c99acb6cc45c39d8cbce3302a73ec.png_pic_center]

C代码运行截图

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

<table> 
 <thead> 
  <tr> 
   <th align="center">No.</th> 
   <th align="center">Python</th> 
   <th align="center">C</th> 
  </tr> 
 </thead> 
 <tbody> 
  <tr> 
   <td align="center">1</td> 
   <td align="center">1635ms</td> 
   <td align="center">2596ms</td> 
  </tr> 
  <tr> 
   <td align="center">2</td> 
   <td align="center">1652ms</td> 
   <td align="center">2599ms</td> 
  </tr> 
  <tr> 
   <td align="center">3</td> 
   <td align="center">1696ms</td> 
   <td align="center">2609ms</td> 
  </tr> 
  <tr> 
   <td align="center">4</td> 
   <td align="center">1673ms</td> 
   <td align="center">2557ms</td> 
  </tr> 
  <tr> 
   <td align="center">5</td> 
   <td align="center">1633ms</td> 
   <td align="center">2550ms</td> 
  </tr> 
  <tr> 
   <td align="center">6</td> 
   <td align="center">1632ms</td> 
   <td align="center">2584ms</td> 
  </tr> 
  <tr> 
   <td align="center">7</td> 
   <td align="center">1626ms</td> 
   <td align="center">2567ms</td> 
  </tr> 
  <tr> 
   <td align="center">8</td> 
   <td align="center">1729ms</td> 
   <td align="center">2603ms</td> 
  </tr> 
  <tr> 
   <td align="center">9</td> 
   <td align="center">1642ms</td> 
   <td align="center">2562ms</td> 
  </tr> 
  <tr> 
   <td align="center">1</td> 
   <td align="center">1691ms</td> 
   <td align="center">2593ms</td> 
  </tr> 
  <tr> 
   <td align="center">平均</td> 
   <td align="center">1660ms</td> 
   <td align="center">2582ms</td> 
  </tr> 
 </tbody> 
</table>

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

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_Q1NETiBA5aSp5YWD5rWq5a2Q_size_118_color_FFFFFF_t_70_g_se_x_16_pic_center]: /images/20220829/9b67ad0710654a65ab16d3b88c25378a.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_Q1NETiBA5aSp5YWD5rWq5a2Q_size_85_color_FFFFFF_t_70_g_se_x_16_pic_center]: /images/20220829/1cbb5dfa3e5644188da4b0e8084e6a23.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_Q1NETiBA5aSp5YWD5rWq5a2Q_size_41_color_FFFFFF_t_70_g_se_x_16_pic_center]: /images/20220829/218539f3169f4a72abb18ab2862838f8.png
[stb_GitHub]: https://github.com/nothings/stb/
[183c99acb6cc45c39d8cbce3302a73ec.png_pic_center]: /images/20220829/42485a550efd43d5be9f885ebc88a11e.png