Just wanted to add that the error I mentioned is caused by the flipping of the sign of the surface normal of just one face by rearranging the order of the vertices in the face matrix
Very useful, except for a bug in the compute_curvature function:
compute_curvature will generate an error on line 75 ("dp = sum( normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );") for SOME surfaces. The error stems from E containing indices that are out of range which is caused by line 48 ("A = sparse(double(i),double(j),s,n,n);") where A's values eventually entirely make up the E matrix. The problem occurs when the i and j vectors create the same ordered pair twice in which case the sparse function adds the two s vector elements together for that matrix location resulting in a value that is too large to be used as an index on line 75. For example, if i = [1 1] and j = [2 2] and s = [3 4] then A(1,2) will equal 3 + 4 = 7.
The i and j vectors are created here:
i = [face(1,:) face(2,:) face(3,:)];
j = [face(2,:) face(3,:) face(1,:)];
The fact that your code seems to depend on the order of the vertices in the faces matrix worries me because the curvature should be the same regardless of the order, obviously. To be fair, I don't completely understand how your code works so perhaps the way it is written it works out to not matter except that it does certainly matter when it results in an index out of bounds error as previously described.
