SwiftによるCGPathの曲線近似の実装（Inspired by paper.js）

PathSimplifier.swift

// PathSimplifier.swift

import UIKit

private let TOLERANCE: CGFloat = 1e-6
private let EPSILON: CGFloat = 1e-12

extension CGPoint {
    func add(p: CGPoint) -> CGPoint {
        return CGPointMake(x+p.x, y+p.y)
    }
    func add(n: CGFloat) -> CGPoint {
        return CGPointMake(x+n, y+n)
    }
    func subtract(p: CGPoint) -> CGPoint {
        return CGPointMake(x-p.x, y-p.y)
    }
    func subtract(n: CGFloat) -> CGPoint {
        return CGPointMake(x-n, y-n)
    }
    func multiply(p: CGPoint) -> CGPoint {
        return CGPointMake(x*p.x, y*p.y)
    }
    func multiply(n: CGFloat) -> CGPoint {
        return CGPointMake(x*n, y*n)
    }
    func divide(p: CGPoint) -> CGPoint {
        return CGPointMake(x/p.x, y/p.y)
    }
    func divide(n: CGFloat) -> CGPoint {
        return CGPointMake(x/n, y/n)
    }
    func negate() -> CGPoint {
        return CGPointMake(-x, -y)
    }
    func normalize() -> CGPoint {
        return normalize(1.0)
    }
    func normalize(d: CGFloat) -> CGPoint {
        let length = sqrt(x*x + y*y)
        let scale = length > 0 ? d / length : 0
        return CGPointMake(x*scale, y*scale)
    }
    func dot(p: CGPoint) -> CGFloat {
        return x*p.x + y*p.y
    }
    func distance(p: CGPoint) -> CGFloat {
        let dx = x - p.x
        let dy = y - p.y
        return sqrt(dx*dx + dy*dy)
    }
}

public class PathSimplifier: NSObject {
    public private(set) var curves: [[CGPoint]]
    
    private var path: CGMutablePathRef = CGPathCreateMutable()
    private var points: [CGPoint]
    private var error: CGFloat
    
    public convenience init(points: [CGPoint]) {
        self.init(points: points, error: 2.5)
    }
    
    public init(points: [CGPoint], error: CGFloat) {
        self.points = points
        self.error = error
        self.curves = [[CGPoint]]()
    }
    
    public func simplify() -> CGMutablePathRef {
        let length = points.count
        
        if (length > 1) {
            self.fitCubic(0, last: length-1, tan1: points[1].subtract(points[0]).normalize(), tan2: points[length-2].subtract(points.last!).normalize())
        }
        
        return path
    }
    
    // Fit a Bezier curve to a (sub)set of digitized points
    private func fitCubic(first: Int, last: Int, tan1: CGPoint, tan2: CGPoint) {
        if (last - first == 1) {
            let pt1 = points[first], pt2 = points[last], dist = pt1.distance(pt2) / 3
            if dist > 0 {
                self.addCurve([pt1, pt1.add(tan1.normalize(dist)), pt2.add(tan2.normalize(dist)), pt2])
            } else {
                // None.
            }
            return
        }
        var uPrime = self.chordLengthParameterize(first, last: last)
        var maxError = max(error, error * error)
        var split = 0
        var parametersInOrder = true
        for _ in 0...4 {
            var curve = self.generateBezier(first, last: last, uPrime: uPrime, tan1: tan1, tan2: tan2)
            //  Find max deviation of points to fitted curve
            let max = self.findMaxError(first, last: last, curve: &curve, u: uPrime)
            if (max.0 < error) {
                self.addCurve(curve)
                return
            }
            split = max.1
            // If error not too large, try reparameterization and iteration
            if (max.0 >= maxError && parametersInOrder) {
                break
            }
            parametersInOrder = self.reparameterize(first, last: last, u: &uPrime, curve: &curve)
            maxError = max.0
        }
        // Fitting failed -- split at max error point and fit recursively
        let V1 = points[split - 1].subtract(points[split])
        let V2 = points[split].subtract(points[split + 1])
        let tanCenter = V1.add(V2).divide(2).normalize()
        self.fitCubic(first, last: split, tan1: tan1, tan2: tanCenter)
        self.fitCubic(split, last: last, tan1: tanCenter.negate(), tan2: tan2)
    }
    
    private func addCurve(curve: [CGPoint]) {
        if (CGPathIsEmpty(path)) {
            CGPathMoveToPoint(path, nil, curve[0].x, curve[0].y)
        }
        CGPathAddCurveToPoint(path, nil, curve[1].x, curve[1].y, curve[2].x, curve[2].y, curve[3].x, curve[3].y)
        curves.append(curve)
    }
    
    // Use least-squares method to find Bezier control points for region.
    private func generateBezier(first: Int, last: Int, uPrime: [CGFloat], tan1: CGPoint, tan2: CGPoint) -> [CGPoint] {
        var epsilon: CGFloat = EPSILON
        let pt1: CGPoint = points[first]
        let pt2: CGPoint = points[last]
        // Create the C and X matrices
        var C: [[CGFloat]] = [[0, 0], [0, 0]]
        var X: [CGFloat] = [0, 0]
        
        for i in 0..<(last - first + 1) {
            let u = uPrime[i]
            let t = 1 - u
            let b = 3 * u * t
            let b0 = t * t * t
            let b1 = b * t
            let b2 = b * u
            let b3 = u * u * u
            let a1 = tan1.normalize(b1)
            let a2 = tan2.normalize(b2)
            let tmp = points[first + i].subtract(pt1.multiply(b0 + b1)).subtract(pt2.multiply(b2 + b3))
            C[0][0] += a1.dot(a1)
            C[0][1] += a1.dot(a2)
            // C[1][0] += a1.dot(a2)
            C[1][0] = C[0][1]
            C[1][1] += a2.dot(a2)
            X[0] += a1.dot(tmp)
            X[1] += a2.dot(tmp)
        }
        
        // Compute the determinants of C and X
        let detC0C1: CGFloat = C[0][0] * C[1][1] - C[1][0] * C[0][1]
        var alpha1: CGFloat = 0, alpha2: CGFloat = 0
        if (abs(detC0C1) > epsilon) {
            // Kramer's rule
            let detC0X = C[0][0] * X[1] - C[1][0] * X[0]
            let detXC1 = X[0] * C[1][1] - X[1] * C[0][1]
            // Derive alpha values
            alpha1 = detXC1 / detC0C1
            alpha2 = detC0X / detC0C1
        } else {
            // Matrix is under-determined, try assuming alpha1 == alpha2
            let c0 = C[0][0] + C[0][1]
            let c1 = C[1][0] + C[1][1]
            if (abs(c0) > epsilon) {
                alpha1 = X[0] / c0
                alpha2 = alpha1
            } else if (abs(c1) > epsilon) {
                alpha1 = X[1] / c1
                alpha2 = alpha1
            } else {
                // Handle below
                alpha1 = 0
                alpha2 = 0
            }
        }
        
        // If alpha negative, use the Wu/Barsky heuristic (see text)
        // (if alpha is 0, you get coincident control points that lead to
        // divide by zero in any subsequent NewtonRaphsonRootFind() call.
        let segLength: CGFloat = pt2.distance(pt1)
        let eps = epsilon * segLength
        var handle1, handle2 : CGPoint?
        if (alpha1 < eps || alpha2 < eps) {
            // fall back on standard (probably inaccurate) formula,
            // and subdivide further if needed.
            alpha1 = segLength / 3
            alpha2 = alpha1
        } else {
            
            // Check if the found control points are in the right order when
            // projected onto the line through pt1 and pt2.
            let line = pt2.subtract(pt1)
            // Control points 1 and 2 are positioned an alpha distance out
            // on the tangent vectors, left and right, respectively
            handle1 = tan1.normalize(alpha1)
            handle2 = tan2.normalize(alpha2)
            if (handle1!.dot(line) - handle2!.dot(line) > segLength * segLength) {
                // Fall back to the Wu/Barsky heuristic above.
                alpha1 = segLength / 3
                alpha2 = alpha1
                handle1 = nil
                handle2 = nil
            }
        }
        
        // First and last control points of the Bezier curve are
        // positioned exactly at the first and last data points
        // Control points 1 and 2 are positioned an alpha distance out
        // on the tangent vectors, left and right, respectively
        return [pt1, pt1.add(handle1 ?? tan1.normalize(alpha1)), pt2.add(handle2 ?? tan2.normalize(alpha2)), pt2]
    }
    
    // Given set of points and their parameterization, try to find
    // a better parameterization.
    private func reparameterize(first: Int, last: Int, inout u: [CGFloat], inout curve: [CGPoint]) -> Bool {
        for i in first...last {
            u[i - first] = self.findRoot(&curve, point: points[i], u: u[i - first])
        }
        // Detect if the new parameterization has reordered the points.
        // In that case, we would fit the points of the path in the wrong order.
        for i in 1..<u.count {
            if (u[i] <= u[i - 1]) {
                return false
            }
        }
        return true
    }
    
    // Use Newton-Raphson iteration to find better root.
    private func findRoot(inout curve: [CGPoint], point: CGPoint, u: CGFloat) -> CGFloat {
        var curve1: [CGPoint] = []
        var curve2: [CGPoint] = []
        // Generate control vertices for Q'
        for i in 0...2 {
            curve1.append(curve[i + 1].subtract(curve[i]).multiply(3))
        }
        // Generate control vertices for Q''
        for i in 0...1 {
            curve2.append(curve1[i + 1].subtract(curve1[i]).multiply(2))
        }
        // Compute Q(u), Q'(u) and Q''(u)
        let pt = self.evaluate(3, curve: curve, t: u)
        let pt1 = self.evaluate(2, curve: curve1, t: u)
        let pt2 = self.evaluate(1, curve: curve2, t: u)
        let diff = pt.subtract(point)
        let df = pt1.dot(pt1) + diff.dot(pt2)
        // Compute f(u) / f'(u)
        if (abs(df) < TOLERANCE) {
            return u
        }
        // u = u - f(u) / f'(u)
        return u - diff.dot(pt1) / df
    }
    
    // Evaluate a bezier curve at a particular parameter value
    private func evaluate(degree: Int, curve: [CGPoint], t: CGFloat) -> CGPoint {
        // Triangle computation
        var tmp: [CGPoint] = curve
        for i in 1...degree {
            for j in 0...degree-i {
                tmp[j] = tmp[j].multiply(1 - t).add(tmp[j + 1].multiply(t))
            }
        }
        return tmp[0]
    }
    
    // Assign parameter values to digitized points
    // using relative distances between points.
    private func chordLengthParameterize(first: Int, last: Int) -> [CGFloat] {
        var u: [CGFloat] = [CGFloat](count: last+1, repeatedValue: 0)
        for i in first+1...last {
            u[i - first] = u[i - first - 1] + points[i].distance(points[i - 1])
        }
        let m = last - first
        for i in 1...m {
            u[i] /= u[m]
        }
        return u
    }
    
    // Find the maximum squared distance of digitized points to fitted curve.
    private func findMaxError(first: Int, last: Int, inout curve: [CGPoint], u: [CGFloat]) -> (CGFloat, Int) {
        let half: Double = Double(last - first + 1) * 0.5
        var index: Int = Int(floor(half))
        var maxDist: CGFloat = 0
        for i in first+1..<last {
            let P = self.evaluate(3, curve: curve, t: u[i - first])
            let v = P.subtract(points[i])
            let dist = v.x * v.x + v.y * v.y // squared
            if (dist >= maxDist) {
                maxDist = dist
                index = i
            }
        }
        return (maxDist, index)
    }
    
}