function distances = biharmonic_distance(baseVertices, coord, triang)

nopts = length(coord) / 3;
notrg = length(triang) / 3;

coord = reshape(coord,3,nopts)';
triang = reshape(triang,3,notrg)';
triang = triang + 1; %1-based

bir = [1 1 1];
trgarea = zeros(1,notrg);
for i = 1:notrg
    trgx = [coord(triang(i,1),1) coord(triang(i,2),1) coord(triang(i,3),1)];
    trgy = [coord(triang(i,1),2) coord(triang(i,2),2) coord(triang(i,3),2)];
    trgz = [coord(triang(i,1),3) coord(triang(i,2),3) coord(triang(i,3),3)];
    aa = [trgx; trgy; bir];
    bb = [trgx; trgz; bir];
    cc = [trgy; trgz; bir];
    area = sqrt(det(aa)^2 + det(bb)^2 + det(cc)^2) / 2;
    trgarea(i) = area;
end

%find the approximate voronoi area of each vertex
AM = zeros(nopts, 1);
for i = 1:notrg
    AM(triang(i,1:3)) = AM(triang(i,1:3)) + trgarea(i)/3;
end

A = sparse(nopts, nopts);
%now construct the cotan laplacian
for i = 1:notrg
    for ii = 1:3
        for jj = (ii+1):3
            kk = 6 - ii - jj; % third vertex no
            v1 = triang(i,ii);
            v2 = triang(i,jj);
            v3 = triang(i,kk);
            e1 = [coord(v1,1) coord(v1,2) coord(v1,3)] - [coord(v2,1) coord(v2,2) coord(v2,3)];
            e2 = [coord(v2,1) coord(v2,2) coord(v2,3)] - [coord(v3,1) coord(v3,2) coord(v3,3)];
            e3 = [coord(v1,1) coord(v1,2) coord(v1,3)] - [coord(v3,1) coord(v3,2) coord(v3,3)];
            cosa = e2* e3'/sqrt(sum(e2.^2)*sum(e3.^2));
            sina = sqrt(1 - cosa^2);
            cota = cosa/sina;
            w = 0.5*cota;
            A(v1, v1) = A(v1, v1) - w;
            A(v1, v2) = A(v1, v2) + w;
            A(v2, v2) = A(v2, v2) - w;
            A(v2, v1) = A(v2, v1) + w;
        end
    end
end
A = -A;
T = sparse([1:nopts], [1:nopts], (AM), nopts, nopts, nopts);
A = A*inv(T)*A;

%make the first row&col of A zeros, this forces the first entry of
%the solution vector to be zero; takes care of inf many solns
A(1, :) = zeros(1,nopts);
A(:, 1) = zeros(nopts,1);
A(1, 1) = 1;

sampleidx = baseVertices + 1; %1-based
samplesize = length(sampleidx);

%make the right hand side vector(s) for the system
h = eye(nopts) - (1/nopts)* ones(nopts,nopts);
h(1,:) = zeros(1, nopts);
h = h(1:nopts,sampleidx);

%solve
g = full(A\h);
g = g - repmat(sum(g)/nopts,nopts,1); %%shift the solution so that add up to zero, but weighted
g = g(sampleidx,:);

%compute dists; this is a vectorization of d^2 = g(i,i)+g(j,j)-2g(i,j)
v1 = diag(g)';
v2 = ones(samplesize,1);
m = v2*v1;
distances = sqrt(m + m' - 2*g);
distances = distances';
distances = distances(:);

end
