From 81cd9a4a99139e428cef2dabe2576add65af3826 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 18 Aug 2026 14:09:42 +0000 Subject: [PATCH] fix(regularisers): correct diag_bending/diag_all corner cross-term The boundary-corrected diagonal in `diag_bending` (field 2d/3d, flow 2d/3d) and `diag_all` (flow 2d/3d) expanded each corner weight as w * (fx0*fy0 + fx1*fy0 + fx1*fy0 + fx1*fy1) double-counting `fx1*fy0` and dropping `fx0*fy1`. The correct expansion of (fx0+fx1)*(fy0+fy1) is w * (fx0*fy0 + fx1*fy0 + fx0*fy1 + fx1*fy1) which is exactly the four corner neighbours `matvec_bending` / `matvec_all` already accumulate a few lines above in each of the same files, one-for-one. The error cancels wherever both axes' one-sided boundary signs agree (fx0==fx1 and fy0==fy1 make both forms 4*fx0*fy0), and only manifests where they differ -- which is why the interior-only diagonal tests never caught it. It corrupts the shipped `field_diag`/`flow_diag` output at boundary voxels under sign-flipping boundary conditions (Zero/DST1/DST2) whenever bending is active, and hence any Jacobi preconditioner built from that diagonal. 12 terms across 6 functions in 4 files; a grep for the doubled-term form under regularisers/ now returns nothing. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_016AjQcY78NgbagPSbPJRr6Z --- regularisers/field/2d.h | 2 +- regularisers/field/3d.h | 6 +++--- regularisers/flow/2d.h | 4 ++-- regularisers/flow/3d.h | 12 ++++++------ 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/regularisers/field/2d.h b/regularisers/field/2d.h index 99bbf6e..27447ae 100755 --- a/regularisers/field/2d.h +++ b/regularisers/field/2d.h @@ -461,7 +461,7 @@ struct Kernels> w11 = kernel[5]; w00 -= w10 * (fx0 + fx1) + w01 * (fy0 + fy1) + w20 * (fx00 + fx11) + w02 * (fy00 + fy11) - + w11 * (fx0*fy0 + fx1*fy0 + fx1*fy0 + fx1*fy1); + + w11 * (fx0*fy0 + fx1*fy0 + fx0*fy1 + fx1*fy1); op(out, w00); }; diff --git a/regularisers/field/3d.h b/regularisers/field/3d.h index b393469..adec109 100755 --- a/regularisers/field/3d.h +++ b/regularisers/field/3d.h @@ -514,9 +514,9 @@ struct Kernels> w110 = kernel[7], w101 = kernel[8], w011 = kernel[9]; w000 -= w100 * (fx0 + fx1) + w010 * (fy0 + fy1) + w001 * (fz0 + fz1) + w200 * (fx00 + fx11) + w020 * (fy00 + fy11) + w002 * (fz00 + fz11) - + w110 * (fx0*fy0 + fx1*fy0 + fx1*fy0 + fx1*fy1) - + w101 * (fx0*fz0 + fx1*fz0 + fx1*fz0 + fx1*fz1) - + w011 * (fy0*fz0 + fy1*fz0 + fy1*fz0 + fy1*fz1); + + w110 * (fx0*fy0 + fx1*fy0 + fx0*fy1 + fx1*fy1) + + w101 * (fx0*fz0 + fx1*fz0 + fx0*fz1 + fx1*fz1) + + w011 * (fy0*fz0 + fy1*fz0 + fy0*fz1 + fy1*fz1); op(out, w000); }; diff --git a/regularisers/flow/2d.h b/regularisers/flow/2d.h index d3d517f..3eee37c 100755 --- a/regularisers/flow/2d.h +++ b/regularisers/flow/2d.h @@ -402,7 +402,7 @@ struct RegFlow { w110 = kernel[5]; w000 -= w100 * (fx0 + fx1) + w010 * (fy0 + fy1) + w200 * (fx00 + fx11) + w020 * (fy00 + fy11) - + w110 * (fx0*fy0 + fx1*fy0 + fx1*fy0 + fx1*fy1); + + w110 * (fx0*fy0 + fx1*fy0 + fx0*fy1 + fx1*fy1); op(out, w000); }; setdiag(out[0], kernel); @@ -681,7 +681,7 @@ struct RegFlow { w110 = kernel[5]; w000 -= w100 * (fx0 + fx1) + w010 * (fy0 + fy1) + w200 * (fx00 + fx11) + w020 * (fy00 + fy11) - + w110 * (fx0*fy0 + fx1*fy0 + fx1*fy0 + fx1*fy1); + + w110 * (fx0*fy0 + fx1*fy0 + fx0*fy1 + fx1*fy1); op(*out, w000); }; diff --git a/regularisers/flow/3d.h b/regularisers/flow/3d.h index 2779583..2c9da38 100755 --- a/regularisers/flow/3d.h +++ b/regularisers/flow/3d.h @@ -523,9 +523,9 @@ struct RegFlow { w110 = kernel[7], w101 = kernel[8], w011 = kernel[9]; w000 -= w100 * (fx0 + fx1) + w010 * (fy0 + fy1) + w001 * (fz0 + fz1) + w200 * (fx00 + fx11) + w020 * (fy00 + fy11) + w002 * (fz00 + fz11) - + w110 * (fx0*fy0 + fx1*fy0 + fx1*fy0 + fx1*fy1) - + w101 * (fx0*fz0 + fx1*fz0 + fx1*fz0 + fx1*fz1) - + w011 * (fy0*fz0 + fy1*fz0 + fy1*fz0 + fy1*fz1); + + w110 * (fx0*fy0 + fx1*fy0 + fx0*fy1 + fx1*fy1) + + w101 * (fx0*fz0 + fx1*fz0 + fx0*fz1 + fx1*fz1) + + w011 * (fy0*fz0 + fy1*fz0 + fy0*fz1 + fy1*fz1); op(out, w000); }; setdiag(out[0], kernel); @@ -942,9 +942,9 @@ struct RegFlow { w110 = kernel[7], w101 = kernel[8], w011 = kernel[9]; w000 -= w100 * (fx0 + fx1) + w010 * (fy0 + fy1) + w001 * (fz0 + fz1) + w200 * (fx00 + fx11) + w020 * (fy00 + fy11) + w002 * (fz00 + fz11) - + w110 * (fx0*fy0 + fx1*fy0 + fx1*fy0 + fx1*fy1) - + w101 * (fx0*fz0 + fx1*fz0 + fx1*fz0 + fx1*fz1) - + w011 * (fy0*fz0 + fy1*fz0 + fy1*fz0 + fy1*fz1); + + w110 * (fx0*fy0 + fx1*fy0 + fx0*fy1 + fx1*fy1) + + w101 * (fx0*fz0 + fx1*fz0 + fx0*fz1 + fx1*fz1) + + w011 * (fy0*fz0 + fy1*fz0 + fy0*fz1 + fy1*fz1); op(*out, w000); };