Continue fucking with trochoid normals

Seems like the determinant of the jacobian is ~0, meaning it's not
invertible. That means the normals look fucked up.

3 files changed, 153 insertions, 87 deletions

diff --git a/math.cginc b/math.cginc
index 758ff52..5b33edf 100644
--- a/math.cginc
+++ b/math.cginc
@@ -247,5 +247,50 @@ bool solveQuadratic(float a, float b, float c, out float x0, out float x1)
   return true;
 }
 
+float determinant(float3x3 m)
+{
+  return m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
+            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
+            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
+}
+
+float3x3 invert(float3x3 m)
+{
+  float det = determinant(m);
+
+  float3x3 adj;
+  adj[0][0] =  (m[1][1] * m[2][2] - m[1][2] * m[2][1]);
+  adj[0][1] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]);
+  adj[0][2] =  (m[0][1] * m[1][2] - m[0][2] * m[1][1]);
+  
+  adj[1][0] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]);
+  adj[1][1] =  (m[0][0] * m[2][2] - m[0][2] * m[2][0]);
+  adj[1][2] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]);
+  
+  adj[2][0] =  (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
+  adj[2][1] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]);
+  adj[2][2] =  (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
+
+  return adj * (1.0 / det);
+}
+
+// Return largest number which divides into both 'a' and 'b'.
+// Uses the Euclidean algorithm: repeatedly divide 'b' by 'a'.
+uint gcd(uint a, uint b)
+{
+  uint tmp = a * b;
+  #define GCD_MAX_ITER 10
+  for (uint i = 0; i < GCD_MAX_ITER; i++) {
+    tmp = b;
+    b = a % b;
+    a = tmp;
+    if (a == b) {
+      return a;
+    }
+  }
+  return 1;
+  
+}
+
 #endif  // __MATH_INC
 
diff --git a/tooner_lighting.cginc b/tooner_lighting.cginc
index adf7890..69b19a4 100644
--- a/tooner_lighting.cginc
+++ b/tooner_lighting.cginc
@@ -138,7 +138,40 @@ v2f vert(appdata v)
 #if defined(_TROCHOID)
   {
     o.objPos_pre_trochoid = v.vertex.xyz;
+//#define TROCHOID_DECOMPOSE
+#if defined(TROCHOID_DECOMPOSE)
     v.vertex.xyz = cyl2_to_troch_map(cyl_to_cyl2_map(cart_to_cyl_map(v.vertex.xyz)));
+#else
+    v.vertex.xyz = cart_to_troch_map(v.vertex.xyz);
+#endif
+  }
+#endif
+
+#if defined(_TROCHOID)
+  {
+#if defined(TROCHOID_DECOMPOSE)
+    // Let h(v) be the trochoid of the cartesian coordinates v.
+    // We evaluate h(v) by first mapping it to cylindrical coordinates, then applying a trochoid function defined on those coordinates:
+    //   h(v) = h_cyl(g(v))
+    // g(v) maps v to cylindrical coordinates.
+    // h_cyl(v) evaluates the trochoid function on cylindrical coordinates.
+    // We want to compute h'(v), i.e. its Jacobian.
+    // By the chain rule:
+    //   h'(v) = h_cyl'(g(v)) * g'(v)
+    // The reality is a little more complex: we also apply a distortion to the cylindrical coordinates:
+    //   h(v) = h_cyl(f(g(v)))
+    // where f(v) is that distortion.
+    // Again, by the chain rule:
+    //   h'(v) = h_cyl'(f(g(v))) * f'(g(v)) * g'(v)
+    float3x3 j0 = cart_to_cyl_jacobian(o.objPos_pre_trochoid);
+    float3x3 j1 = cyl_to_cyl2_jacobian(cart_to_cyl_map(o.objPos_pre_trochoid));
+    float3x3 j2 = cyl2_to_troch_jacobian(cyl_to_cyl2_map(cart_to_cyl_map(o.objPos_pre_trochoid)));
+    float3x3 vector_mover = transpose(invert(mul(mul(j2, j1), j0)));
+#else
+    float3x3 vector_mover = transpose(invert(cart_to_troch_jacobian(o.objPos_pre_trochoid)));
+#endif
+    v.normal = mul(vector_mover, v.normal);
+    v.tangent.xyz = mul(vector_mover, v.tangent.xyz);
   }
 #endif
 #if defined(_FACE_ME_WORLD_Y)
@@ -1595,28 +1628,6 @@ float4 effect(inout v2f i, out float depth)
   i.tangent.xyz = normalize(i.tangent.xyz - i.normal * dot(i.tangent.xyz, i.normal));
   //i.tangent.xyz = normalize(i.tangent.xyz);
 
-#if defined(_TROCHOID)
-  {
-    // Let h(v) be the trochoid of the cartesian coordinates v.
-    // We evaluate h(v) by first mapping it to cylindrical coordinates, then applying a trochoid function defined on those coordinates:
-    //   h(v) = h_cyl(g(v))
-    // g(v) maps v to cylindrical coordinates.
-    // h_cyl(v) evaluates the trochoid function on cylindrical coordinates.
-    // We want to compute h'(v), i.e. its Jacobian.
-    // By the chain rule:
-    //   h'(v) = h_cyl'(g(v)) * g'(v)
-    // The reality is a little more complex: we also apply a distortion to the cylindrical coordinates:
-    //   h(v) = h_cyl(f(g(v)))
-    // where f(v) is that distortion.
-    // Again, by the chain rule:
-    //   h'(v) = h_cyl'(f(g(v))) * f'(g(v)) * g'(v)
-    float3x3 j0 = cart_to_cyl_jacobian(i.objPos_pre_trochoid);
-    float3x3 j1 = cyl_to_cyl2_jacobian(cart_to_cyl_map(i.objPos_pre_trochoid));
-    float3x3 j2 = cyl2_to_troch_jacobian(cyl_to_cyl2_map(cart_to_cyl_map(i.objPos_pre_trochoid)));
-    i.normal = normalize(mul(transpose(invert(mul(j0, mul(j1, j2)))), i.normal));
-  }
-#endif
-
 #if defined(_UVSCROLL)
   float2 orig_uv = i.uv0;
   float uv_scroll_mask = round(_UVScroll_Mask.SampleBias(linear_repeat_s, i.uv0, _Global_Sample_Bias));
@@ -2896,6 +2907,7 @@ fixed4 frag(v2f i
   UNITY_APPLY_DITHER_CROSSFADE(i.pos.xy);
   UNITY_SETUP_INSTANCE_ID(i);
   UNITY_SETUP_STEREO_EYE_INDEX_POST_VERTEX(i);
+
   return effect(i, depth);
 }
 
diff --git a/trochoid_math.cginc b/trochoid_math.cginc
index dd6d69c..3a5107e 100644
--- a/trochoid_math.cginc
+++ b/trochoid_math.cginc
@@ -10,24 +10,30 @@
 #define TROCH_Z_THETA_SCALE 0.0002
 #define TROCH_EPSILON 1e-5
 
-float3 cyl2_to_troch_map(float3 pos)
+float3 cyl2_to_troch_map(float3 v)
 {
-  float theta = pos.x;
-
-  float R = _Trochoid_R;
-  float r = _Trochoid_r;
-  float d = _Trochoid_d;
-
-  float x = ((R - r) * cos(theta * R) + d * cos((R - r) * theta * R / r)) * pos.y * TROCH_POSITION_SCALE;
-  float y = ((R - r) * sin(theta * R) - d * sin((R - r) * theta * R / r)) * pos.y * TROCH_POSITION_SCALE;
-  float z = (pos.z + cos(theta * R * 5) * TROCH_POSITION_SCALE + theta * R * TROCH_Z_THETA_SCALE) * pos.y;
+  const float R = _Trochoid_R;
+  const float r = _Trochoid_r;
+  const float d = _Trochoid_d;
+  const float rrrr = (R - r) * R / r;
+
+  float x =
+      cos(v.x * R) * v.y * (R - r) * TROCH_POSITION_SCALE +
+      cos(v.x * rrrr) * v.y * d * TROCH_POSITION_SCALE;
+  float y =
+      sin(v.x * R) * v.y * (R - r) * TROCH_POSITION_SCALE -
+      sin(v.x * rrrr) * v.y * d * TROCH_POSITION_SCALE;
+  float z =
+      (v.x * v.y * R * TROCH_Z_THETA_SCALE +
+      cos(v.x * R * 5) * v.y * TROCH_POSITION_SCALE +
+      v.y * v.z);
 
   return float3(x, y, z);
 }
 
 float3 cyl_to_cyl2_map(float3 v)
 {
-  return float3(v.x, pow(v.y, _Trochoid_Radius_Power) * _Trochoid_Radius_Scale, v.z * _Trochoid_Height_Scale);
+  return float3(v.x * .5, pow(v.y, _Trochoid_Radius_Power) * _Trochoid_Radius_Scale, v.z * _Trochoid_Height_Scale);
 }
 
 float3 cart_to_cyl_map(float3 v)
@@ -35,33 +41,60 @@ float3 cart_to_cyl_map(float3 v)
   return float3(atan2(v.y, v.x), length(v.xy), v.z);
 }
 
-// Compute partial derivatives of trochoid function with respect to cylindrical coordinates
-float3x3 cyl2_to_troch_jacobian(float3 pos)
+float3 cart_to_troch_map(float3 v)
 {
-  float theta = pos.x;
+  return cyl2_to_troch_map(cyl_to_cyl2_map(cart_to_cyl_map(v)));
+}
 
-  float R = _Trochoid_R;
-  float r = _Trochoid_r;
-  float d = _Trochoid_d;
+float3x3 cart_to_troch_jacobian(float3 v)
+{
+  float epsilon = 1e-5;
+  float3 df_dx = (cart_to_troch_map(v + float3(epsilon, 0, 0)) - cart_to_troch_map(v - float3(epsilon, 0, 0))) / (2 * epsilon);
+  float3 df_dy = (cart_to_troch_map(v + float3(0, epsilon, 0)) - cart_to_troch_map(v - float3(0, epsilon, 0))) / (2 * epsilon);
+  float3 df_dz = (cart_to_troch_map(v + float3(0, 0, epsilon)) - cart_to_troch_map(v - float3(0, 0, epsilon))) / (2 * epsilon);
+  return transpose(float3x3(df_dx, df_dy, df_dz));
+}
+
+// Compute partial derivatives of trochoid function with respect to cylindrical coordinates
+float3x3 cyl2_to_troch_jacobian(float3 v)
+{
+  const float R = _Trochoid_R;
+  const float r = _Trochoid_r;
+  const float d = _Trochoid_d;
+  const float rrrr = (R - r) * R / r;
 
 #if 1
   float3 df_dt = float3(
-    ((R * (r - R) *  (d * sin(R * theta * (R - r) / r) + r * sin(R * theta))) / r) * pos.y * TROCH_POSITION_SCALE,
-    ((R * (r - R) * (-d * cos(R * theta * (R - r) / r) + r * cos(R * theta))) / r) * pos.y * TROCH_POSITION_SCALE,
-    (-R * 5 * sin(theta * R * 5) * TROCH_POSITION_SCALE + R * TROCH_Z_THETA_SCALE) * pos.y);
+    -R * sin(v.x * R) * v.y * (R - r) * TROCH_POSITION_SCALE +
+    -rrrr * sin(v.x * rrrr) * v.y * d * TROCH_POSITION_SCALE,
+
+    R * cos(v.x * R) * v.y * (R - r) * TROCH_POSITION_SCALE -
+    rrrr * cos(v.x * rrrr) * v.y * d * TROCH_POSITION_SCALE,
+
+    v.y * R * TROCH_Z_THETA_SCALE -
+    R * 5 * sin(v.x * R * 5) * v.y * TROCH_POSITION_SCALE);
+#else
+  float3 df_dt = (cyl2_to_troch_map(v + float3(TROCH_EPSILON, 0, 0)) - cyl2_to_troch_map(v - float3(TROCH_EPSILON, 0, 0))) / (2 * TROCH_EPSILON);
+#endif
+
+#if 1
   float3 df_dr = float3(
-    ((R - r) * cos(theta * R) + d * cos((R - r) * theta * R / r)) * TROCH_POSITION_SCALE,
-    ((R - r) * sin(theta * R) - d * sin((R - r) * theta * R / r)) * TROCH_POSITION_SCALE,
-    (pos.z + cos(theta * R * 5) * TROCH_POSITION_SCALE + theta * R * TROCH_Z_THETA_SCALE));
+    ((R - r) * cos(v.x * R) + d * cos((R - r) * v.x * R / r)) * TROCH_POSITION_SCALE,
+    ((R - r) * sin(v.x * R) - d * sin((R - r) * v.x * R / r)) * TROCH_POSITION_SCALE,
+    cos(v.x * R * 5) * TROCH_POSITION_SCALE);
+#else
+  float3 df_dr = (cyl2_to_troch_map(v + float3(0, TROCH_EPSILON, 0)) - cyl2_to_troch_map(v - float3(0, TROCH_EPSILON, 0))) / (2 * TROCH_EPSILON);
+#endif
+
+#if 1
   float3 df_dz = float3(
     0,
     0,
-    pos.y);
+    v.y);
 #else
-  float3 df_dt = (cyl2_to_troch_map(pos + float3(TROCH_EPSILON, 0, 0)) - cyl2_to_troch_map(pos - float3(TROCH_EPSILON, 0, 0))) / (2 * TROCH_EPSILON);
-  float3 df_dr = (cyl2_to_troch_map(pos + float3(0, TROCH_EPSILON, 0)) - cyl2_to_troch_map(pos - float3(0, TROCH_EPSILON, 0))) / (2 * TROCH_EPSILON);
-  float3 df_dz = (cyl2_to_troch_map(pos + float3(0, 0, TROCH_EPSILON)) - cyl2_to_troch_map(pos - float3(0, 0, TROCH_EPSILON))) / (2 * TROCH_EPSILON);
+  float3 df_dz = (cyl2_to_troch_map(v + float3(0, 0, TROCH_EPSILON)) - cyl2_to_troch_map(v - float3(0, 0, TROCH_EPSILON))) / (2 * TROCH_EPSILON);
 #endif
+
   float3x3 jacobian_cyl;
   jacobian_cyl[0] = df_dt;
   jacobian_cyl[1] = df_dr;
@@ -69,17 +102,17 @@ float3x3 cyl2_to_troch_jacobian(float3 pos)
   return transpose(jacobian_cyl);
 }
 
-float3x3 cyl_to_cyl2_jacobian(float3 pos)
+float3x3 cyl_to_cyl2_jacobian(float3 v)
 {
-  // f(x, y, z) = <x, (y^_Trochoid_Radius_Power) * _Trochoid_Radius_Scale, z * _Trochoid_Height_Scale>
+  // f(x, y, z) = <x, (y^_Trochoid_Radiu1_Power) * _Trochoid_Radius_Scale, z * _Trochoid_Height_Scale>
 #if 1
   float3 df_dx = float3(1, 0, 0);
-  float3 df_dy = float3(0, _Trochoid_Radius_Power * pow(pos.y, _Trochoid_Radius_Power - 1) * _Trochoid_Radius_Scale, 0);
+  float3 df_dy = float3(0, _Trochoid_Radius_Power * pow(v.y, _Trochoid_Radius_Power - 1) * _Trochoid_Radius_Scale, 0);
   float3 df_dz = float3(0, 0, _Trochoid_Height_Scale);
 #else
-  float3 df_dx = (cyl_to_cyl2_map(pos + float3(TROCH_EPSILON, 0, 0)) - cyl_to_cyl2_map(pos - float3(TROCH_EPSILON, 0, 0))) / (2 * TROCH_EPSILON);
-  float3 df_dy = (cyl_to_cyl2_map(pos + float3(0, TROCH_EPSILON, 0)) - cyl_to_cyl2_map(pos - float3(0, TROCH_EPSILON, 0))) / (2 * TROCH_EPSILON);
-  float3 df_dz = (cyl_to_cyl2_map(pos + float3(0, 0, TROCH_EPSILON)) - cyl_to_cyl2_map(pos - float3(0, 0, TROCH_EPSILON))) / (2 * TROCH_EPSILON);
+  float3 df_dx = (cyl_to_cyl2_map(v + float3(TROCH_EPSILON, 0, 0)) - cyl_to_cyl2_map(v - float3(TROCH_EPSILON, 0, 0))) / (2 * TROCH_EPSILON);
+  float3 df_dy = (cyl_to_cyl2_map(v + float3(0, TROCH_EPSILON, 0)) - cyl_to_cyl2_map(v - float3(0, TROCH_EPSILON, 0))) / (2 * TROCH_EPSILON);
+  float3 df_dz = (cyl_to_cyl2_map(v + float3(0, 0, TROCH_EPSILON)) - cyl_to_cyl2_map(v - float3(0, 0, TROCH_EPSILON))) / (2 * TROCH_EPSILON);
 #endif
   float3x3 jacobian_cyl_to_cyl2;
   jacobian_cyl_to_cyl2[0] = df_dx;
@@ -89,28 +122,28 @@ float3x3 cyl_to_cyl2_jacobian(float3 pos)
 }
 
 // Compute partial derivatives of transform from cartesian to cylindrical coordinates
-float3x3 cart_to_cyl_jacobian(float3 pos)
+float3x3 cart_to_cyl_jacobian(float3 v)
 {
   // Compute partial derivatives of transform from cartesian to cylindrical coordinates
   // return float3(atan2(v.y, v.x), length(v.xy), v.z);
   // theta = atan2(y, x)
 #if 1
   float3 dtheta_dcart = float3(
-    -pos.y / dot(pos.xy, pos.xy),
-    pos.x / dot(pos.xy, pos.xy),
+    -v.y / dot(v.xy, v.xy),
+    v.x / dot(v.xy, v.xy),
     0);
 #else
-  float3 dtheta_dcart = (cart_to_cyl_map(pos + float3(TROCH_EPSILON, 0, 0)) - cart_to_cyl_map(pos - float3(TROCH_EPSILON, 0, 0))) / (2 * TROCH_EPSILON);
+  float3 dtheta_dcart = (cart_to_cyl_map(v + float3(TROCH_EPSILON, 0, 0)) - cart_to_cyl_map(v - float3(TROCH_EPSILON, 0, 0))) / (2 * TROCH_EPSILON);
 #endif
 
   // radius = (x^2 + y^2)^(1/2)
 #if 1
   float3 dr_dcart = float3(
-    pos.x / sqrt(pos.x * pos.x + pos.y * pos.y),
-    pos.y / sqrt(pos.x * pos.x + pos.y * pos.y),
+    v.x / sqrt(v.x * v.x + v.y * v.y),
+    v.y / sqrt(v.x * v.x + v.y * v.y),
     0);
 #else
-  float3 dr_dcart = (cart_to_cyl_map(pos + float3(0, TROCH_EPSILON, 0)) - cart_to_cyl_map(pos - float3(0, TROCH_EPSILON, 0))) / (2 * TROCH_EPSILON);
+  float3 dr_dcart = (cart_to_cyl_map(v + float3(0, TROCH_EPSILON, 0)) - cart_to_cyl_map(v - float3(0, TROCH_EPSILON, 0))) / (2 * TROCH_EPSILON);
 #endif
 
 #if 1
@@ -119,7 +152,7 @@ float3x3 cart_to_cyl_jacobian(float3 pos)
     0,
     1);
 #else
-  float3 dz_dcart = (cart_to_cyl_map(pos + float3(0, 0, TROCH_EPSILON)) - cart_to_cyl_map(pos - float3(0, 0, TROCH_EPSILON))) / (2 * TROCH_EPSILON);
+  float3 dz_dcart = (cart_to_cyl_map(v + float3(0, 0, TROCH_EPSILON)) - cart_to_cyl_map(v - float3(0, 0, TROCH_EPSILON))) / (2 * TROCH_EPSILON);
 #endif
   float3x3 jacobian_cart_to_cyl;
   jacobian_cart_to_cyl[0] = dtheta_dcart;
@@ -128,30 +161,6 @@ float3x3 cart_to_cyl_jacobian(float3 pos)
   return transpose(jacobian_cart_to_cyl);
 }
 
-float3x3 invert(float3x3 m)
-{
-  float det = m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
-            - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
-            + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
-
-  float3x3 adj;
-  adj[0][0] =  (m[1][1] * m[2][2] - m[1][2] * m[2][1]);
-  adj[0][1] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]);
-  adj[0][2] =  (m[0][1] * m[1][2] - m[0][2] * m[1][1]);
-  
-  adj[1][0] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]);
-  adj[1][1] =  (m[0][0] * m[2][2] - m[0][2] * m[2][0]);
-  adj[1][2] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]);
-  
-  adj[2][0] =  (m[1][0] * m[2][1] - m[1][1] * m[2][0]);
-  adj[2][1] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]);
-  adj[2][2] =  (m[0][0] * m[1][1] - m[0][1] * m[1][0]);
-
-  return adj * (1.0 / det);
-}
-
 #endif  // _TROCHOID
 
 #endif  // __TROCHOID_MATH
-
-