@implementation NSBezierPath (MyBezierPathAdditions) /* Solve cubic equation ax^3 + bx^2 + cx + d = 0 by the Cardano formula */ /* If all three solutions are real numbers, 1 is returned with s[0], s[1], s[2] as the solutions. */ /* If two are conjugated complex numbers, 0 is returned with s[0] as the only real solution, and s[1] and s[2] as the real and imaginary part of one of the complex solutions (so that the complex solutions are s[1] + i * s[2] and s[1] - i * s[2]). */ static int sSolveCubicEquation(double a, double b, double c, double d, double *s) { double p, h, k, w; p = c / (a * 3.0) - b * b / (a * a * 9.0); h = -d / (a * 2.0) - b * b * b / (a * a * a * 27.0) + b * c / (a * a * 6.0); k = h * h + p * p * p; w = b / (-3.0 * a); if (k > 0.0) { double m, n; m = cbrt(h + sqrt(k)); /* cbrt(): cubic root function */ n = cbrt(h - sqrt(k)); s[0] = w + m + n; s[1] = w - (m + n) * 0.5; s[2] = 0.86602540378443859659 * (m - n); /* sqrt(3) * 0.5 */ return 0; } else { double r, th, cs, sn; r = pow(h * h - k, 1.0 / 6.0); th = atan2(sqrt(-k), h) / 3.0; cs = cos(th); sn = sin(th); s[0] = w + 2 * r * cs; s[1] = w + r * (-cs + 1.7320508075688771932 * sn); /* sqrt(3) */ s[2] = w + r * (-cs - 1.7320508075688771932 * sn); return 1; } } /* Get the crossing points of an NSBezierPath and a line segment */ /* The return value is n + a (n is an integer and a is a float with 0 <= a < 1), where n is the index of the bezier path element and a is the parameter indicating the position of the crossing point within the path element. This value is for the first found crossing point. If no crossing points are present, -1.0 is returned. In addition, if cptr and nptr are non-NULL, then *cptr contains a malloc()'ed array of the coordinates of *all* crossing points (as NSPoint's) and *nptr contains the number of returned points. */ - (float)crossLineSegmentWithStartPoint:(NSPoint)spoint endPoint:(NSPoint)epoint crossPoints:(NSPoint **)cptr count:(int *)nptr { NSPoint points[4], *cpoints; int base, i, c, ncpoints, nccount; float s, t, retval; base = 0; c = [self elementCount]; points[0] = NSZeroPoint; if (cptr != NULL && nptr == NULL) cptr = NULL; /* Both cptr and nptr should be non-NULL if crossing points are to be obtained */ nccount = ncpoints = 0; cpoints = NULL; retval = -1.0; for (i = 0; i < c; i++) { int type = [self elementAtIndex:i associatedPoints:points + 1]; if (type == NSMoveToBezierPathElement) { points[0] = points[1]; continue; } else if (type == NSLineToBezierPathElement) { // Line segment float a11, a12, a21, a22, d1, d2, det; d1 = spoint.x - points[0].x; d2 = spoint.y - points[0].y; a11 = points[1].x - points[0].x; a21 = points[1].y - points[0].y; a12 = spoint.x - epoint.x; a22 = spoint.y - epoint.y; det = (a11 * a22 - a12 * a21); s = (a22 * d1 - a12 * d2) / det; t = (-a21 * d1 + a11 * d2) / det; if (s >= 0 && s <= 1 && t >= 0 && t <= 1) { if (cptr == NULL) return (float)base + s; if (retval < 0) retval = (float)base + s; if (nccount >= ncpoints) { ncpoints += 8; cpoints = (NSPoint *)realloc(cpoints, sizeof(NSPoint) * ncpoints); } cpoints[nccount++] = NSMakePoint(spoint.x + t * (epoint.x - spoint.x), spoint.y + t * (epoint.y - spoint.y)); } points[0] = points[1]; base++; } else if (type == NSCurveToBezierPathElement) { // Curve segment double a, b, c, a3, b3, c3, d3; double sol[3]; int j, n, ybase; a = epoint.y - spoint.y; b = spoint.x - epoint.x; c = spoint.y * epoint.x - spoint.x * epoint.y; a3 = a * (-points[0].x + 3 * points[1].x - 3 * points[2].x + points[3].x) + b * (-points[0].y + 3 * points[1].y - 3 * points[2].y + points[3].y); b3 = 3 * a * (points[0].x - 2 * points[1].x + points[2].x) + 3 * b * (points[0].y - 2 * points[1].y + points[2].y); c3 = 3 * a * (-points[0].x + points[1].x) + 3 * b * (-points[0].y + points[1].y); d3 = c + a * points[0].x + b * points[0].y; // solve a*t^3 + b*t^2 + c*t + d = 0 for t if (a3 < -1e-8 || a3 > 1e-8) { n = sSolveCubicEquation(a3, b3, c3, d3, sol); n = (n ? 3 : 1); } else if (b3 < -1e-8 || b3 > 1e-8) { /* Quadratic equation */ s = c3 * c3 - 4.0 * b3 * d3; if (s < 0.0) n = 0; /* No real solution */ else { n = 2; s = sqrt(s); sol[0] = (-c3 + s) / (2.0 * b3); sol[1] = (-c3 - s) / (2.0 * b3); } } else if (c3 < -1e-8 || c3 > 1e-8) { /* Linear equation */ n = 1; sol[0] = -d3 / c3; } else n = 0; ybase = (fabs(a) > fabs(b)); /* Flag to avoid zero division */ for (j = 0; j < n; j++) { NSPoint p; /* Crossing point */ t = sol[j]; if (ybase) { p.y = (-points[0].y + 3 * points[1].y - 3 * points[2].y + points[3].y) * t * t * t + 3 * (points[0].y - 2 * points[1].y + points[2].y) * t * t + 3 * (-points[0].y + points[1].y) * t + points[0].y; s = (p.y - spoint.y) / (epoint.y - spoint.y); p.x = s * (epoint.x - spoint.x) + spoint.x; } else { p.x = (-points[0].x + 3 * points[1].x - 3 * points[2].x + points[3].x) * t * t * t + 3 * (points[0].x - 2 * points[1].x + points[2].x) * t * t + 3 * (-points[0].x + points[1].x) * t + points[0].x; s = (p.x - spoint.x) / (epoint.x - spoint.x); p.y = s * (epoint.y - spoint.y) + spoint.y; } if (s >= 0 && s <= 1 && t >= 0 && t <= 1) { if (cptr == NULL) return (float)base + t; if (retval < 0) retval = (float)base + t; if (nccount >= ncpoints) { ncpoints += 8; cpoints = (NSPoint *)realloc(cpoints, sizeof(NSPoint) * ncpoints); } cpoints[nccount++] = p; } } points[0] = points[3]; base++; } } if (cptr != NULL) { *nptr = nccount; *cptr = cpoints; } return retval; } @end