正射影像

攝影測量 的第二個作業，不只幾何改正， 把傾斜攝影的照片作成正射影像。 我就用上學期的 code， 加上上次作業寫的影像取樣， 組合成影像正射化的程式。

上學期的 bottom up 共線式

上學期 航遙測概論 最後產出是一個 bottom up 共線式， 輸入內外方位參數，能把物空間座標轉到相空間； 或把相空間座標和高程，轉成物空間； 並作透鏡畸變改正。

上學期是用 coffee script 寫的， 寫太爛現在根本不忍直視，就不改了，直接引用。 反正結果是對的就好。

範例

// 引用
const ColinearityEquation = require('colinearityEquation')

// 下面是使用方法

// 外方位參數
const eop = {
    // 平移減少計算誤差
    // xo: 169570.443932,
    xo: 570.443932,
    // yo: 2544089.886767,
    yo: 89.886767,
    zo: 116.355578,
    // 角度轉弧度
    omega: 25.769170/180*Math.PI,
    phi: -0.478151/180*Math.PI,
    kappa:-4.904690/180*Math.PI
}

// 內方位參數
const xp = 0.00797976759482130898
const yp = -0.01632040141493363086
const c = 3.68128070687046982101
const iop = [xp, yp, -c]

// 透鏡畸變
iop.k = [
    -0.00261341937662040030,
    0.00015213497306022675,
    -0.00000127943044334018
]
iop.p = [0.00002392585197592054, 0.00001764654855396589]


// 輸入內外方位參數，創建共線式物件
const ce = new ColinearityEquation(eop, iop)

// tw97 (2544200, 169100, 19)
ce.groundToPhoto([200, 100, 19]) // 輸出像座標 [0.3, 0.1, 0]

// 像空間到物空間未知數不足，
// 要物空間 X Y Z 任一已知，這裡是高程 Z 已知為正高 19m。
ce.photoToGround([0.3, 0.1, 0], [NaN, NaN, 19]) // 輸出物座標 [200, 100, 19]

上次作業的影像處理程式

上次作業寫的程式： 輸入物空間到像空間的轉換函數， 和物空間轉換範圍、取樣距離； 就能把從該範圍自動取樣， 把每個地物座標用輸入的函數轉到像座標， 取該像座標的顏色填到物空間。

但上學期轉出來的結果是像空間座標， 要轉成 pixel 要乘上 pixel size 再平移， 所以要多一層轉換函數。

// 引用
const ReferenceImage = require('reference-image').ReferenceImage
const X = 0
const Y = 1
const Z = 2

// 輸入檔名和物空間到像空間的轉換函數，
// 構造一個 reference image 物件。
const imgp = new ReferenceImage('DJI_0339.jpg', groundToPhoto)

function groundToPhoto(ground) {
    // 補上高程座標
    ground[Z] = 19.7
    // 甪上面的共線式轉換物空間到像空間
    const photo = ce.groundToPhoto(ground)

    // 從像空間轉到 pixel 值
    const scale = 631.8955603017934
    const pixel = point.map((length) => length * scale)
    scalePoint[X] = 1995.5 - scalePoint[X]
    scalePoint[Y] = 1495.5 + scalePoint[Y]
    return scalePoint
}

imgp.then((img) => {
    // 轉換範圍，座標是平移 (-2544000, -169000, 0) 的 twd97
    const start = [100, 300]
    const end = [200, 500]
    img.gsd = 0.1 // 0.1m = 10cm 約要跑半小時
    return img.reference(start, end)
}).then((referenceImage) => {
    // 轉換完存檔
    referenceImage.write('DJI_0339-ortho.png')
})

方向相反

發現之前寫的程式有部份寫錯， 造成轉出來的結果上下相反。所以就改了。 上一次的作業沒有要比較實際座標， 這次才有，才發現數字不太對； 不然都只關注有沒有扭曲而已。

主要是像空間是原點在光軸中心， 而數位影像是用 pixel，原點在左上角， Y 座標方向是相反的。

包裝

之後因為有三張影像，就把共同的部份抽出來，寫成一個函數。 輸入檔名、內外方位、範圍、取樣區間，就能轉換好。 github 原始碼 。

const ortho = require('./ortho')
ortho.run(file, eop, iop, [start, end], gsd)

轉換一張 250x200，gsd 設 0.1m， 也就是 2500x2000=5000000 個點， 約耗時 30 分鐘；一秒轉換 2777 個點。

成果

轉換成果是扇型，因為是傾斜攝影。 有些樹或建築被壓到地上， 因為沒有 dtm，轉換時一律把高程設成 19.6， 差不多是球場的高度。

超出相片範圍我沒有處理， Jimp 的作法還是一樣內插，就用最近的點， 就變成從邊界延伸一條一條的了。

js [flk.fig(339, 'DJI_0339-ortho'), flk.fig(340, 'DJI_0340-ortho'), flk.fig(341, 'DJI_0341-ortho')]

可以看出各張影像有些微不同， 像是從不同角度傾斜拍攝。 可能是因為沒有用 dtm，高程用同一數值的結果。

分析

之後放到 QGIS 裡，點出各控制點在不同照片的座標， 再和真實座標比較：

js flk.gn('vector-clean-map')

圖中向量是從控制點真實座標， 延伸往正射後影像中量測的控制點位置， 為了明顯看出系統誤，再延伸五倍的結果。 而底圖是正射後的 340 影像。

三張正射影像，誤差方向不同， 也分別在北中南為誤差最小處。 東邊的球場，因為正高差較多，所以都偏較多。

系統來說， 預估誤差方向的作法是將共線式的向量投影到地面上， 高程值如果有誤，就會延向量投影方向移動。 移動方向和拍攝傾角有關，因傾角會影響共線式的向量。

js flk.gn('error-dtm')