wip, trying to unwrap the torus

1 files changed, 54 insertions, 23 deletions

diff --git a/vertex_deformation.slang b/vertex_deformation.slang
index 3118394..aea1b5b 100644
--- a/vertex_deformation.slang
+++ b/vertex_deformation.slang
@@ -97,7 +97,8 @@ float2 project_x_onto_y(float2 x, float2 y) {
 // be orthogonal `s_cart`.
 [Differentiable]
 public float3 plane_to_tube(float3 xyz,
-    no_diff float3 p, no_diff float3 r_cart, no_diff float3 s_cart, no_diff float t) {
+    no_diff float3 p, no_diff float3 r_cart, no_diff float3 s_cart,
+    no_diff float t) {
   // Convert from cartesian to (r, s, r x s) space.
   // Ensure orthonormal basis vectors.
   // TODO remove normalize, do at higher level of stack.
@@ -157,39 +158,69 @@ public void plane_to_tube_undeform_normal(float3 xyz, inout float3 normal,
 // Maps a tube with circular cross section on the xz plane to a quad on the xy
 // plane.
 [Differentiable]
-public float3 tube_to_plane(float3 xyz, no_diff float t) {
-  float x0 = xyz.x;
-  float y0 = xyz.y;
-  float z0 = xyz.z;
+public float3 tube_to_plane(float3 xyz,
+    no_diff float3 p, no_diff float3 r_cart, no_diff float3 s_cart,
+    no_diff float t) {
+  // Convert from cartesian to (r, s, r x s) space.
+  // Ensure orthonormal basis vectors.
+  // TODO remove normalize, do at higher level of stack.
+  r_cart = normalize(r_cart);
+  s_cart = normalize(s_cart);
+  float3 rxs_cart = cross(s_cart, r_cart);
+  float3x3 to_rsrxs = float3x3(r_cart, s_cart, rxs_cart);
+  float3x3 to_cart  = inverse(to_rsrxs, determinant(to_rsrxs));
+
+  // Translate origin to `p` then change into (r, s, r x s) basis.
+  xyz = mul(to_rsrxs, xyz - p);
+
+  // Components in pivot basis: n (neutral axis), b (tangential), h (axial s).
+  float epsilon = 1e-4f;
+  float r = xyz.x;     // Lr
+  float s = xyz.y;     // Lv0 x Lr
+  float rxs = xyz.z;   // Lv0
+  float2 v0 = float2(r, rxs);
 
-  // Tube has 2m circumference.
-  // C = pi d = 2 pi r
-  // r = C / (2 pi) = 2 / 2 pi = 1 / pi
-  float ref_r = RCP_PI;
-  float theta = atan2(z0 - ref_r, x0 - ref_r) - PI;
-  float r = length(float2(x0, z0));
+  float theta0 = atan2(rxs, r);
+  float Lr  = length(v0);
+  float Lv0 = theta0 * Lr;
+  float phi = atan2(Lv0, Lr);
+  float rr = length(float2(Lr, Lv0));
 
-  x0 = x0 + r * cos(theta);
-  z0 = z0 - r * sin(theta);
+  // Blend between two vectors:
+  //  v0: vector of length ||Lr + Lv0|| at angle atan2(Lv0, Lr)
+  //  v1: vector of length ||Lr|| at angle ||Lv0|| / ||Lr||
+  // Interpolate in polar coordinates to make it wrap nicely.
+  float r0 = Lr;
+  float r1 = rr;
 
-  // Equivalent to (1 / (|t| + eps)) * sign(t)
-  float radius = 1.0f / (dabs(t) + 1e-4f) * sign(t);
+  float theta1 = phi;
 
-  float x = sin(theta / radius) * (radius * PI_RCP - z0);
-  float z = (1.0f - cos(theta / radius)) * (radius * PI_RCP) + cos(theta / radius) * z0;
+  // Interpolate polar coordinates.
+  float radius = dlerp(r0, r1, t);
+  float theta = dlerp(theta0, theta1, t);
 
-  return float3(x, y0, z);
+  // Map into (r, rxs) basis.
+  float2 nb_t = float2(cos(theta), sin(theta)) * radius;
+
+  // Un-project from (r, rxs) plane to full (r, s, rxs) basis.
+  xyz = float3(nb_t.x, s, nb_t.y);
+
+  // Map back to cartesian basis.
+  xyz = mul(to_cart, xyz) + p;
+
+  return xyz;
 }
 
 public void tube_to_plane_normal(inout float3 xyz, inout float3 normal,
-    inout float3 tangent, float t) {
-  R3R3_NORMALS(xyz, normal, tangent, tube_to_plane, t);
+    inout float3 tangent, float3 p, float3 r, float3 s, float t) {
+  R3R3_NORMALS(xyz, normal, tangent, tube_to_plane, p, r, s, t);
 }
 
-public void tube_to_plane_undeform_normal(float3 xyz, inout float3 normal,
-    inout float3 tangent, float t) {
+public void tube_to_plane_undeform_normal(
+    float3 xyz, inout float3 normal, inout float3 tangent,
+    float3 p, float3 r, float3 s, float t) {
   R3R3_DECLARE_BASIS_VECTORS(xyz);
-  R3R3_AUTODIFF_BASIS_VECTORS(tube_to_plane, t);
+  R3R3_AUTODIFF_BASIS_VECTORS(tube_to_plane, p, r, s, t);
   R3R3_UNDEFORM_NORMAL_AND_TANGENT(normal, tangent);
 }