Stereo Matching using Block SAD

這裡假設左右眼兩張影像已經經過 rectification，在 epipolar line 上限制一個小範圍的search range用以加速，簡單用 SAD 計算最小的左右視差作為 disparity。

使用 middlebury 的 dataset，勉強可以看出深度，但畢竟不到 100 行的小品 code 所以品質很差XD。

 

作為使用不同 window size 計算 SAD的比較，下圖左是用較大 window size 的結果，較為密實但粗略，小的 window size 則是會有較精確的但易破碎的結果；

因此才有了 coarse to fine 的演算法策略，先計算大的 window size 或低解析度的影樣匹配，然後再重複用小的 window size在鄰近 disparity 位置搜尋精確解。

但是 coarse to fine 也有無法解決的地方，太強烈的 smooth 特性反而失去細節，舉例來說，位於半島型的突出小東西或獨立的細小物體 ( 下圖綠框部分 ) ，在這些 window 中，真正要比較差異的區域其實是很少的，只在 window 正中間一小塊，越是 coarse 的範圍看到其他的部份也越多， 結果易忽略細小的目標物，而收斂成適合周邊範圍的 disparity；造成一種「穿越」或「填補」的現象。 coarse layer 定義下來的 disparity 若離目標物深度很遠，在 fine layer 相鄰深度的 refine 就會失敗 ( coarse layer 就消失的東西救不回來 )；許多論文靠著定義特殊的 label 或 segment 來區分這些小區塊做特別處理。Middlebury stereo Evaluation 上目前最好的效果算是這類方法。

後面附上 source code，使用 openCV library讀寫圖

 

#include <time.h>
#include <opencv2/opencv.hpp>

using namespace cv;

// parameter
int winSize     = 35;
int searchRange = 200;
float scale     = 0.25;

// auto parameter
int halfWinSize     = winSize /  2;
int halfSearchRange = searchRange / 2;

Mat1f getPatch(const Mat1f &img, int cx, int cy) {
	Range rangeY(cy - halfWinSize, cy + halfWinSize + 1);
	Range rangeX(cx - halfWinSize, cx + halfWinSize + 1);
	return img(rangeY, rangeX);
}

void stereoMatch(const Mat1f &imgSrc, const Mat1f &imgDst, Mat1f &disparity) {
	disparity = Mat1f::zeros( imgSrc.size() );

	for (int cy = halfWinSize; cy < imgSrc.rows-halfWinSize; ++cy) {
		for (int cxSrc = halfWinSize; cxSrc < imgSrc.cols-halfWinSize; ++cxSrc) {
			// left patch
			Mat1f patchL = getPatch( imgSrc, cxSrc, cy );

			// epipolar line search range
			int cxDstMin = max(cxSrc - halfSearchRange, halfWinSize);
			int cxDstMax = min(cxSrc + halfSearchRange, imgSrc.cols-halfWinSize);

			// find best match disparity
			float minSad        = FLT_MAX;
			float bestDisparity = 0;
			for (int cxDst = cxDstMin; cxDst < cxDstMax; ++cxDst) {
				// right patch
				Mat1f patchR = getPatch( imgDst, cxDst, cy );

				// patch diff: Sum of Absolute Difference
				Mat1f diff = abs( patchL - patchR );
				float sad  = sum(diff)[0];

				// update best SAD and disparity
				if ( sad < minSad ) {
					minSad        = sad;
					bestDisparity = abs( cxDst - cxSrc );
				}
			}

			disparity(cy, cxSrc) = bestDisparity;
		}
	}
}

int main(int argc, char *argv[]) {
	const char *fileNameL = argv[1];
	const char *fileNameR = argv[2];

	// read input image
	Mat1f imgL, imgR;
	imread(fileNameL, 0).convertTo(imgL, CV_32FC1);
	imread(fileNameR, 0).convertTo(imgR, CV_32FC1);

	// scaling down
	resize(imgL, imgL, Size(), scale, scale);
	resize(imgR, imgR, Size(), scale, scale);

	// compute disparity
	clock_t t0 = clock();
	Mat1f disparityL, disparityR;
	stereoMatch(imgL, imgR, disparityL); // disparity Left
	stereoMatch(imgR, imgL, disparityR); // disparity Right
	clock_t t1 = clock();
	printf("time: %f sec\n", (double)  (t1-t0) / CLOCKS_PER_SEC);

	// write file
	normalize(disparityL, disparityL, 0, 1, NORM_MINMAX);
	normalize(disparityR, disparityR, 0, 1, NORM_MINMAX);
	imshow("disparityL", disparityL);
	imshow("disparityR", disparityR);
	imwrite("disparityL.png", disparityL*255);
	imwrite("disparityR.png", disparityR*255);
	waitKey();

	return 0;
}