% Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
% Licensed under the MIT License. See LICENSE.txt in the project root for full license information.

% Computes Gaussian, mean, max and min curvature values given a 2D
% grid/lattice representation of the input mesh
% 2D lattice is needed for the finite difference approximation needed for the differential geometry statistics

% Version History
%   2017/03/28  Marc Bosch       Initial version
function [K,H,Cmax1,Cmin1,normsurf] = fundmoments(X,Y,Z,lambda,M,N)

% Xtemp=X;
% dx=1:size(Xtemp,2);
% dy=1:size(Xtemp,1);
% [X,Y]=meshgrid(dx,dy);


%individual metrics: from "perceptual quality metrics for 3D meshes..."
%metrics surfature
%http://www.seas.upenn.edu/~cis70005/cis700sl7pdf.pdf
%http://www.seas.upenn.edu/~cis70005/cis700sl8pdf.pdf
%gaussian,mean,min,max curvature
%mean=(GL-2FM+EN)/(2(EG-F^2))
%gaussian=(LN-M^2)/(EG0F^2)
partderiv1_h(:,:,1) = imfilter(X,[-1/2 0 1/2],'replicate');%central difference
aux = imfilter(X,[-1 1],'replicate');%for the edges of the lattice
partderiv1_h(:,1,1)=aux(:,1);partderiv1_h(:,end,1)=aux(:,end-1);
partderiv2_h(:,:,1) = imfilter(partderiv1_h(:,:,1),[-1/2 0 1/2],'replicate');
aux = imfilter(partderiv1_h(:,:,1),[-1 1],'replicate');
partderiv2_h(:,1,1)=aux(:,1);partderiv2_h(:,end,1)=aux(:,end-1);
partcrossderiv_v(:,:,1) = imfilter(partderiv1_h(:,:,1),[-1/2 0 1/2]','replicate');
aux = imfilter(partderiv1_h(:,:,1),[-1 1]','replicate');
partcrossderiv_v(1,:,1)=aux(1,:);partcrossderiv_v(end,:,1)=aux(end-1,:);
partderiv1_v(:,:,1) = imfilter(X,[-1/2 0 1/2]','replicate');
aux = imfilter(X,[-1 1]','replicate');
partderiv1_v(1,:,1)=aux(1,:);partderiv1_v(end,:,1)=aux(end-1,:);
partcrossderiv_h(:,:,1) = imfilter(partderiv1_v(:,:,1),[-1/2 0 1/2],'replicate');
aux = imfilter(partderiv1_v(:,:,1),[-1 1],'replicate');
partcrossderiv_h(:,1,1)=aux(:,1);partcrossderiv_h(:,end,1)=aux(:,end-1);
partderiv2_v(:,:,1) = imfilter(partderiv1_v(:,:,1),[-1/2 0 1/2]','replicate');
aux = imfilter(partderiv1_v(:,:,1),[-1 1]','replicate');
partderiv2_v(1,:,1)=aux(1,:);partderiv2_v(end,:,1)=aux(end-1,:);


partderiv1_h(:,:,2) = imfilter(Y,[-1/2 0 1/2],'replicate');%central difference
aux = imfilter(Y,[-1 1],'replicate');%for the edges of the lattice
partderiv1_h(:,1,2)=aux(:,1);partderiv1_h(:,end,2)=aux(:,end-1);
partderiv2_h(:,:,2) = imfilter(partderiv1_h(:,:,2),[-1/2 0 1/2],'replicate');
aux = imfilter(partderiv1_h(:,:,2),[-1 1],'replicate');
partderiv2_h(:,1,2)=aux(:,1);partderiv2_h(:,end,2)=aux(:,end-1);
partcrossderiv_v(:,:,2) = imfilter(partderiv1_h(:,:,2),[-1/2 0 1/2]','replicate');
aux = imfilter(partderiv1_h(:,:,2),[-1 1]','replicate');
partcrossderiv_v(1,:,2)=aux(1,:);partcrossderiv_v(end,:,2)=aux(end-1,:);
partderiv1_v(:,:,2) = imfilter(Y,[-1/2 0 1/2]','replicate');
aux = imfilter(Y,[-1 1]','replicate');
partderiv1_v(1,:,2)=aux(1,:);partderiv1_v(end,:,2)=aux(end-1,:);
partcrossderiv_h(:,:,2) = imfilter(partderiv1_v(:,:,2),[-1/2 0 1/2],'replicate');
aux = imfilter(partderiv1_v(:,:,2),[-1 1],'replicate');
partcrossderiv_h(:,1,2)=aux(:,1);partcrossderiv_h(:,end,2)=aux(:,end-1);
partderiv2_v(:,:,2) = imfilter(partderiv1_v(:,:,2),[-1/2 0 1/2]','replicate');
aux = imfilter(partderiv1_v(:,:,2),[-1 1]','replicate');
partderiv2_v(1,:,2)=aux(1,:);partderiv2_v(end,:,2)=aux(end-1,:);


partderiv1_h(:,:,3) = imfilter(Z,[-1/2 0 1/2],'replicate');%central difference
aux = imfilter(Z,[-1 1],'replicate');%for the edges of the lattice
partderiv1_h(:,1,3)=aux(:,1);partderiv1_h(:,end,3)=aux(:,end-1);
partderiv2_h(:,:,3) = imfilter(partderiv1_h(:,:,3),[-1/2 0 1/2],'replicate');
aux = imfilter(partderiv1_h(:,:,3),[-1 1],'replicate');
partderiv2_h(:,1,3)=aux(:,1);partderiv2_h(:,end,3)=aux(:,end-1);
partcrossderiv_v(:,:,3) = imfilter(partderiv1_h(:,:,3),[-1/2 0 1/2]','replicate');
aux = imfilter(partderiv1_h(:,:,3),[-1 1]','replicate');
partcrossderiv_v(1,:,3)=aux(1,:);partcrossderiv_v(end,:,3)=aux(end-1,:);
partderiv1_v(:,:,3) = imfilter(Z,[-1/2 0 1/2]','replicate');
aux = imfilter(Z,[-1 1]','replicate');
partderiv1_v(1,:,3)=aux(1,:);partderiv1_v(end,:,3)=aux(end-1,:);
partcrossderiv_h(:,:,3) = imfilter(partderiv1_v(:,:,3),[-1/2 0 1/2],'replicate');
aux = imfilter(partderiv1_v(:,:,3),[-1 1],'replicate');
partcrossderiv_h(:,1,3)=aux(:,1);partcrossderiv_h(:,end,3)=aux(:,end-1);
partderiv2_v(:,:,3) = imfilter(partderiv1_v(:,:,3),[-1/2 0 1/2]','replicate');
aux = imfilter(partderiv1_v(:,:,3),[-1 1]','replicate');
partderiv2_v(1,:,3)=aux(1,:);partderiv2_v(end,:,3)=aux(end-1,:);

%https://en.wikipedia.org/wiki/Gaussian_curvature
%https://en.wikipedia.org/wiki/First_fundamental_form
%https://en.wikipedia.org/wiki/Second_fundamental_form
%vector normal to the surface of a face
partderiv1_h_buffer(:,1) = reshape(partderiv1_h(:,:,1),[M*N 1]);
partderiv1_h_buffer(:,2) = reshape(partderiv1_h(:,:,2),[M*N 1]);
partderiv1_h_buffer(:,3) = reshape(partderiv1_h(:,:,3),[M*N 1]);

partderiv2_h_buffer(:,1) = reshape(partderiv2_h(:,:,1),[M*N 1]);
partderiv2_h_buffer(:,2) = reshape(partderiv2_h(:,:,2),[M*N 1]);
partderiv2_h_buffer(:,3) = reshape(partderiv2_h(:,:,3),[M*N 1]);

partcrossderiv_h_buffer(:,1) = reshape(partcrossderiv_h(:,:,1),[M*N 1]);
partcrossderiv_h_buffer(:,2) = reshape(partcrossderiv_h(:,:,2),[M*N 1]);
partcrossderiv_h_buffer(:,3) = reshape(partcrossderiv_h(:,:,3),[M*N 1]);

partderiv1_v_buffer(:,1) = reshape(partderiv1_v(:,:,1),[M*N 1]);
partderiv1_v_buffer(:,2) = reshape(partderiv1_v(:,:,2),[M*N 1]);
partderiv1_v_buffer(:,3) = reshape(partderiv1_v(:,:,3),[M*N 1]);

partderiv2_v_buffer(:,1) = reshape(partderiv2_v(:,:,1),[M*N 1]);
partderiv2_v_buffer(:,2) = reshape(partderiv2_v(:,:,2),[M*N 1]);
partderiv2_v_buffer(:,3) = reshape(partderiv2_v(:,:,3),[M*N 1]);

partcrossderiv_v_buffer(:,1) = reshape(partcrossderiv_v(:,:,1),[M*N 1]);
partcrossderiv_v_buffer(:,2) = reshape(partcrossderiv_v(:,:,2),[M*N 1]);
partcrossderiv_v_buffer(:,3) = reshape(partcrossderiv_v(:,:,3),[M*N 1]);

%https://en.wikipedia.org/wiki/Gaussian_curvature
%https://en.wikipedia.org/wiki/First_fundamental_form
%https://en.wikipedia.org/wiki/Second_fundamental_form
%vector normal to the surface of a face
crossproduct_der=cross(partderiv1_h_buffer,partderiv1_v_buffer);
normsurf = crossproduct_der./repmat(sqrt(sum(crossproduct_der.*crossproduct_der,2)),[1 3]);
%dot products to obtain fundamental coefficients
%first coeffs
E = sum(partderiv1_h_buffer.*partderiv1_h_buffer,2);
F = sum(partderiv1_h_buffer.*partderiv1_v_buffer,2);
G = sum(partderiv1_v_buffer.*partderiv1_v_buffer,2);

%second coeffs
L = sum(partderiv2_h_buffer.*normsurf,2);
Mtemp = sum(partcrossderiv_v_buffer.*normsurf,2);
Ntemp = sum(partderiv2_v_buffer.*normsurf,2);

% Mean Curvature
Htemp = (E.*Ntemp + G.*L - 2.*F.*Mtemp)./(2*(E.*G - F.^2));
H = reshape(Htemp,[M N]);
% Gaussian Curvature
Ktemp = (L.*Ntemp - Mtemp.^2)./(E.*G - F.^2);
K = reshape(Ktemp,[M N]);
% min,max curvatures
Cmax1 = H + sqrt(H.^2 - K);
Cmin1 = H - sqrt(H.^2 - K);
%filter outliers (>lambda*std)

K(abs(K-mean(K(:)))>lambda*std(K(:)))=mean(K(:))+lambda*std(K(:));
H(abs(H-mean(H(:)))>lambda*std(H(:)))=mean(H(:))+lambda*std(H(:));
Cmax1(abs(Cmax1-mean(Cmax1(:)))>lambda*std(Cmax1(:)))=mean(Cmax1(:))+lambda*std(Cmax1(:));
Cmin1(abs(Cmin1-mean(Cmin1(:)))>lambda*std(Cmin1(:)))=mean(Cmin1(:))+lambda*std(Cmin1(:));

