頂點式線性平均內插擬合演算法

最近用Matlab做了各種不同模型的三維數據的模擬實驗 ，不少需要對三維數據進行擬合，matlab 自帶一個擬合工具箱（cftool），確實強大，使用工具箱做數據擬合省不少事，對數據集足夠的原始模型而言，cftool能很準確的擬合出我們想要的結果。當然我們不講Matlab ，本篇通篇跟Matlab毫無關係，~~話講多了，下面開始。

一、數據原型

首先導入一張完整的數據表格，我們可以看到這是一個y=9 x=15 ，大小是9x15的矩陣，矩陣每個元素對應一個數值z，第二張圖是我們在第一張圖裡面的數據對應位置挑出來的數據。大小是3x5的有效矩陣。

這樣我們就在原始數據虛擬出了一個殘缺度高達90%左右的數據，那我們後面的內插擬合就是針對這個殘缺數據而做相應的計算。

1、概念導入

給你兩個數，（129 192）讓你在這兩個數中間線性插入一個值，那我們插入的方式其實有很多，例如在二維平麵線性 B2 =（A1+A2）/2, 或三維空間線性，我們用到如下。

A： 192 129

B ：192 157 129

2、一階擬合

經過上述的概念導入，我們接下對我們殘缺數據進行一階插值，參考下面公式

黑色框為有效的頂點數據，我們用演算法遍歷任意兩個頂點中，只要有一個是有效數據(黑框的，這點很重要，很重要，很重要，重要的事情說三遍)插一個數據(藍色框數據)，如行相鄰的129與192兩個頂點內插157，129與267兩個相鄰頂點內插186。結果如下圖所示。

第一階內插完了嗎？ 不，再看上面的數據。插完第一輪後新增相鄰頂點還需要內插，由於我的數據矩陣不大，第二輪插完後我的數據已經找不到相鄰空白點了。如果數據矩陣非常大，在演算法實現上，這個節點我們的代碼可以寫成遞歸形式，直到插完才開始退棧(效果如下)。

3、二階擬合

我們的二階擬合完全是根據一階擬合的結果來調整的，我的二階設計是通過四個有效頂點（黑色框框為有效頂點，重要的事情說一遍）內插一個綠框數據。AD對角，BC對角，參考如下公式

註：插完之後不執行一階擬合，一階是需要黑框數據支持，任意兩個藍色、綠色數據不需要。

4、三階擬合

三階插值跟二階插值基本是一樣的，唯一不一樣的地方：任意四個頂點（不限定是否黑框），參考上面公式插入。得到第二個圖。

5、刪除輔助點

最後告訴大家，上面除了黑色框框的為有效數據，其他顏色框框的為輔助數據，類似於解幾何題作的輔助線。現在把他刪掉吧。參考刪除之後的結果如下。

7、結果輸出

(此演算法是建立在對稱，規則的模型基礎上滴，暫未試過不規則模型的適應性)。走之前請點贊！

最後附上代碼實現,(僅提供相關的部分，供參考並非一一對應)

註明:Rx=30 TX=18

// 類型1 取模等於0為偶,對應原矩陣"頂點真實數據" if(i%2 == 0 && j%2 == 0){ //頂點不需要插值}// 類型2if(i%2 == 0 && j%2 != 0){ // 邊距8不考慮插值 if(j > (ucRxNumVA - 2) || (i == 0) || i > (ucTxNum - 2)) { // 邊距8不考慮插值 } else { // 插值輔助點1 usDiffTemp1 = ((SINT16)sqrt(FT_Diff1[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff1[((i + 2) / 2) * (ucRxNum / 2) + ((j + 2) / 2)]) + (SINT16)sqrt(FT_Diff1[(i + 2) / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff1[i / 2 * (ucRxNum / 2) + ((j + 2) / 2)])) / 2; // 插值輔助點2 usDiffTemp2 = ((SINT16)sqrt(FT_Diff1[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff1[((i - 2) / 2) * (ucRxNum / 2) + ((j + 2) / 2)]) + (SINT16)sqrt(FT_Diff1[(i - 2) / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff1[i / 2 * (ucRxNum / 2) + ((j + 2) / 2)])) / 2; // 通過輔助點內插 usDiffTemp1 = (SINT16)sqrt(usDiffTemp1 * usDiffTemp2); usDiffTemp2 = (SINT16)sqrt(FT_Diff1[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff1[i / 2 * (ucRxNum / 2) + ((j + 2) / 2)]); // 矩陣2同樣方法內插 usDiffTemp1 = ((SINT16)sqrt(FT_Diff2[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[((i + 2) / 2) * (ucRxNum / 2) + ((j + 2) / 2)]) + (SINT16)sqrt(FT_Diff2[(i + 2) / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[i / 2 * (ucRxNum / 2) + ((j + 2) / 2)])) / 2; usDiffTemp2 = ((SINT16)sqrt(FT_Diff2[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[((i - 2) / 2) * (ucRxNum / 2) + ((j + 2) / 2)]) + (SINT16)sqrt(FT_Diff2[(i - 2) / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[i / 2 * (ucRxNum / 2) + ((j + 2) / 2)])) / 2; usDiffTemp1 = (SINT16)sqrt(usDiffTemp1 * usDiffTemp2); usDiffTemp2 = (SINT16)sqrt(FT_Diff2[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[i / 2 * (ucRxNum / 2) + ((j + 2) / 2)]); }}// 類型3if(i%2 != 0 && j%2 == 0){ // 邊距8不考慮插值 if(i > (ucTxNum - 2) || j == 0 || j > (ucRxNumVA - 2)) { // 邊距8不考慮插值 } else { // 插值輔助點1 usDiffTemp1 = ((SINT16)sqrt(FT_Diff1[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff1[(i + 2) / 2 * (ucRxNum / 2) + (j - 2) / 2]) + (SINT16)sqrt(FT_Diff1[(i + 2) / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff1[i / 2 * (ucRxNum / 2) + ((j - 2) / 2)])) / 2; // 插值輔助點2 usDiffTemp2 = ((SINT16)sqrt(FT_Diff1[i / 2 * (ucRxNum/ 2) + (j / 2)] * FT_Diff1[(i + 2) / 2 * (ucRxNum / 2) + (j + 2) / 2]) + (SINT16)sqrt(FT_Diff1[(i + 2) / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff1[i / 2 * (ucRxNum / 2) + ((j + 2) / 2)])) / 2; // 通過輔助點內插 usDiffTemp1 = (SINT16)sqrt(usDiffTemp1 * usDiffTemp2); usDiffTemp2 = (SINT16)sqrt(FT_Diff1[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff1[(i + 2) / 2 * (ucRxNum / 2) + (j / 2)]); // 矩陣2同樣方法內插 usDiffTemp1 = ((SINT16)sqrt(FT_Diff2[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[(i + 2) / 2 * (ucRxNum / 2) + (j - 2) / 2]) + (SINT16)sqrt(FT_Diff2[(i + 2) / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[i / 2 * (ucRxNum / 2) + ((j - 2) / 2)])) / 2; usDiffTemp2 = ((SINT16)sqrt(FT_Diff2[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[(i + 2) / 2 * (ucRxNum / 2) + (j + 2) / 2]) + (SINT16)sqrt(FT_Diff2[(i + 2) / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[i / 2 * (ucRxNum / 2) + ((j + 2) / 2)])) / 2; usDiffTemp1 = (SINT16)sqrt(usDiffTemp1 * usDiffTemp2); usDiffTemp2 = (SINT16)sqrt(FT_Diff2[i / 2 * (ucRxNum / 2) + (j / 2)] * FT_Diff2[(i + 2) / 2 * (ucRxNum / 2) + (j / 2)]); }}// 類型4if(i%2 != 0 && j%2 != 0){ // 邊距8不考慮插值 if(i > (ucTxNum - 2) || j > (ucRxNumVA - 2)) { // 邊距8不考慮插值 } else { // 通過有效點內插 usDiffTemp1 = ((SINT16)sqrt(FT_Diff1[i / 2 * (ucRxNum / 2) + j / 2] * FT_Diff1[(i + 2) / 2 * (ucRxNum / 2) + (j + 2) / 2]) + (SINT16)sqrt(FT_Diff1[(i + 2) / 2 * (ucRxNum / 2) + j / 2] * FT_Diff1[i / 2 * (ucRxNum / 2) + (j + 2) / 2])) / 2; // 通過有效點內插 usDiffTemp2 = ((SINT16)sqrt(FT_Diff2[i / 2 * (ucRxNum / 2) + j / 2] * FT_Diff2[(i + 2) / 2 * (ucRxNum / 2) + (j + 2) / 2]) + (SINT16)sqrt(FT_Diff2[(i + 2) / 2 * (ucRxNum / 2) + j / 2] * FT_Diff2[i / 2 * (ucRxNum / 2) + (j + 2) / 2])) / 2; }}

推薦閱讀：