Computing gradient of mean curvature on a mesh using gp toolbox, un... (2024)

% 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');

% 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');

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

⋮