function [S,options,RowMatchMat] = SymmetryOperator_01(FiltIm, options);
%function [S,options,RowMatchMat] = SymmetryOperator_01(FiltIm, options);
%

d.HelpTextInfo = GetHelpTextInfo;
d.NormalizationRadius = 4;
d.NormalizationEpsilon = 10;
d.PoolRadius = 4;
d.PoolSampleRate = 8;
if(nargin <2), options=[];, end
options = ResolveMissingOptions(options,d);

FiltIm = LocalNormalization(FiltIm,options);
cccaPoolSupport = CalculatePoolSupport(size(FiltIm),options);
nRows = size(cccaPoolSupport,2);
nCols = size(cccaPoolSupport{1},2);
PoolReps = GetPoolRepresentative(FiltIm, cccaPoolSupport,options);
for iRow = 1:nRows
  for iCol = 1:nCols
    Sample = FiltIm(cccaPoolSupport{iRow}{iCol}{1},cccaPoolSupport{iRow}{iCol}{2},:);
    Sample = nLayerImage2MatrixOfPixels(Sample);
    for jCol = 1:nCols
      Rep = PoolReps(iRow,jCol,:);
      RowMatchMat(iCol,jCol,iRow) = BestMirrorMatch(Rep(:),Sample);
    end
  end
end
keyboard;

function PoolReps = GetPoolRepresentative(FiltIm, cccaPoolSupport,options);
nRows = size(cccaPoolSupport,2);
nCols = size(cccaPoolSupport{1},2);
PoolReps = zeros(nRows, nCols, size(FiltIm,3));
for iRow = 1:nRows
  for iCol = 1:nCols
    Sample = FiltIm(cccaPoolSupport{iRow}{iCol}{1},cccaPoolSupport{iRow}{iCol}{2},:);
    Sample = nLayerImage2MatrixOfPixels(Sample);
    Mag = sum(Sample .^2 ,1);
    [s,si] = sort(-Mag);
    PoolReps(iRow,iCol,:) = Sample(:,si(1));
  end
end

function cccaPoolSupport = CalculatePoolSupport(imsz,options)
nPoolRows = ceil(imsz(1) / options.PoolSampleRate);
nPoolCols = ceil(imsz(2) / options.PoolSampleRate);
for i = 1:nPoolRows
  for j = 1:nPoolCols
     xmin = max(1,(j-1)*options.PoolSampleRate);
     xmax = min(imsz(2),xmin + 2*options.PoolRadius);
     ymin = max(1,(i-1)*options.PoolSampleRate);
     ymax = min(imsz(2),ymin + 2*options.PoolRadius);
     cccaPoolSupport{i}{j}{1} = ymin:ymax;
     cccaPoolSupport{i}{j}{2} = xmin:xmax;
  end
end


function FiltIm = LocalNormalization(FiltIm,options);
p = padimage(FiltIm, options.NormalizationRadius,'symmetric');
Denom = sumfilter(FiltIm, options.NormalizationRadius);
Denom = sqrt(sum(Denom.^2,3));
Denom = max(options.NormalizationEpsilon,Denom);
FiltIm = FiltIm ./ repmat(Denom,[1 1 size(FiltIm,3)]);
FiltIm = unpadimage(FiltIm, options.NormalizationRadius);

function MatchScore = BestMirrorMatch(r,M)
%Assuming vertical mirror, only the last two filters need to be flipped.
M = M([1,4,3,2],:);
MatchScore = (max(r' * M)) / norm(r);