RELATED ARTICLE: Cantarella, J., Deguchi, T. and Shonkwiler, C. (2014), Probability Theory of Random Polygons from the Quaternionic Viewpoint. Commun. Pur. Appl. Math., 67: 1658-1699. https://doi.org/10.1002/cpa.21480. arXiv: arXiv:1206.3161v2 [math. DG].

I have mentioned many times before (1, 2, 3, 4, 5, 6, 7) that there's a correspondence between $n$-gons in the plane and points in either the Stiefel manifold $V_2(\mathbb{R}^n)$ of orthonormal 2-frames in $\mathbb{R}^n$ or the Grassmann manifold $G_2(\mathbb{R}^n)$ of 2-dimensional linear subspaces of $\mathbb{R}^n$. Of course, since geodesics in the Stiefel or Grassmann manifold are well-understood, this makes it easy to find nice parametrizations of morphs from one planar polygon to another.

In fact, the story generalizes: if instead of polygons in the plane you care about polygons in space (as many who are interested in modeling ring polymers in fact do), then you can again pass to a Stiefel or Grassmann manifold: either $V_2(\mathbb{C}^n)$, the Stiefel manifold of Hermitian orthonormal 2-frames in $\mathbb{C}^n$, or $G_2(\mathbb{C}^n)$, the Grassmannian of 2-dimensional complex linear subspaces of $\mathbb{C}^n$. This is first discussed in Hausmann and Knutson's paper and in the context of random polygons in my paper with Cantarella and Deguchi (arXiv version), but here's the short version:

We can represent a point in $V_2(\mathbb{C}^n)$ as an $n \times 2$ complex matrix with orthonormal columns. To get from such a matrix to a polygon, the map is actually rather simples: treat each row $(a_\ell, b_\ell)$, which is a pair of complex numbers, as a quaternion $q_\ell = a_\ell + j b_\ell$. Now, conjugating by the purely imaginary quaternion $i$ gives $\bar{q}_\ell i q_\ell$, which turns out to be a purely imaginary quaternion $x_\ell i + y_\ell j + z_\ell k$ and which we can treat as a three-dimensional vector $e_\ell = (x_\ell,y_\ell,z_\ell)$. But now just think of $e_\ell$ as the $\ell$th edge of a polygon; the magical thing is that the fact that the two columns of the original matrix were (Hermitian) orthogonal and the same length exactly guarantees that $\sum e_\ell = 0$, which is to say that the polygon closes up.

Here's some code for turning a point in the Stiefel manifold into either edges or vertices of a polygon:

<< Quaternions`

QuaternionEdges[StiefelForm_] := 
  Conjugate[#] ** Quaternion[0, 1, 0, 0] ** # & /@ (Quaternion @@ # & /@ Transpose[StiefelForm]);
QEdges[StiefelForm_] := List @@ #[[2 ;;]] & /@ QuaternionEdges[StiefelForm];
QVertices[StiefelForm_] := Accumulate[QEdges[StiefelForm]];

Now, geodesics in the Stiefel manifold are simple: given two such $n \times 2$ matrices, just rotate the first column of the first towards the first column of the second, while simultaneously rotating the second column of the first toward the second column of the second:

StiefelGeodesic[{A_, B_}, {C_, D_}, t_] := Module[{cPerp, dPerp, dist1, dist2},
   {cPerp, dPerp} = {Normalize[C - (C.A)*A], Normalize[D - (D.B)*B]};
   dist1 = ArcCos[A.C];
   dist2 = ArcCos[B.D];
   {Cos[t*dist1]*A + Sin[t*dist1]*cPerp, Cos[t*dist2]*B + Sin[t*dist2]*dPerp}
   ];

Now, in fact, the right way to think of a point in the Stiefel manifold is as a framed polygon, meaning each edge is really a triple of vectors: the edge vector itself (which I like to think of as the tangent vector) together with a normal and binormal vector, each of the same length as the edge vector and all three pairwise perpendicular. In terms of the above map from the Stiefel manifold to polygon space, we actually map $q_\ell$ (the $\ell$th row of the matrix thought of as a quaternion) to the triple $(\bar{q}_\ell i q_\ell, \bar{q}_\ell j q_\ell, \bar{q}_\ell k q_\ell)$.

Therefore, if given a particular polygon in space, you can (non-uniquely) turn this into a point in the Stiefel manifold by first framing the polygon (in this case, by just taking the cross product of consecutive tangents to be the normal, and then the cross product of tangent and normal to get the binormal):

PolygonToFramedPolygon[edges_] := Module[{normals, binormals},
   normals = Table[Norm[edges[[i]]] Normalize[Cross[edges[[i]], RotateLeft[edges[[i]]]]], {i, 1, Length[edges]}];
   binormals = Table[Norm[edges[[i]]] Normalize[Cross[edges[[i]], normals[[i]]]], {i, 1, Length[edges]}];
   Table[{edges[[i]], normals[[i]], binormals[[i]]}, {i, 1, Length[edges]}]
   ];

...and then inverting the $q \mapsto (\bar{q}_\ell i q_\ell, \bar{q}_\ell j q_\ell, \bar{q}_\ell k q_\ell)$ map (this map is 2-to-1, so this involves a choice):

MatrixToQuaternion[M_] := Module[{length, q0, q1, q2, q3},
   length = Norm[Transpose[M][[1]]];
   q0 = 1/2 Sqrt[length + Tr[M]];
   {q1, q2, q3} = 1/(4 q0) {M[[3, 2]] - M[[2, 3]], M[[1, 3]] - M[[3, 1]], M[[2, 1]] - M[[1, 2]]};
   {q0, q1, q2, q3}
   ];

By using a few auxiliary functions, I can create the single function PolygonVerticesToFrame which takes a list of vertices in space as input and outputs a point on the Stiefel manifold which maps to it:

Real4nToComplex2n[M_] := Transpose[Complex @@ # & /@ # & /@ (Partition[#, 2] & /@ Transpose[M])];
Complex2nToReal4n[M_] := Transpose[Flatten /@ (ReIm /@ Transpose[M])];

VerticesToEdges[verts_] := RotateLeft[verts] - verts;

PolygonVerticesToFrame[verts_] := 
  Orthogonalize@
   Real4nToComplex2n[
    Transpose[MatrixToQuaternion /@ PolygonToFramedPolygon[VerticesToEdges[verts]]]];

For this animation, I'm transforming the $(3,4)$ torus knot into the $(4,3)$ torus knot (obviously not by isotopies), so I get each of the two torus knots as a list of vertices and turn them into Stiefel manifold points:

Stereo3D[{x1_, y1_, x2_, y2_}] := {x1/(1 - y2), y1/(1 - y2), x2/(1 - y2)};

torus1 = Table[
   Stereo3D[
    1/Sqrt[2] RotationMatrix[\[Pi]/8, {{1, 0, 0, 0}, {0, 1, 0, 0}}].
       {Cos[3 \[Theta]], Sin[3 \[Theta]], Cos[4 \[Theta]], Sin[4 \[Theta]]}], 
    {\[Theta], 0., 2 \[Pi], 2 \[Pi]/1000}];

torus2 = Table[
   Stereo3D[
    1/Sqrt[2] RotationMatrix[-\[Pi]/72, {{1, 0, 0, 0}, {0, 1, 0, 0}}].
        {Cos[4 \[Theta]], Sin[4 \[Theta]], Cos[3 \[Theta]], Sin[3 \[Theta]]}], 
    {\[Theta], 0., 2 \[Pi], 2 \[Pi]/1000}];

torus1frame = PolygonVerticesToFrame[torus1];

torus2frame = PolygonVerticesToFrame[torus2];

Finally, then, with the help of smootheststep to make the animation look nice, here's the Manipulate for the animation:

smootheststep[t_] := -20 t^7 + 70 t^6 - 84 t^5 + 35 t^4;

DynamicModule[
 {gon1 = torus1frame, gon2 = torus2frame, thickness = .008, t,
  cols = RGBColor /@ {"#FF2E63", "#EAEAEA", "#08D9D6", "#252A34"}},
 Manipulate[
  t = smootheststep[s];
  Graphics3D[{Blend[cols[[;; 3]], t], Specularity[White, 40],
    Tube[Append[#, First[#]] &@(# - ConstantArray[Mean[#], Length[#]] &[
        QVertices[Complex2nToReal4n[Orthogonalize@StiefelGeodesic[gon1, gon2, t]]]]), 
     thickness]},
   ImageSize -> 540, PlotRange -> 1/3, ViewAngle -> \[Pi]/500, 
   ViewPoint -> {0, 0, 100}, ViewVertical -> {-1, 0, 0}, 
   Lighting -> "Neutral", Boxed -> False, Background -> cols[[-1]]],
  {s, 0, 1}]
 ]