雙目立體視差圖計算
clear;
%% 載入2張立體影象
left = imread('left.png');
right = imread('right.png');
sizeI = size(left);
% 顯示覆合影象
zero = zeros(sizeI(1), sizeI(2));
channelRed = left(:,:,1);
channelBlue = right(:,:,3);
composite = cat(3, channelRed, zero, channelBlue);
figure(1);
subplot(2,2,1);
imshow(left);
axis image;
title('Left Image');
subplot(2,2,2);
imshow(right);
axis image;
title('Right Image');
subplot(2,2,3);
imshow(composite);
axis image;
title('Composite Image');
%% 基本的塊匹配
% 通過估計子畫素的塊匹配計算視差
disp('執行基本的塊匹配~');
% 啟動定時器
tic();
% 平均3個顏色通道值將RGB影象轉換為灰度影象
leftI = mean(left, 3);
rightI = mean(right, 3);
% SHD
bitsUint8 = 8;
leftI = im2uint8(leftI./255.0);
rightI = im2uint8(rightI./255.0);
% DbasicSubpixel將儲存塊匹配的結果,元素值為單精度32位浮點數
DbasicSubpixel = zeros(size(leftI), 'single');
% 獲得影象大小
[imgHeight, imgWidth] = size(leftI);
% 視差範圍定義離第1幅影象中的塊位置多少畫素遠來搜尋其它影象中的匹配塊。對於大小為450x375的影象,視差範圍為50是合適的
disparityRange = 50;
% 定義塊匹配的塊大小
halfBlockSize = 5;
blockSize = 2 * halfBlockSize + 1;
% 對於影象中的每行(m)畫素
for (m = 1 : imgHeight)
% 為模板和塊設定最小/最大塊邊界
% 比如:第1行,minr = 1 且 maxr = 4
minr = max(1, m - halfBlockSize);
maxr = min(imgHeight, m + halfBlockSize);
% 對於影象中的每列(n)畫素
for (n = 1 : imgWidth)
% 為模板設定最小/最大邊界
% 比如:第1列,minc = 1 且 maxc = 4
minc = max(1, n - halfBlockSize);
maxc = min(imgWidth, n + halfBlockSize);
% 將模板位置定義為搜尋邊界,限制搜尋使其不會超出影象邊界
% 'mind'為能夠搜尋至左邊的最大畫素數;'maxd'為能夠搜尋至右邊的最大畫素數
% 這裡僅需要向右搜尋,所以mind為0
% 對於要求雙向搜尋的影象,設定mind為max(-disparityRange, 1 - minc)
mind = 0;
maxd = min(disparityRange, imgWidth - maxc);
% 選擇右邊的影象塊用作模板
template = rightI(minr:maxr, minc:maxc);
% 獲得本次搜尋的影象塊數
numBlocks = maxd - mind + 1;
% 建立向量來儲存塊偏差
blockDiffs = zeros(numBlocks, 1);
% 計算模板和每塊的偏差
for (i = mind : maxd)
%選擇左邊影象距離為'i'處的塊
block = leftI(minr:maxr, (minc + i):(maxc + i));
% 計算塊的基於1的索引放進'blockDiffs'向量
blockIndex = i - mind + 1;
%{
% NCC(Normalized Cross Correlation)
ncc = 0;
nccNumerator = 0;
nccDenominator = 0;
nccDenominatorRightWindow = 0;
nccDenominatorLeftWindow = 0;
%}
% 計算模板和塊間差的絕對值的和(SAD)作為結果
for (j = minr : maxr)
for (k = minc : maxc)
%{
% SAD(Sum of Absolute Differences)
blockDiff = abs(rightI(j, k) - leftI(j, k + i));
% SSD(Sum of Squared Differences)
% blockDiff = (rightI(j, k) - leftI(j, k + i)) * (rightI(j, k) - leftI(j, k + i));
blockDiffs(blockIndex, 1) = blockDiffs(blockIndex, 1) + blockDiff;
%}
%{
% NCC
nccNumerator = nccNumerator + (rightI(j, k) * leftI(j, k + i));
nccDenominatorLeftWindow = nccDenominatorLeftWindow + (leftI(j, k + i) * leftI(j, k + i));
nccDenominatorRightWindow = nccDenominatorRightWindow + (rightI(j, k) * rightI(j, k));
%}
end
end
%{
% SAD
blockDiffs(blockIndex, 1) = sum(sum(abs(template - block)));
%}
%{
% NCC
nccDenominator = sqrt(nccDenominatorRightWindow * nccDenominatorLeftWindow);
ncc = nccNumerator / nccDenominator;
blockDiffs(blockIndex, 1) = ncc;
%}
% SHD(Sum of Hamming Distances)
blockXOR = bitxor(template, block);
distance = uint8(zeros(maxr - minr + 1, maxc - minc + 1));
for (k = 1 : bitsUint8)
distance = distance + bitget(blockXOR, k);
end
blockDiffs(blockIndex, 1) = sum(sum(distance));
end
% SAD值排序找到最近匹配(最小偏差),這裡僅需要索引列表
% SAD/SSD/SHD
[temp, sortedIndeces] = sort(blockDiffs, 'ascend');
%{
% NCC
[temp, sortedIndeces] = sort(blockDiffs, 'descend');
%}
% 獲得最近匹配塊的基於1的索引
bestMatchIndex = sortedIndeces(1, 1);
% 將該塊基於1的索引恢復為偏移量
% 這是基本的塊匹配產生的最後的視差結果
d = bestMatchIndex + mind - 1;
%{
% 通過插入計算視差的子畫素估計
% 子畫素估計要求用左右邊的塊, 所以如果最佳匹配塊在搜尋窗的邊緣則忽略估計
if ((bestMatchIndex == 1) || (bestMatchIndex == numBlocks))
% 忽略子畫素估計並儲存初始視差值
DbasicSubpixel(m, n) = d;
else
% 取最近匹配塊(C2)的SAD值和最近的鄰居(C1和C3)
C1 = blockDiffs(bestMatchIndex - 1);
C2 = blockDiffs(bestMatchIndex);
C3 = blockDiffs(bestMatchIndex + 1);
% 調整視差:估計最佳匹配位置的子畫素位置
DbasicSubpixel(m, n) = d - (0.5 * (C3 - C1) / (C1 - (2 * C2) + C3));
end
%}
DbasicSubpixel(m, n) = d;
end
% 每10行更新過程
if (mod(m, 10) == 0)
fprintf('影象行:%d / %d (%.0f%%)\n', m, imgHeight, (m / imgHeight) * 100);
end
end
% 顯示計算時間
elapsed = toc();
fprintf('計算視差圖花費 %.2f min.\n', elapsed / 60.0);
%% 顯示視差圖
fprintf('顯示視差圖~\n');
% 切換到影象4
subplot(2,2,4);
% 第2個引數為空矩陣,從而告訴imshow用資料的最小/最大值,並且對映資料範圍來顯示顏色
imshow(DbasicSubpixel, []);
% 去掉顏色圖會顯示灰度視差圖
colormap('jet');
colorbar;
% 指定視差圖的最小/最大值
%caxis([0 disparityRange]);
%設定顯示的標題
title(strcat('Basic block matching, Sub-px acc., Search left, Block size = ', num2str(blockSize)));
%% 誤匹配畫素的百分比
thresh = 1;
computedDisparityMap = double(DbasicSubpixel);
groundTruthDisparityMap = double(imread('disp6.png'));
scale = max(max(groundTruthDisparityMap)) / max(max(computedDisparityMap));
computedDisparityMap = computedDisparityMap * scale;
numPixels=0;
perBADMatch=0.0;
tic;
for (i = 1 + halfBlockSize : imgHeight - halfBlockSize)
for(j = 1 + halfBlockSize : imgWidth - halfBlockSize - disparityRange)
if(groundTruthDisparityMap(i,j) ~= 0)
if(abs(computedDisparityMap(i,j) - groundTruthDisparityMap(i,j)) > thresh)
perBADMatch = perBADMatch + 1;
end
numPixels = numPixels + 1;
end
end
end
perBADMatch = perBADMatch / numPixels
timeTaken = toc;
fprintf('計算誤匹配畫素的百分比花費 %.2f min.\n', elapsed / 60.0);