function labels = flat(f,K)
%
% Tom's Topographic Primal Sketch: flat finder
% Computes root median square average value 
% On input:
%    f:  2D image function
%    K:  window size for median computation
% On output:
%    labels:  root median square average value at pixel 
%             for KxK window 
% Calling sequence:
%    labels = flat(f,K)

[m,n] = size(f);
labels = zeros(m,n); Karray = ones(K,K);
Wradius = (K-1)/2; Ksquared = K*K;
rmin = (K+1)/2;   cmin = rmin;
rmax = m-rmin+1; cmax = n-cmin+1;

for r = rmin:rmax
  for c = cmin:cmax
    W = f(r-Wradius:r+Wradius,c-Wradius:c+Wradius);
    Wmedian = median(W(1:Ksquared));
    dif = W-Wmedian*Karray; dif2 = dif.*dif;
    Wrms = sqrt(sum(dif2(1:Ksquared)))/Ksquared;
    labels(r,c) = Wrms;
    end
  end
end
