function saMat = SSDB_Parallelism(X,options);
%
%function saMat = SSDB_Parallelism(X,options);

 d.PatchSize = 5;
 % d.FirstPoolSize = 1;
 d.CorrespondanceRange = 19; %MUST BE ODD
% d.CorrespondanceRange = 7;
 d.SpatialReductionSize = [8,8];
 d.RelativeTransformReductionSize = [5,5];
 d.OrigSize = [128 128];
 d.UseLinearVersionOfParallelism = 0;
 HeqOpts.FastApproximation = 0;
 
 if(not(exist('options')))
   options = [];
 end
 options = ResolveMissingOptions(options,d);
 nEx = size(X,2);
 for i = 1:nEx
   fprintf('SSDB_Parallelism image %d of %d\r',i,nEx);
   stim = reshape(X(:,i),options.OrigSize);
   stim = LocalHistogramEqualize(stim,[10 10],50,HeqOpts);
   % stim4 = FourWayFilter(stim);
   stim2 = TwoWayFilter(stim);
   [M] = AllRelativeCorrelations(stim2,options.CorrespondanceRange,options.PatchSize,options);
   [M] = LocalRelativeDistortion(M,options);
   [M] = SpatialConsensus(M,options);
   [M] = BestLocalTransformation(M,options);
   if(i == 1);
     saMat.X = zeros(length(M(:)),nEx);
   end
   saMat.X(:,i) = M(:);
 end
 
function M = BestLocalTransformation(M,options)
t = options.RelativeTransformReductionSize;
s = size(M);
N = zeros(s(1),s(2),t(1),t(2));
iset = 1:s(3);
jset = 1:s(4);
for i = 1:t(1)
  for j = 1:t(2)
    pooli = find( ((iset / s(3)) > (i-1)/t(1)) .* ((iset / s(3)) <= (i)/t(1)));
    poolj = find( ((jset / s(4)) > (j-1)/t(2)) .* ((iset / s(4)) <= (j)/t(2)));
    if(options.UseLinearVersionOfParallelism)
      N(:,:,i,j) = sum(sum(M(:,:,pooli,poolj),3),4);
    else
      N(:,:,i,j) = max(max(M(:,:,pooli,poolj),[],3),[],4);
    end
  end
end
 
function [M] = SpatialConsensus(M,options);
%function [M] = SpatialConsensus(M,options);
sumfiltersize = (ceil(options.OrigSize ./ options.SpatialReductionSize));
s = size(M);
for j = 1:s(4)
  M(:,:,:,j) = sumfilter(M(:,:,:,j),sumfiltersize);
end
r = options.SpatialReductionSize;
N = zeros(r(1),r(2),s(3),s(4));
for i = 1:s(3)
  for j = 1:s(4)
    N(:,:,i,j) = imresize(M(:,:,i,j),options.SpatialReductionSize,'nearest');
  end
end
M = N;

function [M] = AllRelativeCorrelations(I,range,patchsize,options);
 numranges = size(range,1);
 rangerange = 2*options.CorrespondanceRange + 1;
 M = zeros(size(I,1),size(I,2),rangerange,rangerange);
 M = M + inf;
 s = size(I);
 D = zeros(size(I,1), size(I,2), 2);
 % D = zeros(size(I,1), size(I,2), 4);
 n =0;
 for i = 1:rangerange,
   transi = i - options.CorrespondanceRange - 1;
   if(transi == 0), continue;, end; % (null translations are too unstable)
   % fprintf('correlation %.3d of about %.3d \r',n,floor((rangerange^2) / 2));
   for j = floor(rangerange/2):rangerange,
   % for j = 1:rangerange,
     n = n+1;
     transj = j - options.CorrespondanceRange - 1;
     if(transj == 0), continue;, end; % (null translations are too unstable)
     SHIFTI = TranslateIm(I,[transi,transj]);
     % for iOri = 1:4
     %   D(:,:,iOri) = conv2(ones(patchsize,1),ones(1,patchsize),(I(:,:,iOri)-SHIFTI(:,:,iOri)).^2,'same');
     %   % D(:,:,iOri) = sumfilter((I(:,:,iOri)-SHIFTI(:,:,iOri)).^2,floor(patchsize/2));
     %   % D(:,:,iOri) = I(:,:,iOri)+SHIFTI(:,:,iOri);
     % end
     for iOri = 1:2
       D(:,:,iOri) = conv2(ones(patchsize,1),ones(1,patchsize),(I(:,:,iOri)-SHIFTI(:,:,iOri)).^2,'same');
       % D(:,:,iOri) = sumfilter((I(:,:,iOri)-SHIFTI(:,:,iOri)).^2,floor(patchsize/2));
       % D(:,:,iOri) = I(:,:,iOri)+SHIFTI(:,:,iOri);
     end
     M(:,:,i,j) = sum(D,3);
   end
 end
 for i = 1:rangerange,
   transi = i - options.CorrespondanceRange - 1;
   if(transi == 0), continue;, end; % (null translations are too unstable)
   for j = 1:(floor(rangerange/2)-1),
     transj = j - options.CorrespondanceRange - 1;
     if(transj == 0), continue;, end; % (null translations are too unstable)
     M(:,:,i,j) = TranslateIm(M(:,:,rangerange-i+1,rangerange-j+1),[transi, transj],inf);
   end
 end
 M = M.^-1;
 
%  fprintf('correlation %.3d of about %.3d \n',n,floor((rangerange^2) / 2));
 
function [M] = LocalRelativeDistortion(M,options);
%all layers recieve only the maximum relative distorion among their neighbors   
if(options.UseLinearVersionOfParallelism)
 rangerange = 2*options.CorrespondanceRange + 1;
 for i = 2:(rangerange-1)
   M(:,:,i,:) = M(:,:,i,:) + M(:,:,i-1,:);
   M(:,:,i,:) = M(:,:,i,:) + M(:,:,i+1,:);
 end
 for j = 2:(rangerange-1)
   M(:,:,:,j) = M(:,:,:,j) + M(:,:,:,j-1);
   M(:,:,:,j) = M(:,:,:,j) + M(:,:,:,j+1);
 end
 M = M(:,:,2:3:end,2:3:end);
else
 rangerange = 2*options.CorrespondanceRange + 1;
 for i = 2:(rangerange-1)
   M(:,:,i,:) = max(M(:,:,i,:),M(:,:,i-1,:));
   M(:,:,i,:) = max(M(:,:,i,:),M(:,:,i+1,:));
 end
 for j = 2:(rangerange-1)
   M(:,:,:,j) = max(M(:,:,:,j),M(:,:,:,j-1));
   M(:,:,:,j) = max(M(:,:,:,j),M(:,:,:,j+1));
 end
 M = M(:,:,2:3:end,2:3:end);
end 

function outim = FourWayFilter(inim)
filters{1} = [ 1  1  1; 0  0  0;-1 -1 -1];
filters{2} = [-1  0  1;-1  0  1;-1  0  1];
filters{3} = [ 0  1  1;-1  0  1;-1 -1  0];
filters{4} = [ 1  1  0; 1  0 -1; 0 -1 -1];
outim = abs(ApplyFilterBank(inim,filters));

function outim = TwoWayFilter(inim)
filters{1} = [ 1  1  1; 0  0  0;-1 -1 -1];
filters{2} = [-1  0  1;-1  0  1;-1  0  1];
outim = abs(ApplyFilterBank(inim,filters));
