--- title: "Integer-Grid Maps for Reliable Quad Meshing" subtitle: "Guaranteed Valid Parameterization via Mixed-Integer Optimization" --- ## Paper at a Glance ::: {.callout-note} ## Reference **Integer-Grid Maps for Reliable Quad Meshing** — Bommes, Campen, Ebke, Alliez, Kobbelt (2013) ACM Transactions on Graphics (SIGGRAPH), Vol. 32, No. 4, Article 98. [Download the original paper (PDF)](assets/bommes2013.pdf) ::: Parametrization-based quad meshing produces high-quality meshes but suffers from **degeneracies** — the Jacobian determinant of the piecewise linear map can go to zero or become negative, creating flipped triangles, collapsed edges, and local non-injectivities that prevent extracting a valid quad mesh. This paper defines the class of **Integer-Grid Maps (IGMs)** that are *guaranteed* to produce valid quad meshes, and derives a practical **Mixed-Integer Quadratic Program (MIQP)** formulation to find them. Two key optimizations — **complexity reduction** via parametric-space decimation and **singularity separation** via geodesic constraints — make the formulation tractable for real-world meshes. --- ## 1. Introduction ### 1.1 The Reliability Problem > *"The main source of non-reliable behavior are degeneracies in the constructed parametrization function, i.e. the Jacobian determinant of the piecewise linear map is locally less than or equal to zero."* The standard quad meshing pipeline works by mapping a triangle mesh onto a 2D parametric domain, then extracting the quad mesh from the integer iso-lines of that map. The mapping is piecewise linear — one affine map per triangle. When the **Jacobian determinant** of any triangle's map goes to zero or becomes negative, the triangle either collapses to a line or flips its orientation. These degeneracies destroy the essential property that the mapped integer grid forms a valid quad mesh. > *"Geometrically a value less than zero characterizes triangles that flip their orientation in the domain while a zero value results from local non-injectivities where an edge or even a whole triangle is mapped to a single point."* Existing approaches handle this with heuristics — greedy rounding of continuous solutions and stiffening of high-distortion regions — but these **fail for coarse quad layouts** where the target edge length is comparable to the distance between singularities. Drag the triangle vertices to see how the Jacobian determinant changes. When the triangle flips (vertices go clockwise), the determinant becomes negative: ```{ojs} //| label: fig-jacobian-degeneracy //| code-fold: true { const w = 500, h = 420; const svg = d3.create("svg").attr("viewBox", `0 0 ${w} ${h}`).style("width", "100%").style("max-width", w + "px"); // Mutable vertex state stored on the SVG node itself so it persists across redraws const node = svg.node(); if (!node.__verts) { node.__verts = [[150, 320], [350, 320], [250, 120]]; } const verts = node.__verts; function computeMathDet(v) { const ux = v[1][0] - v[0][0], uy = v[1][1] - v[0][1]; const vx = v[2][0] - v[0][0], vy = v[2][1] - v[0][1]; // Negate because SVG y-axis is flipped vs. math convention return -(ux * vy - uy * vx); } function redraw() { svg.selectAll("*").remove(); // Integer grid for (let x = 50; x <= 450; x += 50) svg.append("line").attr("x1", x).attr("y1", 20).attr("x2", x).attr("y2", 400) .attr("stroke", "#eee").attr("stroke-width", 0.5); for (let y = 20; y <= 400; y += 50) svg.append("line").attr("x1", 50).attr("y1", y).attr("x2", 450).attr("y2", y) .attr("stroke", "#eee").attr("stroke-width", 0.5); const mathDet = computeMathDet(verts); // Triangle fill color const fillColor = mathDet > 500 ? "rgba(39, 174, 96, 0.25)" : mathDet > 100 ? "rgba(39, 174, 96, 0.15)" : mathDet > 0 ? "rgba(243, 156, 18, 0.25)" : mathDet === 0 ? "rgba(149, 165, 166, 0.3)" : "rgba(231, 76, 60, 0.25)"; const strokeColor = mathDet > 0 ? "#27ae60" : mathDet === 0 ? "#95a5a6" : "#e74c3c"; svg.append("polygon") .attr("points", verts.map(v => v.join(",")).join(" ")) .attr("fill", fillColor).attr("stroke", strokeColor).attr("stroke-width", 2.5); // Dashed edge vectors from v0 svg.append("line").attr("x1", verts[0][0]).attr("y1", verts[0][1]) .attr("x2", verts[1][0]).attr("y2", verts[1][1]) .attr("stroke", "steelblue").attr("stroke-width", 1.5).attr("stroke-dasharray", "4,2"); svg.append("line").attr("x1", verts[0][0]).attr("y1", verts[0][1]) .attr("x2", verts[2][0]).attr("y2", verts[2][1]) .attr("stroke", "#e67e22").attr("stroke-width", 1.5).attr("stroke-dasharray", "4,2"); // Draggable vertices const labels = ["u\u2080", "u\u2081", "u\u2082"]; const colors = ["#2c3e50", "steelblue", "#e67e22"]; for (let i = 0; i < 3; i++) { svg.append("circle").attr("cx", verts[i][0]).attr("cy", verts[i][1]).attr("r", 10) .attr("fill", colors[i]).attr("cursor", "grab").attr("data-idx", i) .call(d3.drag() .on("start", function() { d3.select(this).attr("cursor", "grabbing"); }) .on("drag", function(event) { const idx = +d3.select(this).attr("data-idx"); verts[idx][0] = Math.max(20, Math.min(w - 20, event.x)); verts[idx][1] = Math.max(20, Math.min(h - 30, event.y)); redraw(); }) .on("end", function() { d3.select(this).attr("cursor", "grab"); }) ); svg.append("text").attr("x", verts[i][0]).attr("y", verts[i][1] - 15) .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", colors[i]) .attr("font-weight", "bold").attr("pointer-events", "none") .text(labels[i]); } // Determinant display const statusText = mathDet > 500 ? "Valid (well-shaped)" : mathDet > 100 ? "Valid" : mathDet > 0 ? "Near-degenerate" : mathDet === 0 ? "DEGENERATE (collapsed)" : "FLIPPED (negative orientation)"; const statusColor = mathDet > 0 ? "#27ae60" : mathDet === 0 ? "#95a5a6" : "#e74c3c"; svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle") .attr("font-size", 15).attr("fill", statusColor).attr("font-weight", "bold") .text(`det = ${mathDet.toFixed(0)} \u2014 ${statusText}`); } redraw(); return node; } ``` ::: {.callout-tip} ## Why the Jacobian matters For a piecewise linear map $f$, each triangle has a constant $2 \times 2$ Jacobian $J_t = \nabla f_t$. The determinant $\det(J_t)$ measures the signed area ratio between the mapped triangle and the original. When $\det(J_t) > 0$, the triangle preserves orientation. When $\det(J_t) \leq 0$, the triangle has collapsed or flipped — the integer iso-lines cross or overlap, making quad extraction impossible. ::: ### 1.2 Contributions The paper makes three main contributions: 1. **A MIQP formulation** where the minimizer is guaranteed to be an Integer-Grid Map (Sec. 3.1) 2. **A complexity reduction algorithm** that accurately approximates the high-dimensional problem with a low-dimensional one (Sec. 3.2) 3. **Singularity separating conditions** that significantly improve the continuous relaxation of the MIQP (Sec. 3.3) --- ## 2. Integer-Grid Maps > *"The main principle of parametrization based quad meshing algorithms is the mapping of the canonical quad mesh formed by the 2D Cartesian grid of integer iso-lines onto a surface embedded in 3D."* ### 2.1 Piecewise Linear Maps Given a triangle mesh $\mathcal{M} = (V, E, T)$, a map $f$ is the union of per-triangle affine maps: $$f_i : (u_i, v_i, w_i) \in \mathbb{R}^{2 \times 3} \mapsto (p_i, q_i, r_i) \in \mathbb{R}^{3 \times 3}$$ Each triangle $t_i$ in 3D is mapped to a triangle in the 2D parametric domain. The intersection of the integer grid $\mathbb{Z} \times \mathbb{Z}$ with the image of $f$ defines the quad mesh edges. ### 2.2 Transition Functions Since each triangle has its own chart, we need **transition functions** to describe how neighboring charts relate. The transition function from the chart of triangle $t_i$ to that of neighboring triangle $t_j$ is: $$g_{i \to j}(\mathbf{a}) = R_{90}^{r_{ij}} \mathbf{a} + \mathbf{t}_{ij} \tag{1}$$ where $r_{ij} \in \{0, 1, 2, 3\}$ is a rotation by multiples of $\pi/2$ and $\mathbf{t}_{ij} \in \mathbb{Z}^2$ is an integer translation. This means neighboring charts can differ by a **90-degree rotation** and an **integer shift** — both of which leave the Cartesian integer grid invariant. This is exactly the definition of an **integer-grid automorphism**. Adjust the rotation $r_{ij}$ and translation $\mathbf{t}_{ij}$ to see how the transition function relates two neighboring charts: ```{ojs} //| label: fig-transition-function viewof r_ij = Inputs.range([0, 3], {value: 0, step: 1, label: "Rotation r_ij (× 90°)"}) ``` ```{ojs} //| label: fig-transition-t viewof t_ij_x = Inputs.range([-2, 2], {value: 0, step: 1, label: "Translation t_x"}) ``` ```{ojs} //| label: fig-transition-ty viewof t_ij_y = Inputs.range([-2, 2], {value: 0, step: 1, label: "Translation t_y"}) ``` ```{ojs} //| code-fold: true { const w = 700, h = 380; const svg = d3.create("svg").attr("viewBox", `0 0 ${w} ${h}`).style("width", "100%").style("max-width", w + "px"); const leftCx = 170, rightCx = 530, cy = 200; const scale = 50; // Draw grids function drawGrid(gx, gy, label) { for (let i = -3; i <= 3; i++) { svg.append("line") .attr("x1", gx + i*scale).attr("y1", gy - 3*scale) .attr("x2", gx + i*scale).attr("y2", gy + 3*scale) .attr("stroke", "#eee").attr("stroke-width", i === 0 ? 1 : 0.5); svg.append("line") .attr("x1", gx - 3*scale).attr("y1", gy + i*scale) .attr("x2", gx + 3*scale).attr("y2", gy + i*scale) .attr("stroke", "#eee").attr("stroke-width", i === 0 ? 1 : 0.5); } // Integer points for (let ix = -2; ix <= 2; ix++) { for (let iy = -2; iy <= 2; iy++) { svg.append("circle") .attr("cx", gx + ix*scale).attr("cy", gy - iy*scale).attr("r", 2) .attr("fill", "#ccc"); } } svg.append("text").attr("x", gx).attr("y", gy + 3.5*scale) .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", "#333").attr("font-weight", "bold") .text(label); } drawGrid(leftCx, cy, "Chart of t_i"); drawGrid(rightCx, cy, "Chart of t_j"); // Triangle in chart i const triI = [[0.3, 0.5], [1.8, 0.2], [1.0, 1.5]]; const triIPts = triI.map(p => [leftCx + p[0]*scale, cy - p[1]*scale]); svg.append("polygon") .attr("points", triIPts.map(p => p.join(",")).join(" ")) .attr("fill", "rgba(70,130,180,0.15)").attr("stroke", "steelblue").attr("stroke-width", 2); // Shared edge (highlighted) svg.append("line") .attr("x1", triIPts[1][0]).attr("y1", triIPts[1][1]) .attr("x2", triIPts[2][0]).attr("y2", triIPts[2][1]) .attr("stroke", "#e74c3c").attr("stroke-width", 3); // Apply transition function to triangle const angle = r_ij * Math.PI / 2; const cosA = Math.round(Math.cos(angle)), sinA = Math.round(Math.sin(angle)); function applyTransition(p) { const rx = cosA * p[0] - sinA * p[1] + t_ij_x; const ry = sinA * p[0] + cosA * p[1] + t_ij_y; return [rx, ry]; } const triJ = triI.map(p => applyTransition(p)); const triJPts = triJ.map(p => [rightCx + p[0]*scale, cy - p[1]*scale]); svg.append("polygon") .attr("points", triJPts.map(p => p.join(",")).join(" ")) .attr("fill", "rgba(230,126,34,0.15)").attr("stroke", "#e67e22").attr("stroke-width", 2); // Shared edge in j's chart svg.append("line") .attr("x1", triJPts[1][0]).attr("y1", triJPts[1][1]) .attr("x2", triJPts[2][0]).attr("y2", triJPts[2][1]) .attr("stroke", "#e74c3c").attr("stroke-width", 3); // Arrow between charts svg.append("line") .attr("x1", leftCx + 3.2*scale).attr("y1", cy) .attr("x2", rightCx - 3.2*scale).attr("y2", cy) .attr("stroke", "#999").attr("stroke-width", 1.5).attr("marker-end", "url(#arrowhead)"); svg.append("defs").append("marker").attr("id", "arrowhead").attr("viewBox", "0 0 10 10") .attr("refX", 10).attr("refY", 5).attr("markerWidth", 8).attr("markerHeight", 8) .attr("orient", "auto") .append("path").attr("d", "M 0 0 L 10 5 L 0 10 z").attr("fill", "#999"); svg.append("text") .attr("x", (leftCx + rightCx) / 2).attr("y", cy - 15) .attr("text-anchor", "middle").attr("font-size", 13).attr("fill", "#666") .text(`g: R${r_ij * 90}° + (${t_ij_x}, ${t_ij_y})`); // Vertex labels for (let i = 0; i < 3; i++) { svg.append("circle").attr("cx", triIPts[i][0]).attr("cy", triIPts[i][1]).attr("r", 4).attr("fill", "steelblue"); svg.append("circle").attr("cx", triJPts[i][0]).attr("cy", triJPts[i][1]).attr("r", 4).attr("fill", "#e67e22"); } return svg.node(); } ``` ::: {.callout-tip} ## Why integer translations? The integer grid $\mathbb{Z} \times \mathbb{Z}$ is invariant under 90-degree rotations and integer translations. So if chart $i$ and chart $j$ are related by such a transformation, the integer iso-lines in one chart **perfectly align** with those in the other. This is what makes the quad mesh well-defined across chart boundaries. ::: ### 2.3 The Three IGM Conditions The class of **Integer-Grid Maps** is the subset of all possible maps that additionally stitch the grid of integer iso-lines into a valid quad mesh. The necessary and sufficient conditions are: **Condition 1 — Transition Functions**: The transition function $g_{i \to j}$ from the chart of triangle $t_i$ to the chart of neighboring triangle $t_j$ must be an integer-grid automorphism of the form $g_{i \to j}(\mathbf{a}) = R_{90}^{r_{ij}} \mathbf{a} + \mathbf{t}_{ij}$ (Eq. 1). **Condition 2 — Singular Points**: All singular vertices must lie on integer locations in the domain: $$\mathbf{f}^{-1}(s_i) \in \mathbb{Z}^2 \quad \forall s_i \in S \tag{2}$$ This ensures the cone singularities (where the cross field has nonzero index) map to grid intersections, which become the irregular vertices of the quad mesh. **Condition 3 — Consistent Orientation**: All domain triangles $(\mathbf{u}, \mathbf{v}, \mathbf{w})$ with $\mathbf{u}, \mathbf{v}, \mathbf{w} \in \mathbb{R}^2$ must have positive orientation: $$\det[\mathbf{v} - \mathbf{u},\; \mathbf{w} - \mathbf{u}] > 0 \tag{3}$$ This prevents triangle flipping — the degeneracy that makes greedy rounding fail. > *"For the remainder of this paper it is very important to keep in mind that an IGM is a bijective map between two piecewise linear 2-manifold meshes."* Toggle between the three conditions to see what each one enforces on a small mesh in parametric space: ```{ojs} //| label: fig-igm-conditions viewof igm_cond = Inputs.radio(["All conditions", "1: Transition functions", "2: Singular vertices on integers", "3: Consistent orientation"], {value: "All conditions", label: "Show condition"}) ``` ```{ojs} //| code-fold: true { const w = 550, h = 450; const svg = d3.create("svg").attr("viewBox", `0 0 ${w} ${h}`).style("width", "100%").style("max-width", w + "px"); const ox = 75, oy = 375, sc = 70; // Integer grid for (let i = 0; i <= 5; i++) { svg.append("line").attr("x1", ox + i*sc).attr("y1", oy - 5*sc).attr("x2", ox + i*sc).attr("y2", oy) .attr("stroke", "#e8e8e8").attr("stroke-width", 1); svg.append("line").attr("x1", ox).attr("y1", oy - i*sc).attr("x2", ox + 5*sc).attr("y2", oy - i*sc) .attr("stroke", "#e8e8e8").attr("stroke-width", 1); // Integer labels svg.append("text").attr("x", ox + i*sc).attr("y", oy + 18) .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#aaa").text(i); if (i > 0) svg.append("text").attr("x", ox - 15).attr("y", oy - i*sc + 4) .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#aaa").text(i); } // Integer points for (let ix = 0; ix <= 5; ix++) { for (let iy = 0; iy <= 5; iy++) { svg.append("circle") .attr("cx", ox + ix*sc).attr("cy", oy - iy*sc).attr("r", 2).attr("fill", "#ddd"); } } // Define a small mesh: 6 triangles forming a patch with 1 singularity // Vertices in parametric space (u, v) const verts = [ [0.5, 0.5], // 0 [2.0, 0.3], // 1 [3.5, 0.8], // 2 [0.3, 2.0], // 3 [2.0, 2.0], // 4: singularity (on integer!) [3.8, 2.2], // 5 [0.8, 3.5], // 6 [2.2, 3.7], // 7 [3.5, 3.3], // 8 ]; const tris = [ [0, 1, 4], [1, 2, 5], [1, 5, 4], [0, 4, 3], [3, 4, 6], [4, 7, 6], [4, 5, 8], [4, 8, 7], ]; const edges = [ [0,1],[1,2],[2,5],[5,8],[8,7],[7,6],[6,3],[3,0], [0,4],[1,4],[1,5],[3,4],[4,5],[4,6],[4,7],[4,8] ]; const singIdx = 4; // vertex 4 is the singularity const showCond1 = igm_cond === "All conditions" || igm_cond === "1: Transition functions"; const showCond2 = igm_cond === "All conditions" || igm_cond === "2: Singular vertices on integers"; const showCond3 = igm_cond === "All conditions" || igm_cond === "3: Consistent orientation"; // Draw triangles for (const tri of tris) { const pts = tri.map(i => [ox + verts[i][0]*sc, oy - verts[i][1]*sc]); // Condition 3: color by orientation let fillC = "rgba(200,200,200,0.08)"; if (showCond3) { const ux = verts[tri[1]][0] - verts[tri[0]][0]; const uy = verts[tri[1]][1] - verts[tri[0]][1]; const vx = verts[tri[2]][0] - verts[tri[0]][0]; const vy = verts[tri[2]][1] - verts[tri[0]][1]; const det = ux * vy - uy * vx; fillC = det > 0 ? "rgba(39, 174, 96, 0.12)" : "rgba(231, 76, 60, 0.2)"; } svg.append("polygon") .attr("points", pts.map(p => p.join(",")).join(" ")) .attr("fill", fillC) .attr("stroke", "#bbb").attr("stroke-width", 1); } // Condition 1: highlight internal edges showing transition compatibility if (showCond1) { const internalEdges = [[0,4],[1,4],[1,5],[3,4],[4,5],[4,6],[4,7],[4,8]]; for (const e of internalEdges) { const p0 = [ox + verts[e[0]][0]*sc, oy - verts[e[0]][1]*sc]; const p1 = [ox + verts[e[1]][0]*sc, oy - verts[e[1]][1]*sc]; svg.append("line") .attr("x1", p0[0]).attr("y1", p0[1]).attr("x2", p1[0]).attr("y2", p1[1]) .attr("stroke", "steelblue").attr("stroke-width", 2.5).attr("opacity", 0.7); } svg.append("text").attr("x", w - 10).attr("y", 25) .attr("text-anchor", "end").attr("font-size", 12).attr("fill", "steelblue") .text("Internal edges: transition functions must be IGM automorphisms"); } // Condition 2: highlight singularity if (showCond2) { const sx = ox + verts[singIdx][0]*sc, sy = oy - verts[singIdx][1]*sc; // Check if on integer const onInt = Math.abs(verts[singIdx][0] - Math.round(verts[singIdx][0])) < 0.01 && Math.abs(verts[singIdx][1] - Math.round(verts[singIdx][1])) < 0.01; svg.append("circle").attr("cx", sx).attr("cy", sy).attr("r", 14) .attr("fill", "none").attr("stroke", onInt ? "#27ae60" : "#e74c3c") .attr("stroke-width", 2.5).attr("stroke-dasharray", "4,2"); svg.append("text").attr("x", sx + 18).attr("y", sy - 8) .attr("font-size", 11).attr("fill", onInt ? "#27ae60" : "#e74c3c") .text(onInt ? "on Z^2" : "NOT on Z^2"); svg.append("text").attr("x", w - 10).attr("y", 45) .attr("text-anchor", "end").attr("font-size", 12).attr("fill", "#e74c3c") .text("Singular vertex must map to integer coordinates"); } // Condition 3: label if (showCond3) { svg.append("text").attr("x", w - 10).attr("y", 65) .attr("text-anchor", "end").attr("font-size", 12).attr("fill", "#27ae60") .text("All triangles must have positive orientation (green)"); } // Draw vertices for (let i = 0; i < verts.length; i++) { const px = ox + verts[i][0]*sc, py = oy - verts[i][1]*sc; const isSing = i === singIdx; svg.append("circle").attr("cx", px).attr("cy", py) .attr("r", isSing ? 6 : 3.5) .attr("fill", isSing ? "#e74c3c" : "#2c3e50"); } svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", "#666") .text("A small mesh patch in parametric space with integer grid overlay"); return svg.node(); } ``` --- ## 3. The Quality Metric ### 3.1 The Energy $E(f)$ In parametrization-based quad remeshing, a variational quality metric $E_q(f)$ penalizes distortion of the resulting quad elements based on the map $f$. The standard form is a Frobenius-norm measure of how far the Jacobian is from the desired sizing and orientation: $$E(f) = \frac{1}{2} \sum_{t \in T} \mathcal{A}_t \| R_t^T \nabla f_t - H_t \|_{\text{Frobenius}}^2 \tag{E1}$$ where: - $H_t = \text{diag}(w_t^{-1}, h_t^{-1})$ is the locally desired anisotropic sizing - $R_t \in \mathbb{R}^{4 \times 2}$ contains the normalized $u$ and $v$ directions of the cross field - $\mathcal{A}_t$ is the triangle area This energy simultaneously measures how well the parameterization aligns with the cross-field directions and reproduces the specified element density. ### 3.2 Generalized Anisotropy Metric To weight the importance of cross direction versus sizing fidelity, the authors generalize with a higher-order norm: $$E_\alpha^k(f) = \frac{1}{2} \sum_{t \in T} \mathcal{A}_t \| R_t^T \nabla f_t - H_t \|_\alpha^k \tag{E2}$$ where $k \in 2\mathbb{Z}^+$, $\alpha \in [0,1]$, and the anisotropy norm is: $$\left\| \begin{pmatrix} a & b \\ c & d \end{pmatrix} \right\|_\alpha^k = 2\alpha(a^k + d^k) + 2(1 - \alpha)(b^k + c^k)$$ For $\alpha < 0.5$, stretching along the cross-field directions is preferred over angular deviation — useful for coarse quad layouts with cross-field-aligned rectangular patches. ### 3.3 The Naive MINLP (Problem P1) Naively searching for the best IGM leads to a **Mixed-Integer Nonlinear Program**: $$\text{minimize } E_q(f) \quad \text{s.t.} \quad (1),\; (2),\; (3) \tag{P1}$$ This has $6|T| + 3|E|$ unknowns including at least $3|E|$ discrete variables. Due to: - (1) containing **nonlinear** constraints in $r_{ij}$ and **linear** in $\mathbf{t}_{ij}$ - (3) imposing $|T|$ **non-convex quadratic inequality** constraints - The resulting number of unknowns being $O(|V|)$ > *"Unfortunately MINLP problems are very hard to solve since they imply all difficulties from continuous as well as discrete optimization."* Modern mixed-integer solvers need three properties: **(i)** convex objective and constraints, **(ii)** low-dimensional search space, **(iii)** continuous relaxation close to the optimum. The next three sections address each. --- ## 4. Convex Consistent Orientation > *"The feasible set of the consistent orientation condition (3) is a non-convex region $\mathcal{N} \subset \mathbb{R}^6$."* ### 4.1 The Non-Convex Problem Condition (3) requires $\det[\mathbf{v}-\mathbf{u}, \mathbf{w}-\mathbf{u}] > 0$ for every triangle $(\mathbf{u}, \mathbf{v}, \mathbf{w})$. This determinant is a **bilinear** (quadratic) function of the vertex positions — the set of positions satisfying it is non-convex. We need a **convex inner approximation** $\mathcal{C} \subset \mathcal{N}$ such that every point in $\mathcal{C}$ satisfies (3). ### 4.2 The Trisector Construction The key geometric insight: partition $\mathbb{R}^2$ into 3 sectors using ccw-ordered rays $\mathbf{s}_0, \mathbf{s}_1, \mathbf{s}_2$ emanating from a center point $\mathbf{m}$: $$S_i = \{ \mathbf{u} \in \mathbb{R}^2 : (\mathbf{u} - \mathbf{m}) \cdot \mathbf{s}_i^\perp > 0 \;\wedge\; (\mathbf{u} - \mathbf{m}) \cdot \mathbf{s}_{i+1}^\perp < 0 \}$$ If vertex $\mathbf{u}_i$ stays strictly within sector $S_i$ and $\mathbf{m}$ is an interior point of the triangle, then the triangle **cannot flip**. The area is guaranteed positive because it decomposes into three sub-triangles $T_i = (\mathbf{u}_i, \mathbf{u}_{i+1}, \mathbf{m})$, each with positive area. ### 4.3 The Fermat Point > *"To enable maximal rotational freedom we choose a trisector origin $\mathbf{m}$ such that angles between $(\hat{\mathbf{u}}_i - \mathbf{m})$ and $(\hat{\mathbf{u}}_{i+1} - \mathbf{m})$ are $120°$ each."* The **Fermat point** (first Torricelli point) is the point inside a triangle where each pair of vertices subtends an angle of exactly $120°$. By placing the trisector origin at the Fermat point, the three sectors are equally sized ($120°$ each), maximizing the space of feasible deformations. The Fermat point is constructed by attaching equilateral triangles to each edge and intersecting the three connecting diagonals. Drag the triangle vertices to explore the Fermat point and trisector sectors. As long as each vertex stays in its colored sector, the triangle orientation is guaranteed positive: ```{ojs} //| label: fig-trisector viewof tri_preset = Inputs.radio(["Equilateral", "Right triangle", "Obtuse (> 100°)"], {value: "Equilateral", label: "Triangle preset"}) ``` ```{ojs} //| code-fold: true { const totalW = 900, panelH = 420, gap = 20; const leftW = 420, rightW = totalW - leftW - gap; const svg = d3.create("svg") .attr("viewBox", `0 0 ${totalW} ${panelH + 30}`) .style("width", "100%").style("max-width", totalW + "px"); let triV; if (tri_preset === "Equilateral") { triV = [[130, 360], [330, 360], [230, 186]]; } else if (tri_preset === "Right triangle") { triV = [[100, 360], [360, 360], [100, 130]]; } else { triV = [[80, 320], [380, 320], [270, 230]]; } // Shared Fermat point computation function fermatPoint(A, B, C) { function angle(p1, p2, p3) { const v1x = p1[0]-p2[0], v1y = p1[1]-p2[1]; const v2x = p3[0]-p2[0], v2y = p3[1]-p2[1]; const dot = v1x*v2x + v1y*v2y; const m1 = Math.sqrt(v1x*v1x + v1y*v1y); const m2 = Math.sqrt(v2x*v2x + v2y*v2y); return Math.acos(Math.max(-1, Math.min(1, dot/(m1*m2)))); } const angA = angle(B, A, C), angB = angle(A, B, C), angC = angle(A, C, B); if (angA >= 2*Math.PI/3) return [...A]; if (angB >= 2*Math.PI/3) return [...B]; if (angC >= 2*Math.PI/3) return [...C]; let mx = (A[0]+B[0]+C[0])/3, my = (A[1]+B[1]+C[1])/3; for (let iter = 0; iter < 200; iter++) { const dA = Math.sqrt((mx-A[0])**2 + (my-A[1])**2) || 1e-10; const dB = Math.sqrt((mx-B[0])**2 + (my-B[1])**2) || 1e-10; const dC = Math.sqrt((mx-C[0])**2 + (my-C[1])**2) || 1e-10; const wA = 1/dA, wB = 1/dB, wC = 1/dC; const wSum = wA + wB + wC; mx = (wA*A[0] + wB*B[0] + wC*C[0]) / wSum; my = (wA*A[1] + wB*B[1] + wC*C[1]) / wSum; } return [mx, my]; } const fp = fermatPoint(triV[0], triV[1], triV[2]); const sectorColors = ["rgba(231,76,60,0.1)", "rgba(39,174,96,0.1)", "rgba(52,152,219,0.1)"]; const sectorStroke = ["#e74c3c", "#27ae60", "#3498db"]; const vLabels = ["u\u2080", "u\u2081", "u\u2082"]; // --- Helper: draw trisector scene into a clipping group --- function drawScene(g, clipX, clipY, clipW, clipH, showZoomBox) { // Clip to panel bounds const clipId = `clip-${clipX}`; const defs = g.append("defs"); defs.append("clipPath").attr("id", clipId) .append("rect").attr("x", clipX).attr("y", clipY).attr("width", clipW).attr("height", clipH); const scene = g.append("g").attr("clip-path", `url(#${clipId})`); // Panel background scene.append("rect").attr("x", clipX).attr("y", clipY) .attr("width", clipW).attr("height", clipH) .attr("fill", "#fafafa").attr("stroke", "#ccc").attr("stroke-width", 1); // Sector rays and fills const rayLen = 800; for (let i = 0; i < 3; i++) { const dx = triV[i][0] - fp[0], dy = triV[i][1] - fp[1]; const len = Math.sqrt(dx*dx + dy*dy) || 1; const endX = fp[0] + dx/len * rayLen; const endY = fp[1] + dy/len * rayLen; scene.append("line") .attr("x1", fp[0]).attr("y1", fp[1]).attr("x2", endX).attr("y2", endY) .attr("stroke", sectorStroke[i]).attr("stroke-width", 1.5) .attr("stroke-dasharray", "6,3").attr("opacity", 0.6); } for (let i = 0; i < 3; i++) { const j = (i + 1) % 3; const dx1 = triV[i][0] - fp[0], dy1 = triV[i][1] - fp[1]; const dx2 = triV[j][0] - fp[0], dy2 = triV[j][1] - fp[1]; const len1 = Math.sqrt(dx1*dx1+dy1*dy1) || 1; const len2 = Math.sqrt(dx2*dx2+dy2*dy2) || 1; const e1x = fp[0] + dx1/len1*rayLen, e1y = fp[1] + dy1/len1*rayLen; const e2x = fp[0] + dx2/len2*rayLen, e2y = fp[1] + dy2/len2*rayLen; scene.append("polygon") .attr("points", `${fp[0]},${fp[1]} ${e1x},${e1y} ${e2x},${e2y}`) .attr("fill", sectorColors[i]).attr("opacity", 0.5); } // Triangle scene.append("polygon") .attr("points", triV.map(v => v.join(",")).join(" ")) .attr("fill", "rgba(255,255,255,0.6)").attr("stroke", "#2c3e50").attr("stroke-width", 2.5); // Fermat point scene.append("circle").attr("cx", fp[0]).attr("cy", fp[1]).attr("r", 6).attr("fill", "#8e44ad"); scene.append("text").attr("x", fp[0] + 12).attr("y", fp[1] - 8) .attr("font-size", 14).attr("fill", "#8e44ad").attr("font-weight", "bold").text("m"); // Vertices for (let i = 0; i < 3; i++) { scene.append("circle").attr("cx", triV[i][0]).attr("cy", triV[i][1]).attr("r", 7) .attr("fill", sectorStroke[i]); scene.append("text").attr("x", triV[i][0]).attr("y", triV[i][1] - 14) .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", sectorStroke[i]).attr("font-weight", "bold") .text(vLabels[i]); } // Angle labels at Fermat point for (let i = 0; i < 3; i++) { const j = (i + 1) % 3; const a1 = Math.atan2(triV[i][1] - fp[1], triV[i][0] - fp[0]); const a2 = Math.atan2(triV[j][1] - fp[1], triV[j][0] - fp[0]); let da = a2 - a1; while (da < -Math.PI) da += 2*Math.PI; while (da > Math.PI) da -= 2*Math.PI; const angleDeg = Math.abs(da) * 180 / Math.PI; const arcR = 30; const mid = a1 + da/2; scene.append("text") .attr("x", fp[0] + (arcR + 15) * Math.cos(mid)) .attr("y", fp[1] + (arcR + 15) * Math.sin(mid)) .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#8e44ad") .text(`${angleDeg.toFixed(0)}\u00B0`); } // Zoom box indicator on the overview panel if (showZoomBox) { const zoomR = 90; scene.append("rect") .attr("x", fp[0] - zoomR).attr("y", fp[1] - zoomR) .attr("width", 2*zoomR).attr("height", 2*zoomR) .attr("fill", "none").attr("stroke", "#8e44ad") .attr("stroke-width", 2).attr("stroke-dasharray", "4,3"); } } // --- Left panel: zoomed-out overview --- const leftG = svg.append("g"); drawScene(leftG, 0, 25, leftW, panelH, true); svg.append("text").attr("x", leftW/2).attr("y", 16).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", "#555").text("Full triangle with trisector"); // --- Right panel: zoomed-in around Fermat point --- const rx = leftW + gap; const zoomR = 90; // radius in scene coords to show around Fermat point const zoomScale = Math.min(rightW, panelH) / (2 * zoomR); const rightG = svg.append("g"); // Background rightG.append("rect").attr("x", rx).attr("y", 25).attr("width", rightW).attr("height", panelH) .attr("fill", "#fafafa").attr("stroke", "#ccc").attr("stroke-width", 1); // Clip for right panel const clipR = "clip-right"; rightG.append("defs").append("clipPath").attr("id", clipR) .append("rect").attr("x", rx).attr("y", 25).attr("width", rightW).attr("height", panelH); // Transform: center Fermat point in the right panel, scale up const rcx = rx + rightW/2, rcy = 25 + panelH/2; const zoomedG = rightG.append("g").attr("clip-path", `url(#${clipR})`) .append("g") .attr("transform", `translate(${rcx}, ${rcy}) scale(${zoomScale}) translate(${-fp[0]}, ${-fp[1]})`); // Redraw the scene geometry into the zoomed group (without clipping, the outer clip handles it) const zRayLen = 800; for (let i = 0; i < 3; i++) { const dx = triV[i][0] - fp[0], dy = triV[i][1] - fp[1]; const len = Math.sqrt(dx*dx + dy*dy) || 1; const endX = fp[0] + dx/len * zRayLen; const endY = fp[1] + dy/len * zRayLen; zoomedG.append("line") .attr("x1", fp[0]).attr("y1", fp[1]).attr("x2", endX).attr("y2", endY) .attr("stroke", sectorStroke[i]).attr("stroke-width", 1.5 / zoomScale) .attr("stroke-dasharray", `${6/zoomScale},${3/zoomScale}`).attr("opacity", 0.6); } for (let i = 0; i < 3; i++) { const j = (i + 1) % 3; const dx1 = triV[i][0] - fp[0], dy1 = triV[i][1] - fp[1]; const dx2 = triV[j][0] - fp[0], dy2 = triV[j][1] - fp[1]; const len1 = Math.sqrt(dx1*dx1+dy1*dy1) || 1; const len2 = Math.sqrt(dx2*dx2+dy2*dy2) || 1; const e1x = fp[0] + dx1/len1*zRayLen, e1y = fp[1] + dy1/len1*zRayLen; const e2x = fp[0] + dx2/len2*zRayLen, e2y = fp[1] + dy2/len2*zRayLen; zoomedG.append("polygon") .attr("points", `${fp[0]},${fp[1]} ${e1x},${e1y} ${e2x},${e2y}`) .attr("fill", sectorColors[i]).attr("opacity", 0.5); } // Triangle in zoomed view zoomedG.append("polygon") .attr("points", triV.map(v => v.join(",")).join(" ")) .attr("fill", "rgba(255,255,255,0.6)").attr("stroke", "#2c3e50") .attr("stroke-width", 2.5 / zoomScale); // Fermat point zoomedG.append("circle").attr("cx", fp[0]).attr("cy", fp[1]) .attr("r", 6 / zoomScale).attr("fill", "#8e44ad"); zoomedG.append("text").attr("x", fp[0] + 12 / zoomScale).attr("y", fp[1] - 8 / zoomScale) .attr("font-size", 14 / zoomScale).attr("fill", "#8e44ad").attr("font-weight", "bold").text("m"); // Vertices in zoomed view for (let i = 0; i < 3; i++) { zoomedG.append("circle").attr("cx", triV[i][0]).attr("cy", triV[i][1]) .attr("r", 7 / zoomScale).attr("fill", sectorStroke[i]); zoomedG.append("text").attr("x", triV[i][0]).attr("y", triV[i][1] - 14 / zoomScale) .attr("text-anchor", "middle").attr("font-size", 14 / zoomScale) .attr("fill", sectorStroke[i]).attr("font-weight", "bold") .text(vLabels[i]); } // Angle labels in zoomed view for (let i = 0; i < 3; i++) { const j = (i + 1) % 3; const a1 = Math.atan2(triV[i][1] - fp[1], triV[i][0] - fp[0]); const a2 = Math.atan2(triV[j][1] - fp[1], triV[j][0] - fp[0]); let da = a2 - a1; while (da < -Math.PI) da += 2*Math.PI; while (da > Math.PI) da -= 2*Math.PI; const angleDeg = Math.abs(da) * 180 / Math.PI; const arcR = 30 / zoomScale; const mid = a1 + da/2; zoomedG.append("text") .attr("x", fp[0] + (arcR + 15 / zoomScale) * Math.cos(mid)) .attr("y", fp[1] + (arcR + 15 / zoomScale) * Math.sin(mid)) .attr("text-anchor", "middle").attr("font-size", 10 / zoomScale).attr("fill", "#8e44ad") .text(`${angleDeg.toFixed(0)}\u00B0`); } svg.append("text").attr("x", rx + rightW/2).attr("y", 16).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", "#555").text("Zoomed in: Fermat point detail"); return svg.node(); } ``` ::: {.callout-note} ## The power of linearization The trisector construction replaces one **non-convex quadratic** constraint ($\det > 0$) with **six linear inequalities** (each vertex must stay in its sector). This is conservative — some valid configurations are excluded — but it guarantees that every solution within the convex feasible region is a valid, non-degenerate triangle. ::: ### 4.4 Virtual Splitting for Obtuse Triangles The Fermat point only exists inside the triangle when **all angles are less than $120°$**. For triangles with an angle larger than $100°$ (a conservative threshold), the paper "virtually" splits the triangle along the altitude at its largest angle: ```{ojs} //| label: fig-virtual-split viewof obtuse_angle = Inputs.range([60, 150], {value: 90, step: 1, label: "Largest angle at c (degrees)"}) ``` ```{ojs} //| code-fold: true { const w = 700, h = 420; const svg = d3.create("svg").attr("viewBox", `0 0 ${w} ${h}`).style("width", "100%").style("max-width", w + "px"); // --- Construct triangle with exact angle at C --- const baseLen = 180; const baseY = 340; const angleRad = obtuse_angle * Math.PI / 180; // Place A and B symmetrically. C is above with the specified angle ACB. // Using the inscribed angle: height from midpoint of AB to C is // h = (baseLen/2) / tan(angle/2) const halfBase = baseLen / 2; const triH = halfBase / Math.tan(angleRad / 2); // Left panel: original triangle with Fermat point and sectors const leftW = w * 0.48; const leftCx = leftW / 2 + 15; const A = [leftCx - halfBase, baseY]; const B = [leftCx + halfBase, baseY]; const C = [leftCx, baseY - Math.max(25, triH)]; const needsSplit = obtuse_angle > 100; // Compute actual angles at all three vertices function vecAngle(p1, vertex, p2) { const v1x = p1[0]-vertex[0], v1y = p1[1]-vertex[1]; const v2x = p2[0]-vertex[0], v2y = p2[1]-vertex[1]; const dot = v1x*v2x + v1y*v2y; const m1 = Math.sqrt(v1x*v1x+v1y*v1y), m2 = Math.sqrt(v2x*v2x+v2y*v2y); return Math.acos(Math.max(-1, Math.min(1, dot/(m1*m2)))); } const angC = vecAngle(A, C, B) * 180 / Math.PI; const angA = vecAngle(B, A, C) * 180 / Math.PI; const angB = vecAngle(A, B, C) * 180 / Math.PI; // Fermat point via Weiszfeld function fermatPt(P, Q, R) { const aP = vecAngle(Q, P, R), aQ = vecAngle(P, Q, R), aR = vecAngle(P, R, Q); if (aP >= 2*Math.PI/3) return [...P]; if (aQ >= 2*Math.PI/3) return [...Q]; if (aR >= 2*Math.PI/3) return [...R]; let mx = (P[0]+Q[0]+R[0])/3, my = (P[1]+Q[1]+R[1])/3; for (let it = 0; it < 300; it++) { const dP = Math.sqrt((mx-P[0])**2+(my-P[1])**2) || 1e-10; const dQ = Math.sqrt((mx-Q[0])**2+(my-Q[1])**2) || 1e-10; const dR = Math.sqrt((mx-R[0])**2+(my-R[1])**2) || 1e-10; const wP=1/dP, wQ=1/dQ, wR=1/dR, ws=wP+wQ+wR; mx = (wP*P[0]+wQ*Q[0]+wR*R[0])/ws; my = (wP*P[1]+wQ*Q[1]+wR*R[1])/ws; } return [mx, my]; } const fp = fermatPt(A, B, C); const sectorColors = ["rgba(231,76,60,0.12)", "rgba(39,174,96,0.12)", "rgba(52,152,219,0.12)"]; const sectorStroke = ["#e74c3c", "#27ae60", "#3498db"]; const vLabels = ["a", "b", "c"]; const triVerts = [A, B, C]; // --- Left panel title --- svg.append("text").attr("x", leftCx).attr("y", 16).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", "#555").text("Fermat point & trisector sectors"); // Draw trisector sectors (rays from Fermat point through each vertex) const rayLen = 500; for (let i = 0; i < 3; i++) { const dx = triVerts[i][0] - fp[0], dy = triVerts[i][1] - fp[1]; const len = Math.sqrt(dx*dx+dy*dy) || 1; const endX = fp[0] + dx/len*rayLen, endY = fp[1] + dy/len*rayLen; svg.append("line").attr("x1", fp[0]).attr("y1", fp[1]).attr("x2", endX).attr("y2", endY) .attr("stroke", sectorStroke[i]).attr("stroke-width", 1.5).attr("stroke-dasharray", "6,3").attr("opacity", 0.6); } // Sector fills for (let i = 0; i < 3; i++) { const j = (i+1)%3; const dx1 = triVerts[i][0]-fp[0], dy1 = triVerts[i][1]-fp[1]; const dx2 = triVerts[j][0]-fp[0], dy2 = triVerts[j][1]-fp[1]; const l1 = Math.sqrt(dx1*dx1+dy1*dy1)||1, l2 = Math.sqrt(dx2*dx2+dy2*dy2)||1; svg.append("polygon") .attr("points", `${fp[0]},${fp[1]} ${fp[0]+dx1/l1*rayLen},${fp[1]+dy1/l1*rayLen} ${fp[0]+dx2/l2*rayLen},${fp[1]+dy2/l2*rayLen}`) .attr("fill", sectorColors[i]).attr("opacity", 0.5); } // Triangle svg.append("polygon") .attr("points", triVerts.map(v=>v.join(",")).join(" ")) .attr("fill", "rgba(255,255,255,0.6)") .attr("stroke", needsSplit ? "#e74c3c" : "#2c3e50").attr("stroke-width", 2); // Fermat point const fpIsAtVertex = (Math.abs(fp[0]-C[0]) < 2 && Math.abs(fp[1]-C[1]) < 2) || (Math.abs(fp[0]-A[0]) < 2 && Math.abs(fp[1]-A[1]) < 2) || (Math.abs(fp[0]-B[0]) < 2 && Math.abs(fp[1]-B[1]) < 2); svg.append("circle").attr("cx", fp[0]).attr("cy", fp[1]).attr("r", 5) .attr("fill", fpIsAtVertex ? "#e74c3c" : "#8e44ad"); svg.append("text").attr("x", fp[0]+10).attr("y", fp[1]-8) .attr("font-size", 12).attr("fill", fpIsAtVertex ? "#e74c3c" : "#8e44ad").attr("font-weight", "bold") .text(fpIsAtVertex ? "m (at vertex!)" : "m"); // Angle arcs at Fermat point — show sector angles for (let i = 0; i < 3; i++) { const j = (i+1)%3; const a1 = Math.atan2(triVerts[i][1]-fp[1], triVerts[i][0]-fp[0]); const a2 = Math.atan2(triVerts[j][1]-fp[1], triVerts[j][0]-fp[0]); let da = a2 - a1; while (da < -Math.PI) da += 2*Math.PI; while (da > Math.PI) da -= 2*Math.PI; const angleDeg = Math.abs(da) * 180 / Math.PI; const ar = 25; const mid = a1 + da/2; svg.append("text") .attr("x", fp[0] + (ar+12)*Math.cos(mid)) .attr("y", fp[1] + (ar+12)*Math.sin(mid)) .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#8e44ad") .text(`${angleDeg.toFixed(0)}\u00B0`); } // Vertex labels and angle labels at vertices for (let i = 0; i < 3; i++) { svg.append("circle").attr("cx", triVerts[i][0]).attr("cy", triVerts[i][1]).attr("r", 5) .attr("fill", sectorStroke[i]); const offY = i < 2 ? 18 : -12; const offX = i === 0 ? -14 : i === 1 ? 14 : 0; svg.append("text").attr("x", triVerts[i][0]+offX).attr("y", triVerts[i][1]+offY) .attr("text-anchor", "middle").attr("font-size", 13).attr("fill", sectorStroke[i]).attr("font-weight", "bold") .text(vLabels[i]); } // Angle arc at C const arcR = 20; const aCx = Math.atan2(A[1]-C[1], A[0]-C[0]); const bCx = Math.atan2(B[1]-C[1], B[0]-C[0]); let startAng = Math.min(aCx, bCx), endAng = Math.max(aCx, bCx); if (endAng - startAng > Math.PI) { const tmp = startAng; startAng = endAng; endAng = tmp + 2*Math.PI; } const arcPts = d3.range(startAng, endAng, 0.02).map(a => [C[0]+arcR*Math.cos(a), C[1]+arcR*Math.sin(a)]); if (arcPts.length > 1) { svg.append("path").attr("d", d3.line()(arcPts)) .attr("fill", "none").attr("stroke", "#f39c12").attr("stroke-width", 1.5); } svg.append("text").attr("x", C[0]).attr("y", C[1]+35) .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#f39c12") .text(`${angC.toFixed(0)}\u00B0`); // --- Right panel: after virtual split (or explanation) --- const rightX = leftW + 30; const rightW2 = w - rightX - 10; const rightCx = rightX + rightW2/2; svg.append("text").attr("x", rightCx).attr("y", 16).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", "#555") .text(needsSplit ? "After virtual split" : "No split needed"); if (needsSplit) { // Split triangle along altitude from C to AB const footX = C[0], footY = baseY; // AB is horizontal // Remap to right panel const rA = [rightCx - halfBase, baseY]; const rB = [rightCx + halfBase, baseY]; const rC = [rightCx, baseY - Math.max(25, triH)]; const rD = [rightCx, baseY]; // foot of altitude // Draw the two sub-triangles // Sub-tri 1: A, D, C (left) svg.append("polygon") .attr("points", `${rA[0]},${rA[1]} ${rD[0]},${rD[1]} ${rC[0]},${rC[1]}`) .attr("fill", "rgba(39,174,96,0.08)").attr("stroke", "#27ae60").attr("stroke-width", 2); // Sub-tri 2: D, B, C (right) svg.append("polygon") .attr("points", `${rD[0]},${rD[1]} ${rB[0]},${rB[1]} ${rC[0]},${rC[1]}`) .attr("fill", "rgba(52,152,219,0.08)").attr("stroke", "#3498db").attr("stroke-width", 2); // Altitude line svg.append("line").attr("x1", rC[0]).attr("y1", rC[1]).attr("x2", rD[0]).attr("y2", rD[1]) .attr("stroke", "#e74c3c").attr("stroke-width", 2).attr("stroke-dasharray", "5,3"); // Compute Fermat points for each sub-triangle const fp1 = fermatPt(rA, rD, rC); const fp2 = fermatPt(rD, rB, rC); // Draw sectors for sub-tri 1 const subTris = [[rA, rD, rC, fp1, "#27ae60"], [rD, rB, rC, fp2, "#3498db"]]; for (const [sA, sB, sC, sFp, col] of subTris) { const sVerts = [sA, sB, sC]; const srLen = 200; for (let i = 0; i < 3; i++) { const dx = sVerts[i][0]-sFp[0], dy = sVerts[i][1]-sFp[1]; const len = Math.sqrt(dx*dx+dy*dy)||1; svg.append("line").attr("x1", sFp[0]).attr("y1", sFp[1]) .attr("x2", sFp[0]+dx/len*srLen).attr("y2", sFp[1]+dy/len*srLen) .attr("stroke", col).attr("stroke-width", 1).attr("stroke-dasharray", "4,3").attr("opacity", 0.4); } svg.append("circle").attr("cx", sFp[0]).attr("cy", sFp[1]).attr("r", 4).attr("fill", "#8e44ad"); } // Compute max angle in each sub-triangle const maxAng1 = Math.max(vecAngle(rD,rA,rC), vecAngle(rA,rD,rC), vecAngle(rA,rC,rD)) * 180/Math.PI; const maxAng2 = Math.max(vecAngle(rB,rD,rC), vecAngle(rD,rB,rC), vecAngle(rD,rC,rB)) * 180/Math.PI; // Labels svg.append("circle").attr("cx", rD[0]).attr("cy", rD[1]).attr("r", 5).attr("fill", "#e74c3c"); svg.append("text").attr("x", rD[0]+10).attr("y", rD[1]+15) .attr("font-size", 12).attr("fill", "#e74c3c").attr("font-weight", "bold").text("d"); for (const [pt, lab, offX, offY] of [[rA,"a",-14,18],[rB,"b",14,18],[rC,"c",0,-12]]) { svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 4).attr("fill", "#2c3e50"); svg.append("text").attr("x", pt[0]+offX).attr("y", pt[1]+offY) .attr("text-anchor", "middle").attr("font-size", 12).attr("fill", "#2c3e50").text(lab); } svg.append("text").attr("x", rightCx).attr("y", baseY + 30).attr("text-anchor", "middle") .attr("font-size", 11).attr("fill", "#666") .text(`max angles: ${maxAng1.toFixed(0)}\u00B0 and ${maxAng2.toFixed(0)}\u00B0 (both \u2264 90\u00B0)`); svg.append("text").attr("x", rightCx).attr("y", baseY + 46).attr("text-anchor", "middle") .attr("font-size", 11).attr("fill", "#8e44ad") .text("each sub-triangle gets its own Fermat point"); } else { // No split needed: show the healthy triangle with annotation const rA = [rightCx - halfBase, baseY]; const rB = [rightCx + halfBase, baseY]; const rC = [rightCx, baseY - Math.max(25, triH)]; svg.append("polygon") .attr("points", `${rA[0]},${rA[1]} ${rB[0]},${rB[1]} ${rC[0]},${rC[1]}`) .attr("fill", "rgba(39,174,96,0.08)").attr("stroke", "#27ae60").attr("stroke-width", 2); const rfp = fermatPt(rA, rB, rC); // Sector angles at Fermat point const rVerts = [rA, rB, rC]; for (let i = 0; i < 3; i++) { const j = (i+1)%3; const a1 = Math.atan2(rVerts[i][1]-rfp[1], rVerts[i][0]-rfp[0]); const a2 = Math.atan2(rVerts[j][1]-rfp[1], rVerts[j][0]-rfp[0]); let da = a2 - a1; while (da < -Math.PI) da += 2*Math.PI; while (da > Math.PI) da -= 2*Math.PI; const angleDeg = Math.abs(da) * 180/Math.PI; const ar2 = 30; const mid2 = a1 + da/2; svg.append("text") .attr("x", rfp[0]+(ar2+12)*Math.cos(mid2)) .attr("y", rfp[1]+(ar2+12)*Math.sin(mid2)) .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#8e44ad") .text(`${angleDeg.toFixed(0)}\u00B0`); } svg.append("circle").attr("cx", rfp[0]).attr("cy", rfp[1]).attr("r", 5).attr("fill", "#8e44ad"); svg.append("text").attr("x", rfp[0]+10).attr("y", rfp[1]-8) .attr("font-size", 12).attr("fill", "#8e44ad").text("m"); for (const [pt, lab, offX, offY] of [[rA,"a",-14,18],[rB,"b",14,18],[rC,"c",0,-12]]) { svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 4).attr("fill", "#2c3e50"); svg.append("text").attr("x", pt[0]+offX).attr("y", pt[1]+offY) .attr("text-anchor", "middle").attr("font-size", 12).attr("fill", "#2c3e50").text(lab); } svg.append("text").attr("x", rightCx).attr("y", baseY + 30).attr("text-anchor", "middle") .attr("font-size", 11).attr("fill", "#27ae60") .text(`all sectors \u2248 120\u00B0 \u2014 rich feasible space`); } // Bottom status const statusColor = needsSplit ? "#e74c3c" : obtuse_angle > 90 ? "#f39c12" : "#27ae60"; const statusText = obtuse_angle >= 120 ? `${angC.toFixed(0)}\u00B0 \u2265 120\u00B0 \u2014 Fermat point collapses to vertex c (degenerate!)` : needsSplit ? `${angC.toFixed(0)}\u00B0 > 100\u00B0 \u2014 sectors too narrow, virtual split required` : obtuse_angle > 90 ? `${angC.toFixed(0)}\u00B0 > 90\u00B0 but \u2264 100\u00B0 \u2014 Fermat point still well inside` : `${angC.toFixed(0)}\u00B0 \u2264 90\u00B0 \u2014 ideal, well-balanced sectors`; svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", statusColor).attr("font-weight", "bold") .text(statusText); return svg.node(); } ``` > *"The split point can be expressed in barycentric coordinates and thus is linear in the original vertices. In practice we virtually split all triangles exhibiting an angle larger than $100°$."* After the split, no angle exceeds $90°$, guaranteeing the Fermat point exists inside each sub-triangle. ### 4.5 The Linearized Orientation Constraints Combining the trisector construction with normalization for floating-point precision gives 6 linear inequality constraints per triangle: $$\delta_i^0 \left( (\mathbf{u}_i - \mathbf{m}) \cdot \frac{\mathbf{s}_i^\perp}{\|\mathbf{s}_i^\perp\|} - \epsilon \right) \geq 0$$ $$\delta_i^1 \left( -(\mathbf{u}_i - \mathbf{m}) \cdot \frac{\mathbf{s}_{i+1}^\perp}{\|\mathbf{s}_{i+1}^\perp\|} - \epsilon \right) \geq 0 \tag{4}$$ where $2\epsilon$ is a lower bound on the distance between two vertices of the same triangle, and $\delta_i^0, \delta_i^1$ are normalization factors. This gives us a **convex formulation** — linear constraints with a quadratic objective — which modern MIQP solvers can handle. --- ## 5. Complexity Reduction ### 5.1 The Problem > *"One reason for the enormous dimension of discrete variables in problem (P1) is that potentially every vertex in the input triangle mesh can represent a singularity in the quad mesh."* A typical input mesh has tens of thousands of vertices, each contributing discrete variables. Branch-and-cut solvers struggle with search spaces this large. ### 5.2 The Key Insight: Decimate in Parametric Space The idea is to exploit the continuous relaxation $f_{\text{relaxed}}$ (the solution ignoring integer constraints). In flat regions far from singularities, the continuous relaxation is already close to an integer solution. Only near singularities does the integer snapping matter. > *"Accordingly the idea is to perform decimation in flat parametric space and thus being able to represent the input surface by an extremely coarse mesh."* The decimation is performed **in the 2D parametric space** of $\mathcal{M}_{\text{relaxed}}$, not in 3D. After decimation, each remaining triangle in the coarse mesh $\mathcal{M}_-$ overlaps with a set of triangles from $\mathcal{M}_{\text{relaxed}}$ in parameter space. ### 5.3 Decimation Operations The mesh decimation uses three standard operations, all performed as an isometric re-tessellation of flat 2D space: 1. **Halfedge collapse** — merge an edge, removing one vertex (sorted by increasing edge length) 2. **Edge flip** — swap the diagonal of a quad formed by two triangles (to restore Delaunay criterion) 3. **Edge split** — subdivide an edge into two Operations that move cone singularities are forbidden. This preserves the critical topological structure while aggressively simplifying flat regions. Click an edge to apply the selected operation. **Collapse** merges the edge's two endpoints into one, **Flip** swaps the shared diagonal between two adjacent triangles, and **Split** inserts a new vertex at the edge midpoint. The star vertex (orange) is a singularity and cannot be moved or removed: ```{ojs} //| label: fig-decimation-ops viewof decim_op = Inputs.radio(["Halfedge collapse", "Edge flip", "Edge split"], {value: "Halfedge collapse", label: "Operation"}) ``` ```{ojs} //| code-fold: true { 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 node = svg.node(); // --- Mesh data structure (stored on SVG node for persistence) --- if (!node.__mesh || decim_op !== node.__lastOp) { // Reinitialize mesh when operation changes, to start fresh // Vertices: a structured grid-ish mesh with one singularity const pts = [ [100, 460], [250, 460], [400, 460], [550, 460], // 0-3 bottom row [100, 320], [250, 310], [400, 320], [550, 320], // 4-7 mid row [175, 180], [350, 170], [525, 190], // 8-10 upper row [350, 60], // 11 top ]; // Triangles as vertex-index triples (CCW in SVG coords = CW math, but that's fine for viz) const tris = [ [0, 1, 4], [1, 5, 4], [1, 2, 5], [2, 6, 5], [2, 3, 6], [3, 7, 6], [4, 5, 8], [5, 9, 8], [5, 6, 9], [6, 10, 9], [6, 7, 10], [8, 9, 11], [9, 10, 11], ]; node.__mesh = { verts: pts.map(p => [...p]), tris: tris.map(t => [...t]), singularity: 9, // vertex 9 is a cone singularity message: "Click an edge to apply the operation", messageColor: "#555", }; node.__lastOp = decim_op; } const mesh = node.__mesh; const verts = mesh.verts; const tris = mesh.tris; const singIdx = mesh.singularity; // --- Edge utilities --- function edgeKey(a, b) { return a < b ? `${a}-${b}` : `${b}-${a}`; } // Build edge-to-triangles map function buildEdgeMap() { const em = new Map(); for (let ti = 0; ti < tris.length; ti++) { const t = tris[ti]; for (let j = 0; j < 3; j++) { const k = edgeKey(t[j], t[(j+1)%3]); if (!em.has(k)) em.set(k, []); em.get(k).push(ti); } } return em; } // --- Operations --- function doCollapse(a, b) { // Collapse edge a-b: merge b into a (keep a, remove b) // If b is the singularity, swap so we keep the singularity if (b === singIdx) { const tmp = a; a = b; b = tmp; } // Cannot collapse if a is the singularity (it would move it) if (a === singIdx) { mesh.message = "Cannot collapse: would move the singularity vertex"; mesh.messageColor = "#e74c3c"; return; } // Move a to midpoint verts[a] = [(verts[a][0]+verts[b][0])/2, (verts[a][1]+verts[b][1])/2]; // Replace all references to b with a for (const t of tris) { for (let j = 0; j < 3; j++) { if (t[j] === b) t[j] = a; } } // Remove degenerate triangles (where two or more indices are the same) for (let i = tris.length - 1; i >= 0; i--) { const t = tris[i]; if (t[0] === t[1] || t[1] === t[2] || t[0] === t[2]) tris.splice(i, 1); } // Update singularity index if needed (it shouldn't change since we kept it) mesh.message = `Collapsed edge: merged vertex ${b} into ${a}`; mesh.messageColor = "#2980b9"; } function doFlip(a, b) { const em = buildEdgeMap(); const k = edgeKey(a, b); const adjTris = em.get(k); if (!adjTris || adjTris.length !== 2) { mesh.message = "Cannot flip: edge is on the boundary"; mesh.messageColor = "#e74c3c"; return; } // Find the two opposite vertices const t0 = tris[adjTris[0]], t1 = tris[adjTris[1]]; let opp0 = -1, opp1 = -1; for (let j = 0; j < 3; j++) { if (t0[j] !== a && t0[j] !== b) opp0 = t0[j]; if (t1[j] !== a && t1[j] !== b) opp1 = t1[j]; } if (opp0 === -1 || opp1 === -1 || opp0 === opp1) { mesh.message = "Cannot flip: degenerate configuration"; mesh.messageColor = "#e74c3c"; return; } // Check convexity of the quad (opp0, a, opp1, b) — flip only if convex function cross(O, A, B) { return (A[0]-O[0])*(B[1]-O[1]) - (A[1]-O[1])*(B[0]-O[0]); } const quad = [opp0, a, opp1, b].map(i => verts[i]); let convex = true; for (let i = 0; i < 4; i++) { const c = cross(quad[i], quad[(i+1)%4], quad[(i+2)%4]); if (c >= 0) { convex = false; break; } // SVG coords: CW is positive } if (!convex) { mesh.message = "Cannot flip: quad is not convex"; mesh.messageColor = "#e74c3c"; return; } // Replace the two triangles with the flipped pair tris[adjTris[0]] = [opp0, opp1, a]; tris[adjTris[1]] = [opp0, b, opp1]; mesh.message = `Flipped edge: (${a},${b}) \u2192 (${opp0},${opp1})`; mesh.messageColor = "#2980b9"; } function doSplit(a, b) { const em = buildEdgeMap(); const k = edgeKey(a, b); const adjTris = em.get(k) || []; // Insert new vertex at midpoint const newIdx = verts.length; verts.push([(verts[a][0]+verts[b][0])/2, (verts[a][1]+verts[b][1])/2]); // Replace each adjacent triangle with two sub-triangles const toRemove = new Set(adjTris); const newTris = []; for (const ti of adjTris) { const t = tris[ti]; // Find the opposite vertex let opp = -1; for (let j = 0; j < 3; j++) { if (t[j] !== a && t[j] !== b) opp = t[j]; } newTris.push([a, newIdx, opp]); newTris.push([newIdx, b, opp]); } // Remove old triangles (in reverse order) and add new ones const sorted = [...toRemove].sort((x, y) => y - x); for (const ti of sorted) tris.splice(ti, 1); for (const nt of newTris) tris.push(nt); mesh.message = `Split edge: inserted vertex ${newIdx} at midpoint of (${a},${b})`; mesh.messageColor = "#2980b9"; } // --- Drawing --- function redraw() { svg.selectAll("*").remove(); svg.append("rect").attr("width", w).attr("height", h).attr("fill", "#fafafa"); const em = buildEdgeMap(); // Draw triangles for (const t of tris) { const pts = t.map(i => verts[i]); svg.append("polygon") .attr("points", pts.map(p => p.join(",")).join(" ")) .attr("fill", "rgba(52,152,219,0.06)").attr("stroke", "#bdc3c7").attr("stroke-width", 0.5); } // Draw edges as clickable hit targets const drawnEdges = new Set(); for (const t of tris) { for (let j = 0; j < 3; j++) { const a = t[j], b = t[(j+1)%3]; const k = edgeKey(a, b); if (drawnEdges.has(k)) continue; drawnEdges.add(k); const isBoundary = (em.get(k) || []).length === 1; // Visible edge svg.append("line") .attr("x1", verts[a][0]).attr("y1", verts[a][1]) .attr("x2", verts[b][0]).attr("y2", verts[b][1]) .attr("stroke", isBoundary ? "#7f8c8d" : "#2c3e50") .attr("stroke-width", isBoundary ? 1.5 : 1); // Wider invisible hit target svg.append("line") .attr("x1", verts[a][0]).attr("y1", verts[a][1]) .attr("x2", verts[b][0]).attr("y2", verts[b][1]) .attr("stroke", "transparent").attr("stroke-width", 14) .attr("cursor", "pointer") .attr("data-a", a).attr("data-b", b) .on("mouseenter", function() { d3.select(this.previousSibling) .attr("stroke", "#e67e22").attr("stroke-width", 3); }) .on("mouseleave", function() { const ib = (em.get(edgeKey(+d3.select(this).attr("data-a"), +d3.select(this).attr("data-b"))) || []).length === 1; d3.select(this.previousSibling) .attr("stroke", ib ? "#7f8c8d" : "#2c3e50") .attr("stroke-width", ib ? 1.5 : 1); }) .on("click", function() { const ea = +d3.select(this).attr("data-a"); const eb = +d3.select(this).attr("data-b"); if (decim_op === "Halfedge collapse") doCollapse(ea, eb); else if (decim_op === "Edge flip") doFlip(ea, eb); else doSplit(ea, eb); redraw(); }); } } // Draw vertices const usedVerts = new Set(); for (const t of tris) for (const v of t) usedVerts.add(v); for (const vi of usedVerts) { const isSing = vi === singIdx; svg.append("circle") .attr("cx", verts[vi][0]).attr("cy", verts[vi][1]) .attr("r", isSing ? 7 : 4) .attr("fill", isSing ? "#e67e22" : "#2c3e50") .attr("stroke", isSing ? "#d35400" : "none") .attr("stroke-width", isSing ? 2 : 0); } // Singularity label if (usedVerts.has(singIdx)) { svg.append("text") .attr("x", verts[singIdx][0] + 12).attr("y", verts[singIdx][1] - 10) .attr("font-size", 12).attr("fill", "#e67e22").attr("font-weight", "bold") .text("singularity"); } // Stats svg.append("text").attr("x", w/2).attr("y", h - 30).attr("text-anchor", "middle") .attr("font-size", 12).attr("fill", "#888") .text(`${usedVerts.size} vertices \u00B7 ${tris.length} triangles \u00B7 ${drawnEdges.size} edges`); // Message svg.append("text").attr("x", w/2).attr("y", h - 8).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", mesh.messageColor).attr("font-weight", "bold") .text(mesh.message); // Legend svg.append("circle").attr("cx", 20).attr("cy", 18).attr("r", 5).attr("fill", "#e67e22").attr("stroke", "#d35400").attr("stroke-width", 1.5); svg.append("text").attr("x", 32).attr("y", 22).attr("font-size", 11).attr("fill", "#555").text("Singularity (protected)"); svg.append("text").attr("x", w - 20).attr("y", 22).attr("text-anchor", "end") .attr("font-size", 11).attr("fill", "#999").text("Hover an edge, then click to apply"); } redraw(); return node; } ``` ### 5.4 Splitting Integer-Critical Edges After decimation, some triangles with large interior angles can trap vertices in degenerate line configurations. An edge is **integer-critical** if it is incident to a triangle with three singular vertices where the angle opposing the edge exceeds $110°$. > *"To not exclude the possibility of arranging the vertices in the integer-grid on a line, after decimation we split all integer-critical edges."* Each split separates two singular vertices by a regular one, and the number of integer-critical edges decreases monotonically. The visualization below shows a triangle with three singular vertices after decimation. When the opposing angle is large, the only integer arrangement that fits all three vertices is **collinear** — a degenerate triangle. Splitting the integer-critical edge inserts a regular vertex that breaks the degeneracy: ```{ojs} //| label: fig-integer-critical viewof ic_opposing_angle = Inputs.range([70, 150], {value: 90, step: 1, label: "Opposing angle (degrees)"}) ``` ```{ojs} //| code-fold: true { const w = 720, h = 460; const svg = d3.create("svg").attr("viewBox", `0 0 ${w} ${h}`).style("width", "100%").style("max-width", w + "px"); const isCritical = ic_opposing_angle > 110; const leftW = w * 0.48, gap = w * 0.04, rightX = leftW + gap; const rightW = w - rightX - 5; // --- Build a triangle with the given opposing angle at vertex C --- // Edge AB is the "integer-critical edge" (the one opposite to the large angle at C) const baseLen = 160; const baseY = 310; const angleRad = ic_opposing_angle * Math.PI / 180; const triH = (baseLen / 2) / Math.tan(angleRad / 2); // --- LEFT PANEL: The triangle in parametric space --- const leftCx = leftW / 2 + 10; const A = [leftCx - baseLen/2, baseY]; const B = [leftCx + baseLen/2, baseY]; const C = [leftCx, baseY - Math.max(20, triH)]; svg.append("text").attr("x", leftCx).attr("y", 16).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", "#555").text("Triangle with 3 singular vertices"); // Integer grid behind const gridSc = 55; const gridOx = leftCx - 3*gridSc, gridOy = baseY + gridSc; for (let i = -1; i <= 6; i++) { svg.append("line").attr("x1", gridOx+i*gridSc).attr("y1", 25).attr("x2", gridOx+i*gridSc).attr("y2", baseY+40) .attr("stroke", "#eee").attr("stroke-width", 0.5); } for (let j = -1; j <= 6; j++) { svg.append("line").attr("x1", gridOx-gridSc).attr("y1", gridOy-j*gridSc).attr("x2", gridOx+6*gridSc).attr("y2", gridOy-j*gridSc) .attr("stroke", "#eee").attr("stroke-width", 0.5); } // Integer lattice dots for (let i = 0; i <= 5; i++) { for (let j = 0; j <= 5; j++) { svg.append("circle").attr("cx", gridOx+i*gridSc).attr("cy", gridOy-j*gridSc).attr("r", 2).attr("fill", "#ddd"); } } // Triangle svg.append("polygon") .attr("points", `${A[0]},${A[1]} ${B[0]},${B[1]} ${C[0]},${C[1]}`) .attr("fill", isCritical ? "rgba(231,76,60,0.08)" : "rgba(39,174,96,0.06)") .attr("stroke", isCritical ? "#e74c3c" : "#27ae60").attr("stroke-width", 2); // Highlight the integer-critical edge (AB) — the one opposite to the large angle svg.append("line").attr("x1", A[0]).attr("y1", A[1]).attr("x2", B[0]).attr("y2", B[1]) .attr("stroke", isCritical ? "#e74c3c" : "#27ae60") .attr("stroke-width", isCritical ? 4 : 2.5); if (isCritical) { const midAB = [(A[0]+B[0])/2, (A[1]+B[1])/2]; svg.append("text").attr("x", midAB[0]).attr("y", midAB[1]+20) .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#e74c3c").attr("font-weight", "bold") .text("integer-critical edge"); } // Angle arc at C function vecAngle(p1, vertex, p2) { const v1x=p1[0]-vertex[0], v1y=p1[1]-vertex[1]; const v2x=p2[0]-vertex[0], v2y=p2[1]-vertex[1]; const dot=v1x*v2x+v1y*v2y; const m1=Math.sqrt(v1x*v1x+v1y*v1y), m2=Math.sqrt(v2x*v2x+v2y*v2y); return Math.acos(Math.max(-1,Math.min(1,dot/(m1*m2)))); } const angC = vecAngle(A, C, B) * 180 / Math.PI; const arcR = 22; const aCx = Math.atan2(A[1]-C[1], A[0]-C[0]); const bCx = Math.atan2(B[1]-C[1], B[0]-C[0]); let startA = Math.min(aCx,bCx), endA = Math.max(aCx,bCx); if (endA - startA > Math.PI) { const tmp=startA; startA=endA; endA=tmp+2*Math.PI; } const arcPts = d3.range(startA, endA, 0.02).map(a => [C[0]+arcR*Math.cos(a), C[1]+arcR*Math.sin(a)]); if (arcPts.length > 1) { svg.append("path").attr("d", d3.line()(arcPts)) .attr("fill", "none").attr("stroke", "#f39c12").attr("stroke-width", 1.5); } svg.append("text").attr("x", C[0]).attr("y", C[1]+38) .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#f39c12") .text(`${angC.toFixed(0)}\u00B0`); // Vertex labels — all are singularities (red dots) for (const [pt, lab, offX, offY] of [[A,"s\u2081",-16,16],[B,"s\u2082",16,16],[C,"s\u2083",0,-14]]) { svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 6).attr("fill", "#e74c3c"); svg.append("text").attr("x", pt[0]+offX).attr("y", pt[1]+offY) .attr("text-anchor", "middle").attr("font-size", 13).attr("fill", "#e74c3c").attr("font-weight", "bold").text(lab); } // --- RIGHT PANEL: Integer grid snapping --- svg.append("text").attr("x", rightX + rightW/2).attr("y", 16).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", "#555") .text(isCritical ? "After splitting the critical edge" : "Integer snapping (no problem)"); // Right panel integer grid const rSc = 70; const rOx = rightX + 30, rOy = baseY + 10; for (let i = 0; i <= 4; i++) { for (let j = 0; j <= 4; j++) { svg.append("circle").attr("cx", rOx+i*rSc).attr("cy", rOy-j*rSc).attr("r", 3).attr("fill", "#ccc"); } } for (let i = 0; i <= 4; i++) { svg.append("line").attr("x1", rOx+i*rSc).attr("y1", rOy+10).attr("x2", rOx+i*rSc).attr("y2", rOy-4*rSc-10) .attr("stroke", "#e8e8e8").attr("stroke-width", 0.5); svg.append("line").attr("x1", rOx-10).attr("y1", rOy-i*rSc).attr("x2", rOx+4*rSc+10).attr("y2", rOy-i*rSc) .attr("stroke", "#e8e8e8").attr("stroke-width", 0.5); } if (isCritical) { // Show: before split (collinear), after split (non-degenerate) // The 3 singularities snapped to integer grid — forced collinear by the flat triangle const colA = [rOx + 0*rSc, rOy - 1*rSc]; const colB = [rOx + 3*rSc, rOy - 1*rSc]; const colC = [rOx + 1*rSc, rOy - 1*rSc]; // on the same line! // Ghost: the degenerate collinear arrangement svg.append("line").attr("x1", colA[0]).attr("y1", colA[1]).attr("x2", colB[0]).attr("y2", colB[1]) .attr("stroke", "rgba(231,76,60,0.3)").attr("stroke-width", 2); svg.append("text").attr("x", (colA[0]+colB[0])/2).attr("y", colA[1]-12) .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "rgba(231,76,60,0.5)") .text("collinear (degenerate!)"); for (const [pt,lab] of [[colA,"s\u2081"],[colB,"s\u2082"],[colC,"s\u2083"]]) { svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 5) .attr("fill", "none").attr("stroke", "rgba(231,76,60,0.4)").attr("stroke-width", 1.5); } // After split: midpoint of AB is now a regular vertex, free to move off the line const splitA = [rOx + 0*rSc, rOy - 2*rSc]; const splitB = [rOx + 3*rSc, rOy - 2*rSc]; const splitC = [rOx + 1*rSc, rOy - 3*rSc]; const splitD = [rOx + 1.5*rSc, rOy - 2*rSc]; // regular vertex (midpoint of AB), NOT on integer grid // Two sub-triangles: ADC and DBC svg.append("polygon") .attr("points", `${splitA[0]},${splitA[1]} ${splitD[0]},${splitD[1]} ${splitC[0]},${splitC[1]}`) .attr("fill", "rgba(39,174,96,0.08)").attr("stroke", "#27ae60").attr("stroke-width", 2); svg.append("polygon") .attr("points", `${splitD[0]},${splitD[1]} ${splitB[0]},${splitB[1]} ${splitC[0]},${splitC[1]}`) .attr("fill", "rgba(52,152,219,0.08)").attr("stroke", "#3498db").attr("stroke-width", 2); // Singular vertices (must be on integers) for (const [pt,lab] of [[splitA,"s\u2081"],[splitB,"s\u2082"],[splitC,"s\u2083"]]) { svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 6).attr("fill", "#e74c3c"); svg.append("text").attr("x", pt[0]).attr("y", pt[1]-12) .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#e74c3c").attr("font-weight", "bold").text(lab); } // Regular vertex (free to be non-integer) svg.append("circle").attr("cx", splitD[0]).attr("cy", splitD[1]).attr("r", 5).attr("fill", "#3498db"); svg.append("text").attr("x", splitD[0]+12).attr("y", splitD[1]+4) .attr("font-size", 11).attr("fill", "#3498db").attr("font-weight", "bold").text("d (regular)"); svg.append("text").attr("x", rightX + rightW/2).attr("y", baseY + 38).attr("text-anchor", "middle") .attr("font-size", 11).attr("fill", "#27ae60") .text("regular vertex d is free to move off the integer grid"); svg.append("text").attr("x", rightX + rightW/2).attr("y", baseY + 54).attr("text-anchor", "middle") .attr("font-size", 11).attr("fill", "#27ae60") .text("\u2192 non-degenerate triangles are now possible"); } else { // Non-critical: show a healthy integer arrangement const snapA = [rOx + 0*rSc, rOy - 1*rSc]; const snapB = [rOx + 3*rSc, rOy - 1*rSc]; const snapC = [rOx + 1*rSc, rOy - 2*rSc]; svg.append("polygon") .attr("points", `${snapA[0]},${snapA[1]} ${snapB[0]},${snapB[1]} ${snapC[0]},${snapC[1]}`) .attr("fill", "rgba(39,174,96,0.08)").attr("stroke", "#27ae60").attr("stroke-width", 2); for (const [pt,lab] of [[snapA,"s\u2081"],[snapB,"s\u2082"],[snapC,"s\u2083"]]) { svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 6).attr("fill", "#e74c3c"); svg.append("text").attr("x", pt[0]).attr("y", pt[1]-12) .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#e74c3c").attr("font-weight", "bold").text(lab); } svg.append("text").attr("x", rightX + rightW/2).attr("y", baseY + 38).attr("text-anchor", "middle") .attr("font-size", 11).attr("fill", "#27ae60") .text("all 3 singularities snap to Z\u00B2 without collapsing"); } // Status bar const statusColor = isCritical ? "#e74c3c" : "#27ae60"; const statusText = isCritical ? `Opposing angle ${angC.toFixed(0)}\u00B0 > 110\u00B0 \u2014 integer-critical! 3 singular vertices forced collinear on Z\u00B2` : `Opposing angle ${angC.toFixed(0)}\u00B0 \u2264 110\u00B0 \u2014 non-collinear integer arrangement exists`; svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", statusColor).attr("font-weight", "bold") .text(statusText); return svg.node(); } ``` ::: {.callout-tip} ## Why 110° is the threshold Three integer-grid points form a non-degenerate triangle only if they are not collinear. When a triangle has all three vertices as singularities (constrained to $\mathbb{Z}^2$) and one angle exceeds $110°$, the triangle is so flat that the nearest valid integer arrangement places all three on a line. Splitting the edge opposite the obtuse angle inserts a **regular** vertex (not constrained to integers) that can sit freely off the integer grid, breaking the collinearity trap. ::: ### 5.5 Maintaining the Map During decimation, a function $\Psi_i : \mathbb{R}^2 \to T_{\mathcal{M}_{\text{relaxed}}} \times \mathbb{R}^2$ tracks the correspondence between the decimated mesh and the original. For each triangle $t_i$ in the coarse mesh, $\Psi_i(\mathbf{m}) = (t_j, \psi_i(\mathbf{m}))$ records which original triangle $t_j$ contains the barycenter and the local coordinates. The following step-through shows the decimation pipeline on a simple example: <div id="decimate-info" style="font-family: monospace; font-size: 13px; padding: 8px 12px; background: var(--bs-tertiary-bg, #f8f9fa); border-radius: 4px; margin-bottom: 8px;"> Step 1/4: Input mesh in parametric space </div> <div style="margin-bottom: 8px;"> <button id="step-input" class="btn btn-sm btn-primary" onclick="window.decimateStep('input')">1. Input Mesh</button> <button id="step-relax" class="btn btn-sm btn-outline-secondary" onclick="window.decimateStep('relax')" disabled>2. Continuous Relaxation</button> <button id="step-decimate" class="btn btn-sm btn-outline-secondary" onclick="window.decimateStep('decimate')" disabled>3. Decimate</button> <button id="step-split" class="btn btn-sm btn-outline-secondary" onclick="window.decimateStep('split')" disabled>4. Split Critical</button> </div> <canvas id="decimate-canvas" style="width: 100%; max-width: 600px; border: 1px solid var(--bs-border-color, #dee2e6); border-radius: 4px;"></canvas> ```{=html} <script> (function() { const canvas = document.getElementById('decimate-canvas'); const ctx = canvas.getContext('2d'); const info = document.getElementById('decimate-info'); function resizeCanvas() { const rect = canvas.parentElement.getBoundingClientRect(); const w = Math.min(rect.width, 600); canvas.style.width = w + 'px'; canvas.style.height = (w * 0.7) + 'px'; canvas.width = w * devicePixelRatio; canvas.height = w * 0.7 * devicePixelRatio; ctx.setTransform(devicePixelRatio, 0, 0, devicePixelRatio, 0, 0); } let currentStep = 0; const steps = ['input', 'relax', 'decimate', 'split']; const stepNames = ['Input mesh in parametric space', 'Continuous relaxation (solve ignoring integers)', 'Decimate flat regions (keep singularities)', 'Split integer-critical edges']; // Generate a simple mesh const gridN = 6; const verts = []; const tris = []; const sings = [14, 22]; // singularity vertex indices for (let j = 0; j <= gridN; j++) { for (let i = 0; i <= gridN; i++) { verts.push([i / gridN, j / gridN]); } } for (let j = 0; j < gridN; j++) { for (let i = 0; i < gridN; i++) { const v0 = j*(gridN+1)+i, v1 = v0+1, v2 = v0+(gridN+1), v3 = v2+1; tris.push([v0, v1, v3]); tris.push([v0, v3, v2]); } } // Relaxed positions (slightly perturbed near singularities) const relaxedVerts = verts.map((v, i) => { let dx = 0, dy = 0; for (const si of sings) { const sv = verts[si]; const dist = Math.sqrt((v[0]-sv[0])**2 + (v[1]-sv[1])**2); if (dist > 0 && dist < 0.3) { const factor = 0.03 * (0.3 - dist) / 0.3; dx += (v[0]-sv[0]) / dist * factor; dy += (v[1]-sv[1]) / dist * factor; } } return [v[0] + dx, v[1] + dy]; }); // Decimated: only keep vertices near singularities + boundary const decimatedKeep = new Set(); for (const si of sings) decimatedKeep.add(si); // Keep boundary vertices for (let i = 0; i <= gridN; i++) { decimatedKeep.add(i); decimatedKeep.add(gridN*(gridN+1)+i); decimatedKeep.add(i*(gridN+1)); decimatedKeep.add(i*(gridN+1)+gridN); } // Keep vertices adjacent to singularities for (const si of sings) { const row = Math.floor(si / (gridN+1)), col = si % (gridN+1); for (let dr = -1; dr <= 1; dr++) { for (let dc = -1; dc <= 1; dc++) { const r2 = row+dr, c2 = col+dc; if (r2 >= 0 && r2 <= gridN && c2 >= 0 && c2 <= gridN) { decimatedKeep.add(r2*(gridN+1)+c2); } } } } function draw() { const cw = canvas.width / devicePixelRatio; const ch = canvas.height / devicePixelRatio; const bg = getComputedStyle(document.documentElement).getPropertyValue('--bs-body-bg') || '#fff'; const fg = getComputedStyle(document.documentElement).getPropertyValue('--bs-body-color') || '#333'; ctx.fillStyle = bg.trim() || '#fff'; ctx.fillRect(0, 0, cw, ch); const margin = 40; const sx = cw - 2*margin, sy = ch - 2*margin; function toScreen(u, v) { return [margin + u * sx, ch - margin - v * sy]; } // Integer grid ctx.strokeStyle = '#e8e8e8'; ctx.lineWidth = 0.5; for (let i = 0; i <= 1; i += 1/gridN) { const [x1, y1] = toScreen(i, 0); const [x2, y2] = toScreen(i, 1); ctx.beginPath(); ctx.moveTo(x1, y1); ctx.lineTo(x2, y2); ctx.stroke(); const [x3, y3] = toScreen(0, i); const [x4, y4] = toScreen(1, i); ctx.beginPath(); ctx.moveTo(x3, y3); ctx.lineTo(x4, y4); ctx.stroke(); } const useVerts = currentStep >= 1 ? relaxedVerts : verts; const showAll = currentStep < 2; // Draw triangles ctx.strokeStyle = 'rgba(150,150,150,0.4)'; ctx.lineWidth = 0.8; for (const tri of tris) { if (!showAll) { // In decimated mode, only show tris where all vertices are kept if (!tri.every(v => decimatedKeep.has(v))) continue; } const [x0, y0] = toScreen(useVerts[tri[0]][0], useVerts[tri[0]][1]); const [x1, y1] = toScreen(useVerts[tri[1]][0], useVerts[tri[1]][1]); const [x2, y2] = toScreen(useVerts[tri[2]][0], useVerts[tri[2]][1]); ctx.beginPath(); ctx.moveTo(x0, y0); ctx.lineTo(x1, y1); ctx.lineTo(x2, y2); ctx.closePath(); ctx.stroke(); } // Draw vertices for (let i = 0; i < useVerts.length; i++) { if (!showAll && !decimatedKeep.has(i)) continue; const [x, y] = toScreen(useVerts[i][0], useVerts[i][1]); const isSing = sings.includes(i); ctx.beginPath(); ctx.arc(x, y, isSing ? 5 : 2, 0, 2*Math.PI); ctx.fillStyle = isSing ? '#e74c3c' : '#2c3e50'; ctx.fill(); } // Split critical edges (step 4): highlight edges near singularities if (currentStep >= 3) { ctx.strokeStyle = '#27ae60'; ctx.lineWidth = 2; for (const si of sings) { const row = Math.floor(si / (gridN+1)), col = si % (gridN+1); // Draw split markers on edges adjacent to singularity const neighbors = []; if (col > 0) neighbors.push(si - 1); if (col < gridN) neighbors.push(si + 1); if (row > 0) neighbors.push(si - (gridN+1)); if (row < gridN) neighbors.push(si + (gridN+1)); for (const ni of neighbors) { const [x1, y1] = toScreen(useVerts[si][0], useVerts[si][1]); const [x2, y2] = toScreen(useVerts[ni][0], useVerts[ni][1]); const mx = (x1+x2)/2, my = (y1+y2)/2; ctx.beginPath(); ctx.arc(mx, my, 3, 0, 2*Math.PI); ctx.fillStyle = '#27ae60'; ctx.fill(); } } } // Vertex count const vertCount = showAll ? useVerts.length : [...decimatedKeep].length; ctx.fillStyle = fg.trim() || '#333'; ctx.font = '12px monospace'; ctx.textAlign = 'right'; ctx.fillText(`Vertices: ${vertCount}`, cw - 10, 18); } const stepBtns = steps.map(s => document.getElementById('step-' + s)); window.decimateStep = function(step) { const idx = steps.indexOf(step); if (idx < 0) return; currentStep = idx; info.textContent = `Step ${idx+1}/${steps.length}: ${stepNames[idx]}`; for (let i = 0; i <= idx; i++) { stepBtns[i].classList.remove('btn-outline-secondary'); stepBtns[i].classList.add('btn-primary'); stepBtns[i].disabled = false; } if (idx + 1 < steps.length) { stepBtns[idx+1].disabled = false; } draw(); }; resizeCanvas(); window.addEventListener('resize', () => { resizeCanvas(); draw(); }); draw(); })(); </script> ``` --- ## 6. Singularity Separation ### 6.1 The Clustering Problem > *"The reason is that singular vertices might get arbitrarily close to each other within the continuous relaxation, while their distance in the integer-grid $\mathbb{Z} \times \mathbb{Z}$ of a valid IGM can never be less than 1."* When the continuous relaxation places two singularities very close together, the branch-and-cut solver must search a huge space to find an integer solution that separates them. This makes the solver spend most of its time exploring infeasible subtrees. ### 6.2 The Geodesic Disc Strategy For each singular vertex $s_i$, grow a geodesic disc $G_i$ of radius $\sqrt{2}$ in the decimated mesh $\mathcal{M}_-$. For all pairs of singularities $(s_i, s_j)$ within each disc ($i < j$), trace the shortest geodesic path connecting them and add a **linear separation constraint**. The geodesic distances are computed using the **Fast Marching method** driven by the circle update of Novotni and Klein (2002), which is exact for obstacle-free planar cases. ### 6.3 Separator Constraints Along each geodesic path between a close pair of singularities, the separation is enforced in the axis with the larger difference: $$\left\{ \text{sgn}(\Delta u)|_{\mathcal{M}_-} \right\} \cdot \Delta u \geq 1$$ $$\text{or} \quad \left\{ \text{sgn}(\Delta v)|_{\mathcal{M}_-} \right\} \cdot \Delta v \geq 1 \tag{5}$$ where separation is performed along the $u$ axis if $|\Delta u| \geq |\Delta v|$ and along $v$ otherwise. The sign is fixed from the continuous relaxation, making the constraint linear. Adjust the target edge length to see how singularity spacing affects the need for separation constraints: ```{ojs} //| label: fig-singularity-separation viewof sep_scale = Inputs.range([0.5, 3.0], {value: 1.5, step: 0.1, label: "Target edge length"}) ``` ```{ojs} //| code-fold: true { const w = 550, h = 450; const svg = d3.create("svg").attr("viewBox", `0 0 ${w} ${h}`).style("width", "100%").style("max-width", w + "px"); const ox = 50, oy = 400, sc = 80; // Integer grid for (let i = 0; i <= 5; i++) { svg.append("line").attr("x1", ox + i*sc).attr("y1", oy).attr("x2", ox + i*sc).attr("y2", oy - 4*sc) .attr("stroke", "#eee").attr("stroke-width", 0.5); svg.append("line").attr("x1", ox).attr("y1", oy - i*sc).attr("x2", ox + 5*sc).attr("y2", oy - i*sc) .attr("stroke", "#eee").attr("stroke-width", 0.5); } // Integer points for (let ix = 0; ix <= 5; ix++) { for (let iy = 0; iy <= 4; iy++) { svg.append("circle").attr("cx", ox + ix*sc).attr("cy", oy - iy*sc).attr("r", 2).attr("fill", "#ddd"); } } // Singularities — positions scale with target edge length const baseSings = [ {x: 1.5, y: 1.8, label: "s\u2081"}, {x: 2.2, y: 2.1, label: "s\u2082"}, {x: 3.5, y: 1.0, label: "s\u2083"}, {x: 3.8, y: 2.8, label: "s\u2084"}, ]; // As edge length decreases, singularities get closer in the parametric domain const sings = baseSings.map(s => ({ x: s.x, y: s.y, label: s.label })); const discRadius = Math.sqrt(2); // Draw geodesic discs for (const s of sings) { const cx = ox + s.x * sc, cy = oy - s.y * sc; svg.append("circle") .attr("cx", cx).attr("cy", cy).attr("r", discRadius * sc) .attr("fill", "rgba(231,76,60,0.06)") .attr("stroke", "rgba(231,76,60,0.3)").attr("stroke-width", 1).attr("stroke-dasharray", "4,3"); } // Check pairs and draw separators let numSeparators = 0; for (let i = 0; i < sings.length; i++) { for (let j = i + 1; j < sings.length; j++) { const dx = (sings[j].x - sings[i].x) * sep_scale; const dy = (sings[j].y - sings[i].y) * sep_scale; const dist = Math.sqrt(dx*dx + dy*dy); const px1 = ox + sings[i].x * sc, py1 = oy - sings[i].y * sc; const px2 = ox + sings[j].x * sc, py2 = oy - sings[j].y * sc; if (dist < discRadius * sep_scale) { // Need separator numSeparators++; svg.append("line") .attr("x1", px1).attr("y1", py1).attr("x2", px2).attr("y2", py2) .attr("stroke", "#e67e22").attr("stroke-width", 2.5); // Separation direction indicator const midX = (px1+px2)/2, midY = (py1+py2)/2; const sepAxis = Math.abs(dx) >= Math.abs(dy) ? "u" : "v"; svg.append("text").attr("x", midX + 8).attr("y", midY - 5) .attr("font-size", 10).attr("fill", "#e67e22").text(`sep ${sepAxis}`); } else { // No separator needed svg.append("line") .attr("x1", px1).attr("y1", py1).attr("x2", px2).attr("y2", py2) .attr("stroke", "#27ae60").attr("stroke-width", 1).attr("stroke-dasharray", "3,3").attr("opacity", 0.5); } } } // Draw singularities for (const s of sings) { const cx = ox + s.x * sc, cy = oy - s.y * sc; svg.append("circle").attr("cx", cx).attr("cy", cy).attr("r", 7).attr("fill", "#e74c3c"); svg.append("text").attr("x", cx).attr("y", cy - 12) .attr("text-anchor", "middle").attr("font-size", 13).attr("fill", "#e74c3c").attr("font-weight", "bold") .text(s.label); } // Info text svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle") .attr("font-size", 13).attr("fill", "#666") .text(`Target edge length: ${sep_scale.toFixed(1)} — ${numSeparators} separator constraint(s) needed`); // Legend svg.append("line").attr("x1", 15).attr("y1", 18).attr("x2", 35).attr("y2", 18) .attr("stroke", "#e67e22").attr("stroke-width", 2.5); svg.append("text").attr("x", 40).attr("y", 22).attr("font-size", 11).attr("fill", "#666").text("Needs separator (too close)"); svg.append("line").attr("x1", 15).attr("y1", 35).attr("x2", 35).attr("y2", 35) .attr("stroke", "#27ae60").attr("stroke-width", 1).attr("stroke-dasharray", "3,3"); svg.append("text").attr("x", 40).attr("y", 39).attr("font-size", 11).attr("fill", "#666").text("Sufficiently separated"); return svg.node(); } ``` ::: {.callout-tip} ## Why separation helps the solver The singularity separation constraints are **cuts** in the optimization literature — they don't remove any optimal integer solutions, but they tighten the continuous relaxation so it's closer to the best integer solution. This dramatically reduces the search tree that the branch-and-cut solver must explore. In experiments, without separation constraints the solver often fails to find any feasible solution within hours. ::: --- ## 7. The Full Algorithm ### 7.1 The Six-Step Pipeline Based on the described components, the algorithm can be understood as a series of bijective maps between meshes: 1. **Globally smooth parametrization** ($\mathcal{M} \to \mathcal{M}_{\text{relaxed}}$) — Solve the continuous relaxation of the quality metric using IPOPT 2. **Decimation in parametric space** ($\mathcal{M}_{\text{relaxed}} \to \mathcal{M}_-$) — Aggressively simplify flat regions, keeping singularities fixed 3. **Addition of singularity separators** ($\mathcal{M}_- \to \mathcal{M}^{\text{sep}}$) — Add geodesic separation constraints between close singularity pairs 4. **MIQP optimization** ($\mathcal{M}^{\text{sep}} \to \mathcal{M}_{\text{int}}$) — Solve the convex MIQP with branch-and-cut (CPLEX) 5. **(Optional) Refinement** ($\mathcal{M}_{\text{int}} \to \mathcal{M}_{\text{IGM}}$) — Split edges to recover DOFs for smoother iso-lines 6. **Quad mesh extraction** ($\mathcal{M}_{\text{IGM}} \to \mathcal{Q}$) — Contour the integer iso-lines to produce the final quad mesh Each step is visualized in the pipeline step-through below: <div id="pipeline-info" style="font-family: monospace; font-size: 13px; padding: 8px 12px; background: var(--bs-tertiary-bg, #f8f9fa); border-radius: 4px; margin-bottom: 8px;"> Step 1/6: Input mesh M with cross-field and sizing </div> <div style="margin-bottom: 8px; display: flex; flex-wrap: wrap; gap: 4px;"> <button id="pipe-step1" class="btn btn-sm btn-primary" onclick="window.pipelineStep('step1')">1. Parametrize</button> <button id="pipe-step2" class="btn btn-sm btn-outline-secondary" onclick="window.pipelineStep('step2')" disabled>2. Decimate</button> <button id="pipe-step3" class="btn btn-sm btn-outline-secondary" onclick="window.pipelineStep('step3')" disabled>3. Separate</button> <button id="pipe-step4" class="btn btn-sm btn-outline-secondary" onclick="window.pipelineStep('step4')" disabled>4. MIQP</button> <button id="pipe-step5" class="btn btn-sm btn-outline-secondary" onclick="window.pipelineStep('step5')" disabled>5. Refine</button> <button id="pipe-step6" class="btn btn-sm btn-outline-secondary" onclick="window.pipelineStep('step6')" disabled>6. Extract</button> </div> <canvas id="pipeline-canvas" style="width: 100%; max-width: 650px; border: 1px solid var(--bs-border-color, #dee2e6); border-radius: 4px;"></canvas> ```{=html} <script> (function() { const canvas = document.getElementById('pipeline-canvas'); const ctx = canvas.getContext('2d'); const info = document.getElementById('pipeline-info'); function resizeCanvas() { const rect = canvas.parentElement.getBoundingClientRect(); const w = Math.min(rect.width, 650); canvas.style.width = w + 'px'; canvas.style.height = (w * 0.65) + 'px'; canvas.width = w * devicePixelRatio; canvas.height = w * 0.65 * devicePixelRatio; ctx.setTransform(devicePixelRatio, 0, 0, devicePixelRatio, 0, 0); } let currentStep = 0; const steps = ['step1', 'step2', 'step3', 'step4', 'step5', 'step6']; const stepLabels = [ 'M → M_relaxed: Continuous parametrization (IPOPT)', 'M_relaxed → M_-: Decimate flat regions', 'M_- → M^sep: Add singularity separators', 'M^sep → M_int: Solve MIQP (CPLEX branch-and-cut)', 'M_int → M_IGM: Refine mesh (recover DOFs)', 'M_IGM → Q: Extract quad mesh from integer iso-lines' ]; const stepMetrics = [ 'Vertices: ~15000 | Variables: ~30000 | Continuous QP', 'Vertices: ~109 | Variables: ~218 | Decimated mesh', 'Vertices: ~109 | Separators: 3 | Linear constraints added', 'Vertices: ~109 | Integer vars: ~50 | Branch-and-cut', 'Vertices: ~2000 | Only continuous DOFs re-optimized', 'Quad faces extracted | Valid IGM guaranteed' ]; // Generate mesh data for each stage const N = 8; function makeGrid(n, perturbation) { const pts = []; for (let j = 0; j <= n; j++) { for (let i = 0; i <= n; i++) { let x = i / n, y = j / n; if (perturbation && i > 0 && i < n && j > 0 && j < n) { x += (Math.random() - 0.5) * perturbation; y += (Math.random() - 0.5) * perturbation; } pts.push([x, y]); } } const tris = []; for (let j = 0; j < n; j++) { for (let i = 0; i < n; i++) { const v0 = j*(n+1)+i, v1 = v0+1, v2 = v0+(n+1), v3 = v2+1; tris.push([v0, v1, v3]); tris.push([v0, v3, v2]); } } return {pts, tris}; } // Seed random for reproducibility let seed = 42; function seededRandom() { seed = (seed * 16807) % 2147483647; return (seed - 1) / 2147483646; } const fullMesh = makeGrid(N, 0); const singVerts = [20, 52]; // two singularities function draw() { const cw = canvas.width / devicePixelRatio; const ch = canvas.height / devicePixelRatio; const bg = getComputedStyle(document.documentElement).getPropertyValue('--bs-body-bg') || '#fff'; const fg = getComputedStyle(document.documentElement).getPropertyValue('--bs-body-color') || '#333'; ctx.fillStyle = bg.trim() || '#fff'; ctx.fillRect(0, 0, cw, ch); const margin = 50; const drawW = cw - 2*margin, drawH = ch - 2*margin - 20; function toScreen(u, v) { return [margin + u * drawW, ch - margin - v * drawH]; } // Integer grid (light) ctx.strokeStyle = '#e8e8e8'; ctx.lineWidth = 0.5; for (let i = 0; i <= N; i++) { const f = i / N; const [x1, y1] = toScreen(f, 0); const [x2, y2] = toScreen(f, 1); ctx.beginPath(); ctx.moveTo(x1, y1); ctx.lineTo(x2, y2); ctx.stroke(); const [x3, y3] = toScreen(0, f); const [x4, y4] = toScreen(1, f); ctx.beginPath(); ctx.moveTo(x3, y3); ctx.lineTo(x4, y4); ctx.stroke(); } const mesh = fullMesh; const showCoarse = currentStep >= 1 && currentStep <= 4; const showQuad = currentStep >= 5; // Draw mesh edges ctx.strokeStyle = showQuad ? 'steelblue' : 'rgba(150,150,150,0.5)'; ctx.lineWidth = showQuad ? 1.5 : 0.8; if (showQuad) { // Draw quad grid lines ctx.strokeStyle = 'steelblue'; ctx.lineWidth = 1.5; for (let i = 0; i <= N; i++) { const f = i / N; const [x1, y1] = toScreen(f, 0); const [x2, y2] = toScreen(f, 1); ctx.beginPath(); ctx.moveTo(x1, y1); ctx.lineTo(x2, y2); ctx.stroke(); const [x3, y3] = toScreen(0, f); const [x4, y4] = toScreen(1, f); ctx.beginPath(); ctx.moveTo(x3, y3); ctx.lineTo(x4, y4); ctx.stroke(); } } else { for (const tri of mesh.tris) { const [x0, y0] = toScreen(mesh.pts[tri[0]][0], mesh.pts[tri[0]][1]); const [x1, y1] = toScreen(mesh.pts[tri[1]][0], mesh.pts[tri[1]][1]); const [x2, y2] = toScreen(mesh.pts[tri[2]][0], mesh.pts[tri[2]][1]); ctx.beginPath(); ctx.moveTo(x0, y0); ctx.lineTo(x1, y1); ctx.lineTo(x2, y2); ctx.closePath(); ctx.stroke(); } } // Separators (step 3+) if (currentStep >= 2 && currentStep < 5) { ctx.strokeStyle = '#e67e22'; ctx.lineWidth = 2.5; const s1 = mesh.pts[singVerts[0]], s2 = mesh.pts[singVerts[1]]; const [sx1, sy1] = toScreen(s1[0], s1[1]); const [sx2, sy2] = toScreen(s2[0], s2[1]); ctx.setLineDash([6, 4]); ctx.beginPath(); ctx.moveTo(sx1, sy1); ctx.lineTo(sx2, sy2); ctx.stroke(); ctx.setLineDash([]); const mx = (sx1+sx2)/2, my = (sy1+sy2)/2; ctx.fillStyle = '#e67e22'; ctx.font = '11px monospace'; ctx.textAlign = 'center'; ctx.fillText('separator', mx, my - 8); } // MIQP highlight (step 4): show integer snapping if (currentStep >= 3 && currentStep < 5) { for (const si of singVerts) { const p = mesh.pts[si]; // Draw arrow from continuous to nearest integer const intX = Math.round(p[0] * N) / N, intY = Math.round(p[1] * N) / N; const [sx, sy] = toScreen(p[0], p[1]); const [ix, iy] = toScreen(intX, intY); ctx.strokeStyle = '#8e44ad'; ctx.lineWidth = 1.5; ctx.setLineDash([3, 3]); ctx.beginPath(); ctx.moveTo(sx, sy); ctx.lineTo(ix, iy); ctx.stroke(); ctx.setLineDash([]); ctx.beginPath(); ctx.arc(ix, iy, 4, 0, 2*Math.PI); ctx.fillStyle = '#8e44ad'; ctx.fill(); } } // Singularity markers for (const si of singVerts) { const [sx, sy] = toScreen(mesh.pts[si][0], mesh.pts[si][1]); ctx.beginPath(); ctx.arc(sx, sy, 5, 0, 2*Math.PI); ctx.fillStyle = '#e74c3c'; ctx.fill(); } // Step label ctx.fillStyle = fg.trim() || '#333'; ctx.font = '11px monospace'; ctx.textAlign = 'left'; ctx.fillText(stepMetrics[currentStep], 10, 16); } const stepBtns = steps.map(s => document.getElementById('pipe-' + s)); window.pipelineStep = function(step) { const idx = steps.indexOf(step); if (idx < 0) return; currentStep = idx; info.textContent = `Step ${idx+1}/${steps.length}: ${stepLabels[idx]}`; for (let i = 0; i <= idx; i++) { stepBtns[i].classList.remove('btn-outline-secondary'); stepBtns[i].classList.add('btn-primary'); stepBtns[i].disabled = false; } if (idx + 1 < steps.length) { stepBtns[idx+1].disabled = false; } draw(); }; resizeCanvas(); window.addEventListener('resize', () => { resizeCanvas(); draw(); }); draw(); })(); </script> ``` ### 7.2 Implementation Details **Optimization of (P2) type problems** is done in steps 1, 3, 4, and 5. Variables in each triangle are merged along a dual spanning tree, and an arbitrary singularity is fixed to $(0,0)$ to obtain a unique minimum. **Features and Boundaries** are preserved by adding linear constraints ensuring feature edges map onto integer-grid lines. The feature graph consists of closed or open chains of mesh edges ending at feature vertices or singularities. **Numerical Solver**: IPOPT for continuous optimization problems, CPLEX for mixed-integer problems. Programs are restricted to quadratic ($k = 2$). **Lazy Constraints**: Orientation constraints (4) are typically violated for only a small fraction of triangles. The algorithm starts with an empty set and iteratively adds violated constraints until a valid solution is found, significantly reducing runtime. ### 7.3 Optional Refinement The refinement of step 5 is the inverse of decimation — it splits edges in $\mathcal{M}_{\text{int}}$ at their midpoints (prioritized by decreasing edge length) to add DOFs that smooth out the iso-lines. Only continuous DOFs are re-optimized; the integer $\mathbf{t}_{ij}$ values remain fixed. This step terminates when all edges are shorter than a user-prescribed tolerance (default: 1). --- ## 8. Results ### 8.1 Importance of Decimation and Singularity Separation > *"To be able to find any solution without our optimizations, we reduced the mesh from 15k to 3k vertices (by keeping the singularities). In the first test we used our decimation strategy without addition of singularity separators and not a single solution could be found within 8h."* The experiments on the BOTIJO model demonstrate that **both** optimizations are essential: | Configuration | Result | |:---|:---| | No decimation, no separation | No solution in 8+ hours | | Decimation only, no separation | No solution in 8+ hours | | Separation only, no decimation | First low-quality solution after 20 min | | **Both decimation + separation** | **High-quality solution in 17 seconds** | ### 8.2 Comparison with Greedy Rounding Methods The paper compares RMIQ (their method) against QuadCover (QC), Mixed-Integer Quadrangulation (MIQ), and MIQ with stiffening on the ROCKERARM model at various target edge lengths $h$: | Algorithm | $h = 10^{-4}$ | $h = 0.025$ | $h = 0.05$ | $h = 0.2$ | |:---|:---:|:---:|:---:|:---:| | QuadCover | 8 degen | 51 degen | 82 degen | 2047 degen | | MIQ | 8 degen | 18 degen | 42 degen | 568 degen | | MIQ Stiff | 0 degen | 0 degen | 1 degen | 541 degen | | **RMIQ** | **0 degen** | **0 degen** | **0 degen** | **0 degen** | RMIQ produces **zero degenerate faces** at all scales, including coarse layouts where greedy methods completely fail. ### 8.3 Coarse Quad Layouts For coarse quad layout generation (target edge length = 20% of bounding box), RMIQ is compared against Dual Loops meshing: | Model | Singularities | Patches (Dual Loops) | Patches (RMIQ) | Patches (Enhanced CF) | |:---|:---:|:---:|:---:|:---:| | BOTIJO | 72 | 221 | 175 | 112 | | BLOCK | 48 | 76 | 76 | — | | GUY | 40 | 168 | 55 | — | | ELK | 52 | 86 | 62 | 58 | | ROCKERARM | 30 | 115 | 74 | 66 | | FERTILITY | 48 | 98 | 117 | 85 | The global optimization approach produces quad layouts with **fewer patches** compared to the greedy Dual Loops method, while maintaining alignment with the cross field. Slight improvements to the cross field (manually repositioning singularities) can further reduce patch counts. --- ## 9. Connection to Our Work This paper addresses the **parameterization stage** of the quad meshing pipeline — the step that converts a cross field and sizing field into a global parameterization from which the quad mesh is extracted. In our MetricFrameField pipeline, this corresponds to: | Pipeline Stage | Bommes et al. 2013 | MetricFrameField | |:---|:---|:---| | Input | Cross field + sizing field on triangle mesh | Cross field from metric-driven solve | | Continuous relaxation | IPOPT solve of $E_\alpha^k$ | MIQ-style parameterization | | Integer constraints | MIQP with branch-and-cut (CPLEX) | Greedy rounding (MIQ) | | Orientation guarantee | Convex trisector constraints | Not yet implemented | | Complexity reduction | Parametric-space decimation | Not yet implemented | | Singularity separation | Geodesic disc constraints | Not yet implemented | | Output | Guaranteed valid IGM | Parameterization (may have degeneracies) | The key takeaway: the RMIQ approach provides a principled replacement for greedy rounding that **guarantees** a valid quad mesh. The trisector linearization (Sec. 4) and complexity reduction (Sec. 5) are the critical ingredients that make this practical. Integrating these ideas into our pipeline would eliminate the degeneracy failures that currently occur with coarse target sizing.

Integer-Grid Maps for Reliable Quad Meshing

Guaranteed Valid Parameterization via Mixed-Integer Optimization

Paper at a Glance

Reference

Integer-Grid Maps for Reliable Quad Meshing — Bommes, Campen, Ebke, Alliez, Kobbelt (2013)

ACM Transactions on Graphics (SIGGRAPH), Vol. 32, No. 4, Article 98.

Download the original paper (PDF)

Parametrization-based quad meshing produces high-quality meshes but suffers from degeneracies — the Jacobian determinant of the piecewise linear map can go to zero or become negative, creating flipped triangles, collapsed edges, and local non-injectivities that prevent extracting a valid quad mesh. This paper defines the class of Integer-Grid Maps (IGMs) that are guaranteed to produce valid quad meshes, and derives a practical Mixed-Integer Quadratic Program (MIQP) formulation to find them. Two key optimizations — complexity reduction via parametric-space decimation and singularity separation via geodesic constraints — make the formulation tractable for real-world meshes.

1. Introduction

1.1 The Reliability Problem

“The main source of non-reliable behavior are degeneracies in the constructed parametrization function, i.e. the Jacobian determinant of the piecewise linear map is locally less than or equal to zero.”

The standard quad meshing pipeline works by mapping a triangle mesh onto a 2D parametric domain, then extracting the quad mesh from the integer iso-lines of that map. The mapping is piecewise linear — one affine map per triangle. When the Jacobian determinant of any triangle’s map goes to zero or becomes negative, the triangle either collapses to a line or flips its orientation. These degeneracies destroy the essential property that the mapped integer grid forms a valid quad mesh.

“Geometrically a value less than zero characterizes triangles that flip their orientation in the domain while a zero value results from local non-injectivities where an edge or even a whole triangle is mapped to a single point.”

Existing approaches handle this with heuristics — greedy rounding of continuous solutions and stiffening of high-distortion regions — but these fail for coarse quad layouts where the target edge length is comparable to the distance between singularities.

Drag the triangle vertices to see how the Jacobian determinant changes. When the triangle flips (vertices go clockwise), the determinant becomes negative:

Code

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

  // Mutable vertex state stored on the SVG node itself so it persists across redraws
  const node = svg.node();
  if (!node.__verts) {
    node.__verts = [[150, 320], [350, 320], [250, 120]];
  }
  const verts = node.__verts;

  function computeMathDet(v) {
    const ux = v[1][0] - v[0][0], uy = v[1][1] - v[0][1];
    const vx = v[2][0] - v[0][0], vy = v[2][1] - v[0][1];
    // Negate because SVG y-axis is flipped vs. math convention
    return -(ux * vy - uy * vx);
  }

  function redraw() {
    svg.selectAll("*").remove();

    // Integer grid
    for (let x = 50; x <= 450; x += 50)
      svg.append("line").attr("x1", x).attr("y1", 20).attr("x2", x).attr("y2", 400)
        .attr("stroke", "#eee").attr("stroke-width", 0.5);
    for (let y = 20; y <= 400; y += 50)
      svg.append("line").attr("x1", 50).attr("y1", y).attr("x2", 450).attr("y2", y)
        .attr("stroke", "#eee").attr("stroke-width", 0.5);

    const mathDet = computeMathDet(verts);

    // Triangle fill color
    const fillColor = mathDet > 500 ? "rgba(39, 174, 96, 0.25)" :
                      mathDet > 100 ? "rgba(39, 174, 96, 0.15)" :
                      mathDet > 0   ? "rgba(243, 156, 18, 0.25)" :
                      mathDet === 0 ? "rgba(149, 165, 166, 0.3)" :
                                      "rgba(231, 76, 60, 0.25)";
    const strokeColor = mathDet > 0 ? "#27ae60" : mathDet === 0 ? "#95a5a6" : "#e74c3c";

    svg.append("polygon")
      .attr("points", verts.map(v => v.join(",")).join(" "))
      .attr("fill", fillColor).attr("stroke", strokeColor).attr("stroke-width", 2.5);

    // Dashed edge vectors from v0
    svg.append("line").attr("x1", verts[0][0]).attr("y1", verts[0][1])
      .attr("x2", verts[1][0]).attr("y2", verts[1][1])
      .attr("stroke", "steelblue").attr("stroke-width", 1.5).attr("stroke-dasharray", "4,2");
    svg.append("line").attr("x1", verts[0][0]).attr("y1", verts[0][1])
      .attr("x2", verts[2][0]).attr("y2", verts[2][1])
      .attr("stroke", "#e67e22").attr("stroke-width", 1.5).attr("stroke-dasharray", "4,2");

    // Draggable vertices
    const labels = ["u\u2080", "u\u2081", "u\u2082"];
    const colors = ["#2c3e50", "steelblue", "#e67e22"];
    for (let i = 0; i < 3; i++) {
      svg.append("circle").attr("cx", verts[i][0]).attr("cy", verts[i][1]).attr("r", 10)
        .attr("fill", colors[i]).attr("cursor", "grab").attr("data-idx", i)
        .call(d3.drag()
          .on("start", function() { d3.select(this).attr("cursor", "grabbing"); })
          .on("drag", function(event) {
            const idx = +d3.select(this).attr("data-idx");
            verts[idx][0] = Math.max(20, Math.min(w - 20, event.x));
            verts[idx][1] = Math.max(20, Math.min(h - 30, event.y));
            redraw();
          })
          .on("end", function() { d3.select(this).attr("cursor", "grab"); })
        );
      svg.append("text").attr("x", verts[i][0]).attr("y", verts[i][1] - 15)
        .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", colors[i])
        .attr("font-weight", "bold").attr("pointer-events", "none")
        .text(labels[i]);
    }

    // Determinant display
    const statusText = mathDet > 500 ? "Valid (well-shaped)" :
                       mathDet > 100 ? "Valid" :
                       mathDet > 0   ? "Near-degenerate" :
                       mathDet === 0 ? "DEGENERATE (collapsed)" :
                                       "FLIPPED (negative orientation)";
    const statusColor = mathDet > 0 ? "#27ae60" : mathDet === 0 ? "#95a5a6" : "#e74c3c";
    svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle")
      .attr("font-size", 15).attr("fill", statusColor).attr("font-weight", "bold")
      .text(`det = ${mathDet.toFixed(0)}  \u2014  ${statusText}`);
  }

  redraw();
  return node;
}

Figure 1

Why the Jacobian matters

For a piecewise linear map $f$, each triangle has a constant $2 \times 2$ Jacobian $J_t = \nabla f_t$. The determinant $\det(J_t)$ measures the signed area ratio between the mapped triangle and the original. When $\det(J_t) > 0$, the triangle preserves orientation. When $\det(J_t) \leq 0$, the triangle has collapsed or flipped — the integer iso-lines cross or overlap, making quad extraction impossible.

1.2 Contributions

The paper makes three main contributions:

A MIQP formulation where the minimizer is guaranteed to be an Integer-Grid Map (Sec. 3.1)
A complexity reduction algorithm that accurately approximates the high-dimensional problem with a low-dimensional one (Sec. 3.2)
Singularity separating conditions that significantly improve the continuous relaxation of the MIQP (Sec. 3.3)

2. Integer-Grid Maps

“The main principle of parametrization based quad meshing algorithms is the mapping of the canonical quad mesh formed by the 2D Cartesian grid of integer iso-lines onto a surface embedded in 3D.”

2.1 Piecewise Linear Maps

Given a triangle mesh $\mathcal{M} = (V, E, T)$, a map $f$ is the union of per-triangle affine maps:

\[f_i : (u_i, v_i, w_i) \in \mathbb{R}^{2 \times 3} \mapsto (p_i, q_i, r_i) \in \mathbb{R}^{3 \times 3}\]

Each triangle $t_i$ in 3D is mapped to a triangle in the 2D parametric domain. The intersection of the integer grid $\mathbb{Z} \times \mathbb{Z}$ with the image of $f$ defines the quad mesh edges.

2.2 Transition Functions

Since each triangle has its own chart, we need transition functions to describe how neighboring charts relate. The transition function from the chart of triangle $t_i$ to that of neighboring triangle $t_j$ is:

\[g_{i \to j}(\mathbf{a}) = R_{90}^{r_{ij}} \mathbf{a} + \mathbf{t}_{ij} \tag{1}\]

where $r_{ij} \in \{0, 1, 2, 3\}$ is a rotation by multiples of $\pi/2$ and $\mathbf{t}_{ij} \in \mathbb{Z}^2$ is an integer translation.

This means neighboring charts can differ by a 90-degree rotation and an integer shift — both of which leave the Cartesian integer grid invariant. This is exactly the definition of an integer-grid automorphism.

Adjust the rotation $r_{ij}$ and translation $\mathbf{t}_{ij}$ to see how the transition function relates two neighboring charts:

Code

viewof r_ij = Inputs.range([0, 3], {value: 0, step: 1, label: "Rotation r_ij (× 90°)"})

Figure 2

Code

viewof t_ij_x = Inputs.range([-2, 2], {value: 0, step: 1, label: "Translation t_x"})

Figure 3

Code

viewof t_ij_y = Inputs.range([-2, 2], {value: 0, step: 1, label: "Translation t_y"})

Figure 4

Code

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

  const leftCx = 170, rightCx = 530, cy = 200;
  const scale = 50;

  // Draw grids
  function drawGrid(gx, gy, label) {
    for (let i = -3; i <= 3; i++) {
      svg.append("line")
        .attr("x1", gx + i*scale).attr("y1", gy - 3*scale)
        .attr("x2", gx + i*scale).attr("y2", gy + 3*scale)
        .attr("stroke", "#eee").attr("stroke-width", i === 0 ? 1 : 0.5);
      svg.append("line")
        .attr("x1", gx - 3*scale).attr("y1", gy + i*scale)
        .attr("x2", gx + 3*scale).attr("y2", gy + i*scale)
        .attr("stroke", "#eee").attr("stroke-width", i === 0 ? 1 : 0.5);
    }
    // Integer points
    for (let ix = -2; ix <= 2; ix++) {
      for (let iy = -2; iy <= 2; iy++) {
        svg.append("circle")
          .attr("cx", gx + ix*scale).attr("cy", gy - iy*scale).attr("r", 2)
          .attr("fill", "#ccc");
      }
    }
    svg.append("text").attr("x", gx).attr("y", gy + 3.5*scale)
      .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", "#333").attr("font-weight", "bold")
      .text(label);
  }

  drawGrid(leftCx, cy, "Chart of t_i");
  drawGrid(rightCx, cy, "Chart of t_j");

  // Triangle in chart i
  const triI = [[0.3, 0.5], [1.8, 0.2], [1.0, 1.5]];
  const triIPts = triI.map(p => [leftCx + p[0]*scale, cy - p[1]*scale]);
  svg.append("polygon")
    .attr("points", triIPts.map(p => p.join(",")).join(" "))
    .attr("fill", "rgba(70,130,180,0.15)").attr("stroke", "steelblue").attr("stroke-width", 2);

  // Shared edge (highlighted)
  svg.append("line")
    .attr("x1", triIPts[1][0]).attr("y1", triIPts[1][1])
    .attr("x2", triIPts[2][0]).attr("y2", triIPts[2][1])
    .attr("stroke", "#e74c3c").attr("stroke-width", 3);

  // Apply transition function to triangle
  const angle = r_ij * Math.PI / 2;
  const cosA = Math.round(Math.cos(angle)), sinA = Math.round(Math.sin(angle));

  function applyTransition(p) {
    const rx = cosA * p[0] - sinA * p[1] + t_ij_x;
    const ry = sinA * p[0] + cosA * p[1] + t_ij_y;
    return [rx, ry];
  }

  const triJ = triI.map(p => applyTransition(p));
  const triJPts = triJ.map(p => [rightCx + p[0]*scale, cy - p[1]*scale]);
  svg.append("polygon")
    .attr("points", triJPts.map(p => p.join(",")).join(" "))
    .attr("fill", "rgba(230,126,34,0.15)").attr("stroke", "#e67e22").attr("stroke-width", 2);

  // Shared edge in j's chart
  svg.append("line")
    .attr("x1", triJPts[1][0]).attr("y1", triJPts[1][1])
    .attr("x2", triJPts[2][0]).attr("y2", triJPts[2][1])
    .attr("stroke", "#e74c3c").attr("stroke-width", 3);

  // Arrow between charts
  svg.append("line")
    .attr("x1", leftCx + 3.2*scale).attr("y1", cy)
    .attr("x2", rightCx - 3.2*scale).attr("y2", cy)
    .attr("stroke", "#999").attr("stroke-width", 1.5).attr("marker-end", "url(#arrowhead)");

  svg.append("defs").append("marker").attr("id", "arrowhead").attr("viewBox", "0 0 10 10")
    .attr("refX", 10).attr("refY", 5).attr("markerWidth", 8).attr("markerHeight", 8)
    .attr("orient", "auto")
    .append("path").attr("d", "M 0 0 L 10 5 L 0 10 z").attr("fill", "#999");

  svg.append("text")
    .attr("x", (leftCx + rightCx) / 2).attr("y", cy - 15)
    .attr("text-anchor", "middle").attr("font-size", 13).attr("fill", "#666")
    .text(`g: R${r_ij * 90}° + (${t_ij_x}, ${t_ij_y})`);

  // Vertex labels
  for (let i = 0; i < 3; i++) {
    svg.append("circle").attr("cx", triIPts[i][0]).attr("cy", triIPts[i][1]).attr("r", 4).attr("fill", "steelblue");
    svg.append("circle").attr("cx", triJPts[i][0]).attr("cy", triJPts[i][1]).attr("r", 4).attr("fill", "#e67e22");
  }

  return svg.node();
}

Why integer translations?

The integer grid $\mathbb{Z} \times \mathbb{Z}$ is invariant under 90-degree rotations and integer translations. So if chart $i$ and chart $j$ are related by such a transformation, the integer iso-lines in one chart perfectly align with those in the other. This is what makes the quad mesh well-defined across chart boundaries.

2.3 The Three IGM Conditions

The class of Integer-Grid Maps is the subset of all possible maps that additionally stitch the grid of integer iso-lines into a valid quad mesh. The necessary and sufficient conditions are:

Condition 1 — Transition Functions: The transition function $g_{i \to j}$ from the chart of triangle $t_i$ to the chart of neighboring triangle $t_j$ must be an integer-grid automorphism of the form $g_{i \to j}(\mathbf{a}) = R_{90}^{r_{ij}} \mathbf{a} + \mathbf{t}_{ij}$ (Eq. 1).

Condition 2 — Singular Points: All singular vertices must lie on integer locations in the domain:

\[\mathbf{f}^{-1}(s_i) \in \mathbb{Z}^2 \quad \forall s_i \in S \tag{2}\]

This ensures the cone singularities (where the cross field has nonzero index) map to grid intersections, which become the irregular vertices of the quad mesh.

Condition 3 — Consistent Orientation: All domain triangles $(\mathbf{u}, \mathbf{v}, \mathbf{w})$ with $\mathbf{u}, \mathbf{v}, \mathbf{w} \in \mathbb{R}^2$ must have positive orientation:

\[\det[\mathbf{v} - \mathbf{u},\; \mathbf{w} - \mathbf{u}] > 0 \tag{3}\]

This prevents triangle flipping — the degeneracy that makes greedy rounding fail.

“For the remainder of this paper it is very important to keep in mind that an IGM is a bijective map between two piecewise linear 2-manifold meshes.”

Toggle between the three conditions to see what each one enforces on a small mesh in parametric space:

Code

viewof igm_cond = Inputs.radio(["All conditions", "1: Transition functions", "2: Singular vertices on integers", "3: Consistent orientation"], {value: "All conditions", label: "Show condition"})

Figure 5

Code

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

  const ox = 75, oy = 375, sc = 70;

  // Integer grid
  for (let i = 0; i <= 5; i++) {
    svg.append("line").attr("x1", ox + i*sc).attr("y1", oy - 5*sc).attr("x2", ox + i*sc).attr("y2", oy)
      .attr("stroke", "#e8e8e8").attr("stroke-width", 1);
    svg.append("line").attr("x1", ox).attr("y1", oy - i*sc).attr("x2", ox + 5*sc).attr("y2", oy - i*sc)
      .attr("stroke", "#e8e8e8").attr("stroke-width", 1);
    // Integer labels
    svg.append("text").attr("x", ox + i*sc).attr("y", oy + 18)
      .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#aaa").text(i);
    if (i > 0) svg.append("text").attr("x", ox - 15).attr("y", oy - i*sc + 4)
      .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#aaa").text(i);
  }

  // Integer points
  for (let ix = 0; ix <= 5; ix++) {
    for (let iy = 0; iy <= 5; iy++) {
      svg.append("circle")
        .attr("cx", ox + ix*sc).attr("cy", oy - iy*sc).attr("r", 2).attr("fill", "#ddd");
    }
  }

  // Define a small mesh: 6 triangles forming a patch with 1 singularity
  // Vertices in parametric space (u, v)
  const verts = [
    [0.5, 0.5],  // 0
    [2.0, 0.3],  // 1
    [3.5, 0.8],  // 2
    [0.3, 2.0],  // 3
    [2.0, 2.0],  // 4: singularity (on integer!)
    [3.8, 2.2],  // 5
    [0.8, 3.5],  // 6
    [2.2, 3.7],  // 7
    [3.5, 3.3],  // 8
  ];

  const tris = [
    [0, 1, 4], [1, 2, 5], [1, 5, 4],
    [0, 4, 3], [3, 4, 6], [4, 7, 6],
    [4, 5, 8], [4, 8, 7],
  ];

  const edges = [
    [0,1],[1,2],[2,5],[5,8],[8,7],[7,6],[6,3],[3,0],
    [0,4],[1,4],[1,5],[3,4],[4,5],[4,6],[4,7],[4,8]
  ];

  const singIdx = 4; // vertex 4 is the singularity

  const showCond1 = igm_cond === "All conditions" || igm_cond === "1: Transition functions";
  const showCond2 = igm_cond === "All conditions" || igm_cond === "2: Singular vertices on integers";
  const showCond3 = igm_cond === "All conditions" || igm_cond === "3: Consistent orientation";

  // Draw triangles
  for (const tri of tris) {
    const pts = tri.map(i => [ox + verts[i][0]*sc, oy - verts[i][1]*sc]);

    // Condition 3: color by orientation
    let fillC = "rgba(200,200,200,0.08)";
    if (showCond3) {
      const ux = verts[tri[1]][0] - verts[tri[0]][0];
      const uy = verts[tri[1]][1] - verts[tri[0]][1];
      const vx = verts[tri[2]][0] - verts[tri[0]][0];
      const vy = verts[tri[2]][1] - verts[tri[0]][1];
      const det = ux * vy - uy * vx;
      fillC = det > 0 ? "rgba(39, 174, 96, 0.12)" : "rgba(231, 76, 60, 0.2)";
    }

    svg.append("polygon")
      .attr("points", pts.map(p => p.join(",")).join(" "))
      .attr("fill", fillC)
      .attr("stroke", "#bbb").attr("stroke-width", 1);
  }

  // Condition 1: highlight internal edges showing transition compatibility
  if (showCond1) {
    const internalEdges = [[0,4],[1,4],[1,5],[3,4],[4,5],[4,6],[4,7],[4,8]];
    for (const e of internalEdges) {
      const p0 = [ox + verts[e[0]][0]*sc, oy - verts[e[0]][1]*sc];
      const p1 = [ox + verts[e[1]][0]*sc, oy - verts[e[1]][1]*sc];
      svg.append("line")
        .attr("x1", p0[0]).attr("y1", p0[1]).attr("x2", p1[0]).attr("y2", p1[1])
        .attr("stroke", "steelblue").attr("stroke-width", 2.5).attr("opacity", 0.7);
    }
    svg.append("text").attr("x", w - 10).attr("y", 25)
      .attr("text-anchor", "end").attr("font-size", 12).attr("fill", "steelblue")
      .text("Internal edges: transition functions must be IGM automorphisms");
  }

  // Condition 2: highlight singularity
  if (showCond2) {
    const sx = ox + verts[singIdx][0]*sc, sy = oy - verts[singIdx][1]*sc;
    // Check if on integer
    const onInt = Math.abs(verts[singIdx][0] - Math.round(verts[singIdx][0])) < 0.01 &&
                  Math.abs(verts[singIdx][1] - Math.round(verts[singIdx][1])) < 0.01;
    svg.append("circle").attr("cx", sx).attr("cy", sy).attr("r", 14)
      .attr("fill", "none").attr("stroke", onInt ? "#27ae60" : "#e74c3c")
      .attr("stroke-width", 2.5).attr("stroke-dasharray", "4,2");
    svg.append("text").attr("x", sx + 18).attr("y", sy - 8)
      .attr("font-size", 11).attr("fill", onInt ? "#27ae60" : "#e74c3c")
      .text(onInt ? "on Z^2" : "NOT on Z^2");
    svg.append("text").attr("x", w - 10).attr("y", 45)
      .attr("text-anchor", "end").attr("font-size", 12).attr("fill", "#e74c3c")
      .text("Singular vertex must map to integer coordinates");
  }

  // Condition 3: label
  if (showCond3) {
    svg.append("text").attr("x", w - 10).attr("y", 65)
      .attr("text-anchor", "end").attr("font-size", 12).attr("fill", "#27ae60")
      .text("All triangles must have positive orientation (green)");
  }

  // Draw vertices
  for (let i = 0; i < verts.length; i++) {
    const px = ox + verts[i][0]*sc, py = oy - verts[i][1]*sc;
    const isSing = i === singIdx;
    svg.append("circle").attr("cx", px).attr("cy", py)
      .attr("r", isSing ? 6 : 3.5)
      .attr("fill", isSing ? "#e74c3c" : "#2c3e50");
  }

  svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#666")
    .text("A small mesh patch in parametric space with integer grid overlay");

  return svg.node();
}

3. The Quality Metric

3.1 The Energy $E(f)$

In parametrization-based quad remeshing, a variational quality metric $E_q(f)$ penalizes distortion of the resulting quad elements based on the map $f$. The standard form is a Frobenius-norm measure of how far the Jacobian is from the desired sizing and orientation:

\[E(f) = \frac{1}{2} \sum_{t \in T} \mathcal{A}_t \| R_t^T \nabla f_t - H_t \|_{\text{Frobenius}}^2 \tag{E1}\]

where:

$H_t = \text{diag}(w_t^{-1}, h_t^{-1})$ is the locally desired anisotropic sizing
$R_t \in \mathbb{R}^{4 \times 2}$ contains the normalized $u$ and $v$ directions of the cross field
$\mathcal{A}_t$ is the triangle area

This energy simultaneously measures how well the parameterization aligns with the cross-field directions and reproduces the specified element density.

3.2 Generalized Anisotropy Metric

To weight the importance of cross direction versus sizing fidelity, the authors generalize with a higher-order norm:

\[E_\alpha^k(f) = \frac{1}{2} \sum_{t \in T} \mathcal{A}_t \| R_t^T \nabla f_t - H_t \|_\alpha^k \tag{E2}\]

where $k \in 2\mathbb{Z}^+$, $\alpha \in [0,1]$, and the anisotropy norm is:

\[\left\| \begin{pmatrix} a & b \\ c & d \end{pmatrix} \right\|_\alpha^k = 2\alpha(a^k + d^k) + 2(1 - \alpha)(b^k + c^k)\]

For $\alpha < 0.5$, stretching along the cross-field directions is preferred over angular deviation — useful for coarse quad layouts with cross-field-aligned rectangular patches.

3.3 The Naive MINLP (Problem P1)

Naively searching for the best IGM leads to a Mixed-Integer Nonlinear Program:

\[\text{minimize } E_q(f) \quad \text{s.t.} \quad (1),\; (2),\; (3) \tag{P1}\]

This has $6|T| + 3|E|$ unknowns including at least $3|E|$ discrete variables. Due to:

1. containing nonlinear constraints in $r_{ij}$ and linear in $\mathbf{t}_{ij}$
1. imposing $|T|$ non-convex quadratic inequality constraints
The resulting number of unknowns being $O(|V|)$

“Unfortunately MINLP problems are very hard to solve since they imply all difficulties from continuous as well as discrete optimization.”

Modern mixed-integer solvers need three properties: (i) convex objective and constraints, (ii) low-dimensional search space, (iii) continuous relaxation close to the optimum. The next three sections address each.

4. Convex Consistent Orientation

“The feasible set of the consistent orientation condition (3) is a non-convex region $\mathcal{N} \subset \mathbb{R}^6$.”

4.1 The Non-Convex Problem

Condition (3) requires $\det[\mathbf{v}-\mathbf{u}, \mathbf{w}-\mathbf{u}] > 0$ for every triangle $(\mathbf{u}, \mathbf{v}, \mathbf{w})$. This determinant is a bilinear (quadratic) function of the vertex positions — the set of positions satisfying it is non-convex. We need a convex inner approximation $\mathcal{C} \subset \mathcal{N}$ such that every point in $\mathcal{C}$ satisfies (3).

4.2 The Trisector Construction

The key geometric insight: partition $\mathbb{R}^2$ into 3 sectors using ccw-ordered rays $\mathbf{s}_0, \mathbf{s}_1, \mathbf{s}_2$ emanating from a center point $\mathbf{m}$:

\[S_i = \{ \mathbf{u} \in \mathbb{R}^2 : (\mathbf{u} - \mathbf{m}) \cdot \mathbf{s}_i^\perp > 0 \;\wedge\; (\mathbf{u} - \mathbf{m}) \cdot \mathbf{s}_{i+1}^\perp < 0 \}\]

If vertex $\mathbf{u}_i$ stays strictly within sector $S_i$ and $\mathbf{m}$ is an interior point of the triangle, then the triangle cannot flip. The area is guaranteed positive because it decomposes into three sub-triangles $T_i = (\mathbf{u}_i, \mathbf{u}_{i+1}, \mathbf{m})$, each with positive area.

4.3 The Fermat Point

“To enable maximal rotational freedom we choose a trisector origin $\mathbf{m}$ such that angles between $(\hat{\mathbf{u}}_i - \mathbf{m})$ and $(\hat{\mathbf{u}}_{i+1} - \mathbf{m})$ are $120°$ each.”

The Fermat point (first Torricelli point) is the point inside a triangle where each pair of vertices subtends an angle of exactly $120°$. By placing the trisector origin at the Fermat point, the three sectors are equally sized ($120°$ each), maximizing the space of feasible deformations.

The Fermat point is constructed by attaching equilateral triangles to each edge and intersecting the three connecting diagonals.

Drag the triangle vertices to explore the Fermat point and trisector sectors. As long as each vertex stays in its colored sector, the triangle orientation is guaranteed positive:

Code

viewof tri_preset = Inputs.radio(["Equilateral", "Right triangle", "Obtuse (> 100°)"], {value: "Equilateral", label: "Triangle preset"})

Figure 6

Code

{
  const totalW = 900, panelH = 420, gap = 20;
  const leftW = 420, rightW = totalW - leftW - gap;
  const svg = d3.create("svg")
    .attr("viewBox", `0 0 ${totalW} ${panelH + 30}`)
    .style("width", "100%").style("max-width", totalW + "px");

  let triV;
  if (tri_preset === "Equilateral") {
    triV = [[130, 360], [330, 360], [230, 186]];
  } else if (tri_preset === "Right triangle") {
    triV = [[100, 360], [360, 360], [100, 130]];
  } else {
    triV = [[80, 320], [380, 320], [270, 230]];
  }

  // Shared Fermat point computation
  function fermatPoint(A, B, C) {
    function angle(p1, p2, p3) {
      const v1x = p1[0]-p2[0], v1y = p1[1]-p2[1];
      const v2x = p3[0]-p2[0], v2y = p3[1]-p2[1];
      const dot = v1x*v2x + v1y*v2y;
      const m1 = Math.sqrt(v1x*v1x + v1y*v1y);
      const m2 = Math.sqrt(v2x*v2x + v2y*v2y);
      return Math.acos(Math.max(-1, Math.min(1, dot/(m1*m2))));
    }
    const angA = angle(B, A, C), angB = angle(A, B, C), angC = angle(A, C, B);
    if (angA >= 2*Math.PI/3) return [...A];
    if (angB >= 2*Math.PI/3) return [...B];
    if (angC >= 2*Math.PI/3) return [...C];
    let mx = (A[0]+B[0]+C[0])/3, my = (A[1]+B[1]+C[1])/3;
    for (let iter = 0; iter < 200; iter++) {
      const dA = Math.sqrt((mx-A[0])**2 + (my-A[1])**2) || 1e-10;
      const dB = Math.sqrt((mx-B[0])**2 + (my-B[1])**2) || 1e-10;
      const dC = Math.sqrt((mx-C[0])**2 + (my-C[1])**2) || 1e-10;
      const wA = 1/dA, wB = 1/dB, wC = 1/dC;
      const wSum = wA + wB + wC;
      mx = (wA*A[0] + wB*B[0] + wC*C[0]) / wSum;
      my = (wA*A[1] + wB*B[1] + wC*C[1]) / wSum;
    }
    return [mx, my];
  }

  const fp = fermatPoint(triV[0], triV[1], triV[2]);
  const sectorColors = ["rgba(231,76,60,0.1)", "rgba(39,174,96,0.1)", "rgba(52,152,219,0.1)"];
  const sectorStroke = ["#e74c3c", "#27ae60", "#3498db"];
  const vLabels = ["u\u2080", "u\u2081", "u\u2082"];

  // --- Helper: draw trisector scene into a clipping group ---
  function drawScene(g, clipX, clipY, clipW, clipH, showZoomBox) {
    // Clip to panel bounds
    const clipId = `clip-${clipX}`;
    const defs = g.append("defs");
    defs.append("clipPath").attr("id", clipId)
      .append("rect").attr("x", clipX).attr("y", clipY).attr("width", clipW).attr("height", clipH);
    const scene = g.append("g").attr("clip-path", `url(#${clipId})`);

    // Panel background
    scene.append("rect").attr("x", clipX).attr("y", clipY)
      .attr("width", clipW).attr("height", clipH)
      .attr("fill", "#fafafa").attr("stroke", "#ccc").attr("stroke-width", 1);

    // Sector rays and fills
    const rayLen = 800;
    for (let i = 0; i < 3; i++) {
      const dx = triV[i][0] - fp[0], dy = triV[i][1] - fp[1];
      const len = Math.sqrt(dx*dx + dy*dy) || 1;
      const endX = fp[0] + dx/len * rayLen;
      const endY = fp[1] + dy/len * rayLen;
      scene.append("line")
        .attr("x1", fp[0]).attr("y1", fp[1]).attr("x2", endX).attr("y2", endY)
        .attr("stroke", sectorStroke[i]).attr("stroke-width", 1.5)
        .attr("stroke-dasharray", "6,3").attr("opacity", 0.6);
    }
    for (let i = 0; i < 3; i++) {
      const j = (i + 1) % 3;
      const dx1 = triV[i][0] - fp[0], dy1 = triV[i][1] - fp[1];
      const dx2 = triV[j][0] - fp[0], dy2 = triV[j][1] - fp[1];
      const len1 = Math.sqrt(dx1*dx1+dy1*dy1) || 1;
      const len2 = Math.sqrt(dx2*dx2+dy2*dy2) || 1;
      const e1x = fp[0] + dx1/len1*rayLen, e1y = fp[1] + dy1/len1*rayLen;
      const e2x = fp[0] + dx2/len2*rayLen, e2y = fp[1] + dy2/len2*rayLen;
      scene.append("polygon")
        .attr("points", `${fp[0]},${fp[1]} ${e1x},${e1y} ${e2x},${e2y}`)
        .attr("fill", sectorColors[i]).attr("opacity", 0.5);
    }

    // Triangle
    scene.append("polygon")
      .attr("points", triV.map(v => v.join(",")).join(" "))
      .attr("fill", "rgba(255,255,255,0.6)").attr("stroke", "#2c3e50").attr("stroke-width", 2.5);

    // Fermat point
    scene.append("circle").attr("cx", fp[0]).attr("cy", fp[1]).attr("r", 6).attr("fill", "#8e44ad");
    scene.append("text").attr("x", fp[0] + 12).attr("y", fp[1] - 8)
      .attr("font-size", 14).attr("fill", "#8e44ad").attr("font-weight", "bold").text("m");

    // Vertices
    for (let i = 0; i < 3; i++) {
      scene.append("circle").attr("cx", triV[i][0]).attr("cy", triV[i][1]).attr("r", 7)
        .attr("fill", sectorStroke[i]);
      scene.append("text").attr("x", triV[i][0]).attr("y", triV[i][1] - 14)
        .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", sectorStroke[i]).attr("font-weight", "bold")
        .text(vLabels[i]);
    }

    // Angle labels at Fermat point
    for (let i = 0; i < 3; i++) {
      const j = (i + 1) % 3;
      const a1 = Math.atan2(triV[i][1] - fp[1], triV[i][0] - fp[0]);
      const a2 = Math.atan2(triV[j][1] - fp[1], triV[j][0] - fp[0]);
      let da = a2 - a1;
      while (da < -Math.PI) da += 2*Math.PI;
      while (da > Math.PI) da -= 2*Math.PI;
      const angleDeg = Math.abs(da) * 180 / Math.PI;
      const arcR = 30;
      const mid = a1 + da/2;
      scene.append("text")
        .attr("x", fp[0] + (arcR + 15) * Math.cos(mid))
        .attr("y", fp[1] + (arcR + 15) * Math.sin(mid))
        .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#8e44ad")
        .text(`${angleDeg.toFixed(0)}\u00B0`);
    }

    // Zoom box indicator on the overview panel
    if (showZoomBox) {
      const zoomR = 90;
      scene.append("rect")
        .attr("x", fp[0] - zoomR).attr("y", fp[1] - zoomR)
        .attr("width", 2*zoomR).attr("height", 2*zoomR)
        .attr("fill", "none").attr("stroke", "#8e44ad")
        .attr("stroke-width", 2).attr("stroke-dasharray", "4,3");
    }
  }

  // --- Left panel: zoomed-out overview ---
  const leftG = svg.append("g");
  drawScene(leftG, 0, 25, leftW, panelH, true);
  svg.append("text").attr("x", leftW/2).attr("y", 16).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#555").text("Full triangle with trisector");

  // --- Right panel: zoomed-in around Fermat point ---
  const rx = leftW + gap;
  const zoomR = 90; // radius in scene coords to show around Fermat point
  const zoomScale = Math.min(rightW, panelH) / (2 * zoomR);
  const rightG = svg.append("g");

  // Background
  rightG.append("rect").attr("x", rx).attr("y", 25).attr("width", rightW).attr("height", panelH)
    .attr("fill", "#fafafa").attr("stroke", "#ccc").attr("stroke-width", 1);

  // Clip for right panel
  const clipR = "clip-right";
  rightG.append("defs").append("clipPath").attr("id", clipR)
    .append("rect").attr("x", rx).attr("y", 25).attr("width", rightW).attr("height", panelH);

  // Transform: center Fermat point in the right panel, scale up
  const rcx = rx + rightW/2, rcy = 25 + panelH/2;
  const zoomedG = rightG.append("g").attr("clip-path", `url(#${clipR})`)
    .append("g")
    .attr("transform", `translate(${rcx}, ${rcy}) scale(${zoomScale}) translate(${-fp[0]}, ${-fp[1]})`);

  // Redraw the scene geometry into the zoomed group (without clipping, the outer clip handles it)
  const zRayLen = 800;
  for (let i = 0; i < 3; i++) {
    const dx = triV[i][0] - fp[0], dy = triV[i][1] - fp[1];
    const len = Math.sqrt(dx*dx + dy*dy) || 1;
    const endX = fp[0] + dx/len * zRayLen;
    const endY = fp[1] + dy/len * zRayLen;
    zoomedG.append("line")
      .attr("x1", fp[0]).attr("y1", fp[1]).attr("x2", endX).attr("y2", endY)
      .attr("stroke", sectorStroke[i]).attr("stroke-width", 1.5 / zoomScale)
      .attr("stroke-dasharray", `${6/zoomScale},${3/zoomScale}`).attr("opacity", 0.6);
  }
  for (let i = 0; i < 3; i++) {
    const j = (i + 1) % 3;
    const dx1 = triV[i][0] - fp[0], dy1 = triV[i][1] - fp[1];
    const dx2 = triV[j][0] - fp[0], dy2 = triV[j][1] - fp[1];
    const len1 = Math.sqrt(dx1*dx1+dy1*dy1) || 1;
    const len2 = Math.sqrt(dx2*dx2+dy2*dy2) || 1;
    const e1x = fp[0] + dx1/len1*zRayLen, e1y = fp[1] + dy1/len1*zRayLen;
    const e2x = fp[0] + dx2/len2*zRayLen, e2y = fp[1] + dy2/len2*zRayLen;
    zoomedG.append("polygon")
      .attr("points", `${fp[0]},${fp[1]} ${e1x},${e1y} ${e2x},${e2y}`)
      .attr("fill", sectorColors[i]).attr("opacity", 0.5);
  }
  // Triangle in zoomed view
  zoomedG.append("polygon")
    .attr("points", triV.map(v => v.join(",")).join(" "))
    .attr("fill", "rgba(255,255,255,0.6)").attr("stroke", "#2c3e50")
    .attr("stroke-width", 2.5 / zoomScale);
  // Fermat point
  zoomedG.append("circle").attr("cx", fp[0]).attr("cy", fp[1])
    .attr("r", 6 / zoomScale).attr("fill", "#8e44ad");
  zoomedG.append("text").attr("x", fp[0] + 12 / zoomScale).attr("y", fp[1] - 8 / zoomScale)
    .attr("font-size", 14 / zoomScale).attr("fill", "#8e44ad").attr("font-weight", "bold").text("m");
  // Vertices in zoomed view
  for (let i = 0; i < 3; i++) {
    zoomedG.append("circle").attr("cx", triV[i][0]).attr("cy", triV[i][1])
      .attr("r", 7 / zoomScale).attr("fill", sectorStroke[i]);
    zoomedG.append("text").attr("x", triV[i][0]).attr("y", triV[i][1] - 14 / zoomScale)
      .attr("text-anchor", "middle").attr("font-size", 14 / zoomScale)
      .attr("fill", sectorStroke[i]).attr("font-weight", "bold")
      .text(vLabels[i]);
  }
  // Angle labels in zoomed view
  for (let i = 0; i < 3; i++) {
    const j = (i + 1) % 3;
    const a1 = Math.atan2(triV[i][1] - fp[1], triV[i][0] - fp[0]);
    const a2 = Math.atan2(triV[j][1] - fp[1], triV[j][0] - fp[0]);
    let da = a2 - a1;
    while (da < -Math.PI) da += 2*Math.PI;
    while (da > Math.PI) da -= 2*Math.PI;
    const angleDeg = Math.abs(da) * 180 / Math.PI;
    const arcR = 30 / zoomScale;
    const mid = a1 + da/2;
    zoomedG.append("text")
      .attr("x", fp[0] + (arcR + 15 / zoomScale) * Math.cos(mid))
      .attr("y", fp[1] + (arcR + 15 / zoomScale) * Math.sin(mid))
      .attr("text-anchor", "middle").attr("font-size", 10 / zoomScale).attr("fill", "#8e44ad")
      .text(`${angleDeg.toFixed(0)}\u00B0`);
  }

  svg.append("text").attr("x", rx + rightW/2).attr("y", 16).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#555").text("Zoomed in: Fermat point detail");

  return svg.node();
}

The power of linearization

The trisector construction replaces one non-convex quadratic constraint ($\det > 0$) with six linear inequalities (each vertex must stay in its sector). This is conservative — some valid configurations are excluded — but it guarantees that every solution within the convex feasible region is a valid, non-degenerate triangle.

4.4 Virtual Splitting for Obtuse Triangles

The Fermat point only exists inside the triangle when all angles are less than $120°$. For triangles with an angle larger than $100°$ (a conservative threshold), the paper “virtually” splits the triangle along the altitude at its largest angle:

Code

viewof obtuse_angle = Inputs.range([60, 150], {value: 90, step: 1, label: "Largest angle at c (degrees)"})

Figure 7

Code

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

  // --- Construct triangle with exact angle at C ---
  const baseLen = 180;
  const baseY = 340;
  const angleRad = obtuse_angle * Math.PI / 180;

  // Place A and B symmetrically. C is above with the specified angle ACB.
  // Using the inscribed angle: height from midpoint of AB to C is
  // h = (baseLen/2) / tan(angle/2)
  const halfBase = baseLen / 2;
  const triH = halfBase / Math.tan(angleRad / 2);

  // Left panel: original triangle with Fermat point and sectors
  const leftW = w * 0.48;
  const leftCx = leftW / 2 + 15;

  const A = [leftCx - halfBase, baseY];
  const B = [leftCx + halfBase, baseY];
  const C = [leftCx, baseY - Math.max(25, triH)];

  const needsSplit = obtuse_angle > 100;

  // Compute actual angles at all three vertices
  function vecAngle(p1, vertex, p2) {
    const v1x = p1[0]-vertex[0], v1y = p1[1]-vertex[1];
    const v2x = p2[0]-vertex[0], v2y = p2[1]-vertex[1];
    const dot = v1x*v2x + v1y*v2y;
    const m1 = Math.sqrt(v1x*v1x+v1y*v1y), m2 = Math.sqrt(v2x*v2x+v2y*v2y);
    return Math.acos(Math.max(-1, Math.min(1, dot/(m1*m2))));
  }
  const angC = vecAngle(A, C, B) * 180 / Math.PI;
  const angA = vecAngle(B, A, C) * 180 / Math.PI;
  const angB = vecAngle(A, B, C) * 180 / Math.PI;

  // Fermat point via Weiszfeld
  function fermatPt(P, Q, R) {
    const aP = vecAngle(Q, P, R), aQ = vecAngle(P, Q, R), aR = vecAngle(P, R, Q);
    if (aP >= 2*Math.PI/3) return [...P];
    if (aQ >= 2*Math.PI/3) return [...Q];
    if (aR >= 2*Math.PI/3) return [...R];
    let mx = (P[0]+Q[0]+R[0])/3, my = (P[1]+Q[1]+R[1])/3;
    for (let it = 0; it < 300; it++) {
      const dP = Math.sqrt((mx-P[0])**2+(my-P[1])**2) || 1e-10;
      const dQ = Math.sqrt((mx-Q[0])**2+(my-Q[1])**2) || 1e-10;
      const dR = Math.sqrt((mx-R[0])**2+(my-R[1])**2) || 1e-10;
      const wP=1/dP, wQ=1/dQ, wR=1/dR, ws=wP+wQ+wR;
      mx = (wP*P[0]+wQ*Q[0]+wR*R[0])/ws;
      my = (wP*P[1]+wQ*Q[1]+wR*R[1])/ws;
    }
    return [mx, my];
  }

  const fp = fermatPt(A, B, C);
  const sectorColors = ["rgba(231,76,60,0.12)", "rgba(39,174,96,0.12)", "rgba(52,152,219,0.12)"];
  const sectorStroke = ["#e74c3c", "#27ae60", "#3498db"];
  const vLabels = ["a", "b", "c"];
  const triVerts = [A, B, C];

  // --- Left panel title ---
  svg.append("text").attr("x", leftCx).attr("y", 16).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#555").text("Fermat point & trisector sectors");

  // Draw trisector sectors (rays from Fermat point through each vertex)
  const rayLen = 500;
  for (let i = 0; i < 3; i++) {
    const dx = triVerts[i][0] - fp[0], dy = triVerts[i][1] - fp[1];
    const len = Math.sqrt(dx*dx+dy*dy) || 1;
    const endX = fp[0] + dx/len*rayLen, endY = fp[1] + dy/len*rayLen;
    svg.append("line").attr("x1", fp[0]).attr("y1", fp[1]).attr("x2", endX).attr("y2", endY)
      .attr("stroke", sectorStroke[i]).attr("stroke-width", 1.5).attr("stroke-dasharray", "6,3").attr("opacity", 0.6);
  }
  // Sector fills
  for (let i = 0; i < 3; i++) {
    const j = (i+1)%3;
    const dx1 = triVerts[i][0]-fp[0], dy1 = triVerts[i][1]-fp[1];
    const dx2 = triVerts[j][0]-fp[0], dy2 = triVerts[j][1]-fp[1];
    const l1 = Math.sqrt(dx1*dx1+dy1*dy1)||1, l2 = Math.sqrt(dx2*dx2+dy2*dy2)||1;
    svg.append("polygon")
      .attr("points", `${fp[0]},${fp[1]} ${fp[0]+dx1/l1*rayLen},${fp[1]+dy1/l1*rayLen} ${fp[0]+dx2/l2*rayLen},${fp[1]+dy2/l2*rayLen}`)
      .attr("fill", sectorColors[i]).attr("opacity", 0.5);
  }

  // Triangle
  svg.append("polygon")
    .attr("points", triVerts.map(v=>v.join(",")).join(" "))
    .attr("fill", "rgba(255,255,255,0.6)")
    .attr("stroke", needsSplit ? "#e74c3c" : "#2c3e50").attr("stroke-width", 2);

  // Fermat point
  const fpIsAtVertex = (Math.abs(fp[0]-C[0]) < 2 && Math.abs(fp[1]-C[1]) < 2) ||
                        (Math.abs(fp[0]-A[0]) < 2 && Math.abs(fp[1]-A[1]) < 2) ||
                        (Math.abs(fp[0]-B[0]) < 2 && Math.abs(fp[1]-B[1]) < 2);
  svg.append("circle").attr("cx", fp[0]).attr("cy", fp[1]).attr("r", 5)
    .attr("fill", fpIsAtVertex ? "#e74c3c" : "#8e44ad");
  svg.append("text").attr("x", fp[0]+10).attr("y", fp[1]-8)
    .attr("font-size", 12).attr("fill", fpIsAtVertex ? "#e74c3c" : "#8e44ad").attr("font-weight", "bold")
    .text(fpIsAtVertex ? "m (at vertex!)" : "m");

  // Angle arcs at Fermat point — show sector angles
  for (let i = 0; i < 3; i++) {
    const j = (i+1)%3;
    const a1 = Math.atan2(triVerts[i][1]-fp[1], triVerts[i][0]-fp[0]);
    const a2 = Math.atan2(triVerts[j][1]-fp[1], triVerts[j][0]-fp[0]);
    let da = a2 - a1;
    while (da < -Math.PI) da += 2*Math.PI;
    while (da > Math.PI) da -= 2*Math.PI;
    const angleDeg = Math.abs(da) * 180 / Math.PI;
    const ar = 25;
    const mid = a1 + da/2;
    svg.append("text")
      .attr("x", fp[0] + (ar+12)*Math.cos(mid))
      .attr("y", fp[1] + (ar+12)*Math.sin(mid))
      .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#8e44ad")
      .text(`${angleDeg.toFixed(0)}\u00B0`);
  }

  // Vertex labels and angle labels at vertices
  for (let i = 0; i < 3; i++) {
    svg.append("circle").attr("cx", triVerts[i][0]).attr("cy", triVerts[i][1]).attr("r", 5)
      .attr("fill", sectorStroke[i]);
    const offY = i < 2 ? 18 : -12;
    const offX = i === 0 ? -14 : i === 1 ? 14 : 0;
    svg.append("text").attr("x", triVerts[i][0]+offX).attr("y", triVerts[i][1]+offY)
      .attr("text-anchor", "middle").attr("font-size", 13).attr("fill", sectorStroke[i]).attr("font-weight", "bold")
      .text(vLabels[i]);
  }
  // Angle arc at C
  const arcR = 20;
  const aCx = Math.atan2(A[1]-C[1], A[0]-C[0]);
  const bCx = Math.atan2(B[1]-C[1], B[0]-C[0]);
  let startAng = Math.min(aCx, bCx), endAng = Math.max(aCx, bCx);
  if (endAng - startAng > Math.PI) { const tmp = startAng; startAng = endAng; endAng = tmp + 2*Math.PI; }
  const arcPts = d3.range(startAng, endAng, 0.02).map(a => [C[0]+arcR*Math.cos(a), C[1]+arcR*Math.sin(a)]);
  if (arcPts.length > 1) {
    svg.append("path").attr("d", d3.line()(arcPts))
      .attr("fill", "none").attr("stroke", "#f39c12").attr("stroke-width", 1.5);
  }
  svg.append("text").attr("x", C[0]).attr("y", C[1]+35)
    .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#f39c12")
    .text(`${angC.toFixed(0)}\u00B0`);

  // --- Right panel: after virtual split (or explanation) ---
  const rightX = leftW + 30;
  const rightW2 = w - rightX - 10;
  const rightCx = rightX + rightW2/2;

  svg.append("text").attr("x", rightCx).attr("y", 16).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#555")
    .text(needsSplit ? "After virtual split" : "No split needed");

  if (needsSplit) {
    // Split triangle along altitude from C to AB
    const footX = C[0], footY = baseY;  // AB is horizontal
    // Remap to right panel
    const rA = [rightCx - halfBase, baseY];
    const rB = [rightCx + halfBase, baseY];
    const rC = [rightCx, baseY - Math.max(25, triH)];
    const rD = [rightCx, baseY]; // foot of altitude

    // Draw the two sub-triangles
    // Sub-tri 1: A, D, C (left)
    svg.append("polygon")
      .attr("points", `${rA[0]},${rA[1]} ${rD[0]},${rD[1]} ${rC[0]},${rC[1]}`)
      .attr("fill", "rgba(39,174,96,0.08)").attr("stroke", "#27ae60").attr("stroke-width", 2);
    // Sub-tri 2: D, B, C (right)
    svg.append("polygon")
      .attr("points", `${rD[0]},${rD[1]} ${rB[0]},${rB[1]} ${rC[0]},${rC[1]}`)
      .attr("fill", "rgba(52,152,219,0.08)").attr("stroke", "#3498db").attr("stroke-width", 2);

    // Altitude line
    svg.append("line").attr("x1", rC[0]).attr("y1", rC[1]).attr("x2", rD[0]).attr("y2", rD[1])
      .attr("stroke", "#e74c3c").attr("stroke-width", 2).attr("stroke-dasharray", "5,3");

    // Compute Fermat points for each sub-triangle
    const fp1 = fermatPt(rA, rD, rC);
    const fp2 = fermatPt(rD, rB, rC);

    // Draw sectors for sub-tri 1
    const subTris = [[rA, rD, rC, fp1, "#27ae60"], [rD, rB, rC, fp2, "#3498db"]];
    for (const [sA, sB, sC, sFp, col] of subTris) {
      const sVerts = [sA, sB, sC];
      const srLen = 200;
      for (let i = 0; i < 3; i++) {
        const dx = sVerts[i][0]-sFp[0], dy = sVerts[i][1]-sFp[1];
        const len = Math.sqrt(dx*dx+dy*dy)||1;
        svg.append("line").attr("x1", sFp[0]).attr("y1", sFp[1])
          .attr("x2", sFp[0]+dx/len*srLen).attr("y2", sFp[1]+dy/len*srLen)
          .attr("stroke", col).attr("stroke-width", 1).attr("stroke-dasharray", "4,3").attr("opacity", 0.4);
      }
      svg.append("circle").attr("cx", sFp[0]).attr("cy", sFp[1]).attr("r", 4).attr("fill", "#8e44ad");
    }

    // Compute max angle in each sub-triangle
    const maxAng1 = Math.max(vecAngle(rD,rA,rC), vecAngle(rA,rD,rC), vecAngle(rA,rC,rD)) * 180/Math.PI;
    const maxAng2 = Math.max(vecAngle(rB,rD,rC), vecAngle(rD,rB,rC), vecAngle(rD,rC,rB)) * 180/Math.PI;

    // Labels
    svg.append("circle").attr("cx", rD[0]).attr("cy", rD[1]).attr("r", 5).attr("fill", "#e74c3c");
    svg.append("text").attr("x", rD[0]+10).attr("y", rD[1]+15)
      .attr("font-size", 12).attr("fill", "#e74c3c").attr("font-weight", "bold").text("d");

    for (const [pt, lab, offX, offY] of [[rA,"a",-14,18],[rB,"b",14,18],[rC,"c",0,-12]]) {
      svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 4).attr("fill", "#2c3e50");
      svg.append("text").attr("x", pt[0]+offX).attr("y", pt[1]+offY)
        .attr("text-anchor", "middle").attr("font-size", 12).attr("fill", "#2c3e50").text(lab);
    }

    svg.append("text").attr("x", rightCx).attr("y", baseY + 30).attr("text-anchor", "middle")
      .attr("font-size", 11).attr("fill", "#666")
      .text(`max angles: ${maxAng1.toFixed(0)}\u00B0 and ${maxAng2.toFixed(0)}\u00B0 (both \u2264 90\u00B0)`);
    svg.append("text").attr("x", rightCx).attr("y", baseY + 46).attr("text-anchor", "middle")
      .attr("font-size", 11).attr("fill", "#8e44ad")
      .text("each sub-triangle gets its own Fermat point");
  } else {
    // No split needed: show the healthy triangle with annotation
    const rA = [rightCx - halfBase, baseY];
    const rB = [rightCx + halfBase, baseY];
    const rC = [rightCx, baseY - Math.max(25, triH)];
    svg.append("polygon")
      .attr("points", `${rA[0]},${rA[1]} ${rB[0]},${rB[1]} ${rC[0]},${rC[1]}`)
      .attr("fill", "rgba(39,174,96,0.08)").attr("stroke", "#27ae60").attr("stroke-width", 2);

    const rfp = fermatPt(rA, rB, rC);
    // Sector angles at Fermat point
    const rVerts = [rA, rB, rC];
    for (let i = 0; i < 3; i++) {
      const j = (i+1)%3;
      const a1 = Math.atan2(rVerts[i][1]-rfp[1], rVerts[i][0]-rfp[0]);
      const a2 = Math.atan2(rVerts[j][1]-rfp[1], rVerts[j][0]-rfp[0]);
      let da = a2 - a1;
      while (da < -Math.PI) da += 2*Math.PI;
      while (da > Math.PI) da -= 2*Math.PI;
      const angleDeg = Math.abs(da) * 180/Math.PI;
      const ar2 = 30;
      const mid2 = a1 + da/2;
      svg.append("text")
        .attr("x", rfp[0]+(ar2+12)*Math.cos(mid2))
        .attr("y", rfp[1]+(ar2+12)*Math.sin(mid2))
        .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "#8e44ad")
        .text(`${angleDeg.toFixed(0)}\u00B0`);
    }
    svg.append("circle").attr("cx", rfp[0]).attr("cy", rfp[1]).attr("r", 5).attr("fill", "#8e44ad");
    svg.append("text").attr("x", rfp[0]+10).attr("y", rfp[1]-8)
      .attr("font-size", 12).attr("fill", "#8e44ad").text("m");

    for (const [pt, lab, offX, offY] of [[rA,"a",-14,18],[rB,"b",14,18],[rC,"c",0,-12]]) {
      svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 4).attr("fill", "#2c3e50");
      svg.append("text").attr("x", pt[0]+offX).attr("y", pt[1]+offY)
        .attr("text-anchor", "middle").attr("font-size", 12).attr("fill", "#2c3e50").text(lab);
    }

    svg.append("text").attr("x", rightCx).attr("y", baseY + 30).attr("text-anchor", "middle")
      .attr("font-size", 11).attr("fill", "#27ae60")
      .text(`all sectors \u2248 120\u00B0 \u2014 rich feasible space`);
  }

  // Bottom status
  const statusColor = needsSplit ? "#e74c3c" : obtuse_angle > 90 ? "#f39c12" : "#27ae60";
  const statusText = obtuse_angle >= 120
    ? `${angC.toFixed(0)}\u00B0 \u2265 120\u00B0 \u2014 Fermat point collapses to vertex c (degenerate!)`
    : needsSplit
    ? `${angC.toFixed(0)}\u00B0 > 100\u00B0 \u2014 sectors too narrow, virtual split required`
    : obtuse_angle > 90
    ? `${angC.toFixed(0)}\u00B0 > 90\u00B0 but \u2264 100\u00B0 \u2014 Fermat point still well inside`
    : `${angC.toFixed(0)}\u00B0 \u2264 90\u00B0 \u2014 ideal, well-balanced sectors`;

  svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", statusColor).attr("font-weight", "bold")
    .text(statusText);

  return svg.node();
}

“The split point can be expressed in barycentric coordinates and thus is linear in the original vertices. In practice we virtually split all triangles exhibiting an angle larger than $100°$.”

After the split, no angle exceeds $90°$, guaranteeing the Fermat point exists inside each sub-triangle.

4.5 The Linearized Orientation Constraints

Combining the trisector construction with normalization for floating-point precision gives 6 linear inequality constraints per triangle:

\[\delta_i^0 \left( (\mathbf{u}_i - \mathbf{m}) \cdot \frac{\mathbf{s}_i^\perp}{\|\mathbf{s}_i^\perp\|} - \epsilon \right) \geq 0\]

\[\delta_i^1 \left( -(\mathbf{u}_i - \mathbf{m}) \cdot \frac{\mathbf{s}_{i+1}^\perp}{\|\mathbf{s}_{i+1}^\perp\|} - \epsilon \right) \geq 0 \tag{4}\]

where $2\epsilon$ is a lower bound on the distance between two vertices of the same triangle, and $\delta_i^0, \delta_i^1$ are normalization factors. This gives us a convex formulation — linear constraints with a quadratic objective — which modern MIQP solvers can handle.

5. Complexity Reduction

5.1 The Problem

“One reason for the enormous dimension of discrete variables in problem (P1) is that potentially every vertex in the input triangle mesh can represent a singularity in the quad mesh.”

A typical input mesh has tens of thousands of vertices, each contributing discrete variables. Branch-and-cut solvers struggle with search spaces this large.

5.2 The Key Insight: Decimate in Parametric Space

The idea is to exploit the continuous relaxation $f_{\text{relaxed}}$ (the solution ignoring integer constraints). In flat regions far from singularities, the continuous relaxation is already close to an integer solution. Only near singularities does the integer snapping matter.

“Accordingly the idea is to perform decimation in flat parametric space and thus being able to represent the input surface by an extremely coarse mesh.”

The decimation is performed in the 2D parametric space of $\mathcal{M}_{\text{relaxed}}$, not in 3D. After decimation, each remaining triangle in the coarse mesh $\mathcal{M}_-$ overlaps with a set of triangles from $\mathcal{M}_{\text{relaxed}}$ in parameter space.

5.3 Decimation Operations

The mesh decimation uses three standard operations, all performed as an isometric re-tessellation of flat 2D space:

Halfedge collapse — merge an edge, removing one vertex (sorted by increasing edge length)
Edge flip — swap the diagonal of a quad formed by two triangles (to restore Delaunay criterion)
Edge split — subdivide an edge into two

Operations that move cone singularities are forbidden. This preserves the critical topological structure while aggressively simplifying flat regions.

Click an edge to apply the selected operation. Collapse merges the edge’s two endpoints into one, Flip swaps the shared diagonal between two adjacent triangles, and Split inserts a new vertex at the edge midpoint. The star vertex (orange) is a singularity and cannot be moved or removed:

Code

viewof decim_op = Inputs.radio(["Halfedge collapse", "Edge flip", "Edge split"], {value: "Halfedge collapse", label: "Operation"})

Figure 8

Code

{
  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 node = svg.node();

  // --- Mesh data structure (stored on SVG node for persistence) ---
  if (!node.__mesh || decim_op !== node.__lastOp) {
    // Reinitialize mesh when operation changes, to start fresh
    // Vertices: a structured grid-ish mesh with one singularity
    const pts = [
      [100, 460], [250, 460], [400, 460], [550, 460],       // 0-3 bottom row
      [100, 320], [250, 310], [400, 320], [550, 320],       // 4-7 mid row
      [175, 180], [350, 170], [525, 190],                    // 8-10 upper row
      [350, 60],                                              // 11 top
    ];
    // Triangles as vertex-index triples (CCW in SVG coords = CW math, but that's fine for viz)
    const tris = [
      [0, 1, 4], [1, 5, 4], [1, 2, 5], [2, 6, 5],
      [2, 3, 6], [3, 7, 6], [4, 5, 8], [5, 9, 8],
      [5, 6, 9], [6, 10, 9], [6, 7, 10], [8, 9, 11],
      [9, 10, 11],
    ];
    node.__mesh = {
      verts: pts.map(p => [...p]),
      tris: tris.map(t => [...t]),
      singularity: 9,  // vertex 9 is a cone singularity
      message: "Click an edge to apply the operation",
      messageColor: "#555",
    };
    node.__lastOp = decim_op;
  }
  const mesh = node.__mesh;
  const verts = mesh.verts;
  const tris = mesh.tris;
  const singIdx = mesh.singularity;

  // --- Edge utilities ---
  function edgeKey(a, b) { return a < b ? `${a}-${b}` : `${b}-${a}`; }

  // Build edge-to-triangles map
  function buildEdgeMap() {
    const em = new Map();
    for (let ti = 0; ti < tris.length; ti++) {
      const t = tris[ti];
      for (let j = 0; j < 3; j++) {
        const k = edgeKey(t[j], t[(j+1)%3]);
        if (!em.has(k)) em.set(k, []);
        em.get(k).push(ti);
      }
    }
    return em;
  }

  // --- Operations ---
  function doCollapse(a, b) {
    // Collapse edge a-b: merge b into a (keep a, remove b)
    // If b is the singularity, swap so we keep the singularity
    if (b === singIdx) { const tmp = a; a = b; b = tmp; }
    // Cannot collapse if a is the singularity (it would move it)
    if (a === singIdx) {
      mesh.message = "Cannot collapse: would move the singularity vertex";
      mesh.messageColor = "#e74c3c";
      return;
    }
    // Move a to midpoint
    verts[a] = [(verts[a][0]+verts[b][0])/2, (verts[a][1]+verts[b][1])/2];
    // Replace all references to b with a
    for (const t of tris) {
      for (let j = 0; j < 3; j++) { if (t[j] === b) t[j] = a; }
    }
    // Remove degenerate triangles (where two or more indices are the same)
    for (let i = tris.length - 1; i >= 0; i--) {
      const t = tris[i];
      if (t[0] === t[1] || t[1] === t[2] || t[0] === t[2]) tris.splice(i, 1);
    }
    // Update singularity index if needed (it shouldn't change since we kept it)
    mesh.message = `Collapsed edge: merged vertex ${b} into ${a}`;
    mesh.messageColor = "#2980b9";
  }

  function doFlip(a, b) {
    const em = buildEdgeMap();
    const k = edgeKey(a, b);
    const adjTris = em.get(k);
    if (!adjTris || adjTris.length !== 2) {
      mesh.message = "Cannot flip: edge is on the boundary";
      mesh.messageColor = "#e74c3c";
      return;
    }
    // Find the two opposite vertices
    const t0 = tris[adjTris[0]], t1 = tris[adjTris[1]];
    let opp0 = -1, opp1 = -1;
    for (let j = 0; j < 3; j++) {
      if (t0[j] !== a && t0[j] !== b) opp0 = t0[j];
      if (t1[j] !== a && t1[j] !== b) opp1 = t1[j];
    }
    if (opp0 === -1 || opp1 === -1 || opp0 === opp1) {
      mesh.message = "Cannot flip: degenerate configuration";
      mesh.messageColor = "#e74c3c";
      return;
    }
    // Check convexity of the quad (opp0, a, opp1, b) — flip only if convex
    function cross(O, A, B) {
      return (A[0]-O[0])*(B[1]-O[1]) - (A[1]-O[1])*(B[0]-O[0]);
    }
    const quad = [opp0, a, opp1, b].map(i => verts[i]);
    let convex = true;
    for (let i = 0; i < 4; i++) {
      const c = cross(quad[i], quad[(i+1)%4], quad[(i+2)%4]);
      if (c >= 0) { convex = false; break; } // SVG coords: CW is positive
    }
    if (!convex) {
      mesh.message = "Cannot flip: quad is not convex";
      mesh.messageColor = "#e74c3c";
      return;
    }
    // Replace the two triangles with the flipped pair
    tris[adjTris[0]] = [opp0, opp1, a];
    tris[adjTris[1]] = [opp0, b, opp1];
    mesh.message = `Flipped edge: (${a},${b}) \u2192 (${opp0},${opp1})`;
    mesh.messageColor = "#2980b9";
  }

  function doSplit(a, b) {
    const em = buildEdgeMap();
    const k = edgeKey(a, b);
    const adjTris = em.get(k) || [];
    // Insert new vertex at midpoint
    const newIdx = verts.length;
    verts.push([(verts[a][0]+verts[b][0])/2, (verts[a][1]+verts[b][1])/2]);
    // Replace each adjacent triangle with two sub-triangles
    const toRemove = new Set(adjTris);
    const newTris = [];
    for (const ti of adjTris) {
      const t = tris[ti];
      // Find the opposite vertex
      let opp = -1;
      for (let j = 0; j < 3; j++) { if (t[j] !== a && t[j] !== b) opp = t[j]; }
      newTris.push([a, newIdx, opp]);
      newTris.push([newIdx, b, opp]);
    }
    // Remove old triangles (in reverse order) and add new ones
    const sorted = [...toRemove].sort((x, y) => y - x);
    for (const ti of sorted) tris.splice(ti, 1);
    for (const nt of newTris) tris.push(nt);
    mesh.message = `Split edge: inserted vertex ${newIdx} at midpoint of (${a},${b})`;
    mesh.messageColor = "#2980b9";
  }

  // --- Drawing ---
  function redraw() {
    svg.selectAll("*").remove();
    svg.append("rect").attr("width", w).attr("height", h).attr("fill", "#fafafa");

    const em = buildEdgeMap();

    // Draw triangles
    for (const t of tris) {
      const pts = t.map(i => verts[i]);
      svg.append("polygon")
        .attr("points", pts.map(p => p.join(",")).join(" "))
        .attr("fill", "rgba(52,152,219,0.06)").attr("stroke", "#bdc3c7").attr("stroke-width", 0.5);
    }

    // Draw edges as clickable hit targets
    const drawnEdges = new Set();
    for (const t of tris) {
      for (let j = 0; j < 3; j++) {
        const a = t[j], b = t[(j+1)%3];
        const k = edgeKey(a, b);
        if (drawnEdges.has(k)) continue;
        drawnEdges.add(k);
        const isBoundary = (em.get(k) || []).length === 1;

        // Visible edge
        svg.append("line")
          .attr("x1", verts[a][0]).attr("y1", verts[a][1])
          .attr("x2", verts[b][0]).attr("y2", verts[b][1])
          .attr("stroke", isBoundary ? "#7f8c8d" : "#2c3e50")
          .attr("stroke-width", isBoundary ? 1.5 : 1);

        // Wider invisible hit target
        svg.append("line")
          .attr("x1", verts[a][0]).attr("y1", verts[a][1])
          .attr("x2", verts[b][0]).attr("y2", verts[b][1])
          .attr("stroke", "transparent").attr("stroke-width", 14)
          .attr("cursor", "pointer")
          .attr("data-a", a).attr("data-b", b)
          .on("mouseenter", function() {
            d3.select(this.previousSibling)
              .attr("stroke", "#e67e22").attr("stroke-width", 3);
          })
          .on("mouseleave", function() {
            const ib = (em.get(edgeKey(+d3.select(this).attr("data-a"), +d3.select(this).attr("data-b"))) || []).length === 1;
            d3.select(this.previousSibling)
              .attr("stroke", ib ? "#7f8c8d" : "#2c3e50")
              .attr("stroke-width", ib ? 1.5 : 1);
          })
          .on("click", function() {
            const ea = +d3.select(this).attr("data-a");
            const eb = +d3.select(this).attr("data-b");
            if (decim_op === "Halfedge collapse") doCollapse(ea, eb);
            else if (decim_op === "Edge flip") doFlip(ea, eb);
            else doSplit(ea, eb);
            redraw();
          });
      }
    }

    // Draw vertices
    const usedVerts = new Set();
    for (const t of tris) for (const v of t) usedVerts.add(v);

    for (const vi of usedVerts) {
      const isSing = vi === singIdx;
      svg.append("circle")
        .attr("cx", verts[vi][0]).attr("cy", verts[vi][1])
        .attr("r", isSing ? 7 : 4)
        .attr("fill", isSing ? "#e67e22" : "#2c3e50")
        .attr("stroke", isSing ? "#d35400" : "none")
        .attr("stroke-width", isSing ? 2 : 0);
    }

    // Singularity label
    if (usedVerts.has(singIdx)) {
      svg.append("text")
        .attr("x", verts[singIdx][0] + 12).attr("y", verts[singIdx][1] - 10)
        .attr("font-size", 12).attr("fill", "#e67e22").attr("font-weight", "bold")
        .text("singularity");
    }

    // Stats
    svg.append("text").attr("x", w/2).attr("y", h - 30).attr("text-anchor", "middle")
      .attr("font-size", 12).attr("fill", "#888")
      .text(`${usedVerts.size} vertices \u00B7 ${tris.length} triangles \u00B7 ${drawnEdges.size} edges`);

    // Message
    svg.append("text").attr("x", w/2).attr("y", h - 8).attr("text-anchor", "middle")
      .attr("font-size", 13).attr("fill", mesh.messageColor).attr("font-weight", "bold")
      .text(mesh.message);

    // Legend
    svg.append("circle").attr("cx", 20).attr("cy", 18).attr("r", 5).attr("fill", "#e67e22").attr("stroke", "#d35400").attr("stroke-width", 1.5);
    svg.append("text").attr("x", 32).attr("y", 22).attr("font-size", 11).attr("fill", "#555").text("Singularity (protected)");
    svg.append("text").attr("x", w - 20).attr("y", 22).attr("text-anchor", "end")
      .attr("font-size", 11).attr("fill", "#999").text("Hover an edge, then click to apply");
  }

  redraw();
  return node;
}

5.4 Splitting Integer-Critical Edges

After decimation, some triangles with large interior angles can trap vertices in degenerate line configurations. An edge is integer-critical if it is incident to a triangle with three singular vertices where the angle opposing the edge exceeds $110°$.

“To not exclude the possibility of arranging the vertices in the integer-grid on a line, after decimation we split all integer-critical edges.”

Each split separates two singular vertices by a regular one, and the number of integer-critical edges decreases monotonically.

The visualization below shows a triangle with three singular vertices after decimation. When the opposing angle is large, the only integer arrangement that fits all three vertices is collinear — a degenerate triangle. Splitting the integer-critical edge inserts a regular vertex that breaks the degeneracy:

Code

viewof ic_opposing_angle = Inputs.range([70, 150], {value: 90, step: 1, label: "Opposing angle (degrees)"})

Figure 9

Code

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

  const isCritical = ic_opposing_angle > 110;
  const leftW = w * 0.48, gap = w * 0.04, rightX = leftW + gap;
  const rightW = w - rightX - 5;

  // --- Build a triangle with the given opposing angle at vertex C ---
  // Edge AB is the "integer-critical edge" (the one opposite to the large angle at C)
  const baseLen = 160;
  const baseY = 310;
  const angleRad = ic_opposing_angle * Math.PI / 180;
  const triH = (baseLen / 2) / Math.tan(angleRad / 2);

  // --- LEFT PANEL: The triangle in parametric space ---
  const leftCx = leftW / 2 + 10;
  const A = [leftCx - baseLen/2, baseY];
  const B = [leftCx + baseLen/2, baseY];
  const C = [leftCx, baseY - Math.max(20, triH)];

  svg.append("text").attr("x", leftCx).attr("y", 16).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#555").text("Triangle with 3 singular vertices");

  // Integer grid behind
  const gridSc = 55;
  const gridOx = leftCx - 3*gridSc, gridOy = baseY + gridSc;
  for (let i = -1; i <= 6; i++) {
    svg.append("line").attr("x1", gridOx+i*gridSc).attr("y1", 25).attr("x2", gridOx+i*gridSc).attr("y2", baseY+40)
      .attr("stroke", "#eee").attr("stroke-width", 0.5);
  }
  for (let j = -1; j <= 6; j++) {
    svg.append("line").attr("x1", gridOx-gridSc).attr("y1", gridOy-j*gridSc).attr("x2", gridOx+6*gridSc).attr("y2", gridOy-j*gridSc)
      .attr("stroke", "#eee").attr("stroke-width", 0.5);
  }
  // Integer lattice dots
  for (let i = 0; i <= 5; i++) {
    for (let j = 0; j <= 5; j++) {
      svg.append("circle").attr("cx", gridOx+i*gridSc).attr("cy", gridOy-j*gridSc).attr("r", 2).attr("fill", "#ddd");
    }
  }

  // Triangle
  svg.append("polygon")
    .attr("points", `${A[0]},${A[1]} ${B[0]},${B[1]} ${C[0]},${C[1]}`)
    .attr("fill", isCritical ? "rgba(231,76,60,0.08)" : "rgba(39,174,96,0.06)")
    .attr("stroke", isCritical ? "#e74c3c" : "#27ae60").attr("stroke-width", 2);

  // Highlight the integer-critical edge (AB) — the one opposite to the large angle
  svg.append("line").attr("x1", A[0]).attr("y1", A[1]).attr("x2", B[0]).attr("y2", B[1])
    .attr("stroke", isCritical ? "#e74c3c" : "#27ae60")
    .attr("stroke-width", isCritical ? 4 : 2.5);

  if (isCritical) {
    const midAB = [(A[0]+B[0])/2, (A[1]+B[1])/2];
    svg.append("text").attr("x", midAB[0]).attr("y", midAB[1]+20)
      .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#e74c3c").attr("font-weight", "bold")
      .text("integer-critical edge");
  }

  // Angle arc at C
  function vecAngle(p1, vertex, p2) {
    const v1x=p1[0]-vertex[0], v1y=p1[1]-vertex[1];
    const v2x=p2[0]-vertex[0], v2y=p2[1]-vertex[1];
    const dot=v1x*v2x+v1y*v2y;
    const m1=Math.sqrt(v1x*v1x+v1y*v1y), m2=Math.sqrt(v2x*v2x+v2y*v2y);
    return Math.acos(Math.max(-1,Math.min(1,dot/(m1*m2))));
  }
  const angC = vecAngle(A, C, B) * 180 / Math.PI;

  const arcR = 22;
  const aCx = Math.atan2(A[1]-C[1], A[0]-C[0]);
  const bCx = Math.atan2(B[1]-C[1], B[0]-C[0]);
  let startA = Math.min(aCx,bCx), endA = Math.max(aCx,bCx);
  if (endA - startA > Math.PI) { const tmp=startA; startA=endA; endA=tmp+2*Math.PI; }
  const arcPts = d3.range(startA, endA, 0.02).map(a => [C[0]+arcR*Math.cos(a), C[1]+arcR*Math.sin(a)]);
  if (arcPts.length > 1) {
    svg.append("path").attr("d", d3.line()(arcPts))
      .attr("fill", "none").attr("stroke", "#f39c12").attr("stroke-width", 1.5);
  }
  svg.append("text").attr("x", C[0]).attr("y", C[1]+38)
    .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#f39c12")
    .text(`${angC.toFixed(0)}\u00B0`);

  // Vertex labels — all are singularities (red dots)
  for (const [pt, lab, offX, offY] of [[A,"s\u2081",-16,16],[B,"s\u2082",16,16],[C,"s\u2083",0,-14]]) {
    svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 6).attr("fill", "#e74c3c");
    svg.append("text").attr("x", pt[0]+offX).attr("y", pt[1]+offY)
      .attr("text-anchor", "middle").attr("font-size", 13).attr("fill", "#e74c3c").attr("font-weight", "bold").text(lab);
  }

  // --- RIGHT PANEL: Integer grid snapping ---
  svg.append("text").attr("x", rightX + rightW/2).attr("y", 16).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#555")
    .text(isCritical ? "After splitting the critical edge" : "Integer snapping (no problem)");

  // Right panel integer grid
  const rSc = 70;
  const rOx = rightX + 30, rOy = baseY + 10;
  for (let i = 0; i <= 4; i++) {
    for (let j = 0; j <= 4; j++) {
      svg.append("circle").attr("cx", rOx+i*rSc).attr("cy", rOy-j*rSc).attr("r", 3).attr("fill", "#ccc");
    }
  }
  for (let i = 0; i <= 4; i++) {
    svg.append("line").attr("x1", rOx+i*rSc).attr("y1", rOy+10).attr("x2", rOx+i*rSc).attr("y2", rOy-4*rSc-10)
      .attr("stroke", "#e8e8e8").attr("stroke-width", 0.5);
    svg.append("line").attr("x1", rOx-10).attr("y1", rOy-i*rSc).attr("x2", rOx+4*rSc+10).attr("y2", rOy-i*rSc)
      .attr("stroke", "#e8e8e8").attr("stroke-width", 0.5);
  }

  if (isCritical) {
    // Show: before split (collinear), after split (non-degenerate)
    // The 3 singularities snapped to integer grid — forced collinear by the flat triangle
    const colA = [rOx + 0*rSc, rOy - 1*rSc];
    const colB = [rOx + 3*rSc, rOy - 1*rSc];
    const colC = [rOx + 1*rSc, rOy - 1*rSc]; // on the same line!

    // Ghost: the degenerate collinear arrangement
    svg.append("line").attr("x1", colA[0]).attr("y1", colA[1]).attr("x2", colB[0]).attr("y2", colB[1])
      .attr("stroke", "rgba(231,76,60,0.3)").attr("stroke-width", 2);
    svg.append("text").attr("x", (colA[0]+colB[0])/2).attr("y", colA[1]-12)
      .attr("text-anchor", "middle").attr("font-size", 10).attr("fill", "rgba(231,76,60,0.5)")
      .text("collinear (degenerate!)");

    for (const [pt,lab] of [[colA,"s\u2081"],[colB,"s\u2082"],[colC,"s\u2083"]]) {
      svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 5)
        .attr("fill", "none").attr("stroke", "rgba(231,76,60,0.4)").attr("stroke-width", 1.5);
    }

    // After split: midpoint of AB is now a regular vertex, free to move off the line
    const splitA = [rOx + 0*rSc, rOy - 2*rSc];
    const splitB = [rOx + 3*rSc, rOy - 2*rSc];
    const splitC = [rOx + 1*rSc, rOy - 3*rSc];
    const splitD = [rOx + 1.5*rSc, rOy - 2*rSc]; // regular vertex (midpoint of AB), NOT on integer grid

    // Two sub-triangles: ADC and DBC
    svg.append("polygon")
      .attr("points", `${splitA[0]},${splitA[1]} ${splitD[0]},${splitD[1]} ${splitC[0]},${splitC[1]}`)
      .attr("fill", "rgba(39,174,96,0.08)").attr("stroke", "#27ae60").attr("stroke-width", 2);
    svg.append("polygon")
      .attr("points", `${splitD[0]},${splitD[1]} ${splitB[0]},${splitB[1]} ${splitC[0]},${splitC[1]}`)
      .attr("fill", "rgba(52,152,219,0.08)").attr("stroke", "#3498db").attr("stroke-width", 2);

    // Singular vertices (must be on integers)
    for (const [pt,lab] of [[splitA,"s\u2081"],[splitB,"s\u2082"],[splitC,"s\u2083"]]) {
      svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 6).attr("fill", "#e74c3c");
      svg.append("text").attr("x", pt[0]).attr("y", pt[1]-12)
        .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#e74c3c").attr("font-weight", "bold").text(lab);
    }
    // Regular vertex (free to be non-integer)
    svg.append("circle").attr("cx", splitD[0]).attr("cy", splitD[1]).attr("r", 5).attr("fill", "#3498db");
    svg.append("text").attr("x", splitD[0]+12).attr("y", splitD[1]+4)
      .attr("font-size", 11).attr("fill", "#3498db").attr("font-weight", "bold").text("d (regular)");

    svg.append("text").attr("x", rightX + rightW/2).attr("y", baseY + 38).attr("text-anchor", "middle")
      .attr("font-size", 11).attr("fill", "#27ae60")
      .text("regular vertex d is free to move off the integer grid");
    svg.append("text").attr("x", rightX + rightW/2).attr("y", baseY + 54).attr("text-anchor", "middle")
      .attr("font-size", 11).attr("fill", "#27ae60")
      .text("\u2192 non-degenerate triangles are now possible");
  } else {
    // Non-critical: show a healthy integer arrangement
    const snapA = [rOx + 0*rSc, rOy - 1*rSc];
    const snapB = [rOx + 3*rSc, rOy - 1*rSc];
    const snapC = [rOx + 1*rSc, rOy - 2*rSc];

    svg.append("polygon")
      .attr("points", `${snapA[0]},${snapA[1]} ${snapB[0]},${snapB[1]} ${snapC[0]},${snapC[1]}`)
      .attr("fill", "rgba(39,174,96,0.08)").attr("stroke", "#27ae60").attr("stroke-width", 2);

    for (const [pt,lab] of [[snapA,"s\u2081"],[snapB,"s\u2082"],[snapC,"s\u2083"]]) {
      svg.append("circle").attr("cx", pt[0]).attr("cy", pt[1]).attr("r", 6).attr("fill", "#e74c3c");
      svg.append("text").attr("x", pt[0]).attr("y", pt[1]-12)
        .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#e74c3c").attr("font-weight", "bold").text(lab);
    }

    svg.append("text").attr("x", rightX + rightW/2).attr("y", baseY + 38).attr("text-anchor", "middle")
      .attr("font-size", 11).attr("fill", "#27ae60")
      .text("all 3 singularities snap to Z\u00B2 without collapsing");
  }

  // Status bar
  const statusColor = isCritical ? "#e74c3c" : "#27ae60";
  const statusText = isCritical
    ? `Opposing angle ${angC.toFixed(0)}\u00B0 > 110\u00B0 \u2014 integer-critical! 3 singular vertices forced collinear on Z\u00B2`
    : `Opposing angle ${angC.toFixed(0)}\u00B0 \u2264 110\u00B0 \u2014 non-collinear integer arrangement exists`;
  svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", statusColor).attr("font-weight", "bold")
    .text(statusText);

  return svg.node();
}

Why 110° is the threshold

Three integer-grid points form a non-degenerate triangle only if they are not collinear. When a triangle has all three vertices as singularities (constrained to $\mathbb{Z}^2$) and one angle exceeds $110°$, the triangle is so flat that the nearest valid integer arrangement places all three on a line. Splitting the edge opposite the obtuse angle inserts a regular vertex (not constrained to integers) that can sit freely off the integer grid, breaking the collinearity trap.

5.5 Maintaining the Map

During decimation, a function $\Psi_i : \mathbb{R}^2 \to T_{\mathcal{M}_{\text{relaxed}}} \times \mathbb{R}^2$ tracks the correspondence between the decimated mesh and the original. For each triangle $t_i$ in the coarse mesh, $\Psi_i(\mathbf{m}) = (t_j, \psi_i(\mathbf{m}))$ records which original triangle $t_j$ contains the barycenter and the local coordinates.

The following step-through shows the decimation pipeline on a simple example:

Step 1/4: Input mesh in parametric space

6. Singularity Separation

6.1 The Clustering Problem

“The reason is that singular vertices might get arbitrarily close to each other within the continuous relaxation, while their distance in the integer-grid $\mathbb{Z} \times \mathbb{Z}$ of a valid IGM can never be less than 1.”

When the continuous relaxation places two singularities very close together, the branch-and-cut solver must search a huge space to find an integer solution that separates them. This makes the solver spend most of its time exploring infeasible subtrees.

6.2 The Geodesic Disc Strategy

For each singular vertex $s_i$, grow a geodesic disc $G_i$ of radius $\sqrt{2}$ in the decimated mesh $\mathcal{M}_-$. For all pairs of singularities $(s_i, s_j)$ within each disc ($i < j$), trace the shortest geodesic path connecting them and add a linear separation constraint.

The geodesic distances are computed using the Fast Marching method driven by the circle update of Novotni and Klein (2002), which is exact for obstacle-free planar cases.

6.3 Separator Constraints

Along each geodesic path between a close pair of singularities, the separation is enforced in the axis with the larger difference:

\[\left\{ \text{sgn}(\Delta u)|_{\mathcal{M}_-} \right\} \cdot \Delta u \geq 1\] \[\text{or} \quad \left\{ \text{sgn}(\Delta v)|_{\mathcal{M}_-} \right\} \cdot \Delta v \geq 1 \tag{5}\]

where separation is performed along the $u$ axis if $|\Delta u| \geq |\Delta v|$ and along $v$ otherwise. The sign is fixed from the continuous relaxation, making the constraint linear.

Adjust the target edge length to see how singularity spacing affects the need for separation constraints:

Code

viewof sep_scale = Inputs.range([0.5, 3.0], {value: 1.5, step: 0.1, label: "Target edge length"})

Figure 10

Code

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

  const ox = 50, oy = 400, sc = 80;

  // Integer grid
  for (let i = 0; i <= 5; i++) {
    svg.append("line").attr("x1", ox + i*sc).attr("y1", oy).attr("x2", ox + i*sc).attr("y2", oy - 4*sc)
      .attr("stroke", "#eee").attr("stroke-width", 0.5);
    svg.append("line").attr("x1", ox).attr("y1", oy - i*sc).attr("x2", ox + 5*sc).attr("y2", oy - i*sc)
      .attr("stroke", "#eee").attr("stroke-width", 0.5);
  }
  // Integer points
  for (let ix = 0; ix <= 5; ix++) {
    for (let iy = 0; iy <= 4; iy++) {
      svg.append("circle").attr("cx", ox + ix*sc).attr("cy", oy - iy*sc).attr("r", 2).attr("fill", "#ddd");
    }
  }

  // Singularities — positions scale with target edge length
  const baseSings = [
    {x: 1.5, y: 1.8, label: "s\u2081"},
    {x: 2.2, y: 2.1, label: "s\u2082"},
    {x: 3.5, y: 1.0, label: "s\u2083"},
    {x: 3.8, y: 2.8, label: "s\u2084"},
  ];

  // As edge length decreases, singularities get closer in the parametric domain
  const sings = baseSings.map(s => ({
    x: s.x,
    y: s.y,
    label: s.label
  }));

  const discRadius = Math.sqrt(2);

  // Draw geodesic discs
  for (const s of sings) {
    const cx = ox + s.x * sc, cy = oy - s.y * sc;
    svg.append("circle")
      .attr("cx", cx).attr("cy", cy).attr("r", discRadius * sc)
      .attr("fill", "rgba(231,76,60,0.06)")
      .attr("stroke", "rgba(231,76,60,0.3)").attr("stroke-width", 1).attr("stroke-dasharray", "4,3");
  }

  // Check pairs and draw separators
  let numSeparators = 0;
  for (let i = 0; i < sings.length; i++) {
    for (let j = i + 1; j < sings.length; j++) {
      const dx = (sings[j].x - sings[i].x) * sep_scale;
      const dy = (sings[j].y - sings[i].y) * sep_scale;
      const dist = Math.sqrt(dx*dx + dy*dy);

      const px1 = ox + sings[i].x * sc, py1 = oy - sings[i].y * sc;
      const px2 = ox + sings[j].x * sc, py2 = oy - sings[j].y * sc;

      if (dist < discRadius * sep_scale) {
        // Need separator
        numSeparators++;
        svg.append("line")
          .attr("x1", px1).attr("y1", py1).attr("x2", px2).attr("y2", py2)
          .attr("stroke", "#e67e22").attr("stroke-width", 2.5);

        // Separation direction indicator
        const midX = (px1+px2)/2, midY = (py1+py2)/2;
        const sepAxis = Math.abs(dx) >= Math.abs(dy) ? "u" : "v";
        svg.append("text").attr("x", midX + 8).attr("y", midY - 5)
          .attr("font-size", 10).attr("fill", "#e67e22").text(`sep ${sepAxis}`);
      } else {
        // No separator needed
        svg.append("line")
          .attr("x1", px1).attr("y1", py1).attr("x2", px2).attr("y2", py2)
          .attr("stroke", "#27ae60").attr("stroke-width", 1).attr("stroke-dasharray", "3,3").attr("opacity", 0.5);
      }
    }
  }

  // Draw singularities
  for (const s of sings) {
    const cx = ox + s.x * sc, cy = oy - s.y * sc;
    svg.append("circle").attr("cx", cx).attr("cy", cy).attr("r", 7).attr("fill", "#e74c3c");
    svg.append("text").attr("x", cx).attr("y", cy - 12)
      .attr("text-anchor", "middle").attr("font-size", 13).attr("fill", "#e74c3c").attr("font-weight", "bold")
      .text(s.label);
  }

  // Info text
  svg.append("text").attr("x", w/2).attr("y", h - 5).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#666")
    .text(`Target edge length: ${sep_scale.toFixed(1)} — ${numSeparators} separator constraint(s) needed`);

  // Legend
  svg.append("line").attr("x1", 15).attr("y1", 18).attr("x2", 35).attr("y2", 18)
    .attr("stroke", "#e67e22").attr("stroke-width", 2.5);
  svg.append("text").attr("x", 40).attr("y", 22).attr("font-size", 11).attr("fill", "#666").text("Needs separator (too close)");

  svg.append("line").attr("x1", 15).attr("y1", 35).attr("x2", 35).attr("y2", 35)
    .attr("stroke", "#27ae60").attr("stroke-width", 1).attr("stroke-dasharray", "3,3");
  svg.append("text").attr("x", 40).attr("y", 39).attr("font-size", 11).attr("fill", "#666").text("Sufficiently separated");

  return svg.node();
}

Why separation helps the solver

The singularity separation constraints are cuts in the optimization literature — they don’t remove any optimal integer solutions, but they tighten the continuous relaxation so it’s closer to the best integer solution. This dramatically reduces the search tree that the branch-and-cut solver must explore. In experiments, without separation constraints the solver often fails to find any feasible solution within hours.

7. The Full Algorithm

7.1 The Six-Step Pipeline

Based on the described components, the algorithm can be understood as a series of bijective maps between meshes:

Globally smooth parametrization ($\mathcal{M} \to \mathcal{M}_{\text{relaxed}}$) — Solve the continuous relaxation of the quality metric using IPOPT
Decimation in parametric space ($\mathcal{M}_{\text{relaxed}} \to \mathcal{M}_-$) — Aggressively simplify flat regions, keeping singularities fixed
Addition of singularity separators ($\mathcal{M}_- \to \mathcal{M}^{\text{sep}}$) — Add geodesic separation constraints between close singularity pairs
MIQP optimization ($\mathcal{M}^{\text{sep}} \to \mathcal{M}_{\text{int}}$) — Solve the convex MIQP with branch-and-cut (CPLEX)
(Optional) Refinement ($\mathcal{M}_{\text{int}} \to \mathcal{M}_{\text{IGM}}$) — Split edges to recover DOFs for smoother iso-lines
Quad mesh extraction ($\mathcal{M}_{\text{IGM}} \to \mathcal{Q}$) — Contour the integer iso-lines to produce the final quad mesh

Each step is visualized in the pipeline step-through below:

Step 1/6: Input mesh M with cross-field and sizing

7.2 Implementation Details

Optimization of (P2) type problems is done in steps 1, 3, 4, and 5. Variables in each triangle are merged along a dual spanning tree, and an arbitrary singularity is fixed to $(0,0)$ to obtain a unique minimum.

Features and Boundaries are preserved by adding linear constraints ensuring feature edges map onto integer-grid lines. The feature graph consists of closed or open chains of mesh edges ending at feature vertices or singularities.

Numerical Solver: IPOPT for continuous optimization problems, CPLEX for mixed-integer problems. Programs are restricted to quadratic ($k = 2$).

Lazy Constraints: Orientation constraints (4) are typically violated for only a small fraction of triangles. The algorithm starts with an empty set and iteratively adds violated constraints until a valid solution is found, significantly reducing runtime.

8. Results

8.1 Importance of Decimation and Singularity Separation

“To be able to find any solution without our optimizations, we reduced the mesh from 15k to 3k vertices (by keeping the singularities). In the first test we used our decimation strategy without addition of singularity separators and not a single solution could be found within 8h.”

The experiments on the BOTIJO model demonstrate that both optimizations are essential:

Configuration	Result
No decimation, no separation	No solution in 8+ hours
Decimation only, no separation	No solution in 8+ hours
Separation only, no decimation	First low-quality solution after 20 min
Both decimation + separation	High-quality solution in 17 seconds

8.2 Comparison with Greedy Rounding Methods

The paper compares RMIQ (their method) against QuadCover (QC), Mixed-Integer Quadrangulation (MIQ), and MIQ with stiffening on the ROCKERARM model at various target edge lengths $h$:

Algorithm	$h = 10^{-4}$	$h = 0.025$	$h = 0.05$	$h = 0.2$
QuadCover	8 degen	51 degen	82 degen	2047 degen
MIQ	8 degen	18 degen	42 degen	568 degen
MIQ Stiff	0 degen	0 degen	1 degen	541 degen
RMIQ	0 degen	0 degen	0 degen	0 degen

RMIQ produces zero degenerate faces at all scales, including coarse layouts where greedy methods completely fail.

8.3 Coarse Quad Layouts

For coarse quad layout generation (target edge length = 20% of bounding box), RMIQ is compared against Dual Loops meshing:

Model	Singularities	Patches (Dual Loops)	Patches (RMIQ)	Patches (Enhanced CF)
BOTIJO	72	221	175	112
BLOCK	48	76	76	—
GUY	40	168	55	—
ELK	52	86	62	58
ROCKERARM	30	115	74	66
FERTILITY	48	98	117	85

The global optimization approach produces quad layouts with fewer patches compared to the greedy Dual Loops method, while maintaining alignment with the cross field. Slight improvements to the cross field (manually repositioning singularities) can further reduce patch counts.

9. Connection to Our Work

This paper addresses the parameterization stage of the quad meshing pipeline — the step that converts a cross field and sizing field into a global parameterization from which the quad mesh is extracted. In our MetricFrameField pipeline, this corresponds to:

Pipeline Stage	Bommes et al. 2013	MetricFrameField
Input	Cross field + sizing field on triangle mesh	Cross field from metric-driven solve
Continuous relaxation	IPOPT solve of $E_\alpha^k$	MIQ-style parameterization
Integer constraints	MIQP with branch-and-cut (CPLEX)	Greedy rounding (MIQ)
Orientation guarantee	Convex trisector constraints	Not yet implemented
Complexity reduction	Parametric-space decimation	Not yet implemented
Singularity separation	Geodesic disc constraints	Not yet implemented
Output	Guaranteed valid IGM	Parameterization (may have degeneracies)

The key takeaway: the RMIQ approach provides a principled replacement for greedy rounding that guarantees a valid quad mesh. The trisector linearization (Sec. 4) and complexity reduction (Sec. 5) are the critical ingredients that make this practical. Integrating these ideas into our pipeline would eliminate the degeneracy failures that currently occur with coarse target sizing.

Paper at a Glance

1. Introduction

1.1 The Reliability Problem

1.2 Contributions

2. Integer-Grid Maps

2.1 Piecewise Linear Maps

2.2 Transition Functions

2.3 The Three IGM Conditions

3. The Quality Metric

3.1 The Energy \(E(f)\)

3.2 Generalized Anisotropy Metric

3.3 The Naive MINLP (Problem P1)

4. Convex Consistent Orientation

4.1 The Non-Convex Problem

4.2 The Trisector Construction

4.3 The Fermat Point

4.4 Virtual Splitting for Obtuse Triangles

4.5 The Linearized Orientation Constraints

5. Complexity Reduction

5.1 The Problem

5.2 The Key Insight: Decimate in Parametric Space

5.3 Decimation Operations

5.4 Splitting Integer-Critical Edges

5.5 Maintaining the Map

6. Singularity Separation

6.1 The Clustering Problem

6.2 The Geodesic Disc Strategy

6.3 Separator Constraints

7. The Full Algorithm

7.1 The Six-Step Pipeline

7.2 Implementation Details

7.3 Optional Refinement

8. Results

8.1 Importance of Decimation and Singularity Separation

8.2 Comparison with Greedy Rounding Methods

8.3 Coarse Quad Layouts

9. Connection to Our Work