function labels = TTPS(f,K,zero_eps,flat_t)
%
% Tom's Topographic Primal Sketch
%    labels = TTPS(f,K,zero_eps,flat_t)
%
[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;
[fc fr] = gradient(f);
fgn = sqrt(fc.*fc + fr.*fr);

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;
    if Wrms<flat_t
      labels(r,c) = 1;       % flat
    else
      theta = atan2(fr,fc);
    end
  end
end
