function Gstar = calc_Gi_stack( Faces, restV )
%CALC_GI_STACK Given our mesh (faces, and REST vertex positions), calculate
% our stack of Gi Matrices. (Sometimes referred to as G* in the notes).

% We know from Equation 10 in our notes that there is a Gradient matrix Gi, which in the
% 2D case is a 4x2n matrix, where n is the number of vertices. 
% We also know from Equation 17 in our notes that Gi only depends on the rest pose matrices. 
% This function returns G* (see Equation 28), which is a 4mx2n matrix,
% where m is the number of triangles. This is a stack of matrices G1, G2... Gm. 
    facesSize = size(Faces);
    restVSize = size(restV);
    MiTildes = zeros(2, 2, facesSize(1)); 
    Dis = zeros(2, restVSize(1), facesSize(1)); %incidence matrices
    Gstar = zeros(4*facesSize(1), 2*restVSize(1)); %contains stack of all Gi matrices
    
    for i = 1:facesSize(1)
        v1 = transpose(restV(Faces(i, 1),:));
        v2 = transpose(restV(Faces(i, 2),:));
        v3 = transpose(restV(Faces(i, 3),:));
        MiTildes(:,:,i) = [v3 - v1, v2 - v1];
        
        %Subtract the first vertex from the other two (for incidence)
        Dis(:,Faces(i, 1),i) = [-1; -1];
        Dis(1,Faces(i, 3),i) = 1;
        Dis(2,Faces(i, 2),i) = 1;
        
        GstarIdx = 4*(i-1) + 1;
            
        Gstar(GstarIdx:GstarIdx+3,:) = kron((inv(MiTildes(:,:,i))')*Dis(:,:,i), eye(2));
    end

end

