[[419970]]

雷達數(shù)據(jù)可視化是雷達數(shù)據(jù)處理的最后階段，通常是將一個二維數(shù)組的數(shù)據(jù)轉(zhuǎn)換為扇形圖像。這個二維數(shù)組的行數(shù)對應(yīng)著雷達的掃描半徑，掃描半徑越大，行數(shù)越多;數(shù)據(jù)的列數(shù)和雷達的掃描角度相關(guān)，掃描角度越大，列數(shù)越多。

雷達掃描數(shù)據(jù)樣例(掃描半徑1km，掃描范圍130°)

比如，上面這張圖就是一個掃描半徑1km、掃描范圍130°的雷達二維數(shù)據(jù)的直觀顯示，下面這張圖則是由這個數(shù)據(jù)轉(zhuǎn)換得到的扇形圖像。

原始數(shù)據(jù)轉(zhuǎn)扇形圖像(順時針掃描，初始相位155°)

二維數(shù)據(jù)轉(zhuǎn)扇形圖像的原理很簡單，就是將二維數(shù)據(jù)的每一列寫到輸出圖像的對應(yīng)像素上。如果輸出圖像的扇形弧長比原始數(shù)據(jù)的列數(shù)多，則需要插值。下圖是二維數(shù)據(jù)轉(zhuǎn)扇形圖像的原理示意圖，原始數(shù)據(jù)共有8列，而輸出圖像的圓弧長24個像素(弧長取決于雷達掃描角度和掃描半徑)，故每一列數(shù)據(jù)被重復(fù)操作了3次，每次的旋轉(zhuǎn)角度各不相同。

二維數(shù)據(jù)轉(zhuǎn)扇形圖像原理示意圖

近日有網(wǎng)友求援，要我?guī)兔?yōu)化一個用于雷達數(shù)據(jù)可視化的Python腳本。略作分析之后，基于二維數(shù)據(jù)轉(zhuǎn)扇形圖像的基本原理，我為求援的網(wǎng)友重寫了一個新的腳本文件，全部代碼大約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): 
    """將矩形圖像轉(zhuǎn)為環(huán)形 
 
    fn_squ      - 輸入文件名 
    fn_fan      - 輸出文件名 
    angle       - 環(huán)形夾角度數(shù) 
    r0          - 環(huán)形內(nèi)圓半徑，默認r0為0，輸出扇形 
    phase       - 初始相位（原點在輸出圖像的中心，以指向右側(cè)的水平線為0°，逆時針方向為正） 
    cw          - 順時針掃描 
    """ 
 
    im = np.array(Image.open(fn_squ)) # 讀圖像文件為NumPy數(shù)組 
    h, w, d = im.shape # 矩形圖像的高度、寬度和通道數(shù) 
 
    r1 = h + r0 # 扇形半徑 
    k = int(np.ceil(np.radians(angle)*r1/w)) # 插值系數(shù)（自動確定，無需修改） 
 
    xs = np.ones((2*r1-1, 2*r1-1), dtype=np.int32) * -1 # 列索引數(shù)組 
    ys = np.ones((2*r1-1, 2*r1-1), dtype=np.int32) * -1 # 行索引數(shù)組 
    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)] # 從原始數(shù)據(jù)得到扇形圖像數(shù)據(jù) 
    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°的雷達二維數(shù)據(jù)(可直接下載圖像文件作為測試數(shù)據(jù))，這段代碼生成扇形圖像大約耗時1.6秒鐘。發(fā)給求援的網(wǎng)友之后，很快傳來了反饋消息：新的腳本不但可以正常運行，速度更是提升了20倍左右。略帶夸張的千恩萬謝之后，這位網(wǎng)友又說，他們原本對優(yōu)化沒有抱多大期望，只想嘗試一下;如果優(yōu)化結(jié)果不理想的話，打算用C替換這個腳本的;現(xiàn)在好了，處理速度足可滿足需求，無需再用C重寫了。

幫忙的事情算是圓滿結(jié)束了，但這位網(wǎng)友的話卻讓我萌生了一個想法：用C來實現(xiàn)同樣的功能，究竟會比Python快多少呢?平時總聽到很多人說，Python如何如何慢，何不借此問題，讓Python和C來一個正面較量呢?

坐而論道，不如起而行之。幾個小時之后，我寫完了下面這段同樣是實現(xiàn)二維數(shù)據(jù)轉(zhuǎn)扇形圖像的C代碼。其中加載圖像文件和保存圖像文件，借用了GitHub上的一個C/C++圖像庫。這個名為stb的圖像庫，并非無名之輩，單是Contributors就有188人之多，持續(xù)開發(fā)近10年之久，圈內(nèi)也算小有名氣。若要運行下面的代碼，請先去stb的GitHub(https://github.com/nothings/stb/)下載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); // 生成保存轉(zhuǎn)換結(jié)果的數(shù)組 
    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]; 
        } 
    } 
 
    // 將轉(zhuǎn)換結(jié)果保存為文件 
    stbi_write_png(outFile, w_out, h_out, chn, fanData, 0); 
 
    QueryPerformanceCounter(&li); 
    stopTime = li.QuadPart; // 記錄結(jié)束時間 
    int costTime =(int)((stopTime-startTime)*1000/freq); 
    printf("Time cost: %d ms\n", costTime); 
 
    return 0; 
} 

激動人心的時刻終于到了，我迫不及待地點擊了“構(gòu)建并運行”按鈕。看起來一切順利，屏幕迅速滾動并最終定格。

C代碼運行截圖

什么?2668毫秒?竟然比Python慢了1000毫秒?不可能!!!直覺告訴我，一定是哪里出現(xiàn)了問題。接下來我又花了幾個小時，反復(fù)檢查驗證，但結(jié)果和過程都沒有發(fā)現(xiàn)問題。下表是10次運行結(jié)果的耗時記錄，結(jié)果顯示，在相同的測試條件下，Python平均耗時1660毫秒，C平均耗時2582毫秒，Python耗時大約是C的64%。

盡管不可思議，但我現(xiàn)在開始嘗試相信這個結(jié)果了。讀者您呢?要是有疑問或建議，歡迎留言。如有更加高效的Python代碼或著C代碼，請發(fā)私信給我，讓我們一起將這場Python和C語言的正面交鋒延續(xù)下去、延伸開來。