function t = tpslab(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
%
% tpslab: compute topographic sketch curvature labels based on bicubic fit
% On input:
%     k1..k10: bicubic coefficients
% On output:
%     t: topographic sketch labels (see Jain, Computer Vision)
% Call:
%      t = tpslab(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
%

[m,n] = size(k1);
fr = k2;
fc = k3;
gradf = sqrt(fr.*fr + fc.*fc);
[l1,l2,v1,v2,v3,v4] = Hessianf(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);

for i1 = 1:m
  for j1 = 1:n
    if gradf(i1,j1)<eps & l1(i1,j1)<0 & l2(i1,j1)<0
	  t(i1,j1) = 3;  % peak
	elseif gradf(i1,j1)<eps & l1(i1,j1)>0 & l2(i1,j1)>0
	  t(i1,j1) = 9;  % pit
	elseif gradf(i1,j1)>=eps & l1(i1,j1)<0 & ...
	       abs(fr(i1,j1)*v1(i1,j1)+fc(i1,j1)*v2(i1,j1))<eps
	  t(i1,j1) = 1;  % ridge
	elseif gradf(i1,j1)>=eps & l2(i1,j1)<0 & ...
	       abs(fr(i1,j1)*v3(i1,j1)+fc(i1,j1)*v4(i1,j1))<eps
	  t(i1,j1) = 1;  % ridge
	elseif gradf(i1,j1)<eps & l1(i1,j1)<0 & abs(l2(i1,j1))<eps
	  t(i1,j1) = 1;  % ridge
	elseif gradf(i1,j1)>=eps & l1(i1,j1)>0 & ...
	       abs(fr(i1,j1)*v1(i1,j1)+fc(i1,j1)*v2(i1,j1))<eps
	  t(i1,j1) = 8;  % valley
	elseif gradf(i1,j1)>=eps & l2(i1,j1)>0 & ...
	       abs(fr(i1,j1)*v3(i1,j1)+fc(i1,j1)*v4(i1,j1))<eps
	  t(i1,j1) = 8;  % valley
	elseif gradf(i1,j1)<eps & l1(i1,j1)>0 & abs(l2(i1,j1))<eps
	  t(i1,j1) = 8;  % valley
	elseif gradf(i1,j1)<eps & l1(i1,j1)*l2(i1,j1)<0
	  t(i1,j1) = 4;  % saddle
	elseif gradf(i1,j1)<eps & abs(l1(i1,j1))<eps & ...
	       abs(l2(i1,j1))<eps
	  t(i1,j1) = 5;    % flat
	else
	  t(i1,j1) = 2;    % anything else!
	 end
  end
end
