function [N,MeanOfOuterRing] = CenterSurroundContrastNormalization(I,options);
%function [N,MeanOfOuterRing] = CenterSurroundContrastNormalization(I,options);
%
%returns a new image, the same size and shape as N, where the individual layers have
%been normalized in energy by the surrounding area.  The center surround region is a square,
% not a circle (doughnut).

d.CSurround.InnerRadius = 8;
d.CSurround.OuterRadius = 16;
d.bUsePctileRegularization = 1;
d.PctileRegularization = .6;
d.NormalizationType = 'DivideByMeanOfOuterRing';

options = ResolveMissingOptions(options,d);
I = padimage(I,options.CSurround.OuterRadius,'symmetric');
SmallSumFilter = MultiSumFilter(I,options.CSurround.InnerRadius);
BigSumFilter = MultiSumFilter(I,options.CSurround.OuterRadius);
SmallSumFilter = unpadimage(SmallSumFilter,options.CSurround.OuterRadius);
BigSumFilter = unpadimage(BigSumFilter,options.CSurround.OuterRadius);

OuterRing = BigSumFilter - SmallSumFilter;
AreaOfOuterRing = (options.CSurround.OuterRadius + 1)^2 - (options.CSurround.InnerRadius + 1)^2;
MeanOfOuterRing = OuterRing / AreaOfOuterRing;

I = unpadimage(I,options.CSurround.OuterRadius);
if(options.bUsePctileRegularization)
   Icols = reshape(I,[size(I,1)*size(I,2),size(I,3)]);
   for i = 1:size(I,3);
      S = sort(Icols(:,i));
      idx = floor(options.PctileRegularization * length(S));
      Reg(1,1,i) = S(idx);
   end
   MeanOfOuterRing = MeanOfOuterRing + repmat(Reg,[size(I,1),size(I,2),1]);
end

N = I ./ MeanOfOuterRing;

function o = MultiSumFilter(i,param1);
o = i;
for j = 1:size(i,3)
  o(:,:,j) = sumfilter(i(:,:,j),param1);
end