立體匹配之動態規劃DP
動態規劃通常應用於最優化問題的求解中，Baker、Ohta 等將動態規劃引入立體匹配中來獲取視差圖。動態規劃匹配的過程可以分為兩個階段，建立視差空間和動態規劃優化，將立體匹配問題轉化為視差空間內尋找最優路徑的問題。
密集匹配通常會充分利用影像間的核線約束條件，對立體像對進行核線糾正，這樣同名像點肯定位於對應的同名核線上，降低了匹配的難度。視差空間影像DSI（Disparity Space Image）是一種反映核線立體像對間同一掃描線視差關係的資料結構，將視差分析轉化為一種類似的譜分析。其數學表達為：
其中w 為影象的寬度， NaN 表示在該視差值為d 時不可能有對應的同名點，L((* )) 是對應點之間的相似性測度函式。在匹配核線影像時，逐行計算畫素的DSI，將每一行的DSI 依次疊加起來構成一個三維空間，稱之為視差空間，視差空間包含了每個畫素所有可能的視差取值。如下圖所示，圖中的x、y 代表影象的橫縱座標，d 表示視差搜尋範圍。
在建立好視差空間後，立體匹配問題就會轉化為在視差空間內尋找最佳路徑的問題，該路徑上的視差累積的能量最小，相似性最大。一般採用從影象左邊到右邊的順序，逐行進行。
實現過程如下：
//////minDeta、maxDeta為最小、大視差，P1、P2為平滑懲罰項
void Dynamic(Mat lImg,Mat rImg,int minDeta,int maxDeta,int P1,int P2)
{
if (lImg.empty() || rImg.empty() || (lImg.type() != CV_8UC3&&lImg.type() != CV_8UC1) || (rImg.type() != CV_8UC3&&rImg.type() != CV_8UC1)) return;
Mat lGrayImg = lImg.type == CV_8UC3 ? cvtColor(lImg, lGrayImg, COLOR_BGR2GRAY) : lImg;
Mat rGrayImg = rImg.type == CV_8UC3 ? cvtColor(rImg, rGrayImg, COLOR_BGR2GRAY) : rImg;
int detaRang=maxDeta-minDeta+1;///視差範圍
Mat dispMap(lGrayImg.size(),CV_8UC1,Scalar::all(0));//////左視差圖
for(int r=0;r<lGrayImg.rows;r++)
{
printf("動態規劃...%d%%\r",(int)(r*100.0/(lGrayImg.rows)));/////方便檢視進度
Mat dif(Size(detaRang, lGrayImg.cols),CV_32FC1,Scalar::all(0));////視差空間
Mat Emat(Size(detaRang, lGrayImg.cols),CV_32FC1,Scalar::all(0));/////能量空間
Mat Dmat(Size(detaRang, lGrayImg.cols),CV_8UC1,Scalar::all(0));/////路徑走向
for(int c=0;c<lGrayImg.cols;c++)
{
for(int d=0;d<detaRang;d++)
{
int rc=c-d-minDeta;///當前視差下左畫素對應對的右畫素座標
if(rc>=0)
{
dif.ptr<float>(c)[d] = SAD(lGrayImg, rGrayImg,Point2i(c, r),Point2i(c - d, r), 0);//////AD匹配代價
}
else
{
dif.ptr<float>(c)[d]=c<minDeta?0:dif.ptr<float>(c)[d-1];/////防止邊界超限處代價對整體路徑的影響
}
}
}
for(int c=1;c<lGrayImg.cols;c++)
{
for(int match=0;match<detaRang;match++)
{
float tmpmin=FLT_MAX;//////當前狀態下路徑代價
int pmin=NULL;////路徑走向
for(int pre=0;pre<=detaRang;pre++)
{
int Edata=dif.ptr<float>(c)[match];////資料項
int dd=abs(match-pre);
float Esmooth=dd>1?P2:(dd==0?0:P1);////視差等於1與大於1時候的平滑懲罰不同
float E=Edata+Esmooth+Emat.ptr<float>(c-1)[pre];/////路徑能量值
if (E<tmpmin)
{
tmpmin=E;
pmin=pre;
}
}
Emat.ptr<float>(c)[match]=tmpmin;
Dmat.ptr<uchar>(c)[match]=pmin;
}
}
int minID=0;
float min=Emat.ptr<float>(lGrayImg.cols-1)[0];
for(int i=0;i<detaRang;++i)
{
if(Emat.ptr<float>(lGrayImg.cols-1)[i]<min)
{
minID =i;
min=Emat.ptr<float>(lGrayImg.cols-1)[i];
}
}
//取得視差
for(int c= lGrayImg.cols-1;c>=0;--c)
{
int d=Dmat.ptr<uchar>(c)[minID];
minID =d;
dispMap.ptr<uchar>(r)[c]=d+minDeta;/////
}
}
}