internal.py

# SPDX-License-Identifier: GPL-2.0-or-later

import bpy, math, cmath
from mathutils import Vector, Matrix
from collections import namedtuple

units = [
    ('-', 'None', '1.0', 0),
    ('px', 'Pixel', '1.0', 1),
    ('m', 'Meter', '1.0', 2),
    ('dm', 'Decimeter', '0.1', 3),
    ('cm', 'Centimeter', '0.01', 4),
    ('mm', 'Millimeter', '0.001', 5),
    ('yd', 'Yard', '0.9144', 6),
    ('ft', 'Foot', '0.3048', 7),
    ('in', 'Inch', '0.0254', 8)
]

param_tolerance = 0.0001
AABB = namedtuple('AxisAlignedBoundingBox', 'center dimensions')
Plane = namedtuple('Plane', 'normal distance')
Circle = namedtuple('Circle', 'orientation center radius')

def circleOfTriangle(a, b, c):
    # https://en.wikipedia.org/wiki/Circumscribed_circle#Cartesian_coordinates_from_cross-_and_dot-products
    dirBA = a-b
    dirCB = b-c
    dirAC = c-a
    normal = dirBA.cross(dirCB)
    lengthBA = dirBA.length
    lengthCB = dirCB.length
    lengthAC = dirAC.length
    lengthN = normal.length
    if lengthN == 0:
        return None
    factor = -1/(2*lengthN*lengthN)
    alpha = (dirBA@dirAC)*(lengthCB*lengthCB*factor)
    beta = (dirBA@dirCB)*(lengthAC*lengthAC*factor)
    gamma = (dirAC@dirCB)*(lengthBA*lengthBA*factor)
    center = a*alpha+b*beta+c*gamma
    radius = (lengthBA*lengthCB*lengthAC)/(2*lengthN)
    tangent = (a-center).normalized()
    orientation = Matrix.Identity(3)
    orientation.col[2] = normal/lengthN
    orientation.col[1] = (a-center).normalized()
    orientation.col[0] = orientation.col[1].xyz.cross(orientation.col[2].xyz)
    return Circle(orientation=orientation, center=center, radius=radius)

def circleOfBezier(points, tolerance=0.000001, samples=16):
    circle = circleOfTriangle(points[0], bezierPointAt(points, 0.5), points[3])
    if circle == None:
        return None
    variance = 0
    for t in range(0, samples):
        variance += ((circle.center-bezierPointAt(points, (t+1)/(samples-1))).length/circle.radius-1) ** 2
    variance /= samples
    return None if variance > tolerance else circle

def areaOfPolygon(vertices):
    area = 0
    for index, current in enumerate(vertices):
        prev = vertices[index-1]
        area += (current[0]+prev[0])*(current[1]-prev[1])
    return area*0.5

def linePointDistance(begin, dir, point):
    return (point-begin).cross(dir.normalized()).length

def linePlaneIntersection(origin, dir, plane):
    det = dir@plane.normal
    return float('nan') if det == 0 else (plane.distance-origin@plane.normal)/det

def nearestPointOfLines(originA, dirA, originB, dirB, tolerance=0.0):
    # https://en.wikipedia.org/wiki/Skew_lines#Nearest_Points
    normal = dirA.cross(dirB)
    normalA = dirA.cross(normal)
    normalB = dirB.cross(normal)
    divisorA = dirA@normalB
    divisorB = dirB@normalA
    if abs(divisorA) <= tolerance or abs(divisorB) <= tolerance:
        return (float('nan'), float('nan'), None, None)
    else:
        paramA = (originB-originA)@normalB/divisorA
        paramB = (originA-originB)@normalA/divisorB
        return (paramA, paramB, originA+dirA*paramA, originB+dirB*paramB)

def lineSegmentLineSegmentIntersection(beginA, endA, beginB, endB, tolerance=0.001):
    dirA = endA-beginA
    dirB = endB-beginB
    paramA, paramB, pointA, pointB = nearestPointOfLines(beginA, dirA, beginB, dirB)
    if math.isnan(paramA) or (pointA-pointB).length > tolerance or \
       paramA < 0 or paramA > 1 or paramB < 0 or paramB > 1:
        return None
    return (paramA, paramB, pointA, pointB)

def aabbOfPoints(points):
    min = Vector(points[0])
    max = Vector(points[0])
    for point in points:
        for i in range(0, 3):
            if min[i] > point[i]:
                min[i] = point[i]
            if max[i] < point[i]:
                max[i] = point[i]
    return AABB(center=(max+min)*0.5, dimensions=(max-min)*0.5)

def aabbIntersectionTest(a, b, tolerance=0.0):
    for i in range(0, 3):
        if abs(a.center[i]-b.center[i]) > a.dimensions[i]+b.dimensions[i]+tolerance:
            return False
    return True

def isPointInAABB(point, aabb, tolerance=0.0, ignore_axis=None):
    for i in range(0, 3):
        if i != ignore_axis and (point[i] < aabb.center[i]-aabb.dimensions[i]-tolerance or point[i] > aabb.center[i]+aabb.dimensions[i]+tolerance):
            return False
    return True

def lineAABBIntersection(lineBegin, lineEnd, aabb):
    intersections = []
    for i in range(0, 3):
        normal = [0, 0, 0]
        normal = Vector(normal[0:i] + [1] + normal[i+1:])
        for j in range(-1, 2, 2):
            plane = Plane(normal=normal, distance=aabb.center[i]+j*aabb.dimensions[i])
            param = linePlaneIntersection(lineBegin, lineEnd-lineBegin, plane)
            if param < 0 or param > 1 or math.isnan(param):
                continue
            point = lineBegin+param*(lineEnd-lineBegin)
            if isPointInAABB(point, aabb, 0.0, i):
                intersections.append((param, point))
    return intersections

def bezierPointAt(points, t):
    s = 1-t
    return s*s*s*points[0] + 3*s*s*t*points[1] + 3*s*t*t*points[2] + t*t*t*points[3]

def bezierTangentAt(points, t):
    s = 1-t
    return s*s*(points[1]-points[0])+2*s*t*(points[2]-points[1])+t*t*(points[3]-points[2])
    # return s*s*points[0] + (s*s-2*s*t)*points[1] + (2*s*t-t*t)*points[2] + t*t*points[3]

def bezierLength(points, beginT=0, endT=1, samples=1024):
    # https://en.wikipedia.org/wiki/Arc_length#Finding_arc_lengths_by_integrating
    vec = [points[1]-points[0], points[2]-points[1], points[3]-points[2]]
    dot = [vec[0]@vec[0], vec[0]@vec[1], vec[0]@vec[2], vec[1]@vec[1], vec[1]@vec[2], vec[2]@vec[2]]
    factors = [
        dot[0],
        4*(dot[1]-dot[0]),
        6*dot[0]+4*dot[3]+2*dot[2]-12*dot[1],
        12*dot[1]+4*(dot[4]-dot[0]-dot[2])-8*dot[3],
        dot[0]+dot[5]+2*dot[2]+4*(dot[3]-dot[1]-dot[4])
    ]
    # https://en.wikipedia.org/wiki/Trapezoidal_rule
    length = 0
    prev_value = math.sqrt(factors[4]+factors[3]+factors[2]+factors[1]+factors[0])
    for index in range(0, samples+1):
        t = beginT+(endT-beginT)*index/samples
        # value = math.sqrt(factors[4]*(t**4)+factors[3]*(t**3)+factors[2]*(t**2)+factors[1]*t+factors[0])
        value = math.sqrt((((factors[4]*t+factors[3])*t+factors[2])*t+factors[1])*t+factors[0])
        length += (prev_value+value)*0.5
        prev_value = value
    return length*3/samples

# https://en.wikipedia.org/wiki/Root_of_unity
# cubic_roots_of_unity = [cmath.rect(1, i/3*2*math.pi) for i in range(0, 3)]
cubic_roots_of_unity = [complex(1, 0), complex(-1, math.sqrt(3))*0.5, complex(-1, -math.sqrt(3))*0.5]
def bezierRoots(dists, tolerance=0.0001):
    # https://en.wikipedia.org/wiki/Cubic_function
    # y(t) = a*t^3 +b*t^2 +c*t^1 +d*t^0
    a = 3*(dists[1]-dists[2])+dists[3]-dists[0]
    b = 3*(dists[0]-2*dists[1]+dists[2])
    c = 3*(dists[1]-dists[0])
    d = dists[0]
    if abs(a) > tolerance: # Cubic
        E2 = a*c
        E3 = a*a*d
        A = (2*b*b-9*E2)*b+27*E3
        B = b*b-3*E2
        C = ((A+cmath.sqrt(A*A-4*B*B*B))*0.5) ** (1/3)
        roots = []
        for root in cubic_roots_of_unity:
            root *= C
            root = -1/(3*a)*(b+root+B/root)
            if abs(root.imag) < tolerance and root.real > -param_tolerance and root.real < 1.0+param_tolerance:
                roots.append(max(0.0, min(root.real, 1.0)))
        # Remove doubles
        roots.sort()
        for index in range(len(roots)-1, 0, -1):
            if abs(roots[index-1]-roots[index]) < param_tolerance:
                roots.pop(index)
        return roots
    elif abs(b) > tolerance: # Quadratic
        disc = c*c-4*b*d
        if disc < 0:
            return []
        disc = math.sqrt(disc)
        return [(-c-disc)/(2*b), (-c+disc)/(2*b)]
    elif abs(c) > tolerance: # Linear
        root = -d/c