How to Calculate Biharmonic Distance in MMA?

Question

There are many kinds of distances. One of them is Biharmonic Distance where I got the the image below from:

The biharmonic examples are on the left, and the author kept his promise in the paper that this distance keeps a fantastic balance between local properties of the underlying space as well global properties. I am slowly working my way through the paper and terms to be able to write a Mathematica program that can calculate biharmonic distance, but I am not a professional Mathematician/Mechanical Engineer. So a lot of the terms I don't understand without great difficulty in the paper, and I thought that this would be nice to share because it would be great if Mathematica had a function like BiharmonicPointDistance[spatialDataPoints,pointSource] which would give the results in the images above. But there also might be shortcuts built into the Wolfram Language that would make recomputing this easier.

Examples of the family of curves that I am interested in are shown below:

which were all created with the code:

populateEllipse[center_, a_, b_, r_, pointQuantity_, maxAttempts_] :=
Module[
  {initialEllipse, ellipse, reg, pts = {}, n = 1, ellipseGraphics, parametricForms = {}, \[Theta], axes},
  ellipse = Disk[center, {a, b}];
  reg = RegionErosion[ellipse, r];
initialEllipse = {center[[1]] + a Cos[t], center[[2]] + b Sin[t]};
While[Length[pts] < pointQuantity && n <= maxAttempts,
    pts = RandomPointConfiguration[HardcorePointProcess[pointQuantity, 2 r, 2], reg];
    n++;
    If[Length[pts["Points"]] > 0, Break[]]; (* Adjusted to exit loop if pts found *)
  ];
pts = If[Length[pts["Points"]] > pointQuantity, Take[pts["Points"], pointQuantity], pts["Points"]];
(* Generate ellipses and store parametric forms *)
  ellipseGraphics = Graphics[{
    Style[ellipse, FaceForm[None], EdgeForm[Black]],
    (axes = r RandomReal[{r, 2r}, 2];
     [Theta] = RandomReal[{0, 2 [Pi]}];
     AppendTo[parametricForms, {#[[1]] + axes[[1]] Cos[t] Cos[[Theta]] - axes[[2]] Sin[t] Sin[[Theta]], #[[2]] + axes[[1]] Cos[t] Sin[[Theta]] + axes[[2]] Sin[t] Cos[[Theta]]}];
     GeometricTransformation[Disk[#, axes], RotationTransform[[Theta], #]]
    ) & /@ pts
  }];
(* Return both initial ellipse and parametric forms for further use *)
  {initialEllipse, parametricForms}
]
divyEllipsePerimeter[pointQuantity_,parametricForm_]:=
Module[{points,plot},
points=Cases[
  Normal@
    (plot = ParametricPlot[
         parametricForm, {t, 0, 2 Pi},
         Mesh -> pointQuantity, MeshFunctions -> {"ArcLength"}
       ]),
    Point[l_] -> l, Infinity]
];
(* Define a function to append a unique negative z-value to each tuple within a sublist *)
modifySublist[sublist_] :=
Module[{z = -RandomReal[]},
  Table[Append[sublist[[i]], z], {i, 1, Length[sublist]}]
];
ellipseQuantity=6;
samplePoints=300;
diskHeight=RandomReal[{.5,2}];
stored=populateEllipse[{0,0},8,5,.55,ellipseQuantity,50];
upperTopLayer=Table[Append[divyEllipsePerimeter[samplePoints,stored[[1]]][[i]], 0], {i, 1, samplePoints}];
lowerTopLayer=Map[modifySublist,Partition[divyEllipsePerimeter[samplePoints,stored[[2]]],samplePoints]];
topLayercentroids=Map[Mean,lowerTopLayer];
topLayer=Partition[Flatten[Join[upperTopLayer,lowerTopLayer,topLayercentroids]],3];
bottomLayer= Map[{#[[1]], #[[2]], #[[3]] -diskHeight} &, topLayer];
surfacePoints=Join[topLayer,bottomLayer];
ListPlot3D[topLayer,Axes->False,Boxed->False,ColorFunction -> "GrayTones"] 

which was motivated by Tad's answer from a previous post.