Stiefel Manifolds and Polygons

This notebook is intended to be a complement to the paper “Stiefel Manifolds and Polygons” by Clayton Shonkwiler. Below are reference implementations of the algorithms discussed in the paper; indeed, this notebook was used to generate all of the figures in the paper.

The Basic Algorithm

First of all, the basic code for translating between points in the Stiefel manifold (Subscript[St, 2])((ℝ^n)) and n-gons in the plane:

StiefelSample[n_]:=Orthogonalize[RandomVariate[NormalDistribution[],{2,n}]];
FrameToEdges[frame_]:=ReIm[(Complex@@#)^2]&/@Transpose[frame];
EdgesToVertices[edges_]:=Chop[Accumulate[edges]];
FrameToVertices[frame_]:=EdgesToVertices[FrameToEdges[frame]];
EdgesToFrame[edges_]:=Transpose[ReIm[Sqrt[Complex@@#]]&/@edges];
VerticesToEdges[vertices_]:=#-RotateRight[#]&[vertices];
VerticesToFrame[vertices_]:=EdgesToFrame[VerticesToEdges[vertices]];

We can generate a random 200-gon very easily:

frame=StiefelSample[200];

And display it to the screen:

Graphics[{FaceForm[None],EdgeForm[Directive[Thickness[.01],JoinForm["Round"]]],Polygon[FrameToVertices[frame]]}]

Here is a function which will generate a specified number of n-gons and save them as PDF files in the directory /figures (which should be created beforehand):

GenerateRandomNGons[n_,numGons_]:=Module[{frames,graphic},
frames=ParallelTable[StiefelSample[n],{numGons}];
Table[
graphic=Graphics[{FaceForm[None],EdgeForm[Directive[Thickness[3/(50Sqrt[n])],JoinForm["Round"]]],Polygon[FrameToVertices[frames[[i]]]]}];
Export[NotebookDirectory[]<>"figures/"<>IntegerString[n]<>"gon"<>IntegerString[i,10,5]<>".pdf",graphic];,
{i,1,numGons}];
];

And now we can generate 20 random 200-gons very easily:

GenerateRandomNGons[200,20]

Random Triangles

The following is a simple grid of brightly-colored random triangles, created in collaboration with Anne Harding:

randomTriangleFrames=Table[StiefelSample[3],{12*17}];

triangleGrid=GraphicsGrid[Table[Graphics[{FaceForm[None],EdgeForm[Directive[AbsoluteThickness[4],JoinForm["Round"],Which[Mod[i+j,3]==0,Yellow,Mod[i+j,3]==1,Pink,Mod[i+j,3]==2,Lighter[Green]]]],Polygon[FrameToVertices[randomTriangleFrames[[12(i-1)+j]]]]}],{i,1,12},{j,1,17}],ImageSize->1080,Background->RGBColor[1/2,0,1]]

Which we can export:

Export[NotebookDirectory[]<>"figures/Triangles.pdf",triangleGrid];

We can determine whether a triangle (given as a set of edge vectors) is obtuse using the following function:

isObtuse[edges_]:=Max[π-MapThread[VectorAngle,{#,RotateLeft[#,1]}]&[edges]]>π/2;

So then we can count how many triangles in the above image are obtuse:

CountsBy[randomTriangleFrames,isObtuse[FrameToEdges[#]]&]

So almost 87%:

177/204.

Compare to the expected probability of being obtuse from our Monthly paper of

(3/2.)-((3 Log[2])/π)

And a very large polygon:

largeframe=StiefelSample[100000];

With[{cols=RGBColor/@{"#fd5f00","#76b39d","#05004e"}},
Graphics[{Line[Append[#,First[#]]&[FrameToVertices[largeframe]],VertexColors->Table[Blend[cols[[;;2]],Haversine[2π (i-1)/(Length[largeframe[[1]]]-1)]],{i,1,Length[largeframe[[1]]]}]]},ImageSize->800,Background->cols[[-1]]]
]

largePolygonGraphic=With[{cols=RGBColor/@{"#fd5f00","#76b39d","#05004e"}},
Graphics[{Thickness[.0012],JoinForm["Round"],Line[Append[#,First[#]]&[FrameToVertices[largeframe]],VertexColors->Table[Blend[cols[[;;2]],Haversine[2π (i-1)/(Length[largeframe[[1]]]-1)]],{i,1,Length[largeframe[[1]]]}]]},ImageSize->2000,Background->cols[[-1]]]
];

Export[NotebookDirectory[]<>"figures/largePolygonGraphic.png",largePolygonGraphic];

Convex n-Gons

The following functions will permute the edges of any planar n-gon to produce a convex n-gon:

RectifyAngle[θ_] := If[θ<0,θ+2π,θ];
PositiveArg[z_] := RectifyAngle @ Arg[z];
Convexify[edges_]:= RotateRight[SortBy[edges,PositiveArg[Complex@@#]&],RandomInteger[{1,Length[edges]}]];
FrameToConvexEdges[frame_]:=Convexify[FrameToEdges[frame]];
FrameToConvexVertices[frame_]:=EdgesToVertices[FrameToConvexEdges[frame]];

As an example, the random 200-gon we generated above becomes the following convex 200-gon when run through this function:

convex200gon=Graphics[{FaceForm[None],EdgeForm[Directive[Thickness[.05],JoinForm["Round"]]],Polygon[FrameToConvexVertices[frame]]}]

Export[NotebookDirectory[]<>"figures/convex200gon.pdf",convex200gon];

Again, we can produce a large number of convex n-gons with the following function:

GenerateRandomConvexNGons[n_,numGons_]:=Module[{frames,graphic},
frames=ParallelTable[StiefelSample[n],{numGons}];
Table[
graphic=Graphics[{FaceForm[None],EdgeForm[Directive[Thickness[3/(50Sqrt[n])],JoinForm["Round"]]],Polygon[FrameToConvexVertices[frames[[i]]]]}];
Export[NotebookDirectory[]<>"figures/convex"<>IntegerString[n]<>"gon"<>IntegerString[i,10,5]<>".pdf",graphic];,
{i,1,numGons}];
];

GenerateRandomConvexNGons[12,3]

Now, we can produce the nice Riemann mapping theorem picture as follows: first, after playing around a bit, the following point in (Subscript[St, 2])((ℝ^5)) leads to a nice pentagon:

frame5={{1.2247163790775448`,0.8132217749967826`,0.6704157895663052`,0.4971399951917043`,1.0074165700767115`},{0.02552484449081336`,1.298852168127096`,0.4454774408827977`,-1.1569644200177216`,-0.805026846780099`}};

Graphics[{FaceForm[None],EdgeForm[Thickness[.006]],Polygon[FrameToConvexVertices[frame5]]}]

So now we’re going to mean-center and extract the vertices as complex numbers:

Complex@@#&/@(-Table[Mean[#],{Length[#]}]+# &[FrameToConvexVertices[frame5]])

Now, plug these vertices into the Schwarz–Christoffel Toolbox for Matlab, which tells us the parameters for the Schwarz–Christoffel transformation taking the unit disk to this convex polygon. Here’s what to type in MATLAB (to automate this, we could use MATLink):

First, define the polygon p:

>> p = polygon([-1.08829+0.598657i -0.721469-1.02334i 0.77781-0.960816i 1.02882-0.363506i 0.00312993+1.749i])

Then compute the map f:

>> f = diskmap(p);

To make things look nicer, it’s better if the map f sends the origin to the origin, so we re-center:

>> f = center(f,0);

Now we can see the angles (in the alpha column) and the pre-vertices by evaluating:

>> f

So now we import these back into Mathematica:

angles={0.6876,0.5575,0.6399,0.7295,0.3855};
prevertices={0.7773+0.6292I,-0.6559+0.7548I,-0.9318-0.3631I,-0.7475-0.6643I,1.0000+0.0000I};

We can use this information to define our Schwarz–Christoffel map (which we’ll define as a compiled function to make it run a lot faster):

sc=Compile[{{z,_Complex}},
NIntegrate[Exp[Sum[((angles[[k]]-1)(Log[I prevertices[[k]]]+Log[(ζ-prevertices[[k]])/(I prevertices[[k]])])),{k,1,5}]],{ζ,0,z}]]

For some reason the above map sends the unit disk to a rotated version of our original polygon. Rather than worry about it, I just rotated back to produce the following image (that’s what the E^(-Iπ/2.3) is doing):

fancypentagon=Show[ParallelTable[Quiet@ParametricPlot[ReIm[E^(-I π/2.3)sc[r E^(I t)]],{t,0.,2π},PlotRange->{{-1.8,1.8},{-1.5,2.1}},ImageSize->540,Axes->False,PlotStyle->Directive[Thickness[.005],ColorData["BrownCyanTones"][r/2+1/2]],Background->GrayLevel[.2],PlotPoints->200],{r,0.,1,1/20}]]

Export[NotebookDirectory[]<>"figures/fancypentagon.png",fancypentagon];

Quad Tilings

Every quadrilateral tiles the plane. This follows from the fact that the sum of the interior angles of any quadrilateral is (360^◦), so four copies of the quadrilateral can share a vertex, and then we can tile the plane by taking appropriate lattice translations.

The following code makes the basic blocks consisting of four quadrilaterals sharing a vertex and then builds the tiling from the blocks:

BasicBlock[edges_]:=Module[{n},
n=Length[edges];
Table[Accumulate[(-1)^i RotateRight[edges,i]],{i,0,n-1}]
];
QuadTilingBasis[edges_]:={edges[[2]]-edges[[4]],edges[[3]]-edges[[1]]};
EdgesToTessellation[edges_,numx_,numy_]:=Module[{basis},
basis=QuadTilingBasis[edges];
Table[Polygon/@(x*basis[[1]]+y*basis[[2]]+#&/@#&/@BasicBlock[edges]),{x,-(numx-1)/2,(numx-1)/2},{y,-(numy-1)/2,(numy-1)/2}]
];

Since we’re representing quadrilaterals by points in the Stiefel manifold (Subscript[St, 2])((ℝ^4)), it’s handy to have a function transforming orthonormal 2-frames in (ℝ^4) into a tiling:

QuadFrameToTessellation[frame_,numx_,numy_]:=EdgesToTessellation[FrameToEdges[frame],numx,numy];

After some playing around, the following Stiefel manifold points produce pleasing tilings:

frame1={{-0.21728090898823724`,-0.08834582559466465`,0.7450486536045761`,0.6244089408802944`},{-0.7224486073877968`,0.44620730863225605`,-0.4229311116667767`,0.3163800281163406`}};
frame2={{0.42870301121644105`,0.38283874197110207`,-0.633804812205533`,-0.5176289074665839`},{-0.6547467745123633`,0.7145404849942669`,0.14903954331161623`,-0.19627982855932777`}};
frame3=frame={{0.7524477450891495`,0.20416350240966882`,0.6031470499258688`,0.16838435604233948`},{0.147256911052663`,0.37436399392169745`,-0.5206619855607935`,0.7530458810684216`}};

tiling1=With[{frame=frame1,cols=RGBColor/@{"#00C3FF","#F46188","#55E9BC","#3B475E"}},Graphics[{FaceForm[None],EdgeForm[Directive[JoinForm["Round"],Thickness[.007],cols[[1]]]],QuadFrameToTessellation[frame,5,5]},PlotRange->1,ImageSize->540,Background->Last[cols]]]

Export[NotebookDirectory[]<>"figures/quadtiling1.pdf",tiling1];

tiling2=With[{frame=frame2,cols=RGBColor/@{"#00C3FF","#F46188","#55E9BC","#3B475E"}},Graphics[{FaceForm[None],EdgeForm[Directive[JoinForm["Round"],Thickness[.007],cols[[2]]]],QuadFrameToTessellation[frame,5,5]},PlotRange->1,ImageSize->540,Background->Last[cols]]]

Export[NotebookDirectory[]<>"figures/quadtiling2.pdf",tiling2];

tiling3=With[{frame=frame3,cols=RGBColor/@{"#00C3FF","#F46188","#55E9BC","#3B475E"}},Graphics[{FaceForm[None],EdgeForm[Directive[JoinForm["Round"],Thickness[.007],cols[[3]]]],QuadFrameToTessellation[frame,7,7]},PlotRange->1,ImageSize->540,Background->Last[cols]]]

Export[NotebookDirectory[]<>"figures/quadtiling3.pdf",tiling3];

A to Z

Now, we want to know how to morph one polygon into another. Since n-gons are represented by pairs of orthonormal vectors ((OverVector[x]),(OverVector[y])), i.e., points in the Stiefel manifold (Subscript[St, 2])((ℝ^n)), we can get nice paths simply by  directly rotating (Subscript[(OverVector[x]), 1]) to (Subscript[(OverVector[x]), 2]), and directly rotating (Subscript[(OverVector[y]), 1]) to (Subscript[(OverVector[y]), 2]) in the orthogonal complement of the (Subscript[(OverVector[x]), 1])(Subscript[(OverVector[x]), 2])-plane, which we accomplish as follows:

StiefelPath[{x1_,y1_},{x2_,y2_},t_]:=Module[{θ1,θ2},
θ1=VectorAngle[x1,x2];
θ2=VectorAngle[y1,y2];
Orthogonalize[{RotationTransform[t*θ1,{x1,x2}][x1],RotationTransform[t*θ2,{y1,y2}][y1]}]
];

I want to illustrate this by giving the geodesic between a polygonal letter A and a polygonal letter Z. To get polygonal versions of the letters, I use the following toy which allows you to place the vertices of a polygon and recover their coordinates:

Manipulate[Grid[{{Graphics[
{{EdgeForm[{Thick}],FaceForm[None],Polygon[u]},{{PointSize[0.025],Point[#]}& /@ u}},PlotRange->{{-9.5,12.5},{-9.5,12.5}},
GridLines->{Range[-15,15],Range[-15,15]},GridLinesStyle->Directive[Opacity[0.2]],ImageSize->{600,600}]},{u}}],Button["Copy coordinates to clipboard",CopyToClipboard[u]],
{{u,{{1,1},{1,2},{2,2},{2,1}}},Locator,LocatorAutoCreate->True}]

Copying and pasting gives me the vertices of my A and Z:

RawAData={{-3.,.9},{1.07,8.18},{1.85,8.2},{5.96,0.86},{4.89,0.86},{1.77,6.6},{1.18,6.6},{.12,4.58},{2.87,4.6},{3.35,3.76},{-.18,3.76},{-1.72,.94}};

RawZData={{0.82,7.26},{7.61,7.32},{7.65,6.04},{1.92,-.06},{7.76,-.02},{7.8,-1.16},{.71,-1.12},{.71,.02},{6.36,6.12},{.85,6.12},{.82,6.12},{.82,7.26}};

Which I can convert to points in the Stiefel manifold:

AFrame=VerticesToFrame[RawAData];
ZFrame=VerticesToFrame[RawZData];

The smootheststep function produces visually pleasing movement:

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

Now, we build a Manipulate showing the path between A and Z,

DynamicModule[{cols=RGBColor/@{"#f77e5e","#3dbd5d","#303030"},centeredpoints,s},
Manipulate[
s=smootheststep[t];
centeredpoints=#-Table[Mean[#],{Length[#]}]&[FrameToVertices[StiefelPath[AFrame,ZFrame,s]]];
Graphics[{FaceForm[None],EdgeForm[Directive[Thickness[.008],JoinForm["Round"],Blend[cols[[;;2]],s]]],Polygon[centeredpoints]},ImageSize->{540,540},PlotRange->{{-.3+.02s,.3+.02s},{-.28-.04s,.32-.04s}},Background->cols[[-1]]],{t,0,1}]
]

And visualize the entire geodesic in a single graphic by translating appropriately:

Module[{cols=RGBColor/@{"#f77e5e","#3dbd5d","#303030"},centeredpoints,s},
azStills=Table[s=smootheststep[t];
centeredpoints=#+ ConstantArray[{.5s,0+.04s}-Mean[#],Length[#]]&[FrameToVertices[StiefelPath[AFrame,ZFrame,s]]];
Graphics[{FaceForm[None],EdgeForm[Directive[Thickness[.008],Opacity[.5],JoinForm["Round"],Blend[cols[[;;2]],s]]],Polygon[centeredpoints]},Background->cols[[-1]]],{t,0,1,.005}]];

AtoZstill=Show[Reverse@azStills,PlotRange->{{-.29,.75},{-.3,.33}},ImageSize->{800,528}]

Export[NotebookDirectory[]<>"/figures/azStill.png",AtoZstill];

And here’s the corresponding geodesic in the Grassmannian:

GrassmannianGeodesic[frame1_,frame2_,t_]:=Block[{k=Length[frame1],U,Σ,V,newframe1,newframe2,dists},
{U,Σ,V}=SingularValueDecomposition[frame1.Transpose[frame2]];
newframe1=Transpose[U].frame1;
newframe2=V.frame2;
dists=ArcCos/@Diagonal[Σ];
Table[RotationTransform[t*dists[[i]],{newframe1[[i]],newframe2[[i]]}][newframe1[[i]]],{i,1,k}]
];

Module[{cols=RGBColor/@{"#f77e5e","#3dbd5d","#303030"},centeredpoints,s},
azStills2=Table[s=smootheststep[t];
centeredpoints=#+ ConstantArray[{.5s,.04s}-Mean[#],Length[#]]&[FrameToVertices[GrassmannianGeodesic[Normalize/@AFrame,Normalize/@ZFrame,s]]];
Graphics[{FaceForm[None],EdgeForm[Directive[Thickness[.008],Opacity[.5],JoinForm["Round"],Blend[cols[[;;2]],s]]],Polygon[centeredpoints]},Background->cols[[-1]]],{t,0,1,.005}]];

AtoZstill2=Show[Reverse@azStills2,PlotRange->{{-.29,.75},{-.3,.33}},ImageSize->{800,528}]

Export[NotebookDirectory[]<>"/figures/azStill2.png",AtoZstill2];

Instead of just taking the positive square root of each edge to turn a polygon into a point in the Stiefel manifold, we can take either square root, leading to a total of (2^n) possible lifts of any n-gon. 

For example, we can take the positive lift of the regular pentagonal star:

star=Table[{Cos[θ],Sin[θ]},{θ,π/2,4π+π/2-4π/5,4π/5}];

starframe1=VerticesToFrame[star];

And then all possible lifts:

starframe2=Flatten[Table[{i,j,k,l,m}*#&/@starframe1,{i,{-1,1}},{j,{-1,1}},{k,{-1,1}},{l,{-1,1}},{m,{-1,1}}],4];

...and, after some experimentation, the 28th lift produces a nice path from the star to itself (here we’re importing one of Peter Kovesi’s very nice isoluminant color maps):

smoothstep[x_]:=3x^2-2x^3;
IsoCols=Import["https://peterkovesi.com/projects/colourmaps/CETperceptual_csv_0_1.zip","CETperceptual_csv_0_1/CET-I2.csv"];

Block[{centeredpoints,s,cols=RGBColor/@IsoCols},
starToStarStills=ParallelTable[s=smoothstep[t];
centeredpoints=#+ ConstantArray[{3/4t,0}-Mean[#],Length[#]]&[FrameToVertices[Quiet@GrassmannianGeodesic[Normalize/@N[starframe1],Normalize/@N[starframe2[[28]]],s]].RotationMatrix[π]];
Graphics[{FaceForm[None],EdgeForm[Directive[Thickness[.004],JoinForm["Round"],Opacity[.8],Blend[cols,t]]],Polygon[centeredpoints]}],{t,0.,1,1/60}];
]

starPath=Show[Reverse@starToStarStills,ImageSize->800,ImagePadding->60,Background->RGBColor["#142d4c"]]

Export[NotebookDirectory[]<>"figures/starPath.png",starPath];

Knot Morph

The following gives code for translating between polygons in space and points in the Stiefel manifold (Subscript[St, 2])((ℂ^n)) of Hermitian-orthonormal 2-frames in (ℂ^n) (in fact, the Stiefel manifold more naturally corresponds to framed polygons, meaning each edge comes equipped with a choice of two normal vectors, so in the below code we pass from a polygon to a Stiefel manifold point by choosing a particular framing on the polygon: in this case, the discrete Bishop frame):

<<Quaternions`
QuaternionEdges[StiefelForm_]:=Conjugate[#]**Quaternion[0,1,0,0]**#&/@(Quaternion@@#&/@Transpose[StiefelForm]);
QEdges[StiefelForm_] :=List@@#[[2;;]]&/@QuaternionEdges[StiefelForm];
QVertices[StiefelForm_]:=Accumulate[QEdges[StiefelForm]];
Real4nToComplex2n[M_]:=Transpose[Complex@@#&/@#&/@(Partition[#,2]&/@Transpose[M])];
Complex2nToReal4n[M_]:=Transpose[Flatten/@(ReIm/@Transpose[M])];
ComplexStiefelPath[{A_,B_},{C_,D_},t_]:=Module[{cPerp,dPerp,dist1,dist2},
{cPerp,dPerp}={Normalize[C-Re[Conjugate[C].A]*A],Normalize[D-Re[Conjugate[D].B]*B]};
dist1=ArcCos[Re[Conjugate[A].C]];
dist2=ArcCos[Re[Conjugate[B].D]];
Orthogonalize[{Cos[t*dist1]*A+Sin[t*dist1]*cPerp,Cos[t*dist2]*B+Sin[t*dist2]*dPerp}]
];
PolygonToFramedPolygon[edges_]:=Module[{normals,binormals},
binormals=Table[Norm[edges[[i]]]Normalize[Cross[edges[[i]],RotateLeft[edges][[i]]]],{i,1,Length[edges]}];
normals=Table[Norm[edges[[i]]]Normalize[Cross[binormals[[i]],edges[[i]]]],{i,1,Length[edges]}];
Table[{edges[[i]],normals[[i]],binormals[[i]]},{i,1,Length[edges]}]
];

This function for transforming a (scaled) rotation matrix to a quaternion is a slight modification of rotm2q from Reichlin’s quaternion package for Octave.

MatrixToQuaternion[A_]:=Module[{M,T,l,s,w,x,y,z},
l=Norm[A[[1]]];
M=1/l A;
T=Tr[M];
Which[T>0,
s=1/(2Sqrt[T+1]);
w=1/(4s);
x=s(M[[2,3]]-M[[3,2]]);
y=s(M[[3,1]]-M[[1,3]]);
z=s(M[[1,2]]-M[[2,1]]);,
M[[1,1]]>M[[2,2]]&&M[[1,1]]>M[[3,3]],
s=2Sqrt[1+M[[1,1]]-M[[2,2]]-M[[3,3]]];
w=(M[[2,3]]-M[[3,2]])/s;
x=s/4;
y=(M[[1,2]]+M[[2,1]])/s;
z=(M[[1,3]]+M[[3,1]])/s;,
M[[2,2]]>M[[3,3]],
s=2Sqrt[1+M[[2,2]]-M[[1,1]]-M[[3,3]]];
w=(M[[3,1]]-M[[1,3]])/s;
x=(M[[1,2]]+M[[2,1]])/s;
y=s/4;
z=(M[[2,3]]+M[[3,2]])/s;,
True,
s=2Sqrt[1+M[[3,3]]-M[[1,1]]-M[[2,2]]];
w=(M[[1,2]]-M[[2,1]])/s;
x=(M[[1,3]]+M[[3,1]])/s;
y=(M[[2,3]]+M[[3,2]])/s;
z=s/4;
];
Sqrt[l]{w,x,y,z}
];
PolygonVerticesToFrame[verts_]:=Orthogonalize@Real4nToComplex2n[Transpose[MatrixToQuaternion[Transpose[#]]&/@PolygonToFramedPolygon[VerticesToEdges[verts]]]];

Now, define a (3,4)-torus knot:

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

torus34=Table[Stereo3D[1/Sqrt[2]RotationMatrix[π/8,{{1,0,0,0},{0,1,0,0}}].{Cos[3θ],Sin[3θ],Cos[4θ],Sin[4θ]}],{θ,0.,2π,2π/1000}];

And a (4,3)-torus knot:

torus43=Table[Stereo3D[1/Sqrt[2]RotationMatrix[-π/72,{{1,0,0,0},{0,1,0,0}}].{Cos[4θ],Sin[4θ],Cos[3θ],Sin[3θ]}],{θ,0.,2π,2π/1000}];

And build their corresponding frames in the Stiefel manifold:

torus34frame=PolygonVerticesToFrame[torus34];
torus43frame=PolygonVerticesToFrame[torus43];

Here’s a Manipulate showing the path between them:

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

And we can export a few frames:

morph=Module[{gon1=torus34frame,gon2=torus43frame,thickness=.008,goodbasis,cols=RGBColor/@{"#FF2E63","#EAEAEA","#08D9D6","#252A34"}},
ParallelTable[
Graphics3D[{Blend[cols[[;;3]],t],Specularity[White,40],Tube[Append[#,First[#]]&@(#-ConstantArray[Mean[#],Length[#]]&[QVertices[Complex2nToReal4n[ComplexStiefelPath[gon1,gon2,t]]]]),thickness]},ImageSize->540,PlotRange->1/3,ViewAngle->π/500,ViewPoint->{0,0,100},ViewVertical->{-1,0,0},Lighting->"Neutral",Boxed->False,Background->cols[[-1]]],{t,0.,1,1/4},"DistributedContexts"->All]
];

Table[Export[NotebookDirectory[]<>"figures/morph"<>IntegerString[i]<>".png",morph[[i]]];,{i,1,Length[morph]}];

Now, we find paths between different framings on the same polygon. First of all, here are two framings on the regular 200-gon:

circleframe=With[{n=200},
1/10 Real4nToComplex2n[Transpose[MatrixToQuaternion/@(Table[Transpose@Prepend[{{-#[[i,2]],#[[i,1]],0},{0,0,1}},#[[i]]],{i,1,n}]&[Append[#,0]&/@N[CirclePoints[n]]])]]
];
circleframe2=With[{n=200},1/10 Real4nToComplex2n[Transpose[MatrixToQuaternion/@(Table[Transpose@Prepend[RotationMatrix[6π i/n].{{-#[[i,2]],#[[i,1]],0},{0,0,1}},#[[i]]],{i,1,n}]&[Append[#,0]&/@N[CirclePoints[n]]])]]
];

And some code for computing the minimizing geodesic between two points in the Grassmannian:

ComplexProjectionBasis[{A_,B_},{C_,D_}]:=Eigenvectors[Transpose[{A,B}].Conjugate[{A,B}].Transpose[{C,D}].Conjugate[{C,D}],2];

ComplexGrassmannianGeodesic[{A_,B_},{C_,D_},t_]:=Module[{a,b,c,d,cPerp,dPerp,dist1,dist2},
{a,b}=Orthogonalize[ComplexProjectionBasis[{C,D},{A,B}]];
{c,d}=Orthogonalize[ComplexProjectionBasis[{A,B},{C,D}]];
{cPerp,dPerp}={Normalize[c-Re[Conjugate[c].a]*a],Normalize[d-Re[Conjugate[d].b]*b]};
dist1=ArcCos[Re[Conjugate[c].a]];
dist2=ArcCos[Re[Conjugate[d].b]];
{Cos[t*dist1]*a+Sin[t*dist1]*cPerp,Cos[t*dist2]*b+Sin[t*dist2]*dPerp}
];

Also, we’ll need to record the full framing on each edge, not just the edge vector:

QuaternionFramedEdges[StiefelForm_]:={Conjugate[#]**Quaternion[0,1,0,0]**#,Conjugate[#]**Quaternion[0,0,1,0]**#,Conjugate[#]**Quaternion[0,0,0,1]**#}&/@(Quaternion@@#&/@Transpose[StiefelForm]);
QFramedEdges[StiefelForm_] :=List@@#[[2;;]]&/@#&/@QuaternionFramedEdges[StiefelForm];

Now, we build a parametrized surface which twists with the framing in order to visualize the framing:

Block[{frame=Transpose[QFramedEdges[Complex2nToReal4n[circleframe2]]],f,g,h},
f=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[#-ConstantArray[Mean[#],Length[#]]&[Accumulate[frame[[1]]]]]];
g=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[frame[[2]]]];
h=
Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[frame[[3]]]];
ParametricPlot3D[f[t]+5(Cos[θ]g[t]+Sin[θ]h[t]),{t,0,1},{θ,0,2π},MeshFunctions->{#5&},Mesh->5,MeshShading->{Red, Yellow},MeshStyle->None,PlotPoints->50]
]

Here's a geodesic between these two frames:

DynamicModule[{frame,n,f,g,h,cols=RGBColor/@{"#f6490d","#1dced8","#161c2e"}},
Manipulate[
frame=Transpose[QFramedEdges[Complex2nToReal4n[ComplexGrassmannianGeodesic[circleframe,circleframe2,t]]]];
f=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[#-ConstantArray[Mean[#],Length[#]]&[Accumulate[frame[[1]]]]]];
g=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[frame[[2]]]];
h=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[frame[[3]]]];
n=Length[frame[[1]]];
ParametricPlot3D[f[t]+3(Cos[θ]g[t]+Sin[θ]h[t]),{t,0,1},{θ,0,2π},MeshFunctions->{#5&},Mesh->5,MeshShading->{cols[[1]], cols[[2]]},MeshStyle->None,PlotPoints->50,Boxed->False,ViewPoint->-100Normalize[Cross[f[0.],f[.25]]],ViewVertical->f[0.],ViewAngle->π/300,PlotRange->1/2,ImageSize->{800,800},Axes->False,Lighting->"Neutral",Background->cols[[-1]]],
{t,0.,1}]
]

Which we can export as an animated GIF:

sausage3=Block[{t,frame,n,f,g,h,cols=RGBColor/@{"#f6490d","#1dced8","#161c2e"}},
Table[
t=smootheststep[s];
frame=Transpose[QFramedEdges[Complex2nToReal4n[ComplexGrassmannianGeodesic[circleframe,circleframe2,t]]]];
f=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[#-ConstantArray[Mean[#],Length[#]]&[Accumulate[frame[[1]]]]]];
g=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[frame[[2]]]];
h=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[frame[[3]]]];
n=Length[frame[[1]]];
ParametricPlot3D[f[t]+3(Cos[θ]g[t]+Sin[θ]h[t]),{t,0,1},{θ,0,2π},MeshFunctions->{#5&},Mesh->5,MeshShading->{cols[[1]], cols[[2]]},MeshStyle->None,PlotPoints->200,Boxed->False,ViewPoint->-100Normalize[Cross[f[0.],f[.25]]],ViewVertical->f[0.],ViewAngle->π/300,PlotRange->1/2,ImageSize->{800,800},Axes->False,Lighting->"Neutral",Background->cols[[-1]]],
{s,0.,1,#}]&[1/50]
];

Export[NotebookDirectory[]<>"sausage3.gif",Join[#,Drop[Reverse[Drop[#,1]],1]]&[sausage3],"AnimationRepetitions"->Infinity,"DisplayDurations"->1/50];

Or we can combine a few points along the way into a single image:

frameGeodesic=Block[{t,frame,n,f,g,h,cols=RGBColor/@{"#f6490d","#1dced8","#161c2e"}},
Show[Table[frame=Transpose[QFramedEdges[Complex2nToReal4n[ComplexGrassmannianGeodesic[circleframe,circleframe2,s]]]];
f=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[#-ConstantArray[Mean[#],Length[#]]&[Accumulate[frame[[1]]]]]];
g=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[frame[[2]]]];
h=Interpolation[Transpose[{N@Range[0,1,1/(Length[#]-1)],#}]]&[Prepend[#,#[[-1]]]&[frame[[3]]]];
n=Length[frame[[1]]];
ParametricPlot3D[
{0,0,-s}+RotationMatrix[{Cross[f[0.],f[.25]],{0,0,1}}].(f[t]+3(Cos[θ]g[t]+Sin[θ]h[t])),{t,0,1},{θ,0,2π},MeshFunctions->{#5&},Mesh->5,MeshShading->{cols[[1]], cols[[2]]},MeshStyle->None,PlotPoints->200]
,{s,0.,1,1/6}],PlotRange->{{-1/2,1/2},{-1/2,1/2},Automatic},Boxed->False,ViewPoint->100{-.9,0,.1},ImageSize->800,ViewVertical->{0,0,1},Axes->False,Lighting->"Neutral",Background->cols[[-1]]]
]

Export[NotebookDirectory[]<>"figures/frameGeodesic.png",frameGeodesic];