function HisteqIm = LocalHistogramEqualize(img,yx_wind,grain,options);
%function HisteqIm = LocalHistogramEqualize(img,yx_wind,grain,options);
%
%each pixel is approximately histogram equalized within the yx_window neighborhood.

d.FastApproximation = 0;
if(nargin < 4), options = [];,end
options = ResolveMissingOptions(options,d);

if(length(unique(img(:))) == 1)
  % fprintf('LHEQ-- No variation in stimulus\n');
  HisteqIm = ones(size(img));
  return
end

img = im2double(img);
img = (img - min(img(:))) / (max(img(:)) - min(img(:)));
img = img * (grain-1);
img = ceil(img) + 1;
IdxImg = padimage(img,[floor(yx_wind),ceil(yx_wind)],'symmetric');

if(not(options.FastApproximation))
  MLI = zeros(size(IdxImg,1), size(IdxImg,2), grain);
  nWindsY = size(IdxImg,1) - yx_wind(1) + 1;
  nWindsX = size(IdxImg,2) - yx_wind(2) + 1;
  T = zeros(size(IdxImg,1), size(IdxImg,2), grain);
  for i = 1:grain
    T(:,:,i) = (IdxImg == i);
    MLI(:,:,i) = sumfilter(T(:,:,i),[floor(yx_wind),ceil(yx_wind)]);
  end
  MLI = unpadimage(MLI,[floor(yx_wind),ceil(yx_wind)]);
  IdxImg = unpadimage(IdxImg,[floor(yx_wind),ceil(yx_wind)]);
  Cum = (cumsum(MLI,3) - MLI/2)/ prod(yx_wind);
  HisteqIm = DepthSelect(Cum,IdxImg);
else
  M = sumfilter(IdxImg,floor(yx_wind(1)/2)) / prod(yx_wind);
  % S = stdfilter(IdxImg,floor(yx_wind(1)/2));
  H = (IdxImg - M);
  % H = (H ./ S);
  H = im2double(H);
  HisteqIm = (H - min(H(:))) / (max(H(:)) - min(H(:)));
  HisteqIm = unpadimage(HisteqIm,[floor(yx_wind),ceil(yx_wind)]);
end
