function [e1,e2] = flat(f,K)
%
% flat: produces error of fit of KxK window at each pixel to a constant value
%     [e1,e2] = flat(f,K)
% On input:
%     f:   image function
%     K:   window size to evaluate
% On output:
%     e1:  RMS error
%     e2:  L1 norm
%          (returns -1 in unevaluated border pixels)

[m,n] = size(f);
e1 = -1*ones(m,n);
e2 = -1*ones(m,n);
imin = (K+1)/2; imax = m-imin+1;
jmin = imin; jmax = n-jmin+1;
Ksquared = K*K;

for i1 = imin:imax
  for j1 = jmin:jmax
    wimin = i1-imin+1; wimax = wimin+K-1;
    wjmin = j1-jmin+1; wjmax = wjmin+K-1;
    w = f(wimin:wimax,wjmin:wjmax);
    wmedian = median(w(1:Ksquared));
    dif = w-wmedian*ones(K,K);
    e1(i1,j1) = sum(sum(sqrt(dif.*dif)))/Ksquared;
    e2(i1,j1) = max(max(abs(dif)));
  end
end
