前言

自從上世紀(jì)50年代以來，計(jì)算機(jī)視覺技術(shù)一直是學(xué)術(shù)界和工業(yè)界的熱點(diǎn)問題，尤其是近十年來，隨著新一輪人工智能技術(shù)的蓬勃發(fā)展，得到AI加持的計(jì)算機(jī)視覺技術(shù)快速走進(jìn)了人們的日常生活，不斷重塑著人們的生活方式。

顧名思義，計(jì)算機(jī)視覺是指讓計(jì)算機(jī)能夠理解、分析和處理圖像和視頻的技術(shù)，通過研究高效的計(jì)算機(jī)算法和系統(tǒng)，使計(jì)算機(jī)獲得像人類一樣理解和使用視覺信息的能力。在這個(gè)過程中，人們需要綜合運(yùn)用計(jì)算機(jī)科學(xué)、數(shù)學(xué)、物理學(xué)、生理學(xué)、心理學(xué)等跨學(xué)科的知識(shí)才能更好地解決一些實(shí)際問題。

隨著視覺應(yīng)用場(chǎng)景的不斷擴(kuò)大，越來越多的人已經(jīng)認(rèn)識(shí)到，相機(jī)成像技術(shù)處在一個(gè)十分關(guān)鍵的地位。在討論視覺算法的效能時(shí)，人們經(jīng)常會(huì)引用一個(gè)經(jīng)典的英語(yǔ)表述“garbage in，garbage out”，這是在說，如果相機(jī)采集的基礎(chǔ)圖像本身就沒有有用的信息，那么再?gòu)?fù)雜精妙的算法也不可能給出有用的結(jié)果。所以，作為愛芯元智的企業(yè)目標(biāo)之一，就是不斷研發(fā)最前沿的AI-ISP技術(shù)，期望為視覺算法提供優(yōu)質(zhì)的輸入圖像，使上層應(yīng)用能夠在各種高挑戰(zhàn)場(chǎng)景下保持穩(wěn)定的性能輸出。

說到圖像質(zhì)量，這其實(shí)是可以討論一千零一夜的話題，不過總的來說，人們會(huì)從分辨率、解像力、幀率、延遲、色彩還原、細(xì)節(jié)還原、噪聲、動(dòng)態(tài)范圍、偽像、拖影、畸變等多個(gè)不同的維度分別對(duì)圖像質(zhì)量進(jìn)行量化分析，最終形成一個(gè)總的評(píng)價(jià)。在一類典型的應(yīng)用場(chǎng)景中，視覺算法的核心任務(wù)是精確測(cè)量物體的形狀和距離，測(cè)得的參數(shù)將作為一系列重要決策的依據(jù)，有些決策可能關(guān)乎企業(yè)盈虧，而另一些決策則可能關(guān)乎人身安全。此時(shí)，輸入圖像的幾何畸變常常是影響系統(tǒng)性能的關(guān)鍵因素之一，為此人們發(fā)展出了一系列方法對(duì)圖像畸變進(jìn)行分析和處理。

在這篇文章里，我們想分享的是人們?cè)谔幚韴D像幾何畸變問題時(shí)常用的思想方法和數(shù)學(xué)工具，希望相關(guān)的討論能夠?qū)Σ糠肿x者有所啟發(fā)。本文所討論的部分技術(shù)已經(jīng)集成在愛芯元智AX620系列芯片的標(biāo)準(zhǔn)軟件發(fā)布包中，可以幫助智慧城市、機(jī)器視覺等行業(yè)的用戶對(duì)廣角鏡頭進(jìn)行畸變校正。另外，在智慧城市和輔助駕駛領(lǐng)域常見的雙攝、多攝拼接產(chǎn)品形態(tài)，其核心任務(wù)是對(duì)多個(gè)圖像進(jìn)行像素級(jí)校正和融合，所以對(duì)鏡頭參數(shù)標(biāo)定和畸變校正算法提出了更高的要求。在研究更高精度和收斂性的算法過程中，本文討論的方法和工具也會(huì)有所裨益。關(guān)于圖像拼接融合的技術(shù)架構(gòu)已經(jīng)在AX620和AX650系列芯片中得到了支持，關(guān)鍵算法和工具還在不斷優(yōu)化完善中，歡迎感興趣的朋友關(guān)注。

1. 數(shù)學(xué)建模

與很多其它學(xué)科一樣，計(jì)算機(jī)視覺從業(yè)者在分析相機(jī)成像問題時(shí)，會(huì)首先建立一個(gè)關(guān)于相機(jī)系統(tǒng)成像原理的數(shù)學(xué)模型，簡(jiǎn)稱相機(jī)模型（camera model）。一個(gè)最為基本且眾所周知的相機(jī)模型是小孔成像模型（pinhole model），其原理如圖1所示。

圖1 小孔成像模型（pinhole camera model）

這個(gè)模型假定光是沿直線傳播的，模型定義了一個(gè)特殊的點(diǎn)即光心，一個(gè)焦距參數(shù)，也就是成像面與光心之間的距離。通常會(huì)以光心為原點(diǎn)建立相機(jī)坐標(biāo)系，并規(guī)定物方空間的Z坐標(biāo)為正，像方空間的Z坐標(biāo)為負(fù)。為了計(jì)算方便，模型還定義了虛擬成像面，虛擬成像面上的像是正立的，可以避免負(fù)號(hào)。

小孔成像模型的優(yōu)點(diǎn)是足夠簡(jiǎn)單，只需要初中數(shù)學(xué)的相似三角形知識(shí)就可以計(jì)算像的大小。當(dāng)相機(jī)的視場(chǎng)角比較小時(shí)，小孔成像模型是相機(jī)成像過程的一個(gè)良好近似。但是當(dāng)相機(jī)使用廣角鏡頭，或者應(yīng)用場(chǎng)景對(duì)模型精度要求很高時(shí)，小孔成像模型就不夠準(zhǔn)確了，此時(shí)必須在模型中引入成像畸變模型，使模型與實(shí)際情況更加符合。

在工程實(shí)踐中，人們經(jīng)常需要用具體的數(shù)字來定義什么是大，什么是小，什么是成功，什么是失敗。在評(píng)價(jià)一個(gè)成像模型的精度時(shí)，人們會(huì)使用重投影誤差（reprojection error）指標(biāo)作為參考。這個(gè)指標(biāo)的含義是，假設(shè)物方空間中存在一個(gè)3D坐標(biāo)點(diǎn)P，它被相機(jī)成像后，實(shí)際從圖像上提取到的2D坐標(biāo)是p，此值包含了一切設(shè)計(jì)因素和隨機(jī)因素導(dǎo)致的誤差，如果再把P值代入相機(jī)成像模型進(jìn)行計(jì)算，模型預(yù)測(cè)的2D坐標(biāo)是，我們把p與之間的像素距離定義為P點(diǎn)的重投影誤差。圖2顯示了重投影誤差的基本思想。

圖2 重投影誤差（reprojection error）

顯然，一張圖像中會(huì)包含成千上萬(wàn)個(gè)與P點(diǎn)類似的點(diǎn)，如果只用一個(gè)點(diǎn)的誤差值來代表整張圖像那明顯是不太合適的。但是我們可以從一張圖像中篩選出一批有代表性的點(diǎn)，把這些點(diǎn)的重投影誤差的某種平均值作為成像模型的重投影誤差，這樣就比較公允了。進(jìn)一步地，如果算法的輸入數(shù)據(jù)是一個(gè)視頻序列，那么只用一張圖像的誤差來代表整個(gè)序列也同樣是不太合適的，此時(shí)需要把每一張圖像的誤差都計(jì)算出來，然后考察整個(gè)序列的平均值和方差情況，這樣得到的數(shù)據(jù)就會(huì)更可信一些。圖3是使用Matlab工具繪制的一個(gè)圖像序列的重投影誤差，從圖中可以看到，重投影誤差的序列平均值是0.18個(gè)像素，對(duì)于很多應(yīng)用而言這是一個(gè)相當(dāng)優(yōu)秀的值。另外，重投影誤差在圖像之間會(huì)存在一定幅度的波動(dòng)，而波動(dòng)的大小也反應(yīng)了成像模型的穩(wěn)定性。

圖3 圖像序列的重投影誤差舉例

當(dāng)我們掌握了重投影誤差這種評(píng)價(jià)方法之后，假設(shè)有一個(gè)新的成像模型擺在面前，我們就可以客觀地判斷這個(gè)模型是好是壞，是有競(jìng)爭(zhēng)力還是缺乏競(jìng)爭(zhēng)力。在很多時(shí)候，這個(gè)能力可以幫助我們規(guī)避錯(cuò)誤決策的風(fēng)險(xiǎn)，讓我們總是站在一個(gè)相對(duì)安全的位置上。

回到剛才的話題，對(duì)于一個(gè)特定的相機(jī)系統(tǒng)，我們把偏離小孔成像模型的行為定義為一種畸變（distortion），當(dāng)我們發(fā)現(xiàn)使用小孔成像模型得到的重投影誤差超出可接受標(biāo)準(zhǔn)時(shí)，此時(shí)就需要設(shè)計(jì)一個(gè)合適的數(shù)學(xué)模型來定量地描述畸變，然后根據(jù)畸變模型對(duì)原圖像進(jìn)行校正，使校正圖像的重投影誤差收斂到可接受水平。

通過考察畸變行為我們能夠發(fā)現(xiàn)一些基本事實(shí)，首先，大部分的成像畸變都可以歸因于鏡頭的設(shè)計(jì)、制造、裝配等環(huán)節(jié)，因此通常可以用鏡頭畸變來指代成像畸變；其次，鏡頭畸變的特性和觀察尺度有關(guān)，在宏觀尺度上觀察時(shí)，廣角鏡頭和魚眼鏡頭通常表現(xiàn)為桶形畸變，長(zhǎng)焦鏡頭經(jīng)常表現(xiàn)為枕形畸變，這些行為是由鏡頭的光學(xué)設(shè)計(jì)決定的，每個(gè)鏡頭都應(yīng)該具有相同的行為，而在微觀尺度上觀察時(shí)，鏡頭的某些局部區(qū)域會(huì)表現(xiàn)出相對(duì)宏觀規(guī)律的偏離，這往往是加工制造環(huán)節(jié)中一些不可控因素導(dǎo)致的，每個(gè)鏡頭的具體行為不一定相同；最后，不同生產(chǎn)批次的鏡頭，由于材料、加工、裝配等因素的變化，也會(huì)導(dǎo)致批次與批次之間呈現(xiàn)出可觀測(cè)的行為差異。

基于以上觀察我們認(rèn)識(shí)到，鏡頭畸變模型可能需要分成兩類：如果應(yīng)用場(chǎng)景對(duì)模型的精度要求不高，比如圖像主要是服務(wù)于人眼觀察，那么畸變模型只要能夠描述鏡頭的宏觀設(shè)計(jì)行為就可以了；如果應(yīng)用場(chǎng)景對(duì)模型的精度要求很高，比如基于圖像進(jìn)行測(cè)量的場(chǎng)景，那么畸變模型必須能夠描述鏡頭的局部特性，此類模型所使用的參數(shù)數(shù)量一定遠(yuǎn)遠(yuǎn)多于前者，獲取這些參數(shù)的過程也會(huì)更復(fù)雜。

對(duì)于描述鏡頭宏觀設(shè)計(jì)行為的畸變模型，我們不妨對(duì)問題做一些合理的簡(jiǎn)化，首先假設(shè)鏡頭的行為是圓對(duì)稱的，只需要考察任一直徑方向上的成像規(guī)律，其次假設(shè)存在一條理想的光滑曲線可以精確地描述鏡頭的成像行為。基于這兩條假設(shè)，我們可以建立理想像點(diǎn)（基于小孔成像模型）與實(shí)際像點(diǎn)的函數(shù)關(guān)系r'=f(r)，在某些情況下，我們可以把r作為理想像點(diǎn)，r'作為實(shí)際像點(diǎn)，在另一些情況下，也可以把r作為實(shí)際像點(diǎn)，r'作為理想像點(diǎn)。在r=0附近對(duì)光滑曲線f(r)做泰勒展開，并忽略階數(shù)高于N的余項(xiàng)，即可得到關(guān)于鏡頭實(shí)際成像公式的多項(xiàng)式模型：

公式（1）

其中，r為理想像點(diǎn)（或?qū)嶋H像點(diǎn)）到鏡頭光心的距離，p1,?,pN是待定參數(shù)，需要通過實(shí)驗(yàn)標(biāo)定的方法獲得具體的值。

在實(shí)踐中，可以對(duì)(1)做進(jìn)一步簡(jiǎn)化，只使用其中的奇數(shù)次多項(xiàng)式或偶數(shù)次多項(xiàng)式，即：

公式（2）

公式（3）

2. 優(yōu)化問題

在前一節(jié)的討論中，我們已經(jīng)能夠建立一個(gè)關(guān)于相機(jī)成像的數(shù)學(xué)模型，其中包含了鏡頭畸變的模型，但是該模型中存在一組待定系數(shù)p1,?,pN目前還不知道具體值，需要找到最優(yōu)解。在數(shù)學(xué)上，這是一類典型的優(yōu)化問題，已經(jīng)有了很成熟的解決方案。實(shí)際上，早在1809年就出現(xiàn)了Gauss-Newton法來處理非線性最小二乘優(yōu)化問題，但是該方法在迭代過程中要求目標(biāo)函數(shù)的雅可比矩陣必須列滿秩，這一條件限制了它的應(yīng)用。Levenberg于1944年提出了一個(gè)改進(jìn)的方法，然后Marquardt在1963年進(jìn)一步發(fā)展出了Levenberg-Marquardt方法，目前L-M方法已經(jīng)成為非線性優(yōu)化的首選方法。關(guān)于L-M方法的實(shí)現(xiàn)原理已有很多文獻(xiàn)給出了很好的描述，在OpenCV、Python、Matlab等工具中也提供了相應(yīng)的庫(kù)可以直接調(diào)用，對(duì)此感興趣的讀者可以深入研究。

回到我們的具體問題，使用L-M方法求解鏡頭畸變模型的待定系數(shù)的過程是，首先要選定一組樣本G作為訓(xùn)練集，假設(shè)G中包含m個(gè)樣本點(diǎn)。然后，定義變量x=(p1,?,pN)T，表示模型的待定參數(shù)；定義函數(shù)，表示待定參數(shù)為x時(shí)，樣本G對(duì)應(yīng)的模型預(yù)測(cè)值，該函數(shù)返回一個(gè)包含m個(gè)2D點(diǎn)坐標(biāo)的向量；定義函數(shù)S(x)，表示從圖像中提取的關(guān)于樣本G的實(shí)測(cè)值，同理，該函數(shù)返回一個(gè)包含m個(gè)2D點(diǎn)坐標(biāo)的向量；

定義重投影誤差函數(shù)：

定義優(yōu)化的目標(biāo)函數(shù)：

當(dāng)x取最優(yōu)解x*時(shí)，F(xiàn)(x*)取最小值。

此外，我們還需要手動(dòng)提供一個(gè)x的初值x0，作為迭代優(yōu)化的起始點(diǎn)。

在每一次迭代過程中，L-M算法會(huì)計(jì)算函數(shù)f(x)的雅可比矩陣J，然后根據(jù)公式(4)計(jì)算下一個(gè)x點(diǎn)距離當(dāng)前點(diǎn)的步長(zhǎng)?xlm，

公式（4）

其中，μ≥0是一個(gè)阻尼因子，用于動(dòng)態(tài)地調(diào)整?xlm的步長(zhǎng)。

為了方便討論，我們先給出LM算法的偽代碼，然后簡(jiǎn)要地討論一下其中涉及到的數(shù)學(xué)工具，期望起到拋磚引玉的作用。

觀察上述偽代碼可知，LM算法的關(guān)鍵在于如何求解線性方程組(A+μI)h=g，其中h是由待定參數(shù)構(gòu)成的向量，系數(shù)矩陣A=J(x)TJ(x)是一個(gè)n階方陣，它的元素是對(duì)傳感器圖像運(yùn)行某種特征提取算法得到的，比如圖4所示的例子，是使用Matlab軟件從相機(jī)拍攝的棋盤格圖像中提取角點(diǎn)的位置坐標(biāo)。

圖4 使用Matlab工具提取棋盤格特征點(diǎn)

由于圖卡不平整、曝光不均勻、軟件四舍五入等隨機(jī)因素的存在，矩陣A的元素值包含了多種誤差的貢獻(xiàn)，因此本質(zhì)上并不存在一個(gè)“絕對(duì)正確的h”，我們只能希望得到一個(gè)“較好的h”。一般而言，參與計(jì)算的樣本越多，離群樣本的貢獻(xiàn)就越小，得到的結(jié)果就越優(yōu)。接下來的問題是，如何求解h向量才能獲得最佳的數(shù)值精度和穩(wěn)定性，因?yàn)椴煌姆椒ㄔ谔匦陨鲜怯幸欢▍^(qū)別的。在實(shí)踐中，我們經(jīng)常會(huì)使用三種方法來求解線性方程組。

我們假定矩陣A的維度為nxn，h為n維列向量，包含n個(gè)待定的參數(shù)，I為n維單位矩陣，J或J(x)為mxn的雅克比矩陣。為了求解n個(gè)未知數(shù)，數(shù)學(xué)上必須有m≥n，否則這是一個(gè)欠定問題，h有無窮多解。當(dāng)然，我們已經(jīng)知道樣本點(diǎn)是越多越好，所以m≥n本來也是在計(jì)劃之內(nèi)的。

進(jìn)一步地，我們記y=f(x)為m維列向量，記為(m+n)維列向量。另外，我們定義是一個(gè)分塊矩陣，其維度為 m’xn，m’ = (m+n)。借助以上定義，我們可以把兩個(gè)矩陣相加等價(jià)地表示成兩個(gè)矩陣相乘，即：

公式（5）

其中，上標(biāo)T代表矩陣的轉(zhuǎn)置操作。

我們馬上會(huì)發(fā)現(xiàn)，使用兩個(gè)矩陣相乘的表示方法會(huì)給我們帶來很多便利。

方法一，Cholesky分解

對(duì)矩陣(A+μI)進(jìn)行Cholesky分解，即：A+μI=RTR，其中I為單位陣，R為上三角矩陣。注意到矩陣A= J(x)TJ(x)+μI為正定矩陣，則R的對(duì)角元素為正。我們得到 RTRh=g，則利用兩次“回代”即可求出h。

Cholesky分解計(jì)算速度很快，所以使用較為廣泛，但該方法存在較大的計(jì)算誤差。深入研究發(fā)現(xiàn)，Cholesky分解的誤差正比于矩陣(A+μI)的條件數(shù)+，而它又等于矩陣Jμ條件數(shù)的平方。因此為減少數(shù)值計(jì)算的誤差，我們的計(jì)算對(duì)象應(yīng)該是Jμ而不是(A+μI)。

+矩陣的條件數(shù)

矩陣的條件數(shù)（condition number）是數(shù)值分析和線性代數(shù)中的一個(gè)重要概念，任意矩陣B的條件數(shù)為σmax /σmin，其中，σmax表示B的最大奇異值，σmin表示B的最小奇異值。條件數(shù)刻畫了矩陣令向量發(fā)生形變的能力。條件數(shù)越大，向量在變換后可能的變化量就越大，而大的條件數(shù)會(huì)放大數(shù)據(jù)中存在的誤差。關(guān)于條件數(shù)的介紹可參考書籍：Matrix Analysis-Roger A. Horn -2nd Edition

方法二 QR分解

我們注意到，方程(A+μI)h=g=-JTy，這是一個(gè)線性最小二乘問題，即min||Jμh+y'||2的正規(guī)方程（關(guān)于h的梯度等于0的方程），則求解該方程等價(jià)于求對(duì)應(yīng)線性二乘問題的解。

我們對(duì)矩陣Jμ進(jìn)行QR分解得到：Jμ?P=QR=Q1R1，其中，Q是維度為m’xm的單位正交矩陣，R為m’xn維的上三角矩陣；Q1為Q的前n列，維度為m’xn；R1為R的前n行，維度為nxn；P為列變換基本矩陣，使得R的對(duì)角線元素的絕對(duì)值遞減(即|Rii|≥|Rjj|，i＜j)，其維度為nxn。

將該QR分解等式帶入線性最小二乘的式中化簡(jiǎn)可知，R1PTh=-Q1Ty'，然后利用回代法，即可得到h。

QR分解法的計(jì)算誤差正比于矩陣Jμ的條件數(shù)，而不是該矩陣條件數(shù)的平方，所以它比Cholesky分解更穩(wěn)定，但會(huì)增加計(jì)算量。當(dāng)優(yōu)化變量的維度較小時(shí)，QR分解法使用比較廣泛，實(shí)際上很多流行的最優(yōu)化庫(kù)（比如minpack和ceres）都是基于QR分解實(shí)現(xiàn)的。

方法三SVD分解

將矩陣Jμ進(jìn)行SVD分解，即：Jμ=USVT=U1S1VT，其中：U為m’xm’維的單位正交陣；S為m’xn維的對(duì)角陣；V為nxn維的單位正交陣；U1為U的前n列，維度為m’xn；S1為S的前n行，維度為nxn；將SVD分解等式帶入線性最小二乘的式中并化簡(jiǎn)，可知：

簡(jiǎn)單整理可得下方公式，即可計(jì)算出h。

SVD法的穩(wěn)定性很高，但其計(jì)算量也最大。

以上我們討論了求解LM算法核心方程組時(shí)三種常用方法的優(yōu)缺點(diǎn)。另外我們也注意到，以L-M算法為代表的非線性優(yōu)化方法并不能保證總是給出收斂的解，其收斂與否往往與我們手動(dòng)提供的初值質(zhì)量有關(guān)。在工程領(lǐng)域，一個(gè)不保證總是有效的方案其可用性是大打折扣的，必須要有其它機(jī)制作為失效后的補(bǔ)充。因此，一個(gè)十分重要的研究課題是尋找一些穩(wěn)定可靠的線性優(yōu)化方法，一方面為主算法提供一個(gè)良好的迭代初值，另一方面可以作為主算法失效時(shí)的備選方案。限于篇幅，本文對(duì)此就不做展開討論了，如果讀者喜歡，后續(xù)可以專題討論一下這個(gè)問題。

小結(jié)

在本文中，我們主要討論了計(jì)算機(jī)視覺領(lǐng)域經(jīng)常遇到的一個(gè)基礎(chǔ)性問題，即如何定量地評(píng)價(jià)圖像的幾何畸變，以及如何通過數(shù)學(xué)建模的方法為鏡頭建立成像模型，從而能夠定量地描述成像畸變并加以校正。當(dāng)鏡頭模型建立后，為了求解模型中的待定參數(shù)，我們還介紹了L-M方法，這是一種非常實(shí)用的最小二乘非線性優(yōu)化方法，在工程技術(shù)領(lǐng)域有著非常廣泛的應(yīng)用。雖然很多文獻(xiàn)和資料都會(huì)有涉及LM算法的應(yīng)用，但是討論LM算法實(shí)現(xiàn)細(xì)節(jié)的資料卻頗為少見，所以我們借本文簡(jiǎn)要討論了求解LM算法可以使用的數(shù)學(xué)工具以及工具背后的數(shù)理邏輯，希望我們所討論的這些思想方法和具體工具能對(duì)讀者帶來一定的啟發(fā)，非常感謝廣大讀者對(duì)愛芯元智的關(guān)注和支持，我們下期再會(huì)！

審核編輯：湯梓紅