function [row_out,col_out,SAD,adds]=fast_SAD(BSP,template,initial_row,initial_col,search_region,adds)
%[row_out,col_out,SAD,adds]=fast_SAD(BSP,template,initial_row,initial_col,search_region,adds)

%
%This function performs the "brute force" SAD block matching, given:
%BSP: which is the block sum pyramid of the search image.
%
%template: the template BSP
%
%initial_row/initial_col: the initial row and column (upper leftmost point) of the search area in the
%search image
%
%search_region: the number of pixels above/below and left/right of the initial location to search for a match
%adds: takes in the total number of adds to date
%
%Outputs:
%row_out/col_out: location of match
%adds: number of integer adds required to perform the matching.


%To take advantage of the MATLAB 7.0 Speedups in integer math perform the following conversions
%if needed (note this algorithm will not work in previous MATLAB versions).
BSP=uint16(BSP);
template=uint16(template);


%Template location is defined by its upper leftmost coordinate (rather than
%its center coordinate.
[row_size,col_size,dim]=size(BSP);
[template_row,template_col,dim1]=size(template);
boundary_region=template_row;
%Determine the boundaries for the current block to match
row_min=max(1,initial_row-search_region);
row_max=min(row_size-boundary_region-log2(boundary_region),initial_row+search_region-1);
col_min=max(1,initial_col-search_region);
col_max=min(col_size-boundary_region-log2(boundary_region),initial_col+search_region-1);
row_max=uint16(row_max);
row_min=uint16(row_min);
col_min=uint16(col_min);
col_max=uint16(col_max);
initial_row=uint16(initial_row);
initial_col=int16(initial_col);

%The rest of the algo is implemented here.
%The steps are to first check all the different blocks at this location
%Then check only the lowest block at each other pixel.  The lowest block
%corresponds to figure 6 in the paper.
%We have to keep track of the depth we've gone in and the current lower
%bound at each depth.
%After looking at all the pixels we look at the lowest lower bound.  If the
%lowest lower bound is not at the level K-1, we compute the bound at the
%next highest level for this block.  Then we look at the lowest bound
%again.  If this one is not at the K-1 level we look at the next higher
%block....
%We repeat this until the lower bound of a K-1 block is the lowest bound
%for the group.

%To check the sums at a particular level we use:
counter=1;
for i=0:(row_max-row_min)
    for j =0:(col_max-col_min)
        P(counter,:)=[i+row_min j+col_min sum(sum(abs(template(1:2^(dim-1):template_row,1:2^(dim-1):template_col,1)-BSP(row_min+i:2^(dim-1):row_min+i+template_row-1,col_min+j:2^(dim-1):col_min+j+template_col-1)))) 1];
        counter=counter+1;
    end
end

K=1;
[values,index]=min(P(:,3));
check=P(index,:);
while K~=dim    
i=check(1);
j=check(2);
K_temp=P(index,4)+1;
adds=adds+sum(sum(size(template(1:2^(dim-1-K):template_row,1:2^(dim-1-K):template_col,K))));
lower_bound=sum(sum(abs(template(1:2^(dim-1-K):template_row,1:2^(dim-1-K):template_col,K)-BSP(i:2^(dim-1-K):i+template_row-1,j:2^(dim-1-K):j+template_col-1,K))));
P(index,:)=[i j lower_bound K_temp];
[values,index]=min(P(:,3));
check=P(index,:);
K=check(4);
end

row_out=check(1);
col_out=check(2);
SAD=check(3);
