stereoGlobalEnergyMin(leftImage, rightImage, dmax, matching, Alfa)

%*******************************************************************
% Region Based Stereo Matching Algorithm by Global Error Energy
% Minimization by Smoothing Functions method explanied in  the 
% "Obtaining Depth Maps From Color Images By Region Based Stereo
% Matching Algorithms"
% 
% It uses stereo color image pairs. Right camera image is loaded into XR array
% and left camera image is loaded into XL array. You can set Disparity search 
% range by DisparityRange.
% 
%
% Programmer: B. Baykant ALAG�Z
% Ver.1.0 date: 09.2008
%
% Adapted by: Raul Correal (Mar 2012) 
% ([e-mail])
%
% Input:
% - leftImage: left image
% - rightImage: right image
% - dmax: disparity range
% - matching: 1=for point, 2=for line, 3=for 3x3 square
% - Alfa: tolerance coefficient for eliminating unreliable estimation
%
% Output:
% - disparityx: disparity map
% - disparityReliable: disparity map with reliable disparities
% - disparityF: Median Filtered Disparity Map with Reliable Disparities
%
%******************************************************************

function [disparityx disparityReliable disparityF] = stereoGlobalEnergyMin(leftImage, rightImage, dmax, matching, Alfa)


% auto-presettings for the program
[m n p]=size(rightImage);
Edis=1000000*ones(m,n);
disparity=zeros(m,n);
XR=double(rightImage);
XL=double(leftImage);

% process by increasing disparity
for d=0:dmax
    fprintf ('Computing for disparity: %d\n',d);
% composing error energy matrix for every disparity.(Algorithm step:1)
for j=3+d:n-2
    for i=2:m-1
         if p==3
           %kareselfark(i,j-d)=(XR(i,j,1)-XL(i,j-d,1))^2+(XR(i,j,2)-XL(i,j-d,2))^2+(XR(i,j,3)-XL(i,j-d,3))^2;  
           if matching==1
           %point matching 
            ErrorEnergy(i,j-d)=(1/3)*[(XL(i,j,1)-XR(i,j-d,1))^2+(XL(i,j,2)-XR(i,j-d,2))^2+(XL(i,j,3)-XR(i,j-d,3))^2];
           elseif matching==2 
           %block matching with line type window
            ErrorEnergy(i,j-d)=(1/15)*[(XL(i,j,1)-XR(i,j-d,1))^2+(XL(i,j,2)-XR(i,j-d,2))^2+(XL(i,j,3)-XR(i,j-d,3))^2+(XL(i,j-1,1)-XR(i,j-1-d,1))^2+(XL(i,j-1,2)-XR(i,j-1-d,2))^2+(XL(i,j-1,3)-XR(i,j-1-d,3))^2+(XL(i,j+1,1)-XR(i,j+1-d,1))^2+(XL(i,j+1,2)-XR(i,j+1-d,2))^2+(XL(i,j+1,3)-XR(i,j+1-d,3))^2+(XL(i,j-2,1)-XR(i,j-2-d,1))^2+(XL(i,j-2,2)-XR(i,j-2-d,2))^2+(XL(i,j-2,3)-XR(i,j-2-d,3))^2+(XL(i,j+2,1)-XR(i,j+2-d,1))^2+(XL(i,j+2,2)-XR(i,j+2-d,2))^2+(XL(i,j+2,3)-XR(i,j+2-d,3))^2];
           else
               top=0;
               for k=i-1:i+1
                   for l=j-1:j+1
                     top=top+(XL(k,l,1)-XR(k,l-d,1))^2+(XL(k,l,2)-XR(k,l-d,2))^2+(XL(k,l,3)-XR(k,l-d,3))^2;  
                   end
               end
               ErrorEnergy(k,l-d)=(1/27)*top;
           end 
       else
           Disp('ERROR WARNING: Use RGB color image for stereo pair');
        end
    end
end
% applying smooting on error energy surfaces by iterative averaging
% filtering. (Algorithm step:2)
ErrorEnergyFilt=IterativeAveragingFilter(ErrorEnergy,5,[4 4]);

% selecting disparity which has minimum error energy.(Algorithm step:3)
[m1 n1]=size(ErrorEnergyFilt);
for k=1:m1
    for l=1:n1
       if Edis(k,l)>ErrorEnergyFilt(k,l)
            disparity(k,l)=d;
            Edis(k,l)=ErrorEnergyFilt(k,l);
        end        
    end
end
end
% clear 1000000 pre-setting in Edis
for k=1:m
    for l=1:n
       if Edis(k,l)==1000000
            Edis(k,l)=0;
        end        
    end
end

% extracting calculated zone
nx=n-dmax;
for k=2:m-1
    for l=2:nx-1
        disparityx(k,l)=disparity(k,l);
        %Edisx(k,l)=Edis(k,l);
        %regMapx(k,l)=regMap(k,l);
        XLx(k,l)=XL(k,l);
        XRx(k,l)=XR(k,l);
        top=0;
        for x=k-1:k+1
            for y=l-1:l+1
                top=top+(XL(x,y+disparity(k,l),1)-XR(x,y,1))^2+(XL(x,y+disparity(k,l),2)-XR(x,y,2))^2+(XL(x,y+disparity(k,l),3)-XR(x,y,3))^2;  
            end
        end
        Ed(k,l)=(1/27)*top;
    end
end

%calculates error energy treshold for reliablity of disparity
Toplam=0;
for k=1:m-1
    for l=1:nx-1
      Toplam=Toplam+Ed(k,l);       
    end
end
% Error threshold Ve
Ve=Alfa*(Toplam/((m-1)*(nx-1)));

EdReliable=Ed;
disparityReliable=disparityx;
Ne=zeros(m,nx);
for k=1:m-1
    for l=1:nx-1
       if Ed(k,l)>Ve
          % sets unreliable disparity to zero
          disparityReliable(k,l)=0;
          EdReliable(k,l)=0;
          Ne(k,l)=1; % indicates no-estimated state
        end        
    end
end

% calculating reliablities both raw disparity and filtered disparity
TopE=0;
TopER=0;
Sd=0;
for k=1:m-1
    for l=1:nx-1
          TopE=TopE+Ed(k,l);
          if Ne(k,l)==0          
             TopER=TopER+EdReliable(k,l);
             Sd=Sd+1;
          end          
    end
end

ReliablityE=((nx-1)*(m-1))/(TopE);
ReliablityER=(Sd)/(TopER);

% median filtering for repairment of occulations
%disparityF=IterativeAveragingFilter(disparity,5,[4 4]);
disparityF=medfilt2(disparityReliable,[5 5]);

fprintf ('******** Reliablity Report  ********** \n')
fprintf ('Reliablity of the disparity map: %f \n',ReliablityE)
fprintf ('Reliablity of the disparity map filtered: %f \n',ReliablityER)