I find this question extremely interesting, but at the same time it's evident that it's also very complicated. On the other hand, I believe that the winning strategy is always to start from something simple and then, once consolidated, generalize it at will until you reach your goals.
So, assuming that I'm not an expert in the field, what I can say is that it often happens that we've a set of points to which we want to associate a function in one or more variables that describes them in the best possible way. A first question to ask is whether we want a function that passes absolutely through these points (interpolation) or if we want a function that passes us as close as possible, not necessarily that touches all the points (fitting).
In both cases, then, it's necessary to establish the basic functions to be used, usually polynomial, rational or trigonometric. Furthermore, it must be established whether this function is to be defined elementary or piecewise and in the latter case it's necessary to establish how to connect these sections, usually something smooth is desired, etc etc.
The answer to this list of questions depends on many factors, first of all the context in which one operates and the objectives we have set. Specifically, since the starting data may be affected by measurement errors, by noise, I would say that a curve fitting is preferable and in particular I thought that it was advantageous to rely on the Discrete Fourier Transform to quickly find a least squares trigonometric fit to the data.
In light of all this, by defining the following routine in Mathematica:
TrigFit[v_, m_, {x_, x0_, x1_}] := Module[{fc, fs, i, imax, j, n, p, t},
n = Length[v]; imax = Min[m, Floor[(n - 1) / 2]]; t = 2 π (x - x0)/(x1 - x0);
fc = Table[Sum[v[[j + 1]] Cos[2 π i j / n] / n, {j, 0, n - 1}], {i, 0, imax}];
fs = Table[Sum[v[[j + 1]] Sin[2 π i j / n] / n, {j, 0, n - 1}], {i, 0, imax}];
p = fc[[1]] + 2 Sum[fc[[i + 1]] Cos[i t] + fs[[i + 1]] Sin[i t], {i, 1, imax}];
Expand@Rationalize[p, 0.01]]
we can begin with an equally spaced distribution of points of which we wish to determine a function in a variable with at most $m$ modes and defined on the period considered:
pts = Table[N[Cos[x] / (5/4 + Sin[x]) + Sin[x] / (5/4 - Cos[x])],
{x, 0, 2 π - 2 π / 11, 2 π / 11}];
f[x_] = TrigFit[pts, 4, {x, 0, 2 π}];
Show[Plot[f[x], {x, 0, 2 π}],
ListPlot[pts, DataRange -> {0, 2 π - 2 π / 11}, PlotStyle -> Red],
AxesLabel -> {"x", "y"}, GridLines -> Automatic]
TeXForm@f[x]
$\quad\quad\quad\quad\quad\quad$
$\sin (x)+\frac{1}{4} \sin (3 x)+\frac{1}{4} \sin (4 x)+\cos (x)-\frac{1}{4} \cos (3 x)$
Therefore, we can raise the bar by considering a distribution of points in the $\langle x,\,y \rangle$ plane of which we wish to determine a $m$ modes parameterization of the closed curve that best describes them:
pts = Table[N[{JacobiP[5, 1, 2, Cos[t]], JacobiP[5, 1, 2, Sin[t]]}],
{t, 0, 2 π, π / 25}];
{x[t_], y[t_]} = TrigFit[#, 3, {t, 0, 2 π}]& /@ Transpose[pts];
Show[ParametricPlot[{x[t], y[t]}, {t, 0, 2 π}],
ListPlot[pts, PlotStyle -> Red],
AxesLabel -> {"x", "y"}, GridLines -> Automatic]
TeXForm@{x[t], y[t]}
$\quad\quad\quad\quad\quad\quad$
$\left\{-\frac{2 \sin (t)}{5}+\frac{2}{5} \sin (2 t)-\frac{8}{9} \sin (3 t)+\frac{25 \cos (t)}{4}-\frac{13}{4} \cos (2 t)+\frac{34}{7} \cos (3 t)-\frac{17}{9},\frac{584 \sin (t)}{97}-\frac{3}{4} \sin (2 t)-\frac{60}{13} \sin (3 t)+\frac{6 \cos (t)}{13}+\frac{44}{13} \cos (2 t)-\frac{8}{7} \cos (3 t)-\frac{101}{50}\right\}$
It's not all. The power of this approach is that it extends even to closed three-dimensional curves, everything always depends on having a set of sufficiently accurate coordinates. For example, in Mathematica we can find the coordinates of some knots:
pts = First@@KnotData[{8, 3}, "ImageData"];
{x[t_], y[t_], z[t_]} = TrigFit[#, 10, {t, 0, 2 π}]& /@ Transpose[pts];
ParametricPlot3D[{x[t], y[t], z[t]}, {t, 0, 2 π}, Axes -> False,
Boxed -> False, ImageSize -> 400, PlotRange -> All,
PlotStyle -> Directive[Orange, Specularity[White, 70]]
] /. Line[pts_, rest___] :> Tube[pts, 1/25, rest]
TeXForm@{x[t], y[t], z[t]}
$\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$
$\left\{-\frac{\sin (t)}{9}+\frac{2}{9} \sin (3 t)+\frac{1}{17} \sin (5 t)+\frac{1}{14} \sin (7 t)-\frac{\cos (t)}{8}-\frac{2}{81} \cos (2
t)-\frac{2}{11} \cos (3 t)+\frac{2}{23} \cos (5 t)-\frac{2}{43} \cos (7 t)+\frac{2}{71} \cos (10 t),\frac{\sin (t)}{9}-\frac{4}{11} \sin (2
t)-\frac{1}{8} \sin (3 t)+\frac{2}{57} \sin (5 t)-\frac{2}{57} \sin (6 t)-\frac{2 \cos (t)}{19}-\frac{2}{75} \cos (2 t)-\frac{2}{13} \cos (3 t)-\frac{2}{21} \cos (4 t)-\frac{1}{39} \cos (5 t)+\frac{1}{8},-\frac{1}{44} \sin (5 t)-\frac{1}{11} \sin (8 t)-\frac{1}{28} \cos (5 t)-\frac{2}{77} \cos (8 t)\right\}$
At this point, once we understand the academic examples, we can finally put them into practice by considering something similar to what you mentioned in the question. In particular, considering a two-dimensional figure with sufficiently well defined and closed edges, i.e. something like this:
$\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$
we can import it into Mathematica to extract the edge points:
img = Binarize@Import["https://i.stack.imgur.com/RpSbo.png"];
lines = Cases[Normal@ListContourPlot[Reverse@ImageData[img], Contours -> 1], _Line, -1];
pts = Reverse@SortBy[lines[[All, 1]], Length];
and consequently, by selecting the points of the three main edges, we obtain what is desired:
p = Table[TrigFit[#, 10, {t, 0, 2 π}]& /@ Transpose[pts[[i]]], {i, 3}];
ParametricPlot[p, {t, 0, 2 π}, AxesLabel -> {"x", "y"}, GridLines -> Automatic]
$\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$
where, as an example, a parameterization of the orange edge is:
TeXForm@p[[2]]
$\left\{-\frac{432 \sin (t)}{13}-\frac{29}{5} \sin (2 t)-\frac{14}{9} \sin (3 t)-\frac{1}{6} \sin (4 t)+\frac{8}{5} \sin (5 t)-\frac{76}{25} \sin (6 t)+\frac{13}{5} \sin (7 t)+\frac{19}{6} \sin (8 t)-\frac{3}{2} \sin (9 t)-\frac{2}{11} \sin (10 t)+\frac{48 \cos (t)}{5}-\frac{13}{4} \cos (2 t)+\frac{38}{13} \cos (3 t)-\frac{8}{5} \cos (4 t)-\frac{44}{15} \cos (5 t)-\frac{50}{11} \cos (6 t)+\frac{8}{13} \cos (7 t)+\frac{12}{13} \cos (8 t)+\frac{11}{7} \cos (9 t)-\frac{5}{6} \cos (10 t)+73,-\frac{41 \sin (t)}{4}+\frac{36}{7} \sin (2 t)-\frac{4}{3} \sin (3 t)-\frac{6}{5} \sin (4 t)-\frac{4}{7} \sin (5 t)+\frac{2}{3} \sin (6 t)+\frac{1}{7} \sin (7 t)+\frac{3}{5} \sin (8 t)-\frac{14}{11} \sin (9 t)+\frac{22}{21} \sin (10 t)-\frac{155 \cos (t)}{3}-\frac{48}{11} \cos (2 t)-\frac{22}{3} \cos (3 t)-\frac{26}{9} \cos (4 t)-\frac{4}{3} \cos (5 t)-\frac{2}{5} \cos (6 t)-\frac{52}{11} \cos (7 t)+\frac{5}{4} \cos (8 t)-\frac{4}{9} \cos (9 t)-\frac{1}{3} \cos (10 t)+\frac{1168}{9}\right\}$
Once acceptable parameterizations have been obtained for each edge of interest, we can proceed to work on the change of shape and in this regard I had spent a fairly detailed answer. Of course, this will haven't solved all the problems, but I hope it's a good starting point. ^_^
