雷达数据可视化是雷达数据处理的最后阶段, 通常是将一个二维数组的数据转换为扇形图像. 这个二维数组的行数对应着雷达的扫描半径, 扫描半径越大, 行数越多; 数据的列数和雷达的扫描角度相关, 扫描角度越大, 列数越多.
雷达扫描数据样例 (扫描半径 1km, 扫描范围 130°)
比如, 上面这张图就是一个扫描半径 1km, 扫描范围 130° 的雷达二维数据的直观显示, 下面这张图则是由这个数据转换得到的扇形图像.
原始数据转扇形图像 (顺时针扫描, 初始相位 155°)
二维数据转扇形图像的原理很简单, 就是将二维数据的每一列写到输出图像的对应像素上. 如果输出图像的扇形弧长比原始数据的列数多, 则需要插值. 下图是二维数据转扇形图像的原理示意图, 原始数据共有 8 列, 而输出图像的圆弧长 24 个像素 (弧长取决于雷达扫描角度和扫描半径), 故每一列数据被重复操作了 3 次, 每次的旋转角度各不相同.
二维数据转扇形图像原理示意图
近日有网友求援, 要我帮忙优化一个用于雷达数据可视化的 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(https://github.com/nothings/stb/) 下载 stb_image.h 和 stb_image_write.h 两个头文件.
- #include
- #include
- #include
- #define _USE_MATH_DEFINES
- #include
- #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
- fanData[i] = 0; // 初始化像素, 全部透明
- }
- double theta, sinv, cosv;
- int x, y, col_raw, pos_raw, pos_out;
- for (int i=0; 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
- 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;
- }
激动人心的时刻终于到了, 我迫不及待地点击了 "构建并运行" 按钮. 看起来一切顺利, 屏幕迅速滚动并最终定格.
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 语言的正面交锋延续下去, 延伸开来.
来源: http://developer.51cto.com/art/202108/679539.htm