function K = Gcurvb(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
% Gaussian curvature based on bicubic approximation
%    K = Gcurvb(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
%
%   k array's generated by bicubicm

[m,n] = size(k1);
onemat = ones(m,n);
fr = bicubicfr(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);
fr2 = fr.*fr;
fc = bicubicfc(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);
fc2 = fc.*fc;
frr = bicubicfrr(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);
frc = bicubicfrc(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);
fcc = bicubicfcc(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10);

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