function K = Gcurv(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
%   Finds the Gaussian curvature from pre-computed bicubic coef's
%
%    K = Gcurv(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)

[m,n] = size(k1);
onemat = ones(m,n);

%fr = k2 + 2*k4*x + k4*x + k5*y + 3*k7*x*x + 2*k8*x*y + k9*y*y;

fr = k2;
fr2 = fr.*fr;

%fc = k3 + k5.*x + 2.*k6.*y + k8.*x.*x + 2.*k9.*x.*y + 3.*k10.*y.*y;

fc = k3;
fc2 = fc.*fc;

%frr = 2.*k4 + 6.*k7.*x + 2.*k8.*y;
%frc = k5 + 2.*k8.*x + 2.*k9.*y;
%fcc = 2.*k6 + 2.*k9.*x + 6.*k10.*y;

frr = k4;
frc = k5;
fcc = k6;


Knum = (frr.*fcc - frc.*frc);
Kdenom = (onemat+fr2+fc2).*(onemat+fr2+fc2);
K = Knum./Kdenom;
K = K.*finite(K);


