Computing gradient of mean curvature on a mesh using gp toolbox, unexpected error. What am I missing?

I am trying to compute this quantity here H grad(H) where H is the mean curvture

The code breaks at the last line because there is a dimension mismatch between H and grad(H). size(grad(H))= [960 3] and size(V)=[162 3] and size(mean_curvture)=[162 3]. Why is the grad(H) has that dimension? where is the mistake?

% Load a mesh
[V, F] = load_mesh('sphere.off');  

% Compute Mean Curvature and Normal Vector

% 2. Compute Cotangent Laplacian (L) and Mass Matrix (M)
L = cotmatrix(V, F);
M = massmatrix(V, F, 'barycentric');

% 3. Compute Mean Curvature Vector (H)
H = -inv(M) * (L * V);

% 4. Compute the Magnitude of Mean Curvature (mean curvature at each vertex)
mean_curvature = sqrt(sum(H.^2, 2));

% Compute the Normal Vector
normals = per_vertex_normals(V, F); 

% Compute Gaussian Curvature
gaussian_curvature = discrete_gaussian_curvature(V,F);

% Compute Laplacian of Mean Curvature
laplacian_H =  L * mean_curvature;  

% Compute the Gradient of Mean Curvature
% Use gptoolbox's gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;

% Compute the new term: newthing
%H_squared = mean_curvature .^ 2;
%newthing = laplacian_H + ((H_squared - 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;
% Compute the new term: newthing
H_squared = mean_curvature .^ 2;
newthing = laplacian_H + ((H_squared - 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;





% Define small arc equation
h = 0.1;  % Small parameter h
vertex = V;  % V contains the vertices

% Compute the small arc
f_h = vertex + (mean_curvature .* normals) * h + newthing * (h^2 / 2);

% Visualization of small arc
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), 'FaceColor', 'cyan', 'EdgeColor', 'none');
hold on;
plot3(f_h(:,1), f_h(:,2), f_h(:,3), 'r.', 'MarkerSize', 10);
axis equal;
title('Small Arc Visualization');
xlabel('X');
ylabel('Y');
zlabel('Z');

the grad(H) script

% Step 1: Load a Mesh
[V, F] = load_mesh('sphere.off');

% Step 2: Compute Laplace-Beltrami Operator
L = cotmatrix(V, F);  % cotangent Laplace-Beltrami operator

% Step 3: Compute the Mean Curvature Vector (H)
H = -L * V;  % Mean curvature normal vector

% Step 4: Compute the Mean Curvature Magnitude
mean_curvature = sqrt(sum(H.^2, 2));

% Step 5: Compute the Gradient of Mean Curvature
% Use gptoolbox's gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;

% Step 6: Visualization or Further Processing
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), mean_curvature, 'EdgeColor', 'none');
axis equal;
lighting gouraud;
camlight;
colorbar;
title('Mean Curvature Gradient');
xlabel('X');
ylabel('Y');
zlabel('Z');

grad(V, F) from gp toolbox

function [G] = grad(V,F)
  % GRAD Compute the numerical gradient operator for triangle or tet meshes
  %
  % G = grad(V,F)
  %
  % Inputs:
  %   V  #vertices by dim list of mesh vertex positions
  %   F  #faces by simplex-size list of mesh face indices
  % Outputs:
  %   G  #faces*dim by #V Gradient operator
  %
  % Example:
  %   L = cotmatrix(V,F);
  %   G  = grad(V,F);
  %   dblA = doublearea(V,F);
  %   GMG = -G'*repdiag(diag(sparse(dblA)/2),size(V,2))*G;
  %
  %   % Columns of W are scalar fields
  %   G = grad(V,F);
  %   % Compute gradient magnitude for each column in W
  %   GM = squeeze(sqrt(sum(reshape(G*W,size(F,1),size(V,2),size(W,2)).^2,2)));
  %

  dim = size(V,2);
  ss = size(F,2);
  switch ss
  case 2
    % Edge lengths
    len = normrow(V(F(:,2),:)-V(F(:,1),:));
    % Gradient is just staggered grid finite difference
    G = sparse(repmat(1:size(F,1),2,1)',F,[1 -1]./len, size(F,1),size(V,1));
  case 3
    % append with 0s for convenience
    if size(V,2) == 2
      V = [V zeros(size(V,1),1)];
    end
  
    % Gradient of a scalar function defined on piecewise linear elements (mesh)
    % is constant on each triangle i,j,k:
    % grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
    % grad(Xijk) = Xj * (Vi - Vk)^R90 / 2A + Xk * (Vj - Vi)^R90 / 2A + 
    %             -Xi * (Vi - Vk)^R90 / 2A - Xi * (Vj - Vi)^R90 / 2A
    % where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
    % i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of 
    % 90 degrees
    %
    % renaming indices of vertices of triangles for convenience
    i1 = F(:,1); i2 = F(:,2); i3 = F(:,3); 
    % #F x 3 matrices of triangle edge vectors, named after opposite vertices
    v32 = V(i3,:) - V(i2,:);  v13 = V(i1,:) - V(i3,:); v21 = V(i2,:) - V(i1,:);
  
    % area of parallelogram is twice area of triangle
    % area of parallelogram is || v1 x v2 || 
    n  = cross(v32,v13,2); 
    
    % This does correct l2 norm of rows, so that it contains #F list of twice
    % triangle areas
    dblA = normrow(n);
    
    % now normalize normals to get unit normals
    u = normalizerow(n);
  
    % rotate each vector 90 degrees around normal
    %eperp21 = bsxfun(@times,normalizerow(cross(u,v21)),normrow(v21)./dblA);
    %eperp13 = bsxfun(@times,normalizerow(cross(u,v13)),normrow(v13)./dblA);
    eperp21 = bsxfun(@times,cross(u,v21),1./dblA);
    eperp13 = bsxfun(@times,cross(u,v13),1./dblA);
  
    %g =                                              ...
    %  (                                              ...
    %    repmat(X(F(:,2)) - X(F(:,1)),1,3).*eperp13 + ...
    %    repmat(X(F(:,3)) - X(F(:,1)),1,3).*eperp21   ...
    %  );
  
    GI =  ...
      [0*size(F,1)+repmat(1:size(F,1),1,4) ...
      1*size(F,1)+repmat(1:size(F,1),1,4) ...
      2*size(F,1)+repmat(1:size(F,1),1,4)]';
    GJ = repmat([F(:,2);F(:,1);F(:,3);F(:,1)],3,1);
    GV = [eperp13(:,1);-eperp13(:,1);eperp21(:,1);-eperp21(:,1); ...
       eperp13(:,2);-eperp13(:,2);eperp21(:,2);-eperp21(:,2); ...
       eperp13(:,3);-eperp13(:,3);eperp21(:,3);-eperp21(:,3)];
    G = sparse(GI,GJ,GV, 3*size(F,1), size(V,1));
  
    %% Alternatively
    %%
    %% f(x) is piecewise-linear function:
    %%
    %% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk
    %% ∇f(x) = ...                 = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk 
    %%                             = ∇φi fi + ∇φj fj + ∇φk) fk 
    %%
    %% ∇φi = 1/hjk ((Vj-Vk)/||Vj-Vk||)^perp = 
    %%     = ||Vj-Vk|| /(2 Aijk) * ((Vj-Vk)/||Vj-Vk||)^perp 
    %%     = 1/(2 Aijk) * (Vj-Vk)^perp 
    %% 
    %m = size(F,1);
    %eperp32 = bsxfun(@times,cross(u,v32),1./dblA);
    %G = sparse( ...
    %  [0*m + repmat(1:m,1,3) ...
    %   1*m + repmat(1:m,1,3) ...
    %   2*m + repmat(1:m,1,3)]', ...
    %  repmat([F(:,1);F(:,2);F(:,3)],3,1), ...
    %  [eperp32(:,1);eperp13(:,1);eperp21(:,1); ...
    %   eperp32(:,2);eperp13(:,2);eperp21(:,2); ...
    %   eperp32(:,3);eperp13(:,3);eperp21(:,3)], ...
    %  3*m,size(V,1));
  
    if dim == 2
      G = G(1:(size(F,1)*dim),:);
    end
  
  
    % Should be the same as:
    % g = ... 
    %   bsxfun(@times,X(F(:,1)),cross(u,v32)) + ...
    %   bsxfun(@times,X(F(:,2)),cross(u,v13)) + ...
    %   bsxfun(@times,X(F(:,3)),cross(u,v21));
    % g = bsxfun(@rdivide,g,dblA);

  case 4
    % really dealing with tets
    T = F;
    % number of dimensions
    assert(dim == 3);
    % number of vertices
    n = size(V,1);
    % number of elements
    m = size(T,1);
    % simplex size
    assert(size(T,2) == 4);

    if m == 1 && ~isnumeric(V)
      simple_volume = @(ad,r) -sum(ad.*r,2)./6;
      simple_volume = @(ad,bd,cd) simple_volume(ad, ...
        [bd(:,2).*cd(:,3)-bd(:,3).*cd(:,2), ...
        bd(:,3).*cd(:,1)-bd(:,1).*cd(:,3), ...
        bd(:,1).*cd(:,2)-bd(:,2).*cd(:,1)]);
      P = sym('P',[1 3]);
      V1 = V(T(:,1),:);
      V2 = V(T(:,2),:);
      V3 = V(T(:,3),:);
      V4 = V(T(:,4),:);
      V1P = V1-P;
      V2P = V2-P;
      V3P = V3-P;
      V4P = V4-P;
      A1 = simple_volume(V2P,V4P,V3P);
      A2 = simple_volume(V1P,V3P,V4P);
      A3 = simple_volume(V1P,V4P,V2P);
      A4 = simple_volume(V1P,V2P,V3P);
      B = [A1 A2 A3 A4]./(A1+A2+A3+A4);
      G = simplify(jacobian(B,P).');
      return
    end
  
    % f(x) is piecewise-linear function:
    %
    % f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk + φl(x) fl
    % ∇f(x) = ...                 = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk + ∇φl(x) fl
    %                             = ∇φi fi + ∇φj fj + ∇φk fk + ∇φl fl
    %
    % ∇φi = 1/hjk = Ajkl / 3V * (Facejkl)^perp
    %     = Ajkl / 3V * (Vj-Vk)x(Vl-Vk)
    %     = Ajkl / 3V * Njkl / ||Njkl||
    % 
  
    % get all faces
    F = [ ...
      T(:,1) T(:,2) T(:,3); ...
      T(:,1) T(:,3) T(:,4); ...
      T(:,1) T(:,4) T(:,2); ...
      T(:,2) T(:,4) T(:,3)];
    % compute areas of each face
    A = doublearea(V,F)/2;
    N = normalizerow(normals(V,F));
  
    % compute volume of each tet
    vol = volume(V,T);
  
    GI = ...
      [0*m + repmat(1:m,1,4) ...
       1*m + repmat(1:m,1,4) ...
       2*m + repmat(1:m,1,4)];
    GJ = repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1);
    GV = repmat(A./(3*repmat(vol,4,1)),3,1).*N(:);
    G = sparse(GI,GJ,GV, 3*m,n);
  
  end
end

the sphere.off file (mesh data)

[Umar]

Hi @ ehsan,

After analyzing your code and addressing your query, “Why is the grad(H) has that dimension? where is the mistake?”

When you are trying to solve grad_H, you multiply it by mean_curvature, which may lead to dimension mismatches unless mean_curvature is appropriately reshaped or indexed and the output of grad(V,F) is structured in such a way that it may not directly align with your expectations for vertex-wise operations if you do not not handle them correctly. So, to resolve the issue you should solve the gradient of mean curvature, so you can properly handle its application across vertices. Here is the snippet code,

% Load a mesh

[V, F] = load_mesh('sphere.off');  

% Compute Cotangent Laplacian (L) and Mass Matrix (M)

L = cotmatrix(V, F);

M = massmatrix(V, F, 'barycentric');

% Compute Mean Curvature Vector (H)

H = -inv(M) * (L * V);  % This should be [162, 3] for vertices

% Compute the Magnitude of Mean Curvature

mean_curvature = sqrt(sum(H.^2, 2));  % This will be [162, 1]

% Compute the Normal Vector

normals = per_vertex_normals(V, F); 

% Compute Gaussian Curvature

gaussian_curvature = discrete_gaussian_curvature(V,F);

% Compute Laplacian of Mean Curvature

laplacian_H = L * mean_curvature;  

% Correctly compute the Gradient of Mean Curvature

G = grad(V,F);  % G will be [#faces*dim, #vertices]

mean_curvature_vector = repmat(mean_curvature, 1, size(G, 1)/size(V, 

1)); % Reshape as necessary

grad_H = G * mean_curvature_vector'; % Ensure proper multiplication

% Compute new term: newthing

H_squared = mean_curvature .^ 2;

newthing = laplacian_H + ((H_squared - 2 * gaussian_curvature) .*   

mean_curvature) .* normals + mean_curvature .* grad_H;

% Define small arc equation

h = 0.1;  

vertex = V;  

f_h = vertex + (mean_curvature .* normals) * h + newthing * (h^2 / 2);

After implementing these changes, verify that your outputs align with expected physical behaviors and characteristics of curvature in your mesh.Also, consider adding unit tests for smaller meshes to validate your calculations independently before scaling up to larger datasets. Hope this helps, please let me know if you any further questions.

FYI, I couldn’t download “ sphere.off” since I am using MATLAB mobile.