Finding distance point to bezier curve (illustration purpose)

bezier.py

import torch as th


def bezier( t, P1, P2, P3, P4 ):
    return P1*(1-t**3) + P2*3*(1-t**2)*t + P3*3*(1-t)*t*t + P4*t**3

def dist2( a, b):
    return th.sum( (a-b)**2 )

P = th.tensor((1,2))

th.manual_seed(42)

C1 = th.randn((2,))
C2 = th.randn((2,))
C3 = th.randn((2,))
C4 = th.randn((2,))

Pb = bezier( 0.5, C1,C2,C3,C4)
print(Pb)


ut0 = th.nn.Parameter(th.randn((1,))) #unconstrained t0

#Computing P0
opt = th.optim.Adam([ut0],lr=0.1) #You should use newton method here instead
for i in range(200): #with newton we expect to need less than 3 iterations
    opt.zero_grad()
    t0 = th.sigmoid(ut0)#we constrain u to be in [0,1]
    P0 = bezier( t0, C1,C2,C3,C4 )
    err = dist2( P, P0)
    err.backward()
    opt.step()

t0 = t0.detach()

print("ut0.grad") #should be close to 0 at optimum
print(ut0.grad) 

d0 = th.sqrt( err ).detach()

print("d0")
print(d0)
print("t0:")
print( t0 )
print("P0:")
print( P0.detach().numpy() )

ut1 = th.nn.Parameter(-3.0*th.ones((1,))) #we try to start close to the first control point 

P1ur = th.nn.Parameter(th.randn((1,)) ) #the unconstrained radius
P1normed = th.nn.Parameter(th.randn((2,)) ) #the overparametrized direction of P1


#Computing P1 
opt = th.optim.Adam([ut1,P1ur,P1normed],lr=0.1) #You should use newton method here instead
for i in range(200):
    opt.zero_grad()
    t1 = th.sigmoid(ut1)*t0 #we constrain u to be in [0,t0]
    #print("t1 : ")
    #print(t1)
    #we should start in a feasible region where d(P1) < d0
    #A feasible starting point can be computed by the intersection with the bezier curve
    #with the circle centered on P and radius d0 - epsilon
    #The circle in question is already tangent to the curve at t0
    #I don't think bezier curve can cross a circle tangentially (aka P0 be a saddle point) but I may be wrong
    
    #This intersection can be solve by solving a quartic polynomial or
    #This intersection can be approximated by slicing the curve into piecewise segment
    
    #P1 is constrained to be inside the disk centered on P, and of radius < d0
    P1 = P + d0* th.sigmoid( P1ur) * P1normed / th.sqrt(th.sum(P1normed**2))
    P1b = bezier( t1, C1,C2,C3,C4 )

    slackerror = dist2( P1,P1b)
    #a constraint solver will use lagrange multiplier instead of the 100
    err = dist2( P, P1) + 100* slackerror
    #for the solution to be valid in the feasible set we must have dist2(P1,P1b) = 0

    err.backward()
    opt.step()

d1 = th.sqrt(dist2(P,P1b))

print("slackerror : ")
print(slackerror)

print("t1")
print(t1) #if we are close to t0

print("ut1.grad") #should be 0 at optimum
print(ut1.grad)

print("d1")
print(d1) 

print("P1:")
print( P1.detach().numpy() )


#Computing P2 : similar than P1 but on [t0,1]
ut2 = th.nn.Parameter(3.0*th.ones((1,))) #we try to start close to the second control point 
#we can use constrain on the radius min( d0,d1 ) around P

#and so on
#d2 = ...

#d = min( d0,d1,d2)