Skip to content

Fix sign convention in BFF::dirichletToNeumann - #264

Open
joaoajmatos wants to merge 1 commit into
nmwsharp:masterfrom
joaoajmatos:fix/bff-dtn-sign
Open

joaoajmatos wants to merge 1 commit into
nmwsharp:masterfrom
joaoajmatos:fix/bff-dtn-sign

Conversation

@joaoajmatos

Copy link
Copy Markdown

Fixes #263.

dirichletToNeumann() returned the exterior angles for -u while
computeBoundaryPositions() lays out with exp(+u), so any non-constant scale
factors (cones) came out with the wrong sign: a 0.5 rad notch became a 0.5 rad
overlap with inverted triangles at the cone, and the cut sides came out short. flatten() and
flattenFromExteriorAngles() were not affected. Algebra and a before/after
figure are in the issue.

Change: two signs in dirichletToNeumann(), plus a comment stating the
convention. No API change.

Tests: new test/src/boundary_first_flattening_test.cpp (registered in
test/CMakeLists.txt), three cases: a DtN/NtD round trip, scale-factor vs
exterior-angle layouts on the same u, and a flat unit disk with one interior
cone of angle 0.5 cut to the boundary (notch angle, cut-side lengths, outline).
All three fail on v1.1.0 and pass with this change; the full suite goes from
198 to 201 passing.

dirichletToNeumann() returned the exterior angles for -u, while
computeBoundaryPositions() and neumannToDirichlet() use lengths scaled by
exp(+u). flattenFromScaleFactors() therefore produced inconsistent layouts
for any non-constant u, e.g. when placing cone singularities and cutting to
the boundary before flattening. flatten() (u = 0) was unaffected.

Adds test/src/boundary_first_flattening_test.cpp: DtN/NtD round trip,
scale-factor vs exterior-angle layouts, and a flat disk with one cone cut to
the boundary. All three fail on v1.1.0 and pass with this change.

Fixes nmwsharp#263
@joaoajmatos

Copy link
Copy Markdown
Author

Sign convention error in BFF::dirichletToNeumann

geometry-central v1.1.0, src/surface/boundary_first_flattening.cpp. Companion to #263 / #264.

BFF::flattenFromScaleFactors(u) lays out a mesh with the sign of the conformal factor flipped between its two halves: dirichletToNeumann() returns the boundary exterior angles that belong to $-u$, while computeBoundaryPositions() scales edge lengths by $e^{+u}$. The default path flatten() passes $u = 0$ and is unaffected, which is why the error has gone unnoticed. With a cone singularity of angle $\alpha$ the layout behaves as if the cone had angle $-\alpha$: the cut that should open into a notch overlaps itself instead, and the first ring of triangles at the cone inverts. The fix is two signs in one function. Three regression tests fail on v1.1.0 and pass with the fix; the existing suite is unchanged (198 → 201 passing).

1. Setting

Boundary First Flattening (Sawhney & Crane 2017) computes a conformal flattening of a disk-topology mesh from boundary data. Write the target metric as $\tilde g = e^{2u} g$, so lengths scale by $e^{u}$. On a surface with Gaussian curvature $K$ the log conformal factor satisfies the Yamabe equation

$$\Delta u = K - \tilde K,$$

with $\tilde K = 0$ in the interior for a flattening. Discretely, with angle defects $\Omega$ in place of integrated curvature and the cotan Laplacian $L$ (positive semidefinite in geometry-central, so $L = -\Delta$), this reads

$$L,u = -(\Omega - \tilde\Omega).$$

Along the boundary the role of $\tilde K$ is played by the target exterior angles $k$, and the two boundary problems BFF offers are the Dirichlet-to-Neumann map $u_b \mapsto k$ and its inverse $k \mapsto u_b$.

In the code, Omegai holds $2\pi - \theta_v$ at interior vertices, Omegab holds $\pi - \theta_v$ at boundary vertices, and $L$ is partitioned into Lii, Lib, Lbb.

2. What the three functions compute

Layout, computeBoundaryPositions:

targetLength(iV) = geo.edgeLengths[e] * exp(0.5 * (uBdy(tail) + uBdy(tip)));

Edge lengths are scaled by $e^{+u}$. This fixes the convention: $u$ is the log scale factor of the target metric.

Neumann to Dirichlet, neumannToDirichlet:

rhs          = [Omegai ; Omegab - kBdy]
fullSolution = -Lsolver->solve(rhs)

That is $u = -L^{-1}(\Omega - \tilde\Omega)$, i.e. $L u = -(\Omega - \tilde\Omega)$, consistent with $e^{+u}$ above.

Dirichlet to Neumann, dirichletToNeumann, as shipped:

return Omegab - (Lib.transpose() * Liisolver->solve(Omegai - Lib * uBdy)) - Lbb * uBdy;

3. Derivation of the correct map

Split $L u = -(\Omega - \tilde\Omega)$ into interior rows (where $\tilde\Omega = 0$) and boundary rows (where $\tilde\Omega = k$):

$$L_{ii} u_i + L_{ib} u_b = -\Omega_i, \qquad L_{bi} u_i + L_{bb} u_b = -(\Omega_b - k).$$

Eliminating $u_i$:

$$u_i = -L_{ii}^{-1}(\Omega_i + L_{ib} u_b), \qquad k = \Omega_b - L_{bi} L_{ii}^{-1}(\Omega_i + L_{ib} u_b) + L_{bb} u_b.$$

The shipped expression is

$$k_{\text{shipped}}(u_b) = \Omega_b - L_{bi} L_{ii}^{-1}(\Omega_i - L_{ib} u_b) - L_{bb} u_b = k(-u_b).$$

So the function returns the exterior angles of the metric $e^{-2u} g$ while the layout uses $e^{+2u} g$. Two special cases hide the error:

  • $u = 0$: both expressions coincide. This is the flatten() path.
  • $u$ constant: $L\mathbf{1} = 0$, so every $u_b$-dependent term vanishes in both expressions. A uniform rescale is invisible to the exterior angles, as it should be.

flattenFromExteriorAngles() goes through neumannToDirichlet() and is correct. Only flattenFromScaleFactors() is affected, and it is the entry point used for cone singularities.

4. Consequence for cones

Take a flat unit disk and prescribe one cone of angle $\alpha = 0.5$ at the centre by solving the Yamabe equation for $u$ with $u = 0$ on the outer circle. Cut from the centre to the boundary along one ray and flatten. The correct result is a disk with a wedge of angle 0.5 removed; the two cut sides have equal length

$$\ell = \frac{2\pi}{2\pi - 0.5} \approx 1.08646.$$

With the sign flipped the cone is laid out as if $\alpha = -0.5$: instead of removing a wedge the layout gains 0.5 rad of material, so the cut sides cross and the layout overlaps itself. The first ring of 16 triangles around the cone inverts. Both cut sides come out equal to each other but short (1.0805 on the unit disk). On a 5k-triangle sphere cap with a single cone the cut sides are about 20% short.

Flat unit disk with one cone of angle 0.5 cut to the boundary. Left: input. Middle: v1.1.0, the cut sides cross and the layout overlaps itself by 0.5 rad (darker wedge), with inverted triangles at the cone. Right: patched, a 0.5 rad notch with both cut sides at 1.08612 (analytic 1.08646 on the continuum sector).

Quantity Analytic v1.1.0 Patched
Cut side a 1.08646 1.08050 1.08612
Cut side b 1.08646 1.08050 1.08612
Angle at the cut notch 0.5 overlap 0.5 notch 0.5
Inverted triangles at the cone 0 16 0

The residual 0.03% between the patched side length and the analytic sector comes from the mesh (20 rings × 64 sectors) approximating the continuum cone; it shrinks with refinement.

5. Fix

Vector<double> BFF::dirichletToNeumann(const Vector<double>& uBdy) {
  // Interior: Lii u_i = -(Omegai + Lib u_b) (zero target curvature, lengths scale by e^{+u}
  // as in computeBoundaryPositions and neumannToDirichlet). Boundary: k = Omegab + Lbi u_i + Lbb u_b.
  return Omegab - (Lib.transpose() * Liisolver->solve(Omegai + Lib * uBdy)) + Lbb * uBdy;
}

Two signs change. No API change, and no change to flatten() results.

6. Regression tests

New file test/src/boundary_first_flattening_test.cpp, registered in test/CMakeLists.txt. The mesh is a flat unit disk built as a polar grid: a centre vertex, 20 rings of 64 vertices. A "slit" variant duplicates the ring vertices on the ray $\theta = 0$, so the centre becomes a boundary vertex and the mesh stays a topological disk. Scale factors for one interior cone of angle 0.5 come from the Yamabe equation on the uncut disk with $u = 0$ on the outer circle.

Test What it checks v1.1.0 Patched
DirichletNeumannRoundTrip $u_b \to k \to u_b$ returns the input (mean-zero), tolerance 1e-8 FAIL (0.7275) PASS
ScaleFactorsMatchExteriorAngles flattenFromScaleFactors(u) equals flattenFromExteriorAngles(k(u)) up to rigid motion and the uniform scale $e^{\bar u}$, tolerance 1e-8 FAIL (0.1803) PASS
ConeFromScaleFactorsIsIsometricOnTheCut notch angle equals the cone angle; cut sides equal; outline length within 0.00628 of 3D FAIL (notch off by 2.3555) PASS

Full suite:

Tree Tests Result
v1.1.0 + new tests (source unpatched) 201 198 pass, 3 fail (the new ones)
v1.1.0 + patch + new tests 201 201 pass

7. Reproducibility

# pristine v1.1.0 with the submodule (tarballs omit it)
git clone --branch v1.1.0 https://github.com/nmwsharp/geometry-central.git
cd geometry-central && git submodule update --init

# A: tests only            B: full patch
git apply --include='test/*' bff-dtn-sign.patch      # or: git apply bff-dtn-sign.patch

cmake -S test -B build -G Ninja -DCMAKE_BUILD_TYPE=Release
ninja -C build geometry-central-test
./build/bin/geometry-central-test --gtest_filter='BFFSuite.*'
./build/bin/geometry-central-test

Toolchain used: Apple clang (Xcode 27.2), Eigen 3.4.0, Ninja. The BFF results were identical on two independent builds (v1.1.0 and current master).

8. Workaround for v1.1.0 users

Either call dirichletToNeumann(-u), or compute $k$ from the formula in section 3 and call flattenFromBoth(u, k) directly. The second form does not depend on the library version.

References

R. Sawhney and K. Crane, "Boundary First Flattening", ACM Transactions on Graphics 37(1), 2017. N. Sharp et al., geometry-central, https://geometry-central.net.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

BFF::flattenFromScaleFactors gives an overlap where it should give a notch (sign bug in dirichletToNeumann)

1 participant