Frame Field Initialization

Functional Representations of Frames in 2D

Paper at a Glance

NoteReference

On Smooth Frame Field Design — Ray, Sokolov, Reberol (2015)

Download the original paper (PDF)

This page walks through Section 2 of the paper: Functional Representation of Frames in 2D. Each subsection below mirrors a subsection of the paper, with direct quotes followed by interactive toy examples.

View as interactive slides

---

2.1 Problem Settings

“Given a 2D shape, frame field design in 2D consists of finding a smooth frame field aligned with the boundary of the shape. We formulate it as minimizing the field curvature, based on the following definitions.”

The paper defines the key objects:

• “A frame is a set of 4 unit vectors \(f = \{f_i\}, i \in [0,3]\) that is invariant by a rotation of \(\pi/2\). It can be represented by the angle \(\theta\) such that \(\forall i\), \(f_i = (\cos(\theta + i\pi/2), \sin(\theta + i\pi/2))\).”
• “A frame field is a frame per vertex of a 2D shape triangulation.”
• “The boundary constraint: a frame located on a boundary vertex must have one of its member vectors equal to the normal on the boundary.”
• “The rotation angle between two frames is the angle \(\Delta\theta\) of the rotation that transforms one frame into the other. This angle being defined modulo \(\pi/2\), we choose the \(\Delta\theta\) with minimum absolute value.”
• “The curvature of a frame field is the sum over each edge of the squared rotation angle between frames that are defined at the edge extremities.”
• “A triangle is said to be singular if the sum of the rotation angles over his boundary is not equal to 0.”

Drag the slider to rotate the frame — four unit vectors separated by \(\pi/2\):

viewof theta_deg = Inputs.range([0, 90], {value: 30, step: 1, label: "θ (degrees)"})

Figure 1

{
  const w = 380, h = 380, cx = w/2, cy = h/2, r = 140;
  const theta = theta_deg * Math.PI / 180;

  const svg = d3.create("svg").attr("viewBox", `0 0 ${w} ${h}`).style("width", "100%").style("max-width", w + "px");

  // Unit circle
  const circle_pts = d3.range(0, 2*Math.PI + 0.01, 0.02).map(a => [cx + r*Math.cos(a), cy - r*Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circle_pts)).attr("fill", "none").attr("stroke", "#ccc").attr("stroke-width", 1);

  // Axes
  svg.append("line").attr("x1", cx - r - 20).attr("y1", cy).attr("x2", cx + r + 20).attr("y2", cy).attr("stroke", "#ddd");
  svg.append("line").attr("x1", cx).attr("y1", cy - r - 20).attr("x2", cx).attr("y2", cy + r + 20).attr("stroke", "#ddd");

  const colors = ["steelblue", "firebrick", "steelblue", "firebrick"];

  for (let i = 0; i < 4; i++) {
    const angle = theta + i * Math.PI / 2;
    const dx = r * Math.cos(angle), dy = -r * Math.sin(angle);

    // Arrow line
    svg.append("line")
      .attr("x1", cx).attr("y1", cy)
      .attr("x2", cx + dx).attr("y2", cy + dy)
      .attr("stroke", colors[i]).attr("stroke-width", 3);

    // Arrowhead
    const headLen = 12, headAngle = 0.4;
    const a1 = Math.atan2(-dy, dx) + Math.PI + headAngle;
    const a2 = Math.atan2(-dy, dx) + Math.PI - headAngle;
    svg.append("polygon")
      .attr("points", `${cx+dx},${cy+dy} ${cx+dx+headLen*Math.cos(a1)},${cy+dy-headLen*Math.sin(a1)} ${cx+dx+headLen*Math.cos(a2)},${cy+dy-headLen*Math.sin(a2)}`)
      .attr("fill", colors[i]);

    // Label
    svg.append("text")
      .attr("x", cx + 1.15*dx).attr("y", cy + 1.15*dy)
      .attr("text-anchor", "middle").attr("dominant-baseline", "middle")
      .attr("font-size", 16).attr("fill", colors[i])
      .text(`f${i}`);
  }

  // Angle arc
  if (theta > 0.01) {
    const arc_pts = d3.range(0, theta + 0.01, 0.02).map(a => [cx + 40*Math.cos(a), cy - 40*Math.sin(a)]);
    svg.append("path").attr("d", d3.line()(arc_pts)).attr("fill", "none").attr("stroke", "orange").attr("stroke-width", 2);
    svg.append("text")
      .attr("x", cx + 55*Math.cos(theta/2)).attr("y", cy - 55*Math.sin(theta/2))
      .attr("text-anchor", "middle").attr("font-size", 16).attr("fill", "orange")
      .text("θ");
  }

  // Title
  svg.append("text").attr("x", cx).attr("y", 20).attr("text-anchor", "middle").attr("font-size", 16).attr("fill", "#333")
    .text(`2D Frame (θ = ${theta_deg}°)`);

  return svg.node();
}

---

2.2 Functional Approach: Frame Representation

The key insight of the paper is to represent frames not by angles, but by functions.

“The frame aligned with coordinate axes is called the reference frame \(\tilde{f} = \{(1,0),(0,1),(-1,0),(0,-1)\}\). Instead of using the angle approach, we represent it by the function \(\tilde{F}(\alpha) = \cos(4\alpha)\) with \(\alpha \in [0, 2\pi[\), that exhibits the same \(\pi/2\) rotation invariance as the frame.”

“Any other frame \(f\) can be obtained by a rotation of \(\tilde{f}\) by angle \(\theta\). The functional counterpart is to rotate the graph of the function \(\tilde{F}\), namely any frame \(f\) can be represented by a function \(F(\alpha) = \tilde{F}(\alpha - \theta) = \cos(4(\alpha - \theta))\) with \(\alpha \in [0, 2\pi[\).”

“A compact representation of these functions is given by the trigonometric relation \(f(\alpha) = \cos(4\alpha - 4\theta) = \cos(4\theta)\cos(4\alpha) + \sin(4\theta)\sin(4\alpha)\): we see that a frame function \(F\) can be represented by a 2D vector of coefficients \(a = (a_0, a_1)^\top = (\cos(4\theta), \sin(4\theta))^\top\) in Fourier basis \(B = (\cos(4\alpha), \sin(4\alpha))\) i.e. \(F = Ba\).”

TipWhy \(F = Ba\) — unpacking the Fourier basis step by step

The key move is just the angle addition formula for cosine. Start from the rotated frame function:

\[F(\alpha) = \cos(4(\alpha - \theta)) = \cos(4\alpha - 4\theta)\]

Apply the identity \(\cos(A - B) = \cos A \cos B + \sin A \sin B\) with \(A = 4\alpha\) and \(B = 4\theta\):

\[F(\alpha) = \underbrace{\cos(4\theta)}_{a_0} \cdot \underbrace{\cos(4\alpha)}_{B_0(\alpha)} + \underbrace{\sin(4\theta)}_{a_1} \cdot \underbrace{\sin(4\alpha)}_{B_1(\alpha)}\]

Now notice the structure: \(F(\alpha)\) is a linear combination of two fixed basis functions \(B_0(\alpha) = \cos(4\alpha)\) and \(B_1(\alpha) = \sin(4\alpha)\), weighted by coefficients that depend only on \(\theta\). Written as a dot product:

\[F(\alpha) = \begin{pmatrix} \cos(4\alpha) & \sin(4\alpha) \end{pmatrix} \begin{pmatrix} \cos(4\theta) \\ \sin(4\theta) \end{pmatrix} = B(\alpha) \cdot a\]

This is exactly a Fourier series truncated to a single frequency (the 4th harmonic). The functions \(\cos(4\alpha)\) and \(\sin(4\alpha)\) form an orthogonal basis for all signals with period \(\pi/2\) — which is precisely the symmetry of a frame. Any frame function, regardless of \(\theta\), lives in this 2D function space, and the coefficients \(a = (\cos 4\theta, \sin 4\theta)^\top\) are its coordinates in that basis.

Why frequency 4? Because a frame has 4-fold symmetry (\(\pi/2\) periodicity). The function \(F(\alpha)\) must satisfy \(F(\alpha + \pi/2) = F(\alpha)\), so only Fourier modes with period dividing \(\pi/2\) survive — that means frequencies \(4, 8, 12, \ldots\) Since we’re representing a single “bump” per frame vector, frequency 4 alone is sufficient.

“A coefficient vector \(a\) is feasible if and only if there exists \(\theta\) such that \(a = (\cos(4\theta), \sin(4\theta))^\top\). Geometrically, \(a\) is constrained to live on a curve parameterized by \(\theta\). This curve represents, in coefficient space, all possible rotations of the reference frame. In 2D, this curve is the unit circle, so the constraint on \(a\) is simply: \(a^\top a = 1\).”

Drag θ to see how the frame, its function representation, and its coefficient vector are all linked:

viewof theta2_deg = Inputs.range([0, 90], {value: 30, step: 1, label: "θ (degrees)"})

Figure 2

{
  const theta = theta2_deg * Math.PI / 180;
  const W = 750, H = 280;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  // --- Panel 1: Frame vectors ---
  const p1x = 5, p1w = 200, p1cx = p1x + p1w/2, p1cy = H/2, p1r = 80;

  svg.append("text").attr("x", p1cx).attr("y", 18).attr("text-anchor", "middle").attr("font-size", 13).attr("fill", "#333").text("Frame in physical space");

  // circle
  const circ = d3.range(0, 2*Math.PI+0.01, 0.02).map(a => [p1cx+p1r*Math.cos(a), p1cy-p1r*Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circ)).attr("fill", "none").attr("stroke", "#ddd");

  const frameColors = ["steelblue","firebrick","steelblue","firebrick"];
  for (let i = 0; i < 4; i++) {
    const a = theta + i*Math.PI/2;
    svg.append("line").attr("x1",p1cx).attr("y1",p1cy).attr("x2",p1cx+p1r*Math.cos(a)).attr("y2",p1cy-p1r*Math.sin(a))
      .attr("stroke",frameColors[i]).attr("stroke-width",2.5);
  }

  // --- Panel 2: Function F(α) ---
  const p2x = 220, p2w = 280, p2h = 220, p2top = 40;
  const plotL = p2x+40, plotR = p2x+p2w-10, plotT = p2top+10, plotB = p2top+p2h-20;
  const xScale = d3.scaleLinear().domain([0, 2*Math.PI]).range([plotL, plotR]);
  const yScale = d3.scaleLinear().domain([-1.2, 1.2]).range([plotB, plotT]);

  svg.append("text").attr("x",(plotL+plotR)/2).attr("y",18).attr("text-anchor","middle").attr("font-size",13).attr("fill","#333")
    .text(`F(α) = cos(4(α − ${theta2_deg}°))`);

  // axes
  svg.append("line").attr("x1",plotL).attr("y1",yScale(0)).attr("x2",plotR).attr("y2",yScale(0)).attr("stroke","#ccc");
  svg.append("line").attr("x1",plotL).attr("y1",plotT).attr("x2",plotL).attr("y2",plotB).attr("stroke","#ccc");

  // reference function (gray)
  const refPts = d3.range(0, 2*Math.PI, 0.02).map(a => [xScale(a), yScale(Math.cos(4*a))]);
  svg.append("path").attr("d", d3.line()(refPts)).attr("fill","none").attr("stroke","#bbb").attr("stroke-width",1.5).attr("stroke-dasharray","4,3");

  // rotated function
  const rotPts = d3.range(0, 2*Math.PI, 0.02).map(a => [xScale(a), yScale(Math.cos(4*(a - theta)))]);
  svg.append("path").attr("d", d3.line()(rotPts)).attr("fill","none").attr("stroke","firebrick").attr("stroke-width",2);

  // mark peaks
  for (let i = 0; i < 4; i++) {
    const peakAlpha = theta + i*Math.PI/2;
    if (peakAlpha <= 2*Math.PI) {
      svg.append("line").attr("x1",xScale(peakAlpha)).attr("y1",plotT).attr("x2",xScale(peakAlpha)).attr("y2",plotB)
        .attr("stroke","red").attr("stroke-width",0.8).attr("stroke-dasharray","3,3");
      svg.append("text").attr("x",xScale(peakAlpha)).attr("y",plotT-3).attr("text-anchor","middle").attr("font-size",10).attr("fill","red").text(`f${i}`);
    }
  }

  // axis labels
  svg.append("text").attr("x",plotR).attr("y",yScale(0)+15).attr("text-anchor","end").attr("font-size",11).attr("fill","#666").text("α");
  svg.append("text").attr("x",plotL-5).attr("y",plotT).attr("text-anchor","end").attr("font-size",11).attr("fill","#666").text("F(α)");
  // legend
  svg.append("line").attr("x1",plotL+10).attr("y1",plotB+15).attr("x2",plotL+30).attr("y2",plotB+15).attr("stroke","#bbb").attr("stroke-dasharray","4,3");
  svg.append("text").attr("x",plotL+34).attr("y",plotB+19).attr("font-size",10).attr("fill","#999").text("F̃ (ref)");
  svg.append("line").attr("x1",plotL+90).attr("y1",plotB+15).attr("x2",plotL+110).attr("y2",plotB+15).attr("stroke","firebrick").attr("stroke-width",2);
  svg.append("text").attr("x",plotL+114).attr("y",plotB+19).attr("font-size",10).attr("fill","firebrick").text("F (rotated)");

  // --- Panel 3: Coefficient space ---
  const p3x = 520, p3w = 220, p3cx = p3x+p3w/2, p3cy = H/2, p3r = 85;

  svg.append("text").attr("x",p3cx).attr("y",18).attr("text-anchor","middle").attr("font-size",13).attr("fill","#333").text("Coefficient space");

  // unit circle
  const circ3 = d3.range(0, 2*Math.PI+0.01, 0.02).map(a => [p3cx+p3r*Math.cos(a), p3cy-p3r*Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circ3)).attr("fill","none").attr("stroke","#ccc").attr("stroke-width",1.5);

  // axes
  svg.append("line").attr("x1",p3cx-p3r-15).attr("y1",p3cy).attr("x2",p3cx+p3r+15).attr("y2",p3cy).attr("stroke","#ddd");
  svg.append("line").attr("x1",p3cx).attr("y1",p3cy-p3r-15).attr("x2",p3cx).attr("y2",p3cy+p3r+15).attr("stroke","#ddd");

  // coefficient point
  const a0 = Math.cos(4*theta), a1 = Math.sin(4*theta);
  svg.append("line").attr("x1",p3cx).attr("y1",p3cy).attr("x2",p3cx+p3r*a0).attr("y2",p3cy-p3r*a1).attr("stroke","firebrick").attr("stroke-width",1.5).attr("stroke-dasharray","3,3");
  svg.append("circle").attr("cx",p3cx+p3r*a0).attr("cy",p3cy-p3r*a1).attr("r",6).attr("fill","firebrick");
  svg.append("text").attr("x",p3cx+p3r*a0+10).attr("y",p3cy-p3r*a1-8).attr("font-size",12).attr("fill","firebrick")
    .text(`a = (${a0.toFixed(2)}, ${a1.toFixed(2)})`);

  // axis labels
  svg.append("text").attr("x",p3cx+p3r+15).attr("y",p3cy-8).attr("font-size",11).attr("fill","#666").text("a₀");
  svg.append("text").attr("x",p3cx+6).attr("y",p3cy-p3r-10).attr("font-size",11).attr("fill","#666").text("a₁");

  return svg.node();
}

Parametric polar plot (cf. paper Fig. 3)

The paper visualizes frame functions as polar plots using the parametric curve \(x(\alpha) = (1 + F(\alpha))\cos(\alpha)\), \(y(\alpha) = (1 + F(\alpha))\sin(\alpha)\). The resulting “clover” shape has lobes pointing in the direction of each frame vector — the maxima of \(F(\alpha)\) become the tips of the lobes.

“The plot is made using \(x(\alpha) = (1 + F(\alpha))\cos(\alpha)\) and \(y(\alpha) = (1 + F(\alpha))\sin(\alpha)\) for \(\alpha \in [0, 2\pi[\). It is easy to see that corresponding frames are aligned with maxima of the representation functions.” (Fig. 3 caption)

viewof theta_polar_deg = Inputs.range([0, 90], {value: 0, step: 1, label: "θ (degrees)"})

Figure 3

{
  const theta = theta_polar_deg * Math.PI / 180;
  const W = 650, H = 320;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  const panels = [
    {cx: 160, cy: H/2, label: "Reference F̃(α) = cos(4α)", th: 0, color: "#888"},
    {cx: 490, cy: H/2, label: `Rotated F(α) = cos(4(α−${theta_polar_deg}°))`, th: theta, color: "firebrick"}
  ];

  const R = 120; // scale for the polar plot

  for (const p of panels) {
    svg.append("text").attr("x", p.cx).attr("y", 18).attr("text-anchor", "middle").attr("font-size", 12).attr("fill", "#333").text(p.label);

    // Light unit circle for reference
    const unitCirc = d3.range(0, 2*Math.PI+0.01, 0.02).map(a => [p.cx + R*0.5*Math.cos(a), p.cy - R*0.5*Math.sin(a)]);
    svg.append("path").attr("d", d3.line()(unitCirc)).attr("fill", "none").attr("stroke", "#eee").attr("stroke-width", 0.5);

    // Axes
    svg.append("line").attr("x1", p.cx - R - 15).attr("y1", p.cy).attr("x2", p.cx + R + 15).attr("y2", p.cy).attr("stroke", "#eee");
    svg.append("line").attr("x1", p.cx).attr("y1", p.cy - R - 15).attr("x2", p.cx).attr("y2", p.cy + R + 15).attr("stroke", "#eee");

    // Parametric polar plot: r(α) = 1 + F(α) = 1 + cos(4(α - θ))
    const pts = d3.range(0, 2*Math.PI + 0.01, 0.015).map(a => {
      const r = (1 + Math.cos(4*(a - p.th))) * R * 0.5;
      return [p.cx + r*Math.cos(a), p.cy - r*Math.sin(a)];
    });
    svg.append("path").attr("d", d3.line()(pts)).attr("fill", p.color).attr("fill-opacity", 0.1)
      .attr("stroke", p.color).attr("stroke-width", 2);

    // Frame vectors at the tips
    const frameColors = ["steelblue", "#d35", "steelblue", "#d35"];
    for (let i = 0; i < 4; i++) {
      const a = p.th + i * Math.PI/2;
      const tipR = R; // r = (1+1)*R*0.5 = R at maxima
      const dx = tipR * Math.cos(a), dy = -tipR * Math.sin(a);
      svg.append("line").attr("x1", p.cx).attr("y1", p.cy).attr("x2", p.cx + dx).attr("y2", p.cy + dy)
        .attr("stroke", frameColors[i]).attr("stroke-width", 1.5).attr("stroke-dasharray", "4,3");
      svg.append("circle").attr("cx", p.cx + dx).attr("cy", p.cy + dy).attr("r", 3).attr("fill", frameColors[i]);
    }

    // Center dot
    svg.append("circle").attr("cx", p.cx).attr("cy", p.cy).attr("r", 2.5).attr("fill", "#333");
  }

  // Arrow between panels
  svg.append("line").attr("x1", 280).attr("y1", H/2).attr("x2", 360).attr("y2", H/2).attr("stroke", "#aaa").attr("stroke-width", 1.5);
  svg.append("polygon").attr("points", "360," + H/2 + " 352," + (H/2-5) + " 352," + (H/2+5)).attr("fill", "#aaa");
  svg.append("text").attr("x", 320).attr("y", H/2 - 10).attr("text-anchor", "middle").attr("font-size", 12).attr("fill", "#888").text("rotate by θ");

  return svg.node();
}

The four-lobed “clover” shape is the signature of a frame: the lobes point exactly where the frame vectors point. Rotating the frame by θ rotates the entire clover. At θ = 0° (reference frame), the lobes align with the coordinate axes. Drag the slider and watch the lobes follow the frame vectors.

---

2.3 Functional Approach: Objective Function

“We have defined the field curvature as the sum over each edge of the squared difference between frames at the edges extremities. In our formalization, the difference between two frames (at vertices \(i\) and \(j\)) is no longer the rotation angle, but the \(L^2\) norm of the difference between the corresponding functions: \(\int_0^{2\pi} (F^j(\alpha) - F^i(\alpha))^2 d\alpha\).”

The paper derives the energy:

\[E = \sum_{ij} \int_0^{2\pi} (F^j(\alpha) - F^i(\alpha))^2 \, d\alpha\] \[= \sum_{ij} \int_0^{2\pi} (Ba^j - Ba^i)^2 \, d\alpha\] \[= \sum_{ij} (a^j - a^i)^\top \left(\int_0^{2\pi} B^\top B \, d\alpha \right) (a^j - a^i)\]

“As the function basis \(B\) is orthogonal, and all functions are of norm \(\sqrt{\pi}\), the expression simplifies to:”

\[E = \pi \sum_{ij} \|a^j - a^i\|^2 \quad \text{(Eq. 1)}\]

This is remarkable: the smoothness energy in function space is simply the sum of squared Euclidean distances between coefficient vectors!

Drag the two frame angles to compare arc length (field curvature) vs chord length (our energy \(E\)):

viewof thetaI_deg = Inputs.range([0, 90], {value: 0, step: 1, label: "θⁱ (degrees)"})
viewof thetaJ_deg = Inputs.range([0, 90], {value: 30, step: 1, label: "θʲ (degrees)"})

(a)

(b)

Figure 4

{
  const tI = thetaI_deg * Math.PI / 180, tJ = thetaJ_deg * Math.PI / 180;
  const W = 500, H = 500, cx = W/2, cy = H/2, r = 180;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  // Unit circle
  const circ = d3.range(0, 2*Math.PI+0.01, 0.02).map(a => [cx+r*Math.cos(a), cy-r*Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circ)).attr("fill","none").attr("stroke","#ddd").attr("stroke-width",1.5);

  // Coefficient vectors (4θ maps to coefficient space)
  const ai = [Math.cos(4*tI), Math.sin(4*tI)];
  const aj = [Math.cos(4*tJ), Math.sin(4*tJ)];

  // Arc between them on unit circle
  let arcStart = 4*tI, arcEnd = 4*tJ;
  if (arcEnd < arcStart) [arcStart, arcEnd] = [arcEnd, arcStart];
  const arcPts = d3.range(arcStart, arcEnd+0.01, 0.02).map(a => [cx+r*Math.cos(a), cy-r*Math.sin(a)]);
  if (arcPts.length > 1) {
    svg.append("path").attr("d", d3.line()(arcPts)).attr("fill","none").attr("stroke","orange").attr("stroke-width",4).attr("opacity",0.8);
  }

  // Chord
  svg.append("line").attr("x1",cx+r*ai[0]).attr("y1",cy-r*ai[1]).attr("x2",cx+r*aj[0]).attr("y2",cy-r*aj[1])
    .attr("stroke","steelblue").attr("stroke-width",3).attr("stroke-dasharray","6,4");

  // Points
  svg.append("circle").attr("cx",cx+r*ai[0]).attr("cy",cy-r*ai[1]).attr("r",7).attr("fill","firebrick");
  svg.append("text").attr("x",cx+r*ai[0]-15).attr("y",cy-r*ai[1]-12).attr("font-size",15).attr("fill","firebrick").attr("font-weight","bold").text("aⁱ");

  svg.append("circle").attr("cx",cx+r*aj[0]).attr("cy",cy-r*aj[1]).attr("r",7).attr("fill","forestgreen");
  svg.append("text").attr("x",cx+r*aj[0]+8).attr("y",cy-r*aj[1]-12).attr("font-size",15).attr("fill","forestgreen").attr("font-weight","bold").text("aʲ");

  // Compute values
  const arcLen = Math.abs(arcEnd - arcStart);
  const chordLen = Math.sqrt((aj[0]-ai[0])**2 + (aj[1]-ai[1])**2);

  // Info box
  svg.append("text").attr("x",20).attr("y",H-50).attr("font-size",14).attr("fill","orange").text(`Arc length (4Δθ): ${arcLen.toFixed(3)}`);
  svg.append("text").attr("x",20).attr("y",H-30).attr("font-size",14).attr("fill","steelblue").text(`Chord length ‖aʲ−aⁱ‖: ${chordLen.toFixed(3)}`);
  svg.append("text").attr("x",20).attr("y",H-10).attr("font-size",13).attr("fill","#666").text(`Ratio (chord/arc): ${arcLen > 0.001 ? (chordLen/arcLen).toFixed(3) : "—"}`);

  // Legend
  svg.append("line").attr("x1",20).attr("y1",25).attr("x2",50).attr("y2",25).attr("stroke","orange").attr("stroke-width",3);
  svg.append("text").attr("x",55).attr("y",29).attr("font-size",12).attr("fill","#666").text("Arc (rotation angle)");
  svg.append("line").attr("x1",220).attr("y1",25).attr("x2",250).attr("y2",25).attr("stroke","steelblue").attr("stroke-width",3).attr("stroke-dasharray","6,4");
  svg.append("text").attr("x",255).attr("y",29).attr("font-size",12).attr("fill","#666").text("Chord (‖aʲ−aⁱ‖)");

  return svg.node();
}

Notice that for small rotations, arc \(\approx\) chord. For large rotations, the chord is shorter — the energy \(E\) is smoother than the field curvature near singularities.

---

2.4 Functional Approach: Constraints

There are two constraints:

Feasibility: Each coefficient vector must represent a valid frame.

“As discussed in §2.2 the first constraint is clearly that the variables \(a^i\) must be feasible (i.e. there exists a frame represented by \(a^i\)).”

This means \(a^\top a = 1\) — the coefficient vector must lie on the unit circle.

Boundary alignment: Frames at boundary vertices must align with the boundary normal.

“The second constraint is that frames of boundary vertices \(i\) must have one member vector equals to the normal direction. If \(\theta^i\) denotes the normal direction, the frame can be directly fixed by satisfying two equations:”

\[a_0^i = \cos(4\theta^i) \quad \quad a_1^i = \sin(4\theta^i) \quad \text{(Eq. 2)}\]

“However, as we already have the feasibility constraint \(a^{i\top} a^i = 1\), enforcing only one equation has the same effect:”

\[a_0^i \cos(4\theta^i) + a_1^i \sin(4\theta^i) = 1 \quad \text{(Eq. 3)}\]

Drag to change the boundary normal direction and see how it pins the coefficient vector:

viewof normal_deg = Inputs.range([0, 360], {value: 90, step: 1, label: "Normal angle θⁿ (degrees)"})

Figure 5

{
  const theta_n = normal_deg * Math.PI / 180;
  const W = 700, H = 360;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  // --- Left panel: physical space ---
  const lcx = 150, lcy = H/2 + 20, lr = 100;

  svg.append("text").attr("x",lcx).attr("y",20).attr("text-anchor","middle").attr("font-size",14).attr("fill","#333").text("Physical space: boundary frame");

  // Boundary line (perpendicular to normal)
  const bx = Math.cos(theta_n + Math.PI/2), by = Math.sin(theta_n + Math.PI/2);
  svg.append("line").attr("x1",lcx-150*bx).attr("y1",lcy+150*by).attr("x2",lcx+150*bx).attr("y2",lcy-150*by)
    .attr("stroke","#888").attr("stroke-width",3);

  // Normal arrow
  const nx = 80*Math.cos(theta_n), ny = -80*Math.sin(theta_n);
  svg.append("line").attr("x1",lcx).attr("y1",lcy).attr("x2",lcx+nx).attr("y2",lcy+ny).attr("stroke","orange").attr("stroke-width",3);
  svg.append("text").attr("x",lcx+nx*1.2).attr("y",lcy+ny*1.2).attr("text-anchor","middle").attr("font-size",16).attr("fill","orange").attr("font-weight","bold").text("n̂");

  // Frame aligned with normal
  const fColors = ["steelblue","#4a9","steelblue","#4a9"];
  for (let i = 0; i < 4; i++) {
    const a = theta_n + i*Math.PI/2;
    const fdx = 60*Math.cos(a), fdy = -60*Math.sin(a);
    svg.append("line").attr("x1",lcx).attr("y1",lcy).attr("x2",lcx+fdx).attr("y2",lcy+fdy)
      .attr("stroke",fColors[i]).attr("stroke-width",2);
  }

  // --- Right panel: coefficient space ---
  const rcx = 500, rcy = H/2 + 20, rr = 130;

  svg.append("text").attr("x",rcx).attr("y",20).attr("text-anchor","middle").attr("font-size",14).attr("fill","#333").text("Coefficient space: boundary constraint");

  // Unit circle
  const circ = d3.range(0, 2*Math.PI+0.01, 0.02).map(a => [rcx+rr*Math.cos(a), rcy-rr*Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circ)).attr("fill","none").attr("stroke","#ccc").attr("stroke-width",2);

  // Axes
  svg.append("line").attr("x1",rcx-rr-20).attr("y1",rcy).attr("x2",rcx+rr+20).attr("y2",rcy).attr("stroke","#eee");
  svg.append("line").attr("x1",rcx).attr("y1",rcy-rr-20).attr("x2",rcx).attr("y2",rcy+rr+20).attr("stroke","#eee");
  svg.append("text").attr("x",rcx+rr+22).attr("y",rcy-5).attr("font-size",11).attr("fill","#999").text("a₀");
  svg.append("text").attr("x",rcx+5).attr("y",rcy-rr-12).attr("font-size",11).attr("fill","#999").text("a₁");

  // Constraint line: a₀ cos(4θⁿ) + a₁ sin(4θⁿ) = 1
  const c0 = Math.cos(4*theta_n), c1 = Math.sin(4*theta_n);
  // Draw constraint line across the panel
  const linePts = [];
  for (let t = -2; t <= 2; t += 0.01) {
    // parametric: point on line = c + t * perp
    const px = c0 + t*(-c1), py = c1 + t*c0;
    const sx = rcx + rr*px, sy = rcy - rr*py;
    linePts.push([sx, sy]);
  }
  // filter to viewport
  const filtered = linePts.filter(([x,y]) => x > rcx-rr-20 && x < rcx+rr+20 && y > rcy-rr-20 && y < rcy+rr+20);
  if (filtered.length > 1)
    svg.append("path").attr("d", d3.line()(filtered)).attr("fill","none").attr("stroke","red").attr("stroke-width",1.5).attr("stroke-dasharray","5,4");

  // The constrained point
  const a0 = Math.cos(4*theta_n), a1 = Math.sin(4*theta_n);
  svg.append("circle").attr("cx",rcx+rr*a0).attr("cy",rcy-rr*a1).attr("r",7).attr("fill","orange");
  svg.append("text").attr("x",rcx+rr*a0+10).attr("y",rcy-rr*a1-10).attr("font-size",12).attr("fill","orange")
    .text(`a = (${a0.toFixed(2)}, ${a1.toFixed(2)})`);

  // Legend
  svg.append("line").attr("x1",rcx-rr).attr("y1",rcy+rr+20).attr("x2",rcx-rr+25).attr("y2",rcy+rr+20).attr("stroke","red").attr("stroke-dasharray","5,4");
  svg.append("text").attr("x",rcx-rr+30).attr("y",rcy+rr+24).attr("font-size",11).attr("fill","#666").text("a₀cos4θⁿ + a₁sin4θⁿ = 1 (Eq. 3)");

  return svg.node();
}

---

2.5 Toy Example

“It is natural to ask the question: ‘Does minimizing our energy minimizes the field curvature as well?’”

“Two frames \(f^i\) and \(f^j\) are represented by \(a^i\) and \(a^j\), both located on the unit circle. The field curvature measures the circle arc length between them, whereas our \(L^2\) norm is the chord length between \(a^i\) and \(a^j\).”

“From a practical point of view, we want to produce smooth fields, so most edges will have low curvature. In this case the objective function \(E\) is almost proportional to the field curvature. If, however, two adjacent frames are not similar (e.g. they are close to singularities), then the function \(E\) is smoother than the field curvature, making it easier for the optimization algorithm to move singularities.”

The paper introduces a toy problem: a chain of 4 vertices with the two endpoints locked.

“Let us illustrate our intuition on a very simple interpolation example: a chain of four vertices having its extremities locked. The toy problem is therefore to find two frames interpolating the extremity frames.”

Drag θ₁ and θ₂ to manually interpolate between the locked endpoints. Watch the energy values update:

viewof theta0_lock = Inputs.range([0, 90], {value: 0, step: 1, label: "θ₀ locked (degrees)"})
viewof theta3_lock = Inputs.range([0, 90], {value: 60, step: 1, label: "θ₃ locked (degrees)"})
viewof theta1_free = Inputs.range([-45, 90], {value: 20, step: 1, label: "θ₁ free (degrees)"})
viewof theta2_free = Inputs.range([-45, 90], {value: 40, step: 1, label: "θ₂ free (degrees)"})

(a)

(b)

(c)

(d)

Figure 6

{
  const t0 = theta0_lock*Math.PI/180, t1 = theta1_free*Math.PI/180;
  const t2 = theta2_free*Math.PI/180, t3 = theta3_lock*Math.PI/180;
  const thetas = [t0, t1, t2, t3];

  const W = 700, H = 400;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  // --- Left: frame chain ---
  const lw = 350;
  svg.append("text").attr("x", lw/2).attr("y", 20).attr("text-anchor","middle").attr("font-size",14).attr("fill","#333").text("Frame chain");

  const chainY = 100;
  const xs = [40, 130, 220, 310];
  const labels = ["P₀","P₁","P₂","P₃"];
  const isLocked = [true, false, false, true];
  const colrs = ["orange","steelblue","steelblue","orange"];

  // edges
  svg.append("line").attr("x1",xs[0]).attr("y1",chainY).attr("x2",xs[3]).attr("y2",chainY).attr("stroke","#bbb").attr("stroke-width",2);

  for (let i = 0; i < 4; i++) {
    // frame vectors
    for (let k = 0; k < 4; k++) {
      const a = thetas[i] + k*Math.PI/2;
      svg.append("line").attr("x1",xs[i]).attr("y1",chainY).attr("x2",xs[i]+28*Math.cos(a)).attr("y2",chainY-28*Math.sin(a))
        .attr("stroke",colrs[i]).attr("stroke-width",2);
    }
    svg.append("circle").attr("cx",xs[i]).attr("cy",chainY).attr("r",4).attr("fill",colrs[i]);
    svg.append("text").attr("x",xs[i]).attr("y",chainY+45).attr("text-anchor","middle").attr("font-size",12).attr("fill",colrs[i])
      .text(labels[i] + (isLocked[i] ? " 🔒" : ""));
  }

  // --- Energy computations ---
  function rotAngle(a, b) {
    let d = b - a;
    d = ((d + Math.PI/4) % (Math.PI/2) + Math.PI/2) % (Math.PI/2) - Math.PI/4;
    return d;
  }

  const fieldCurv = rotAngle(t0,t1)**2 + rotAngle(t1,t2)**2 + rotAngle(t2,t3)**2;

  const coeffs = thetas.map(t => [Math.cos(4*t), Math.sin(4*t)]);
  let objE = 0;
  for (let e = 0; e < 3; e++) {
    objE += (coeffs[e+1][0]-coeffs[e][0])**2 + (coeffs[e+1][1]-coeffs[e][1])**2;
  }
  objE *= Math.PI;

  // Energy display
  const ey = 180;
  svg.append("text").attr("x",20).attr("y",ey).attr("font-size",13).attr("fill","#333").attr("font-weight","bold").text("Energy comparison:");
  svg.append("text").attr("x",20).attr("y",ey+22).attr("font-size",13).attr("fill","purple").text(`Field curvature Σ(Δθ)² = ${fieldCurv.toFixed(4)}`);
  svg.append("text").attr("x",20).attr("y",ey+44).attr("font-size",13).attr("fill","teal").text(`Objective E = π Σ‖aʲ−aⁱ‖² = ${objE.toFixed(4)}`);

  // --- Right: coefficient space ---
  const rcx = 530, rcy = H/2, rr = 140;
  svg.append("text").attr("x",rcx).attr("y",20).attr("text-anchor","middle").attr("font-size",14).attr("fill","#333").text("Coefficient space");

  const circ = d3.range(0, 2*Math.PI+0.01, 0.02).map(a => [rcx+rr*Math.cos(a), rcy-rr*Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circ)).attr("fill","none").attr("stroke","#ddd").attr("stroke-width",1.5);

  // axes
  svg.append("line").attr("x1",rcx-rr-10).attr("y1",rcy).attr("x2",rcx+rr+10).attr("y2",rcy).attr("stroke","#eee");
  svg.append("line").attr("x1",rcx).attr("y1",rcy-rr-10).attr("x2",rcx).attr("y2",rcy+rr+10).attr("stroke","#eee");

  // edges in coefficient space (chords)
  for (let e = 0; e < 3; e++) {
    svg.append("line")
      .attr("x1",rcx+rr*coeffs[e][0]).attr("y1",rcy-rr*coeffs[e][1])
      .attr("x2",rcx+rr*coeffs[e+1][0]).attr("y2",rcy-rr*coeffs[e+1][1])
      .attr("stroke","#aaa").attr("stroke-width",1.5);
  }

  // points
  for (let i = 0; i < 4; i++) {
    svg.append("circle").attr("cx",rcx+rr*coeffs[i][0]).attr("cy",rcy-rr*coeffs[i][1]).attr("r",6).attr("fill",colrs[i]);
    svg.append("text").attr("x",rcx+rr*coeffs[i][0]+10).attr("y",rcy-rr*coeffs[i][1]-8).attr("font-size",12).attr("fill",colrs[i]).text(labels[i]);
  }

  return svg.node();
}

Energy landscape and local minima (cf. paper Fig. 4 & 5)

“The field curvature is not smooth (it is piecewise quadratic) and we can observe that there are two local minima. Our objective function is smooth, and has exactly the same minima on this example.”

To see this, we plot both energies as functions of θ₁ and θ₂ (the free angles). The locked endpoints define two natural interpolation paths — counterclockwise (short way) and clockwise (long way around the π/2-periodic circle). Each path yields a local minimum.

viewof land_t0 = Inputs.range([0, 90], {value: 0, step: 1, label: "θ₀ locked (degrees)"})
viewof land_t3 = Inputs.range([0, 90], {value: 60, step: 1, label: "θ₃ locked (degrees)"})

(a)

(b)

Figure 7

{
  const t0 = land_t0 * Math.PI / 180, t3 = land_t3 * Math.PI / 180;
  const N = 80;
  const lo = -Math.PI/4, hi = Math.PI/2;
  const step = (hi - lo) / (N - 1);

  // Rotation angle mod π/2, in [−π/4, π/4]
  function rotAngle(a, b) {
    let d = b - a;
    d = ((d + Math.PI/4) % (Math.PI/2) + Math.PI/2) % (Math.PI/2) - Math.PI/4;
    return d;
  }

  // Compute energy grids
  const curvGrid = [], objGrid = [];
  let curvMin = Infinity, objMin = Infinity, curvMax = 0, objMax = 0;
  for (let j = 0; j < N; j++) {
    curvGrid[j] = []; objGrid[j] = [];
    const t2 = lo + j * step;
    for (let i = 0; i < N; i++) {
      const t1 = lo + i * step;
      const curv = rotAngle(t0,t1)**2 + rotAngle(t1,t2)**2 + rotAngle(t2,t3)**2;
      const c0 = [Math.cos(4*t0), Math.sin(4*t0)];
      const c1 = [Math.cos(4*t1), Math.sin(4*t1)];
      const c2 = [Math.cos(4*t2), Math.sin(4*t2)];
      const c3 = [Math.cos(4*t3), Math.sin(4*t3)];
      const obj = ((c1[0]-c0[0])**2+(c1[1]-c0[1])**2 + (c2[0]-c1[0])**2+(c2[1]-c1[1])**2 + (c3[0]-c2[0])**2+(c3[1]-c2[1])**2) * Math.PI;
      curvGrid[j][i] = curv;
      objGrid[j][i] = obj;
      curvMin = Math.min(curvMin, curv); curvMax = Math.max(curvMax, curv);
      objMin = Math.min(objMin, obj); objMax = Math.max(objMax, obj);
    }
  }

  // Find minima by brute force
  function findMinima(grid) {
    let mins = [];
    for (let j = 1; j < N-1; j++) {
      for (let i = 1; i < N-1; i++) {
        const v = grid[j][i];
        if (v < grid[j-1][i] && v < grid[j+1][i] && v < grid[j][i-1] && v < grid[j][i+1]
            && v < grid[j-1][i-1] && v < grid[j-1][i+1] && v < grid[j+1][i-1] && v < grid[j+1][i+1]) {
          mins.push({i, j, v, t1: lo + i*step, t2: lo + j*step});
        }
      }
    }
    mins.sort((a,b) => a.v - b.v);
    return mins;
  }

  // Helper to draw a frame chain
  function drawFrameChain(svg, thetas, x0, y0, spacing, frameR, label, labelColor) {
    const chainColors = ["orange","steelblue","steelblue","orange"];
    svg.append("line").attr("x1", x0).attr("y1", y0).attr("x2", x0 + 3*spacing).attr("y2", y0)
      .attr("stroke","#ccc").attr("stroke-width",1);
    for (let i = 0; i < 4; i++) {
      const cx = x0 + i*spacing;
      for (let k = 0; k < 4; k++) {
        const ang = thetas[i] + k*Math.PI/2;
        svg.append("line").attr("x1", cx).attr("y1", y0)
          .attr("x2", cx + frameR*Math.cos(ang)).attr("y2", y0 - frameR*Math.sin(ang))
          .attr("stroke", chainColors[i]).attr("stroke-width", 1.5);
      }
      svg.append("circle").attr("cx", cx).attr("cy", y0).attr("r", 2).attr("fill", chainColors[i]);
    }
    if (label) {
      svg.append("text").attr("x", x0 - 5).attr("y", y0 + 5).attr("text-anchor","end")
        .attr("font-size", 11).attr("fill", labelColor || "red").attr("font-weight","bold").text(label);
    }
  }

  const W = 700, H = 580;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  const panels = [
    {x: 10, w: 310, grid: curvGrid, min: curvMin, max: curvMax, title: "Field curvature Σ(Δθ)²"},
    {x: 370, w: 310, grid: objGrid, min: objMin, max: objMax, title: "Objective E = π Σ‖aʲ−aⁱ‖²"}
  ];

  const cSize = 280, cTop = 40;

  // Viridis-ish color scale
  function viridis(t) {
    const r = Math.max(0, Math.min(255, Math.round(68 + t * (253-68))));
    const g = Math.max(0, Math.min(255, Math.round(1 + t * (231-1))));
    const b = Math.max(0, Math.min(255, Math.round(84 + (t < 0.5 ? t*2*(168-84) : (168 + (t-0.5)*2*(37-168))))));
    return `rgb(${r},${g},${b})`;
  }

  // Collect all minima across both panels (use objGrid for the canonical minima)
  const objMins = findMinima(objGrid);

  for (const p of panels) {
    svg.append("text").attr("x", p.x + cSize/2 + 5).attr("y", 16).attr("text-anchor","middle").attr("font-size",12).attr("fill","#333").text(p.title);

    // Draw heatmap
    const cellW = cSize / N, cellH = cSize / N;
    const range = p.max - p.min || 1;
    for (let j = 0; j < N; j++) {
      for (let i = 0; i < N; i++) {
        const t = (p.grid[j][i] - p.min) / range;
        svg.append("rect")
          .attr("x", p.x + i * cellW)
          .attr("y", cTop + (N - 1 - j) * cellH)
          .attr("width", cellW + 0.5)
          .attr("height", cellH + 0.5)
          .attr("fill", viridis(t));
      }
    }

    // Axis labels
    svg.append("text").attr("x", p.x + cSize/2).attr("y", cTop + cSize + 18).attr("text-anchor","middle").attr("font-size",11).attr("fill","#666").text("θ₁");
    svg.append("text").attr("x", p.x - 6).attr("y", cTop + cSize/2).attr("text-anchor","middle").attr("font-size",11).attr("fill","#666")
      .attr("transform", `rotate(-90, ${p.x - 6}, ${cTop + cSize/2})`).text("θ₂");

    // Tick labels
    const ticks = [-45, 0, 45, 90];
    for (const deg of ticks) {
      const rad = deg * Math.PI / 180;
      const frac = (rad - lo) / (hi - lo);
      svg.append("text").attr("x", p.x + frac * cSize).attr("y", cTop + cSize + 30).attr("text-anchor","middle").attr("font-size",9).attr("fill","#999").text(deg + "°");
    }

    // Mark minima
    const mins = findMinima(p.grid);
    const minLabels = ["P₀", "P₁"];
    const minColors = ["teal", "purple"];
    for (let k = 0; k < Math.min(mins.length, 2); k++) {
      const m = mins[k];
      const mx = p.x + (m.i / (N-1)) * cSize;
      const my = cTop + (1 - m.j / (N-1)) * cSize;
      svg.append("circle").attr("cx", mx).attr("cy", my).attr("r", 6).attr("fill","none").attr("stroke",minColors[k]).attr("stroke-width",2.5);
      svg.append("text").attr("x", mx + 8).attr("y", my - 6).attr("font-size", 13).attr("fill",minColors[k]).attr("font-weight","bold").text(minLabels[k]);
    }
  }

  // --- Frame chain insets below the heatmaps ---
  const insetTop = cTop + cSize + 50;
  const minLabels2 = ["P₀", "P₁"];
  const minColors2 = ["teal", "purple"];

  svg.append("text").attr("x", W/2).attr("y", insetTop - 8).attr("text-anchor","middle").attr("font-size", 13).attr("fill","#333").attr("font-weight","bold")
    .text("Frames at each minimum");

  // Compute the two minima frame chains
  // P₀: CCW interpolation
  const ccw_thetas = [t0, t0 + (t3-t0)/3, t0 + 2*(t3-t0)/3, t3];
  // P₁: CW interpolation (wrap through π/2)
  const gap_cw = (Math.PI/2 - (t3 - t0));
  const cw_thetas = [t0, t0 - gap_cw/3, t0 - 2*gap_cw/3, t3];

  const minThetas = [ccw_thetas, cw_thetas];

  for (let k = 0; k < 2; k++) {
    const rowY = insetTop + 20 + k * 80;
    const chainX = 100;
    const spacing = 110;

    drawFrameChain(svg, minThetas[k], chainX, rowY, spacing, 22, minLabels2[k], minColors2[k]);

    // Angle labels under each vertex
    const vertLabels = ["θ₀","θ₁","θ₂","θ₃"];
    for (let i = 0; i < 4; i++) {
      const cx = chainX + i * spacing;
      const deg = (minThetas[k][i] * 180 / Math.PI).toFixed(0);
      svg.append("text").attr("x", cx).attr("y", rowY + 32).attr("text-anchor","middle").attr("font-size",9).attr("fill","#888").text(`${deg}°`);
    }

    // Δθ labels between vertices
    for (let i = 0; i < 3; i++) {
      const x1 = chainX + i*spacing + 25, x2 = chainX + (i+1)*spacing - 25;
      const diff = minThetas[k][i+1] - minThetas[k][i];
      const deg = (diff * 180 / Math.PI).toFixed(1);
      svg.append("text").attr("x", (x1+x2)/2).attr("y", rowY - 26).attr("text-anchor","middle").attr("font-size",9).attr("fill",minColors2[k]).text(`Δθ=${deg}°`);
    }

    // Energy value
    const coeffs = minThetas[k].map(t => [Math.cos(4*t), Math.sin(4*t)]);
    let E = 0;
    for (let e = 0; e < 3; e++) E += (coeffs[e+1][0]-coeffs[e][0])**2 + (coeffs[e+1][1]-coeffs[e][1])**2;
    E *= Math.PI;
    svg.append("text").attr("x", chainX + 3*spacing + 30).attr("y", rowY + 5).attr("font-size", 10).attr("fill", minColors2[k])
      .text(`E = ${E.toFixed(2)}`);
  }

  return svg.node();
}

Key observations:

• The field curvature (left) has sharp creases — it is piecewise quadratic because the rotation angle wraps modulo π/2. A gradient descent optimizer can get stuck on the creases.
• The objective function E (right) is smooth (globally quadratic in the coefficient vectors), making it much easier to optimize.
• Both share the same local minima P₀ and P₁, so minimizing E finds the same frames.
• P₀ (teal) interpolates counterclockwise with small, uniform rotation steps — this is the global minimum.
• P₁ (purple) takes the long way around clockwise with larger jumps — only a local minimum.

The smoothness of E makes it easier for gradient-based methods to find P₀ and escape from P₁.

---

2.6 Implementation

“We have to minimize our objective function \(E\) (eq. 1) with linear equality constraints on boundary vertices (eq. 2 or eq. 3) and quadratic equality constraints \(a^{i\top} a^i = 1\) for each vertex \(i\).”

“We use [Kowalski et al. 2012]’s algorithm to solve this problem. It finds an initial solution by relaxing the quadratic feasibility constraints on \(a^i\) and finding the nearest feasible solution. Then it performs a number of smoothing iterations to ameliorate the solution.”

Initialization

“Here, we relax the feasibility constraints, so we need to minimize the quadratic function \(E\) subject to linear boundary constraints. To do it, we simply replace the linear constraints by a strong penalty term in the objective function, leading to a new quadratic function to optimize.”

“More precisely, the new quadratic function is expressed in the form \(\|AX - b\|^2\) where \(A\) is a matrix, \(X\) our variable vector (\(X_{2i} = a_0^i\) and \(X_{2i+1} = a_1^i\)) and \(b\) is a vector.”

The system \(AX - b = 0\) is built line-by-line:

• “initial objective function \(E\): for each edge \(ij\), we add two equations: \(\sqrt{\pi}(X_{2i} - X_{2j}) = 0\) and \(\sqrt{\pi}(X_{2i+1} - X_{2j+1}) = 0\)“*

• “boundary constraints: for each boundary vertex \(i\), we add two equations: \(\lambda X_{2i} = \lambda \cos 4\theta^i\) and \(\lambda X_{2i+1} = \lambda \sin 4\theta^i\), where we set \(\lambda = 100\) in our experiments.”*

“From \(A\) and \(b\), we just need to solve the linear system \(A^\top A X = A^\top b\) to obtain a minimizer of \(\|AX - b\|^2\). Then from \(X\) we can obtain an initialization of \(a^i\) by projecting corresponding vectors on the set of feasible coefficients:

\(a^i \leftarrow (X_{2i}, X_{2i+1})^\top / \|(X_{2i}, X_{2i+1})\|\).”*

Here is the actual \(A\) matrix and \(b\) vector for our toy problem (4 vertices, 3 edges, vertices 0 and 3 locked). The numerical values update with the locked angles:

viewof mat_t0_deg = Inputs.range([0, 90], {value: 0, step: 1, label: "θ₀ locked (degrees)"})
viewof mat_t3_deg = Inputs.range([0, 90], {value: 60, step: 1, label: "θ₃ locked (degrees)"})
viewof mat_lambda = Inputs.range([1, 200], {value: 100, step: 1, label: "λ (boundary penalty)"})

(a)

(b)

(c)

Figure 8

{
  const t0 = mat_t0_deg * Math.PI / 180;
  const t3 = mat_t3_deg * Math.PI / 180;
  const sp = Math.sqrt(Math.PI);
  const lam = mat_lambda;

  // Computed boundary values
  const c4t0 = Math.cos(4*t0), s4t0 = Math.sin(4*t0);
  const c4t3 = Math.cos(4*t3), s4t3 = Math.sin(4*t3);

  const W = 700, H = 520;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  const smoothColor = "steelblue";
  const bdryColor = "#d97706";
  const feasColor = "#8b5cf6";

  // Layout
  const matL = 95, matT = 70;
  const cellW = 55, cellH = 26;
  const nCols = 8; // X = [a₀⁰, a₁⁰, a₀¹, a₁¹, a₀², a₁², a₀³, a₁³]
  const colLabels = ["a₀⁰","a₁⁰","a₀¹","a₁¹","a₀²","a₁²","a₀³","a₁³"];

  const sp_s = sp.toFixed(2);
  const msp_s = (-sp).toFixed(2);

  // Format number concisely
  function fmt(v) {
    if (Math.abs(v) < 0.005) return "0";
    if (Math.abs(v - Math.round(v)) < 0.005) return Math.round(v).toString();
    return v.toFixed(2);
  }

  // Rows: A values (as numbers) and b values
  const rows = [
    // --- Smoothness: edge 0-1 ---
    {An: [sp,0,-sp,0,0,0,0,0], bn: 0, color: smoothColor, src: "E", note: "edge 0–1, a₀"},
    {An: [0,sp,0,-sp,0,0,0,0], bn: 0, color: smoothColor, src: "E", note: "edge 0–1, a₁"},
    // --- Smoothness: edge 1-2 ---
    {An: [0,0,sp,0,-sp,0,0,0], bn: 0, color: smoothColor, src: "E", note: "edge 1–2, a₀"},
    {An: [0,0,0,sp,0,-sp,0,0], bn: 0, color: smoothColor, src: "E", note: "edge 1–2, a₁"},
    // --- Smoothness: edge 2-3 ---
    {An: [0,0,0,0,sp,0,-sp,0], bn: 0, color: smoothColor, src: "E", note: "edge 2–3, a₀"},
    {An: [0,0,0,0,0,sp,0,-sp], bn: 0, color: smoothColor, src: "E", note: "edge 2–3, a₁"},
    // --- Boundary: vertex 0 ---
    {An: [lam,0,0,0,0,0,0,0], bn: lam*c4t0, color: bdryColor, src: "BC", note: "v₀, a₀"},
    {An: [0,lam,0,0,0,0,0,0], bn: lam*s4t0, color: bdryColor, src: "BC", note: "v₀, a₁"},
    // --- Boundary: vertex 3 ---
    {An: [0,0,0,0,0,0,lam,0], bn: lam*c4t3, color: bdryColor, src: "BC", note: "v₃, a₀"},
    {An: [0,0,0,0,0,0,0,lam], bn: lam*s4t3, color: bdryColor, src: "BC", note: "v₃, a₁"},
  ];
  const nRows = rows.length;

  // Title
  svg.append("text").attr("x", W/2).attr("y", 20).attr("text-anchor","middle").attr("font-size",14).attr("fill","#333")
    .text("Initialization system ‖AX − b‖² for the 4-vertex toy problem");
  svg.append("text").attr("x", W/2).attr("y", 38).attr("text-anchor","middle").attr("font-size",12).attr("fill","#888")
    .text(`θ₀ = ${mat_t0_deg}°, θ₃ = ${mat_t3_deg}°, λ = ${mat_lambda}`);

  // X vector label
  svg.append("text").attr("x", matL + nCols*cellW/2).attr("y", matT - 22)
    .attr("text-anchor","middle").attr("font-size",12).attr("fill","#333").attr("font-weight","bold")
    .text("A");

  // Column headers
  for (let c = 0; c < nCols; c++) {
    const isLocked = c < 2 || c >= 6;
    svg.append("text")
      .attr("x", matL + c*cellW + cellW/2).attr("y", matT - 8)
      .attr("text-anchor","middle").attr("font-size", 10).attr("fill", isLocked ? bdryColor : "#555")
      .attr("font-weight", isLocked ? "bold" : "normal")
      .text(colLabels[c]);
  }

  // b label
  const bCol = matL + nCols*cellW + 25;
  svg.append("text").attr("x", bCol + 28).attr("y", matT - 22)
    .attr("text-anchor","middle").attr("font-size",12).attr("fill","#333").attr("font-weight","bold")
    .text("b");

  // Note column
  const noteCol = bCol + 65;

  // Brackets for A
  const bkL = matL - 6, bkR = matL + nCols*cellW + 4;
  const bkT = matT - 2, bkB = matT + nRows*cellH + 2;
  svg.append("path").attr("d", `M${bkL+4},${bkT} L${bkL},${bkT} L${bkL},${bkB} L${bkL+4},${bkB}`)
    .attr("fill","none").attr("stroke","#999").attr("stroke-width",1.5);
  svg.append("path").attr("d", `M${bkR-4},${bkT} L${bkR},${bkT} L${bkR},${bkB} L${bkR-4},${bkB}`)
    .attr("fill","none").attr("stroke","#999").attr("stroke-width",1.5);

  // Brackets for b
  const bbL = bCol - 4, bbR = bCol + 58;
  svg.append("path").attr("d", `M${bbL+4},${bkT} L${bbL},${bkT} L${bbL},${bkB} L${bbL+4},${bkB}`)
    .attr("fill","none").attr("stroke","#999").attr("stroke-width",1.5);
  svg.append("path").attr("d", `M${bbR-4},${bkT} L${bbR},${bkT} L${bbR},${bkB} L${bbR-4},${bkB}`)
    .attr("fill","none").attr("stroke","#999").attr("stroke-width",1.5);

  // Source column header
  svg.append("text").attr("x", matL - 30).attr("y", matT - 8).attr("text-anchor","middle").attr("font-size",10).attr("fill","#999").text("src");

  // Draw rows
  for (let r = 0; r < nRows; r++) {
    const row = rows[r];
    const y = matT + r*cellH;

    // Row background stripe
    svg.append("rect").attr("x", matL).attr("y", y).attr("width", nCols*cellW).attr("height", cellH)
      .attr("fill", row.color).attr("opacity", 0.06);

    // Source tag
    svg.append("text").attr("x", matL - 30).attr("y", y + cellH/2 + 4)
      .attr("text-anchor","middle").attr("font-size", 9).attr("fill", row.color).attr("font-weight","bold")
      .text(row.src);

    // A cells
    for (let c = 0; c < nCols; c++) {
      const v = row.An[c];
      svg.append("rect").attr("x", matL + c*cellW).attr("y", y).attr("width", cellW).attr("height", cellH)
        .attr("fill","none").attr("stroke","#e8e8e8").attr("stroke-width",0.5);
      const txt = fmt(v);
      svg.append("text")
        .attr("x", matL + c*cellW + cellW/2).attr("y", y + cellH/2 + 4)
        .attr("text-anchor","middle")
        .attr("font-size", txt.length > 5 ? 8 : (txt.length > 3 ? 9 : 10))
        .attr("fill", v !== 0 ? row.color : "#e0e0e0")
        .text(txt);
    }

    // b cell
    svg.append("rect").attr("x", bCol).attr("y", y).attr("width", 55).attr("height", cellH)
      .attr("fill", row.color).attr("opacity", 0.06);
    svg.append("rect").attr("x", bCol).attr("y", y).attr("width", 55).attr("height", cellH)
      .attr("fill","none").attr("stroke","#e8e8e8").attr("stroke-width",0.5);
    const btxt = fmt(row.bn);
    svg.append("text").attr("x", bCol + 27).attr("y", y + cellH/2 + 4)
      .attr("text-anchor","middle").attr("font-size", btxt.length > 5 ? 8 : 10)
      .attr("fill", row.color)
      .text(btxt);

    // Note
    svg.append("text").attr("x", noteCol).attr("y", y + cellH/2 + 4)
      .attr("font-size", 9).attr("fill","#aaa").text(row.note);
  }

  // Legend
  const legY = matT + nRows*cellH + 25;
  const items = [
    {color: smoothColor, tag: "E", label: `Smoothness: √π · (Xᵢ − Xⱼ) = 0 for each edge, each component (√π ≈ ${sp.toFixed(3)})`},
    {color: bdryColor, tag: "BC", label: `Boundary (Eq. 2): λ · Xᵢ = λ · cos(4θⁿ) and λ · Xᵢ₊₁ = λ · sin(4θⁿ) for locked vertices`},
    {color: feasColor, tag: "F", label: "Feasibility (smoothing only, not shown): λ · (ã₀Xᵢ + ã₁Xᵢ₊₁) = λ, using prev. iteration's ã"},
  ];
  for (let k = 0; k < items.length; k++) {
    const iy = legY + k * 22;
    svg.append("rect").attr("x", 15).attr("y", iy - 8).attr("width", 14).attr("height", 14)
      .attr("fill", items[k].color).attr("opacity", 0.15).attr("stroke", items[k].color).attr("stroke-width",1);
    svg.append("text").attr("x", 22).attr("y", iy + 2).attr("text-anchor","middle").attr("font-size",9)
      .attr("fill", items[k].color).attr("font-weight","bold").text(items[k].tag);
    svg.append("text").attr("x", 38).attr("y", iy + 2).attr("font-size", 10).attr("fill","#555").text(items[k].label);
  }

  // Annotation: what changes with θ
  svg.append("text").attr("x", W/2).attr("y", legY + 80).attr("text-anchor","middle").attr("font-size", 11).attr("fill","#666")
    .text("Drag the sliders — θ changes the orange b values, λ scales the orange rows vs. the blue rows.");

  return svg.node();
}

Solve \(A^\top A X = A^\top b\), then project: \(a^i \leftarrow (X_{2i}, X_{2i+1})^\top / \|(X_{2i}, X_{2i+1})\|\).

Smoothing Iterations

“Each smoothing iteration is similar to the initialization problem, except that we add to the objective function a new quadratic penalty term that corresponds to a linear approximation of the feasibility constraint as done in [Kowalski et al. 2012, p. 6].”

“As before it is expressed by a new set of linear equations when constructing \(A\) and \(b\): for each vertex \(i\), we add one equation \(\lambda(X_{2i}a_0^i + X_{2i+1}a_1^i - 1) = 0\), where \(a^i\) denotes the solution obtained at the previous iteration.”

Below is the full solve for the toy problem. Use the slider to step through smoothing iterations and watch the coefficient vectors converge on the unit circle:

viewof impl_t0 = Inputs.range([0, 90], {value: 0, step: 1, label: "θ₀ locked (degrees)"})
viewof impl_t3 = Inputs.range([0, 90], {value: 60, step: 1, label: "θ₃ locked (degrees)"})
viewof impl_lambda_bdry = Inputs.range([1, 200], {value: 100, step: 1, label: "λ boundary"})
viewof impl_lambda_feas = Inputs.range([1, 200], {value: 100, step: 1, label: "λ feasibility"})
viewof n_smooth = Inputs.range([0, 10], {value: 0, step: 1, label: "Smoothing iterations"})

(a)

(b)

(c)

(d)

(e)

Figure 9

{
  const t0 = impl_t0*Math.PI/180, t3 = impl_t3*Math.PI/180;
  const nv = 4, lam = impl_lambda_bdry, lam_f = impl_lambda_feas;
  const edgeList = [[0,1],[1,2],[2,3]];
  const bdry = {0: t0, 3: t3};

  // Simple 8x8 linear system solver (toy size, use direct Gaussian elimination)
  function solveLS(n, getRowsFn) {
    // Build AᵀA and Aᵀb from list of (row_coeffs, rhs)
    const AtA = Array.from({length: n}, () => new Float64Array(n));
    const Atb = new Float64Array(n);
    for (const {coeffs, rhs} of getRowsFn()) {
      for (let i = 0; i < n; i++) {
        for (let j = 0; j < n; j++) {
          AtA[i][j] += coeffs[i] * coeffs[j];
        }
        Atb[i] += coeffs[i] * rhs;
      }
    }
    // Gaussian elimination
    const A = AtA.map(r => [...r]);
    const b = [...Atb];
    for (let col = 0; col < n; col++) {
      let maxR = col;
      for (let r = col+1; r < n; r++) if (Math.abs(A[r][col]) > Math.abs(A[maxR][col])) maxR = r;
      [A[col], A[maxR]] = [A[maxR], A[col]];
      [b[col], b[maxR]] = [b[maxR], b[col]];
      if (Math.abs(A[col][col]) < 1e-15) continue;
      for (let r = col+1; r < n; r++) {
        const f = A[r][col]/A[col][col];
        for (let c = col; c < n; c++) A[r][c] -= f*A[col][c];
        b[r] -= f*b[col];
      }
    }
    const x = new Float64Array(n);
    for (let i = n-1; i >= 0; i--) {
      let s = b[i];
      for (let j = i+1; j < n; j++) s -= A[i][j]*x[j];
      x[i] = s / A[i][i];
    }
    return x;
  }

  function buildRows(a_prev) {
    return function*() {
      const sqpi = Math.sqrt(Math.PI);
      // Smoothness
      for (const [vi,vj] of edgeList) {
        for (let comp = 0; comp < 2; comp++) {
          const c = new Float64Array(2*nv);
          c[2*vi+comp] = sqpi; c[2*vj+comp] = -sqpi;
          yield {coeffs: c, rhs: 0};
        }
      }
      // Boundary
      for (const [vi_s, ti] of Object.entries(bdry)) {
        const vi = parseInt(vi_s);
        const c0 = new Float64Array(2*nv); c0[2*vi] = lam;
        yield {coeffs: c0, rhs: lam*Math.cos(4*ti)};
        const c1 = new Float64Array(2*nv); c1[2*vi+1] = lam;
        yield {coeffs: c1, rhs: lam*Math.sin(4*ti)};
      }
      // Feasibility linearization (if a_prev provided)
      if (a_prev) {
        for (let vi = 0; vi < nv; vi++) {
          const c = new Float64Array(2*nv);
          c[2*vi] = lam_f*a_prev[vi][0]; c[2*vi+1] = lam_f*a_prev[vi][1];
          yield {coeffs: c, rhs: lam_f};
        }
      }
    };
  }

  // Initialization (no feasibility constraint)
  const X0 = solveLS(2*nv, buildRows(null));
  const a_relaxed = [];
  for (let i = 0; i < nv; i++) a_relaxed.push([X0[2*i], X0[2*i+1]]);

  // Project
  let a_current = a_relaxed.map(a => {
    const n = Math.sqrt(a[0]**2+a[1]**2);
    return n > 1e-10 ? [a[0]/n, a[1]/n] : [1, 0];
  });
  const a_after_proj = a_current.map(a => [...a]);

  // Smoothing iterations
  for (let it = 0; it < n_smooth; it++) {
    const X = solveLS(2*nv, buildRows(a_current));
    a_current = [];
    for (let i = 0; i < nv; i++) {
      const v = [X[2*i], X[2*i+1]];
      const n = Math.sqrt(v[0]**2+v[1]**2);
      a_current.push(n > 1e-10 ? [v[0]/n, v[1]/n] : [1,0]);
    }
  }

  // --- Draw ---
  const W = 700, H = 400;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");
  const colrs = ["orange","steelblue","forestgreen","orange"];
  const lbls = ["P₀","P₁","P₂","P₃"];

  // --- Left: coefficient space ---
  const lcx = 180, lcy = H/2, lr = 150;
  svg.append("text").attr("x",lcx).attr("y",18).attr("text-anchor","middle").attr("font-size",14).attr("fill","#333")
    .text(n_smooth === 0 ? "After initialization + projection" : `After ${n_smooth} smoothing iteration(s)`);

  const circ = d3.range(0, 2*Math.PI+0.01, 0.02).map(a => [lcx+lr*Math.cos(a), lcy-lr*Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circ)).attr("fill","none").attr("stroke","#ddd").attr("stroke-width",2);

  // Show relaxed points (diamonds) if iter=0
  if (n_smooth === 0) {
    for (let i = 0; i < nv; i++) {
      const rx = lcx+lr*a_relaxed[i][0], ry = lcy-lr*a_relaxed[i][1];
      svg.append("rect").attr("x",rx-5).attr("y",ry-5).attr("width",10).attr("height",10)
        .attr("fill","none").attr("stroke",colrs[i]).attr("stroke-width",1.5)
        .attr("transform",`rotate(45,${rx},${ry})`);
    }
    // chord line
    svg.append("line")
      .attr("x1",lcx+lr*a_relaxed[0][0]).attr("y1",lcy-lr*a_relaxed[0][1])
      .attr("x2",lcx+lr*a_relaxed[3][0]).attr("y2",lcy-lr*a_relaxed[3][1])
      .attr("stroke","red").attr("stroke-width",1).attr("stroke-dasharray","4,3");
  }

  // Edges in coeff space
  for (const [vi,vj] of edgeList) {
    svg.append("line")
      .attr("x1",lcx+lr*a_current[vi][0]).attr("y1",lcy-lr*a_current[vi][1])
      .attr("x2",lcx+lr*a_current[vj][0]).attr("y2",lcy-lr*a_current[vj][1])
      .attr("stroke","#aaa").attr("stroke-width",1.5);
  }

  // Points
  for (let i = 0; i < nv; i++) {
    svg.append("circle").attr("cx",lcx+lr*a_current[i][0]).attr("cy",lcy-lr*a_current[i][1]).attr("r",7).attr("fill",colrs[i]);
    svg.append("text").attr("x",lcx+lr*a_current[i][0]+10).attr("y",lcy-lr*a_current[i][1]-10).attr("font-size",13).attr("fill",colrs[i]).text(lbls[i]);
  }

  // --- Right: frame chain ---
  const chainX0 = 390, chainSpacing = 80, chainY = H/2;
  svg.append("text").attr("x",chainX0+1.5*chainSpacing).attr("y",18).attr("text-anchor","middle").attr("font-size",14).attr("fill","#333").text("Resulting frames");

  svg.append("line").attr("x1",chainX0).attr("y1",chainY).attr("x2",chainX0+3*chainSpacing).attr("y2",chainY).attr("stroke","#bbb").attr("stroke-width",2);

  for (let i = 0; i < nv; i++) {
    const cx = chainX0 + i*chainSpacing;
    const th = Math.atan2(a_current[i][1], a_current[i][0]) / 4;
    for (let k = 0; k < 4; k++) {
      const a = th + k*Math.PI/2;
      svg.append("line").attr("x1",cx).attr("y1",chainY).attr("x2",cx+30*Math.cos(a)).attr("y2",chainY-30*Math.sin(a))
        .attr("stroke",colrs[i]).attr("stroke-width",2);
    }
    svg.append("circle").attr("cx",cx).attr("cy",chainY).attr("r",4).attr("fill",colrs[i]);
    svg.append("text").attr("x",cx).attr("y",chainY+50).attr("text-anchor","middle").attr("font-size",12).attr("fill",colrs[i]).text(lbls[i]);
  }

  // Energy display
  let objE = 0;
  for (const [vi,vj] of edgeList) {
    objE += (a_current[vj][0]-a_current[vi][0])**2 + (a_current[vj][1]-a_current[vi][1])**2;
  }
  objE *= Math.PI;
  svg.append("text").attr("x",chainX0).attr("y",H-30).attr("font-size",13).attr("fill","teal").text(`E = ${objE.toFixed(4)}`);

  // Legend for initialization
  if (n_smooth === 0) {
    svg.append("rect").attr("x",20).attr("y",H-35).attr("width",8).attr("height",8).attr("fill","none").attr("stroke","#888").attr("transform","rotate(45,24," + (H-31) + ")");
    svg.append("text").attr("x",35).attr("y",H-27).attr("font-size",11).attr("fill","#888").text("Relaxed (off circle)");
    svg.append("circle").attr("cx",170).attr("cy",H-31).attr("r",5).attr("fill","#888");
    svg.append("text").attr("x",180).attr("y",H-27).attr("font-size",11).attr("fill","#888").text("Projected (on circle)");
  }

  return svg.node();
}

---

2.7 Toy Problem Revisited

“This section explains our optimization approach on the toy problem already presented in § 2.5. As we mentioned before, at the initialization step we relax the constraints of feasibility of \(a^i\).”

“Top row of Figure 6 shows the solution of the initialization stage. Intuitively, we allow the points \(a^i\) not to be on the unit circle. Hence \(a^1\) and \(a^2\) lie on the chord between locked points \(a^0\) and \(a^3\).”

“Note that the initialization stage produces the correct sense of rotation (Figure 5). However it does not directly produce the optimum point \(P_0\). The reason being that the initialization stage produces points \(a^1\) and \(a^2\) (before normalization) in the way that all three chord segments are equal: \(\|a_0 - a_1\| = \|a_1 - a_2\| = \|a_2 - a_3\|\). But after projecting the points onto the feasible circle, the corresponding arc lengths are not equal. Therefore, we need a few smoothing iterations to reach the optimum.”

The interactive figure in §2.6 above demonstrates this: set smoothing iterations to 0 to see the initialization, then increase to watch convergence.

---

Interactive Sandbox: Variable-Length Chain

The toy problem in §2.5 fixed 4 vertices with 2 free. Here you can control how many interior frames there are, and lock or unlock any frame — including the endpoints.

• Click any frame cross to toggle its lock. Locked frames show a dashed orange ring and 🔒.
• Drag a locked frame to rotate it — the cross follows your cursor, and the solver re-runs when you release.
• Watch how the solver redistributes curvature when you pin frames in unexpected orientations.

viewof sb_n_inner = Inputs.range([1, 12], {value: 4, step: 1, label: "Interior frames"})
viewof sb_smooth = Inputs.range([0, 15], {value: 3, step: 1, label: "Smoothing iterations"})
viewof sb_lambda = Inputs.range([1, 200], {value: 100, step: 1, label: "λ (penalty weight)"})

(a)

(b)

(c)

Figure 10

// Persistent state: vertex index → angle in degrees (locked), or absent (free).
// Keys "L" and "R" control the two endpoints. Integer keys control interior vertices.
// Endpoints start locked; interior vertices start free.
mutable sb_locks = ({L: 0, R: 60})

Figure 11

{
  const nInner = sb_n_inner;
  const nv = nInner + 2;
  const lam = sb_lambda;
  const nSmooth = sb_smooth;
  const pi = Math.PI;
  const locks = sb_locks;

  // Edge list
  const edgeList = [];
  for (let i = 0; i < nv - 1; i++) edgeList.push([i, i + 1]);

  // Helper: lock-state key for vertex index
  function lockKey(i) { return i === 0 ? "L" : i === nv - 1 ? "R" : i; }

  // --- Solver (reusable — takes a boundary map) ---
  function solveLS(n, rowIter) {
    const AtA = Array.from({length: n}, () => new Float64Array(n));
    const Atb = new Float64Array(n);
    for (const {coeffs, rhs} of rowIter) {
      const nz = [];
      for (let i = 0; i < n; i++) if (coeffs[i] !== 0) nz.push(i);
      for (const i of nz) {
        for (const j of nz) AtA[i][j] += coeffs[i] * coeffs[j];
        Atb[i] += coeffs[i] * rhs;
      }
    }
    const A = AtA.map(r => [...r]);
    const b = [...Atb];
    for (let col = 0; col < n; col++) {
      let maxR = col;
      for (let r = col + 1; r < n; r++) if (Math.abs(A[r][col]) > Math.abs(A[maxR][col])) maxR = r;
      [A[col], A[maxR]] = [A[maxR], A[col]];
      [b[col], b[maxR]] = [b[maxR], b[col]];
      if (Math.abs(A[col][col]) < 1e-15) continue;
      for (let r = col + 1; r < n; r++) {
        const f = A[r][col] / A[col][col];
        for (let c = col; c < n; c++) A[r][c] -= f * A[col][c];
        b[r] -= f * b[col];
      }
    }
    const x = new Float64Array(n);
    for (let i = n - 1; i >= 0; i--) {
      let s = b[i];
      for (let j = i + 1; j < n; j++) s -= A[i][j] * x[j];
      x[i] = s / A[i][i];
    }
    return x;
  }

  function* makeRows(bdryMap, a_prev) {
    const sqpi = Math.sqrt(pi);
    for (const [vi, vj] of edgeList) {
      for (let comp = 0; comp < 2; comp++) {
        const c = new Float64Array(2 * nv);
        c[2 * vi + comp] = sqpi; c[2 * vj + comp] = -sqpi;
        yield {coeffs: c, rhs: 0};
      }
    }
    for (const [vi, ti] of bdryMap) {
      const c0 = new Float64Array(2 * nv); c0[2 * vi] = lam;
      yield {coeffs: c0, rhs: lam * Math.cos(4 * ti)};
      const c1 = new Float64Array(2 * nv); c1[2 * vi + 1] = lam;
      yield {coeffs: c1, rhs: lam * Math.sin(4 * ti)};
    }
    if (a_prev) {
      for (let vi = 0; vi < nv; vi++) {
        const c = new Float64Array(2 * nv);
        c[2 * vi] = lam * a_prev[vi][0]; c[2 * vi + 1] = lam * a_prev[vi][1];
        yield {coeffs: c, rhs: lam};
      }
    }
  }

  function solveField(bdryMap, iters) {
    const X0 = solveLS(2 * nv, makeRows(bdryMap, null));
    let a = [];
    for (let i = 0; i < nv; i++) {
      const v0 = X0[2 * i], v1 = X0[2 * i + 1];
      const nm = Math.sqrt(v0 ** 2 + v1 ** 2);
      a.push(nm > 1e-10 ? [v0 / nm, v1 / nm] : [1, 0]);
    }
    for (let it = 0; it < iters; it++) {
      const X = solveLS(2 * nv, makeRows(bdryMap, a));
      a = [];
      for (let i = 0; i < nv; i++) {
        const v0 = X[2 * i], v1 = X[2 * i + 1];
        const nm = Math.sqrt(v0 ** 2 + v1 ** 2);
        a.push(nm > 1e-10 ? [v0 / nm, v1 / nm] : [1, 0]);
      }
    }
    return a;
  }

  function buildBdry(locksObj) {
    const b = new Map();
    if (typeof locksObj.L === "number") b.set(0, locksObj.L * pi / 180);
    if (typeof locksObj.R === "number") b.set(nv - 1, locksObj.R * pi / 180);
    for (let i = 1; i <= nInner; i++) {
      if (typeof locksObj[i] === "number") b.set(i, locksObj[i] * pi / 180);
    }
    return b;
  }

  function computeEnergy(a) {
    let E = 0;
    for (const [vi, vj] of edgeList) {
      E += (a[vj][0] - a[vi][0]) ** 2 + (a[vj][1] - a[vi][1]) ** 2;
    }
    return E * pi;
  }

  // Initial solve
  const bdry = buildBdry(locks);
  let a_current = solveField(bdry, nSmooth);

  // --- Visualization ---
  const W = 720, H = 380;
  const chainW = 420, coeffW = 260;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`)
    .style("width", "100%").style("max-width", W + "px");

  // --- Left: frame chain ---
  const chainMarginL = 40, chainMarginR = 30;
  const chainSpacing = (chainW - chainMarginL - chainMarginR) / Math.max(1, nv - 1);
  const chainY = H / 2;
  const frameR = Math.min(28, chainSpacing * 0.4);

  svg.append("text").attr("x", chainW / 2).attr("y", 18)
    .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", "#333")
    .text("Frame chain — click to lock/unlock, drag to rotate");

  svg.append("line")
    .attr("x1", chainMarginL).attr("y1", chainY)
    .attr("x2", chainMarginL + (nv - 1) * chainSpacing).attr("y2", chainY)
    .attr("stroke", "#bbb").attr("stroke-width", 2);

  // --- Right: coefficient space (static scaffold) ---
  const rcx = chainW + coeffW / 2 + 10, rcy = H / 2, rr = Math.min(140, coeffW / 2 - 20);
  svg.append("text").attr("x", rcx).attr("y", 18)
    .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", "#333").text("Coefficient space");
  const circ = d3.range(0, 2 * pi + 0.01, 0.02).map(a => [rcx + rr * Math.cos(a), rcy - rr * Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circ)).attr("fill", "none").attr("stroke", "#ddd").attr("stroke-width", 1.5);
  svg.append("line").attr("x1", rcx - rr - 10).attr("y1", rcy).attr("x2", rcx + rr + 10).attr("y2", rcy).attr("stroke", "#eee");
  svg.append("line").attr("x1", rcx).attr("y1", rcy - rr - 10).attr("x2", rcx).attr("y2", rcy + rr + 10).attr("stroke", "#eee");

  // --- Store element references for live updates ---
  const frameGs = [];
  const coeffEdgeLines = [];
  const coeffDots = [];
  const coeffLabels = [];

  // Coeff space edge lines
  for (let e = 0; e < edgeList.length; e++) {
    const [vi, vj] = edgeList[e];
    coeffEdgeLines.push(svg.append("line")
      .attr("x1", rcx + rr * a_current[vi][0]).attr("y1", rcy - rr * a_current[vi][1])
      .attr("x2", rcx + rr * a_current[vj][0]).attr("y2", rcy - rr * a_current[vj][1])
      .attr("stroke", "#aaa").attr("stroke-width", 1.5));
  }

  // Draw each vertex (chain frame + coeff dot)
  for (let i = 0; i < nv; i++) {
    const cx = chainMarginL + i * chainSpacing;
    const th = Math.atan2(a_current[i][1], a_current[i][0]) / 4;
    const isLocked = bdry.has(i);
    const color = isLocked ? "#d97706" : "steelblue";

    // Frame cross group
    const g = svg.append("g");
    for (let k = 0; k < 4; k++) {
      const ang = th + k * pi / 2;
      g.append("line")
        .attr("x1", cx).attr("y1", chainY)
        .attr("x2", cx + frameR * Math.cos(ang)).attr("y2", chainY - frameR * Math.sin(ang))
        .attr("stroke", color).attr("stroke-width", 2);
    }
    frameGs.push(g);

    // Lock ring
    if (isLocked) {
      svg.append("circle").attr("cx", cx).attr("cy", chainY).attr("r", frameR + 5)
        .attr("fill", "none").attr("stroke", "#d97706").attr("stroke-width", 1.5)
        .attr("stroke-dasharray", "4,3");
    }

    // Center dot
    svg.append("circle").attr("cx", cx).attr("cy", chainY).attr("r", 4).attr("fill", color);

    // Label
    const label = i === 0 ? "P₀" : i === nv - 1 ? `P${nv - 1}` : `P${i}`;
    const angleDeg = isLocked ? locks[lockKey(i)] : null;
    const labelText = isLocked ? `${label} 🔒 ${angleDeg}°` : label;
    svg.append("text").attr("x", cx).attr("y", chainY + frameR + 18)
      .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", color)
      .text(labelText);

    // Coeff space dot + label
    coeffDots.push(svg.append("circle")
      .attr("cx", rcx + rr * a_current[i][0]).attr("cy", rcy - rr * a_current[i][1])
      .attr("r", 6).attr("fill", color));
    coeffLabels.push(svg.append("text")
      .attr("x", rcx + rr * a_current[i][0] + 9).attr("y", rcy - rr * a_current[i][1] - 9)
      .attr("font-size", 11).attr("fill", color).text(label));

    // --- Interaction: click to toggle lock, drag to rotate ---
    const hitTarget = svg.append("circle").attr("cx", cx).attr("cy", chainY).attr("r", frameR + 8)
      .attr("fill", "transparent").attr("cursor", "pointer");

    let isDragging = false;

    const drag = d3.drag()
      .on("start", () => { isDragging = false; })
      .on("drag", (event) => {
        if (!isLocked) return;
        isDragging = true;

        // Compute dragged angle
        const dx = event.x - cx, dy = -(event.y - chainY);
        const dragAngle = Math.atan2(dy, dx);
        const deg90 = ((dragAngle * 180 / pi) % 90 + 90) % 90;

        // Build temp boundary and re-solve
        const tempLocks = {...locks};
        tempLocks[lockKey(i)] = deg90;
        const tempBdry = buildBdry(tempLocks);
        const a_solved = solveField(tempBdry, nSmooth);

        // Update ALL frame crosses
        for (let v = 0; v < nv; v++) {
          const vx = chainMarginL + v * chainSpacing;
          const vth = Math.atan2(a_solved[v][1], a_solved[v][0]) / 4;
          const vc = tempBdry.has(v) ? "#d97706" : "steelblue";
          frameGs[v].selectAll("line").remove();
          for (let k = 0; k < 4; k++) {
            const ang = vth + k * pi / 2;
            frameGs[v].append("line")
              .attr("x1", vx).attr("y1", chainY)
              .attr("x2", vx + frameR * Math.cos(ang)).attr("y2", chainY - frameR * Math.sin(ang))
              .attr("stroke", vc).attr("stroke-width", 2);
          }
          // Update coeff space
          coeffDots[v].attr("cx", rcx + rr * a_solved[v][0]).attr("cy", rcy - rr * a_solved[v][1]);
          coeffLabels[v].attr("x", rcx + rr * a_solved[v][0] + 9).attr("y", rcy - rr * a_solved[v][1] - 9);
        }
        for (let e = 0; e < edgeList.length; e++) {
          const [vi, vj] = edgeList[e];
          coeffEdgeLines[e]
            .attr("x1", rcx + rr * a_solved[vi][0]).attr("y1", rcy - rr * a_solved[vi][1])
            .attr("x2", rcx + rr * a_solved[vj][0]).attr("y2", rcy - rr * a_solved[vj][1]);
        }
        energyText.text(`E = ${computeEnergy(a_solved).toFixed(4)}`);
      })
      .on("end", (event) => {
        if (isDragging && isLocked) {
          const dx = event.x - cx, dy = -(event.y - chainY);
          const dragAngle = Math.atan2(dy, dx);
          const deg90 = Math.round(((dragAngle * 180 / pi) % 90 + 90) % 90);
          const newLocks = {...mutable sb_locks};
          newLocks[lockKey(i)] = deg90;
          mutable sb_locks = newLocks;
        } else if (!isDragging) {
          const newLocks = {...mutable sb_locks};
          const key = lockKey(i);
          if (typeof newLocks[key] === "number") {
            delete newLocks[key];
          } else {
            const angle = Math.atan2(a_current[i][1], a_current[i][0]) / 4 * 180 / pi;
            newLocks[key] = Math.round(((angle % 90) + 90) % 90);
          }
          mutable sb_locks = newLocks;
        }
      });

    hitTarget.call(drag);
  }

  // --- Stats ---
  const nLocked = bdry.size, nFree = nv - nLocked;
  const energyText = svg.append("text").attr("x", chainW + 15).attr("y", H - 35)
    .attr("font-size", 12).attr("fill", "teal")
    .text(`E = ${computeEnergy(a_current).toFixed(4)}`);
  svg.append("text").attr("x", chainW + 15).attr("y", H - 18)
    .attr("font-size", 11).attr("fill", "#888")
    .text(`${nv} vertices (${nLocked} locked, ${nFree} free)`);

  return svg.node();
}

Figure 12

What to try:

• Start with all interior frames free — the solver evenly distributes the rotation between endpoints.
• Unlock an endpoint (click it) — with no boundary conditions on that side, the solver extends the field freely.
• Lock a middle frame at 0° while the right endpoint is at 60° — curvature concentrates near the locked frame.
• Drag a locked frame to rotate it and watch the energy change in real time.
• Unlock everything — with zero constraints the solver has no preferred orientation (the energy is trivially zero).

---

2.8 Boundary constraints: the parallelogram trap (cf. paper Fig. 8)

“Boundary constraint. Using it for the initialization step is wrong: for example, a normal constraint of angle \(\theta = 0\) forces \(a_0 = 1\) but let \(a_1\) free. As illustrated in Figure 8, it can produce very bad frame fields. Therefore we must use two separate equations (2).”

“Boundary constraints for global optimization: if we only use one constraint (eq. 3), the initialization step finds a constant frame field on a parallelogram (left). It is therefore mandatory to lock the frames to get the proper boundary constraints (right).” (Fig. 8 caption)

This is a subtle but critical point. During the initialization step, we relax the feasibility constraint \(a^\top a = 1\). Without feasibility, Eq. 3 (\(a_0 \cos 4\theta^n + a_1 \sin 4\theta^n = 1\)) is a line in coefficient space, not a point. The optimizer is free to slide along that line — and it will, choosing whatever minimizes the smoothness energy.

On a rectangle this is invisible (all normals are axis-aligned, so everything maps to the same point anyway). A parallelogram makes the problem obvious: the slanted sides have normals that genuinely differ from horizontal/vertical, and the frames must rotate to follow them. With Eq. 3, the optimizer ignores this and produces a nearly constant field.

The fix: use both equations from Eq. 2 (\(a_0 = \cos 4\theta^n\), \(a_1 = \sin 4\theta^n\)) during initialization. This pins each boundary vertex to an exact point in coefficient space.

Drag the skew angle to tilt the parallelogram. Toggle between constraint types to see the difference:

viewof skew_deg = Inputs.range([0, 40], {value: 25, step: 1, label: "Skew angle (degrees)"})
viewof grid_nx = Inputs.range([3, 12], {value: 6, step: 1, label: "Grid columns"})
viewof grid_ny = Inputs.range([3, 10], {value: 5, step: 1, label: "Grid rows"})
viewof constraint_mode = Inputs.radio(["Eq. 2: both components (correct)", "Eq. 3: single equation (wrong)"], {value: "Eq. 2: both components (correct)", label: "Boundary constraint"})
viewof para_lambda_bdry = Inputs.range([1, 200], {value: 100, step: 1, label: "λ boundary"})
viewof para_lambda_feas = Inputs.range([1, 200], {value: 100, step: 1, label: "λ feasibility"})
viewof para_smooth = Inputs.range([0, 10], {value: 0, step: 1, label: "Smoothing iterations"})

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 13

{
  const skew = skew_deg * Math.PI / 180;

  // Build a parallelogram mesh: nx×ny grid, sheared by skew angle
  const nx = grid_nx, ny = grid_ny;
  const nv = nx * ny;
  const baseW = 3, baseH = 2;
  const verts = [];
  for (let j = 0; j < ny; j++) {
    const ty = j / (ny - 1);
    const shear = ty * baseH * Math.tan(skew); // horizontal shift increases with row
    for (let i = 0; i < nx; i++) {
      const tx = i / (nx - 1);
      verts.push([tx * baseW + shear, ty * baseH]);
    }
  }

  // Edges: horizontal + vertical (grid connectivity)
  const edges = [];
  for (let j = 0; j < ny; j++)
    for (let i = 0; i < nx - 1; i++)
      edges.push([j * nx + i, j * nx + i + 1]);
  for (let j = 0; j < ny - 1; j++)
    for (let i = 0; i < nx; i++)
      edges.push([j * nx + i, (j + 1) * nx + i]);

  // Boundary normals: computed from edge directions
  // Bottom: normal = -ŷ = 3π/2
  // Top: normal = +ŷ = π/2
  // Right: edge goes from (baseW, 0) to (baseW + baseH·tan(skew), baseH)
  //   direction = (tan(skew), 1), outward normal = rotate -90° = (1, -tan(skew)) → angle
  // Left: edge goes from (0,0) to (baseH·tan(skew), baseH)
  //   direction = (tan(skew), 1), outward normal = rotate +90° = (-1, tan(skew)) → angle
  const rightNormal = Math.atan2(-Math.tan(skew), 1); // points right-and-down
  const leftNormal = Math.atan2(Math.tan(skew), -1);   // points left-and-up

  const bdryMap = {};
  // Bottom row
  for (let i = 0; i < nx; i++) bdryMap[i] = 3 * Math.PI / 2;
  // Top row
  for (let i = 0; i < nx; i++) bdryMap[(ny - 1) * nx + i] = Math.PI / 2;
  // Left column (not corners, already set)
  for (let j = 1; j < ny - 1; j++) bdryMap[j * nx] = leftNormal;
  // Right column
  for (let j = 1; j < ny - 1; j++) bdryMap[j * nx + nx - 1] = rightNormal;
  // Corners: at a skewed parallelogram, the two adjacent normals don't
  // differ by exactly π/2, so no single frame can perfectly align with both.
  // We assign each corner the normal of its horizontal edge (top/bottom).
  // The smoothness energy will handle the transition to the slanted sides.
  bdryMap[0] = 3 * Math.PI / 2;                     // bottom-left → bottom normal
  bdryMap[nx-1] = 3 * Math.PI / 2;                  // bottom-right → bottom normal
  bdryMap[(ny-1)*nx] = Math.PI / 2;                  // top-left → top normal
  bdryMap[(ny-1)*nx + nx-1] = Math.PI / 2;           // top-right → top normal

  const useEq2 = constraint_mode.startsWith("Eq. 2");
  const lam = para_lambda_bdry, lam_f = para_lambda_feas, sqpi = Math.sqrt(Math.PI);

  // Gaussian elimination solver
  function solveSystem(n, rowGen) {
    const AtA = Array.from({length: n}, () => new Float64Array(n));
    const Atb = new Float64Array(n);
    for (const {coeffs, rhs} of rowGen()) {
      for (let ii = 0; ii < n; ii++) {
        for (let jj = 0; jj < n; jj++) AtA[ii][jj] += coeffs[ii] * coeffs[jj];
        Atb[ii] += coeffs[ii] * rhs;
      }
    }
    const A = AtA.map(r => [...r]);
    const b = [...Atb];
    for (let col = 0; col < n; col++) {
      let maxR = col;
      for (let r = col+1; r < n; r++) if (Math.abs(A[r][col]) > Math.abs(A[maxR][col])) maxR = r;
      [A[col], A[maxR]] = [A[maxR], A[col]];
      [b[col], b[maxR]] = [b[maxR], b[col]];
      if (Math.abs(A[col][col]) < 1e-15) continue;
      for (let r = col+1; r < n; r++) {
        const f = A[r][col] / A[col][col];
        for (let c = col; c < n; c++) A[r][c] -= f * A[col][c];
        b[r] -= f * b[col];
      }
    }
    const x = new Float64Array(n);
    for (let ii = n-1; ii >= 0; ii--) {
      let s = b[ii];
      for (let jj = ii+1; jj < n; jj++) s -= A[ii][jj] * x[jj];
      x[ii] = s / A[ii][ii];
    }
    return x;
  }

  function buildRows(a_prev) {
    return function*() {
      // Smoothness
      for (const [vi, vj] of edges) {
        for (let comp = 0; comp < 2; comp++) {
          const c = new Float64Array(2*nv);
          c[2*vi+comp] = sqpi; c[2*vj+comp] = -sqpi;
          yield {coeffs: c, rhs: 0};
        }
      }
      // Boundary constraints
      for (const [vi_s, theta_n] of Object.entries(bdryMap)) {
        const vi = parseInt(vi_s);
        if (useEq2) {
          const c0 = new Float64Array(2*nv); c0[2*vi] = lam;
          yield {coeffs: c0, rhs: lam * Math.cos(4*theta_n)};
          const c1 = new Float64Array(2*nv); c1[2*vi+1] = lam;
          yield {coeffs: c1, rhs: lam * Math.sin(4*theta_n)};
        } else {
          const c = new Float64Array(2*nv);
          c[2*vi] = lam * Math.cos(4*theta_n);
          c[2*vi+1] = lam * Math.sin(4*theta_n);
          yield {coeffs: c, rhs: lam};
        }
      }
      if (a_prev) {
        for (let vi = 0; vi < nv; vi++) {
          const c = new Float64Array(2*nv);
          c[2*vi] = lam_f * a_prev[vi][0];
          c[2*vi+1] = lam_f * a_prev[vi][1];
          yield {coeffs: c, rhs: lam_f};
        }
      }
    };
  }

  // Solve: initialization + projection
  const X0 = solveSystem(2*nv, buildRows(null));
  let a_current = [];
  for (let i = 0; i < nv; i++) {
    const v = [X0[2*i], X0[2*i+1]];
    const nm = Math.sqrt(v[0]**2 + v[1]**2);
    a_current.push(nm > 1e-10 ? [v[0]/nm, v[1]/nm] : [1, 0]);
  }

  // Smoothing iterations
  for (let it = 0; it < para_smooth; it++) {
    const X = solveSystem(2*nv, buildRows(a_current));
    a_current = [];
    for (let i = 0; i < nv; i++) {
      const v = [X[2*i], X[2*i+1]];
      const nm = Math.sqrt(v[0]**2 + v[1]**2);
      a_current.push(nm > 1e-10 ? [v[0]/nm, v[1]/nm] : [1, 0]);
    }
  }

  // --- Draw ---
  const W = 700, H = 440;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  // Compute bounding box of vertices for scaling
  let minX = Infinity, maxX = -Infinity, minY = Infinity, maxY = -Infinity;
  for (const v of verts) {
    minX = Math.min(minX, v[0]); maxX = Math.max(maxX, v[0]);
    minY = Math.min(minY, v[1]); maxY = Math.max(maxY, v[1]);
  }
  const meshL = 40, meshR = 350, meshT = 50, meshB = 360;
  const meshW = meshR - meshL, meshH = meshB - meshT;
  const dataW = maxX - minX || 1, dataH = maxY - minY || 1;
  const scale = Math.min(meshW / dataW, meshH / dataH) * 0.85;
  const offX = meshL + (meshW - dataW * scale) / 2;
  const offY = meshT + (meshH - dataH * scale) / 2;
  function tx(v) { return offX + (v[0] - minX) * scale; }
  function ty(v) { return meshB - (offY - meshT) - (v[1] - minY) * scale; }

  svg.append("text").attr("x", (meshL+meshR)/2).attr("y", 20).attr("text-anchor","middle").attr("font-size",13).attr("fill","#333")
    .text(useEq2 ? "Eq. 2: frames follow the boundary" : "Eq. 3: frames ignore the slant!");

  // Draw edges
  for (const [vi, vj] of edges) {
    svg.append("line")
      .attr("x1", tx(verts[vi])).attr("y1", ty(verts[vi]))
      .attr("x2", tx(verts[vj])).attr("y2", ty(verts[vj]))
      .attr("stroke", "#ddd").attr("stroke-width", 0.8);
  }

  // Draw frames — scale cross size to half the minimum edge length so they never overlap
  let minEdgeLen = Infinity;
  for (const [vi, vj] of edges) {
    const dx = tx(verts[vi]) - tx(verts[vj]);
    const dy = ty(verts[vi]) - ty(verts[vj]);
    minEdgeLen = Math.min(minEdgeLen, Math.sqrt(dx*dx + dy*dy));
  }
  const frameLen = Math.max(4, minEdgeLen * 0.42);
  for (let vi = 0; vi < nv; vi++) {
    const px = tx(verts[vi]), py = ty(verts[vi]);
    const a = a_current[vi];
    const th = Math.atan2(a[1], a[0]) / 4;
    const isBdry = vi in bdryMap;

    for (let k = 0; k < 4; k++) {
      const ang = th + k * Math.PI/2;
      svg.append("line")
        .attr("x1", px).attr("y1", py)
        .attr("x2", px + frameLen*Math.cos(ang)).attr("y2", py - frameLen*Math.sin(ang))
        .attr("stroke", isBdry ? "orange" : "steelblue").attr("stroke-width", 1.5);
    }
    svg.append("circle").attr("cx", px).attr("cy", py).attr("r", 2).attr("fill", isBdry ? "orange" : "steelblue");
  }

  // Draw boundary normals
  for (const [vi_s, theta_n] of Object.entries(bdryMap)) {
    const vi = parseInt(vi_s);
    const px = tx(verts[vi]), py = ty(verts[vi]);
    svg.append("line")
      .attr("x1", px).attr("y1", py)
      .attr("x2", px + 12*Math.cos(theta_n)).attr("y2", py - 12*Math.sin(theta_n))
      .attr("stroke", "red").attr("stroke-width", 1).attr("stroke-dasharray", "2,2");
  }

  // --- Right panel: coefficient space ---
  const rcx = 540, rcy = H/2, rr = 150;
  svg.append("text").attr("x", rcx).attr("y", 20).attr("text-anchor","middle").attr("font-size",13).attr("fill","#333").text("Coefficient space");

  const circ = d3.range(0, 2*Math.PI+0.01, 0.02).map(a => [rcx+rr*Math.cos(a), rcy-rr*Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circ)).attr("fill","none").attr("stroke","#ddd").attr("stroke-width",1.5);
  svg.append("line").attr("x1",rcx-rr-10).attr("y1",rcy).attr("x2",rcx+rr+10).attr("y2",rcy).attr("stroke","#eee");
  svg.append("line").attr("x1",rcx).attr("y1",rcy-rr-10).attr("x2",rcx).attr("y2",rcy+rr+10).attr("stroke","#eee");
  svg.append("text").attr("x",rcx+rr+12).attr("y",rcy-5).attr("font-size",10).attr("fill","#999").text("a₀");
  svg.append("text").attr("x",rcx+5).attr("y",rcy-rr-8).attr("font-size",10).attr("fill","#999").text("a₁");

  // Mark where each distinct boundary normal maps on the circle
  const uniqueNormals = new Map();
  for (const [vi_s, theta_n] of Object.entries(bdryMap)) {
    const key = theta_n.toFixed(4);
    if (!uniqueNormals.has(key)) uniqueNormals.set(key, theta_n);
  }
  for (const [, theta_n] of uniqueNormals) {
    const ca = Math.cos(4*theta_n), sa = Math.sin(4*theta_n);
    svg.append("circle").attr("cx", rcx+rr*ca).attr("cy", rcy-rr*sa).attr("r", 5)
      .attr("fill","none").attr("stroke","red").attr("stroke-width",1.5);
    const ndeg = (theta_n * 180 / Math.PI).toFixed(0);
    svg.append("text").attr("x", rcx+rr*ca + 8).attr("y", rcy-rr*sa - 8)
      .attr("font-size", 9).attr("fill","red").text(`θⁿ=${ndeg}°`);

    // In Eq3 mode, draw the constraint line through this point
    if (!useEq2) {
      // Line: ca·a₀ + sa·a₁ = 1
      // Parametric: (ca, sa) + t·(-sa, ca)
      const pts = [];
      for (let t = -2; t <= 2; t += 0.01) {
        const px = ca + t*(-sa), py = sa + t*ca;
        const sx2 = rcx + rr*px, sy2 = rcy - rr*py;
        if (sx2 > rcx-rr-15 && sx2 < rcx+rr+15 && sy2 > rcy-rr-15 && sy2 < rcy+rr+15) pts.push([sx2, sy2]);
      }
      if (pts.length > 1)
        svg.append("path").attr("d", d3.line()(pts)).attr("fill","none").attr("stroke","red").attr("stroke-width",1).attr("stroke-dasharray","4,3").attr("opacity",0.6);
    }
  }

  // Plot all coefficient vectors
  for (let vi = 0; vi < nv; vi++) {
    const a = a_current[vi];
    const isBdry = vi in bdryMap;
    svg.append("circle")
      .attr("cx", rcx + rr*a[0]).attr("cy", rcy - rr*a[1]).attr("r", 3)
      .attr("fill", isBdry ? "orange" : "steelblue").attr("opacity", 0.8);
  }

  if (!useEq2) {
    svg.append("text").attr("x", rcx).attr("y", H - 8).attr("text-anchor","middle").attr("font-size",11).attr("fill","red")
      .text("Dashed lines = Eq. 3 constraints (a₁ free along each line)");
  }

  // Legend
  svg.append("circle").attr("cx", meshL).attr("cy", H-15).attr("r", 4).attr("fill","orange");
  svg.append("text").attr("x", meshL+8).attr("y", H-11).attr("font-size",10).attr("fill","#888").text("boundary");
  svg.append("circle").attr("cx", meshL+80).attr("cy", H-15).attr("r", 4).attr("fill","steelblue");
  svg.append("text").attr("x", meshL+88).attr("y", H-11).attr("font-size",10).attr("fill","#888").text("interior");
  svg.append("line").attr("x1", meshL+140).attr("y1", H-15).attr("x2", meshL+155).attr("y2", H-15).attr("stroke","red").attr("stroke-dasharray","2,2");
  svg.append("text").attr("x", meshL+159).attr("y", H-11).attr("font-size",10).attr("fill","#888").text("normals");

  return svg.node();
}

WarningWhy does Eq. 3 fail on a parallelogram?

On a rectangle (skew = 0°), all four boundary normals (0°, 90°, 180°, 270°) map to the same point in coefficient space due to the \(\cos(4\theta)\) frequency — the problem is invisible.

A parallelogram breaks this degeneracy. The slanted sides have normals at angles \(\theta^n \neq k \cdot 90°\), which map to different points on the unit circle. Eq. 2 correctly pins each boundary vertex to its target point. But Eq. 3 only constrains each vertex to lie on a line through that point — and the optimizer slides along those lines to find a compromise that minimizes smoothness energy, producing a nearly constant field that ignores the slant.

Try it: set skew to 25°, toggle to Eq. 3, and observe how all frames point nearly the same direction. Then switch to Eq. 2 and watch the boundary frames properly rotate to align with the slanted sides. Add smoothing iterations to see how the interior frames interpolate.

---

Interactive Sandbox: Rectangle Grid

Now apply the same idea to a 2D grid. Boundary frames start locked to their outward normals. Click any frame to toggle its lock, drag a locked frame to rotate it.

viewof rect_nx = Inputs.range([3, 12], {value: 6, step: 1, label: "Columns"})
viewof rect_ny = Inputs.range([3, 10], {value: 5, step: 1, label: "Rows"})
viewof rect_smooth = Inputs.range([0, 15], {value: 3, step: 1, label: "Smoothing iterations"})
viewof rect_lambda = Inputs.range([1, 200], {value: 100, step: 1, label: "λ (penalty weight)"})

(a)

(b)

(c)

(d)

Figure 14

// Lock state for rectangle grid.
// Key = vertex index (string). Value = angle in degrees (locked) or absent (free).
// "bdry_default" flag: when true, boundary verts use their default normal.
// We store explicit overrides only.
mutable rect_locks = ({})

Figure 15

{
  const nx = rect_nx, ny = rect_ny;
  const nv = nx * ny;
  const lam = rect_lambda;
  const nSmooth = rect_smooth;
  const pi = Math.PI;
  const locks = rect_locks;

  // --- Grid vertex positions ---
  const baseW = 3, baseH = 2;
  const verts = [];
  for (let j = 0; j < ny; j++) {
    for (let i = 0; i < nx; i++) {
      verts.push([i / (nx - 1) * baseW, j / (ny - 1) * baseH]);
    }
  }

  // --- Edges: horizontal + vertical ---
  const edges = [];
  for (let j = 0; j < ny; j++)
    for (let i = 0; i < nx - 1; i++)
      edges.push([j * nx + i, j * nx + i + 1]);
  for (let j = 0; j < ny - 1; j++)
    for (let i = 0; i < nx; i++)
      edges.push([j * nx + i, (j + 1) * nx + i]);

  // --- Default boundary normals ---
  const defaultBdryAngle = {};
  for (let i = 0; i < nx; i++) defaultBdryAngle[i] = 3 * pi / 2;               // bottom
  for (let i = 0; i < nx; i++) defaultBdryAngle[(ny - 1) * nx + i] = pi / 2;    // top
  for (let j = 1; j < ny - 1; j++) defaultBdryAngle[j * nx] = pi;               // left
  for (let j = 1; j < ny - 1; j++) defaultBdryAngle[j * nx + nx - 1] = 0;       // right

  // --- Solver (sparse AtA assembly for performance during drag) ---
  const sqpi = Math.sqrt(pi);

  function solveSystem(n, rowIter) {
    const AtA = Array.from({length: n}, () => new Float64Array(n));
    const Atb = new Float64Array(n);
    for (const {coeffs, rhs} of rowIter) {
      const nz = [];
      for (let ii = 0; ii < n; ii++) if (coeffs[ii] !== 0) nz.push(ii);
      for (const ii of nz) {
        for (const jj of nz) AtA[ii][jj] += coeffs[ii] * coeffs[jj];
        Atb[ii] += coeffs[ii] * rhs;
      }
    }
    const A = AtA.map(r => [...r]);
    const b = [...Atb];
    for (let col = 0; col < n; col++) {
      let maxR = col;
      for (let r = col + 1; r < n; r++) if (Math.abs(A[r][col]) > Math.abs(A[maxR][col])) maxR = r;
      [A[col], A[maxR]] = [A[maxR], A[col]];
      [b[col], b[maxR]] = [b[maxR], b[col]];
      if (Math.abs(A[col][col]) < 1e-15) continue;
      for (let r = col + 1; r < n; r++) {
        const f = A[r][col] / A[col][col];
        for (let c = col; c < n; c++) A[r][c] -= f * A[col][c];
        b[r] -= f * b[col];
      }
    }
    const x = new Float64Array(n);
    for (let ii = n - 1; ii >= 0; ii--) {
      let s = b[ii];
      for (let jj = ii + 1; jj < n; jj++) s -= A[ii][jj] * x[jj];
      x[ii] = s / A[ii][ii];
    }
    return x;
  }

  function* makeRows(bdryMap, a_prev) {
    for (const [vi, vj] of edges) {
      for (let comp = 0; comp < 2; comp++) {
        const c = new Float64Array(2 * nv);
        c[2 * vi + comp] = sqpi; c[2 * vj + comp] = -sqpi;
        yield {coeffs: c, rhs: 0};
      }
    }
    for (const [vi, ti] of bdryMap) {
      const c0 = new Float64Array(2 * nv); c0[2 * vi] = lam;
      yield {coeffs: c0, rhs: lam * Math.cos(4 * ti)};
      const c1 = new Float64Array(2 * nv); c1[2 * vi + 1] = lam;
      yield {coeffs: c1, rhs: lam * Math.sin(4 * ti)};
    }
    if (a_prev) {
      for (let vi = 0; vi < nv; vi++) {
        const c = new Float64Array(2 * nv);
        c[2 * vi] = lam * a_prev[vi][0]; c[2 * vi + 1] = lam * a_prev[vi][1];
        yield {coeffs: c, rhs: lam};
      }
    }
  }

  function solveField(bdryMap, iters) {
    const X0 = solveSystem(2 * nv, makeRows(bdryMap, null));
    let a = [];
    for (let i = 0; i < nv; i++) {
      const v0 = X0[2 * i], v1 = X0[2 * i + 1];
      const nm = Math.sqrt(v0 ** 2 + v1 ** 2);
      a.push(nm > 1e-10 ? [v0 / nm, v1 / nm] : [1, 0]);
    }
    for (let it = 0; it < iters; it++) {
      const X = solveSystem(2 * nv, makeRows(bdryMap, a));
      a = [];
      for (let i = 0; i < nv; i++) {
        const v0 = X[2 * i], v1 = X[2 * i + 1];
        const nm = Math.sqrt(v0 ** 2 + v1 ** 2);
        a.push(nm > 1e-10 ? [v0 / nm, v1 / nm] : [1, 0]);
      }
    }
    return a;
  }

  function buildBdry(locksObj) {
    const b = new Map();
    for (let vi = 0; vi < nv; vi++) {
      const key = String(vi);
      const isBdryDefault = vi in defaultBdryAngle;
      if (locksObj[key] === "free") continue;
      else if (typeof locksObj[key] === "number") b.set(vi, locksObj[key] * pi / 180);
      else if (isBdryDefault) b.set(vi, defaultBdryAngle[vi]);
    }
    return b;
  }

  function computeEnergy(a) {
    let E = 0;
    for (const [vi, vj] of edges) {
      E += (a[vj][0] - a[vi][0]) ** 2 + (a[vj][1] - a[vi][1]) ** 2;
    }
    return E * pi;
  }

  // --- Initial solve ---
  const bdry = buildBdry(locks);
  let a_current = solveField(bdry, nSmooth);

  // Cap smoothing iterations during drag for responsiveness on large grids
  const dragSmooth = Math.min(nSmooth, 1);

  // --- Visualization ---
  const W = 700, H = 480;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`)
    .style("width", "100%").style("max-width", W + "px");

  // Scale grid to SVG
  const meshL = 40, meshR = W - 40, meshT = 50, meshB = H - 60;
  const meshW = meshR - meshL, meshH = meshB - meshT;
  const dataW = baseW, dataH = baseH;
  const scale = Math.min(meshW / dataW, meshH / dataH) * 0.85;
  const offX = meshL + (meshW - dataW * scale) / 2;
  const offY = meshT + (meshH + dataH * scale) / 2;
  function tx(v) { return offX + v[0] * scale; }
  function ty(v) { return offY - v[1] * scale; }

  svg.append("text").attr("x", W / 2).attr("y", 20)
    .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", "#333")
    .text("Rectangle grid — click to lock/unlock, drag to rotate");

  // Draw edges
  for (const [vi, vj] of edges) {
    svg.append("line")
      .attr("x1", tx(verts[vi])).attr("y1", ty(verts[vi]))
      .attr("x2", tx(verts[vj])).attr("y2", ty(verts[vj]))
      .attr("stroke", "#e8e8e8").attr("stroke-width", 0.8);
  }

  // Frame cross size
  let minEdgeLen = Infinity;
  for (const [vi, vj] of edges) {
    const dx = tx(verts[vi]) - tx(verts[vj]);
    const dy = ty(verts[vi]) - ty(verts[vj]);
    minEdgeLen = Math.min(minEdgeLen, Math.sqrt(dx * dx + dy * dy));
  }
  const frameLen = Math.max(4, minEdgeLen * 0.40);

  // --- Store frame group refs for live drag updates ---
  const frameGs = [];
  // Precompute screen positions
  const screenX = [], screenY = [];
  for (let vi = 0; vi < nv; vi++) {
    screenX.push(tx(verts[vi]));
    screenY.push(ty(verts[vi]));
  }

  // Draw frames with interaction
  for (let vi = 0; vi < nv; vi++) {
    const px = screenX[vi], py = screenY[vi];
    const a = a_current[vi];
    const th = Math.atan2(a[1], a[0]) / 4;
    const isLocked = bdry.has(vi);
    const color = isLocked ? "#d97706" : "steelblue";

    // Frame cross in a group
    const frameG = svg.append("g");
    for (let k = 0; k < 4; k++) {
      const ang = th + k * pi / 2;
      frameG.append("line")
        .attr("x1", px).attr("y1", py)
        .attr("x2", px + frameLen * Math.cos(ang)).attr("y2", py - frameLen * Math.sin(ang))
        .attr("stroke", color).attr("stroke-width", 1.5);
    }
    frameGs.push(frameG);

    // Lock ring for locked frames
    if (isLocked) {
      svg.append("circle").attr("cx", px).attr("cy", py).attr("r", frameLen + 3)
        .attr("fill", "none").attr("stroke", "#d97706").attr("stroke-width", 1)
        .attr("stroke-dasharray", "3,2").attr("opacity", 0.6);
    }

    // Center dot
    svg.append("circle").attr("cx", px).attr("cy", py).attr("r", 2).attr("fill", color);

    // --- Interaction: click to toggle, drag to rotate ---
    const hitR = Math.max(frameLen + 4, 12);
    const hitTarget = svg.append("circle").attr("cx", px).attr("cy", py).attr("r", hitR)
      .attr("fill", "transparent").attr("cursor", "pointer");

    let isDragging = false;

    const drag = d3.drag()
      .on("start", () => { isDragging = false; })
      .on("drag", (event) => {
        if (!isLocked) return;
        isDragging = true;
        const dx = event.x - px, dy = -(event.y - py);
        const dragAngle = Math.atan2(dy, dx);
        const deg90 = ((dragAngle * 180 / pi) % 90 + 90) % 90;

        // Build temp boundary and re-solve
        const tempLocks = {...locks};
        tempLocks[String(vi)] = deg90;
        const tempBdry = buildBdry(tempLocks);
        const a_solved = solveField(tempBdry, dragSmooth);

        // Update ALL frame crosses
        for (let v = 0; v < nv; v++) {
          const vth = Math.atan2(a_solved[v][1], a_solved[v][0]) / 4;
          const vc = tempBdry.has(v) ? "#d97706" : "steelblue";
          frameGs[v].selectAll("line").remove();
          for (let k = 0; k < 4; k++) {
            const ang = vth + k * pi / 2;
            frameGs[v].append("line")
              .attr("x1", screenX[v]).attr("y1", screenY[v])
              .attr("x2", screenX[v] + frameLen * Math.cos(ang))
              .attr("y2", screenY[v] - frameLen * Math.sin(ang))
              .attr("stroke", vc).attr("stroke-width", 1.5);
          }
        }
        energyText.text(`E = ${computeEnergy(a_solved).toFixed(4)}    ${nv} vertices    ${nSmooth} smoothing iter.`);
      })
      .on("end", (event) => {
        const key = String(vi);
        if (isDragging && isLocked) {
          const dx = event.x - px, dy = -(event.y - py);
          const dragAngle = Math.atan2(dy, dx);
          const deg90 = Math.round(((dragAngle * 180 / pi) % 90 + 90) % 90);
          const newLocks = {...mutable rect_locks};
          newLocks[key] = deg90;
          mutable rect_locks = newLocks;
        } else if (!isDragging) {
          const newLocks = {...mutable rect_locks};
          const isBdryDefault = vi in defaultBdryAngle;

          if (isLocked) {
            if (isBdryDefault && typeof newLocks[key] !== "number") {
              newLocks[key] = "free";
            } else {
              delete newLocks[key];
              if (isBdryDefault) newLocks[key] = "free";
            }
          } else {
            const angle = Math.atan2(a_current[vi][1], a_current[vi][0]) / 4 * 180 / pi;
            newLocks[key] = Math.round(((angle % 90) + 90) % 90);
          }
          mutable rect_locks = newLocks;
        }
      });

    hitTarget.call(drag);
  }

  // --- Stats ---
  const nLocked = bdry.size, nFree = nv - nLocked;
  const energyText = svg.append("text").attr("x", 40).attr("y", H - 15)
    .attr("font-size", 12).attr("fill", "teal")
    .text(`E = ${computeEnergy(a_current).toFixed(4)}    ${nv} vertices (${nLocked} locked, ${nFree} free)    ${nSmooth} smoothing iter.`);

  // Legend
  svg.append("circle").attr("cx", W - 220).attr("cy", H - 15).attr("r", 5).attr("fill", "#d97706");
  svg.append("text").attr("x", W - 210).attr("y", H - 11).attr("font-size", 10).attr("fill", "#888").text("locked");
  svg.append("circle").attr("cx", W - 160).attr("cy", H - 15).attr("r", 5).attr("fill", "steelblue");
  svg.append("text").attr("x", W - 150).attr("y", H - 11).attr("font-size", 10).attr("fill", "#888").text("free");

  // Reset hint
  svg.append("text").attr("x", W / 2).attr("y", H - 35)
    .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#bbb")
    .text("Change grid size to reset all locks");

  return svg.node();
}

Figure 16

What to try:

• The default rectangle has all boundary frames locked to outward normals. On a rectangle these are all axis-aligned, so the entire field is trivially uniform (\(E \approx 0\)).
• Unlock a boundary vertex — the solver smoothly fills in the gap. Unlock an entire side to see the field relax away from the boundary.
• Lock an interior vertex at 45° — this forces the surrounding field to rotate, creating a visible “ripple” of curvature around the locked frame.
• Lock two interior vertices at opposite angles — the solver must interpolate between them, concentrating curvature along the shortest path.
• Unlock all boundary vertices (click them all) and lock a single interior frame — the entire field aligns to that one constraint.

---

Summary

The key takeaways from Section 2:

Concept	Angle representation	Functional representation
Frame	Angle \(\theta\)	Coefficient vector \(a = (\cos 4\theta, \sin 4\theta)^\top\)
Feasibility	Always feasible	\(a^\top a = 1\) (unit circle)
Smoothness energy	\(\sum_{ij} (\Delta\theta_{ij})^2\) (piecewise quadratic)	\(\pi \sum_{ij} \|a^j - a^i\|^2\) (quadratic)
Boundary constraint	\(\theta^i = \theta^{\text{normal}}\)	\(a_0 \cos 4\theta^n + a_1 \sin 4\theta^n = 1\)
Optimization	Nonlinear (angle wrapping)	Linear system + projection

The functional approach converts a nonlinear optimization over angles into a sequence of linear least-squares problems followed by projection onto the unit circle.

---