function flabel = TPS(f,l,zero_eps,flat_t,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
%
% Topographic Primal Sketch
%       flabel = TPS(f,l,zero_eps,flat_t,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)

[m,n] = size(f);
table3 = zeros(11,11);
table3(7,7) = 7; table3(7,8) = 7;
table3(8,8) = 8; table3(8,9) = 9; table3(8,10) = 9;
table3(9,9) = 9;
table3(10,10) = 10; table3(10,9) = 9;
table3(11,11) = 11; table3(11,10) = 11;

flabel = zeros(m,n);
%[k1,k2,k3,k4,k5,k6,k7,k8,k9,k10] = bicubicm(f,l);
gn = gradnorm(f);
lhalf = (l-1)/2;
rmin = (l+1)/2;   cmin = rmin;
rmax = m-rmin+1; cmax = n-cmin+1;
for r = rmin:rmax
  for c = cmin:cmax
   H = [2*k4(r,c) k5(r,c); k5(r,c) 2*k6(r,c)];
   [fc fr] = gradient(f);
   [W,L] = eig(H);
   if abs(gn(r,c))<zero_eps & abs(L(1,1))<zero_eps & abs(L(2,2))<zero_eps
     flabel(r,c) = 1;
   else 
    if abs(abs(L(1,1))-abs(L(2,2)))<zero_eps & abs(L(1,1))>zero_eps
      Wt = -inv(H)*[1 0]';
      W(:,1) = [Wt(1) Wt(2)]';
      W(:,2) = [-W(2,1) W(1,1)]';
    end;
    theta(1) = atan2(W(2,1),W(1,1)); s1 = sin(theta(1)); c1 = cos(theta(1));
    theta(2) = atan2(W(2,2),W(1,2)); s2 = sin(theta(2)); c2 = cos(theta(2));
    A1 = k7(r,c)*(s1^3)+k8(r,c)*(s1^2)*c1+k9(r,c)*(c1^2)*s1+k10(r,c)*(c1^3);
    B1 = k4(r,c)*(s1^2)+k5(r,c)*s1*c1+k6(r,c)*(c1^2);
    C1 = k2(r,c)*s1+k3(r,c)*c1;
    A2 = k7(r,c)*(s2^3)+k8(r,c)*(s2^2)*c2+k9(r,c)*(c2^2)*s2+k10(r,c)*(c2^3);
    B2 = k4(r,c)*(s2^2)+k5(r,c)*s2*c2+k6(r,c)*(c2^2);
    C2 = k2(r,c)*s2+k3(r,c)*c2;
    maxrho1 = abs(1+c1*(1-c1));
    maxrho2 = abs(1+c2*(1-c2));
    rt1 = (-B1+sqrt(B1*B1-3*A1*C1))/(3*A1);
    rt2 = (-B1-sqrt(B1*B1-3*A1*C1))/(3*A1);
    rt3 = (-B2+sqrt(B2*B2-3*A2*C2))/(3*A2);
    rt4 = (-B2-sqrt(B2*B2-3*A2*C2))/(3*A2);
    rts1(1,1) = real(rt1); rts1(1,2) = imag(rt1);
    rts1(2,1) = real(rt2); rts1(2,2) = imag(rt2);
    rts2(1,1) = real(rt3); rts2(1,2) = imag(rt3);
    rts2(2,1) = real(rt4); rts2(2,2) = imag(rt4);
    zcs = [];
    if (abs(rts1(1,2))<zero_eps) & (abs(rts1(1,1))<maxrho1)
      zcs = [zcs; rts1(1,1) theta(1)];
    end;
    if (abs(rts1(2,2))<zero_eps) & (abs(rts1(2,1))<maxrho1)
      zcs = [zcs; rts1(2,1) theta(1)];
    end;
    if (abs(rts2(1,2))<zero_eps) & (abs(rts2(1,1))<maxrho1)
      zcs = [zcs; rts2(1,1) theta(2)];
    end;
    if (abs(rts2(2,2))<zero_eps) & (abs(rts2(2,1))<maxrho1)
      zcs = [zcs; rts2(2,1) theta(2)];
    end;
    k = [k1(r,c),k2(r,c),k3(r,c),k4(r,c),k5(r,c),k6(r,c),k7(r,c), ...
         k8(r,c),k9(r,c),k10(r,c)];
    [zcnt,two] = size(zcs);
    if zcnt==0   % no zero crossing in pixel
      label = zcnt0(gn(r,c),fc(r,c),fr(r,c),k,L(2,2),L(1,1),zero_eps);
    elseif zcnt == 1                 % larger eigenvalue
      label = zcnt1(zcs(1,1),zcs(1,2),k,zero_eps);
    elseif zcnt == 2
      label1 = zcnt1(zcs(1,1),zcs(1,2),k,zero_eps);
      label2 = zcnt1(zcs(2,1),zcs(2,2),k,zero_eps);
      if label1 == 0
        label = label2;
      elseif label2 == 0
        label = label1;
      else 
        label = table3(label1,label2);
      end;
    else label = 0;
    end;
    flabel(r,c) = label;
   end;
  end;
end;
