Calculating the exact length of an SVG Arc in Python?

Question

I want to be able to calculate the exact length of an SVG Arc. I can do all the manipulations rather easily. But, I am unsure of whether there is a solution at all or the exact implementation of the solution.

Here's the exact solution for the circumference of the ellipse. Using popular libraries is fine. I fully grasp that there are no easy solution as they will all require hypergeometric functions to be exact.

from scipy import pi, sqrt
from scipy.special import hyp2f1

def exact(a, b):
    t = ((a - b) / (a + b)) ** 2
    return pi * (a + b) * hyp2f1(-0.5, -0.5, 1, t)

a = 2.667950e9
b = 6.782819e8
print(exact(a, b))

My idea is to have this as opt-in code if you happen to have scipy installed it'll use the exact super-fancy solution, else it'll fall back to the weaker approximation code (progressively smaller line segments until error is small). The problem is the math level here is above me. And I don't know if there's some ways to specify a start and stop point for that.

Most of the approximation solutions are for ellipses, but I only want the arc. There may also be a solution unknown to me, for calculating the length of an arc on an ellipse but since the start and end position can be anywhere. It doesn't seem to be instantly viable to say a sweep angle is 15% of the total possible angle therefore it's 15% of the ellipse circumference.

A more effective less fancy arc approximation might also be nice. There are progressively better ellipse approximations but I can't go from ellipse circumference to arc length, so those are currently not helpful.

---

Let's say the arc parameterization is the start and end points on the ellipse. Since that's how SVG is parameterized. But, anything that isn't tautological like arc_length parameterization is a correct answer.

Answer 1

If you want to calculate this with your bare hands and the std lib, you can base your calculation on the following formula. This is only valid for two points on the upper half of the ellipse because of the acos but we're going to use it with the angles directly.

The calculation consists in these steps:

1. Start with the SVG data: start point, a ,b rotation, long arc, sweep, end point
2. Rotate the coordinate system to match the horizontal axis of the ellipse.
3. Solve a system of 4 equations with 4 unknowns to get the centre point and the angles corresponding to the start and end point
4. Approximate the integral by a discreet sum over small segments. This is where you could use scipy.special.ellipeinc, as suggested in the comments.

Step 2 is easy, just use a rotation matrix (note the angle rot is positive in the clockwise direction):

 m = [
        [math.cos(rot), math.sin(rot)], 
        [-math.sin(rot), math.cos(rot)]
   ]

Step 3 is very well explained in this answer. Note the value obtained for a1 is modulo pi because it is obtained with atan. That means that you need to calculate the centre points for the two angles t1 and t2 and check they match. If they don't, add pi to a1 and check again.

Step 4 is quite straightforward. Divide the interval [t1, t2] into n segments, get the value of the function at the end of each segment, time it by the segment length and sum all this up. You can try refining this by taking the value of the function at the mid-point of each segment, but I'm not sure there is much benefit to that. The number of segments is likely to have more effect on the precision.

Here is a very rough Python version of the above (please bear with the ugly coding style, I was doing this on my mobile whilst traveling )

import math

PREC = 1E-6

# matrix vector multiplication
def transform(m, p):
    return ((sum(x * y for x, y in zip(m_r, p))) for m_r in m)

# the partial integral function        
def ellipse_part_integral(t1, t2, a, b, n=100):

    # function to integrate
    def f(t):
        return math.sqrt(1 - (1 - a**2 / b**2) * math.sin(t)**2)


    start = min(t1, t2)
    seg_len = abs(t1 - t2) / n
    return - b * sum(f(start + seg_len * (i + 1)) * seg_len for i in range(n))


def ellipse_arc_length(x1, y1, a, b, rot, large_arc, sweep, x2, y2):
    if abs(x1 - x2) < PREC and abs(y1 - y2) < PREC:
        return 0

    # get rot in radians
    rot = math.pi / 180 * rot
    # get the coordinates in the rotated coordinate system
    m = [
        [math.cos(rot), math.sin(rot)], 
        [- math.sin(rot), math.cos(rot)]
    ]
    x1_loc, y1_loc, x2_loc, y2_loc = *transform(m, (x1,y1)), *transform(m, (x2,y2))

    r1 = (x1_loc - x2_loc) / (2 * a)
    r2 = (y2_loc - y1_loc) / (2 * b)

    # avoid division by 0 if both points have same y coord
    if abs(r2) > PREC:
        a1 = math.atan(r1 / r2)
    else:
        a1 = r1 / abs(r1) * math.pi / 2

    if abs(math.cos(a1)) > PREC:
        a2 = math.asin(r2 / math.cos(a1))
    else:
        a2 = math.asin(r1 / math.sin(a1))

    # calculate the angle of start and end point
    t1 = a1 + a2
    t2 = a1 - a2

    # calculate centre point coords
    x0 = x1_loc - a * math.cos(t1)
    y0 = y1_loc - b * math.sin(t1)

    x0s = x2_loc - a * math.cos(t2)
    y0s = y2_loc - b * math.sin(t2)


    # a1 value is mod pi so the centres may not match
    # if they don't, check a1 + pi
    if abs(x0 - x0s) > PREC or abs(y0 - y0s) > PREC:
        a1 = a1 + math.pi
        t1 = a1 + a2
        t2 = a1 - a2

        x0 = x1_loc - a * math.cos(t1)
        y0 = y1_loc - b * math.sin(t1)

        x0s = x2_loc - a * math.cos(t2)
        y0s = y2_loc - b * math.sin(t2)

    # get the angles in the range [0, 2 * pi]
    if t1 < 0:
        t1 += 2 * math.pi
    if t2 < 0:
        t2 += 2 * math.pi

    # increase minimum by 2 * pi for a large arc
    if large_arc:
        if t1 < t2:
            t1 += 2 * math.pi 
        else:
            t2 += 2 * math.pi

    return ellipse_part_integral(t1, t2, a, b)

print(ellipse_arc_length(0, 0, 40, 40, 0, False, True, 80, 0))

The good news is that the sweep flag doesn't matter as long as you're just looking for the length of the arc.

I'm not 100% sure the modulo pi problem is handled correctly and the implementation above may have a few bugs. Nevertheless, it gave me a good approximation of the length in the simple case of a half circle, so I dare calling it WIP. Let me know if this is worth pursuing, I can have a further look when I'll be seated at a computer. Or maybe someone can come up with a clean way of doing this in the meantime?