基于DSP的無人機遙感影像SIFT算法設計與實現

孫 鵬，肖 經，趙海盟，劉 帆，晏 磊，3，趙紅穎

（1. 桂林電子科技大學機電工程學院，廣西桂林541004；2. 空間信息集成與3S工程應用北京市重點實驗室（北京大學），北京100871；3. 廣西高校無人機遙測重點實驗室（桂林航天工業學院），廣西桂林541004）

（?通信作者電子郵箱zhaohaimeng@guat.edu.cn）

0 引言

尺度不變特征變換（Scale-Invariant Feature Transform，SIFT）算法［1］在無人機（Unmanned Aerial Vehicle，UAV）遙感影像特征點提取和匹配領域內扮演著重要角色，隨著時間推移得到不斷的發展和完善。

數字信號處理器（Digital Signal Processor，DSP）是一種適用于密集型數據運算與實時信號處理的微處理器［2］?，F階段，在高性能DSP 平臺上實現部分復雜圖像處理算法成為了DSP 應用研究的一個熱點。在高性能DSP 平臺上實現完整的SIFT 算法，可以減輕對主計算機中央處理器（Central Processing Unit，CPU）、圖形處理器（Graphics Processing Unit，GPU）和高速緩存等資源的占用，特別是對圖像算法處理過程中的定點、浮點迭代運算相比計算機平臺處理速度更快。同時，DSP 平臺具有體積小、功耗低和便于集成等優點［3］，采用DSP 實現無人機遙感影像臨場數據快速處理具有很大的優勢。

前期，已經有一些人在DSP 上進行了SIFT 算法的部分研究。許飛等［4］在EVM6670L 平臺上運行SIFT 算法時，基于進程間通信（Inter-Process Communication，IPC）模塊實現4 核處理器協同處理；劉顏開等［5］在TMS320C6678 平臺上實現8 核協同處理SIFT 算法。以上研究，既沒有使用本文采用的DSP內核的硬件計算單元直接處理單精度浮點型像素數據的乘法計算，也沒有充分發揮DSP的軟件流水特長，對計算過程并未進行高質量的重新編排；此外，研究中處理的圖像尺寸太小，遠遠小于本文實驗采用的大尺寸無人機遙感圖像，因此無法適應于無人機組網遙感影像的臨場快速處理。

1 Rob-hess SIFT算法

SIFT算法是第一個通過穩健的描述子將一定程度的不變量與尺度、旋轉、放射變換和光照聯系起來描述局部特征的方法［6］。它被認為是計算機視覺研究的一個熱點［7-8］，也具有計算密集的特征［9-10］。

本文的參考對象為Rob-hess 基于OpenCV（Open source Computer Vision library）和C語言編譯環境編寫的通用SIFT算法程序。Rob-hess SIFT 算法實現特征點提取的主要過程如圖1所示。

圖1 Rob-hess SIFT算法的特征點提取流程Fig. 1 Feature point extraction flowchart by Rob-hess SIFT algorithm

2 SIFT算法在DSP平臺上的實現

2.1 方案設計

德州儀器C66x系列DSP內核硬件的組成如圖2所示。

圖2 C66x DSP內核Fig. 2 C66x DSP kernel

主要數據計算模塊包括支持單精度浮點型數據運算的硬件乘法器（Hardware Multiplier，MPY）、累加器和雙數據通道［11］。

為了充分發揮DSP 內核計算性能，本文重新設計了SIFT算法的圖像數據結構、軟件流水和動態數據存儲。首先，重構圖像數據結構和圖像函數使得硬件乘法器可以直接參與單精度浮點型像素數據的乘法計算；其次，采用軟件流水技術重新編排算法的循環計算，以增強算法計算的并行能力；最后，遷移算法產生的動態數據至片外第三代雙倍速率同步動態隨機存儲器（Double Data Rate 3 synchronous dynamic random access memory，DDR3），以提升動態數據的存儲空間。具體實現流程如圖3所示。

圖3 DSP平臺上SIFT算法的總體設計Fig. 3 Overall design of SIFT algorithm on DSP platform

2.2 圖像數據類型轉換

硬件計算單元直接參與數學計算可以提升計算的速度。但是，原始SIFT 算法使用char 型地址空間存儲圖像像素數據，使得DSP 內核的單精度浮點型硬件乘法器不能夠直接計算算法單精度浮點型像素數據的乘法。

為了使DSP內核的硬件乘法器能夠直接參與單精度浮點型像素數據的乘法計算，需要根據硬件乘法器的輸入、輸出特性對圖像數據結構和圖像函數進行重構，具體實現流程如圖4所示。

1）IplImage圖像數據結構的重構。

在SIFT 算法原始IplImage 圖像數據結構中，char 型、float型像素數據都存儲于char*型指針imageData 指定的地址空間。一個char 型圖像像素數據存儲于一個char 型存儲空間內，一個float 型像素數據需要存儲于4 個char 型存儲空間內。由于，DSP 內核單精度浮點型硬件乘法器不能夠參與char 型存儲空間內數據的計算，無法發揮DSP 優良的單精度浮點型數據計算能力。

為了使單精度浮點型硬件乘法器可以直接計算SIFT 算法float 型像素數據的乘法，本文對算法的IplImage 圖像數據結構進行重構。

首先，在IplImage結構中新增float*型指針imageData1，并使用該指針為float型像素數據分配存儲空間；其次，在算法浮點運算段啟用imageData1指針直接存儲像素點的float型像素數據，實現float 型像素數據的存儲空間類型由char 型轉變為float 型；最后，對圖像浮點型像素數據的計算過程進行調整，并在編譯軟件中啟用硬件乘法器處理該型像素數據的乘法計算。新IplImage 結構中圖像單精度浮點型像素數據訪問形式如下。

新IplImage 結構中圖像gray32 第row 行、第col 列像素點的浮點型像素數據訪問形式如式（1）所示:

float_val=*（gray32-＞imageData1+

gray32-＞widthstep?row+col） （1）

新IplImage 結構中圖像gray32 第n 個浮點型像素數據的訪問形式如式（2）所示:

float_val=*（gray32-＞imageData1+n） （2）其中:gray32 為一個IplImage 型結構體圖像的起始地址；row、col 和n 為被定義的int 型變量；gray32-＞widthstep 為gray32 圖像中每一行像素數據存儲時占用的存儲空間；gray32-＞imageData1為新數據結構中圖像單精度浮點型像素數據的起始地址。

2）圖像函數的重構。

原始SIFT 算法程序調用cvCreateImage 函數創建圖像結構的首地址并根據圖像頭信息為結構體的像素數據分配適當存儲空間。在圖像創建過程中，調用的子函數主要包括創建圖像頭信息的函數cvCreateImageHeader 和為像素數據分配存儲空間的函數cvCreateData。

原始圖像像素數據存儲空間的首地址由像素數據空間分配函數cvCreateData 提供。為了使新增加的imageData1 指針能夠指向合法的存儲空間，需要對函數cvCreateData 進行重構。當輸入的圖像像素數據類型為32 位float 型時，函數cvCreateData 的子函數ialloc 開始分配存儲空間并將返回的地址強制轉換為float*型；然后，將得到的float*型地址傳遞到指針imageData1。具體實現方式如式（3）所示:

img-＞imageData1=（float*）ialloc

（（size_t）img-＞imageSize） （3）

其中:img-＞imageSize 為圖像像素數據需要占據空間的大小（單位:字節），由圖像的行、列和像素點數據類型決定；img-＞imageData1 為圖像img 的單精度浮點型像素數據在新IplImage 結構中存儲后返回的起始地址；iallco 為空間分配函數，在完成空間分配后會返回void*型起始地址（該地址可以被強制轉換為其他類型）。

2.3 軟件流水編排

2.3.1 編排設計

德州儀器C66x DSP 內核具有8 路超長指令字（Very Long Instruction Word，VLIW）、浮點數據路徑［12］和8 組單精度浮點型硬件乘法器，內核主頻高達1.25 GHz（Giga Hertz），浮點性能高達20 GFLOP（Giga Floating-point Operations per Second）。由于，其指令集架構（Instruction Set Architecture，ISA）是基于VLIW 架構，每個指令加載包的長度為256 位［13］。因此，在C66x 內核DSP 平臺上使用的單指令多數據流（Single Instruction Multiple Data，SIMD）指令可以同時執行最多8 條32位指令。

使用軟件流水技術對算法程序重新進行編排，使得編譯器可以更加均衡地使用芯片內核的硬件資源；同時，部分無關聯性指令的不同執行階段可以被同時執行，縮短各指令執行時的等待間隔，可以更充分地發揮系統的并行計算能力。

軟件流水編排的實質是:降低函數指針之間的關聯性以增強計算的并行性能；將算法的子函數內嵌入主程序；對多層循環體進行簡化與展開；使用軟件流水技術將數據計算段編排到預設的pipeline 中；同時，結合CCS（Code Composer Studio）5.5 集成開發環境的DSP 軟件流水優化功能，優化流水編排的層級，最大化使用硬件資源。采用四級流水編排前后的指令執行情況如圖5所示。

圖5 指令執行對比Fig. 5 Comparison of instruction execution

在圖5 中，指令甲、乙與指令1、2、3、4、5 為處理器的相同執行指令。沒有進行軟件流水編排時，指令乙只能在指令甲完成寫回操作后才可以得到執行。對算法程序運行重新編排后，程序中指令1 的寫回、指令2 的執行、指令3 的譯碼和指令4 的預取操作可以被同時執行，縮短了指令的執行間隔。該處為流水編排后循環核的開始；指令2的寫回、指令3的執行、指令4 的譯碼、指令5 的預取也會同時執行，為循環核的一部分；循環核左側部分為流水循環填充（pipeline loop prolog）；循環核右側部分為流水循環排空（pipeline loop epilog）。圖5中，同樣的時鐘周期內，原始程序僅可以執行甲、乙兩條完整指令，而流水編排后的程序可以執行5 條完整指令，算法程序的計算速度得以提升。

流水編排的具體實現流程如圖6所示。

圖6 軟件流水編排流程Fig. 6 Flowchart of software pipeline arrangement

2.3.2 具體實現

為了減少計算過程中輸入、輸出數據指針的關聯性，在函數中使用restrict 和const 關鍵字以聲明函數非關聯性指針指向不同內存塊。當函數輸入、輸出指針存在關聯性且在計算獨立時，需要對程序的輸出指針進行調整；建立過渡內存塊，并將輸出指針指向過渡區塊；完成計算后，再將過渡區塊的內容拷貝到原始程序指定的內存塊中。鄰近插值函數變化如下。

原始源代碼:

void resizeImg（IplImage*gray，IplImage*Big）

插入關鍵字后源代碼:

void resizeImg（const IplImage*gray，IplImage*Big）

以上源代碼中，gray 為輸入指針，Big 為輸出指針。在源代碼中使用const關鍵字以聲明輸入、輸出指針為非關聯性指針且指向不同的內存塊。

由于含有子函數的循環體無法通過優化器編排為一個pipeline。因此，需要將迭代計算調用的子函數內嵌入主程序循環體內。

此外，DSP 集成開發環境的編譯器優化循環計算時，只在循環計算的內層中形成一個pipeline。因此，需要對多重循環計算進行簡化和展開，使得計算可以更加充分地被編排入預設的pipeline。圖像歸一化的計算變化如下。

原始源代碼:

以上源代碼中，row3 和col3 是被定義的int 型變量；height1 和width1 為圖像的每列、行像素點數目；alpha 為一浮點數；gray8-＞imageData 為歸一化前圖像像素數據起始地址；gray-＞imageData1 為歸一化后圖像像素數據起始地址。原始的算法計算程序只能將圖像每一行像素數據的計算納入pipeline，而簡化后的計算程序可以將圖像所有像素數據的計算都納入pipeline，從而縮短了指令間等待時間。

軟件編譯器優化設置可以改善軟件流水性能［14］。編譯器優化選項使用如下:

1）啟用Assume no irregular alias or loop be-havior，以聲明程序中沒有使用alasing技術。

2）設置Optimization level 為3，即選擇-o3以進行文件級別的優化。

3）啟用Program mode compilation，即使能-pm 以配合-o3使算法實現程序級別優化。

4）設置Specify call assumptions when opti-mizing 為3，選擇-op3以控制程序的優化級。

5）設置Optimize for code size 為2，選擇-ms2 以縮小代碼的部分尺寸。

6）設置Generate optimizer information file at level 為2，選擇-on2以生成優化信息文件。

2.4 數據存儲遷移

SIFT 算法具有占用存儲空間大的特征，處理1 000×750彩色正射影像時，主要內存占用情況如圖7所示。分析圖7可知:SIFT 算法處理1000 × 750 彩色影像時，構建高斯差分金字塔過程中需要的存儲空間最大，約為168 MB；考慮到SIFT算法計算過程中的其他變量和特征點情況，處理1000 × 750彩色影像所需要的實際存儲空間應大于168 MB。芯片內部的SHRAM（shared memory）存儲空間為64 MB，不能夠滿足算法運行過程中的存儲需求。

圖7 主要內存分布Fig. 7 Main memory distribution

為了使DSP 平臺SIFT 算法可以處理1 000×750 彩色正射影像，需要將算法的動態數據從芯片內部的SHRAM 轉移至片外DDR3存儲器。由于原始DDR3存儲器內各子存儲空間，也無法支持SIFT 算法的數據存儲，因此需要將DDR3 存儲器中不連續的原始子存儲空間進行合并，以拓展存儲器中的連續存儲空間。

為了完成上述的存儲遷移，需要調整存放鏈接器的配置信息CMD（算法工程中后綴為.cmd）文件。具體調整細節為:根據平臺板載存儲空間重新編寫CMD 文件中的SECTIONS（目標存儲器模型段）和MEMORY（硬件資源描述段）。

1）CMD文件的SECTIONS編寫。

德州儀器C66x DSP 內核的片內內存與CPU 保持相同的時鐘頻率。但是，片內存儲空間過小，不適應大量數據的存放。在本文研究的過程中，使用DSP 平臺搭載的DDR3 存儲器存儲SIFT 算法計算產生的動態數據，改變動態數據存儲空間的具體方式如下。

在新創建圖像結構體分配內存時，直接采用了ialloc 函數，同時，需要調整CMD 文件中SECTIONS（目標存儲器模型段）的描述代碼，詳細源代碼變化如下。

原始CMD文件中動態內存的存儲器模型描述源代碼:

.sysmem ＞SHRAM

重新編寫CMD 文件中動態內存的存儲器模型描述源代碼:

.sysmem ＞DDR3

2）CMD文件的MEMORY編寫。

CMD 文件的MEMORY 代碼段用于描述系統實際硬件資源配置。為了充分使用DDR3 存儲器的存儲空間，需要將存儲器內不連續的硬件描述區域合并為一個以DDR3 為名稱的連續存儲空間，具體的硬件描述源代碼變換如下。

原始DDR3存儲器地址空間的硬件描述源代碼:

DDR3:origin=0x80000000 length=0x10000000

DDR3A:origin=0x90000000 length=0x10000000

DDR3B:origin=0xA0000000 length=0x10000000

合并后DDR3存儲器地址空間的硬件描述源代碼:

DDR3:origin=0x80000000 length=0x7FFFFFFF

在以上代碼中，DDR3、DDR3A、DDR3B 分別為DDR3 存儲器內子存儲空間的名稱；origin 表示該子存儲空間的起始地址；length 表示該子存儲空間的有效長度；DSP 平臺實際搭載的DDR3 存儲器起始地址為0x80000000，有效長度為0x7FFFFFFF（容量為2 GB），該容量遠大于SHRAM 內64 MB存儲空間。

通過重新編寫CMD 文件SECTIONS 和MEMORY，將DSP硬件平臺DDR3 存儲器內不連續的存儲空間轉變為連續的存儲空間，并將SIFT 算法的動態數據遷移至該段連續的存儲空間內，為算法處理較大尺寸的圖像打下了基礎。

3 實驗與結果分析

3.1 硬件平臺設計

根據上述方案，選擇德州儀器66AKH12 處理器作為DSP硬件平臺主處理器，平臺如圖8 所示。66AKH12 是基于KeyStone Ⅱ架構的高性能處理器，主處理器內部集成了8 核的C66x 內核組和4 核的ARM（Advanced RISC Machine）A15內核組［15］。XDS（Extended Debugging Simulator）200mini 模塊為仿真器；MiniUSB（Universal Serial Bus）接口固定在仿真器上，用于與外部連接。

圖8 DSP測試平臺Fig. 8 DSP test platform

在實驗中，設置C66x 內核組的0 核作為處理器核，配置DSP 的時鐘頻率為1.2 GHz。測試平臺通過硬件Board to Board 連接器連接至XDS200mini 模塊，并使用MiniUSB 電纜將XDS200mini 連接至PC。PC 平臺的CCS 5.5 軟件集成開發環境用于程序的調試、加載和中斷等。運行設計如圖9 所示。測試圖像為一無人機遙感任務采集的正射彩色影像，原始影像尺寸為1 000 × 750，如圖10所示。

圖9 硬件運行設計Fig. 9 Hardware running design

圖10 測試圖像Fig. 10 Test image

3.2 結果分析

1）算法精度分析。

DSP為32位處理器，原始SIFT算法程序的處理平臺為64位處理器。由于，32位處理器與64位處理器數據總線位數的不同，在數據長度為64 位double 型數據的計算過程中會產生相應的計算誤差；而且，在不修改計算數據的數據類型時，該部分的計算誤差會一直存在。由于在構建高斯尺度空間的計算過程中進行了大量的double 型數據的乘法、加法和減法迭代計算，因此DSP 平臺SIFT 算法程序與原始SIFT 程序所提取特征點坐標、尺度、方向的小數部分和128 維特征描述子也會存在相應的差異。不同平臺SIFT 算法的實際處理情況如圖11所示。

圖11 不同平臺的特征點處理結果Fig. 11 Processing results of feature points on different platforms

圖11 中A、B 列分別為DSP 平臺SIFT 算法程序和原始SIFT 算法程序提取的特征點信息，C 列為DSP 平臺SIFT 算法程序提取的特征點數據與原始算法提取數據之差。在A、B列中feat1 為1 000×750 測試圖像（如圖10 所示）的名字。x、y 為提取特征點的坐標信息，scl 為特征點的尺度，ori 為特征點的方向角，ii 為特征點的編號；從第8 行到第135 行數據為特征點0 的128 維特征描述子；分析C 列可以發現，算法在不同平臺上提取特征點0 的描述子信息誤差為0，特征點0 坐標、尺度的小數部分存在細小的誤差（此誤差由處理器之間硬件差異造成）。由于EXCLE 表格中第2 行到第7 行的特征信息數據不是單一數值，因此，在進行不同平臺特征點信息逐差計算的過程中，此處系統自動使用“#VALUE！”以表示“引用單元格錯誤”。由于圖11 無法展示全體特征點的特征信息，因此使用表1對該信息進行統計。

表1 特征點信息（1000 × 750）Tab. 1 Information of feature points（1000 × 750）

分析表1 可以發現，處理測試的1000 × 750 圖像時，本研究的DSP 平臺SIFT 算法程序與原始算法程序的處理結果保持高度一致:兩個平臺提取的特征點數量和特征點坐標信息保持一致，在3 324 個相同坐標特征點所對應的特征描述子中，相同描述子比例的高達99.97%；可見DSP 平臺SIFT 算法程序對圖像特征點的識別、描述能力與原始算法程序基本保持一致。由于使用特征點信息進行圖像匹配時，處理過程對圖像特征點的坐標、尺度、方向角和128 維特征描述子的誤差要求遠大于本研究算法程序的提取誤差，因此，該算法程序完全滿足臨場實時快速處理無人機組網遙感影像的精度需求。

2）耗時分析。

DSP 平臺SIFT 算法的運行時間由CCS 5.5 軟件自帶Count Event模塊統計。詳細實驗結果如圖12所示。

圖12 DSP平臺運行的SIFT算法Fig. 12 SIFT algorithm running on DSP platform

實驗中，處理器為C66x 內核組的0 核，主頻為1.2 GHz。當程序運行至圖12 第42 行時，受到間斷點的限制，板載仿真在此處暫停；此時，Count Event模塊統計的時鐘顯示圖像預處理耗時，此處耗時不計入SIFT 算法的總耗時。點擊繼續運行后，Count Event 模塊自動清零并開始重新統計CPU 運行時鐘數據，SIFT 算法程序開始被執行；當程序運行至第46行時，板載仿真進入暫停，此時Count Event 模塊統計的時鐘顯示SIFT算法處理1000 × 750 影像所消耗的CPU 時鐘數據。由圖12分析可得，在特征點提取過程中，DSP平臺SIFT算法程序共消耗2 525 842 766 個CPU 時鐘周期，約為2.53 GHz。結合DSP主頻配置和表1 分析可得:在確保提取到高質量特征點坐標、尺度、方向和特征描述子前提下，DSP 平臺的0 核運行本文的SIFT 算法程序處理測試的1000 × 750圖像需要運行2.108 s，滿足在單核處理器平臺上臨場實時快速處理無人機組網遙感影像的時間要求。

4 結語

本文基于DSP平臺的軟硬件資源實現了具有完整功能的Rob-hess SIFT算法。通過對圖像數據結構、圖像函數的重構、軟件流水的重新編排和動態數據存儲空間的遷移，實現了對大尺寸無人機遙感圖像的高精度快速處理。實驗結果表明，本文算法的精度、運行耗時和處理圖像的尺寸滿足臨場處理需求。

下一步，可以對數據的輸入輸出接口進行深入設計，實現無人機遙感影像數據在PC 和DSP 間的快速交互。同時，在DSP 開發的過程中，進行多任務調度實現可以進一步提升處理速度。