Cacher is the code snippet organizer for pro developers

function [ deformedSkull ] = fitSkull(skullMesh, skinMesh)
% fitSkull
%   Detailed explanation goes here

%% Parameter preparation

% forehead:     1217
% between eyes:  857
% nose bridge:   946
% head right:    866
% head left:    5053
anchorIndices = [ 1217, 857, 946, 866, 5053 ];
% forehead:     4.5 mm
% between eyes: 7.0 mm
% nose bridge:  2.0 mm
% head right:   3.5 mm
% head left:    3.5 mm
anchorThickness = [ 0.45, 0.70, 0.20, 0.35, 0.35 ];
minThickness = 0.2;
maxThickness = 0.7;

numSkullVertices = numel(skullMesh.vidx);
freeVarIndices = 1:numSkullVertices;
freeVarIndices(anchorIndices) = [];

% Precomputed mesh tree to accelerate ray-mesh collision detection
skullMeshTree = opcodemesh(skullMesh.v', skullMesh.f');

%% Get default Euler-Lagrangian for shell minimization
% Euler-Lagrangian solution to the shell energy minimization
%
%    -Δx + Δ²x = 0
% 
% with replacing Laplacian(Δ) with discrete Laplace Beltrami operator

L = discreteLaplaceBeltrami(skullMesh);
stretchingPenalty = 1.0;
bendingPenalty = 1.0;
lagrangian = -(stretchingPenalty * L) + (bendingPenalty * L^2);

% Eliminate hard constraint terms from Lagrangian

anchorPositions = skullMesh.v(anchorIndices, :) + ...
    (skullMesh.n(anchorIndices, :) .* anchorThickness);

hardConstRhs = -lagrangian(:, anchorIndices) * anchorPositions;
hardConstLagrangian = lagrangian(:, freeVarIndices);

%% Iterations

% Initial skull deformation is given by hard constraint only
displacement = hardConstRhs \ hardConstLagrangian;

thickness = zeros(1, numel(freeVarIndices));
softConstWeights = zeros(1, numel(freeVarIndices));


end

function [ L ] = discreteLaplaceBeltrami(mesh)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

N = numel(mesh.vidx);

% Create cotangent angle weight with diagonal entires sums up to the 
% others of the same row ...
L = makeOneRingWeights(mesh, 'cotan', 0) / 2;
L(1:N+1:end) = -sum(L, 2);

% ... divided by voronoi-"mixed" vertex area.
area = vertexArea(mesh, [], 'mixed');
for i = 1:N
    L(i, :) = L(i, :) / area(i);
end

end