%%%%%%%%%%%%%%belief propagation%%%%%%%%%%%%%%%%% %基于可信传播的立体匹配算法 %代价函数:AD %优化策略:max-product belief propagation,accelerated updating %参考论文:Comparison of Graph Cuts with Belief Propagation for Stereo, using %Identical MRF Parameters %参考论文:Stereo Matching Using Belief Propagation % %dlut %wzb clear; clc; T=4;%Plots平滑模型的梯度阈值 Truncation=20;%代价函数的截断阈值 P=4;%Plots模型参数 s=10;%Plos模型参数 D=50;%parameter for converting cost into compatibility exp(-C/D) iter=50;%belief propagation 算法的迭代次数 %公共参数 T = input('Input the smooth parameter of Plots model, T[4]: '); if isempty(T)%Plots平滑模型的梯度阈值 T = 4; end s = input('Input the smooth parameter of Plots model, s[10]: '); if isempty(s)%Plots模型参数 s = 10; end P = input('Input the smooth parameter of Plots model, P[4]: '); if isempty(P)%Plots模型参数 P = 4; end Truncation = input('Input the truncation value of cost, Truncation[20]: '); if isempty(Truncation)%代价函数的截断阈值 Truncation = 20; end iter = input('Input the iteration times, iter[50]: '); if isempty(iter) iter = 50; end %图像名 imagename = input('Input the image pairs name, imagename[tsukuba]: '); if isempty(imagename) imagename = 'tsukuba'; end reference = imread([imagename 'left.png']);%参考图像左图像 disp(['The reference image is: ' imagename 'left.png']); target = imread([imagename 'right.png']);%目标图像 disp(['The target image is: ' imagename 'right.png']); dmax = input('Input the disparity range(the max disparity)[16]:');%最大视差值,候选视差范围从0-dmax if isempty(dmax) dmax=16;%19 end sdmax='The disparity range is: 0--'; eval(['disp(''[' sdmax num2str(dmax) ']'')']); scalar = input('Input the disparity scalar when saved, scalar[16]:'); if isempty(scalar) scalar=16;%19 end disp('running....'); %%%%%%%%%%%%%belief propagation accelerated left image is the reference%%% rgbreference = reference; reference = rgb2gray(reference);reference = double(reference); target = rgb2gray(target);target = double(target); [nr,nc] = size(reference);%图像大小 disp = zeros(nr,nc); leftgradient = zeros(nr,nc); leftgradient(:,2:end) = abs(reference(:,2:end)-reference(:,1:nc-1));%左梯度 leftgradient = leftgradient>T; rightgradient = zeros(nr,nc); rightgradient(:,1:nc-1) = abs(reference(:,1:nc-1)-reference(:,2:end));%右梯度 rightgradient = rightgradient>T; upgradient = zeros(nr,nc); upgradient(2:end,:) = abs(reference(2:end,:)-reference(1:nr-1,:));%上梯度 upgradient = upgradient>T; downgradient = zeros(nr,nc); downgradient(1:nr-1,:) = abs(reference(1:nr-1,:)-reference(2:end,:));%下梯度 downgradient = downgradient>T; %prawcost = Ad(reference,target,dmax); prawcost=zeros(nr,nc,dmax+1); for d=0:dmax%d为循环变量 prawcost(:,d+1:nc,d+1)=abs(reference(:,d+1:nc)-target(:,1:nc-d)); for temp=1:d%图像边界处理 prawcost(:,temp,d+1)=prawcost(:,d+1,d+1); end end prawcost(prawcost>Truncation) = Truncation;%代价函数 prawmessage = exp(-prawcost/D); clear reference; clear target; leftmessage = ones(nr,nc,dmax+1);%从左邻域像素传递来的信息存储矩阵 rightmessage = leftmessage;%从右邻域像素传递来的信息存储矩阵 upmessage = leftmessage;%从上邻域像素传来的信息存储矩阵 downmessage = leftmessage;%从下邻域像素传来的信息存储矩阵 smooth1 = exp(-s/D);%平滑项 smooth2 = exp(-P*s/D); E = []; tic for iteriter = 1:iter %beleif propagation 迭代过程 %,由于参考图像为左图像,参考图像左边界像素在右图像中没有匹配点的概率很大,所以选择信息先从右向左传播的策略 tempupdownpraw = upmessage.*downmessage.*prawmessage;%上信息,下信息,源信息的乘积 for row = 2:nr-1%计算右信息,信息从右向左传播 for col = 2:nc-1 currentcol = nc-col+1;%当前像素列坐标 tempvector = rightmessage(row,currentcol+1,:).*tempupdownpraw(row,currentcol+1,:);%右信息,上信息,下信息,源信息的乘积 tempvector = tempvector(:); if rightgradient(row,currentcol) smooth = smooth1; else smooth = smooth2; end tempvector = tempvector*smooth; maxvectorelement = max(tempvector);%加入平滑项的最大传递信息 for d = 1:dmax+1 tempvectorelement = tempvector(d)/smooth; rightmessage(row,currentcol,d) = max(maxvectorelement,tempvectorelement); %去掉平滑项的最大传递信息 end rightmessage(row,currentcol,:) = rightmessage(row,currentcol,:)*(dmax+1)/sum(rightmessage(row,currentcol,:));%归一化 end end%计算右信息,信息从右向左传播 for row = 2:nr-1%计算左信息,信息从左向右传播 for col = 2:nc-1 tempvector = leftmessage(row,col-1,:).*tempupdownpraw(row,col-1,:);%左信息,上信息,下信息,源信息的乘积 tempvector = tempvector(:); if leftgradient(row,col) smooth = smooth1; else smooth = smooth2; end tempvector = tempvector*smooth; maxvectorelement = max(tempvector);%加入平滑项的最大传递信息 for d = 1:dmax+1 tempvectorelement = tempvector(d)/smooth; leftmessage(row,col,d) = max(maxvectorelement,tempvectorelement); %去掉平滑项的最大传递信息 end leftmessage(row,col,:) = leftmessage(row,col,:)*(dmax+1)/sum(leftmessage(row,col,:));%归一化 end end%计算左信息,信息从左向右传播 templeftrightpraw = leftmessage.*rightmessage.*prawmessage;%左信息,右信息,源信息的乘积 for col = 2:nc-1%计算上信息,信息从上向下传播 for row = 2:nr-1 tempvector = upmessage(row-1,col,:).*templeftrightpraw(row-1,col,:);%左信息,右信息,上信息,源信息的乘积 tempvector = tempvector(:); if upgradient(row,col) smooth = smooth1; else smooth = smooth2; end tempvector = tempvector*smooth; maxvectorelement = max(tempvector);%加入平滑项的最大传递信息 for d = 1:dmax+1 tempvectorelement = tempvector(d)/smooth; upmessage(row,col,d) = max(maxvectorelement,tempvectorelement); %去掉平滑项的最大传递信息 end upmessage(row,col,:) = upmessage(row,col,:)*(dmax+1)/sum(upmessage(row,col,:));%归一化 end end%计算上信息,信息从上向下传播 for col = 2:nc-1%计算下信息,信息从下向上传播 for row = 2:nr-1 currentrow = nr-row+1; tempvector = downmessage(currentrow+1,col,:).*templeftrightpraw(currentrow+1,col,:);%左信息,右信息,上信息,源信息的乘积 tempvector = tempvector(:); if downgradient(currentrow,col) smooth = smooth1; else smooth = smooth2; end tempvector = tempvector*smooth; maxvectorelement = max(tempvector);%加入平滑项的最大传递信息 for d = 1:dmax+1 tempvectorelement = tempvector(d)/smooth; downmessage(currentrow,col,d) = max(maxvectorelement,tempvectorelement); %去掉平滑项的最大传递信息 end downmessage(currentrow,col,:) = downmessage(currentrow,col,:)*(dmax+1)/sum(downmessage(currentrow,col,:));%归一化 end end%计算下信息,信息从下向上传播 belief = downmessage.*upmessage.*leftmessage.*rightmessage.*prawmessage; ED = 0; [temp disp] = max(belief,[],3);%选择可信度最大的作为最终视差求取结果 for row = 1:nr%data term for col = 1:nc ED = ED+prawcost(row,col,disp(row,col)); end end disp = disp-1;%数组下表从1开始,而视差从0开始 ES = 0; temp1 = zeros(size(disp)); temp1(1:nr-1,:) = disp(1:nr-1,:)-disp(2:nr,:); tmp1 = temp1~=0&downgradient;%列方向上的大梯度平滑项 ES = ES+length(find(tmp1))*s; tmp1 = temp1~=0&~downgradient;%列方向上的小梯度平滑项 ES = ES+length(find(tmp1))*s*P; temp2 = zeros(size(disp)); temp2(:,1:nc-1) = disp(:,1:nc-1)-disp(:,2:nc); tmp2 = temp2~=0&rightgradient;%行方向上的大梯度平滑项 ES = ES+length(find(tmp2))*s; tmp2 = temp2~=0&~rightgradient;%行方向上的小梯度平滑项 ES = ES+length(find(tmp2))*s*P; E = [E;[ED ES]]; eval(['imwrite(uint8(disp*scalar),''bpadisp' num2str(iteriter) 's' num2str(s) 'P' num2str(P) 'T' num2str(T) 'AD' num2str(Truncation) imagename '.png'')']);%将视差结果存储为图像 end %beleif propagation 迭代过程 eval(['save bpas' num2str(s) 'P' num2str(P) 'T' num2str(T) 'AD' num2str(Truncation) imagename 'energy.mat E']); toc %%%%%%%%%%%%%belief propagation accelerated left image is the reference%