This problem is just begging for Newton approximation like recursions.
How to find distance of arc up to precision required?
> First iteration: Find out the distance between both endpoints. Create
> new points at 1/4 of a distance from 1/4 first control point to 3/4 of
> second control point, plus 3/4 of a distance from 1st point, plus 1/2
> the distance. Measure length of a polyline comprising these five
> points. Is second distance minus first distance smaller than our
> tolerance?
> - yes => return distance
> - no => recursion Other iterations: Create points at 1/8, 3/8, 5/8 and 7/8 (or 1/16th, 3/16ths etc., and so on in subsequent iterations) Is
> new distance - previous distance < tolerance? Yes/no, return
> distance/recursion.
In this way not only can you tell what the distance is, but also what is maximum imprecision in your calculation, plus you never calculate more than what is absolutely necessary.
To find the point at t I would strongly suggest using pure Newton's approximation, as this is what it calculates, and within defined toleration.
Pyhon code, created and tested in Blender 2.83 LTS.
```python
from math import sqrt
def lengthOfCurve(points, H1, H2, toleration):
A = points[0]
E = points[-1]
if len(points) == 2: # first round
Ba = [(H1[0] - A[0])*0.25 + A[0], (H1[1] - A[1])*0.25 + A[1], (H1[2] - A[2])*0.25 + A[2]]
Bb = [(E[0] - H2[0])*0.25 + H2[0], (E[1] - H2[1])*0.25 + H2[1], (E[2] - H2[2])*0.25 + H2[2]]
B = [(Bb[0] - Ba[0])*0.25 + Ba[0], (Bb[1] - Ba[1])*0.25 + Ba[1], (Bb[2] - Ba[2])*0.25 + Ba[2]]
Ca = [(H1[0] - A[0])*0.5 + A[0], (H1[1] - A[1])*0.5 + A[1], (H1[2] - A[2])*0.5 + A[2]]
Cb = [(E[0] - H2[0])*0.5 + H2[0], (E[1] - H2[1])*0.5 + H2[1], (E[2] - H2[2])*0.5 + H2[2]]
C = [(Cb[0] - Ca[0])*0.5 + Ca[0], (Cb[1] - Ca[1])*0.5 + Ca[1], (Cb[2] - Ca[2])*0.5 + Ca[2]]
Da = [(H1[0] - A[0])*0.75 + A[0], (H1[1] - A[1])*0.75 + A[1], (H1[2] - A[2])*0.75 + A[2]]
Db = [(E[0] - H2[0])*0.75 + H2[0], (E[1] - H2[1])*0.75 + H2[1], (E[2] - H2[2])*0.75 + H2[2]]
D = [(Db[0] - Da[0])*0.75 + Da[0], (Db[1] - Da[1])*0.75 + Da[1], (Db[2] - Da[2])*0.75 + Da[2]]
points = [A, B, C, D, E]
originalLength = sizeOfVector(A, E)
newLength = sizeOfVector(A, B) + sizeOfVector(B, C) + sizeOfVector(C, D) + sizeOfVector(D, E)
if newLength - toleration < originalLength:
return(round(newLength/toleration)*toleration)
return(lengthOfCurve(points, H1, H2, toleration))
else: # following rounds
numberOfSegments = len(points) - 1
points2 = []
i = numberOfSegments
while i > 0:
t = 2*i - 1
pointA = [(H1[0] - A[0])*(t/(2*numberOfSegments)) + A[0], (H1[1] - A[1])*(t/(2*numberOfSegments)) + A[1], (H1[2] - A[2])*(t/(2*numberOfSegments)) + A[2]]
pointB = [(E[0] - H2[0])*(t/(2*numberOfSegments)) + H2[0], (E[1] - H2[1])*(t/(2*numberOfSegments)) + H2[1], (E[2] - H2[2])*(t/(2*numberOfSegments)) + H2[2]]
point = [(pointB[0] - pointA[0])*(t/(2*numberOfSegments)) + pointA[0], (pointB[1] - pointA[1])*(t/(2*numberOfSegments)) + pointA[1], (pointB[2] - pointA[2])*(t/(2*numberOfSegments)) + pointA[2]]
points2.insert(0, point) # I am prepending, because I am advancing bacwards
i -= 1
newPoints = []
i = 0
while i < numberOfSegments: # join points and points2 using zipper
newPoints.append(points[i])
newPoints.append(points2[i])
i += 1
newPoints.append(points[numberOfSegments])
newLength = 0
i = numberOfSegments
while i > 0:
newLength += sizeOfVector(newPoints[2*i], newPoints[2*i - 1])
newLength += sizeOfVector(newPoints[2*i - 1], newPoints[2*i - 2])
i -= 1
originalLength = 0
i = int(numberOfSegments/2)
while i > 0:
originalLength += sizeOfVector(points[2*i], points[2*i - 1])
originalLength += sizeOfVector(points[2*i - 1], points[2*i - 2])
i -= 1
if newLength - toleration < originalLength:
return(round(newLength/toleration)*toleration)
return(lengthOfCurve(newPoints, H1, H2, toleration))
def sizeOfVector(P, Q):
return(sqrt((P[0] - Q[0])*(P[0] - Q[0]) + (P[1] - Q[1])*(P[1] - Q[1]) + (P[2] - Q[2])*(P[2] - Q[2])))
point1 = [0, 0, 0]
point2 = [1, 0, 0]
handle1 = [0, 1, 0]
handle2 = [1, 1, 0]
toleration = 0.01
print(lengthOfCurve([point1, point2], handle1, handle2, toleration))
```
EDIT:
Ok, I got the t working as well, within toleration, as expected. It was far from simple approximation. Let me know if there is dire need for it in English. I did all of this in a single day and I feel a bit exhausted to translate the second chunk from Czech at the moment.
```python
from math import sqrt
def dylkaKrivky(body, H1, H2, tolerance):
A = body[0]
E = body[-1]
if len(body) == 2: # prvni kolo
Ba = [(H1[0] - A[0])*0.25 + A[0], (H1[1] - A[1])*0.25 + A[1], (H1[2] - A[2])*0.25 + A[2]]
Bb = [(E[0] - H2[0])*0.25 + H2[0], (E[1] - H2[1])*0.25 + H2[1], (E[2] - H2[2])*0.25 + H2[2]]
B = [(Bb[0] - Ba[0])*0.25 + Ba[0], (Bb[1] - Ba[1])*0.25 + Ba[1], (Bb[2] - Ba[2])*0.25 + Ba[2]]
Ca = [(H1[0] - A[0])*0.5 + A[0], (H1[1] - A[1])*0.5 + A[1], (H1[2] - A[2])*0.5 + A[2]]
Cb = [(E[0] - H2[0])*0.5 + H2[0], (E[1] - H2[1])*0.5 + H2[1], (E[2] - H2[2])*0.5 + H2[2]]
C = [(Cb[0] - Ca[0])*0.5 + Ca[0], (Cb[1] - Ca[1])*0.5 + Ca[1], (Cb[2] - Ca[2])*0.5 + Ca[2]]
Da = [(H1[0] - A[0])*0.75 + A[0], (H1[1] - A[1])*0.75 + A[1], (H1[2] - A[2])*0.75 + A[2]]
Db = [(E[0] - H2[0])*0.75 + H2[0], (E[1] - H2[1])*0.75 + H2[1], (E[2] - H2[2])*0.75 + H2[2]]
D = [(Db[0] - Da[0])*0.75 + Da[0], (Db[1] - Da[1])*0.75 + Da[1], (Db[2] - Da[2])*0.75 + Da[2]]
body = [A, B, C, D, E]
dylkaPuvodni = velikostVektoru(A, E)
dylkaNova = velikostVektoru(A, B) + velikostVektoru(B, C) + velikostVektoru(C, D) + velikostVektoru(D, E)
if dylkaNova - tolerance < dylkaPuvodni:
return(round(dylkaNova/tolerance)*tolerance)
return(dylkaKrivky(body, H1, H2, tolerance))
else: # nasledujici kola
usecek = len(body) - 1
body2 = []
i = usecek
while i > 0:
t = 2*i - 1
bodA = [(H1[0] - A[0])*(t/(2*usecek)) + A[0], (H1[1] - A[1])*(t/(2*usecek)) + A[1], (H1[2] - A[2])*(t/(2*usecek)) + A[2]]
bodB = [(E[0] - H2[0])*(t/(2*usecek)) + H2[0], (E[1] - H2[1])*(t/(2*usecek)) + H2[1], (E[2] - H2[2])*(t/(2*usecek)) + H2[2]]
bod = [(bodB[0] - bodA[0])*(t/(2*usecek)) + bodA[0], (bodB[1] - bodA[1])*(t/(2*usecek)) + bodA[1], (bodB[2] - bodA[2])*(t/(2*usecek)) + bodA[2]]
body2.insert(0, bod) # predpojuju, protoze postupju odzadu
i -= 1
noveBody = []
i = 0
while i < usecek: # spojit body a body2 zipem
noveBody.append(body[i])
noveBody.append(body2[i])
i += 1
noveBody.append(body[usecek])
dylkaNova = 0
i = usecek
while i > 0:
dylkaNova += velikostVektoru(noveBody[2*i], noveBody[2*i - 1])
dylkaNova += velikostVektoru(noveBody[2*i - 1], noveBody[2*i - 2])
i -= 1
dylkaPuvodni = 0
i = int(usecek/2)
while i > 0:
dylkaPuvodni += velikostVektoru(body[2*i], body[2*i - 1])
dylkaPuvodni += velikostVektoru(body[2*i - 1], body[2*i - 2])
i -= 1
if dylkaNova - tolerance < dylkaPuvodni:
return(round(dylkaNova/tolerance)*tolerance)
return(dylkaKrivky(noveBody, H1, H2, tolerance))
def velikostVektoru(P, Q):
return(sqrt((P[0] - Q[0])*(P[0] - Q[0]) + (P[1] - Q[1])*(P[1] - Q[1]) + (P[2] - Q[2])*(P[2] - Q[2])))
bod1 = [0, 0, 0]
bod2 = [1, 0, 0]
madlo1 = [0, 1, 0]
madlo2 = [1, 1, 0]
tolerance = 0.01
print("dylka krivky = ", dylkaKrivky([bod1, bod2], madlo1, madlo2, tolerance))
def pomocna(body, H1, H2, tolerance):
A = body[0]
E = body[-1]
if len(body) == 2: # prvni kolo
Ba = [(H1[0] - A[0])*0.25 + A[0], (H1[1] - A[1])*0.25 + A[1], (H1[2] - A[2])*0.25 + A[2]]
Bb = [(E[0] - H2[0])*0.25 + H2[0], (E[1] - H2[1])*0.25 + H2[1], (E[2] - H2[2])*0.25 + H2[2]]
B = [(Bb[0] - Ba[0])*0.25 + Ba[0], (Bb[1] - Ba[1])*0.25 + Ba[1], (Bb[2] - Ba[2])*0.25 + Ba[2]]
Ca = [(H1[0] - A[0])*0.5 + A[0], (H1[1] - A[1])*0.5 + A[1], (H1[2] - A[2])*0.5 + A[2]]
Cb = [(E[0] - H2[0])*0.5 + H2[0], (E[1] - H2[1])*0.5 + H2[1], (E[2] - H2[2])*0.5 + H2[2]]
C = [(Cb[0] - Ca[0])*0.5 + Ca[0], (Cb[1] - Ca[1])*0.5 + Ca[1], (Cb[2] - Ca[2])*0.5 + Ca[2]]
Da = [(H1[0] - A[0])*0.75 + A[0], (H1[1] - A[1])*0.75 + A[1], (H1[2] - A[2])*0.75 + A[2]]
Db = [(E[0] - H2[0])*0.75 + H2[0], (E[1] - H2[1])*0.75 + H2[1], (E[2] - H2[2])*0.75 + H2[2]]
D = [(Db[0] - Da[0])*0.75 + Da[0], (Db[1] - Da[1])*0.75 + Da[1], (Db[2] - Da[2])*0.75 + Da[2]]
body = [A, B, C, D, E]
dylkaPuvodni = velikostVektoru(A, E)
dylkaNova = velikostVektoru(A, B) + velikostVektoru(B, C) + velikostVektoru(C, D) + velikostVektoru(D, E)
if dylkaNova - tolerance < dylkaPuvodni:
return(body)
return(pomocna(body, H1, H2, tolerance))
else: # nasledujici kola
usecek = len(body) - 1
body2 = []
i = usecek
while i > 0:
t = 2*i - 1
bodA = [(H1[0] - A[0])*(t/(2*usecek)) + A[0], (H1[1] - A[1])*(t/(2*usecek)) + A[1], (H1[2] - A[2])*(t/(2*usecek)) + A[2]]
bodB = [(E[0] - H2[0])*(t/(2*usecek)) + H2[0], (E[1] - H2[1])*(t/(2*usecek)) + H2[1], (E[2] - H2[2])*(t/(2*usecek)) + H2[2]]
bod = [(bodB[0] - bodA[0])*(t/(2*usecek)) + bodA[0], (bodB[1] - bodA[1])*(t/(2*usecek)) + bodA[1], (bodB[2] - bodA[2])*(t/(2*usecek)) + bodA[2]]
body2.insert(0, bod) # predpojuju, protoze postupju odzadu
i -= 1
noveBody = []
i = 0
while i < usecek: # spojit body a body2 zipem
noveBody.append(body[i])
noveBody.append(body2[i])
i += 1
noveBody.append(body[usecek])
dylkaNova = 0
i = usecek
while i > 0:
dylkaNova += velikostVektoru(noveBody[2*i], noveBody[2*i - 1])
dylkaNova += velikostVektoru(noveBody[2*i - 1], noveBody[2*i - 2])
i -= 1
dylkaPuvodni = 0
i = int(usecek/2)
while i > 0:
dylkaPuvodni += velikostVektoru(body[2*i], body[2*i - 1])
dylkaPuvodni += velikostVektoru(body[2*i - 1], body[2*i - 2])
i -= 1
if dylkaNova - tolerance < dylkaPuvodni:
return(noveBody)
return(pomocna(noveBody, H1, H2, tolerance))
def bodNaKrivce(bod1, bod2, madlo1, madlo2, kde, tolerance):
dylkaCelkem = dylkaKrivky([bod1, bod2], madlo1, madlo2, tolerance)
dylkaBodu = dylkaCelkem*kde
sekvence = pomocna([bod1, bod2], madlo1, madlo2, tolerance)
celkovaDylka = 0
for i, s in enumerate(sekvence):
print("dylka pred ", celkovaDylka)
celkovaDylka += velikostVektoru(sekvence[i], sekvence[i + 1])
if celkovaDylka >= dylkaBodu:
print("skutecna dylka ", celkovaDylka)
return(sekvence[i + 1])
bod1 = [0, 0, 0]
bod2 = [1, 0, 0]
madlo1 = [0, 1, 0]
madlo2 = [1, 1, 0]
tolerance = 0.01
kde = 0.3
print("bod na krivce ", bodNaKrivce(bod1, bod2, madlo1, madlo2, kde, tolerance))
```