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.
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.
2nd day edit 2: Fixed precision inaccuracy, made as foolproof as possible, translated into EN
from math import sqrt def lengthOfCurve(points, H1, H2, toleration): A = points[0] E = points[-1] if len(points) == 2: # round one 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 = lengthOfVector(A, E) newLength = lengthOfVector(A, B) + lengthOfVector(B, C) + lengthOfVector(C, D) + lengthOfVector(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 advance back first i -= 1 newPoints = [] i = 0 while i < numberOfSegments: # join points and points2 using zipper method newPoints.append(points[i]) newPoints.append(points2[i]) i += 1 newPoints.append(points[numberOfSegments]) newLength = 0 i = numberOfSegments while i > 0: newLength += lengthOfVector(newPoints[2*i], newPoints[2*i - 1]) newLength += lengthOfVector(newPoints[2*i - 1], newPoints[2*i - 2]) i -= 1 originalLength = 0 i = int(numberOfSegments/2) while i > 0: originalLength += lengthOfVector(points[2*i], points[2*i - 1]) originalLength += lengthOfVector(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 lengthOfVector(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("length of curve = ", lengthOfCurve([point1, point2], handle1, handle2, toleration)) def findPoints(points, H1, H2, toleration, lengthOfCurv): A = points[0] E = points[-1] if len(points) == 2: # round one 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] longestSegment = lengthOfVector(A, B) lenBC = lengthOfVector(B, C) lenCD = lengthOfVector(C, D) lenDE = lengthOfVector(D, E) if lenBC > longestSegment: longestSegment = lenBC if lenCD > longestSegment: longestSegment = lenBC if lenDE > longestSegment: longestSegment = lenBC if longestSegment < toleration: return(points) return(findPoints(points, H1, H2, toleration, lengthOfCurv)) 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 advance back first i -= 1 newPoints = [] i = 0 while i < numberOfSegments: # join points and points2 using zipper method newPoints.append(points[i]) newPoints.append(points2[i]) i += 1 newPoints.append(points[numberOfSegments]) newLength = 0 i = numberOfSegments while i > 0: newLength += lengthOfVector(newPoints[2*i], newPoints[2*i - 1]) newLength += lengthOfVector(newPoints[2*i - 1], newPoints[2*i - 2]) i -= 1 originalLength = 0 i = int(numberOfSegments/2) longestSegment = 0 while i > 0: byHowMuch = lengthOfVector(points[2*i], points[2*i - 1]) if byHowMuch > longestSegment: longestSegment = byHowMuch originalLength += byHowMuch byHowMuch = lengthOfVector(points[2*i - 1], points[2*i - 2]) if byHowMuch > longestSegment: longestSegment = byHowMuch originalLength += byHowMuch i -= 1 #print(longestSegment) if longestSegment < toleration: return(newPoints) return(findPoints(newPoints, H1, H2, toleration, lengthOfCurv)) def pointOnCurve(point1, point2, handle1, handle2, placement, toleration, lengthOfCurv): aggregateLength = lengthOfCurve([point1, point2], handle1, handle2, toleration) lengthOfUnfinishedCurve = aggregateLength*placement sequence = findPoints([point1, point2], handle1, handle2, toleration, lengthOfCurv) aggregateLength = 0 for i, s in enumerate(sequence): aggregateLength += lengthOfVector([sequence[i][0], sequence[i][1], sequence[i][2]], [sequence[i + 1][0], sequence[i + 1][1], sequence[i + 1][2]]) if aggregateLength >= lengthOfUnfinishedCurve: return([round(sequence[i + 1][0]/toleration)*toleration, round(sequence[i + 1][1]/toleration)*toleration, round(sequence[i + 1][2]/toleration)*toleration]) point1 = [0, 0, 0] point2 = [1, 0, 0] handle1 = [0, 1, 0] handle2 = [1, 1, 0] toleration = 0.01 placement = 0.3 lengthOfCurv = lengthOfCurve([point1, point2], handle1, handle2, toleration) print("point on curve = ", pointOnCurve(point1, point2, handle1, handle2, placement, toleration, lengthOfCurv))