From 4e8d518a61055c945834da17f74c6839c869ebf0 Mon Sep 17 00:00:00 2001
From: Alex Fink <000024@gmail.com>
Date: Sat, 14 Dec 2024 16:52:50 +0000
Subject: [PATCH] rel curvature computations

---
 renderer/pyscript-offline/public/main.py    |   8 +-
 renderer/pyscript-offline/public/relaxer.py | 111 +++++++++++++++++++-
 2 files changed, 113 insertions(+), 6 deletions(-)

diff --git a/renderer/pyscript-offline/public/main.py b/renderer/pyscript-offline/public/main.py
index fc2b3f2..c07bf19 100644
--- a/renderer/pyscript-offline/public/main.py
+++ b/renderer/pyscript-offline/public/main.py
@@ -16,7 +16,7 @@
 firstsg = relaxer.DifferentialGlyph.from_glyph(firstsg)
 firstsg.x = -2.
 firstsg.y = 0.
-firstsg.angle = -math.pi/2#math.pi/6#-math.pi/2
+firstsg.angle = math.pi/6#-math.pi/2
 text.add_glyph(firstsg)
 
 cat = dictionary.glyph_by_id("cat")
@@ -31,6 +31,6 @@
 
 canvas.render(text)
 
-document.body.append(f"Non-constant-velocity penalty: {str(relaxer.velocity_penalty(rel))}. ")
-cat.dx = 1.
-document.body.append(f"Derivative: {str(relaxer.deriv_velocity_penalty(rel))}. ")
+document.body.append(f"Max curvature: {str(relaxer.curvature_penalty(rel))}. ")
+firstsg.dangle = 1.
+document.body.append(f"Derivative: {str(relaxer.deriv_curvature_penalty(rel))}. ")
diff --git a/renderer/pyscript-offline/public/relaxer.py b/renderer/pyscript-offline/public/relaxer.py
index 93c631f..52201fa 100644
--- a/renderer/pyscript-offline/public/relaxer.py
+++ b/renderer/pyscript-offline/public/relaxer.py
@@ -93,11 +93,13 @@ def dpoly(rel):
   bezier = svgpathtools.CubicBezier(rel.bp0.dz, rel.bp0.handledz, rel.bp1.handledz, rel.bp1.dz)
   return bezier.poly()
 
+
+
 # The particular penalties below are not necessarily the ones we'll want to
 # go with in the end.
 
 def velocity_penalty(rel):
-  "Penalty score for this rel not going at constant velocity."
+  "Penalty for this rel not going at constant velocity."
   bezier = rel.svgpathtools_bezier()
   dd = bezier.poly().deriv().deriv()
   ddx = [float(dd[i].real) for i in range(0,2)]
@@ -108,7 +110,7 @@ def velocity_penalty(rel):
     ddx[0]*ddx[0]+ddy[0]*ddy[0]
 
 def deriv_velocity_penalty(rel):
-  """Derivative of penalty score for this rel not going at constant velocity.
+  """Derivative of penalty for this rel not going at constant velocity.
   
   The glyphs which rel is attached to should be DifferentialGlyph objects."""
   bezier = rel.svgpathtools_bezier()
@@ -124,3 +126,108 @@ def deriv_velocity_penalty(rel):
     ddx[0]*du_ddx[1]+ddy[0]*du_ddy[1] +\
     du_ddx[0]*ddx[1]+du_ddy[0]*ddy[1] +\
     2*(ddx[0]*du_ddx[0]+ddy[0]*du_ddy[0])
+
+def curvature_penalty(rel):
+  "Maximum unsigned curvature attained at any point of this rel."
+  return _curvature_penalty(rel, False)
+
+def deriv_curvature_penalty(rel):
+  "Derivative of maximum unsigned curvature attained."
+  return _curvature_penalty(rel, True)
+
+def _curvature_penalty(rel, differentiate):
+  """Computations for maximum curvature on this rel.
+  
+  Parameters:
+  differentiate: boolean, whether to return the derivative."""
+  # Let x, y be the coordinates of the cubic curve.
+  # Let ' denote derivative in the parameter of the curve (which we call t).
+  # Signed curvature is k = (x'y'' - y'x'')/((x'^2+y^'2)^(3/2)).
+  # k' = [ (x''y''+x'y'''-y''x''-y'x''')(x'^2+y^'2)^(3/2)
+  # - (x'y''-y'x'')*3/2(x'^2+y'^2)^(1/2)*2(x'x''+y'y'') ]
+  # / ((x^2+y^2)^3)
+  # = [ (x'y'''-y'x''')(x'^2+y^'2) - (x'y''-y'x'')*3(x'x''+y'y'') ]
+  # / ((x^2+y^2)^(5/2))).
+  # For a cubic curve the degree of the numerator is generically 5
+  # because of cancellations in the leading coefficient.
+  bezier = rel.svgpathtools_bezier()
+  d = bezier.poly().deriv()
+  dx = numpy.poly1d([float(c.real) for c in d.c])
+  dy = numpy.poly1d([float(c.imag) for c in d.c])
+  ddx = dx.deriv()
+  ddy = dy.deriv()
+  dddx = ddx.deriv()
+  dddy = ddy.deriv()
+  
+  def curvature(t):
+    "Evaluation of the signed curvature at curve parameter t."
+    dxt = dx(t)
+    dyt = dy(t)
+    return (dxt*ddy(t)-dyt*ddx(t)) / pow(dxt*dxt+dyt*dyt, 1.5)
+  
+  # numerator of k'
+  ndk = numpy.polysub(numpy.polymul(
+      numpy.polysub(numpy.polymul(dx, dddy), numpy.polymul(dy, dddx)),
+      numpy.polyadd(numpy.polymul(dx, dx), numpy.polymul(dy, dy)) ),
+    numpy.polymul(
+      numpy.polysub(numpy.polymul(dx, ddy), numpy.polymul(dy, ddx)),
+      numpy.polyadd(numpy.polymul(dx, ddx), numpy.polymul(dy, ddy)) )*3.
+    )
+  criticals = numpy.roots(ndk)
+  criticals = [0., 1.] + [t.real for t in criticals if t.imag==0. and t.real >= 0. and t.real <= 1.]
+  k_at_criticals = [curvature(t) for t in criticals]
+  abs_kmax, tmax, imax = max((abs(k_at_criticals[i]),criticals[i],i) for i in range(len(criticals)))
+  if not differentiate:
+    return abs_kmax
+  
+  sign_max = -1. if k_at_criticals[imax] < 0. else 1.
+  
+  du_d = dpoly(rel).deriv()
+  du_dx = numpy.poly1d([float(c.real) for c in du_d.c])
+  du_dy = numpy.poly1d([float(c.imag) for c in du_d.c])
+  du_ddx = du_dx.deriv()
+  du_ddy = du_dy.deriv()
+  
+  dxt = dx(tmax)
+  dyt = dy(tmax)
+  ddxt = ddx(tmax)
+  ddyt = ddy(tmax)
+  du_dxt = du_dx(tmax)
+  du_dyt = du_dy(tmax)
+  du_ddxt = du_ddx(tmax)
+  du_ddyt = du_ddy(tmax)
+  
+  # If the maximum curvature is attained at an endpoint,
+  # we just want the u-derivative of curvature there.
+  # The evaluations du_dxt, etc. above are correct for this purpose.
+  #
+  # But if the maximum curvature is attained at an interior critical point tmax,
+  # work out a=dt/du on the locus of critical points by implicit differentiation,
+  # and then differentiate the max curvature not in direction du but du+a dt.
+  if imax >= 2:
+    dddxt = dddx(tmax)
+    dddyt = dddy(tmax)
+    du_dddx = du_ddx.deriv()
+    du_dddy = du_ddy.deriv()
+    du_dddxt = du_dddx(tmax)
+    du_dddyt = du_dddy(tmax)
+    # Consider p = ndk(t) + u du_ndk(t) and find (dp/du)/(dp/dt)
+    # at (t, u) = (tmax, 0).
+    dp = ndk.deriv()(tmax) # ndk.deriv() = (dp/dt)|u=0
+    # And below is dp/du, which is independent of u, at tmax.
+    # recall, ndk = (x'y'''-y'x''')(x'^2+y^'2) - (x'y''-y'x'')*3(x'x''+y'y'')
+    du_p = (du_dxt*dddyt+dxt*du_dddyt-du_dyt*dddxt-dyt*du_dddxt)*(dxt*dxt+dyt*dyt) + \
+      2*(dxt*dddyt-dyt*dddxt)*(du_dxt*dxt+du_dyt*dyt) - \
+      3*(du_dxt*ddyt+dxt*du_ddyt-du_dyt*ddxt-dyt*du_ddxt)*(dxt*ddxt+dyt*ddyt) - \
+      3*(dxt*ddyt-dyt*ddxt)*(du_dxt*ddxt+dxt*du_ddxt+du_dyt*ddyt+dyt*du_ddyt)
+    a = du_p/dp
+    # Replace the four du_ variables used below by derivatives in direction du+a dt.
+    du_dxt += a*ddxt
+    du_dyt += a*ddyt
+    du_ddxt += a*dddxt
+    du_ddyt += a*dddyt
+  
+  du_k = ((du_dxt*ddyt+dxt*du_ddyt-du_dyt*ddxt-dyt*du_ddxt)*(dxt*dxt+dyt*dyt) - \
+    3*(dxt*ddyt-dyt*ddxt)*(du_dxt*dxt+du_dyt*dyt)) / \
+    pow(dxt*dxt+dyt*dyt, 2.5)
+  return sign_max * du_k